THREE-DIMENSIONAL SPONTANEOUS MAGNETIC RECONNECTION

Current sheets are abundant in magnetized plasmas. Similar to thin vortices of hydrodynamics, they are naturally created by the nonlinear evolution of the conductive fluid (Parker 1994; Biskamp 2000; Priest & Forbes 2000). Magnetic X-points naturally evolve into current sheets due to currents' mutual attraction, creating the so-called Y-point configuration (Figure 1). Perhaps one of the most conspicuous phenomena associated with current sheets in plasmas are solar flares, bursts of radiation of up to $6\times {10}^{33}$ erg in X-rays. Following a big solar flare, coronal mass ejections (CME) occur, hinting to the global rearrangement of the magnetic field, which is called magnetic reconnection. Another well-known process is a magnetospheric storm, which is a perturbation of magnetosphere, associated with reconnection in the magnetotail. While CME demonstrates that there was a topological rearrangement of the magnetic field, careful observations near the flare site typically estimate the rate of inflow of magnetic field lines, called the reconnection rate, to a 0.001–0.1 fraction of the Alfvén speed ${v}_{{\rm{A}}}=B/\sqrt{4\pi \rho }$ (see, e.g., Dere 1996).

Zoom In Zoom Out Reset image size

In a well-conductive plasma, one might expect that current sheets are non-dissipative and, therefore, invisible. Indeed, the Sweet–Parker (SP) model (Parker 1957; Sweet 1969) predicts a very low reconnection rate for most astrophysical and space magnetic configurations. The dimensionless number characterizing plasma conductivity is the Lundquist number $S={{Lv}}_{{\rm{A}}}/\eta$, where L is the length of the layer and η is magnetic diffusivity. High Lundquist number means the magnetic resistive decay time, ${L}^{2}/\eta$, is much larger than the Alfvén crossing time, $L/{v}_{{\rm{A}}}$. For laminar thin current sheets, the SP model (Parker 1957; Sweet 1969) predicts a reconnection rate of ${v}_{{\rm{A}}}/\sqrt{S}$. This enhancement compared to resistive diffusion is due to the fact that magnetic field diffuses only through a thin width of the current sheet $L/\sqrt{S}$. This speed, however, is extremely low for most astrophysical and space magnetic configurations. The SP model, in the limit of very high S, becomes consistent with the so-called frozen-in condition of ideally conductive fluids. The SP prediction, therefore, contradicts the idea that a discontinuity in the magnetic field may result in an arbitrary reconnection rate independent of resistivity, as in Syrovatskii's model (Syrovatskii 1971). The search for fast reconnection has shifted toward microscopic effects beyond magnetohydrodynamics (MHD), e.g., effects in collisionless plasmas (Drake et al. 2006; Che et al. 2011; Daughton et al. 2011). Alternative approaches suggest that in the presence of turbulence the magnetic field lines will be stochastic (Lazarian & Vishniac 1999; Kowal et al. 2009; Eyink et al. 2011, 2013) which would lead to fast reconnection, having implications for particle acceleration (Lazarian et al. 2008; Beresnyak & Lazarian 2015; Beresnyak & Li 2016). Studies of the resistive tearing instability in a thin current sheet (Biskamp 1986) surprisingly resulted in a conclusion that it becomes faster and not slower with decreasing resistivity (Loureiro et al. 2007; Huang & Bhattacharjee 2010; Loureiro et al. 2012) at sufficiently high S. It has become clear that the SP model is problematic at high S because thin SP current sheets are unstable above the critical Lundquist number $S={{Lv}}_{{\rm{A}}}/\eta \sim {10}^{4}$. In a two-dimensional (2D) resistive reconnection scenario, a secondary instability of the current sheet between magnetic islands will result in a resistivity-independent reconnection (Uzdensky et al. 2010). 2D MHD simulations measured reconnection speeds around $0.01\div0.03{v}_{{\rm{A}}}$ (Huang & Bhattacharjee 2010; Loureiro et al. 2012) and observed hierarchical formation and ejection of plasmoids. Plasma simulations demonstrated that collisionless thin current layers are also unstable (Daughton et al. 2009).

In this paper we conduct 3D simulations of a thin current sheet with significant imposed mean field, starting with oblique tearing and developing into nonlinear phase which we call spontaneous turbulent reconnection.

In what follows, Section 2 describes the simulation setup, Section 3 offers an overview of the results of the simulations in terms of bulk average quantities, such as total energy budget and its evolution, and Sections 4 and 5 describe local properties of turbulence, spectrum, and anisotropy. Section 6 comments on the global nature of the perturbation of the magnitude of the magnetic field (slow mode). Section 7 discusses in some detail the differences between the reconnection picture in 2D and 3D cases, Section 8 proposes a phenomenological model for the reconnection rate, Section 9 discusses implications for particle acceleration, and Section 10 points out that a thin current layer could be viewed from an electromagnetic viewpoint as having non-zero resistance per unit length, even in the limit of vanishing resistivity. Section 11 is a discussion.

One of the simplest setups to study nonlinear development of tearing is a periodic setup with the mean field B_z0 threading the box and the reconnecting field $\pm {B}_{y0}$ changing sign in the x direction. We also consider the incompressible case, in which the problem has only two defining dimensionless numbers: (1) the Lundquist number S; and (2) the ratio ${B}_{y0}/{B}_{z0}$.

This very simple geometry physically corresponds to the initial (albeit already nonlinear) development of tearing before the outflow becomes important. The cartoon in Figure 1 illustrates the typical progression of high-S spontaneous reconnection by showing a cut perpendicular to the current and global mean field direction: the magnetic configuration with the X-point (I) may develop into a thin current sheet (II), the latter develops instability, the instability becomes nonlinear and produces a turbulent current layer (III) which later expands and produces the classic picture with inflow and outflow (IV). We study the initial regime before the outflow becomes important, designated as regime III in Figure 1, which is especially interesting because it has the highest volumetric dissipation rate.

We use all-periodic setups, so we actually simulated two current layers, and all the results are the average between the properties of these two layers. Figure 2 demonstrates the setup and shows a simulation snapshot during the development of the nonlinear phase, with the magnitude of B shown as grayscale on the surface of the box. One of the numerical challenges in studying 3D spontaneous reconnection was to break through the critical Lundquist number barrier of 10⁴, which require sufficiently big boxes. In simulations with imposed large-scale perturbation, such as Daughton et al. (2011), Oishi et al. (2015), and Huang & Bhattacharjee (2016), the Lundquist number is directly estimated using the size of the perturbation, which is typically the box size. In this case simulations try to reproduce the whole current layer in an SP configuration in regime IV of Figure 1. We, on the other hand, try to simulate a zoom-in of the middle of the current layer in regime III which initially looks like a planar current sheet but at later times will develop large-scale perturbations. If we define the Lundquist number using the box size, as $S={v}_{{Ay}}L/\eta$, and the system size L_S is actually bigger than the zoom-in box size, the Lundquist number of the whole system ${S}_{S}={v}_{{Ay}}{L}_{S}/\eta$ is larger than S which is quoted in Table 1. We can, therefore, safely assume that the larger system is unstable to tearing, just as our planar current sheet is unstable. A subtle difference of these two types of setup is that the simulations with global large-scale initial perturbations aim to describe the stationary regime IV at later times, $t\gg {L}_{S}/{V}_{{\rm{A}}}$, with a finite, albeit large, S, while our planar current sheet setup aims to simulate earlier times, $t\lt {L}_{S}/{V}_{{\rm{A}}}$, when the global outflow did not have time to form yet. We also assume that ${L}_{S}\gg L$ i.e., the global Lundquist number S_S is asymptotically large, so we can ignore gradients from the large-scale setup of the system and only impose small-scale perturbations. The end time of our simulations is determined by the development of large-scale structures, as long as these structures become comparable with the box size; the box size deemed not sufficient, similarly to simulations of the nonlinear Kelvin–Helmholtz instability.

Zoom In Zoom Out Reset image size

Table 1.  Spontaneous Reconnection MHD Experiments

Run	N3	Dissipation	S or ${S}_{4}^{0.4}$	${B}_{y0}/{B}_{z0}$	${v}_{r}/{v}_{{\rm{A}}}$
N1	5763	$-3.6\cdot {10}^{-4}{k}^{2}$	$1.7\times {10}^{4}$	1.0	0.0124
H1B1	5763	$-2.4\cdot {10}^{-9}{k}^{4}$	$2.5\times {10}^{4}$	0.5	0.0214
H1B2	5763	$-2.4\cdot {10}^{-9}{k}^{4}$	$2.5\times {10}^{4}$	1.0	0.0210
H1B3	5763	$-2.4\cdot {10}^{-9}{k}^{4}$	$2.5\times {10}^{4}$	2.0	0.0187
N2	7683	$-2.5\cdot {10}^{-4}{k}^{2}$	$2.5\times {10}^{4}$	1.0	0.0117
H2	7683	$-9.4\cdot {10}^{-10}{k}^{4}$	$3.7\times {10}^{4}$	1.0	0.0183
N3	11523	$-1.4\cdot {10}^{-4}{k}^{2}$	$4.4\times {10}^{4}$	1.0	0.0155
H3	11523	$-9.7\cdot {10}^{-10}{k}^{4}$	$3.6\times {10}^{4}$	1.0	0.0146
N4	15363	$-9.8\cdot {10}^{-5}{k}^{2}$	$6.4\times {10}^{4}$	1.0	0.0154
H4	15363	$-3.7\cdot {10}^{-10}{k}^{4}$	$5.4\times {10}^{4}$	1.0	0.0144

Download table as:  ASCIITypeset image

We used a pseudospectral code with explicit dissipation coefficients—viscosity and magnetic diffusivity equal to each other. The code solves incompressible resistive-viscous MHD equations, does not have inherent dissipation or dispersion grid errors, and is mostly described in Beresnyak (2014), except in the present paper we used full MHD, not reduced MHD. These simulations are DNS, i.e., they are well-resolved with the dimensionless maximum wavenumber satisfying ${k}_{\max }{l}_{\eta }\gt 1$, where ${l}_{\eta }={({\eta }^{3}/\epsilon )}^{1/4}$ is the dissipation (Kolmogorov) scale.¹

Normal viscosity and magnetic diffusivity $\eta {{\rm{\nabla }}}^{2}B$ was used in one series of simulations and hyper-diffusivity of the forth order, ${\eta }_{4}{{\rm{\nabla }}}^{4}B$ in the other. The Lundquist number was defined in terms of the box size $L=2\pi$ and the reconnecting field ${v}_{{Ay}}=1$ as $S={v}_{{Ay}}L/\eta$ and the hyper-Lundquist number as ${S}_{4}={v}_{{Ay}}{L}^{3}/{\eta }_{4}$. The Lundquist number S and the hyper-Lundquist number S₄ that give the same dynamical range of scales between the dissipation scale and outer scale, e.g., the ratio $L/{l}_{\eta }$, are related by $S={S}_{4}^{4/10}$, assuming Kolmogorov scaling. The rationale behind the two dissipation schemes was to check the influence not only on S, but also on the dissipation functional form to the bulk quantities, such as reconnection and dissipation rate. The magnetic Prandtl number was unity: ${{\Pr }}_{m}=\nu /\eta =1$, which was motivated by the desire to reach the highest possible S while staying well-resolved. As we will show below, the amount of magnetic and kinetic dissipation were close to each other, which is typical for an ordinary MHD cascade.

The simulations were set up with a thin current sheet with Harris profile and seeded with small initial perturbations, $\sim {10}^{-6}$ of the magnetic energy. These perturbations subsequently evolved due to oblique tearing near the current sheet itself, the bulk of the volume was almost undisturbed. At later times this evolved into a fully nonlinear turbulent current layer. The boundary between the undisturbed volume, which had nearly zero current, and the turbulent current layer was fairly well-defined, as shown in Figure 2. The inside of the current layer was determined as opposite to the undisturbed fluid where current density is always small. The point at which the current density exceeds a certain threshold in magnitude was regarded as the beginning of the current layer. We varied current density thresholds to study the dependence on the measurement of the layer width Δ. The procedure to obtain the errors for that measurement was to vary the lower and upper thresholds for current by a factor of two. The list of all simulations is presented in Table 1.

The evolution of the current layer width and the inferred reconnection rate are shown in Figure 3. The system initially contained the energy density of the opposing field ${B}_{y0}^{2}/8\pi$, which was free to dissipate, and the energy density of the mean imposed field ${B}_{z0}^{2}/8\pi$, which had to be conserved due to conservation of total flux through the x–y plane. After subtracting the latter contribution, we designate a dimensionless free energy density as

$$\begin{eqnarray}&&w=(4\pi \rho {v}^{2}+{B}^{2}-{B}_{z0}^{2})/{B}_{y0}^{2},\end{eqnarray} \tag{ 1 }$$

which is unity in the undisturbed fluid. After $t\approx 2$ turbulence in the layer fully develops, and the average w within the layer ${w}_{t}\approx 0.6$, while the undisturbed part still has w = 1. The dissipation of ${w}_{d}\equiv 1-{w}_{t}=0.4$ fraction of energy happens during the development of turbulence and stays approximately constant.

Zoom In Zoom Out Reset image size

We inferred the reconnection rate to be the growth speed of the current layer width ${v}_{r}=d{\rm{\Delta }}/{dt}$. This is different from the conventional definition which is an inflow speed in stage IV (Figure 1). In stage III, however, it is a meaningful definition in terms of how much free magnetic energy is available to the system per unit time per unit area of the current layer. From this energetic viewpoint, the inflow definition and our definition are similar.

V_r was around $0.015{v}_{{Ay}}$ for high Lundquist numbers and appears to be only weakly dependent on the imposed mean field B_z0 (Figure 3). The dissipation rate per unit area from both sides of the current sheet (note a factor of two) can be calculated from w_d and v_r as

$$\begin{eqnarray}&&{\epsilon }_{S}=2{w}_{d}{v}_{r}(1/2)\rho {v}_{{Ay}}^{2}\approx 0.006\rho {v}_{{Ay}}^{3},\end{eqnarray} \tag{ 2 }$$

where we note that conventional dissipation rate per unit mass, traditionally used in the theory of incompressible turbulence, will be expressed using current layer width Δ as

$$\begin{eqnarray}&&\epsilon =(1/\rho ){\epsilon }_{S}/{\rm{\Delta }}={w}_{d}{v}_{r}{v}_{{Ay}}^{2}/{\rm{\Delta }},\end{eqnarray} \tag{ 3 }$$

and will depend on time. The turbulence in the expanding current layer is not stationary turbulence in the sense that it grows in volumes and produces turbulent energy as well as dissipates energy. The outer scale of this turbulence also grows in time.

The contribution to turbulent fraction of the energy density ${w}_{t}\approx 0.6$ was partitioned to ∼0.55 in x and y magnetic components, ∼0.02–0.04 in the $\delta {B}_{z}$ component, and ∼0.01–0.02 in kinetic energy. The turbulence in the current layer was strongly anisotropic with respect to the B_z0 direction, with wavevector predominantly perpendicular to z. The B_x and B_y components, carrying most of the energy, therefore represented Alfvénic perturbations, while the sub-dominant $\delta {B}_{z}$ component was the slow-mode (pseudo-Alfvén) perturbation. We discuss anisotropy in more detail in Section 6, noticing that the reconnection rate depending only weakly on B_z0 is not surprising, since the anisotropic turbulence of Alfvénic perturbations, also known as Alfvénic turbulence or reduced MHD turbulence, possesses rescaling symmetry with respect to B_z0 (Beresnyak 2012), which we actually confirm in Section 6. The domination of Alfvénic perturbations in reconnection with strong mean field is extremely important as it sheds light on fluid-like behavior in plasma simulations, e.g., (Daughton et al. 2009), in spite of significant collisionless effects. The explanation for this is that reduced MHD is well-applicable to collisionless plasmas on scales above the ion Larmor radius r_L, and that plasma does not require significant collisional terms to behave like reduced MHD fluid (Schekochihin et al. 2009).

The structure of the perturbed current layer looked rather turbulent (Figure 4) with only a small fraction of the initial current sheet structure being retained.

Zoom In Zoom Out Reset image size

We define the power spectra of turbulent perturbations in the y–z plane as

$$\begin{eqnarray}&&E({k}_{l})={L}^{-1}\int \hat{f}({{\boldsymbol{k}}}_{l}){\hat{f}}^{* }({{\boldsymbol{k}}}_{l})d\phi \,{dx}\end{eqnarray} \tag{ 4 }$$

where ${k}_{l}=({k}_{y},{k}_{z})$—a wavevector in the y–z plane, $\hat{f}({{\boldsymbol{k}}}_{l})$—are Fourier transforms of either v or B. Neglecting k_x in this spectrum is necessary to get rid of the contribution from the $\pm {B}_{y0}$ jump across the current layer which happens in the x direction. The spectrum presented in Figure 5 has the magnetic contribution dominating over the kinetic one on large scales, but tends to approximate equipartition on smaller scales. This is not surprising, since turbulence is driven by magnetic energy. This spectral picture, qualitatively, is characteristic of decaying magnetic turbulence, including cases when the initial field was completely random (see, e.g., Biskamp 2003; Brandenburg et al. 2015).

Zoom In Zoom Out Reset image size

The total energy spectral slope was around $-1.5\div-1.7$, roughly consistent with Goldreich–Sridhar (Goldreich & Sridhar 1995) scaling. The slopes between −1 and −3 are indicative of local-in-scale turbulence. Precise measurement of the spectral slope in these simulations was difficult due to limited inertial range. However, the scale-locality would imply that at sufficiently high S the inertial range scaling and anisotropy will be the same as in the homogeneous driven MHD turbulence. In the next section we test the anisotropy component of this conjecture.

The scale-locality is a key ingredient in theories of turbulent reconnection. Indeed full scale-locality will imply that large-scale quantities, such as reconnection rate and dissipation rate per unit area, should be independent of any microphysics. Speaking in practical terms, if both ion Larmor radius r_L and ion skip depth d_i are much smaller than the minimum size of the problem, the layer width Δ, reconnection rate will be independent of microphysics. Note that in regime IV, stationary reconnection, Δ will become constant around $0.015L$.² This allows us to estimate the range of applicability of this particular mechanism of fast reconnection.

We  have calculated second-order structure functions of the turbulent v and B fields inside the current layer, e.g., for velocity:

$$\begin{eqnarray}&&{\mathrm{SF}}_{v}^{2}{({l}_{\parallel },{l}_{\perp })=\langle (v({\boldsymbol{r}}-{\boldsymbol{l}})-v({\boldsymbol{r}}))}^{2}{\rangle }_{{\boldsymbol{r}}}.\end{eqnarray} \tag{ 5 }$$

Note that we assume that it depends only on the component of ${\boldsymbol{l}}$ parallel and perpendicular to the magnetic field. Two types of such measurement are possible: when the parallel direction is determined by the global mean magnetic field, in our case the z direction, or local magnetic field ${\boldsymbol{B}}$. Scale-dependent anisotropy of the Goldreich & Sridhar (1995) model is observed with the local measurement (see, e.g., Cho & Vishniac 2000; Beresnyak & Lazarian 2009).

Using these structure functions for v and B we built correspondence between ${\lambda }_{\parallel }$ and ${\lambda }_{\perp }$ by equating SF values in parallel and perpendicular directions. More details on this type of measurement can be found in Beresnyak & Lazarian (2009). Figure 6 shows anisotropy ${\lambda }_{\parallel }/{\lambda }_{\perp }$ as a function of ${\lambda }_{\perp }$. One thing to notice is that the value of anisotropy in this particular simulation, with ${B}_{y}/{B}_{z}=1$, is around 20. At the same time the rms value of velocity perturbation is around $0.08\mbox{--}0.11{v}_{{Ay}}$. The interaction strength parameter $\xi =\delta v{\lambda }_{\parallel }/{v}_{{\rm{A}}}{\lambda }_{\perp }$ will be around unity, i.e., these perturbations are, approximately, "critically balanced." This means that from an MHD perspective we are dealing with "strong" turbulence, i.e., nonlinear interaction terms have the same contribution as the tension of the mean field B_z. Also note that this anisotropy corresponds to the angle of the field line bending ∼1/20 which is much smaller than the angle of the initial stripes of developing oblique tearing (45^o in the case of ${B}_{y}/{B}_{z}=1$). So the turbulence self-organizes itself into being strong and forgets properties of the oblique tearing that initiated it.

Zoom In Zoom Out Reset image size

The evidence for scale-dependent anisotropy is only tentative. Considering rather short inertial range, the expected law of scale-dependency from Goldreich & Sridhar (1995) is ${\lambda }_{\parallel }/{\lambda }_{\perp }\sim {\lambda }_{\perp }^{-1/3}$, see Figure 6.

Another important indicator is how anisotropy varies with the value of the mean field B_z0, this is presented on Figure 7. We plotted both local and global anisotropy measurements. Note how ${\lambda }_{\parallel }$ measurements have the same shape and are increasing with increasing B_z0. In purely Alfvénic dynamics, also called reduced MHD, ${\lambda }_{\parallel }$ is strictly proportional to B_z0 (see, e.g., Beresnyak 2012). We noticed that this scaling is almost perfect between ${B}_{z0}=1$ and ${B}_{z0}=2$ cases, which further confirms that Alfvénic dynamics dominates in the current layer turbulence.

Zoom In Zoom Out Reset image size

In the ${B}_{y0}/{B}_{z0}=1$ case, apart from Alfvén mode, which contains 93% of the total energy of turbulent motions, the other 7% are perturbations in B_z. These perturbations are strongly anisotropic and this component represents the pseudo-Alfvén or slow mode. The perturbations are energetically dominated by large scales and have a well-defined global structure: namely, the perturbation is negative (decreasing B_z) on the edge of the layer and positive (increasing B_z) in the middle of the layer, see Figure 8. Note that total B_z flux must be conserved. Why does turbulence create large-scale structure in B_z so that it is larger in the center? This could be due to the diamagnetism of turbulence, which is stronger on the edge, where turbulence is more intense, so that the diamagnetism of turbulence pushes B_z flux toward the center. This conjecture will require future research. The peculiar structure of the mean flux through the layer can have consequences for particle acceleration in turbulent current layers, as particles are more likely to be trapped in the low-B_z regions. When the mean field increases, the effect becomes negligible due to weaker coupling between the Alfvén and slow mode.

Zoom In Zoom Out Reset image size

The results reported above appear to be qualitatively different from previous 2D results, which was to be expected: the geometrical constraints in the 2D magnetic configuration naturally feature magnetic separatrixes, X-points, and magnetic islands, which are normally absent in 3D. Also, the dynamic influence of the global mean field, which is present in a generic reconnection geometry, is completely ignored by the 2D treatment. Another important point is that if B_z tension plays no important role in 2D, the perturbations that govern the 2D case are not Alfvénic and the arguments for the reduced MHD analogy we used in the above section are also not applicable.

The 3D spontaneous reconnection that we studied here proceeded in a different way than the 2D case, which, in Loureiro et al. (2012) was dominated by the ejection of plasmoids and had significant time-dependence. In the 3D case considered here, we observed a very steady rate (Figure 3) with the turbulent current layer slowly eating through the mostly undisturbed fluid and turbulence being fueled by the free energy of the oppositely directed magnetic fields. The 3D case was also different from the 2D case in that the memory of the initial conditions, i.e., the location of the original current sheet, was largely forgotten. In the 2D case the current sheet remains precisely where it was, up to very high S, generating and ejecting plasmoids along the same line. In the 3D case only small pieces of the original current sheet are visible after t = 10 and the layer otherwise looks turbulent. Few large-scale structures in 3D may be called flux ropes, however, unlike 2D, they are turbulent inside (Figure 4). In addition, the number of these structures does not depend on S as it does in the 2D case described in Uzdensky et al. (2010). Another difference with 2D is the asymmetry of emerging turbulence with respect to the original current sheet—in 3D often the upper or a lower part of the layer dominates.

The classic X-point inflow/outflow picture is usually preserved in 2D in each X-point between plasmoids. In our simulations such a simple picture was not observed, see Figure 9.

Zoom In Zoom Out Reset image size

The physical reason for the resistivity-independent rate also appears to be different: in the 2D case it relies on the hierarchical formation and ejection of plasmoids (Uzdensky et al. 2010), while in 3D we hypothesize this to be a consequence of turbulence locality, similar to models of reconnection due to ambient turbulence (Lazarian & Vishniac 1999). The difference between 3D reconnection with ambient turbulence and spontaneous reconnection is that in the spontaneous case there is no external agent that drives reconnection and there is no parameter of the amplitude of ambient turbulence as in, e.g., Kowal et al. (2009).

Resistively independent turbulent reconnection has been argued to be the case of reconnection due to ambient turbulence (Lazarian & Vishniac 1999; Eyink et al. 2011). In this paper we extend this result to the case where reconnection develops spontaneously without an external agent. We demonstrate this by showing that turbulence in the current layer resembles ordinary turbulence and is likely local in scale. At sufficiently large Lundquist numbers, the dynamics of large scales, which determines such properties as bulk reconnection and dissipation rates, will be disconnected from the dynamics on plasma scales and dissipation parameters. It is therefore natural to expect reconnection and dissipation rates to go to asymptotic universal values. Whether the locality argument can be applied in 2D is unclear because large plasmoids may couple large and small scales directly.

The reconnection speed in ambient turbulence was argued to be proportional to the kinetic energy density (Lazarian & Vishniac 1999). Our ${v}_{r}=0.015{v}_{{Ay}}$ measurement obtained in the absence of ambient turbulence could be seen as a lower limit on fluid reconnection speed in 3D.

2D theory explaining observed reconnection rate is based on a hierarchical plasmoid model (Uzdensky et al. 2010). This model predicts that the number of plasmoids scales linearly with S and this hierarchy truncates on the scale of the critical layer. Thus the reconnection rate is just equal to the SP rate at the critical value of S: ${v}_{r}={v}_{{\rm{A}}}{S}_{\mathrm{crit}}^{-1/2}$. In our picture reconnection is not due to large-scale tearing, but due to turbulence on the edge of the layer. This turbulence, as we showed above is strong and cannot be considered a linear stage of any instability.

One of the interesting empirical facts about spontaneous reconnection that we observed in simulations is that the reconnection rate is constant in time and the level of velocity and magnetic perturbations remains approximately on the same level as well. Can this be reconciled with the turbulent picture, despite the volumetric dissipation rate inside the layer depending on time? The basic scaling for turbulent velocity as a function of turbulent cascade rate is

$$\begin{eqnarray}&&\delta {v}_{l}^{2}={C}_{K,v}{\epsilon }^{2/3}{l}^{2/3},\end{eqnarray} \tag{ 6 }$$

where l is a scale of interest and we introduced the Kolmogorov constant ${C}_{K,v}$ that refers to the velocity perturbation, not the total energy spectrum. If we argue that turbulence is driven on the scale of the current layer thickness, i.e., $l={\rm{\Delta }}$, and use Equation (3) for the dissipation rate, we calculate that $\delta {v}_{l}^{2}$ indeed does not depend on Δ and, therefore, on time:

$$\begin{eqnarray}&&\delta {v}_{l}^{2}={C}_{K,v}{w}_{d}^{2/3}{v}_{r}^{2/3}{v}_{{Ay}}^{4/3}.\end{eqnarray} \tag{ 7 }$$

The fact that $\delta {v}_{l}$ is constant in time is consistent with our numerical measurement. How reconnection rate depends on $\delta {v}_{l}$ is not immediately obvious. One possibility is to use the expression for the turbulent reconnection rate from Lazarian & Vishniac (1999), but it is not clear how the layer width relates to the time-dependent injection scale l. Assuming that v_r depends only on the rms velocity $\delta {v}_{l}$  in the current layer, however, in the manner ${v}_{r}={v}_{{Ay}}f(\delta v/{v}_{{Ay}})$ we can obtain the time-independent rate by substituting Equation (7) and solving the nonlinear equation

$$\begin{eqnarray}&&{v}_{r}={v}_{{Ay}}f({C}_{K,v}^{1/2}{w}_{d}^{1/3}{v}_{r}^{1/3}{v}_{{Ay}}^{-1/3}).\end{eqnarray} \tag{ 8 }$$

In particular, using the ${v}_{r}\sim {M}_{{\rm{A}}}^{2}$ dependence from Lazarian & Vishniac (1999), e.g., choosing ${v}_{r}={C}_{\mathrm{LV}}\delta {v}_{l}^{2}/{v}_{{Ay}}$ we obtain

$$\begin{eqnarray}&&{v}_{r}={C}_{\mathrm{LV}}^{3}{C}_{K,v}^{3}{w}_{d}^{2}{v}_{{Ay}}.\end{eqnarray} \tag{ 9 }$$

The constant C_LV has not yet been precisely measured (c.f. Kowal et al. 2009). The constant ${C}_{K,v}$ can be obtained in our simulations and refers to the ordinary Kolmogorov constant, as well as the fraction of total cascade energy that resides in its kinetic part. Introducing the ratio of kinetic to magnetic energy as r_A, which is around 0.13 in our simulations, we can estimate ${C}_{K,v}={C}_{K}{r}_{{\rm{A}}}/(1+{r}_{{\rm{A}}})\approx 0.48$ using C_K = 4.2 from Beresnyak (2011). This gives the reconnection rate of $0.018{C}_{\mathrm{LV}}^{3}{v}_{{Ay}}$, which corresponds to the measured value if ${C}_{\mathrm{LV}}=0.94$.

Interestingly, this expression depends only on the basic properties of well-developed turbulence, dimensionless numbers C_K and r_A, and not on the properties of the instability.

The actual dissipation mechanisms of reconnection which will result in observable phenomena are still debated. It is plausible that dissipation in plasmas sometimes results in heating, and sometimes in the acceleration of fast particles. For example, the acceleration on shock fronts starts with particles being pulled out of the thermal pool due to an extremely high velocity gradient at the shock front itself. Similarly, solar X-ray flares, which produce accelerated electrons, have been brought up as proof that current layers must have microscopic widths to allow for plasma effects, including parallel electric field, and electron acceleration. The logic is the following. Suppose that the turbulent reconnection picture is true and current layers are wide compared with plasma scales and electrons only "feel" local turbulent perturbations, in which case electron acceleration will be, basically, stochastic turbulent acceleration. Turbulent acceleration, specifically from quasilinear theory (Schlickeiser 2002), was calculated to be second order in v/c, too slow in many practical cases, while acceleration by the electric field $E\sim {{\boldsymbol{v}}}_{r}\times {\boldsymbol{B}}/c$ is first order in v/c and should dominate.

We recently found that the very basis of the above argument, the claim that turbulent acceleration must be second order, is in fact not true. In particular, we found analytically that if turbulence is fueled by magnetic energy, like in the case of spontaneous turbulent reconnection, the structure of magnetic and electric fields in this turbulence is such that the average acceleration by curvature drift is positive. This is due to a mathematical relation between the MHD term which is responsible for energy transfer between kinetic and magnetic energy and the term responsible for the curvature drift acceleration (Beresnyak & Li 2016). This does not require extra assumptions such as that particles have to be trapped for a considerable time in magnetic islands, in fact the whole volume of turbulence will be the first order accelerator. So, as long as the particle gyro-radius is smaller than the current layer width, the acceleration of these particles will be efficient. The expression for the acceleration rate we derived in Beresnyak & Li (2016)

$$\begin{eqnarray}&&\displaystyle \frac{d{ \mathcal E }}{{dt}}={{ \mathcal E }}_{\parallel }\displaystyle \frac{8\pi }{{B}^{2}}{ \mathcal D },\end{eqnarray} \tag{ 10 }$$

relates it to the energy transfer from magnetic to kinetic energy ${ \mathcal D }$ and the acceleration rate. In ordinary driven turbulence this term is zero, while in spontaneous reconnection it is equal to half of the total dissipation rate . Since this mechanism results in average acceleration for all particles, all electrons are predicted to be accelerated to approximately the same energy, $0.35T{(L/{d}_{i})}^{1/4}$ (Beresnyak & Li 2016), where T is the thermal energy. The subsequent transition to regime IV will result in an outflow, particle escape, and additional acceleration due to converging magnetic field lines, which will result in the formation of a power-law tail. It is interesting that X-ray emission during solar flares indeed features a thermal component and a power-law tail.

Our result suggests that in the high-S limit all macroscopic properties of reconnection are expressed in terms of macroscopic plasma properties ρ, v_A, and L, and independent of microscopic dissipation. This result is quite spectacular considering that individual field line reconnection does depend on microphysics.

Let us think of the turbulent current layer as a conductor. The electromotive force between two points separated by a large distance in the z direction will be expressed as

$$\begin{eqnarray}&&{ \mathcal E }=\int \displaystyle \frac{{\boldsymbol{v}}\times {\boldsymbol{B}}}{c}\,{dz}\approx \displaystyle \frac{{v}_{r}}{c}{B}_{y0}{L}_{z},\end{eqnarray} \tag{ 11 }$$

and should be independent of whether we take a point inside the current sheet or outside of it. At the same time, the total current flowing through the width L_y of the current layer will be $I=(c/4\pi )2{B}_{y0}{L}_{y}$. Taking the ratio of the two, we obtain the effective resistance of the current layer

$$\begin{eqnarray}&&R=\displaystyle \frac{1}{2}\displaystyle \frac{{v}_{r}}{c}\displaystyle \frac{{L}_{z}}{{L}_{y}}\left(\displaystyle \frac{4\pi }{c}\right)\approx 0.0075\displaystyle \frac{{v}_{{\rm{A}}}}{c}\displaystyle \frac{{L}_{z}}{{L}_{y}}\left(\displaystyle \frac{4\pi }{c}\right).\end{eqnarray} \tag{ 12 }$$

Note that ${R}_{0}=4\pi /c$ ($\sim 376.73{\rm{\Omega }}$ in SI units) is known as the impedance of free space, and our final result is independent of microscopic resistivity. Such a resistance will dissipate energy in conductive fluids in spite of the fact that microscopic resistivity of plasma in most astrophysical objects can be considered negligible. In the limit of infinitely heavy plasma, ${v}_{{\rm{A}}}/c\to 0$, this will result in zero resistance, consistent with the ordinary resistance expression proportional to resistivity. For a very light plasma, i.e., the relativistic force-free magnetically dominated limit,  ${v}_{{\rm{A}}}/c=1$ and the resistance is a sizable fraction of the impedance of free space.

For example, jets in active galactic nuclei are self-contained electromagnetic structures carrying large-scale poloidal current, with at least part of the return poloidal current flowing in a layer separating the magnetic pressure-dominated jet and the outside medium (Begelman et al. 1984). A Poynting-dominated jet has an impedance of ∼90Ω (Lovelace & Kronberg 2013), and since ${v}_{{\rm{A}}}/c\sim 1$ for rarefied electron–positron plasma, we can estimate that jets with aspect ratios of ${L}_{z}/{L}_{y}\gt 650$ will dissipate a sizable fraction of their energy in an outer current layer, due to this layer's fundamental fluid resistance. Whether this will result in an outer layer's visibility is an open question, but considering the result of the previous section, the first order acceleration of particles and non-thermal emission from this layer is highly likely.

Another example is dissipation in pulsar magnetospheres, which feature the return current layer in the equatorial plane beyond the light cylinder (Uzdensky & Spitkovsky 2012). The current layer separating open and closed field lines within the light cylinder (see, e.g., Arons 2011) can also result in acceleration.

Our simulations clearly demonstrate that turbulence must be a part of high-Lundquist reconnection, e.g., astrophysical reconnection. Quite importantly, this turbulence is not random, but contains non-trivial correlations that come from the fact that energy is transferred from magnetic to kinetic, these correlations likely to result in efficient particle acceleration (Beresnyak & Li 2016), which can help explain why reconnection on the Sun results in powerful X-ray flares. Observational evidence favoring spontaneous turbulent reconnection includes magnetospheric observations that show an enhanced level of turbulence inside current sheets (Matsumoto et al. 2003; Cattell et al. 2005).

Many astrophysical objects, such as the interstellar medium in our Galaxy, feature relatively high levels of ambient turbulence and the reconnection is argued to be fast (Lazarian & Vishniac 1999) due to existing magnetic field stochasticity. In highly magnetized environments, such the solar surface or pulsar wind nebulae, the velocity of ambient turbulence may be tiny compared to the local Alfvén speed, and so turbulence spontaneously generated by the current sheet and fueled by the reconnecting field itself is more important. More qualitatively, if ${M}_{{\rm{A}}}=\delta v/{v}_{{\rm{A}}}\lt \sqrt{0.015}\approx 0.12$, spontaneous reconnection will dominate. Future parameter studies in simulations with the Lundquist number above the critical and varying level of ambient turbulence should clarify the transition between turbulent reconnection due to ambient turbulence and spontaneous reconnection. So far, simulations with driven turbulence have $S\lt {10}^{4}$, so the non-driven case is still consistent with the SP picture (Kowal et al. 2009), with resistive rate higher than 0.015 that we measured in this paper.

Currently, two completely opposite ways to explain fast spontaneous reconnection exist. First is the model that relates reconnection rate to the critical Lundquist number, i.e., the claim that reconnection rate depends on stability to linear resistive tearing. This gives a dimensionless rate of ${S}_{\mathrm{crit}}^{-1/2}$ (Uzdensky et al. 2010). The second relates reconnection to the inherent properties of strong turbulence and gives the reconnection rate of ${C}_{K}^{3}{r}_{{\rm{A}}}^{3}/{(1+{r}_{{\rm{A}}})}^{3}$ (this paper). Both pictures reasonably agree with the measurement, but are hard to reconcile with each other. One way to connect these pictures is to imagine that turbulence produces constant anomalous resistivity which brings the effective Lundquist number to at or below the critical value, just to barely suppress tearing. The counterargument to this is that linear stability studies have only been performed in the laminar SP regime and the large value of S_crit is the result of a rather non-trivial interplay between resistive tearing rate $\sim {\eta }^{1/2}$ and thinning of the current layer. In a non-laminar case it is easy to argue for ${S}_{\mathrm{crit}}\sim 1$, but hard to argue for ${S}_{\mathrm{crit}}\sim {10}^{4}$.

In stationary reconnection with speed v_r that develops an outflow with speed ${v}_{\mathrm{out}}$ and reaches the quasi-stationary state as in Figure 1 panel (IV) the inflow speed is balanced by the outflow ${v}_{\mathrm{in}}={v}_{r}$. In this case the current layer width also reaches the stationary value of ${\rm{\Delta }}={{Lv}}_{r}/{v}_{\mathrm{out}}$, neglecting compressibility. Similar to the cascade locality argument, e.g., the argument why v_r should be independent of S, the outflow speed also must be a fraction of v_A. It is not necessarily equal to v_A as is often assumed. The outflow develops under the force of magnetic tension from the B_y component but we showed that a sizable fraction of this energy is dissipated and not converted into kinetic motion. We expect the outflow speed to be a sizable fraction of v_Ay, this fraction being lower by a factor of $\sim \sqrt{1-0.4}\approx 0.77$, where ${w}_{d}\approx 0.4$ is a dissipation factor that we found in this paper.

The reconnection rate could be affected by the presence of an outflow in regime IV. Increasing the box size will allow larger-scale fluid motions which could emulate the outflow effect locally; this also corresponds to higher S. We showed that with increasing S the reconnection rate goes to a constant. This is probably associated with the fact that most of the activity which results in growth of the current layer happens on the boundary. Likewise, the simulation naturally included the effects of the local outflows that develop on larger and larger scales as reconnection progresses, while the reconnection rate stays relatively stable (Figure 3).

It is also worth mentioning that the fraction of dissipated energy that we measured, ${w}_{d}\approx 0.4$, will limit the compression ratio of the current layer. Previously it was thought that most of the magnetic energy during reconnection could be spent to accelerate the outflow jet, which is why the outflow speed was always estimated as being equal to the Alfvén speed (see above). In this case the compression ratio may be arbitrarily high for low beta plasmas, going as $1/\beta$. In our case the compression ratio will be limited to $1/{w}_{t}(\gamma -1)\approx 3.8$ for monoatomic gas.

I am grateful to Alex Schekochihin, Nuno Loureiro, Fan Guo, Alex Lazarian, Ethan Vishniac, Homa Karimabadi, Bill Daughton, and Hui Li for fruitful discussions. Computations were performed on NICS Kraken with XSEDE allocation TG-AST110057 and on LANL institutional computing resources.