Analytic scaling relations are derived for a phenomenological model of the plasmoid instability in an evolving current sheet, including the effects of reconnection outflow. Two scenarios are considered, where the plasmoid instability can be triggered either by an injected initial perturbation or by the natural noise of the system (here referred to as the system noise). The two scenarios lead to different scaling relations because the initial noise decays when the linear growth of the plasmoid instability is not sufficiently fast to overcome the advection loss caused by the reconnection outflow, whereas the system noise represents the lowest level of fluctuations in the system. The leading order approximation for the current sheet width at disruption takes the form of a power law multiplied by a logarithmic factor, and from that, the scaling relations for the wavenumber and the linear growth rate of the dominant mode are obtained. When the effects of the outflow are neglected, the scaling relations agree, up to the leading order approximation, with previously derived scaling relations based on a principle of least time. The analytical scaling relations are validated with numerical solutions of the model.

## I. INTRODUCTION

Thin current sheets, which appear to be ubiquitous in laboratory, space, and astrophysical plasmas, are known to be unstable to the plasmoid instability,^{1–3} which disrupts current sheets to form smaller structures such as plasmoids (or flux ropes) and secondary current sheets. Plasmoid-instability-mediated disruption of reconnecting current sheets plays a crucial role in triggering the transition from slow to fast magnetic reconnection^{3–9} and also modifies the energy cascade in magnetohydrodynamic (MHD) turbulence, making the energy power spectrum steeper than that predicted by traditional MHD turbulence theories.^{10–17}

In early theoretical studies of the plasmoid instability, it was commonly assumed that the aspect ratio of the current sheet follows the Sweet-Parker scaling $ L / a \u2243 S $,^{18,19} where *S* ≡ *LV _{A}*/

*η*is the Lundquist number. Here,

*L*and

*a*denote the half-length and the half-width of the current sheet, respectively,

*V*is the Alfvén speed, and

_{A}*η*is the magnetic diffusivity. For the resistive magnetohydrodynamics (MHD) model, this assumption yields the scaling relations

*γ*∼

*S*

^{1/4}for the linear growth rate and

*k*∼

*S*

^{3/8}for the wavenumber of the fastest growing mode.

^{2,3,20}Following the same assumption, scaling relations have also been derived for other models, including Hall MHD

^{21}and viscoresistive MHD.

^{22}

However, a current sheet in nature typically forms dynamically, starting from a broader sheet that gradually thins down. The fact that the Sweet-Parker-based scaling *γ* ∼ *S*^{1/4} diverges in the asymptotic limit of *S* → *∞* indicates that for a high-*S* current sheet, the disruption may occur before the current sheet realizes the Sweet-Parker aspect ratio.^{23} The precise conditions for current sheet disruption have been a topic of active research in recent literature.^{9,23–29} By using a principle of least time, Comisso *et al.* showed that the current sheet width, the linear growth rate, and the dominant mode wavenumber at disruption do not follow power-law scaling relations.^{27,28} Huang *et al.* proposed a phenomenological model incorporating the effects of reconnection outflow.^{9} Numerical solutions of this model have been extensively tested with results from direct numerical simulations.

The objective of this work is to develop a methodology to obtain analytical scaling relations for the phenomenological model of Huang *et al.* for a special class of evolving current sheets, which is modeled by a Harris sheet thinning down exponentially in time, where the upstream magnetic field remains constant. This paper is organized as follows. Section II gives an overview of the model. Section III gives a formal solution of the model equation based on the method of characteristics. The solution is then applied in Sec. IV to obtain scaling relations for key features at the disruption time, such as the current sheet width, the linear growth rate, and the dominant mode wavenumber. Here, we consider two possible scenarios. In the first scenario, an initial noise is injected into the system to trigger the plasmoid instability. In the second scenario, the plasmoid instability evolves from the natural noise of the system, referred to as the “system noise” in this paper.^{30} The two scenarios lead to different scaling relations because the initial noise tends to decay during an early time when the linear growth of the plasmoid instability is not sufficiently fast to overcome the advection loss caused by the reconnection outflow; on the other hand, the system noise represents the lowest level of fluctuations pertaining to the system even when no external noise is injected. In Sec. V, we consider the case where the effects of the outflow are neglected. In this case, we show that the scaling relations agree with that of Comisso *et al.* up to the leading order approximation. We also validate the analytical scaling relations by comparing them with numerical solutions. We conclude in Sec. VI.

## II. PHENOMENOLOGICAL MODEL

In this section, we give a brief outline of the phenomenological model. The reader is referred to Ref. 9 for the details of the derivation. Figure 1 shows a schematic diagram for the system we consider here. The following equation governs the linear phase of the plasmoid instability in an evolving current sheet:

Here, $ f ( k , t ) \u2261 | B \u0302 z ( k , t ) | / B 0 L 0 $ is the Fourier amplitude of the magnetic field fluctuation normal to the current sheet as a function of the wavenumber *k* and the time *t*, normalized to the characteristic magnetic field *B*_{0} and length *L*_{0} of the system, *v _{o}* is the outflow speed,

*γ*(

*k*,

*t*) is the linear growth rate of the tearing instability, and

*L*is the half-length of the current sheet. Here, the linear growth rate

*γ*explicitly depends on time because the current sheet can evolve in time. The differential operator in the left-hand-side incorporates the stretching effect on wavelength due to the reconnection outflow; the right-hand-side represents the effects of linear growth, advection loss due to the outflow, and evolution of the current sheet length. The domain of wavenumber

*k*is limited to

*k*≥

*π*/

*L*because the wavelength must be smaller than the current sheet length 2

*L*.

The linear growth rate *γ*(*k*, *t*) depends on the profile of magnetic field *B _{x}*(

*z*,

*t*), which can be a function of time depending on the system under consideration. The model also allows the current sheet half-length

*L*to change in time, as represented by the (1/2

*L*)

*dL*/

*dt*term. The outflow speed

*v*

_{0}can also be time-dependent. To integrate the model equation, we need

*γ*(

*k*,

*t*),

*v*(

_{o}*t*),

*L*(

*t*), and an initial condition

*f*=

*f*

_{0}(

*k*) as inputs. With these conditions provided, we can integrate equation (1) in time until the plasmoid instability enters the nonlinear regime, which occurs when the typical half-width

*w*of plasmoids exceeds the inner layer half-width

*δ*. At that moment, the fluctuating part of the current density $ J \u0303 $ is of the same order of the background current density

*J*, implying that the current sheet has lost its integrity. Therefore, we take

*w*=

*δ*as the condition for current sheet disruption.

This condition for current sheet disruption depends critically on the half-width *w* of plasmoids. However, the precise definition of *w* is nuanced. The nuance comes from the fact that while *w* is well-defined when the magnetic fluctuation $ B \u0303 $ is a single-mode Fourier harmonic (as is often assumed in tearing mode analyses), it is not so when the fluctuation is a superposition of a continuum of Fourier modes. In this model problem, fortunately, a meaningful half-width *w* of plasmoids can be defined because, as we will see, the solution *f* typically becomes localized around a dominant wavenumber *k _{d}* as time progresses. For a spectrum

*f*(

*k*) localized around

*k*=

*k*, we define the plasmoid half-width

_{d}*w*as

where the prime denotes *d*/*dz* evaluated at *z* = 0. In Eq. (2), the fluctuation amplitude $ B \u0303 $ must take into account the contribution from all the neighboring modes of the dominant mode. Precisely, a superposition of all the modes in the range *k* ∈ [*k _{d}*/

*ξ*,

*k*] gives a fluctuation amplitude

_{d}ξHere, the *O*(1) parameter *ξ* sets the range of superposition. As the Fourier spectrum *f*(*k*) is usually well-localized at the disruption time, the amplitude $ B \u0303 $ is insensitive to the choice of *ξ*, provided that *ξ* is not too close to unity (e.g., *ξ* = 1.5 was used in Ref. 9). Note that for a single-mode perturbation, definition (2) reduces to the usual expression for the magnetic island half-width.^{31}

The plasmoid half-width *w* must be compared with the inner layer half-width *δ* to determine the condition for current sheet disruption. As the inner layer half-width *δ* depends on the wavenumber, we adopt the inner layer half-width at the dominant wavenumber *k _{d}* for this comparison. For resistive tearing modes we consider in this study, the inner layer half-width for

*k*=

*k*is given by

_{d}where *η* is the magnetic diffusivity. Here, the prime again denotes *d*/*dz* evaluated at *z* = 0, and *V _{A}* is the Alfvén speed of the reconnecting component

*B*. The criterion

_{x}*w*=

*δ*, with

*w*and

*δ*defined by Eqs. (2)–(4), then determines the condition for current sheet disruption. This criterion has been extensively tested and validated with numerical simulations.

^{9}

## III. INTEGRAL SOLUTION OF THE MODEL EQUATION FOR A THINNING HARRIS SHEET

Until now, the model has been general regarding the time evolution of both the current sheet profile and the reconnection outflow. To fix ideas and make the model analytically tractable, we assume that the current sheet can be modeled by a Harris sheet and the time evolution of the width *a* obeys the form

In this form, assuming *a _{∞}* ≪

*a*

_{0}, the current sheet width

*a*decreases exponentially at early times with $ a \u2243 a 0 e \u2212 t / \tau $, and

*a*→

*a*as

_{∞}*t*→

*∞*. We further assume that the upstream magnetic field

*B*and the current sheet half-length

_{x}*L*remain fixed in time and that the outflow speed

*v*is the upstream Alfvén speed

_{o}*V*. Under these conditions, if the plasmoid instability does not disrupt the current sheet, eventually the current sheet will approach the Sweet-Parker current sheet, which is a steady-state solution of reconnection governed by resistive MHD. Therefore, we set the asymptotic half-width

_{A}*a*to be the Sweet-Parker width, i.e., $ a \u221e = L / S $, where

_{∞}*S*≡

*LV*/

_{A}*η*is the Lundquist number. These assumptions are consistent with direct numerical simulations reported in Ref. 9. Heuristic derivations for Eq. (5) can be found in Refs. 32 and 9.

Under the assumption that the upstream *B _{x}* and the current sheet length

*L*are fixed in time, it is natural to choose

*B*

_{0}=

*B*and

_{x}*L*

_{0}=

*L*as the characteristic magnetic field and length scale. Hereafter, we normalize lengths and times with respect to the current sheet half-length

*L*and the Alfvén time

*τ*≡

_{A}*L*/

*V*, as follows: $ t \u0302 = t / \tau A , \u2009 \tau \u0302 = \tau / \tau A , k \u0302 = kL , \u2009 a \u0302 = a / L $, and $ \gamma \u0302 = \gamma \tau A $. The model equation in normalized variables is

_{A}where *dL*/*dt* = 0 is assumed.

We can solve Eq. (6) by the method of characteristics. The characteristic starting from $ k \u0302 = k \u0302 0 $ at $ t \u0302 = t \u0302 0 $ is given by

which represents the mode-stretching effect due to the reconnection outflow. Now the differential operator on the left-hand-side of Eq. (6) is simply the time derivative along characteristics and hence the solution can be obtained by integrating along them, yielding

To perform the integration, we need the linear growth rate $ \gamma \u0302 $ as a function of the wavenumber and time. The linear growth rate of the tearing mode is governed by the tearing stability index $ \Delta \u2032 $,^{33} which depends entirely on the solution of linearized ideal MHD force-free equation in the outer region away from the resonant surface (where **k**⋅**B **=** **0). The tearing instability requires the condition $ \Delta \u2032 > 0 $ being satisfied. While the complete dispersion relation for resistive tearing modes can be expressed in terms of $ \Delta \u2032 $, the expression is involved and requires solving a transcendental equation.^{34} Analytically tractable approximations can be obtained in two limits: The small- $ \Delta \u2032 $ regime for short-wavelength modes and the large- $ \Delta \u2032 $ regime for long-wavelength modes. In the small- $ \Delta \u2032 $ regime (large *k* and $ \Delta \u2032 \delta \u226a 1 $), the linear growth rate is given by^{33}

where $ C \Gamma = ( \pi \u2212 1 \Gamma ( 1 / 4 ) / \Gamma ( 3 / 4 ) ) 4 / 5 \u2248 0.953 $. In the large- $ \Delta \u2032 $ regime (small *k* and $ \Delta \u2032 \delta \u226b 1 $), the linear growth rate is approximately given by^{34}

Here, the derivative $ ( k V A ) \u2032 $ is evaluated at the resonant surface. For a Harris sheet profile with $ V A ( z ) = V A tanh ( z / a ) $, the tearing stability index $ \Delta \u2032 $ is given by

From Eqs. (9)–(11), we obtain the normalized growth rates for a Harris sheet in the small- $ \Delta \u2032 $ regime

and in the large- $ \Delta \u2032 $ regime

for an arbitrary ζ. In this work, we take the value ζ = 3/2, which gives a nearly exact approximation of the true dispersion relation. Furthermore, we will neglect the factor $ ( 1 \u2212 k \u0302 2 a \u0302 2 ) $ in Eq. (12) for our analytical derivations. This approximation can be justified *a posteriori* by noting that the dominant mode wavenumber $ k \u0302 d $ follows the same scaling as the fastest growing mode wavenumber $ k \u0302 max $ [see Eq. (35) and the discussion thereafter] and that $ k \u0302 max a \u0302 \u2243 ( S a \u0302 ) \u2212 1 / 4 $ (see below); therefore, in the high-*S* regime, the condition $ k \u0302 a \u0302 \u226a 1 $ is satisfied in the neighborhood (in the *k*-space) of the dominant mode, which is our primary interest. Note that, from the scaling of the disruption current sheet width, Eq. (47), the condition $ S a \u0302 \u226b 1 $ is satisfied. Putting these approximations together yields an approximate expression for the linear growth rate that is valid for $ k \u0302 a \u0302 \u226a 1 $, i.e.,

The fastest growing mode occurs at the transition between the small- $ \Delta \u2032 $ and the large- $ \Delta \u2032 $ regimes, i.e., at $ k \u0302 a \u0302 \u2243 ( S a \u0302 ) \u2212 1 / 4 $ where $ \gamma \u0302 s \u2243 \gamma \u0302 l $. More precisely, the fastest growing wavenumber $ k \u0302 max $ is given by (see, e.g., Ref. 35)

and the corresponding growth rate is

## IV. SCALING RELATIONS FOR CURRENT SHEET DISRUPTION

In this section, we derive scaling relations for key features at current sheet disruption in the high-*S* regime, under the assumptions detailed in Sec. III. Specifically, the disruption time $ t \u0302 d $, the current sheet width $ a \u0302 d $, the dominant mode wavenumber $ k \u0302 d $, and the linear growth rate $ \gamma \u0302 d $ of the dominant mode are derived analytically. Here, we consider two cases. In the first case, an initial noise is injected into the system. The initial noise is assumed to be much larger than the system noise; therefore, the latter is ignored and set to zero in the analysis. In the second case, there is no initial noise, and the system noise is the seed from which the plasmoid instability develops.

### A. Case I: Initial perturbation

Let the initial perturbation be $ f 0 ( k \u0302 0 ) $ at $ t \u0302 0 = 0 $. To integrate the linear growth rate along a characteristic in Eq. (8), it is convenient to define a new variable

In the high-*S* regime we consider here, it can be shown *a posteriori* that current sheet disruption occurs when $ a \u0302 d \u226b a \u0302 \u221e $, because $ a \u0302 d $ follows a weaker scaling with respect to *S* [see Eq. (47)] than $ a \u0302 \u221e $, which is assumed to be the Sweet-Parker width $ a \u0302 SP = S \u2212 1 / 2 $. Hence, we will ignore *a _{∞}* in Eq. (5) and assume $ a \u0302 = a \u0302 0 e \u2212 t \u0302 / \tau \u0302 $ in the following analysis. The combined effect of mode stretching ( $ k \u0302 = k \u0302 0 e \u2212 ( t \u0302 \u2212 t \u0302 0 ) $) and current sheet thinning $ ( a \u0302 = a \u0302 0 e \u2212 t \u0302 / \tau \u0302 ) $ implies that the variable

*φ*decreases exponentially in time as

where the time scale $ \tau \u0302 * $ is defined as

The linear growth rate, Eq. (15), can be expressed in terms of *φ* instead of $ k \u0302 $ as

We can now integrate the growth rate along the characteristic with a change of variable from $ t \u0302 $ to *φ*, yielding

where

Here, _{2}*F*_{1} is the hypergeometric function^{36} and

To proceed from Eq. (22), we further make the approximation that in the high-*S* regime, the contribution from the lower boundary, $ G ( \phi 0 ; \tau \u0302 ) $, is negligible compared to the contribution from the upper boundary, $ G ( \phi ; \tau \u0302 ) $, at the disruption time. This approximation can be justified *a posteriori* using the scaling relation of the disruption current sheet width $ a \u0302 d $ in Eq. (47). Because $ a \u0302 d $ decreases as *S* increases, it follows from the relation $ \phi / \phi 0 = ( a \u0302 / a \u0302 0 ) \tau \u0302 / \tau \u0302 * $ that in the high-*S* regime, the condition *φ* ≪ *φ*_{0} is satisfied at the disruption time, where the variable *φ* is typically an *O*(1) quantity in the neighborhood of the dominant mode [see the discussion in the paragraph below Eq. (34)]. The function $ G ( \phi \u2032 ; \tau \u0302 ) $ is monotonically decreasing with respect to $ \phi \u2032 $. Its overall behavior is well captured by two asymptotic expressions for small and large $ \phi \u2032 $, which can be obtained using the asymptotic expansion of the hypergeometric function near *x* → *∞*

and the Taylor expansion near *x* = 0

for the special case with *b* =* c*, Eq. (25) is not valid and must be replaced by

It follows that the function $ G ( \phi \u2032 ; \tau \u0302 ) $ is approximately

when $ \phi \u2032 \u2009 \u2272 \u2009 O ( 1 ) $, and

when $ \phi \u2032 \u2273 O ( 1 ) $. For the special case $ \tau \u0302 = 1 $ that yields $ \tau \u0302 * / \tau \u0302 \u266f = 2 / 3 $, Eq. (28) is replaced by

Because *φ*_{0} ≫ *φ* ∼ *O*(1), it follows from Eqs. (28)–(30) that $ G ( \phi 0 ; \tau \u0302 ) \u226a G ( \phi ; \tau \u0302 ) $.

Ignoring $ G ( \phi 0 ; \tau \u0302 ) $, we can simplify Eq. (22) as

where

and

The solution of the model equation can now be written as

with $ k \u0302 0 = k \u0302 e t \u0302 , \u2009 h = \tau \u0302 \u266f C \Gamma 5 / 8 S \u2212 1 / 2 a \u0302 \u2212 3 / 2 / 2 $, and $ \phi = C \Gamma \u2212 3 / 4 S 1 / 5 k \u0302 4 / 5 a \u0302 $ being substituted into the right-hand side.

The dominant wavenumber $ k \u0302 d $ at the disruption time can be obtained by locating the maximum of *f*. Because the exponential factor in Eq. (34) is highly localized, whereas most random fluctuations have a broadband spectrum, we further assume that $ f 0 ( k \u0302 0 ) $ is slowly varying compared to the exponential factor. Under this assumption, the dominant wavenumber is approximately the one that maximizes the exponent in Eq. (34).^{37} Let *φ* = *φ _{d}* correspond to the maximum of $ g ( \phi ; \tau \u0302 ) $. The value of

*φ*depends on $ \tau \u0302 $ and, in general, has to be obtained numerically; typically,

_{d}*φ*is an

_{d}*O*(1) quantity. With the value of

*φ*obtained, the dominant wavenumber $ k \u0302 d $ and linear growth rate $ \gamma \u0302 d $ are given by

_{d}and

where the constants *c _{k}* and

*c*are defined as

_{γ}and

Equation (35) shows that the dominant wavenumber and the fastest growing wavenumber follow the same scaling [cf. Eq. (16)]. However, as we will see in Table I, typically the constant *c _{k}* is smaller than the coefficient in Eq. (16). Therefore, the dominant wavenumber is smaller than the fastest-growing wavenumber. This discrepancy between the dominant wavenumber and the fastest-growing wavenumber is due to the fact that the mode amplitude at a given wavenumber $ k \u0302 $ is determined by the entire history of the linear growth rate following the characteristic. Even though the dominant mode is not the fastest growing mode at the moment, it has higher linear growth rates at earlier times.

^{9}

. | $ \tau \u0302 = log \u2009 2 , $ . | $ \tau \u0302 \u2192 0 , $ . | $ \tau \u0302 = log \u2009 2 , $ . | $ \tau \u0302 = log \u2009 2 , $ . |
---|---|---|---|---|

. | χ = 0
. | χ = 0
. | χ = 1
. | χ = 3
. |

c _{a} | 0.873 | 0.838 | 0.873 | 0.873 |

c _{k} | 0.563 | 0.705 | 0.563 | 0.563 |

c _{γ} | 0.538 | 0.576 | 0.538 | 0.538 |

c _{χ} | 0.0349 | 0.0313 | 0.0111 | 0.00111 |

. | $ \tau \u0302 = log \u2009 2 , $ . | $ \tau \u0302 \u2192 0 , $ . | $ \tau \u0302 = log \u2009 2 , $ . | $ \tau \u0302 = log \u2009 2 , $ . |
---|---|---|---|---|

. | χ = 0
. | χ = 0
. | χ = 1
. | χ = 3
. |

c _{a} | 0.873 | 0.838 | 0.873 | 0.873 |

c _{k} | 0.563 | 0.705 | 0.563 | 0.563 |

c _{γ} | 0.538 | 0.576 | 0.538 | 0.538 |

c _{χ} | 0.0349 | 0.0313 | 0.0111 | 0.00111 |

Equation (35) gives the relation between $ k \u0302 d $ and $ a \u0302 d $ but we still need to determine the latter. This final step is accomplished by using the disruption condition $ w \u0302 = \delta \u0302 $. From Eqs. (2) and (4), the disruption condition is equivalent to

where the magnetic fluctuation $ B \u0303 $ is calculated by integrating the fluctuation spectrum *f* over a neighborhood of the dominant wavenumber as prescribed in Eq. (3). Because the fluctuation spectrum *f* is localized around the dominant wavenumber, we can first expand the exponent in Eq. (34) near $ k \u0302 d $ and then approximate the integration by a Gaussian integral. Expanding *g* in the neighborhood of $ k \u0302 d $ yields

Hence, in the neighborhood of the dominant mode,

where

and $ t \u0302 d $ is the disruption time. Now we can evaluate the magnetic fluctuation $ B \u0303 $ by using Eq. (41) and extending the domain of integration to (−*∞*, *∞*) (see Fig. 2 for an illustration of this approximation), yielding

Using Eqs. (43) and (36) in the disruption condition (39) yields the condition to determine the disruption width $ a \u0302 d $

In deriving Eq. (44), we have used Eq. (35) for $ k \u0302 d $, Eq. (32) for *h*, Eq. (42) for *A*, and the relation $ e t \u0302 d = ( a \u0302 0 / a \u0302 d ) \tau \u0302 $.

Equation (44) is the fundamental equation to determine the disruption current sheet width $ a \u0302 d $ for a general initial condition $ f 0 ( k \u0302 0 ) $; its validity only requires that $ f 0 ( k \u0302 0 ) $ is slowly varying such that the dominant wavenumber in the spectrum, given by Eq. (34), is determined by the peak in the exponent. However, Eq. (44) is difficult to solve for a general $ f 0 ( k \u0302 0 ) $ without resorting to numerical methods. To make further progress, we assume that the initial condition takes a power-law form $ f 0 ( k \u0302 0 ) = \u03f5 k \u0302 0 \u2212 \chi $. Then, Eq. (44) becomes

where

Equation (45) can be solved for $ a \u0302 d $ in terms of the lower branch *W*_{–1} of the Lambert *W* function.^{38} The solution is

where

and

An *a posteriori* consistency check for the Gaussian integral approximation in Eq. (43) is in order. Note that *g*(*φ _{d}*) and $ g \u2033 ( \phi d ) $ are typically

*O*(1) quantities; therefore, the validity of the Gaussian integral approximation requires the condition $ h ( a \u0302 d ) \u226b 1 $ being satisfied. Using Eq. (47) in Eq. (32), the condition is

This condition is satisfied only if $ | \Xi | \u226a 1 $. Therefore, the noise amplitude *ϵ* must be sufficiently small for the approximation to be valid.

Because the condition $ | \Xi | \u226a 1 $ is satisfied, we can use the asymptotic expansion of *W*_{–1}(*z*) as *z* → 0^{−}, i.e.,

to obtain the leading order approximation for $ a \u0302 d $

Here, we have ignored an *O*(1) factor $ c \chi ( \theta c a 3 / 2 ) \u2212 1 / \theta $ within the logarithm, which can be attributed to the next order correction of the form $ log \u2009 ( log \u2009 ( \u2026 ) ) $.

Now we have obtained the current sheet width $ a \u0302 d $ when disruption occurs. The disruption time $ t \u0302 d $ can be obtained by $ t \u0302 d = \tau \u0302 \u2009 log \u2009 ( a \u0302 0 / a \u0302 d ) $. The dominant wavenumber $ k \u0302 d $ and the linear growth rate $ \gamma \u0302 d $ can be obtained by substituting $ a \u0302 d $ into Eqs. (35) and (36).

### B. Case II: System noise

The calculation in Sec. IV A assumes that an initial perturbation is applied at $ t \u0302 = 0 $. If the current sheet width is sufficiently broad such that the linear growth rate $ \gamma \u0302 ( k \u0302 ) < 1 / 2 $, and then the perturbation amplitude at wavenumber $ k \u0302 $ decreases in time because of the advection loss caused by the reconnection outflow [Eq. (6)]. If the linear growth rate never rises above 1/2, then the perturbation amplitude will decrease monotonically in time and asymptotically approach zero as $ t \u0302 \u2192 \u221e $. However, since noise is present in any natural system, there will always be some fluctuations in the system. Noise can be introduced into the model by explicitly adding a source term or by setting a lower bound (as a function of $ k \u0302 $) to the fluctuation amplitude. Here, we adopt the second approach: Instead of an initial perturbation, we now use $ f 0 ( k \u0302 ) $ to describe the system noise, i.e., $ f 0 ( k \u0302 ) $ represents the lower bound of $ f ( k \u0302 $). In this section, we consider the situation where the system noise provides the seed for the plasmoid instability.

Because $ \gamma \u0302 $ must be greater than 1/2 for a mode to grow when we integrate along a characteristic in Eq. (8), the starting time $ t \u0302 0 $ is determined by the condition $ \gamma \u0302 = 1 / 2 $. For neighboring modes of the dominant mode at the disruption time, it can be shown *a posteriori* that those modes start to grow when they are in the small- $ \Delta \u2032 $ regime.^{39} Using the relations $ k \u0302 0 = k \u0302 e t \u0302 \u2212 t \u0302 0 $ and $ a \u0302 ( t \u0302 0 ) = a \u0302 e ( t \u0302 \u2212 t \u0302 0 ) / \tau \u0302 $ in the condition using the small- $ \Delta \u2032 $ dispersion relation

yields

hence, the wavenumber $ k \u0302 0 $ when the mode starts to grow is

Following the same calculation leading to Eq. (34), the solution is

with Eqs. (55) and (56), $ h = \tau \u0302 \u266f C \Gamma 5 / 8 S \u2212 1 / 2 a \u0302 \u2212 3 / 2 / 2 $, and $ \phi = C \Gamma \u2212 3 / 4 S 1 / 5 k \u0302 4 / 5 a \u0302 $ being substituted into the right-hand-side.

The disruption condition (39) yields the equation to determine $ a \u0302 d $

If we assume a power-law system noise $ f 0 ( k \u0302 ) = \u03f5 k \u0302 \u2212 \chi $, Eq. (59) becomes

where

The solution for $ a \u0302 d $ can be expressed in terms of the lower-branch Lambert function *W*_{−1} as

where

and *c _{a}* is defined in Eq. (48). Using the asymptotic expansion (52), the leading order approximation of $ a \u0302 d $ is

Here, we again ignore an *O*(1) factor within the logarithm, which can be attributed to the next order correction of the form log(log(…)).

## V. DISCUSSION

### A. Effects of the reconnection outflow

Our phenomenological model includes the effects of the reconnection outflow. It is instructive to investigate how the outflow affects the scaling relations. The effects of the outflow can be turned off by setting *v _{o}* to zero in Eq. (1). While we have repeated the calculation for cases without the outflow, here we present an easier way that yields the same results. Neglecting the effects of the outflow can be formally achieved by assuming that the current sheet thins down almost instantaneous, i.e., the current sheet thinning time scale

*τ*≪

*τ*, such that the mode-stretching and advective loss due to the outflow have no time to take effect. Taking the limit $ \tau \u0302 = \tau / \tau A \u2192 0 $ is thus equivalent to neglecting the effects of the outflow. In this limit, the scaling relations obtained in Secs. IV A and IV B reduce to the same results; i.e., it makes no difference whether the plasmoid instability evolves from an initial noise or a system noise if the effects of the outflow are ignored.

_{A}In the $ \tau \u0302 \u2192 0 $ limit, the two time scales $ \tau \u0302 * $ and $ \tau \u0302 \u266f $ become identical to $ \tau \u0302 $. The variables *θ* and Ξ become

and

The disruption current sheet width is given by

with the leading order approximation

Therefore, ignoring the outflow changes the time scale from $ \tau \u0302 \u266f $ to $ \tau \u0302 $ in the scaling relation of $ a \u0302 d $; the power indices of $ \tau \u0302 $ and *S* within the logarithmic factor are also different. Otherwise, the scaling relation (69) is similar to Eqs. (53) and (65).

### B. Comparison with the scaling relations of Comisso *et al.*

Now we are in a position to compare the analytical scaling relations in this work with that obtained previously by Comisso *et al.*,^{27,28} because the latter also ignore the effects of outflow. However, as Comisso *et al.* used a different way to describe the fluctuations, we need a translation of notations before a comparison can be made.

Comisso *et al.* described the fluctuations by the island size $ w \u0302 = 2 ( \psi \u0302 a \u0302 ) 1 / 2 $ as a function of $ k \u0302 $ and $ t \u0302 $, where $ \psi \u0302 $ is the magnetic flux of the island. They propose a “principle of least time” that determines the dominant mode to disrupt the current sheet by the wavenumber that satisfies the condition $ w \u0302 ( k \u0302 , t \u0302 ) = \delta \u0302 ( k \u0302 , t \u0302 ) $ with the shortest time. For an initial condition $ \psi \u0302 0 ( k \u0302 ) = \u03f5 k \u0302 \u2212 \alpha $, the leading order approximation of the current sheet width at disruption is given by Eq. (26) of Ref. 28

The present work describes the fluctuations by a spectrum $ f ( k \u0302 , t \u0302 ) $ of the magnetic field. To compare Eq. (70) with Eq. (69), we need to find a correspondence between *f* and $ \psi \u0302 $. This can be obtained by noting that the magnetic field fluctuation $ B \u0303 \u223c k \u0302 1 / 2 f B 0 $ [from Eq. (3)] and $ B \u0303 \u223c k \u0302 \psi \u0302 B 0 $ [from Eq. (2) and the relation $ w \u0302 = 2 ( \psi \u0302 a \u0302 ) 1 / 2 $]; hence, the correspondence is $ f \u223c k \u0302 1 / 2 \psi \u0302 $. For the assumed power-law initial perturbations $ f = \u03f5 k \u0302 \u2212 \chi $ and $ \psi \u0302 = \u03f5 k \u0302 \u2212 \alpha $, we have the correspondence $ \alpha \u2194 \chi + 1 / 2 $. Substituting *α* → *χ* + 1/2 into Eq. (69), we recover the scaling relation (70). Therefore, we conclude that when the effects of the outflow are ignored, the phenomenological model of this work yields the same scaling relations as the relations obtained using the principle of least time. However, the agreement is only up to the leading order approximation. It can be shown that the two approaches give different results if the next order corrections of the Lambert *W* function of the form log(log(…)) are included.

### C. Comparison with numerical solutions of the model equation

We have employed several approximations when deriving the analytical scaling relations. These approximations include: (i) ignoring the asymptotic current sheet width *a _{∞}* in Eq. (5); (ii) using an approximate linear growth rate Eq. (15) that ignores the $ k \u0302 a \u0302 = 1 $ stability threshold; (iii) ignoring the contribution from the lower bound

*φ*

_{0}of the integral in Eq. (22); (iv) assuming the initial or system noise spectrum $ f 0 ( k \u0302 ) $ to be slowly varying compared to the exponential factor in Eq. (34); and (v) replacing the integral in (43) by a Gaussian integral and pushing the bounds to ±

*∞*. Among these approximations, (iv) is an assumption that must be satisfied by the noise spectrum $ f ( k \u0302 0 ) $; the other approximations can be justified

*a posteriori*in the limit of high-

*S*and low-noise.

We now test the analytical scaling relations with the results from numerical solutions of the model equation. For the numerical solutions, we do not employ any approximation except the linear growth rate, where Eq. (14) with *ζ* = 3/2 is used. Using the approximate linear growth rate is not a significant compromise because the expression is nearly exact; note that the analytical calculation employs a further simplified expression for the linear growth rate, Eq. (15).

To fix ideas, we employ the same setup as in Ref. 9, i.e., with $ a \u0302 0 = 1 / \pi $ and $ \tau \u0302 = log \u2009 2 $. We vary the parameters *ϵ*, *χ*, and *S* for cases with either an initial noise or a system noise, with or without the outflow. We integrate the normalized model equation (6) forward in time along characteristics. Because the wavenumber $ k \u0302 $ decreases exponentially in time along a characteristic, it is convenient to set the mesh in the wavenumber $ k \u0302 $ equispaced with respect to $ log \u2009 k \u0302 $ and set the time step $ \Delta t \u0302 = \Delta \u2009 log \u2009 k \u0302 $. In this way, the wavenumber $ k \u0302 $ shifts exactly one grid space in one time step, which significantly simplifies the time integration along characteristics. We perform the time integration with a second-order trapezoidal scheme. For cases with system noise, in each time step the solution $ f ( k \u0302 ) $ is compared with $ f 0 ( k \u0302 ) $; we set $ f ( k \u0302 ) $ to $ f 0 ( k \u0302 ) $ for those $ k \u0302 $ where $ f ( k \u0302 ) < f 0 ( k \u0302 ) $. In each time step, we identify the dominant wavenumber and calculate the plasmoid size $ w \u0302 $ with Eq. (3), where *ξ* = 1.5 has been used. The time integration continues until the disruption condition $ \delta \u0302 = w \u0302 $ is met. Then we record $ a \u0302 d , \u2009 k \u0302 d $, and $ \gamma \u0302 d $ at the disruption time $ t \u0302 d $.

The analytical scaling relations involve constants *c _{a}*,

*c*,

_{k}*c*, and

_{γ}*c*. These constants depend on $ \tau \u0302 $ because they depend on

_{χ}*φ*that maximizes the function $ g ( \phi ; \tau \u0302 ) $; additionally,

_{d}*c*and

_{a}*c*also depend on

_{χ}*g*(

*φ*) and $ g \u2033 ( \phi d ) $. For $ \tau \u0302 = log \u2009 2 $, we numerically obtain

_{d}*φ*≈ 0.655,

_{d}*g*(

*φ*) ≈ 0.841, and $ g \u2033 ( \phi d ) \u2248 \u2212 0.820 $. If the outflow is neglected, i.e., taking the limit $ \tau \u0302 \u2192 0 $, they become

_{d}*φ*≈ 0.784,

_{d}*g*(

*φ*) ≈ 0.791, and $ g \u2033 ( \phi d ) \u2248 \u2212 0.629 $. These values are used in Eqs. (37), (38), (46), and (48) to obtain the constants

_{d}*c*,

_{a}*c*,

_{k}*c*, and

_{γ}*c*. The values of these constants for different settings reported in this work are listed in Table I. In the following comparisons, we use the precise analytical scaling relations involving the Lambert function, but the leading order approximate scaling relations are almost as good.

_{χ}We start with a setup having the initial noise $ f 0 ( k \u0302 ) = \u03f5 $. The scalings for $ t \u0302 d , \u2009 a \u0302 d , \u2009 k \u0302 d $, and $ \gamma \u0302 d $ with respect to *S* are shown in Fig. 3 for two cases *ϵ* = 10^{−10} and *ϵ* = 10^{−20}. Here, the symbols denote the scalings obtained from numerical solutions, the dash-dotted lines are the analytical scalings, and the black dashed lines are scalings based on the fastest growing mode of the Sweet-Parker current sheet [i.e., by setting $ a \u0302 = a \u0302 SP = S \u2212 1 / 2 $ in Eqs. (16) and (17)]. We can see that the analytical scalings accurately agree with the numerically results in the high-*S* regime. On the other hand, while the Sweet-Parker scaling $ a \u0302 SP = S \u2212 1 / 2 $ agrees with $ a \u0302 d $ in the low to moderately-high *S* regimes, the Sweet-Parker-based scalings for $ k \u0302 max $ and $ \gamma \u0302 max $ considerably depart from the numerical results of $ k \u0302 d $ and $ \gamma \u0302 d $; the difference is especially substantial when *S* is low. The discrepancies are due to the stretching effect of outflow that makes the wavelength of the dominant mode significantly longer than that of the fastest growing mode.

Next, we test the analytical scalings for the other two cases: (1) system noise instead of initial noise and (2) ignoring the effects of outflow. Here, we set $ f 0 ( k \u0302 ) = \u03f5 $ with *ϵ* = 10^{−20}. We only present the results for $ k \u0302 d $ and $ \gamma \u0302 d $ since they provide more stringent tests than $ t \u0302 d $ and $ a \u0302 d $. The scalings of $ k \u0302 d $ and $ \gamma \u0302 d $ with respect to *S* are shown in Fig. 4; for comparison, we also show results that have been presented in Fig. 3 for cases of initial noise. We again see that the analytical scalings are in excellent agreement with the numerical results in the high-*S* regime. In that regime, the scalings of $ k \u0302 d $ are quite similar for all three cases; on the other hand, ignoring the effect of the outflow only slightly affects the scaling of $ \gamma \u0302 s $ compared to the case with a system noise, whereas $ \gamma \u0302 s $ is higher for the case with an initial noise. The linear growth rate $ \gamma \u0302 s $ is higher for the initial noise case because the noise decays during the early period when all the modes are stable; therefore, the disruption occurs at a later time compared to the other two cases. As the current sheet is thinner at a later time, the linear growth rate is higher. Note that the Sweet-Parker-based scalings agree quite well with numerical results at the low-*S* regime when the outflow is ignored, confirming our earlier statement that the discrepancies are due to the effect of the outflow.

The analytical scaling relations are derived under the assumption that the noise spectrum $ f 0 ( k \u0302 ) $ is slowly varying compared to the exponential factor in Eq. (34). For the cases we have tested so far, this assumption is trivially satisfied with $ f 0 ( k \u0302 ) = \u03f5 $. Now we test the analytical scalings for a power-law system noise $ f 0 ( k \u0302 ) = \u03f5 k \u0302 \u2212 \chi $, with *ϵ* = 10^{−20}. Figure 5 shows the scalings of $ k \u0302 d $ and $ \gamma \u0302 d $ for the cases *χ* = 1 and *χ* = 3. We can see that even for a relatively steep noise spectrum with *χ* = 3, the analytical scalings remain quite accurate.

## VI. CONCLUSION

In conclusion, we have derived analytical scaling relations for the phenomenological model of plasmoid-mediated disruption of evolving current sheets, including the effects of the reconnection outflow. The analytical scaling relations are derived for two scenarios, one with an initial noise and another with a system noise. The former is useful for comparison with numerical simulations or laboratory experiments where an initial perturbation can be injected; the latter is more suitable for describing the plasmoid instability in natural systems. If the effects of outflow are neglected, and with a proper translation of notations, the scalings obtained from the present model agree with that from previous works based on a principle of least time,^{27,28} up to the leading order approximation.

The phenomenological model is valid in the linear regime of the plasmoid instability up to the point of current sheet disruption. The time scale *t _{d}* for current sheet disruption is typically on the order of a few Alfvén time

*τ*. Here, we estimate the Alfvén times for various systems of interest. For a large-scale coronal mass ejection (CME), the post-CME current sheet can extend to a length beyond several times that of the solar radius $ R \u2299 \u2243 7 \xd7 10 8 m $; the Alfvén speed can be estimated from the speeds of moving blobs in the current sheet, typically in the range of 2 × 10

_{A}^{5}–10

^{6}m/s.

^{40,41}Assuming

*L*∼ 10

^{9}m and

*V*∼ 10

_{A}^{6}m/s, the Alfvén time

*τ*is on the order of 1000 s. For ultraviolet bursts in the solar transition region,

_{A}^{42,43}we estimate the current sheet length

*L*∼ 10

^{6}m from forward modeling

^{44}and the Alfvén speed

*V*∼ 10

_{A}^{5}m/s from the Doppler broadening of emission line profiles, yielding

*τ*∼ 10 s. For laboratory magnetic reconnection experiments,

_{A}^{45,46}the Alfvén time

*τ*is typically on the order of microseconds, which is within the temporal resolution of

_{A}*in situ*measurements.

^{47}Thus, it may be possible to test some of the predictions from the phenomenological model with future generations of laboratory experiments.

^{48}

Although the present study focuses on resistive MHD, the methodology is quite general. The analysis could be generalized to include other effects such as viscosity, the Hall effect, or ambipolar diffusion from partially ionized plasmas. It is also important to further test the predictions of the phenomenological model with full resistive MHD simulations, especially in the high-*S* regime. While the present analytical calculation is limited to the special class of evolving current sheets where the upstream magnetic field remains constant, the phenomenological model could also describe cases where the upstream magnetic field evolves in time, e.g., as in Refs. 49 and 50. These directions will be pursued in the future.

## ACKNOWLEDGMENTS

We thank the anonymous referee for many constructive suggestions. This work was supported by the National Science Foundation, Grant Nos. AGS-1338944 and AGS-1460169; the Department of Energy, Grant No. DE-SC0016470; and the National Aeronautics and Space Administration, Grant No. 80NSSC18K1285. A. Bhattacharjee would like to acknowledge the hospitality of the Flatiron Institute, supported by the Simons Foundation.

## References

It always requires some “seed” noise to trigger the plasmoid instability. In MHD simulations that tend to be clean (i.e., with very low system noise due to discretization and round-off errors), a seed noise is usually provided by adding a random noise in the initial condition; this is an example of initial noise. In contrast, particle-in-cell simulations typically do not require an initial noise to trigger the plasmoid instability because they are inherently noisy, but an initial noise can be seeded as well. The random noise in particle-in-cell simulations is an example of system noise. The same applies to laboratory and natural systems. The plasmoid instability in such systems is usually triggered by the system noise, but an initial noise may also be introduced in laboratory experiments.

This assumption is not valid if $f0(k\u03020)$ is narrowband. In that case, we have to locate the maximum of the full $f(k\u0302,t\u0302)$.

This can be seen as follows. Note that the dominant wavenumber $k\u0302d$ scales in the same way as the fastest-growing wavenumber $k\u0302max\u223cS\u22121/4a\u0302d\u22125/4$ at the disruption time, where the latter is at the transition between the small-$\Delta \u2032$ and the large-$\Delta \u2032$ regimes. Because $k\u0302=k\u0302det\u0302d\u2212t\u0302$ along a characteristic and $a\u0302\u221de\u2212t\u0302/\tau \u0302$, the wavenumber $k\u0302$ is larger than $k\u0302max\u223cS\u22121/4a\u0302de\u22125(t\u0302d\u2212t\u0302)/4\tau \u0302$ at an earlier time $t\u0302$ with $t\u0302d\u2212t\u0302>O(1)$, which places the mode in the small-$\Delta \u2032$ regime.