# Self-generated turbulence in magnetic reconnection

## Abstract

Classical Sweet-Parker models of reconnection predict that reconnection rates depend inversely on the resistivity, usually parameterized using the dimensionless Lundquist number (). We describe magnetohydrodynamic (MHD) simulations using a static, nested grid that show the development of a three-dimensional instability in the plane of a current sheet between reversing field lines without a guide field. The instability leads to rapid reconnection of magnetic field lines at a rate independent of over at least the range resolved by the simulations. We find that this instability occurs even for cases with that in our models appear stable to the recently described, two-dimensional, plasmoid instability. Our results suggest that three-dimensional, MHD processes alone produce fast (resistivity independent) reconnection without recourse to kinetic effects or external turbulence. The unstable reconnection layers provide a self-consistent environment in which the extensively studied turbulent reconnection process can occur.

Department of Physics, Farmingdale State College, Farmingdale, NY 11735 \altaffiliationalso Department of Astrophysics, American Museum of Natural History, New York, NY \affiliationDepartment of Astrophysics, American Museum of Natural History, New York, NY \affiliationDepartment of Physics, Florida State University, Tallahassee, FL \affiliationDepartment of Physics, Barnard College, New York, NY

## 1. Introduction

During magnetic reconnection, magnetic field lines change topology, resulting in the conversion of magnetic energy into both thermal energy and kinetic energy of bulk flows and non-thermal particles. The rate at which this process occurs in the classical Sweet-Parker picture (Sweet, 1958; Parker, 1957) depends on the Lundquist or magnetic Reynolds number , where is the Alfvén speed, a characteristic length of the system, and the resistivity. The Sweet-Parker rate is orders of magnitude too slow to explain the fast reconnection seen in low resistivity plasmas during solar flares and sawtooth crashes in tokamaks (Yamada et al., 2010). Because it is a fundamental plasma process, reconnection is thought to be important in astrophysical environments as diverse as the heliosphere (e.g. Edmondson et al., 2010) and microquasars (Khiali et al., 2015).

The identification of the 2D plasmoid instability 1(Biskamp, 1986; Loureiro et al., 2007; Huang & Bhattacharjee, 2013), a super-Alfvénic, small-scale instability, has provided a mechanism to greatly speed up Sweet-Parker reconnection. However, this instability has primarily been studied in two dimensions (2D), assuming symmetry in the plane of the current sheet. The reason for this dimensional reduction is that the reconnection process is inherently multi-scale, with a large separation between the global scale of the reconnection layer and the resistive length where the instability grows. Even 2D simulations tax state of the art computational resources if uniform grids are used.

Lazarian & Vishniac (1999) argued that reconnection in the presence of any sort of turbulence would be fast, because the turbulence would drive many points of contact between the opposed field lines. This idea has been put on a rigorous mathematical basis (Eyink et al., 2011) as reviewed by Lazarian et al. (2015a) and Lazarian et al. (2015). Indeed, recent modelling suggests that turbulent reconnection may be responsible for the radio and gamma-ray emission from accreting black holes (Singh et al., 2015). However, these ideas all require that reconnection proceed at a large fraction of . Numerical models examining reconnection in forced turbulence support this theory, starting with (Kowal et al., 2009). In this work, we demonstrate that turbulent reconnection proceeds in a very similar fashion when the turbulence is self-generated from an instability of the reconnection layer itself.

Here, we describe a set of nested grid simulations that model the reconnection layer in 3D over a broad range of , without any forcing or guide field. These simulations show that a startlingly fast, 3D, instability occurs in the plane of the current sheet, which was assumed uniform in the 2D simulations. This instability drives a large increase in the rate of reconnection, that we show remains independent of over two orders of magnitude of variation in the resistivity.

Boozer (2012b, 2013) has argued that reconnection requires a point geometry to proceed, so that it is an inherently three-dimensional (3D) process. The work we describe here demonstrates that such a 3D geometry naturally arises even from 2D initial conditions, resulting in fast reconnection apparently independent of .

Previous work in this field has shown 3D instability, but has not provided a clear demonstration of independence of reconnection rate from . Dahlburg et al. (2003, 2005) focused on the case of a current sheet with a strong guide field, and found a 3D instability set in for a weak enough guide field, which they called a secondary instability. However, they did not measure the scaling of the reconnection rate with . Lapenta & Bettarini (2011) reported the breakdown of an initially 2D Harris sheet into a fully 3D reconnection region with greatly enhanced reconnection rate. An MHD kink instability on a central plasmoid was followed by a Rayleigh-Taylor instability driven by the reconnection jet interacting with the plasmoids at the ends of the layer. However, again, no test of the dependence on was performed. Another numerical experiment has shown that thin, 3D, current layers are unstable to infinitesimal perturbations and reconnect at a rate apparently independent of Lundquist number (Beresnyak, 2013), but only a factor of three variation in was explored. Edmondson et al. (2010) studied the formation of coronal current sheets due to photospheric forcing in a global, 3D, AMR simulation. They concluded that the dynamics of the current sheet were 3D, allowing a steady rather than the bursty reconnection rate found by 2D models of the plasmoid instability. Finally, 3D reconnection in the collisionless limit has been explored by Daughton et al. (2011) and Pritchett (2013). That work largely focused on particular kinetic effects that drive dissipation at the smallest scales.

## 2. Methods

We use the mesh refinement code Enzo, which solves the compressible, adiabatic, resistive, MHD equations (Bryan et al., 2014). We use static refinement to focus computational effort on the current sheets. Our computations are performed within a cubic, 3D volume, with periodic boundary conditions in all three dimensions. From the available algorithmic options, we choose piecewise linear reconstruction, the HLLD Riemann shock-capturing solver, and constrained transport (Gardiner & Stone, 2005) to ensure to machine precision (Collins et al., 2010). We performed all analysis using the yt toolkit (Turk et al., 2011).

Our initial condition is a pair of oppositely directed, parallel current sheets to accommodate the periodic boundary conditions, each perturbed following the GEM Reconnection Challenge (Birn et al., 2001) to initiate Sweet-Parker reconnection. All but two of our runs are initialized with low-amplitude, 3D velocity perturbations with mean Alfvèn Mach number . These perturbations have a spectrum with wavenumbers ranging from . We choose and , except for run C+, which has and . We do not continue to force the velocity field during the simulation. We normalize all lengths to the size of the box , densities to the initial density at the center of the sheets , and times to the Alfvén crossing time of each sheet where is the Alfvén speed of the upstream plasma. is the scale length of the initial current sheet. Table 1 lists the parameters of our runs.

A minimum resolution requirement for reconnection is the proper
resolution of the Sweet-Parker current layer, whose width , where is the length of the layer^{1}

Figure 1 shows of the initial Sweet-Parker current sheet at , long before any unstable perturbations have grown to significant amplitudes. All runs with have current sheet widths that agree well with the Sweet-Parker prediction, because they are resolved by zones across the sheets. The run with demonstrates the effects of marginal resolution, while the run with is only resolved by zones, and is a factor of four too thick. We do not use this last run (run J) in our subsequent analysis, although it serves as an important limit on the numerical resistivity of our code.

run | notes | |||

A | ||||

A* | double resolution | |||

B | ||||

C | ||||

C+ | perturbation | |||

D | ||||

E | ||||

F | – | stable to 3D instability | ||

G | – | no initial perturbations | ||

H | – | no initial perturbations | ||

J | – | underresolved, unanalysed |

Lundquist number Resistivity in code units Decay rate of magnetic energy (see text)

## 3. Field dynamics

Reconnection in our models begins at the Sweet-Parker rate expected for a stable field reversal, as shown by the width of the current sheet. This leads to the initial slow decline of the volume integrated magnetic energy for all simulations (Figure 2), as well as the low values of integrated kinetic energy. Instability along the plane of the current sheet then sets in, driving far faster reconnection, and transferring energy from the magnetic field into the flow, as shown by the sudden drop in magnetic energy and the corresponding rise in kinetic energy. The morphology of the onset and growth of the instability is shown in the middle panels of Figure 3, while the final panel shows its saturated state.

The
development of the 3D instability results in the buckling of the current sheet in the
– plane, with a characteristic wavenumber (third
panel of Fig. 3). The
simulations of Lapenta & Bettarini (2011) can be seen to show similar
behavior, though it is not emphasized in their paper. They used a
thin box in the third dimension, so they only had two wavelengths in
that direction. Ours is a factor of six deeper in the direction than
theirs^{2}

We find that resistivity stabilizes the 3D instability for . Run D with shows the instability clearly through the growth of kinetic energy, although the total amount of reconnection driven by the turbulence is small, because the large resistivity has already allowed significant laminar reconnection to occur.

When , the transition from laminar to turbulent reconnection begins with the rapid growth of kinetic energy apparently driven by a kink-type instability along the plasmoids in the -direction. However, this is not a classical kink instability, as the interior of the plasmoids is a demagnetized reconnection region, rather than a column of current-carrying plasma surrounded by vacuum as would be true in the classical case. The growth of kinetic energy and decay of magnetic energy must be due to reconnection occurring where field lines are driven together by these instabilities rather than a simple rearrangement of the horizontal field by them.

It appears from our models that the 3D instability may actually grow independently of the 2D plasmoid instability. Run E with lacks evidence for the growth of the 2D plasmoid instability seen at higher Lundquist numbers by ourselves and previous authors (Loureiro et al., 2007; Samtaney et al., 2009; Uzdensky et al., 2010; Huang & Bhattacharjee, 2010; Loureiro et al., 2012; Huang & Bhattacharjee, 2013), but nonetheless shows the growth of the 3D instability, albeit with a delay in onset of rapid growth (Fig. 2).

This delay occurs because growth of the 2D plasmoid instability triggers secondary Richtmyer-Meshkov instability (Richtmyer, 1960; Meshkov, 1969) that accelerates onset of the 3D instability, but is not required for 3D instability to occur. The acceleration of a flow across the density contrast between the plasmoid and the surrounding current sheet along the axis drives the Richtmyer-Meshkov instability (analogous to the Rayleigh-Taylor instability that occurs in a gravitational field). This instability begins when weak transient shocks from the formation of the plasmoids at the sides of the domain () pass over the density gradient produced by the first plasmoid to form around the initial X line at . The secondary instability drives weak, initial mixing along the X line.

We examined three numerical issues with further runs. First, to determine if resolution affects our major result, we run our marginally resolved run A at twice the resolution (run A*). This run results in essentially identical growth rate .

Second, we checked that the instability is not purely numerical by performing Runs G & H, identical to the unstable Runs E & B respectively, except without any initial perturbations. The instability did not grow in either of these runs. In Run G, which has and is thus stable to the 2D plasmoid instability, the entire 3D volume settled down into a steady, laminar reconnection at the Sweet-Parker rate, with a sheet width . Run H, on the other hand, is unperturbed in the – plane, but is unstable to the plasmoid instability. In this case, we find a vigorous plasmoid instability that is entirely symmetric along the axis, demonstrating that our code is sufficiently stable to recover the 2D results in 3D if there are no explicit 3D perturbations.

Third, run C+ was performed at standard resolution with smaller scale perturbations (). We find the growth rate is higher by a factor for the lower perturbations, suggesting a wavenumber dependence for the underlying instability. We will pursue a formal stability analysis of the instability in a separate paper, and this wavenumber dependence represents an important test for that work.

## 4. Scaling

Figures 2 and 3 show three phases of reconnection: the slow, Sweet-Parker phase while linear instabilities grow, a rapid exponential phase in which 3D effects dominate reconnection, and finally a saturated, MHD turbulent phase. The Sweet-Parker reconnection rate is , implying far slower than observed reconnection at the large values of typical of Solar and space plasmas. Figure 4 shows the decay rate during the rapid reconnection phase as a function of . While runs over two orders of magnitude, varies by a factor of only about three, with no discernible functional relationship to . Thus, the 3D instability offers a fast reconnection mechanism that occurs at a rate apparently independent of , without appeal to either kinetic effects or anomalous resistivity.

Once the instability saturates, the current sheet thickens considerably (see the right panel of Fig. 3), and the picture of a steady flow of fresh field from upstream (i.e. from the direction above and below the layer) no longer holds in our simulations. At the end of our simulation, there is still plenty of field left to reconnect. The turbulence self-consistently driven in the reconnection layer allows stochastic reconnection to occur (Lazarian & Vishniac, 1999; Eyink et al., 2011, 2013). The thickening of the reconnection layer we see is consistent with their model: the diffusion of the large-scale magnetic fields controls the ultimate reconnection rate. Our periodic boundary conditions, similar to those of Beresnyak (2013), do not allow for a self-consistent steady state to occur. However, during the rapid, resistivity-independent phase of the instability, the reconnection region grows.

## 5. Discussion

We have shown that for Lundquist numbers , current sheets become turbulent in the direction perpendicular to the field along the sheet, leading to fast, 3D, magnetic reconnection (i.e. decay rate of magnetic energy independent of resistivity ). At high , individual 2D plasmoids rapidly lose their identities, as the current sheet splits into filaments parallel to the field direction. These 3D effects, driven by rapidly growing instabilities along current sheets, appear essential to understanding reconnection. This provides strong support to the geometric ideas advanced by Boozer (2012b, a, 2013) as well as the turbulent reconnection model developed by Lazarian, Vishniac, and coworkers (Lazarian & Vishniac, 1999; Lazarian et al., 2015, 2015a). To demonstrate the latter point, in the bottom panel of Figure 2 we show the reconnection layer width as a function of time. During the rapid growth phase, , in pleasing agreement with the theory described in Lazarian et al. (2015a). We find that the layer expansion velocity is , roughly a factor of 4–5 smaller than that reported by Beresnyak (2013). We suspect the discrepancy is due to the fact that his simulations include a guide field, although the lower diffusivity of his pseudo-spectral code could also play a role. Nevertheless, the agreement on the linear form of the growth rate despite our different setups and codes supports the turbulent reconnection model.

As a result of the 3D instability, the initial current sheet develops into a thick region of MHD turbulence. Lazarian et al. (2015a) speculate that the growth in the plane perpendicular to the reconnection (i.e. along the direction in our simulations) could be due to Kelvin-Helmholz instability. Once the turbulent state is reached, the decay of magnetic energy in our model slows dramatically, back to a rate comparable to the initial Sweet-Parker rate. However, we stress that the state of the system after the rapidly reconnecting, unstable phase is radically different from the state before it, and the slow subsequent reconnection may depend on the geometry we chose for our simulations, which does not continue to force the system on large scales, unlike, for example Solar flares Dudík et al. (2014). However, this turbulence naturally produces the conditions required for fast stochastic reconnection (Lazarian & Vishniac, 1999; Eyink et al., 2011, 2013).

Future work employing deeper grid hierarchies, adaptive resolution elements placed in regions of high , and simulations with large scale forcing through inflow boundary conditions will clarify the outcomes of the instability and its role in understanding reconnection in astrophysical environments. Specifically, by removing periodic boundary conditions (or isolating them via deep AMR hierarchies), we will be able to make a more detailed study of a key prediction of the Lazarian & Vishniac (1999) model, the rate of broadening of reconnection layer.

## 6. Acknowledgments

JSO & M-MML acknowledge support from NSF grant AST10-09802. MT
participated in the REU program at AMNH, funded by NSF grant
AST10-04591. We thank the anonymous referees for comments that
significantly improved the paper. We also thank A. Boozer,
C. R. DeVore, and A. Hubbard for useful discussions. This work used
the Extreme Science and Engineering Discovery Environment (XSEDE),
which is supported by National Science Foundation grant number
OCI-1053575. The computations presented here were performed on
*Kraken* at the National Institute for Computational Science,
under grant TG-AST120045, and *Stampede* at the Texas Advanced
Computing Center under grants TG-MCA99S024 and TG-AST140008.

### Footnotes

- A
popular alternative is to use to represent the
*half*-width of the current sheet, but in that case, is the half-length as well. - Note that the Lapenta & Bettarini (2011) box is oriented so their axis corresponds to our axis

### References

- Beresnyak, A. 2013, ArXiv e-prints, arXiv:1301.7424
- Birn, J., Drake, J. F., Shay, M. A., et al. 2001, J. Geophys. Res., 106, 3715
- Biskamp, D. 1986, Physics of Fluids, 29, 1520
- Boozer, A. H. 2012a, Physics of Plasmas, 19, 092902
- —. 2012b, Physics of Plasmas, 19, 112901
- —. 2013, Physics of Plasmas, 20, 032903
- Bryan, G. L., Norman, M. L., O’Shea, B. W., et al. 2014, ApJS, 211, 19
- Collins, D. C., Xu, H., Norman, M. L., Li, H., & Li, S. 2010, The Astrophysical Journal Supplement Series, 186, 308
- Dahlburg, R. B., Klimchuk, J. A., & Antiochos, S. K. 2003, Advances in Space Research, 32, 1029
- —. 2005, ApJ, 622, 1191
- Daughton, W., Roytershteyn, V., Karimabadi, H., et al. 2011, Nat. Phys., 7, 539
- Dudík, J., Janvier, M., Aulanier, G., et al. 2014, ApJ, 784, 144
- Edmondson, J. K., Antiochos, S. K., DeVore, C. R., & Zurbuchen, T. H. 2010, ApJ, 718, 72
- Eyink, G., Vishniac, E., Lalescu, C., et al. 2013, Nature, 497, 466
- Eyink, G. L., Lazarian, A., & Vishniac, E. T. 2011, ApJ, 743, 51
- Gardiner, T. A. & Stone, J. M. 2005, J. Comput. Phys, 205, 509
- Huang, Y.-M., & Bhattacharjee, A. 2010, Phys. Plasmas, 17, 062104
- —. 2013, Phys. Plasmas, 20, 055702
- Khiali, B., de Gouveia Dal Pino, E. M., & del Valle, M. V. 2015, MNRAS, 449, 34
- Kowal, G., Lazarian, A., Vishniac, E. T., & Otmianowska-Mazur, K. 2009, ApJ, 700, 63
- Lapenta, G., & Bettarini, L. 2011, Europhys. Lett., 93, 65001
- Lazarian, A., Eyink, G. L., Vishniac, E. T., & Kowal, G. 2015a, in Astrophysics and Space Science Library, Vol. 407, Astrophysics and Space Science Library, ed. A. Lazarian, E. M. de Gouveia Dal Pino, & C. Melioli, 311
- Lazarian, A., Eyink, G. L., Vishniac, E. T., & Kowal, G. 2015b, ArXiv e-prints, arXiv:1502.01396
- Lazarian, A., & Vishniac, E. T. 1999, ApJ, 517, 700
- Loureiro, N. F., Samtaney, R., Schekochihin, A. A., & Uzdensky, D. A. 2012, Phys. Plasmas, 19, 042303
- Loureiro, N. F., Schekochihin, A. A., & Cowley, S. C. 2007, Phys. Plasmas, 14, 100703
- Meshkov, E. E. 1969, Fluid Dynamics, 4, 101
- Parker, E. N. 1957, J. Geophys. Res., 62, 509
- Pritchett, P. L. 2013, Physics of Plasmas, 20, 080703
- Richtmyer, R. D. 1960, Comm. Pure Appl. Math., 13, 297
- Samtaney, R., Loureiro, N., Uzdensky, D., Schekochihin, A., & Cowley, S. 2009, Phys. Rev. Lett., 103, 4
- Singh, C. B., de Gouveia Dal Pino, E. M., & Kadowaki, L. H. S. 2015, ApJ, 799, L20
- Sweet, P. A. 1958, in Electromagnetic Phenomena in Cosmical Physics (Cambridge, UK: Cambridge University Press), 123
- Turk, M. J., Smith, B. D., Oishi, J. S., et al. 2011, ApJS, 192, 9
- Uzdensky, D., Loureiro, N., & Schekochihin, A. 2010, Phys. Rev. Lett., 105, 235002
- Yamada, M., Kulsrud, R., & Ji, H. 2010, Rev. Mod. Phys., 82, 603