Marginal Stability of Sweet-Parker Type Current Sheets at Low Lundquist Numbers
aa r X i v : . [ phy s i c s . p l a s m - ph ] A p r Draft version April 11, 2018
Typeset using L A TEX default style in AASTeX62
Marginal Stability of Sweet-Parker Type Current Sheets at Low Lundquist Numbers
Chen Shi ( 时 辰 ), Marco Velli, and Anna Tenerani EPSS, UCLA, Los Angeles, CA 90095, USA (Received; Revised; Accepted)
Submitted to ApJABSTRACTMagnetohydrodynamic simulations have shown that a non-unique critical Lundquist number S c exists, hovering around S c ∼ , above which threshold Sweet-Parker type stationary reconnectingconfigurations become unstable to a fast tearing mode dominated by plasmoid generation. It is knownthat the flow along the sheet plays a stabilizing role, though a satisfactory explanation of the non-universality and variable critical Lundquist numbers observed is still lacking. Here we discuss thisquestion using 2D linear MHD simulations and linear stability analyses of Sweet-Parker type currentsheets in the presence of background stationary inflows and outflows at low Lundquist numbers ( S ≤ ). Simulations show that the inhomogeneous outflow stabilizes the current sheet by stretching thegrowing magnetic islands and at the same time evacuating the magnetic islands out of the currentsheet. This limits the time during which fluctuations which begin at any given wave-length can remainunstable, rendering the instability non-exponential. We find that the linear theory based on theexpanding-wavelength assumption works well for S larger than ∼ Keywords: magnetic reconnection, magnetohydrodynamics INTRODUCTIONMagnetic reconnection is a phenomenon where magnetic energy is converted into the thermal and kinetic energy ofthe plasma via an induction electric field, leading to a reconfigured magnetic field topology. The classic stationaryreconnection model was developed by Sweet and Parker (SP) (Parker 1957; Sweet 1958), who studied magnetic fieldannihilation in a two-dimensional current layer with a macroscopic half-length L . It predicts a reconnection rate R ∼ S − / where S is the Lundquist number defined as S = LV Au /η ( V Au and η are the upstream Alfv´en speed andthe magnetic diffusivity). For current sheets in space plasmas, S is usually very large (e.g. ∼ in the solar corona),and the typical time-scales are much longer than the duration of explosive events such as flares: thus, even though theinverse aspect ratio of the current sheet is small, a/L ≃ S − / ( a is the half-thickness of the current sheet), the SPmodel cannot explain explosive reconnection properly.However, simulations by (Biskamp 1986) showed that the SP stationary reconnecting current sheet becomes unstableto an extremely fast super-tearing, or plasmoid instability, once a critical value for the Lundquist number 10 ≤ S c ≤ is exceeded. Loureiro et al. (2007) pointed out, (after (Tajima & Shibata 2002)), that the maximum growth rateof the instability scales as ∼ S / . The positive scaling of the instability growth rate with Lundquist number, and itsdivergence in the ideal limit, led Pucci & Velli (2014) to point out that the only conclusion of such stability studiesshould be that large- S SP current sheets should not form in nature in the first place. They identified a critical aspectratio scaling
L/a ∼ S / above which current sheets become unstable even in the limit of ideal MHD. The criticalaspect ratio scaling separates current sheets for which the tearing instability diverges, from those in which the growthrate of the tearing mode goes to zero as S → ∞ . At the critical scaling the tearing instability becomes ideal – the Corresponding author: Chen [email protected]
Shi et al. growth rate becomes independent of the Lundquist number. Resistive MHD simulations by Tenerani et al. (2015) andLandi et al. (2015) have confirmed the existence and relevance of the critical aspect ratio even for nonlinear simulations.The works mentioned above all focus on large Lundquist numbers ( S ≥ ), though as mentioned above theearly MHD simulations suggested that at S ≤ the current sheet is actually stable and laminar (Biskamp 1986).Bulanov et al. (1978) proposed that, unlike static current sheets which are always unstable (albeit to the very slow,for macroscopic current layers, tearing instability), the inhomogeneous outflow of the SP configuration can have astabilizing effect leading to the Lundquist number threshold for the current sheet stability. Biskamp (2005) furtheranalyzed Bulanov’s model and estimated roughly the threshold for the aspect ratio L/a to be ∼ , correspondingto a Lundquist number S ∼ . However, Loureiro et al. (2007) claimed that the outflow only plays a mildlystabilizing role. Indeed it was further suggested (Loureiro et al. (2013)) that stabilization might be provided by amore mathematical criterion, namely that in order for the tearing instability to develop in a SP current sheet, theinner singular layer which tearing develops should be significantly narrower, i.e. by about a factor of 3, than theoverall current sheet thickness. Remark that this criterion might be important in invalidating the classical asymptoticexpansions (Furth et al. 1963) with which the analytical approximation of the dispersion relation for resistive modeswas originally derived, but arguably has no real physical meaning: indeed, as we will discuss further in the conclusions,it can not justify the observed SP stabilization.More recently the 2D linear MHD simulation by Ni et al. (2010) confirmed that, though the effect of flows is negli-gible when S is large, its importance increases with decreasing S . In their simulations, initial random perturbationsgrow linearly at first but then saturate. They determine a critical Lundquist number S ∼ S ≤ ). We use open boundary conditions in both the inflow and outflow directions in thesimulations. Our result shows that the inhomogeneous outflow stretches the magnetic islands forming and growing inthe current sheet while at the same time evacuating them from the simulation domain. Both stretching and evacuationplay a stabilizing role. In addition, the stretching effect changes the linear instability from the exponential behaviorof the static instability, to one in which growth is only exponential over a limited time-interval. We confirm that theoverall growth of fluctuations is limited, the reason being that the effective wave-number of unstable perturbationsdecreases over time limiting the overall period of instability.The following section describes our numerical method; Section (3) describes the numerical results. Sections (4)and (5) then describe SP stability while Section (6) shows numerical results for a different background configuration.Section (7) describes the relevant outlook of our results. NUMERICAL METHOD AND SIMULATION SETUPThe governing equations for the resistive stability of SP current sheets are given by the linearized resistive MHDequations (and the adiabatic equation for closure - this decouples the dynamic instability from potential thermodynamiceffects): ∂ρ ∂t + u · ∇ ρ + ρ ∇ · u + u · ∇ ρ + ρ ∇ · u = 0 (1a) ∂ u ∂t + u · ∇ u + u · ∇ u + 14 πρ [ ∇ b · B − B · ∇ b + ∇ B · b − b · ∇ B ] + ρ ρ u · ∇ u = − ρ ∇ p (1b) ∂φ ∂t − u × b − u × B − η ∇ φ = 0 (1c) ∂T ∂t + u · ∇ T + ( κ − ∇ · u ) T + u · ∇ T + ( κ − ∇ · u ) T = 0 . (1d)The variables with subscript “0” are the stationary background fields and those with subscript “1” are the perturbedfields which are evolved by the code. p = ρ T + ρ T is the perturbed pressure, ρ, T are density and temperature earing instability in 2D current sheet u is the velocity and the adiabatic index κ = 5 / φ is the perturbed magnetic flux function from whichthe perturbed magnetic field b is obtained b = ∂φ ∂y ˆ e x − ∂φ ∂x ˆ e y The simulation domain is a 256 ×
256 rectangular box whose x / y axis is along/across the current sheet. As mentionedin the introduction, the Lundquist number is defined with the half-length of the box L : S = LV Au η (2)In our simulations, the length of the box is fixed to be 2 ( − ≤ x ≤
1) which means all the lengths are normalized tothe half length of the current sheet. As we are studying the Sweet-Parker type current sheets, the half width of thecurrent sheet scales with S as: aL = S − (3)so a changes with the Lundquist number in our simulations. As a result, the width of the simulation domain isalso changing according to the Lundquist number in order to maintain proper range and resolution in y direction(max( | y | ) ∼ a ). Derivatives are calculated using the 6th order Compact Finite Difference (CFD) scheme (Lele1992). At the boundaries, projected characteristics are applied to keep the boundaries open (Landi et al. 2005) (seeAppendix A). We use a 4th order Runge-Kutta method to advance in time.The background fields are given by: B = B ( y )ˆ e x (4a) u = Γ x ˆ e x − Γ y ˆ e y (4b)where Γ is a constant controlling how fast the outflow is accelerated. For Sweet-Parker type current sheets, the outflowspeed is the upstream Alfv´en speed so we have Γ = V Au L = τ − A (5)where τ A is the Alfv´en crossing time of the current sheet. The plot of u is shown in the left panel of Fig (1). ρ isuniform and set as 1 / π so the local Alfv´en speed is simply V A ( y ) = B ( y ) (6) T ( x, y ) is calculated by T ( x, y ) = ¯ T −
12 ( u x + u y + B ) (7)so that the 0th order momentum equation is satisfied. ¯ T is a constant and set to be 2 in the simulations. Note thatself-consistency requires that ∇ × ( u · ∇ u − πρ B · ∇ B ) = 0, which is satisfied by Eq (4).The self-consistent equilibrium magnetic field B ( y ) is B ( y ) = ¯ B .
54 exp[ − ( y √ a ) ] Z y/ √ a e s ds (8)shown as the blue curve in the right panel of Fig (1). For this equilibrium, the amplitude of B has a maximumand then decreases toward 0 far from the center of the sheet. In other words, because of the background inflow, themagnetic field is concentrated near the midplane of the current sheet and we use the peak value ¯ B (set as 1) as thecharacteristic (“upstream”) magnetic field in this case. In Section (3) – (5), we adopt this equilibrium field for thesimulations and the linear stability analysis. In Section (6), we also show simulation results based on the widely-usedHarris current sheet field B ( y ) = ¯ B tanh( ya ) (9)shown as the orange curve in the right panel of Fig (1) where ¯ B (set as 1) is the asymptotic value of B .Note that both L and V Au are fixed to be 1 thus the Lundquist number is determined solely by the resistivity η inthe code. Except for the runs shown in Section (5), we initiate the simulations with white noise, i.e. perturbations ofthe magnetic flux φ uniformly distributed in k space with random phases peaking near the midplane of the currentsheet. Shi et al. −1.00 −0.75 −0.50 −0.25 0.00 0.25 0.50 0.75 1.00 x −0.100−0.075−0.050−0.0250.0000.0250.0500.0750.100 y Sp ee d −0.100 −0.075 −0.050 −0.025 0.000 0.025 0.050 0.075 0.100 y −1.00−0.75−0.50−0.250.000.250.500.751.00 B ( y ) / ̄ B self̄consistentHarris Figure 1.
The background fields in the simulations with S = 10 . Left: Streamlines of the background flow u . Right: Twotypes of B ( y ) used in the simulations. −0.10−0.050.000.050.10 y t = 0.400 −0.0040−0.0032−0.0024−0.0016−0.00080.00000.00080.00160.00240.0032−0.10−0.050.000.050.10 y t = 1.400 −4.0−3.2−2.4−1.6−0.80.00.81.62.43.2−1.00 −0.75 −0.50 −0.25 0.00 0.25 0.50 0.75 1.00 x −0.10−0.050.000.050.10 y t = 3.000 −1600−1200−800−400040080012001600 0.000.250.500.751.001.251.501.75 b y ( x = ) t = 1.400 −0.10 −0.05 0.00 0.05 0.10 y −20−1001020 u y ( x = ) Figure 2.
Left column: Snapshots of the perturbed magnetic flux function φ in the S = 10 run. Right column: Profiles of b y and u y along y -axis at t = 1 .
4. 3.
SIMULATION RESULTSIn this section we describe simulation results of the runs with the self-consistent background fields. Fig (2) leftcolumn shows the evolution of the perturbed magnetic flux function φ in the run S = 10 ( a = 0 . b y and u y along the y -axis at t = 1 .
4. As will be shownin Section (4), these are essentially identical to the eigenfunctions which may be calculated from the linear theory byincorporating the expanding wavelength assumption k ( t ) = k exp( − Γ t ); (10) earing instability in 2D current sheet k for the tearing mode in a current sheet with a linearly accelerating outflow u x = Γ x . Note, however,that in the 2D current sheet, it is impossible to strictly define a purely monochromatic wave because of the lack ofperiodicity along x . In any case, the 1D assumption takes the stretching effect of the background flow properly intoaccount: imagine two points separated by one wavelength λ ( t ) along the flow x ( t ) = x ( t ) + λ ( t ) (11)After time dt , the distance between the two points becomes λ ( t + dt ) = x − x + [ u x ( x ) − u x ( x )] dt = λ ( t ) + du x dx λ ( t ) dt (12)which means 1 λ dλdt = du x dx (13)and we easily get Eq (10). The result shown in Fig (2) is qualitatively consistent with the time-dependent wave numberassumption: we will discuss this in more detail in Section (4).We calculate the perturbed magnetic and kinetic energies inside the simulation domain E k = Z
12 ( u x + u y ) , E b = Z
12 ( b x + b y ) (14)as functions of time. These are shown in Fig (3): from top-left to bottom-right the panels display runs for S =10 , ,
300 and 100 respectively. One immediate result from Fig (3) is that the growth rate of the energy in theperturbed fields decreases with decreasing Lundquist number. At S ∼ x -component of the fields dominates over the y -component energy. This is very clear for runs S = 10 , and 300 where the ratios E kx /E ky and E bx /E by tend toincrease with time. This phenomenon can also be explained qualitatively by the stretching effect of the inhomogeneousoutflow. Consider the 1st order magnetic field b , which must be divergence free: ∂b x ∂x + ∂b y ∂y = 0 (15)If the perturbations are approximated by a mode with a time-dependent wave number (Eq (10)), we can estimate ∂b x ∂x ∼ k b x exp( − Γ t ) (16)and ∂b y ∂y ∼ b y δ (17)where δ is the thickness of the inner layer of the tearing mode (Loureiro et al. 2007) and is approximately constantwith time. We then get b y b x ∼ ( k δ ) exp( − Γ t ) (18)which means that the growth of b y is slower than that of b x . We can apply the same analysis to u assuming thatthe incompressible modes dominate.In practice, b y , as the reconnecting component of the field, is the more interesting quantity. Fig (4) shows thegrowth of b y in several runs. The upper panel shows the amplitude of b y as function of time for runs with different S .Here we use the maximum of | b y | inside the whole simulation domain to represent the amplitude of b y . For S < | b y | does not grow with time. For the runs with S > | b y | grows after an initial transient phase and then thegrowth gradually slows down until reaching a maximum. The two dots on each curve mark the start and the end ofthe growth phase. The saturation time corresponds roughly to the time when most of the magnetic islands are ejectedout. We use the ratio between the amplitudes of b y at the two points to measure the total growth of | b y | during Shi et al. t E n e r g y S = 10000 E k E b t −3 −2 −1 E n e r g y S = 1000 E kx E bx t −3 −2 −1 E n e r g y S = 300 E ky E by t −6 −5 −4 −3 E n e r g y S = 100
Figure 3.
Energy of the perturbed velocity and magnetic fields as functions of time for different S . Blue lines are magneticenergy and red lines are kinetic energy. Dashed lines are the energy in x components, dashed-dotted lines are the energy in y components and solid lines are the total energy. a simulation, which is shown as the blue dots in the lower panel of Fig (4) as a function of Lundquist number. Theorange dots in the same panel are the total growth of b y estimated through the linear theory which will be discussedin next section. At S = 1000, | b y | grows by a factor of 17 while at S = 10 it grows by a factor of ∼ . × . If wearbitrarily choose the criterion that | b y | grows by a factor of 10 in order to perturb the current sheet strongly, thecritical Lundquist number needs to be ∼ LINEAR STABILITY ANALYSIS OF THE 2D CURRENT SHEETSIf we use the same background flow as used in the simulations (Eq (4b)) and include the expanding effect of theoutflow on wavelengths by using a time-dependent wave vector as given by Eq (10), assuming that v and b are functionsof y only (rigorously they should be functions of both y and t ), of the form u y b y ! = − iv ( y ) b ( y ) ! exp (cid:2) ik ( t ) x + Z t γ ( k ( t ′ )) dt ′ (cid:3) (19)we can reduce Eq (1) to a set of ordinary differential equations y ( v ′′′ − k v ′ ) − ( γ + 1) v ′′ + ( γ − k v = k La [ B ( y )( b ′′ − k b ) − B ′′ ( y ) b ] (20a)( γ + 1) b − yb ′ = k La B ( y ) v + 1 S L a ( b ′′ − k b ) (20b)In deriving Eq (20) we have also assumed incompressibility and adopted the following normalizations: length to thehalf thickness of the current sheet a , magnetic field to the upstream magnetic field ¯ B , speed to the upstream Alfv´en earing instability in 2D current sheet t −4 −2 m a x ( | b y | )( t ) S=10000S=6000S=3000S=1000S=500S=300S=1005001000 3000 6000 10000 S g r o w t h o f | b y | SimulationLinear theory
Figure 4.
Upper panel: Maximum amplitude of b y in the simulation domain as functions of time for several runs. The twodots on each curve mark the start and the end of the growth phase. Lower panel: Blue dots: growth of | b y | during the growthphase of the simulations (the ratio between max( | b y | ) values at the two dots on each curve in the upper panel) as a functionof S . Orange dots: growth of | b y | predicted by the linear theory. See G ( S ) (Eq (26)) and the discussion in Section (4). speed V Au , time to Alfv´en time τ A = L/V Au where L is the half length of the current sheet. S is defined as LV Au /η and we study Sweet-Parker type current sheets so that a/L = S − / . Eq (20) is a 1D eigenvalue problem with theboundary conditions that v and b vanish at infinity. We numerically solve this eigenvalue problem (see Appendix B)and present the result in Fig (5). The upper panel of Fig (5) is the dispersion relation γ ( k ). The shape of the γ ( k )curves is similar to that of the dispersion relation of classic tearing mode: there is a fastest growing mode k f for eachLundquist number and the growth rate decreases for larger or smaller k . Unlike the tearing mode in static currentsheets whose growth rate is never negative, now part of the γ ( k ) curve goes below the γ = 0 line. For S .
70 thewhole curve is below γ = 0 which means that in this configuration of the background fields, the tearing mode does notgrow at all for S . γ is ( γ + 1) v ′′ . In the same equation for the tearing mode inside a static current sheet, the termbecomes γv ′′ . That is to say γ ∼ γ s − γ is the growth rate with the presence of flow and γ s is the growth rate without the flow. The difference 1(Γ)comes from the time-derivative of k in deriving Eq (20). This means that in the current sheet with the backgroundflow, stretching of the magnetic islands due to the inhomogeneous outflow slows down their growth. Compared withthe equations derived by (Bulanov et al. 1978) (who takes into account only the outflow), Eq (20) displays two majordifferences: (1) the coefficient in front of v ′′ changes from γ +2 to γ +1. (2) y -convection terms appear ( y ( v ′′′ − k v ′ ) and Shi et al. yb ′ ). In Bulanov’s equation, the constant 2(Γ) consists of one Γ which comes from dk/dt (expansion of the wavelength)and another one Γ coming from du x /dx (expansion of the fluid volume) so their model is even more stable than Eq(20). However, when the inflow is included, the du x /dx is cancelled out by du y /dy because of incompressibility. Infact, the background configuration considered by (Bulanov et al. 1978) is not self-consistent because their u = Γ x ˆ e x is compressible while they still assumed a uniform density. The newly-included y -convection terms, especially the term y ( v ′′′ − k v ′ ), though not important at large S , affects the calculated γ at low S ( . ) when y ( v ′′′ − k v ′ ) becomescomparable to ( γ + 1) v ′′ . We also solved Eq (20) without the two y -convection terms (not shown here) and it turnsout that the growth rate becomes larger, i.e., the y -convection terms somewhat contribute to the stabilization of thecurrent sheet. But as will be discussed further below, at very low S the linear analysis method becomes invalid so weare unable to evaluate exactly how much stabilization the y -convection brings.In the middle panel of Fig (5) we compare the eigenfunctions calculated through Eq (20) with the profiles takenfrom the simulation for S = 10 . The blue lines are the calculated eigenfunctions at ka = 0 .
15, i.e., the wavelength λ ≈ . u y and b y along the y axis taken from the S = 10 run at t = 1 . ≈ .
4. The amplitudes of the eigenfunctions and the cuts fromthe simulation are adjusted to be the same. We can see that the shapes of the eigenfunctions by the linear theory arevery close to the profiles taken from the linear simulation, confirming the validity of the linear theory in this case.Based on the linear theory, we are able to explain the saturation phenomenon observed in Fig (4). For a mode withwave number k ( t ) at time t , the instantaneous growth rate is γ ( k ( t )). As time increases, k ( t ) decreases and thus theinstantaneous growth rate moves along the γ − k curve toward left. After k ( t ) passes the fastest-growing wave number k f , the growth rate starts to decrease and eventually goes to 0 (even below 0). To verify this theory, we first estimatethe dominant mode, i.e. the mode whose amplitude is the largest, at a certain time and compare it with the simulationresult. We assume that the initial perturbations are uniformly distributed in k -space. At time t , the wave number k = k ( t ) corresponds to the initial wave number k = k ( t ) e t (in our normalization Γ = 1). At any time 0 ≤ t ′ ≤ t thismode corresponds to k ′ = k e − t ′ = k ( t ) e t − t ′ (22)and the instantaneous growth rate at t ′ is γ ′ = γ ( k ′ ) = γ [ k ( t ) e t − t ′ ] (23)In order to get the amplitude of the mode k ( t ), we integrate along the γ − k curve from k to k ( t ): | b y | ( t, k ( t )) = exp (cid:2) Z t γ ( k e t − t ′ ) dt ′ (cid:3) (24)and then we are able to find the wavenumber k d ( t ) whose amplitude is the largest. In the left panel of Fig (6),the orange curve shows the estimated dominant wave number through Eq (24) for S = 10 and the blue dots showthe dominant wave number calculated from the simulation S = 10 . In analyzing the simulation result, we first cut b y ( x, y, t ) along x axis, i.e. the midplane of the current sheet, and then multiply it by a Tukey window to removethe effect of non-periodicity along x before we apply Fourier transform to it. The steps in the blue dots are due tothe resolution in k space when we do Fourier transform. We shift the orange curve by ∆ t = 0 . t & .
0. The large differencebetween them for t . . b y as a function of time from the simulation, we are able to calculate the instantaneous growth rate at a given timeby γ ( t ) = 1 | b y | d | b y | dt (25)and then relate it with the calculated dominant wave number at the same time to get a γ − k curve, which is shown inthe right panel of Fig (6). Again the blue dots are the simulation result and the orange curve is the dispersion relationgiven by the linear theory (the curve taken from Fig (5)). The theory and the simulations agree with each other verywell, verifying the physical picture that γ moves along the γ − k curve.We estimate the total growth of | b y | using the dispersion relation shown in Fig (5) by doing the integral along the γ − k curves: G ( S ) = exp (cid:20) Z t c γ ( k e − t ′ | S ) dt ′ (cid:21) (26) earing instability in 2D current sheet ka −101234567 γ τ A S = 70S = 500S = 1000S = 3000S = 6000S = 10000−10 −5 0 5 10 y/a b −10 −5 0 5 10 y/a −1.0−0.50.00.51.0 v −10 −5 0 5 10 y/a b −10 −5 0 5 10 y/a −1.0−0.50.00.51.0 v Linear TheorySimulation
Figure 5.
Linear stability analysis based on Eq (20). Upper panel: Dispersion relation γ − k for different S . Middle panel:Blue solid line: eigenfunctions calculated through linear analysis for S = 10 and ka = 0 .
15. Orange dashed line: profiles of u y and b y taken from the simulation when ka ≈ .
15 (right column of Fig (2), t = 1 . S = 500 and ka = 0 .
20 ( t ≈ . Shi et al. t k a SimulationLinear theory 0.0 0.1 0.2 0.3 0.4 0.5 ka γ τ A SimulationLinear theory
Figure 6.
Left panel: Dominant wave number as a function of time. Blue dots are calculated from the simulation S = 10 andthe orange curve is the estimate through the linear theory. The steps in the blue dots are due to the finite resolution in k whenwe apply Fourier transform. The orange curve is shifted in time by ∆ t = 0 . ∼ .
1. Right panel: Dispersion relation γ − k . Blue dots are calculated from the simulation and theorange curve is the linear theory. where t = 0 is at the right zero point of the γ − k curve and t c is when k exp( − t c ) = π/L , i.e., when the wavelengthreaches the length of the current sheet (generally very close to the left zero point of the curve). Note that Eq (26) isactually an estimate of the upper limit of the total growth of | b y | because we integrate through the whole positivepart of each γ − k curve, thus the estimate is larger than the real growth. The calculated G ( S ) is shown as the orangedots in the lower panel of Fig (4). We can see that it is in accordance with the simulation though a little bit largerand in general the theory result and the simulation result differ by a factor of smaller than 2.Although the linear stability analysis based on the k = k exp( − Γ t ) assumption works quite well for S & , wemust point out that, for small Lundquist numbers, this assumption becomes invalid because of the small growth rate( γ . τ − A ), i.e., when the growth time scale for a magnetic island becomes larger than the time scale for it to beejected out of the current sheet. Actually, when γ < Γ(= τ − A ), we cannot assume that the functions v and b in Eq(20) are independent of time anymore. Note that in Eq (20) k is a function of time so mathematically the assumptionEq (19) is not rigorous. Stated otherwise, we are unable to pose a 1D eigenvalue problem for a 2D Sweet-Parker typecurrent sheet at very small Lundquist numbers. In the lower panel of Fig (5), we compare the simulation result withthe eigenfunctions for S = 500 and ka = 0 .
20. Clearly the simulation result deviates from the eigenfunctions. There isno doubt that as S becomes smaller the difference between the simulation and the approximate 1D linear theory willbecome larger. The linear theory predicts a threshold S c ≈
70 below which b y does not grow while from the upperpanel of Fig (4) we see that the threshold is at least ≈ THE EFFECT OF THE INITIAL PERTURBATIONIn this section we discuss the stability of the current sheet in terms of initial conditions, i.e. in terms of the initialperturbation. Our simulation is linear so that the amplitude of the initial noise is arbitrary. In Section (3) we chosethe criterion for an unstable current sheet to be that | b y | grows by a factor of 10 which gives a critical Lundquistnumber ∼ | b y | to grow sufficiently that the sheet is strongly affected, i.e. nonlinear evolution is triggered, at times comparable totheir evacuation from the domain.Because the configuration of the current sheet in this study is fully 2D and the simulation is open-boundary, it isnatural to propose that the shape of the initial perturbation also affects the stability. To be more specific, where theinitial perturbation is excited on the current sheet is important because the generation location of a magnetic islanddetermines how long the island grows before it is ejected out of the current sheet. Here the location refers to thedistance to x = 0, i.e. the location along the outflow. Based on the S = 10 run, we carried out three more runswith identical background configuration but for each run we modulate the initial perturbation by multiplying it with earing instability in 2D current sheet t −5 −3 −1 m a x ( | b y | )( t ) x c =0.0x c =0.4x c =0.8no modulation−1.00 −0.75 −0.50 −0.25 0.00 0.25 0.50 0.75 1.00 x m o d u l a t i o n x c =0.0x c =0.4x c =0.8 Figure 7.
Left: Modulation functions multiplied to the initial perturbations in the three runs with the identical configurationsas the run S = 10 shown in Fig (2). Right: Amplitudes of b y as functions of time in the three runs with the modulationfunctions shown in the left panel. The red dashed line is the original run without modulation, i.e. the S = 10 curve in theupper panel of Fig (4). a gaussian function: f ( x | x c ) = 12 (cid:2) exp( − ( x − x c ) σ ) + exp( − ( x + x c ) σ ) (cid:3) (27)where x c = 0 .
0, 0 . . σ = 0 .
1. The 3 modulation functions are plotted in the left panel of Fig (7). Throughthe modulation functions we are able to localize the initial perturbation around ± x c . We calculate the amplitudes b y as functions of time for each of the three runs and present the results in the right panel of Fig (7). The red dashedline shows the result of the original run without modulation, i.e. the S = 10 curve in the upper panel of Fig (4).We see that the curve of the x c = 0 . x c = 0 . | b y | is almost the same as the x c = 0 . t . .
0) but later the growth rapidlystops. As a result, the total growth of | b y | in this run is a factor of ∼ . × rather than the factor 3 . × ofthe x c = 0 . x c = 0 .
8, there is nearly no growth phase for | b y | . These results agree with the idea that theevacuation time of the magnetic island is a key factor in determining how much the island can grow. Assuming thata magnetic island is generated at x = x c , we can estimate the evacuation time of the magnetic island τ e ( x c ) = Z x c dxu x ( x ) = Z x c dxx = ln( 1 x c ) (28)We then get τ e (0 .
4) = 0 .
92 and τ e (0 .
8) = 0 .
22 which agree with the simulation results. SIMULATION RESULTS WITH THE HARRIS CURRENT SHEET MODELWe carry out the same simulations as those in Section (3) but change B to the Harris current sheet field (Eq (9))and present the results in Fig (8) and (9). We also carry out the simplified 1D linear stability analysis with the Harriscurrent sheet field, shown in Fig (10) where the dashed line is the result for a static Harris current sheet and is plottedto show the stabilizing effect of the background flow (see Eq (21) and the discussion in Section (4)). Note that weare still assuming the background fields are stationary though technically speaking the velocity field Eq (4b) and theHarris magnetic field Eq (9) do not satisfy the induction equation. Indeed the 0th order fields should vary with timeand their evolution will hence affect the evolution of the 1st order quantities. In this section, however, we ignore thenon-self-consistency of the background fields to illustrate the dependence of the current sheet stability on differentbackground magnetic fields. In reality, this implies that our results are physically significant only for growth rateslarge enough ( γτ A >
1) where the evolution of the background fields is negeligible.2
Shi et al. −0.10−0.050.000.050.10 y t = 0.400 −0.0008−0.0006−0.0004−0.00020.00000.00020.00040.0006−0.10−0.050.000.050.10 y t = 1.400 −0.040−0.032−0.024−0.016−0.0080.0000.0080.0160.0240.032−1.00 −0.75 −0.50 −0.25 0.00 0.25 0.50 0.75 1.00 x −0.10−0.050.000.050.10 y t = 3.000 −1.6−1.2−0.8−0.40.00.40.81.2 0.020.030.040.05 b y ( x = ) t = 1.400 −0.10 −0.05 0.00 0.05 0.10 y −0.15−0.10−0.050.000.050.100.15 u y ( x = ) Figure 8.
Left column: Snapshots of the perturbed magnetic flux function φ in the S = 10 run with a Harris current sheetfield. Right column: Profiles of b y and u y along y -axis at t = 1 . Compared with the self-consistent B case, the growth rate with B ( y ) = ¯ B tanh( y/a ) is smaller. By comparingFig (8) with Fig (2), we see that the magnetic islands are much thinner with the self-consistent field because of themore concentrated B and that is why the growth rate with the self-consistent B is larger. The growth of | b y | is alsomuch smaller than before: at S = 1000, | b y | grows by a factor of 1 . S = 10 it grows by a factor of 400 (3 . × for self-consistent field). If we adopt the same criterion for an unstable currentsheet as in Section (3) that | b y | grows by a factor of 10 , the critical Lundquist number is at least ∼ CONCLUSIONIn this study, we have explored the tearing instability inside 2D Sweet-Parker type current sheets at low Lundquistnumbers using both 2D linear open-boundary MHD simulations and a simplified 1D linear stability analysis. Wefind that in the linear stability analysis, the assumption of a time-dependent wave number k = k e − Γ t proposed by(Bulanov et al. 1978) works quite well for S & . As predicted by the linear theory, the simulations show that theinhomogeneous background outflow stretches the growing magnetic islands and the stretching slows down the islandgrowth. On the other hand, because of the finite size of the current sheet, the magnetic islands are evacuated out ofthe current sheet within a finite time (several Alfv´en times) and thus they cannot grow infinitely (as would happen ina periodic current sheet). In the simulations initialized with white noise we clearly observe a saturation phase afterthe growing phase, corresponding to the time when most of the magnetic islands are ejected out. As a consequence,the excitation location of the initial perturbation is an important factor in determining how much perturbations cangrow.We estimate through the simulations the critical Lundquist number above which the current sheet can be viewed as“strongly perturbed”. The result is dependent on the shape of the background magnetic field B ( y ) we use. Generally,the Harris current sheet field is more stable than the self-consistent background field (Eq (8)). Based on the criterionthat | b y | grows by a factor 10 during the growth phase, the Lundquist number threshold can be ∼ ∼ B ( y ) is used. We may say that under this specific criterion the threshold S c is several thousands. earing instability in 2D current sheet t −6 −5 −4 −3 −2 −1 m a x ( | b y | )( t ) S=10000S=6000S=3000S=2000S=1000S=500S=300S=1001000 2000 3000 6000 10000 S g r o w t h o f | b y | SimulationLinear theory
Figure 9.
Same as Fig (4) except for the use of the Harris current sheet field.
However, no universal criterion which is independent of context can be defined, since amplitude and localization offluctuations in the initial current are fundamental to the subsequent evolution.We find that the flow structure of the SP sheet adequately accounts for the low S stabilization without the need ofother factors, such as proposed by (Loureiro et al. 2013). Their suggested instability threshold, given by δa = f (29)where δ is the thickness of the inner layer, a is the thickness of the current sheet, and f ≃ / S compatible with findings of some simulations. On the other hand, the scaling of the stabilization criterionwith current sheet aspect ratio appears to go in the wrong direction. Indeed, consider a current sheet whose aspectratio is aL ∼ S − α . (30)The asymptotic scaling of δ then becomes (see e.g. (Tenerani et al. 2016)) δa ∼ S − (1 − α ) (31)and the threshold Lundquist number S c from the criterion above would become S c = ( 1 f ) − α (32)For a Sweet-Parker type current sheet ( α = 1 /
2) and f = 1 /
3, Eq (32) gives S c ≈ . × , close to 10 . However, S c given by Eq (32) is an increasing function of α , which means that the thinner a current sheet the more stable it would4 Shi et al. ka −1012345 γ τ A S = 330S = 600S = 1000S = 2000S = 6000S = 10000S = 10000, no flow−10 −5 0 5 10 y/a b −10 −5 0 5 10 y/a −1.0−0.50.00.51.0 v −10 −5 0 5 10 y/a b −10 −5 0 5 10 y/a −1.0−0.50.00.51.0 v Linear TheorySimulation
Figure 10.
Linear stability analysis based on Eq (20) and B = tanh( y/a ). Upper panel: Dispersion relation γ − k for different S . Middle panel: Blue solid line: eigenfunctions calculated through linear analysis for S = 10 and ka = 0 .
15. Orange dashedline: profiles of u y and b y taken from the simulation when ka ≈ .
15 (right column of Fig (8), t = 1 . S = 500 and ka = 0 .
20 ( t ≈ . earing instability in 2D current sheet < ) asymptoticscalings really do not work in the first place. On the other hand, we have demonstrated that velocity fields indeedprovide the right criteria for stability.In terms of astrophysical structures, such as for example coronal magnetic fields or prominence eruptions, fluctuationsare always present, propagating from other coronal regions or from the photosphere. As shown by Tenerani et al.(2016), at very high Lundquist numbers 2D current sheets tend to disrupt in a quasi self-similar way, in which caseflurrys of smaller scale current sheets are generated iteratively. At each step, the Lundquist number decreases, andfinally they will cross the critical S at which further disruption is quenched. However, at those scales where S becomessmall enough, kinetic effects, explored in (Pucci et al. 2017), become fundamental. We plan to address this questionin subsequent works.We would like to thank Fulvia Pucci for many useful discussions. This research was supported by the NSF-DOEPartnership in Basic Plasma Science and Engineering award n. 1619611 and the NASA Parker Solar Probe ObservatoryScientist grant NNX15AF34G. This work used the Extreme Science and Engineering Discovery Environment (XSEDE)(Towns et al. 2014) Comet at the San Diego Supercomputer Center through allocation TG-AST160007. XSEDE issupported by National Science Foundation grant number ACI-1548562.APPENDIX A. CHARACTERISTIC FORM OF THE LINEAR MHD EQUATIONSThe linear ideal MHD equation set is a partial differential equation system and can be written in the form: ∂ U ∂t + A ∂ U ∂x + C ∂ U ∂y + D = 0 (A1)where U is the vector of the 1st order variables: U = ( ρ , u x , u y , u z , b x , b y , b z , T ) (A2)The coefficient matrices A and C include only the 0th order quantities and the matrix D = FU contains otherterms which do not include any derivatives of the 1st order variables. Compared with the case of the nonlinear MHDequations, the first 3 terms on the L.H.S. of Eq (A1) are unchanged except that the derivatives are for 1st ordervariables only and the coefficients are 0th order. The matrix D now contains more terms than the nonlinear casebecause all the derivatives of the 0th order quantities are included in it, e.g. [ u · ∇ ρ + ( ∇ · u ) ρ ] in Eq (1a). Notethat D does not affect the projection of characteristics because it does not have any “waves” of U .The calculation of the characteristics is almost the same as shown in (Landi et al. 2005). For example, along x , wefirst diagonalize the matrix A A = SΛS − (A3)where the diagonal matrixΛ = diag { u x + f x , u x + a x , u x + s x , u x , u x − s x , u x − a x , u x − f x } (A4)consists of the speeds of the 7 MHD modes (projected along the x direction). u x , s x , a x , f x are the local entropymode speed, slow magnetosonic wave speed, Alfv´en speed and fast magnetosonic wave speed projected along x . Notethat when we calculate the characteristics along x , there is no need to consider b x because of the divergence-freecondition. Then we are able to decompose A ∂ U /∂x into the linear combinations of the 7 characteristics. The resultis quite lengthy and is not shown here, but note that one can write down the result easily by adding subscript “0” toall the variables in the coefficients and subscript “1” to all the variables in the derivatives for the equations shown inthe appendix of (Landi et al. 2005).In order to impose open-boundary condition on a certain side of the simulation box, we decompose the correspondingderivatives into the linear combination of the characteristics, e.g., decompose the x -derivative terms at the right/leftboundaries. For each of the 7 characteristics, we decide whether it is propagating inward or outward of the domainby looking at its mode speed. If it is outgoing, we keep this characteristic; otherwise we set it as 0. By doing this weensure that no information is going into the simulation domain and the boundary is “open”.6 Shi et al.
When the resistivity is non-zero, the 2nd order derivative term will cause some reflections at the boundaries. Butpractically, as the resistivity is quite small, the resistive term does not influence the simulation significantly at theboundaries. B. SERIES-EXPANSION METHOD FOR SOLVING THE ORDINARY DIFFERENTIAL EQUATION SETEq (20) is a 5th order homogeneous boundary value problem (BVP) with a regular singular point y = 0. Thewidely used methods for BVPs, e.g. the shooting method, usually involve the integrals of the ODEs over the domain.However, for ODEs with singular points such as Eq (20) it is hard to evaluate the integrals. Thus we adopt the methodof series-expansion to solve Eq (20). For simplicity, we write Eq (20) in the form L ( v, b | γ, k, S ) = 0 (B5)where L is the linear operator corresponding to the ODE set with γ, k, S (and B ( y )) as parameters. Our goal is,by imposing proper boundary conditions, to determine the dispersion relation γ ( k ) and the solution v ( y ) and b ( y )simultaneously.Before discussing the series expansion method, we would like to first clarify the boundary conditions. The equationis 5th order so 5 boundary conditions are needed. We set the domain to be y ∈ [0 , L y ] where L y = 17 so that the rightboundary is far enough to impose the asymptotic boundary conditions. We expect v to be odd and b to be even in y so that v (0) = 0 , b ′ (0) = 0 (B6)For large enough y , B ( y ) (Eq (8)) can be approximated by B ( y ) ≈ . × y + O (1 /y ) (B7)Far from the center of the current sheet, the perturbations asymptote to 0. Plugging Eq (B7) in Eq (20) and eliminating v , we get a 5th order ODE for b , which by keeping only the 5th and 4th order derivative terms a priori, becomes b (5) + ( yb ) (4) = 0 (B8)which gives the solution b ( y ) ∝ exp( − y ) (B9)and we further get v ( y ) = 11 . × S − k + γ + 2 k yb ( y ) (B10)The 3 boundary conditions at the right boundary are then given by b ′ + L y b = 0 (B11a) v ′ + ( L y − L y ) v = 0 (B11b) v − . S − k + γ + 2 k L y b = 0 (B11c)Eq (B6) and (B11) are the 5 boundary conditions for the BVP. Note that Eq (B11) are not equivalent to open-boundaryconditions adopted in the simulations unless L y → ∞ where all 1st order quantities vanish. But as the simulation boxis large enough max( | y | ) ∼ a and the initial perturbations are localized around y = 0, the difference in the boundaryconditions is negligible and we are still allowed to make comparisons between the linear theory and the simulations.We expand v ( y ), b ( y ) and B ( y ) in power series of y around a point y c v ( y ) = N X n =0 v n ( y − y c ) n , b ( y ) = N X n =0 b n ( y − y c ) n , B ( y ) = N X n =0 B n ( y − y c ) n (B12) earing instability in 2D current sheet v n and b n are quantities to be determined and N = 100 is the cut-off order of the series. B n are given by theTaylor expansion of B ( y ) around y c . Once we get the values of v n and b n , the solution v ( y ) and b ( y ) can be recovered.We plug Eq (B12) into Eq (20) and get the recurrence relations by matching the coefficients in front of ( y − y c ) n : y c ( n + 3)( n + 2)( n + 1) v n +3 = ( n + 2)( n + 1)( γ + 1 − n ) v n +2 + k y c ( n + 1) v n +1 − ( γ − − n ) k v n + kS n X l =0 (cid:2) ( l + 2)( l + 1) B n − l b l +2 − k B n − l b l − ( n − l + 2)( n − l + 1) B n − l +2 b l (cid:3) (B13a)( n + 2)( n + 1) b n +2 = ( γ + 1 − n + k ) b n − y c ( n + 1) b n +1 − kS n X l =0 B n − l v l (B13b)From Eq (B13) it is obvious that there are only 5 free parameters ( v , v , v , b , b ): all the v n for n ≥ b n for n ≥ v ( y ) and b ( y ): v ( y ) = v ( y | v , v , v , b , b ) , b ( y ) = b ( y | v , v , v , b , b ) (B14)Then the 5 boundary conditions can be written as a linear algebraic equation set M ( γ, k, S ) V = 0 (B15)where V = ( v , v , v , b , b ) T and M is a matrix which is a function of γ , k , S . In order to get a nontrivial solution,we need det( M ) = 0 (B16)which gives the dispersion relation γ ( k, S ) and then we get v ( y ) and b ( y ).In practice, we must do the expansion at multiple points ( y c ) to ensure convergence of the series. If we use N c pointsto expand the functions, we need 5 N c equations to determine the solution. Apart from the 5 boundary conditions Eq(B6) and (B11), the other 5( N c −
1) equations come from the continuities of v ( y ), v ′ ( y ), v ′′ ( y ), b ( y ), b ′ ( y ) at N c − y m ), each of which lies between two neighbouring expansion points. In solving Eq (20) we use 11expansion points. REFERENCES), each of which lies between two neighbouring expansion points. In solving Eq (20) we use 11expansion points. REFERENCES