General theory of the plasmoid instability

The rapid conversion of magnetic energy into plasma particle energy through the process of magnetic reconnection is of great importance in the realm of plasma physics and astrophysics.^1–4 Sawtooth crashes, magnetospheric substorms, stellar and gamma-ray flares are just a few examples of phenomena in which magnetic reconnection plays an essential role.

In large systems, such as those found in space and astrophysical environments, the potential formation of highly elongated current sheets would result in extremely low reconnection rates, which fail to account for the observed fast energy release rates.^5–7 However, such current sheets are subject to a violent linear instability that leads to their breakup, giving rise to a tremendous increase in the reconnection rate that appears to be very weakly dependent on the Lundquist number of the system in the nonlinear regime.^8–17 This crucial instability, which serves as a trigger of fast reconnection, is the plasmoid instability,² thus dubbed as it leads to the formation of plasmoids.

In the widely studied Sweet-Parker current sheets, which are characterized by an inverse aspect ratio $a/L∼S−1/2$⁠, Tajima and Shibata,¹ as well as Loureiro et al.,¹⁸ have found that the growth rate γ and the wavenumber k of the plasmoid instability obey $γτA∼S1/4$ and kL ∼ S^3∕8, where τ_A is the Alfvénic timescale based on the length of the current sheet. Since the Lundquist number S is extremely large in most space and astrophysical plasmas,¹⁹ the linear growth of the instability turns out to be surprisingly fast, and the number of plasmoids produced is also very high. Other notable works have since followed, which have verified and extended the work on the plasmoid instability in different contexts.^20–24 

Despite the success of the theory, its limitations soon became evident. For sufficiently high growth rates, Sweet-Parker current sheets cannot be attained as current layers are linearly unstable and disrupt before this state is achieved. In order to bypass this limitation, Pucci and Velli²⁵ conjectured that current sheets break up when γτ_A ∼ 1. Later, Uzdensky and Loureiro²⁶ considered a similar criterion (γτ = 1) as the end-point of the linear stage of the instability, presenting an appealing but heuristic discussion for the case of a current sheet evolving on the timescale τ. However, as demonstrated below, the most interesting physics occurs in the regime γτ > 1. Yet, a quantitative theory that encompasses all of the intricate physics underlying the onset and development of the plasmoid instability has been elusive. It is the purpose of this letter to develop a quantitative theory of the plasmoid instability in time evolving current sheets based on a principle of least time. We obtain new and surprising results, such as scalings that are not power laws, and provide a detailed picture of how the plasmoid instability arises and develops.

In a time evolving current sheet, tearing modes become unstable at different times and exhibit different instantaneous growth rates γ(k, t). Their amplitude changes in time according to $ψ(k,t)=ψ0 exp(∫t0tγ(t′)dt′)$⁠, where ψ represents the tearing eigenfunction and $ψ0:=ψ(k,t0)$⁠. During the linear evolution, the amplitudes of the different modes are small, thereby ensuring that they do not affect each other and the current sheet evolution. Their linear evolution ends when the plasmoid half-width² 

grows to the same order of the inner resistive layer width

Note that ψ is now taken to be implicitly evaluated at the resonant surface, B₀ is the reconnecting magnetic field upstream of the current sheet, and a is the current sheet half-width. Here, w may be regarded as a label of the perturbation amplitude and represents the plasmoid width only if the associated mode is dominant with respect to the others (or, less restrictively, for multiple dominant modes, as long as they are few in number and sufficiently localized in the spectrum). This does turn out to be the case at the end of the linear phase.

In such a complex scenario, we intend to determine the mode (k_*) that emerges from the linear phase and transitions into the nonlinear phase, and the time (t_*) at which this transition occurs. The most general way to address this question arguably relies on the formulation of a principle of least time for the plasmoid instability, i.e., the mode of the plasmoid instability that emerges from the linear phase is the one that traverses it in the least time. Mathematically, to implement our formulation, we introduce the function

It is self-evident that $w(k,t0)≪δin(k,t0)$ is necessary for the linear evolution to exist, viz., the initial mode amplitude has to be sufficiently small in large-S plasmas. We assume a continuous spectrum of wavenumbers k. Then, the least time principle yields the mathematical relations

The second relation is equivalent to the condition $dt*/dk*=0$⁠, where $t*=t*(k*)$ implicitly follows from the first relation. This allows us to interpret t_* as the variable that is extremized. In addition, it can also be shown a posteriori that $|(k/γ)(∂γ/∂k)|≪1$⁠, and consequently, $δin∝k−1/2$ holds true in the neighborhood of k_*. Since w is localized and has a much stronger dependence than δ_in on k, the mode that completes the linear phase in the least time is the dominant (larger amplitude) one that enters the nonlinear phase.

In our subsequent discussion, it is convenient to work with dimensionless quantities. In particular, we normalize all the lengths to the current sheet half-length L, the magnetic field to the upstream field B₀, and the time to the Alfvén time τ_A = L/v_A. We use carets to indicate the dimensionless quantities – $â=a/L, δ̂in=δin/L, k̂=kL, ψ̂=ψ/LB0, t̂=t/τA$⁠, and $γ̂=γτA$⁠. Therefore, the normalized magnetic diffusivity corresponds to the inverse of the Lundquist number, i.e., $η̂−1=S:=vAL/η$⁠.

Although our framework is altogether general, we are interested in the case where L and B₀ remain approximately constant, while the current sheet width decreases in time via $â(t̂)=â0f(t̂)$⁠. Here, $f(t̂)$ is a function that must obey $f(t̂0)=1$ and $limt̂→∞f(t̂)=S−1/2$⁠. Indeed, $â=S−1/2$ is the natural lower limit to the thickness of a reconnection layer, due to the increase in the Ohmic heating when B₀/a increases.³ 

Now, we can explicitly rewrite the first of Eq. (4) as

where $ŵ0:=ŵ(k̂,t̂0)=2(ψ̂0â0)1/2$ is the label for the initial perturbation amplitude. We assume that the initial perturbation is the same for all wavelengths, but other possibilities can be easily handled with our treatment. Then, we combine the second of Eq. (4) with Eq. (5) to obtain the least time equation

This equation leads us to $k̂*, γ̂(k̂*,t̂*):=γ̂*$ and $δ̂in(k̂*,t̂*):=δ̂in*$ as a function of the inverse aspect ratio $â(t̂*):=â*$⁠.

The procedure delineated above is fairly general, and can be applied to tearing modes with arbitrary degrees of collisionality. However, to proceed further, we must specify an expression for $γ̂$⁠. Here, for the sake of definiteness, we focus on resistive tearing modes. When the evolution of the current sheet width is slow, i.e., $â−1dâ/dt̂<γ̂$⁠, the growth rate of the tearing mode can be computed using the instantaneous value of $â$⁠. For resistive tearing modes, two simple algebraic relations exist, which are valid in two different regimes, depending on the value of the tearing stability parameter $Δ̂′$⁠.²⁷ In the small-$Δ̂′$ regime, defined as $Δ̂′δ̂in≪1$⁠, the classic analysis by Furth et al.²⁷ demonstrated that

where $cΓ=[(2π)−1Γ(1/4)/Γ(3/4)]4/5≈0.55$⁠. On the other hand, in the large-$Δ̂′$ regime, defined as $Δ̂′δ̂in≳1$⁠, Coppi et al.²⁸ showed that the growth rate becomes independent of $Δ̂′$⁠, resulting in

As we are interested in the entire domain of $Δ̂′$⁠, we seek an expression for $γ̂$ that (i) is a reasonable approximation of the exact growth rate,²⁸ (ii) reduces to (7) and (8) in the appropriate limits, and (iii) is simple enough to be analytically tractable. For this purpose, we consider the half-harmonic mean of this two relations, i.e.,

Since the harmonic mean is a Schur-concave function,²⁹ in addition to being quite simple, we have verified that it fulfills all of the criteria described above. Furthermore, choosing a different approximation for $γ̂$ such as the simpler one employed in Refs. 1, 8, 15, 23, and 25 leads to the same scaling relations presented below, albeit with slightly different numerical factors.

Considering a common Harris-like current sheet, it is known that²⁷ $Δ̂′â=2[(k̂â)−1−k̂â]$⁠. As we are not interested in the very slow-growing part of the mode evolution, we consider the regime $k̂â≪1$⁠. Then, without any further simplification, it is possible to obtain the following expression from Eq. (6):

A careful consideration of this equation reveals that the two terms on the right-hand-side must approximately balance each other. Therefore, we end up with

where c_k is a $O(1)$ coefficient. Then, using Eq. (11), we obtain the expressions

where c_γ and c_δ are also $O(1)$ coefficients. These relations show that the dominant mode at the end of the linear phase exhibits the same scaling properties of the fastest growing mode,²⁷ the latter of which satisfies the equation $∂γ̂/∂k̂|k̂f*,t̂*=0$⁠. This property arises because the additional contributions in Eq. (6) can be shown to be negligible if $ŵ0≪δ̂in(k̂*,t̂0)$⁠. If the current layer has the time to evolve until the Sweet-Parker inverse aspect ratio (⁠$â*→S−1/2$⁠) is attained, then we obtain

which exactly matches previous studies of the plasmoid instability^1,8,18,22,24 that were undertaken assuming a fixed Sweet-Parker current sheet.

However, for very high Lundquist numbers, the plasmoids complete their linear evolution well before the Sweet-Parker aspect ratio is reached. Therefore, we have to evaluate $â*$ for a more general case. This can be done by substituting the relations (11) and (12) into Eq. (5), which give us the following inverse-aspect-ratio equation:

This expression yields the final inverse aspect ratio $â*$ for a general current sheet evolution $â(t̂)$⁠. It shows clearly that $â*$⁠, and consequently the scaling relations of $t̂*, γ̂*, k̂*$ and $δ̂in*$⁠, cannot be universal, as they must depend on the specific form of the function $â(t̂)$⁠.

We proceed further by considering, arguably, the most common case of current sheet thinning – the exponential thinning – which is typically the result of an instability-driven current sheet. Subsequently, we generalize the exponential thinning to encompass even algebraic cases. Other, more uncommon, possibilities could also be investigated, since the developed framework is fairly general.

We consider a current sheet width that shrinks in time as

where $â∞=S−1/2$⁠. Note that the plasmoid width (1) starts to grow only when $γ̂>1/τ$⁠. For $t̂*>τ ln(1+â0S1/2)$⁠, the current sheet approaches $â*≃S−1/2$⁠. Therefore, one can easily recover the relations (13). However, for the most interesting case $t̂*<τ ln(1+â0S1/2)$⁠,³⁰ which occurs for very large S-values, one has to solve Eq. (14). In this case $t̂*≃τ ln(â0/â*)$⁠, and, using $â*≪â0$⁠, we obtain

where c_a ≈ 0.3. This implicit equation for $â*$ can be solved by iteration using $â*(0)≃caτ2/3S−1/3$ as the lowest order solution. Then, to the next order, we find

The logarithmic contribution to this expression is extremely important, since $γ̂*$ and $k̂*$ exhibit a super-linear dependence on $1/â*$⁠, implying that the former two are strongly affected by changes in $â*$⁠. Equation (17) reveals that $â*$ is smaller than $â*(0)$⁠. Its value depends not only on S but also on the perturbation amplitude $ŵ0=2(ψ̂0â0)1/2$ and the time scale of the driving process τ.

Substituting Eq. (17) into Eqs. (11) and (12), we duly obtain

These relations are clearly very different from relations (13). It is instructive to compare our analytical solutions with the numerical solutions of the full principle of least time computed from Eq. (4). The solutions for $k̂*$ and $γ̂*$ are shown in Figs. 1 and 2 for two different values of $ŵ0$⁠. The analytical Sweet-Parker-limit solution is accurate only for moderately high S-values, while Eqs. (18) and (19) are increasingly exact for larger values of S. The behavior of $γ̂*$ is non-monotonic in S, while $k̂*$ (proportional to the number of plasmoids) does exhibit a monotonic behavior, although it is lower when compared to the Sweet-Parker-based solution for large S.

It is important to consider the time scale of the plasmoid instability. From Eq. (17), it follows that

Naturally, this time depends on the calibration of the “clock,” i.e., on the starting time that is set by the initial inverse aspect ratio $â0$⁠. While this constitutes a ready observable, it is not the intrinsic time scale of the instability (τ_p). This is because the mode $k̂*$ remains quiescent for an extended period of time, while it is subject to rapid growth only over a small fraction of $t̂*$⁠.

The emergent mode effectively starts growing when its instantaneous growth rate is such that $γ̂(k̂*,t̂on)>1/τ$⁠. Therefore, upon using Eq. (9) and retaining the dominant terms, we obtain $âon≃k̂*S−1/2τ3/2$⁠. Consequently, the intrinsic time scale of the plasmoid instability, defined via $τp=τ ln(âon/â*)$⁠, becomes

Because of the very weak dependence on S (of the form $ln lnS$⁠) and the perturbation amplitude, it is manifest that the intrinsic time scale of the plasmoid instability is near-universal for exponentially thinning current sheets.

Hitherto, we have focused our attention on the exponential case on account of its commonality. However, other cases may occur, in which the current sheet thinning depends algebraically on time, as is possible for forced reconnection in a stable plasma.^16,31–34 Therefore, we consider a generalized current shrinking function of the form

which recovers the exponential thinning in the limit n → ∞. Here, we proceed in a more straightforward way by exploiting the fact that $γ̂*≈γ̂f(t̂*)$⁠, where $γ̂f$ is the instantaneous growth rate of the fastest growing mode. This implies that we can approximate $γ̂(â)$ with $γ̂s(â)$ in Eq. (14). In the large $t̂*$ regime, we recover the Sweet-Parker-limit solution. Otherwise, for $t̂*<nτ[(1+â0S1/2)1/n−1]$⁠, we obtain a generalized version of Eq. (16), using again $Δ̂′â≃2(k̂â)−1$ and $â*≪â0$⁠. Evaluating the implicit equation along the same lines as Eq. (16), we obtain

where $ν:=2n/(2+3n)$ and $λ:=ck−2/5n/(2+4n)$⁠. For n → ∞, we recover the same scaling as Eq. (17), while it is evident that different choices of n lead to different expressions for $â*$⁠. The n-dependence of $â*$ is illustrated in Fig. 3. We can see that $â*$ decreases for larger values of n, and that the domain of validity of the Sweet-Parker-limit solution reduces with decreasing n-values.

Equation (24) enables us to obtain the scaling relations of the plasmoid instability for different cases of current sheet thinning. Let us focus on a current sheet thinning that is inversely proportional in time (n = 1), which has been studied widely, especially in analyzing forced reconnection in Taylor's model.^{16,31–33,35} We find that

Therefore, substituting this relation into Eqs. (11) and (12), we obtain $k̂*∼L11/2/(â0τ)1/2, γ̂*∼S−1/5L13/5/(â0τ)3/5$⁠, and $ŵ*∼S−2/5L1−3/10(â0τ)3/10$⁠, where $L1:=ln(τ1/10â03/5/S3/10ŵ0)$⁠. The final time is $t̂*∼τâ0/â*$⁠, much higher than the exponential thinning case. These scaling laws are considerably different when compared to the latter – the plasmoid instability is less violent and a low number of plasmoids emerge in the nonlinear phase.

In this letter, we have generalized the previous treatments of the plasmoid instability by formulating a principle of least time for plasmoids in time evolving current sheets. We have shown that the scaling relations of all relevant parameters are dependent not only on S but also on the perturbation amplitude and the characteristic rate of current sheet thinning. We also presented a detailed explanation of the dynamics of the plasmoid instability – the system remains quiescent for a certain period of time until the thinning reaches a critical value. Once this value has been attained, the plasmoid instability occurs on a short time scale, leading to explosive growth of the plasmoids. Thus, the developed theory sheds new light on the onset problem of fast magnetic reconnection, and also highlights the role of the current sheet thinning in determining the onset time. Direct numerical simulations supporting the theory will appear in a forthcoming paper.

A final remark concerning a universal feature of our results is in order. It is common in all realms of science to seek the existence of power laws, despite the fact that they are, sometimes, intrinsically simplistic.³⁶ In contrast, we find that the scaling relations of the plasmoid instability are not true power laws – a result that has never been derived or predicted before.

It is a pleasure to acknowledge fruitful discussions with Eero Hirvijoki, Hantao Ji, Russell Kulsrud, Roscoe White, and Yao Zhou. This research was supported by the NSF Grant Nos. AGS-1338944 and AGS-1460169, and by the DOE Grant No. DE-AC02-09CH-11466.