Scalings pertaining to current sheet disruption mediated by the plasmoid instability

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 ≃ 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_A is the Alfvén speed, and η 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²¹ and viscoresistive MHD.²² 

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.²³ 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.⁹ 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.³⁰ 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.

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 ) ≡ | B ̂ 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₀ and length L₀ 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 2L.

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/2L)dL/dt term. The outflow speed v₀ 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₀(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 ̃$ 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 ̃$ 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_d, we define the plasmoid half-width w as

where the prime denotes d/dz evaluated at z = 0. In Eq. (2), the fluctuation amplitude $B ̃$ 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_dξ] gives a fluctuation amplitude

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 ̃$ 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.³¹ 

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_d is given by

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_x. The criterion 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.⁹ 

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₀, the current sheet width a decreases exponentially at early times with $a ≃ a 0 e − t / τ$⁠, and a → a_∞ as t → ∞. We further assume that the upstream magnetic field B_x and the current sheet half-length L remain fixed in time and that the outflow speed v_o is the upstream Alfvén speed V_A. 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_∞ to be the Sweet-Parker width, i.e., $a ∞ = 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₀ = B_x and L₀ = 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_A, as follows: $t ̂ = t / τ A , τ ̂ = τ / τ A , k ̂ = kL , a ̂ = a / L$⁠, and $γ ̂ = γ τ A$⁠. The model equation in normalized variables is

where dL/dt = 0 is assumed.

We can solve Eq. (6) by the method of characteristics. The characteristic starting from $k ̂ = k ̂ 0$ at $t ̂ = t ̂ 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 $γ ̂$ as a function of the wavenumber and time. The linear growth rate of the tearing mode is governed by the tearing stability index $Δ ′$⁠,³³ 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 $Δ ′ > 0$ being satisfied. While the complete dispersion relation for resistive tearing modes can be expressed in terms of $Δ ′$⁠, the expression is involved and requires solving a transcendental equation.³⁴ Analytically tractable approximations can be obtained in two limits: The small- $Δ ′$ regime for short-wavelength modes and the large- $Δ ′$ regime for long-wavelength modes. In the small- $Δ ′$ regime (large k and $Δ ′ δ ≪ 1$⁠), the linear growth rate is given by³³ 

where $C Γ = ( π − 1 Γ ( 1 / 4 ) / Γ ( 3 / 4 ) ) 4 / 5 ≈ 0.953$⁠. In the large- $Δ ′$ regime (small k and $Δ ′ δ ≫ 1$⁠), the linear growth rate is approximately given by³⁴ 

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

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

and in the large- $Δ ′$ regime

The two asymptotic limits (12) and (13) can be smoothly connected by a function

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 − k ̂ 2 a ̂ 2 )$ in Eq. (12) for our analytical derivations. This approximation can be justified a posteriori by noting that the dominant mode wavenumber $k ̂ d$ follows the same scaling as the fastest growing mode wavenumber $k ̂ max$ [see Eq. (35) and the discussion thereafter] and that $k ̂ max a ̂ ≃ ( S a ̂ ) − 1 / 4$ (see below); therefore, in the high-S regime, the condition $k ̂ a ̂ ≪ 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 ̂ ≫ 1$ is satisfied. Putting these approximations together yields an approximate expression for the linear growth rate that is valid for $k ̂ a ̂ ≪ 1$⁠, i.e.,

We will apply this expression to Eq. (8) to derive scaling relations in Sec. IV.

The fastest growing mode occurs at the transition between the small- $Δ ′$ and the large- $Δ ′$ regimes, i.e., at $k ̂ a ̂ ≃ ( S a ̂ ) − 1 / 4$ where $γ ̂ s ≃ γ ̂ l$⁠. More precisely, the fastest growing wavenumber $k ̂ max$ is given by (see, e.g., Ref. 35)

and the corresponding growth rate is

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 ̂ d$⁠, the current sheet width $a ̂ d$⁠, the dominant mode wavenumber $k ̂ d$⁠, and the linear growth rate $γ ̂ 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.

Let the initial perturbation be $f 0 ( k ̂ 0 )$ at $t ̂ 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 ̂ d ≫ a ̂ ∞$⁠, because $a ̂ d$ follows a weaker scaling with respect to S [see Eq. (47)] than $a ̂ ∞$⁠, which is assumed to be the Sweet-Parker width $a ̂ SP = S − 1 / 2$⁠. Hence, we will ignore a_∞ in Eq. (5) and assume $a ̂ = a ̂ 0 e − t ̂ / τ ̂$ in the following analysis. The combined effect of mode stretching (⁠ $k ̂ = k ̂ 0 e − ( t ̂ − t ̂ 0 )$⁠) and current sheet thinning $( a ̂ = a ̂ 0 e − t ̂ / τ ̂ )$ implies that the variable φ decreases exponentially in time as

where the time scale $τ ̂ *$ is defined as

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

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

where

Here, ₂F₁ is the hypergeometric function³⁶ and

To proceed from Eq. (22), we further make the approximation that in the high-S regime, the contribution from the lower boundary, $G ( φ 0 ; τ ̂ )$⁠, is negligible compared to the contribution from the upper boundary, $G ( φ ; τ ̂ )$⁠, at the disruption time. This approximation can be justified a posteriori using the scaling relation of the disruption current sheet width $a ̂ d$ in Eq. (47). Because $a ̂ d$ decreases as S increases, it follows from the relation $φ / φ 0 = ( a ̂ / a ̂ 0 ) τ ̂ / τ ̂ *$ that in the high-S regime, the condition φ ≪ φ₀ 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 ( φ ′ ; τ ̂ )$ is monotonically decreasing with respect to $φ ′$⁠. Its overall behavior is well captured by two asymptotic expressions for small and large $φ ′$⁠, 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 ( φ ′ ; τ ̂ )$ is approximately

when $φ ′ ≲ O ( 1 )$⁠, and

when $φ ′ ≳ O ( 1 )$⁠. For the special case $τ ̂ = 1$ that yields $τ ̂ * / τ ̂ ♯ = 2 / 3$⁠, Eq. (28) is replaced by

Because φ₀ ≫ φ ∼ O(1), it follows from Eqs. (28)–(30) that $G ( φ 0 ; τ ̂ ) ≪ G ( φ ; τ ̂ )$⁠.

Ignoring $G ( φ 0 ; τ ̂ )$⁠, we can simplify Eq. (22) as

where

and

The solution of the model equation can now be written as

with $k ̂ 0 = k ̂ e t ̂ , h = τ ̂ ♯ C Γ 5 / 8 S − 1 / 2 a ̂ − 3 / 2 / 2$⁠, and $φ = C Γ − 3 / 4 S 1 / 5 k ̂ 4 / 5 a ̂$ being substituted into the right-hand side.

The dominant wavenumber $k ̂ 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 ̂ 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).³⁷ Let φ = φ_d correspond to the maximum of $g ( φ ; τ ̂ )$⁠. The value of φ_d depends on $τ ̂$ and, in general, has to be obtained numerically; typically, φ_d is an O(1) quantity. With the value of φ_d obtained, the dominant wavenumber $k ̂ d$ and linear growth rate $γ ̂ d$ are given by

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 ̂$ 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.⁹ 

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

where the magnetic fluctuation $B ̃$ 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 ̂ d$ and then approximate the integration by a Gaussian integral. Expanding g in the neighborhood of $k ̂ d$ yields

Hence, in the neighborhood of the dominant mode,

where

and $t ̂ d$ is the disruption time. Now we can evaluate the magnetic fluctuation $B ̃$ 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 ̂ d$

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

Equation (44) is the fundamental equation to determine the disruption current sheet width $a ̂ d$ for a general initial condition $f 0 ( k ̂ 0 )$⁠; its validity only requires that $f 0 ( k ̂ 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 ̂ 0 )$ without resorting to numerical methods. To make further progress, we assume that the initial condition takes a power-law form $f 0 ( k ̂ 0 ) = ϵ k ̂ 0 − χ$⁠. Then, Eq. (44) becomes

where

Equation (45) can be solved for $a ̂ d$ in terms of the lower branch W_–1 of the Lambert W function.³⁸ 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 ″ ( φ d )$ are typically O(1) quantities; therefore, the validity of the Gaussian integral approximation requires the condition $h ( a ̂ d ) ≫ 1$ being satisfied. Using Eq. (47) in Eq. (32), the condition is

This condition is satisfied only if $| Ξ | ≪ 1$⁠. Therefore, the noise amplitude ϵ must be sufficiently small for the approximation to be valid.

Because the condition $| Ξ | ≪ 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 ̂ d$

Here, we have ignored an O(1) factor $c χ ( θ c a 3 / 2 ) − 1 / θ$ within the logarithm, which can be attributed to the next order correction of the form $log ( log ( … ) )$⁠.

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

The calculation in Sec. IV A assumes that an initial perturbation is applied at $t ̂ = 0$⁠. If the current sheet width is sufficiently broad such that the linear growth rate $γ ̂ ( k ̂ ) < 1 / 2$⁠, and then the perturbation amplitude at wavenumber $k ̂$ 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 ̂ → ∞$⁠. 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 ̂$⁠) to the fluctuation amplitude. Here, we adopt the second approach: Instead of an initial perturbation, we now use $f 0 ( k ̂ )$ to describe the system noise, i.e., $f 0 ( k ̂ )$ represents the lower bound of $f ( k ̂$⁠). In this section, we consider the situation where the system noise provides the seed for the plasmoid instability.

Because $γ ̂$ must be greater than 1/2 for a mode to grow when we integrate along a characteristic in Eq. (8), the starting time $t ̂ 0$ is determined by the condition $γ ̂ = 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- $Δ ′$ regime.³⁹ Using the relations $k ̂ 0 = k ̂ e t ̂ − t ̂ 0$ and $a ̂ ( t ̂ 0 ) = a ̂ e ( t ̂ − t ̂ 0 ) / τ ̂$ in the condition using the small- $Δ ′$ dispersion relation

yields

hence, the wavenumber $k ̂ 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 = τ ̂ ♯ C Γ 5 / 8 S − 1 / 2 a ̂ − 3 / 2 / 2$⁠, and $φ = C Γ − 3 / 4 S 1 / 5 k ̂ 4 / 5 a ̂$ being substituted into the right-hand-side.

Equations (35)–(40) remain valid, but A in Eqs. (41)–(43) must be replaced by

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

If we assume a power-law system noise $f 0 ( k ̂ ) = ϵ k ̂ − χ$⁠, Eq. (59) becomes

where

The solution for $a ̂ d$ can be expressed in terms of the lower-branch Lambert function W₋₁ as

where

and c_a is defined in Eq. (48). Using the asymptotic expansion (52), the leading order approximation of $a ̂ 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(…)).

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 τ ≪ τ_A, such that the mode-stretching and advective loss due to the outflow have no time to take effect. Taking the limit $τ ̂ = τ / τ A → 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.

In the $τ ̂ → 0$ limit, the two time scales $τ ̂ *$ and $τ ̂ ♯$ become identical to $τ ̂$⁠. The variables θ and Ξ become

and

The disruption current sheet width is given by

with the leading order approximation