A general theory of the onset and development of the plasmoid instability is formulated by means of a principle of least time. The scaling relations for the final aspect ratio, transition time to rapid onset, growth rate, and number of plasmoids are derived and shown to depend on the initial perturbation amplitude $(w\u03020)$, the characteristic rate of current sheet evolution (1/*τ*), and the Lundquist number (*S*). They are *not* simple power laws, and are proportional to $S\alpha \tau \beta [ln\u2009f(S,\tau ,w\u03020)]\sigma $. The detailed dynamics of the instability is also elucidated, and shown to comprise of a period of quiescence followed by sudden growth over a short time scale.

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,^{2} 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\u223cS\u22121/2$, Tajima and Shibata,^{1} as well as Loureiro *et al.*,^{18} have found that the growth rate *γ* and the wavenumber *k* of the plasmoid instability obey $\gamma \tau A\u223cS1/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,

^{19}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^{25} conjectured that current sheets break up when *γτ _{A}* ∼ 1. Later, Uzdensky and Loureiro

^{26}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 $\psi (k,t)=\psi 0\u2009exp(\u222bt0t\gamma (t\u2032)dt\u2032)$, where *ψ* represents the tearing eigenfunction and $\psi 0:=\psi (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^{2}

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*_{0} 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)\u226a\delta 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/\gamma )(\u2202\gamma /\u2202k)|\u226a1$, and consequently, $\delta in\u221dk\u22121/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*_{0}, and the time to the Alfvén time *τ _{A}* =

*L*/

*v*. We use carets to indicate the dimensionless quantities – $a\u0302=a/L,\u2009\delta \u0302in=\delta in/L,\u2009k\u0302=kL,\u2009\psi \u0302=\psi /LB0,\u2009t\u0302=t/\tau A$, and $\gamma \u0302=\gamma \tau A$. Therefore, the normalized magnetic diffusivity corresponds to the inverse of the Lundquist number, i.e., $\eta \u0302\u22121=S:=vAL/\eta $.

_{A}Although our framework is altogether general, we are interested in the case where *L* and *B*_{0} remain approximately constant, while the current sheet width decreases in time via $a\u0302(t\u0302)=a\u03020f(t\u0302)$. Here, $f(t\u0302)$ is a function that must obey $f(t\u03020)=1$ and $limt\u0302\u2192\u221ef(t\u0302)=S\u22121/2$. Indeed, $a\u0302=S\u22121/2$ is the natural lower limit to the thickness of a reconnection layer, due to the increase in the Ohmic heating when *B*_{0}/*a* increases.^{3}

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

where $w\u03020:=w\u0302(k\u0302,t\u03020)=2(\psi \u03020a\u03020)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\u0302*,\u2009\gamma \u0302(k\u0302*,t\u0302*):=\gamma \u0302*$ and $\delta \u0302in(k\u0302*,t\u0302*):=\delta \u0302in*$ as a function of the inverse aspect ratio $a\u0302(t\u0302*):=a\u0302*$.

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 $\gamma \u0302$. Here, for the sake of definiteness, we focus on resistive tearing modes. When the evolution of the current sheet width is slow, i.e., $a\u0302\u22121da\u0302/dt\u0302<\gamma \u0302$, the growth rate of the tearing mode can be computed using the instantaneous value of $a\u0302$. 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 $\Delta \u0302\u2032$.^{27} In the small-$\Delta \u0302\u2032$ regime, defined as $\Delta \u0302\u2032\delta \u0302in\u226a1$, the classic analysis by Furth *et al.*^{27} demonstrated that

where $c\Gamma =[(2\pi )\u22121\Gamma (1/4)/\Gamma (3/4)]4/5\u22480.55$. On the other hand, in the large-$\Delta \u0302\u2032$ regime, defined as $\Delta \u0302\u2032\delta \u0302in\u22731$, Coppi *et al.*^{28} showed that the growth rate becomes independent of $\Delta \u0302\u2032$, resulting in

As we are interested in the entire domain of $\Delta \u0302\u2032$, we seek an expression for $\gamma \u0302$ that (i) is a reasonable approximation of the exact growth rate,^{28} (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,^{29} in addition to being quite simple, we have verified that it fulfills all of the criteria described above. Furthermore, choosing a different approximation for $\gamma \u0302$ 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^{27} $\Delta \u0302\u2032a\u0302=2[(k\u0302a\u0302)\u22121\u2212k\u0302a\u0302]$. As we are not interested in the very slow-growing part of the mode evolution, we consider the regime $k\u0302a\u0302\u226a1$. 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,

^{27}the latter of which satisfies the equation $\u2202\gamma \u0302/\u2202k\u0302|k\u0302f*,t\u0302*=0$. This property arises because the additional contributions in Eq. (6) can be shown to be negligible if $w\u03020\u226a\delta \u0302in(k\u0302*,t\u03020)$. If the current layer has the time to evolve until the Sweet-Parker inverse aspect ratio ($a\u0302*\u2192S\u22121/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 $a\u0302*$ 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 $a\u0302*$ for a general current sheet evolution $a\u0302(t\u0302)$. It shows clearly that $a\u0302*$, and consequently the scaling relations of $t\u0302*,\u2009\gamma \u0302*,\u2009k\u0302*$ and $\delta \u0302in*$, cannot be universal, as they must depend on the specific form of the function $a\u0302(t\u0302)$.

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 $a\u0302\u221e=S\u22121/2$. Note that the plasmoid width (1) starts to grow only when $\gamma \u0302>1/\tau $. For $t\u0302*>\tau \u2009ln(1+a\u03020S1/2)$, the current sheet approaches $a\u0302*\u2243S\u22121/2$. Therefore, one can easily recover the relations (13). However, for the most interesting case $t\u0302*<\tau \u2009ln(1+a\u03020S1/2)$,^{30} which occurs for very large *S*-values, one has to solve Eq. (14). In this case $t\u0302*\u2243\tau \u2009ln(a\u03020/a\u0302*)$, and, using $a\u0302*\u226aa\u03020$, we obtain

where *c _{a}* ≈ 0.3. This implicit equation for $a\u0302*$ can be solved by iteration using $a\u0302*(0)\u2243ca\tau 2/3S\u22121/3$ as the lowest order solution. Then, to the next order, we find

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

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\u0302*$ and $\gamma \u0302*$ are shown in Figs. 1 and 2 for two different values of $w\u03020$. 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 $\gamma \u0302*$ is *non-monotonic* in *S*, while $k\u0302*$ (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 $a\u03020$. While this constitutes a ready observable, it is not the *intrinsic* time scale of the instability (*τ _{p}*). This is because the mode $k\u0302*$ remains quiescent for an extended period of time, while it is subject to rapid growth only over a small fraction of $t\u0302*$.

The emergent mode effectively starts growing when its instantaneous growth rate is such that $\gamma \u0302(k\u0302*,t\u0302on)>1/\tau $. Therefore, upon using Eq. (9) and retaining the dominant terms, we obtain $a\u0302on\u2243k\u0302*S\u22121/2\tau 3/2$. Consequently, the intrinsic time scale of the plasmoid instability, defined via $\tau p=\tau \u2009ln(a\u0302on/a\u0302*)$, becomes

Because of the very weak dependence on *S* (of the form $ln\u2009lnS$) 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 $\gamma \u0302*\u2248\gamma \u0302f(t\u0302*)$, where $\gamma \u0302f$ is the instantaneous growth rate of the fastest growing mode. This implies that we can approximate $\gamma \u0302(a\u0302)$ with $\gamma \u0302s(a\u0302)$ in Eq. (14). In the large $t\u0302*$ regime, we recover the Sweet-Parker-limit solution. Otherwise, for $t\u0302*<n\tau [(1+a\u03020S1/2)1/n\u22121]$, we obtain a generalized version of Eq. (16), using again $\Delta \u0302\u2032a\u0302\u22432(k\u0302a\u0302)\u22121$ and $a\u0302*\u226aa\u03020$. Evaluating the implicit equation along the same lines as Eq. (16), we obtain

where $\nu :=2n/(2+3n)$ and $\lambda :=ck\u22122/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 $a\u0302*$. The *n*-dependence of $a\u0302*$ is illustrated in Fig. 3. We can see that $a\u0302*$ 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\u0302*\u223cL11/2/(a\u03020\tau )1/2,\u2009\gamma \u0302*\u223cS\u22121/5L13/5/(a\u03020\tau )3/5$, and $w\u0302*\u223cS\u22122/5L1\u22123/10(a\u03020\tau )3/10$, where $L1:=ln(\tau 1/10a\u030203/5/S3/10w\u03020)$. The final time is $t\u0302*\u223c\tau a\u03020/a\u0302*$, 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.^{36} 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.

## References

This criterion can be reformulated as $S1/4>ca\u22123/2\tau \u22121ln(\tau 1/6S\u22121/3a\u030201/2w\u03020)$ by using the expression for $t\u0302*$ that is computed later in Eq. (21).