Magnetic reconnection, topological changes in magnetic fields, is a fundamental process in magnetized plasmas. It is associated with energy release in regions of magnetic field annihilation, but this is only one facet of this process. Astrophysical fluid flows normally have very large Reynolds numbers and are expected to be turbulent, in agreement with observations. In strong turbulence, magnetic field lines constantly reconnect everywhere and on all scales, thus making magnetic reconnection an intrinsic part of the turbulent cascade. We note in particular that this is inconsistent with the usual practice of magnetic field lines as persistent dynamical elements. A number of theoretical, numerical, and observational studies starting with the paper done by Lazarian and Vishniac [Astrophys. J. **517**, 700–718 (1999)] proposed that 3D turbulence makes magnetic reconnection fast and that magnetic reconnection and turbulence are intrinsically connected. In particular, we discuss the dramatic violation of the textbook concept of magnetic flux-freezing in the presence of turbulence. We demonstrate that in the presence of turbulence, the plasma effects are subdominant to turbulence as far as the magnetic reconnection is concerned. The latter fact justifies a magnetohydrodynamiclike treatment of magnetic reconnection on all scales much larger than the relevant plasma scales. We discuss the numerical and observational evidence supporting the turbulent reconnection model. In particular, we demonstrate that the tearing reconnection is suppressed in 3D, and unlike the 2D settings, 3D reconnection induces turbulence that makes magnetic reconnection independent of resistivity. We show that turbulent reconnection dramatically affects key astrophysical processes, e.g., star formation, turbulent dynamo, and acceleration of cosmic rays. We provide criticism of the concept of “reconnection-mediated turbulence” and explain why turbulent reconnection is very different from enhanced turbulent resistivity and hyper-resistivity and why the latter have fatal conceptual flaws.

## I. PROBLEM OF MAGNETIC RECONNECTION IN ASTROPHYSICS

Magnetic fields are ubiquitous in astrophysical systems and critically affect the dynamics and properties of magnetized plasmas over an extended range of scales. These scales are typically much larger than any relevant plasma scales, e.g., the ion inertial length. The magnetic field of the Earth's magnetosphere, the case-study of many *in situ* measurements, is an exception in the astrophysical context in terms of the involved scales. In what follows, we will discuss magnetic fields at much larger scales.

One textbook concept related to magnetic fields, widely used in astrophysical studies, is that magnetic fields in highly conducting plasmas are nearly perfectly frozen into the fluid and retain their topology for all time (Alfvén, 1942; Parker, 1979). This concept of flux-freezing is at the heart of many theories, e.g., the theory of star formation in the magnetized interstellar medium (ISM).

At the same time, there is ample evidence that magnetic fields in astrophysical systems do change their topology, with solar flares being a classical example of such a process (see, e.g., Parker, 1970; Lovelace, 1976; Priest and Forbes, 2002). The studies on the process enabling such a topology change, i.e., magnetic reconnection, were historically motivated by observations of the solar corona (Innes *et al.*, 1997; Yokoyama and Shibata, 1995; Masuda *et al.*, 1994), which also created a number of misconceptions. First of all, magnetic reconnection is frequently associated only with the processes of conversion of magnetic energy into heat. In addition, this understanding influenced attempts to find very peculiar conditions for flux conservation violation. For instance, in a fundamental study by Priest and Forbes (2002), a number of examples of magnetic configurations that produce fast reconnection are discussed. The attempts to accelerate the reconnection rate for special plasma conditions, e.g., in plasmas with very small collision rates (see, e.g., Shay *et al.*, 1998; Drake, 2001, 2006; Daughton *et al.*, 2006; Uzdensky and Kulsrud, 2006; Bhattacharjee *et al.*, 2003; Zweibel and Yamada, 2009; Yamada *et al.*, 2010), follow the same logic of seeking special circumstances where the magnetic flux-freezing is violated.

We feel that the above treatment of magnetic reconnection is not generic or widely applicable. Solar flares (Sturrock, 1966) are just one vivid example of reconnection activity. In reality, magnetic reconnection is a ubiquitous process taking place in various astrophysical environments, both collisionless and collisional (see Shibata and Magara, 2011). For instance, magnetic reconnection is a part of large-scale dynamo acting in stellar interiors (Parker, 1993; Ossendrijver, 2003). Magnetic reconnection is also required to make the theory of magnetohydrodynamic (MHD) turbulence (Goldreich and Sridhar, 1995) self-consistent (Jafari and Vishniac, 2019b).

The problem of magnetic reconnection is not limited to explaining its typically fast rates. In fact, it is often overlooked that the observations of solar activity indicate that the reconnection rates can significantly vary. This presents serious problems for theories that are based on relating the rate of magnetic reconnection to the peculiar properties of plasmas. These properties do not change, for instance, as the flux gets accumulated prior to a solar flare and gets annihilated during the flare. This may suggest that magnetic reconnection should have a sort of trigger.

It is also worth noting that for decades, a paradoxical situation existed with the treatment of astrophysical magnetic reconnection. On the one hand, the core of the reconnection community claimed numerous restrictions either on magnetic field configurations or the properties of reconnecting plasmas. On the other hand, the practitioners of numerical astrophysics used their MHD codes to simulate various astrophysical settings without worrying about the aforementioned restrictions. At the same time, the question of the actual effect of magnetic reconnection on the validity of numerical studies was somehow avoided.

The absence of an effective communication between the communities is surprising. For instance, for many years, it was considered that in order for reconnection to be fast, it must take place in collisionless environments. However, if collisionless reconnection were the only way to make reconnection rapid, then numerical MHD simulations of many astrophysical processes including those of the interstellar medium (ISM), which is collisional, are in error. Fortunately for numerical practitioners, the observations of the collisional solar photosphere indicate that the reconnection is fast in these environments as well (see Sec. VII A).

In this review, we provide evidence that, in fact, turbulence determines the rate of reconnection in realistic 3D astrophysical systems. The concept of 3D turbulent reconnection was introduced by (Lazarian and Vishniac, 1999, henceforth LV99). This model was followed by several subsequent theoretical studies, for instance, Eyink *et al.*, (2011, henceforth ELV11), Eyink (2015), and Jafari *et al.* (2018). The LV99 theory has been supported by numerical simulations in nonrelativistic (see Kowal *et al.*, 2009, 2012; Eyink *et al.*, 2013; Beresnyak, 2017; Oishi *et al.*, 2015; Kowal *et al.*, 2017) and relativistic settings (see Takamoto and Lazarian, 2016; Takamoto, 2018). In addition, different pieces of observational evidence, discussed in the review, support the predictions of the 3D turbulent reconnection theory of LV99.

To understand the essence of the LV99 model, one can recall that the problem of magnetic reconnection in astrophysical systems can be viewed as the problem related to the scale disparity. Reconnection occurs on very large scales, while the dissipation processes take place at the smallest plasma scales, which are set by, e.g., resistivity. LV99 theory shows how turbulence, as a multiscale phenomenon, can indeed resolve the problem of scale disparity and other problems that plague the traditional reconnection models.

While the idea of turbulent reconnection is nearly coeval with the ideas of 2D collisionless Hall-MHD reconnection, the latter has been substituted more recently by the ideas of tearing/plasmoid reconnection (Shibata and Tanuma, 2001). In fact, 2D numerical work has already shown that magnetic reconnection can be accelerated significantly by the tearing instability (Loureiro *et al.*, 2007; Lapenta, 2008; Daughton *et al.*, 2009a, 2009b; Bhattacharjee *et al.*, 2009; Cassak *et al.*, 2009).

In a sense, this model has some similarities to the turbulent reconnection. Indeed, it departs from the Sweet-Parker Y-point reconnection but allows instabilities to evolve within the current sheet. We will discuss how in realistic 3D configurations this type of instability transfers to turbulent reconnection.

One should keep in mind that we deal with 3D turbulent reconnection on large scales at which an MHD-like approximation is applicable. In particular, we consider scales much larger than the ion Larmor radius or skin depth, where an ideal Ohm's law is generally agreed to be valid. At the same time, usual tests of reconnection are based on the *in situ* measurements of magnetospheric reconnection. This type of reconnection occurs on scales comparable to the ion inertial length, and therefore, it is atypical for the large-scale reconnection that occurs in most astrophysical systems. For the reconnection explored at larger scales, e.g., in the solar wind (Gosling, 2012; Lalescu *et al.*, 2015), the measured properties of reconnection are consistent with the general expectations based on the LV99 theory (see Sec. VII B).

Occasionally, the process of magnetic reconnection is understood in a very narrow sense, e.g., magnetic reconnection necessarily in low-*β* plasmas with the magnetic energy much larger than the thermal energy. Here, we consider turbulent reconnection as a generic process, which is applicable to plasmas of arbitrary *β* and for arbitrary arrangements of reconnecting magnetic fluxes.

Our review is devoted to the astrophysical reconnection, although the magnetic reconnection is a fundamental physical process, and therefore, our conclusions can be applied also to the laboratory reconnection. We discuss the theory and observations of astrophysical turbulence in Sec. II. The LV99 model of 3D turbulent reconnection and its extensions are described in Secs. III and IV. There we explain how turbulent reconnection and flux-freezing violation in turbulent fluids are inter-related. The numerical tests of turbulent reconnection and the testing of flux freezing violation in turbulent fluids are provided in Sec. V. Testing with both MHD codes and Hall-MHD codes are discussed. The initial tests of selected ideas of turbulent reconnection with 3D kinetic simulations are presented in Sec. VI. The observational evidence supporting turbulent reconnection is briefly discussed in Sec. VII. The magnetic reconnection in a high-Prandtl number medium, reconnection in relativistic fluids, and reconnection in the presence of whistler turbulence are discussed in Sec. VIII. Astrophysical implications of reconnection for turbulent dynamo, star formation, as well as for particle acceleration and gamma ray bursts are covered in Sec. IX. In Sec. X, we discuss the differences between our model of turbulent reconnection and other ideas on the effects of turbulence on reconnection and resistivity. We explain the reasons for our disagreement with (i) the currently popular reconnection-mediated turbulence idea and (ii) the concepts of turbulent resistivity and turbulent ambipolar diffusion. There we also outline the problems of the mean field approach to study reconnection, compare turbulent and tearing reconnection, and explain why reconnection rates are not expected to have a “universal” value of 0.1 Alfvén velocity. We present our summary in Sec. XI.

## II. TURBULENCE AS A NATURAL STATE OF ASTROPHYSICAL FLUIDS AND ITS MHD DESCRIPTION

### A. Observational evidence

Observations of the diffuse warm ISM reveal a Kolmogorov spectrum of electron density fluctuations (see Armstrong *et al.*, 1995; Chepurnov and Lazarian, 2010). Similar spectral slopes of supersonic velocity fluctuations are measured using Doppler shifted spectral lines (see Lazarian, 2009, for a review) and through the *in situ* measurements of the solar wind fluctuations (Leamon *et al.*, 1998). The evidence of turbulence being ubiquitous comes from the nonthermal broadening of spectral lines and measures obtained by other techniques (see Burkhart *et al.*, 2010). This fact is not surprising as magnetized astrophysical plasmas generally have very large Reynolds numbers. Indeed, the length scales involved are large, while the motions of charged particles in the direction perpendicular to magnetic fields are of the order of the Larmor radius. Plasma flows at these high Reynolds numbers are prey to numerous linear and finite-amplitude instabilities. This induces turbulent motions readily.

The precise origin of the plasma turbulence differs from instance to instance. It is sometimes driven by an external energy source, such as supernova explosions in the ISM (Norman and Ferrara, 1996; Ferrière, 2001), merger events and active galactic nuclei (AGN) outflows in the intercluster medium (ICM) (Subramanian *et al.*, 2006; Enßlin and Vogt, 2006; Chandran, 2005), and baroclinic forcing behind shock waves in interstellar clouds. In other cases, turbulence is spontaneous, with available energy released by a rich array of instabilities, such as magnetorotational instability (MRI) in accretion disks (Balbus and Hawley, 1998; Jafari and Vishniac, 2018b) and kink instability of twisted flux tubes in the solar corona (Galsgaard and Nordlund, 1997; Gerrard and Hood, 2003). In all these mentioned cases, turbulence is not driven by reconnection. Nevertheless, we would like to stress that the driving of turbulence through the energy release in the reconnection zone is important for the transfer from laminar to turbulent reconnection. Altogether, whatever its origin is, the signatures of plasma turbulence are seen throughout astrophysical media.

Indeed, observations show that turbulence is ubiquitous in all astrophysical plasmas. The spectrum of electron density fluctuations in the Milky Way is presented in Fig. 1, as one dramatic piece of evidence for cosmic turbulence. Similar spectra are reported in Leamon *et al.* (1998) and Bale *et al.* (2005) for solar wind, Padoan *et al.* (2006) for molecular clouds, and Vogt and Enßlin (2005) for the ICM. New techniques for studying turbulence, e.g., the Velocity Channel Analysis (VCA) and Velocity Coordinate Spectrum (VCS) techniques (Lazarian and Pogosyan, 2000, 2004, 2006), have provided important insights into the velocity spectra of turbulence in molecular clouds (see Padoan *et al.*, 2006, 2010), Galactic, and extragalactic atomic hydrogen [Stanimirović and Lazarian, 2001; Chepurnov *et al.*, 2010, 2015, see also the review by Lazarian (2009), where a compilation of velocity and density spectra obtained with contemporary HI and CO data are presented].

### B. Basic phenomenology of MHD turbulence

It is generally accepted that an MHD-like description of turbulence is applicable to the motions of magnetized plasmas at sufficiently large scales (e.g., see discussion in Schekochihin *et al.*, 2009). Standard MHD or a multispecies variant is valid at length-scales much larger than the collisional mean-free-path of the plasma. However, even in a nearly collisionless but magnetized plasma, “kinetic MHD” or “kinetic reduced-MHD” equations are satisfied at scales much larger than the ion gyroradius (Kulsrud, 1983; Kunz *et al.*, 2015). Since collisions in those plasmas are too infrequent to bring the particle distributions to Maxwellian, the pressure tensor is anisotropic and given not by a thermodynamic equation of state but instead by 1D kinetic equations along particle guiding line-centers. This subtlety is irrelevant for the incompressible shear-Alfvén wave modes, which nearly uncouple from the kinetic/compressible modes and satisfy the well-known “reduced MHD” (RMHD) equations [for relativistic turbulence (see Takamoto and Lazarian, 2017) the transfer of energy to fast modes is more important]. Shear-Alfvén waves are often energetically dominant, as in the solar wind where they contain about 90% of the total energy of plasma and fields (see Dobrowolny *et al.*, 1980). For this reason, we shall generally assume that MHD is an accurate model at sufficiently large scales or at least a good working approximation within which to explore the effects of turbulence on large-scale reconnection.

We shall invoke below in our discussion of strong MHD turbulence the (Goldreich and Sridhar, 1995, henceforth GS95) theory of Alfvénic turbulence. We believe that the empirical evidence largely supports GS95 over its competitors and gives a good description of MHD turbulence (apart from the phenomenon of small-scale intermittency). However, we wish to avoid here any controversy regarding the “correct” phenomenology. It may be shown following LV99, Appendix D, that the predictions on turbulent reconnection are qualitatively independent of the particular phenomenology adopted and, in particular, do not depend sensitively upon the predicted spectral slopes. Therefore, our use of GS95 theory may be regarded by its persistent skeptics as a mere convenience since it is the simplest model that incorporates the small-scale anisotropy due to a local magnetic field.

#### 1. Derivation of the MHD turbulence relations

MHD turbulence was traditionally viewed as the process dominated by the Alfvén wave interactions. Below, we follow this path to show how to obtain GS95 relations and go beyond them (see also Cho *et al.*, 2003a).

For Alfvénic perturbations, the relative perturbations of velocities and magnetic fields are related in the following way:

where *B _{l}* is the perturbation of the magnetic field

*B*at the scale

*l*,

*B*is the perturbation of the magnetic field at the injection scale

_{L}*L*, while

*u*is the velocity fluctuation at the scale

_{l}*l*in the turbulent flow with energy injected with velocity

*u*.

_{L}To understand Alfvén wave damping, it is advantageous to present MHD turbulence as a superposition of colliding Alfvén wave packets. Consider such packets with parallel scales $l||$ and perpendicular scales $l\u22a5$. The system of reference with respect to the local magnetic field is what is relevant when dealing with MHD turbulence. In the rest of the discussion, we operate with $l\Vert $ and $l\u22a5$ which are defined in the frame of wavepackets moving along the local magnetic field. The wave packet collision induces a change

where the first term is the change of the energy of a packet as it interacts with the oppositely moving wave packet. The time of the interaction should be identified with the time of the passage of the given wave packet through the oppositely moving packet of the size $l||$. Thus, the interaction time is given by $\Delta t\u223cl||/VA$.

The characteristic rate of the cascade is due to the change in the structure of the oppositely moving wave packet. The latter has the rate $ul/l\u22a5$. As a result, Eq. (2) gives

The fractional energy change per collision is $\Delta E/E$, which measures the strength of the nonlinear interaction

In Eq. (4), *f* measures the ratio of the rate of shearing of the wave packet $ul/l\u22a5$ to the rate of the propagation of the wave packet $VA/l||$. Two cases can be identified. For $f\u226a1$, the shearing rate is much smaller than the propagation rate, and the cascade process is a random walk. This means that

which is the number of steps required for the O(1) cascade of energy, and therefore, the cascade time is

For $\u2135\u226b1$, the turbulence corresponds to the “weak MHD turbulence.” Naturally, $\u2135$ cannot be smaller than unity. The limiting case of $\u2135\u22481$ corresponds to “strong MHD turbulence.”

Consider first the case of weak MHD turbulence, which requires $VA\u226auL$ and very high-frequency oscillations. This cascade thus proceeds by resonant 3-wave Alfvénic interactions that satisfy $\omega 1+\omega 2=\omega 3$ in addition to $k1+k2=k3,$ where $k$ and *ω* denote the wavevector and angular frequency, respectively. Because of the dispersion relation $\omega =sVAk\Vert $ with $s=\xb11$ for shear-Alfvén waves traveling parallel and antiparallel to magnetic fields and because only counterpropagating waves interact, one finds that $k2,||=0$ and $k1,||=k3,||.$ Thus, resonantly interacting wavepackets have unchanging $k||$, but $k\u22a5$ increases by distortions through collisions. The decrease in $l\u22a5$ while $l||$ does not change signifies the increase in the energy change per collision. This eventually makes $\u2135$ of the order of unity and the cascade transitions to strong MHD turbulence. In this case, one gets

which signifies the cascade time being equal to the wave period $\u223c\Delta t$. It was hypothesized by Goldreich and Sridhar (1995) that the further decrease in $l\u22a5$ entails the corresponding decrease in $l||$ to keep Eq. (7) satisfied, which is the “critical balance condition.” The ability to change $l||$ by the nonresonant strong cascade means that the frequencies of interacting waves increase. However, this increase has a natural limit because a too high frequency would restore the resonance condition $\omega 1+\omega 2=\omega 3$, and this would shut off the parallel cascade. Thus, the strong turbulent cascade has a tendency to operate at marginal resonance, where the critical balance relation (7) is satisfied.

The cascade of turbulent energy satisfies the relation (Batchelor, 1953)

which for the hydrodynamic cascade provides

where the relation $tcas\u2248l/ul$ is used.

For the weak cascade $\u2135\u226b1$

where Eqs. (8) and (6) are used (see LV99). The isotropic turbulence injection at scale $\u2113\Vert \u2243\u2113\u22a5\u2243L$ results in the third relation in Eq. (10). Taking into account that $l||$ is constant for weak turbulence, it is easy to see that Eq. (10) provides

which is different from the hydrodynamic $\u223cl1/3$ scaling [using the relation $kE(k)\u223cuk2$ it is easy to show that the spectrum of weak turbulence is $Ek,weak\u223ck\u22a5\u22122$ (LV99; Galtier *et al.*, 2000)].

For “sub-Alfvénic turbulence” injected isotropically at length scale *L* with *u _{L}* <

*V*, the transition to the strong regime, i.e., $\u2135\u22481,$ happens at the scale (LV99)

_{A}Therefore, weak turbulence has a limited, i.e., $[L,LMA2]$, inertial range and transits at *l _{trans}* into strong turbulence. The velocity at the transition follows from the $\u2135\u22481$ condition given by Eqs. (5) and (4)

The scaling for the strong turbulence range at $l<ltrans$ can be easily obtained. Indeed, the turbulence cascades over one wave period, which according to Eq. (23), are equal to $l\u22a5/ul$. Substituting the latter in Eq. (8), one gets

The latter is analogous to the hydrodynamic Kolmogorov cascade in the direction perpendicular to the local direction of the magnetic field. This strong MHD turbulence cascade starts at *l _{trans}*, and it has the injection velocity given by Eq. (13). This provides the another way to get the scaling relations for velocity in strong Alfvénic turbulence (LV99)

which can be rewritten in terms of the injection velocity *u _{L}* as

Substituting this in Eq. (7), we get the relationship between the parallel and perpendicular scales of the eddies (LV99)

The relations Eqs. (17) and (15) reduce to the GS95 scaling for trans-Alfvénic turbulence if $MA\u22611$.

In the opposite case of “super-Alfvénic turbulence,” i.e., for $uL>VA,$ the scales close to the injection scale are essentially hydrodynamic as the influence of magnetic forces is minimal. Thus, the velocity is Kolmogorov

The cascade changes its nature at the scale

for which the turbulent velocity is equal to the Alfvén velocity (Lazarian and Pogosyan, 2016). The cascade rate for $l<lA$ is

This case likewise reduces to trans-Alfvénic turbulence but with *l _{A}* acting as the injection scale. At scales $\u2113<lA$

and

The wave description of strong MHD turbulence is very productive in the case when the turbulence is imbalanced, i.e., the flow of wave energy from one direction along field lines significantly exceeds the flow from the other direction. An example of imbalanced MHD theory that agrees with numerical simulation is the theory by Beresnyak and Lazarian (2008) (see also Beresnyak and Lazarian, 2010).

#### 2. Eddy description of MHD turbulence

In our above discussion, we revealed the analogy between the Alfvénic turbulence and Kolmogorov turbulence. This analogy can be extended in the case of strong Alfvénic turbulence. In fact, the GS95 relations can be naturally obtained if one assumes that such turbulence presents a collection of eddies that are mixing the conducting fluid perpendicular to the direction of the magnetic field present at the location of the eddies. This description can be justified in the presence of fast reconnection that we are describing in the review, and this supports our claim that the MHD turbulence and turbulent magnetic reconnection are intrinsically connected.

As it was described in LV99, the reconnection of magnetic fields associated with an eddy takes place within the eddy turnover time (see also Sec. III). In this situation, the random driving preferentially excites the motions that induce minimal bending of magnetic field lines, i.e., the hydrodynamic type motions of eddies that mix magnetized fluid perpendicular to the direction of the magnetic field. This is the picture of MHD turbulence that follows from the LV99 description. The key concept in the eddy description is that magnetic perturbations are measured in terms of the local direction of the magnetic field, not the mean field. This concept is missing in GS95, and this causes a lot of confusion with different authors attempting to test the GS95 scale-dependent anisotropy in the frame of the mean magnetic field. In fact, the latter anisotropy only exists in the local frame of the eddies, while the anisotropy in the reference frame of the mean magnetic field is scale-independent and is determined by the anisotropy of the largest eddies (see Cho *et al.*, 2003a).

Incidentally, within the description suggested in LV99, it is clear that the gradients of both velocities and magnetic field are expected to be perpendicular to the local magnetic field direction. Therefore, by tracing the directions of the gradients, one can trace the magnetic field direction in turbulent media. This idea is at the foundations of the new technique of magnetic field tracing that employs gradients of intensities within spectroscopic channel maps (Lazarian and Yuen, 2018b), gradients of synchrotron intensity (Lazarian *et al.*, 2018), and gradients of synchrotron polarization (Lazarian and Yuen, 2018a). The magnetic field tracing by the new technique provides results in good agreement with polarization data (see Hu *et al.*, 2019), which is an indirect way of confirming both the ubiquitous presence of the MHD turbulence in the interstellar medium and a successful confirmation of the theoretical predictions.

Adopting the picture of eddies perpendicular to the local magnetic field of the eddies, one can claim that in terms of perpendicular motions, we expect to have the Komogorov picture, i.e., $u\u2113\u223cl\u22a51/3$. The critical balance relations in the picture of the eddies also follow naturally. Indeed, as the eddy provides the mixing of magnetic field lines in the perpendicular direction over the time scale $l\u22a5/u\u2113$, it sends an Alfvén wave with the same period. This is exactly

where $l||$ and $l\u22a5$ are eddy axes parallel and perpendicular to the “local” direction of the magnetic field and the magnetic field of the eddy. Note that this requires a more sophisticated way to analyze numerical simulations. For instance, a way of determining the local magnetic field in numerical simulations using structure functions is discussed in the study by Cho and Vishniac (2000). In fact, the notion of the local magnetic field is the essential part of the modern understanding of Alfvénic turbulence, and it was added to the MHD turbulence theory by the later studies (LV99, Cho and Vishniac, 2000; Maron and Goldreich, 2001). This notion is natural within the picture of LV99 turbulent reconnection. Indeed, if one considers the aforementioned eddy picture, the use of the local magnetic field means that the eddies are only affected by the magnetic field around them and do not know anything about the mean magnetic field obtained by averaging over the entire turbulent volume.

It is assumed in GS95 that the energy is injected at the scale *L* with the injection velocity equal to the Alfvén velocity. If the injection velocity *u _{L}* is smaller than $VA,$ then the eddies at a small scale get more elongated. For $MA<1$, the eddy description is valid only for the strong turbulence, i.e., starting with the scales less than $ltrans=LMA2$. Taking

*l*as the injection scale of the turbulence and the velocity at this scale $VAMA2$ [see Eq. (13)] as the injection velocity, one can repeat the arguments about the correspondence of the time scales of parallel and perpendicular motions to obtain the relations given by Eqs. (15) and (17).

_{trans}#### 3. Compressible MHD turbulence

The compressible MHD turbulence can be approximated by the superposition of 3 cascades of basic MHD modes, namely, Alfvénic, slow, and fast. The original idea of the decomposition of low amplitude MHD perturbations into the modes can be traced to (Dobrowolny *et al.*, 1980), while the practical realization for the actual trans-Alfvénic turbulence was first performed in (Cho and Lazarian, 2002). We, however, will concentrate on the Alfvénic part of the MHD turbulent cascade, as this component is most important for the theory of reconnection that we consider in the review. In nonrelativistic MHD turbulence, the backreaction of slow and fast magnetosonic modes (Cho and Lazarian, 2002, 2003; Kowal and Lazarian, 2010) on the Alfvénic cascade is insignificant (Cho and Lazarian, 2002; Goldreich and Sridhar, 1995; Lithwick and Goldreich, 2001). Therefore, in nonrelativistic MHD turbulence, the above relations obtained for Alfvénic turbulence are not much affected by the presence of compressible modes.

In the rest of this section, we discuss how the partial ionization and relativistic effects change the picture of Alfvénic turbulence.

### C. MHD turbulence in a partially ionized medium

In partially ionized astrophysical plasmas in, e.g., the early universe, cold interstellar phases, protoplanetary disks, and solar atmosphere, the MHD turbulence is subjected to the damping effects due to the presence of neutrals. Since we are interested in the turbulent reconnection associated with Alfvénic turbulence, here we focus on the damping of Alfvénic turbulence.

The damping effects in a partially ionized medium mainly arise from the collisional friction between ions and neutrals, i.e., ion-neutral collisional damping (IN) and the viscosity in neutral fluid, i.e., neutral viscous damping (NV) (Lithwick and Goldreich, 2001; Lazarian *et al.*, 2004). Their relative importance depends on the turbulence and plasma parameters (Xu *et al.*, 2015, 2016; Xu and Lazarian, 2017b). In super-Alfvénic turbulence with the turbulent energy larger than the magnetic energy, IN dominates over NV if the condition

is satisfied, where *ξ _{n}* is the neutral fraction,

*ν*is the neutral-ion collisional frequency,

_{ni}*ν*is the neutral viscosity,

_{n}*L*is the driving scale of turbulence, and

*u*is the turbulent velocity at

_{L}*L*. In sub-Alfvénic turbulence with the turbulent energy smaller than the magnetic energy, under the condition

IN is more important than NV. Here, $MA=uL/VA$ is the Alfvén Mach number and *V _{A}* is the Alfvén speed. An illustration of the above condition in sub-Alfvénic turbulence can be found in Fig. 2, where the typical parameters of the interstellar turbulence are used (Xu and Lazarian, 2017b). Xu

*et al.*(2016) found that in the warm neutral medium of the ISM and the solar chromosphere, NV plays a significant role in damping MHD turbulence.

When the damping effects are significant, MHD turbulence is suppressed. The damping scale where the MHD turbulence is dissipated can be determined by comparing the turbulent energy cascade rate with the damping rate. Different from the damping scale of linear MHD waves, the damping scale of MHD turbulence has parallel and perpendicular components, $kdam,||$ and $kdam,\u22a5$, due to the anisotropy of MHD turbulence. In the frame of local magnetic field, i.e., the magnetic field averaged over the length scale of interest, magnetic fields are mixed by turbulent motions in the perpendicular direction, and the perturbed magnetic fields have wavelike motions in the parallel direction. For strong MHD turbulence, the lifetime of waves is equal to the eddy turnover time. As the turbulence anisotropy increases toward small length scales, when the damping scale of MHD turbulence is much smaller than *L*, it can be approximated by its perpendicular component in terms of magnitude.

Under the damping effects of both IN and NV, a general expression of the damping scale is (Xu and Lazarian, 2017b)

for super-Alfvénic turbulence and

for sub-Alfvénic turbulence, where *V _{Ai}* is the Alfvén speed of the ionized gas,

*ν*is the ion-neutral collisional frequency, and $lA=LMA\u22123$ in super-Alfvénic turbulence.

_{in}When IN is the dominant damping effect, the damping is related to the coupling state between ions and neutrals. The most severe damping arises in the weakly coupled regime in a weakly ionized medium, where neutrals decouple from ions, but ions are still coupled with neutrals due to their frequent collisions with surrounding neutrals. In this situation, the (perpendicular) damping scale for super- and sub-Alfvénic turbulence can be simplified as (Xu and Lazarian, 2017b)

and

Its relation to the neutral-ion decoupling scale is

for both super- and sub-Alfvénic turbulence. The difference between the decoupling and damping scales has been confirmed by two-fluid MHD simulations (Burkhart *et al.*, 2015). The above relation shows that the decoupling of neutrals from ions plays a key role in the IN damping of MHD turbulence.

After neutrals decouple from ions, the hydrodynamic turbulence carried by neutrals is only subject to viscous damping, and its NV damping scale is smaller than the IN damping scale of MHD turbulence in ions. This differential damping of neutrals and ions can account for the observed linewidth difference of the coexistent neutrals and ions in molecular clouds (Li and Houde, 2008; Xu *et al.*, 2015).

When NV dominates over IN, the damping scale is approximately (Xu and Lazarian, 2017b)

for super-Alfvénic turbulence and

for sub-Alfvénic turbulence. MHD turbulence is damped at the NV damping scale. On the length scales smaller than the NV damping scale and larger than the IN damping scale, magnetic fluctuations persist, and the magnetic energy has a spectral form $k\u22121$. This “new regime of MHD turbulence” in the subviscous range was first analytically predicted by Lazarian *et al.* (2004) and further confirmed by MHD simulations (Cho *et al.*, 2002a, 2003b). The $k\u22121$ spectrum of magnetic energy in the subviscous range has also been observed in dynamo simulations in the nonlinear stage of dynamo (Haugen *et al.*, 2004; Schekochihin *et al.*, 2004).

### D. Relativistic MHD turbulence

Due to its numerical and theoretical simplicity, MHD turbulence in the relativistic force-free regime has been studied first. Relativistic force-free formalism can be used for a system, such as the magnetosphere of a pulsar or a black hole, in which the magnetic energy density is much larger than that of matter.

Scaling relations for relativistic Alfvénic MHD turbulence were first derived by Thompson and Blaes (1998). The predictions of these theories in terms of the spectral slope and anisotropy coincide with those in GS95 theory. These relations were first numerically tested by Cho (2005) who performed a numerical simulation of a decaying relativistic force-free MHD turbulence with a numerical resolution of 512^{3} and calculated the energy spectrum and anisotropy of eddy structures. After *t *=* *3, the energy spectrum with $E\u223ck5/3$ was shown to decrease in amplitude without changing its slope. The anisotropy was also measured and was shown to be consistent with the GS95 expectations, i.e., $k||\u223ck\u22a52/3$. Thus, in terms of Alfvénic turbulence, the correspondence between the nonrelativistic theory and its relativistic counterpart is excellent. Interestingly enough, the relativistic simulations of imbalanced Alfvénic turbulence, i.e., the turbulence with the excess of the energy flux moving in one direction, performed by Cho and Lazarian (2014) provided results in agreement with the theory of nonrelativistic imbalanced turbulence in Beresnyak and Lazarian (2008).

The situation is somewhat different in the case of relativistic compressible theory. The initial studies of relativistic MHD turbulence (see Zhang *et al.*, 2009; Inoue *et al.*, 2011; Beckwith and Stone, 2011; Zrake and MacFadyen, 2012, 2013; Garrison and Nguyen, 2015) have not revealed much differences between the nonrelativistic turbulence and relativistic turbulence. However, more recent studies by Takamoto and Lazarian (2016) and (2017) provided the decomposition of the turbulent motions into Alfvén, slow, and fast modes. The comparison of the obtained results with the analogous decomposition into modes in the nonrelativistic case in the study by Cho and Lazarian (2002) and (2003) and Kowal and Lazarian (2010) revealed significant differences in terms of the coupling of the compressible and incompressible motions in the two cases. These differences are illustrated in Fig. 3 where the transfer of energy between Alfvén and fast modes is shown in the relativistic case. Note that the corresponding transfer of energy is suppressed for nonrelativistic compressible turbulence. In terms of anisotropy, this energy transfer makes fast modes more similar to Alfvén ones. We find that this effect of energy transfer is important for describing turbulent relativistic reconnection.

### E. Turbulence: Summary of basic facts

For our discussion of turbulent reconnection, it is important that turbulence is ubiquitous in astrophysical environments and that Alfvénic modes usually represent the dominant component of the MHD turbulent cascade. These modes are present in collisional and collisionless environments. The turbulence in most astrophysical settings has the sources different from the magnetic reconnection. However, it is difficult to understand the properties of turbulent fluids if magnetic reconnection does not happen on the dynamical time of the fluid motions.

The partial ionization of plasmas changes the properties of MHD turbulence increasing the scale at which MHD turbulence is damped. The relativistic MHD turbulence and nonrelativistic MHD turbulences have significant similarities, although they show a difference in the coupling of fast and Alfvén modes. A more detailed description of the theory and implications of MHD turbulence can be found in a book by Beresnyak and Lazarian (2019).

## III. ANALYTICAL MODEL OF TURBULENT RECONNECTION

### A. Derivation of the LV99 reconnection rate

The 3D model of turbulent reconnection proposed in LV99 provides a natural generalization of the classical Sweet-Parker model (Parker, 1957; Sweet, 1958). Thus, we remind the reader first the basic properties of this model.

In the Sweet-Parker model, two regions with uniform “laminar” magnetic fields are separated by a thin current sheet. The magnetic fields are frozen in through the two regions, and the violation of the flux freezing takes place over a thin slot *δ*, with the latter being determined by Ohmic diffusion. Thus, the speed of reconnection is determined roughly by the resistivity divided by the sheet thickness, i.e.,

Obviously, the smaller the Δ value, the larger the reconnection speed. Small Δ, however, prevents plasmas from leaving the reconnection region.

The Sweet-Parker model deals with the “steady state reconnection,” and therefore, the plasma in the diffusion region must be ejected from the edge of the current sheet at the Alfvén speed, *V _{A}*. As a result, the reconnection speed is limited by the conservation of the mass condition

where *L _{x}* is the lateral extend of the current sheet.

The conditions imposed by Eqs. (33) and (34) act in the opposing directions as far as the thickness of the reconnection region Δ is concerned. As a result, the Sweet-Parker reconnection cannot be fast as it faces contradictory requirements. Therefore, the compromise speed of reconnection is the Sweet-Parker reconnection speed that is reduced from the Alfvén speed by the square root of the Lundquist number, $S\u2261LxVA/\eta $, i.e.,

The Lundquist number in astrophysical conditions is enormous, and it can be 10^{16} or larger. In this situation, the Sweet-Parker reconnection rate is absolutely negligible.

Due to the disparity between *L _{x}* and Δ, the Sweet-Parker model cannot explain most of the observed astrophysical reconnection events, e.g., solar flares. A way to deal with this problem was suggested by Petschek (1964), who speculated that in special conditions, the reconnection configurations may have magnetic field lines converging to the reconnection zone at a sharp angle making

*L*and Δ comparable. The corresponding configuration of the magnetic field is known as X-point reconnection as opposed to the Sweet-Parker extended current sheet Y-type reconnection. The research exploring X-point reconnection was the dominant trend in the attempts of obtaining fast reconnection for decades. However, the major problem of such configurations is to preserve the X-point configuration over astrophysically large scales.

_{x}The direction taken by LV99 was to modify the Y-type reconnection in order to make Δ macroscopic and even astrophysically large. Naturally, this is impossible with plasma effects. However, LV99 proposed that turbulence can induce magnetic field wandering that can do the job. The corresponding model of magnetic reconnection is illustrated in Fig. 4.

Similar to the Sweet-Parker model, the turbulent model developed in LV99 deals with a generic configuration, which arises naturally as magnetic flux tubes crossing each other make their way one through another. The outflow thickness Δ is determined by the large-scale magnetic field wandering. This wandering depends on the level of turbulence. As a result, the rate of reconnection in the LV99 model is determined by the level of MHD turbulence. The LV99 reconnection rate is fast, i.e., its rate does not depend on the fluid/plasma resistivity. However, the rates of turbulent reconnection can change significantly depending on the level of turbulence. This is an important prediction of LV99 theory which makes it different from the theories that, for instance, appeal to plasma effects for describing fast reconnection.

To quantify the reconnection rate in the LV99 model, one should use the scaling relations for Alfvénic turbulence from Sec. II. In fact, it is reconnection at low Alfvén Mach numbers which corresponds to the classical representation of magnetic reconnection flux tubes. This in fact required LV99 to generalize the GS95 MHD theory for arbitrary Alfvén Mach numbers.

As the first step consider the problem of magnetic field wandering, as this magnetic wandering determines the thickness of outflow Δ. We may note that magnetic wandering determined in LV99 is not only important for magnetic reconnection. In fact, it is a fundamental property of the turbulent magnetic field, and it is important for different problems, e.g., problems of heat transfer (Narayan and Medvedev, 2001, Lazarian, 2006) and cosmic ray propagation (Yan and Lazarian, 2008, Lazarian and Yan, 2014).

Using the MHD turbulence relations, one can describe the diffusion of magnetic field lines. Indeed, a bundle of field lines confined within a region of width *y* at some particular point spreads out perpendicular to the mean magnetic field direction as magnetic field lines are traced in either direction. The rate of field line spread is given by (see more in Sec. III B 3)

where $\lambda ||\u22121\u2248\u2113||\u22121,\u2009\u2113||$ is the parallel scale. The corresponding transversal scale, $\u2113\u22a5$, is $\u223c\u27e8y2\u27e91/2$, and *x* is the distance measured along the magnetic field direction. As a result, using Eq. (21), one gets

where the substitution $\u27e8y2\u27e91/2$ for $\u2113\u22a5$ is used. This expression for the magnetic field dispersion is applicable only when *y* is small enough for us to use the strong turbulence scaling relations, i.e., when $\u27e8y2\u27e9<Li2(uL/VA)4$.

For scales less than *L _{i}*, it follows from Eq. (37) that the mean squared magnetic field wandering is equal to

which provides a superdiffusive dependence of $\u27e8y2\u27e9\u223cs3$, which was one of the major findings in LV99 as far as the properties of the magnetic field in turbulence are concerned. Due to the above scaling, the probability of the wandering magnetic field line to return back to the reconnection layer is very small.

The dependence of $\u27e8y2\u27e9\u223cMA4$ is significant as it shows that if the turbulent kinetic energy is a small portion of magnetic energy of the fluid, the magnetic field lines are nearly straight. This provides a smooth transition to the slow reconnection in the absence of turbulence and, as we discuss later, to the flux freezing being a fair approximation to the low *M _{A}* fluids.

For describing magnetic field wandering over scales larger than the turbulence injection scale, i.e., $s\u226bLi$, magnetic field wandering obeys the usual random walk scaling. The corresponding random walk step is $s/Li$ with the mean squared displacement per step equal to $Li2(uL/VA)4$. As a result

which provides the other limit for the magnetic field wandering. While the diffusive regime of magnetic field lines given by Eq. (39) was well known in the cosmic ray literature (see Schlickeiser, 2002), the superdiffusive one was introduced in LV99 [using the relation $kE(k)\u223cuk2$ it is easy to show that the spectrum of weak turbulence is $Ek,weak\u223ck\u22a5\u22122$ (LV99; Galtier *et al.*, 2000)].

Using Eqs. (38) and (39), one can describe the thickness of the outflow for the reconnection of the antiparallel magnetic fields in turbulent media, i.e., in the setting generalizing Sweet-Parker reconnection without the guide field for turbulent media. Measuring *s* along the horizontal *x*-axis along the magnetic field reversal and combining Eqs. (38) and (39), it is possible to derive the thickness of the outflow Δ in the aforementioned two regimes. Then, using the mass conservation Eq. (52), it is straightforward to obtain

where $VAMA2$ is proportional to the turbulent eddy speed. As we discussed earlier, the obtained reconnection rate varies depending on Alfvén Mach number *M _{A}* and for $MA\u223c1$ can represent a large fraction of the Alfvén speed in the case when

*L*and

_{i}*L*are not too comparable.

_{x}Similar to the Sweet-Parker picture, the contribution of the guide field does not change the dynamics of magnetic reconnection. Figure 4 presents a *XY* plane of the cross section of the reconnection region perpendicular to the guide field *B _{z}* shared by two fluxes. The guide field should be ejected together with plasma from the reconnection region. Indeed, both the random velocity and magnetic field can be decomposed into components in the direction of the guide field and perpendicular to it. The field wondering in the direction perpendicular to the guide field is affected by the velocity perturbations in the same direction. The guide field does not oppose to the corresponding bending of the magnetic field. The outflow velocity is, however, determined only by the Alfvén velocity corresponding to the perpendicular component of the magnetic field. In fact, for the sake of clarity, it could be right to change

*M*for $MA,XY$, where the latter denotes the Alfvén Mach number corresponding to the Alfvén velocity in the plane perpendicular to the guide field. This notation of the Alfvén Mach number is not common, and we will not use it therefore.

_{A}It was communicated to us by Andrey Beresnyak that the independence of the reconnection rate on the guide field can be also justified from the point of the Reduced MHD (RMHD) approach (see Beresnyak, 2012). Formally, the RMHD is applicable for describing incompressible fluids when the perturbations are much smaller than *V _{A}*. However, in practical cases, RMHD is well applicable for many collisionless systems, e.g., solar wind (see Schekochihin

*et al.*, 2009), and its practical validity is good even for trans-Alfvénic turbulence (Beresnyak

*et al.*, 2009). Within the RMHD, the dynamics of turbulent motions is not changed if the ratio $B0/l||$ stays the same. If no constraints are imposed on the parallel scale of perturbations $l||$, they increase with the increase in the guide field

*B*

_{0}as was demonstrated in numerical simulations by Beresnyak (2017). The dynamics in terms of perpendicular perturbations $l\u22a5$ and field wondering that these perturbations entailed increasing the reconnection outflow stayed the same.

For distances along the magnetic field less than the injection scale *L _{i}* using the definition of the cascading rate for strong MHD turbulence given by Eq. (14), we can rewrite the expression for magnetic field wandering [see Eq. (38)] as

from which it is clear that the increase in *V _{A}* can be naturally compensated by the increase in the parallel scale along the magnetic field.

For the scales larger than the injection scale, Eq. (39) can be rewritten as

where the expression for the cascading of the weak MHD cascade given by Eq. (10) is used. It is also obvious that the simultaneous increase in the parallel scale of both injection and turbulent motions and the Alfvén velocity does not change the magnetic field wandering.

For the formal applicability of RMHD, one should assume that the guide field *B _{z}* is much stronger than the reconnecting magnetic field

*B*. The relevant setting corresponds to the reconnection of two flux tubes intersecting at a small angle

_{x}*α*. The corresponding setting was discussed in LV99 for the reconnection induced by individual eddies within MHD turbulence, and it was shown there that the reconnection happens within one eddy turnover time. Here, we consider a more generic situation of reconnecting flux tubes with an arbitrary turbulence driving.

In the above setting, $Bx\u2248Bz\u2009sin\u2009\theta /2$, and the perturbations at the injection scale *L _{i}* result in the perturbations in the

*XY*plane with the scale $Li,XY\u223cLi\theta /2$. Similarly, the path

*x*measured along the magnetic flux tubes projects into the path $xXY\u2248x\theta /2$ in the

*XY*plane. If one considers the dynamics of magnetic field lines in the

*XY*-plane, the speed of propagation of Alfvénic perturbations along the

*x*-axis is $VA\theta /2$, and therefore, in this approximation, it is obvious that the growth of the outflow region Δ as seen in the

*XY*-plane does not depend on the angle

*θ*between reconnecting fluxes, i.e., magnetic reconnection proceeds independently of the value of guide field but depends only on the

*B*field and the properties of turbulence in the

_{x}*XY*-plane.

In numerical experiments described in Sec. V, the properties of driving and the value of the reconnecting component of the magnetic field were fixed in the *XY*-plane. Therefore, the magnetic field wandering for such settings can be described by Eqs. (41) and (42) with *V _{A}* arising only from

*B*and

_{x}*L*being the injection scale of turbulence in the

_{i}*XY*-plane. Correspondingly, Eq. (40) should use $MA\u2248\delta Bx/Bx$, where $\delta Bx$ is a turbulent perturbation of

*B*at the injection scale

_{x}*L*.

_{i}The LV99 theory of turbulent reconnection has provided a number of predictions very different from the contemporary theories of fast magnetic reconnection. First of all, the turbulent reconnection was identified as a generic process that takes place everywhere in the magnetized conducting fluid. A violation of flux freezing implicitly follows from that. Second, LV99 does not prescribe the given reconnection rate but provides an expression that shows that the speed of reconnection varies depending on the level of turbulence. This can explain the variability of the energy release in reconnection events observed in, e.g., solar flares. Indeed, many of the observation processes require that the rate of reconnection changes. In LV99 theory, this is controlled by the level of turbulence that can vary. For instance, the low rate of turbulent reconnection can explain how magnetic flux can be accumulated prior to a solar flare. Third, the independence on the plasma parameters makes the turbulent reconnection a generic reconnection process in the astrophysical settings.

While the earlier theories of turbulent reconnection were focused on X-point reconnection, LV99 showed that Y-type reconnection can be fast. The tearing reconnection that was developed later (see Sec. X C) presented another version of the fast Y-type reconnection. Similar to LV99, it departed from the regular structures prescribed by Sweet-Parker and Petschek reconnection models.

Being Y-type fast reconnection, LV99 necessarily involved a significant volume of magnetized fluid in the reconnection process. This volume reconnection changes the nature of the energetic particle acceleration as it is discussed in Sec. IX C.

### B. Flux freezing violation: Richardson dispersion

#### 1. Alfvén flux freezing

Since the seminal work by Alfvén (1942), the “flux freezing principle” has become a powerful tool to estimate the solutions of MHD equations in a diverse range of problems (Parker, 1979; Kulsrud, 2005). Besides its utility in plasma physics and in astrophysics, it has found many applications, for example, in explaining the low angular momentum of stars and the spiral structure of magnetic field lines in the solar wind. Despite its success in addressing diverse problems, however, there have been also fatal failures (see Burlaga *et al.*, 1982; Khabarova and Obridko, 2012; Richardson *et al.*, 2013; Eyink, 2015). For instance, assuming flux freezing, the magnetic topology in a high-conductivity magnetized fluid should not change rapidly, whereas the observations of solar flares and coronal mass ejections (CMEs) indicate otherwise. Also, in star formation, the magnetic pressure of the in-falling magnetized material should be strong enough to prevent gravitational collapse altogether, if magnetic flux-freezing held with good accuracy. Finally, if the magnetic field were frozen into the fluid in small-scale astrophysical dynamos, the tangled structure of the generated magnetic field would quench its further growth.

In order to account for these difficulties with magnetic flux-freezing, one might find it tempting to appeal to additional effects besides collisional resistivity. For instance, one might consider additional terms, such as Hall effect and electron pressure anisotropy, in the generalized Ohm's law (electron momentum equation) (Burch *et al.*, 2016; Cassak, 2016). This does not answer the question, unfortunately, how magnetic reconnection events can be fast in fluids with high rates of collisions. Indeed, both observations and theory imply fast magnetic reconnection in the solar photosphere and other environments such as collisional parts of accretion disks, which cannot be explained by appealing to collisionless mechanisms (Chae *et al.*, 2010; Singh *et al.*, 2011). The problem persists even for weakly collisional systems where an ideal MHD-type description should be valid at length-scales much larger than the ion gyroradius, above which flux-freezing is often assumed to hold. However, observations indicate that magnetic structures with length scales much larger than ion gyroradius reconnect very rapidly, e.g., at the heliospheric current sheet (HCS) in the solar wind (Gosling, 2012). It is unclear how collisionless mechanisms at scales below the ion cyclotron radius can lead to fast magnetic reconnection on much larger scales. This is the essence of the “scale problem” in magnetic reconnection (Ji *et al.*, 2019).

While the necessity of fast magnetic reconnection was widely accepted, most of the studies assumed that the fast reconnection happens for very special configurations of magnetic fields. This way the community could reconcile at least some of the observational evidence of fast reconnection and the concept of magnetic flux freezing. LV99, as discussed previously, identified magnetic reconnection as an intrinsic property of the MHD turbulent cascade, implying the violation of the flux freezing principle (see Vishniac and Lazarian, 1999, for an explicit statement about the flux freezing in turbulent fluids).

Below we provide a description of the Richardson dispersion, the phenomenon that manifests a gross violation of the flux freezing in magnetized turbulent fluids. This numerically confirmed phenomenon (see Sec. V B) provides another way of rederiving the LV99 rate of turbulent reconnection. A more rigorous approach to the flux freezing violation is discussed in Sec. IV.

#### 2. Time dependent Richardson dispersion

Another way to understand the violation of flux freezing in turbulent fluids is related to the concept of Richardson dispersion. This process initially introduced in hydrodynamic turbulence carries over to MHD turbulence. Richardson (1926) empirically analyzed the dispersion of the volcanic ashes during the eruption and formulated the law that can be derived using hydrodynamic turbulence theory. Note, however, that the corresponding theory was formulated much later by Kolmogorov (1941). Using this theory, one can predict that the separation between two particles obeys the equation $d/dt[l(t)]\u223cv(l)\u223c\alpha l1/3$, where *α* is proportional to a cubical root of the energy cascading rate. The solution of this equation $l2\u223ct3$ describes the Richardson dispersion, i.e., the mean square separation between particles increases in proportion to the cube of time. The accelerated growth of the separation is due to the fact that the larger the separation, the larger the eddies that induce the separation. In Kolmogorov turbulence, the larger eddies have larger velocity dispersion. An interesting feature of this solution is that it provides the type of fast separation even if the initial separation of particles is zero. Formally, this means the violation of Laplacian determinism. Mathematically, the above paradox is resolved by accounting to the fact that the turbulent field is not differentiable (the Kolmogorov velocity field is Holder continuous, i.e., $|v(r1)\u2212v(r2)|\u2272C|r1\u2212r2|1/3.$), and therefore, it is not unexpected that the initial value problem does not have a unique solution. In any realistic physical settings, the turbulence is damped at a nonzero scale *l _{min}*, which makes the paradox less vivid. At the scales less than

*l*, the separation of points is not stochastic, and it follows the Lyapunov growth (see Lazarian, 2006). However, at the scales larger

_{min}*l*but smaller than the injection scale of turbulence

_{min}*L*, i.e., over the inertial range of turbulence, the explosive growth of separation is still present.

_{max}In terms of motions perpendicular to the local direction of the magnetic field, the scaling of velocities in the GS95 picture corresponds to the Kolmogorov one. Therefore, we do expect the plasma particles to separate in accordance with the Richardson law over the inertial range. The situation is illustrated in Fig. 5. The magnetic field is shown to be perpendicular to the plane of the picture. A cloud of test particles undergoes an ordinary diffusion for $a\u226bLmax$, undergoing Lyapunov separation of field lines for $a<lmin$ and exhibiting Richardson dispersion over the inertial range $[Lmax,lmin]$. Obviously, the Richardson-type dynamics is not possible without magnetic field lines disentangling themselves via fast reconnection. This corresponds to the notion in LV99 that turbulent reconnection is required for turbulent magnetized fluid to preserve the fluid-type behavior. We note that if reconnection were slow, the intersecting magnetic field flux tubes would form unresolved knots. This inevitably would arrest the fluid motions, making the fluid with a dynamically important magnetic field Jello or feltlike substance. The spectrum of magnetic fields would be very modified in this case with a significant increase in energy at small scales. Only acoustic type waves can propagate in this substance and definitely no present-day numerical simulations of MHD type can describe it. This picture of slow reconnection in turbulent media contradicts both numerical simulations and observations of interstellar turbulence (see Secs. II A and II B). The turbulent reconnection, on the contrary, predicts the free evolution of eddies that produce the Kolmogorov-type cascade for eddies perpendicular to the local direction of the magnetic field as discussed in Sec. II B.

The Richardson dispersion is a numerically proven phenomenon (see Sec. V B), and it was used in ELV11 to rederive the LV99 rates of reconnection.

Consider first the Sweet-Parker reconnection. The Ohmic diffusion induces the mean square displacement across the reconnection layer in a time *t* given by

where $\lambda =c2/4\pi \sigma $ is the magnetic diffusivity. As in our previous treatment, the magnetized plasma is ejected out of the sides of the reconnection layer of length *L _{x}* at a velocity of order

*V*. Thus, the time that the lines can spend in the resistive layer is the Alfvén crossing time $tA=Lx/VA$. Therefore, the reconnected field lines are separated by a distance

_{A}where *S* is the Lundquist number. Combining Eqs. (34) and (44), one obtains the Sweet-Parker result, $vrec=VA/S$.

The difference with the turbulent case is that instead of Ohmic diffusion, the Richardson dispersion dominates (ELV11). In this case, the mean squared separation of particles follows the Richardson law, i.e., $\u27e8|x1(t)\u2212x2(t)|2\u27e9\u2248\u03f5t3$, where *t* is the time and *ϵ* is the energy cascading rate. For sub-Alfvénic turbulence, $\u03f5\u2248uL4/(VALi)$ (see LV99), and therefore, one can write

in the limit of $L<Li$. Similar considerations allowed ELV11 to rederive the LV99 expression for $L>Li$, which differs from Eq. (46) by the change in the power 1/2 to $\u22121/2$.

#### 3. Field line wandering

Magnetic fields produce a new effect compared to the Richardson dispersion in hydrodynamics. In magnetized fluids, one can trace magnetic line separation. This effect has been discussed for decades in the cosmic ray literature as the cause of the perpendicular diffusion of cosmic rays in the galactic magnetic field (Jokipii, 1973). It is worth mentioning that the correct quantitative description of the effect was obtained only in LV99. Earlier papers assumed the diffusive separation of magnetic field lines, which is the effect taking place only if magnetic field lines are separated by more than the turbulence injection scale *L*. It was shown in LV99 that at scales less than *L*, the magnetic field lines separate superdiffusively. This changes a lot both the propagation and the acceleration of cosmic rays (see Lazarian and Yan, 2014).

Dealing with magnetic reconnection in Sec. III, we have discussed the separation of magnetic field lines. Figure 6 illustrates the spread of magnetic field lines in the perpendicular direction in the situation when magnetic field lines are traced by particles moving along them. This figure confirms the prediction for the magnetic field line separation in Eq. (38), i.e., the squared separation between the magnetic field lines increases in proportion to the cube of the distance measured along the magnetic field lines (see Fig. 6).

The other important prediction in Eq. (38) is the change in the square of the separation between the magnetic field lines, which increases with the forth power of the Alfvén Mach number, i.e., $\u223cMA4$. The successful testing of this prediction is demonstrated in Fig. 7, which is taken from the study by Xu and Yan (2013). The asterisks correspond to the calculations obtained using the parallel particle velocity. This way of calculating the separation provides a better description of magnetic field wandering. The fitted dependence is $MA\u22123.8$, which is very different from the $MA2$ dependence traditionally accepted in the cosmic ray literature.

We note that the magnetic field tracing with particles is not the only way to numerically explore the magnetic field separation. In fact, the analytical predictions of magnetic field line superdiffusion were also explored numerically by the direct tracing of magnetic field lines first with low resolution by Lazarian *et al.* (2004) and with higher resolution in Beresnyak (2013).

Taken together, these numerical testing provides a solid confirmation of the predictions given by Eq. (38). It is significant that the turbulent magnetic field wandering induced by strong MHD turbulence provides the separation of magnetic field lines $\delta lines2\u223cs3$, where the distance *s* along the instantaneous realization of a magnetic field line. Thus, for sufficiently large *s*, all parts of the volume of magnetized fluid get connected. In other words, the entire volume becomes accessible to particles moving along magnetic field lines. This is important for a wide variety of astrophysical processes, e.g., the cosmic ray propagation and acceleration (Yan and Lazarian, 2008; Lazarian and Yan, 2014) and heat propagation (Narayan and Medvedev, 2001; Lazarian, 2006).

We would like to stress again that magnetic wandering cannot be understood using the notion of magnetic field lines preserving their identify in turbulent flows. Magnetic field lines wandering in the turbulent magnetized flow constantly reconnect, inducing the exchange of plasmas and magnetic field. In view of this, it is interesting to recall the considerations on the nontrivial nature of astrophysical magnetic field lines that can be found in a prophetic book by Parker (1979). There it was mentioned that a magnetic field line may be traced through a significant part of the galaxy. We add that this is only possible through constant reforming of the magnetic field lines via reconnection.

We note that there is a formal analogy between the Richardson diffusion that predicts that the mean squared separation between particles $\delta particles2$ that grow with time as *t*^{3} and the mean squared separation between the magnetic field lines $\delta lines2$ that grow with the distance *s* as *s*^{3}. Therefore, rather than using a new term to describe the latter dependence theoretically predicted in LV99, it was decided in ELV11 to term this new type of dispersion as “the Richardson diffusion in space,” while keeping the term “Richardson dispersion” for the original time-dependent process empirically discovered in Richardson (1926).

### C. Applicability of the LV99 approach

The derivation of the LV99 reconnection rate assumes that all magnetic energy is converted to the kinetic energy of the outflow. This may not be true when *M _{A}* is large. In this case, the energy dissipation in the reconnection layer becomes non-negligible, and the decrease in the outflow velocity is expected (see ELV11).

The corresponding effect of turbulent dissipation can be estimated from steady-state energy balance in the reconnection layer

where the kinetic energy of the outflow is balanced against the magnetic energy transported into the layer minus the energy dissipated by turbulence. The turbulent dissipation is given by $\epsilon =uL4/VALi$ for sub-Alfvénic turbulence (Kraichnan, 1965). Dividing (47) by $\Delta =Lxvrec/vout$, one obtains

The corresponding solutions are easy to get by introducing the ratios $f=vout/VA$ and $r=2MA4(Lx/Li)$ which are the measures of the outflow speed and the energy dissipated by turbulence in units of the available magnetic energy. This gives

When *r *=* *0, the only solution of (49) with *f *>* *0 is $f=1,$ which is the LV99 estimate $vout=VA$ for $MA\u226a1$. For the larger values of *r*, $f\u22431\u2212(r/2)$, in agreement with the formula $f=(1\u2212r)1/2$ that follows from Eq. (65) in ELV11. This implies a decrease in *v _{out}* as compared to

*V*. Note that formula (49) can be used when

_{A}*r*is not too large. The largest value of

*r*for which a positive, real

*f*exists is $rmax=2/27\u22480.385$. For

*r*, one gets $fmin=1/3\u22480.577$. This implies that the LV99 approach is limited to sufficiently small

_{max}*M*. For $Lx\u2243Li$, one may consider values of

_{A}*M*up to 0.662. Therefore, the original LV99 approach is accurate for $MA\u22721$.

_{A}While the accuracy of LV99 expressions can be decreased due to the additional dissipation of the turbulent energy, the predicted phenomenon of fast reconnection arising from the turbulence persists for any *M _{A}*. In fact, one can see that the effect of the reduced outflow velocity “increases” the reconnection rate. The reason is that field-lines now spend a time $Lx/vout$ exiting from the reconnection layer, greater than that assumed in LV99 by a factor of $1/f$. This implies a thicker reconnection layer Δ due to the longer time-interval of Richardson dispersion in the layer, greater than the LV99 estimate by a factor of $(1/f)3/2$. The net reconnection speed $vrec=vout\Delta /Lx$ is thus larger by a factor of $(1/f)1/2$. The increased width Δ more than offsets the reduced outflow velocity

*v*. At even larger

_{out}*M*, magnetic fields at a large scale are weak that they are easily bent and twisted by the fluid turbulent motions. However, at the scales smaller than

_{A}*l*, we expect the original LV99 treatment to be accurate again.

_{A}### D. LV99 model: Summary of facts

The model of turbulent reconnection suggested in LV99 departed from the earlier attempts to explain astrophysical magnetic reconnection as something that required special conditions. On the contrary, it suggested that reconnection takes place everywhere in realistically turbulent fluids. Moreover, the transition to turbulent reconnection is inevitable as the 3D outflow from the reconnection region gets turbulent for sufficiently large Lundquist numbers.

The model identified turbulent magnetic field wandering as the process that opens up the reconnection region outflow making reconnection fast, i.e., independent of fluid resistivity or microscopic plasma effects. The derived dependences for the field wandering were successfully tested with the numerical simulations.

The model provided analytical expressions relating the reconnection rate with the level of turbulence. It also provided the relationships between the thickness of the reconnection layer and the level of turbulence. This explains that the rate of magnetic reconnection may vary significantly and also provides ways for quantitative numerical and observational testing of the model.

## IV. VIOLATION OF FLUX-FREEZING, RENORMALIZATION, AND THE SCALE PROBLEM

In Sec. III B, we discussed the issue of flux freezing on the basis of the simplified treatment of Richardson dispersion, treating the latter more like an empirical fact and using the corresponding scaling relations to rederive the LV99 reconnection rates. However, the phenomenon of Richardson dispersion presents a much more fundamental property of magnetic fields in turbulent fluids. Below we discuss the violation of flux freezing on the basis of a more rigorous theoretical treatment.

### A. Spontaneous stochasticity

The conventional flux-freezing principle has not only problems with explaining observations. It also encounters serious theoretical difficulties. This can be exemplified by naively extending the concept to turbulent MHD flows with both small resistivity and small viscosity. To understand the subtlety of such singular limits, consider the problem of incompressible fluid turbulence. As is well known, the kinetic energy dissipation rate, in a fluid with viscosity *ν* and velocity field $u$, does not vanish in the limit of vanishing viscosity [e.g., see Sreenivasan (1984, 1998)]. On the contrary, $lim\nu \u21920\nu |\u2207u|2>0$, which implies a “dissipative anomaly.” This in turn requires that velocity gradients must blow up, with $\u2207u\u2192\u221e$ in the limit $\nu \u21920$. In fact, as Onsager pointed out in a prescient work (see Sec. IV B below), the velocity field in a turbulent fluid must develop Hölder-type singularities in the limit as viscosity tends to zero (Onsager, 1949; Eyink and Sreenivasan, 2006; Eyink, 2018). These singularities render all hydrodynamic equations containing velocity gradients naively undefined. Extending Onsager's analysis to MHD plasmas, turbulent solutions are found to suffer several dissipative anomalies that require diverging gradients of all MHD fields, e.g., velocity field, magnetic field, and density (see, e.g., Mininni and Pouquet, 2009; Caflisch *et al.*, 1997; Eyink and Drivas, 2018). For example, magnetic-field gradients blow up in the limit as magnetic diffusivity *η* tends to zero, with $\u2207B\u2192\u221e$ as $\eta \u21920$. Because velocity fields in particular become nonsmooth and singular as $\nu \u21920$, the magnetic field lines cannot freeze into the fluid in a standard deterministic sense. Instead, a nondeterministic frozen-in property known as “stochastic flux freezing” (Sec. IV C) has been established by Eyink (2009) and (2011), which is intrinsically random in nature due to the roughness of the advecting velocity field (Eyink, 2010) and is associated with the “spontaneous stochasticity” of Lagrangian particle trajectories (Bernard *et al.*, 1998; Eyink *et al.*, 2011, 2013).

The latter concept is crucial and must be briefly explained. In essence, spontaneous stochasticity means that stochastic Lagrangian particle trajectories remain nonunique and indeterministic (random) for turbulent flows even in the limit of vanishing viscosity and molecular diffusivity (Bernard *et al.*, 1998). This phenomenon is directly related to the Richardson superballistic behavior (Richardson, 1926) (see Sec. III B). In a turbulent flow, therefore, infinitely many magnetic field lines are advected to each point in space even in the ideal limit and must be averaged to obtain the resultant magnetic field (Eyink, 2011; Eyink *et al.*, 2013). Due to Richardson dispersion, in fact, even a tiny breakdown of magnetic flux freezing at small scales will be explosively amplified into larger scales in the turbulent inertial range (Eyink, 2015, see also Sec. IV D below). The rates of such rapid amplification are nonvanishing in the ideal limit (Eyink, 2011; Eyink *et al.*, 2011). This breakdown of conventional magnetic flux freezing has already been verified in numerical simulations of resistive MHD turbulence at high conductivity (e.g., Eyink *et al.*, 2013).

To appreciate the importance of (asymptotic) singularities in MHD, we may note that without the differentiability of fields, even fundamental concepts such as the notion of field lines needs to be reconsidered. Taking magnetic field $B(x)$ as an example, its “field line” through point $x$ is conventionally defined as the integral curve $\xi (s;x)$ whose tangent vector at any point is parallel to the vector field at that point, so that, parameterized with arc length *s*

If the magnetic field satisfies the following inequality for some constant *C *>* *0

with *h *=* *1 (Lipschitz condition), then there exists a unique field-line $\xi (s;x)$ passing through point $x$; however, if this condition holds merely with $0<h<1$ (Hölder continuity), then there are generally infinitely many such integral curves (see Jafari and Vishniac, 2019a, 2019b; Jafari *et al.*, 2019b). If the Lipschitz condition is not satisfied and hence the uniqueness theorem cannot be applied for a given magnetic field, what would a magnetic field line mean after all? When resistivity is tiny but nonzero, then $B$ is differentiable, and indeed, Eq. (50) has unique solutions. However, these solutions exhibit such extreme sensitivity to initial conditions that the distance $|\xi (s;x\u2032)\u2212\xi (s;x)|$ between field lines through neighboring points becomes “independent of” $x\u2032\u2212x$ at distances $s\u226b|x\u2032\u2212x|$ along the lines (LV99) (we provide more detail on this effect in Sec. III B). Such explosive separation of magnetic field lines is not observed if the field lines exhibit standard “deterministic chaos,” where the memory of the initial separation is preserved. Exactly analogous difficulties plague the concepts of “Lagrangian fluid particle” or “fluid trajectory,” which are ill-conditioned and asymptotically nonunique if the fluid velocity $v(x,t)$ becomes spatially non-Lipschitz as $\nu \u21920.$

The singularities associated with dissipative anomalies lead to subtle effects whose consistent description requires great care; otherwise, contradictory conclusions are obtained. Similar difficulties were encountered already many decades ago in quantum field theory and condensed matter physics, which could be alleviated by mathematical methods collectively known as regularization and renormalization. As we review in Sec. IV B, such methods can be applied also to plasma turbulence by exploiting a spatial coarse-graining operation, which removes divergences and leads to meaningful equations. Physically speaking, “coarse-graining” simply means averaging quantities in a volume in space. Indeed, what is measured in the laboratory as magnetic or velocity fields at a point in space and time is actually always an average over a small volume around that point during a time interval. Paradoxes are avoided by developing a theory based upon such quantities that can be measured in principle rather than idealized pointwise fields that are strictly unobservable. Such renormalization introduces an “observation scale” or “resolution scale” which is entirely arbitrary. However, obviously, no objective physical fact can depend upon this arbitrary scale. In the language of quantum field theory and statistical physics, this principle is called “renormalization-group invariance” (Gross, 1976).

### B. Renormalization of MHD equations and singularities

We illustrate Onsager's renormalization-type arguments in one of the simplest models with realistic astrophysical applications and the standard MHD equations for a compressible, magnetized, single-component fluid with mass density $\rho ,$ velocity $v,$ magnetic field $B,$ and internal energy density *u*

Here, $I$ is the identity tensor and $p=p(\rho ,u)$ is the thermodynamic equation of state for the pressure and

with shear viscosity $\mu ,$ bulk viscosity $\zeta ,$ thermal conductivity $\kappa ,$ and electrical conductivity *σ* (or magnetic diffusivity $\eta =c2/4\pi \sigma $). See Landau *et al.* (2013), Chap. VIII, Secs. 65 and 66.

The local balance of mechanical energy for the above model is given by the following equation:

Just as for the case of incompressible, neutral fluids studied by Onsager, turbulent solutions are characterized by “dissipative anomalies” (see, e.g., Mininni and Pouquet, 2009) in which

Similar anomalous dissipation must then appear also in the balances of internal energy and entropy. Such anomalies obviously require diverging space-gradients $|\u2207v|,\u2009|\u2207B|\u2192\u221e$ for $\mu ,\u2009\zeta ,\u2009\kappa ,\u2009\eta \u21920$ or a “violet catastrophe” in the words of Onsager (1945).

To obtain a well-defined dynamical description in this singular ideal limit, one must regularize these UV divergences. For example, one may employ a coarse-graining or “block-spin” regularization, defined as

and likewise for $v,$$\rho ,$*u*, etc., with $G\u2113(r)=\u2113\u22123G(r/\u2113)$ for some positive, smooth, rapidly decaying kernel $G(\rho )$ normalized as $\u222bd3\rho \u2009G(\rho )=1.$ Coarse-grained gradients such as $\u2207B\xaf\u2113(x,t)$ necessarily remain finite as $\mu ,\zeta ,\kappa ,\eta \u21920,$ and their UV divergences are removed. It should be emphasized again that such a coarse-graining is purely passive and corresponds to observing the flow “with spectacles off” so that eddies of size $<\u2113$ cannot be resolved. The basic principle of renormalization-group invariance is that this regularization scale $\u2113$ is completely arbitrary, and no objective physics can depend upon it. Because the coarse-graining operation commutes with space-time derivatives, the effective equations of motion for eddies of size $>\u2113$ are exactly the same as the original MHD equations with just an additional overbar $()\xaf\u2113$.

However, even when there are dissipative anomalies in the microscopic or “bare” balance equations for conserved quantities, the dissipative transport terms such as $\u2207\xd7(\eta \u2207\xd7B)\xaf\u2113$ in the coarse-grained MHD equations can be shown to vanish for $\mu ,\zeta ,\kappa ,\eta \u226a1$ at fixed $\u2113$ or for $\u2113\u226b1$ at fixed $\mu ,\zeta ,\kappa ,\eta $. These dissipative terms are thus “irrelevant” in the renormalization-group sense. It follows, when $\mu ,\zeta ,\kappa ,\eta \u226a1,$ that there exists an “inertial-range” of length-scales $\u2113$ at which the ideal MHD equations hold in the “coarse-grained sense,” i.e.,

Taking the limit $\mu ,\zeta ,\kappa ,\eta \u21920$, these coarse-grained equations hold for all $\u2113>0,$ which is equivalent to the standard mathematical notion of a “weak”/“distributional” solution or what Onsager called a “more general description” than standard differential equations.

The validity of the ideal MHD equations in this generalized sense in not sufficient, however, to guarantee that naive balances of mechanical energy, internal energy, and entropy will hold. To show this clearly, the above ideal equation can be cosmetically rewritten by introducing the “density-weighted Favre-average”

and “pth-order cumulants” $\tau \xaf\u2113(f1,f2,\u2026,fp)$ and $\tau \u0303\u2113(f1,f2,\u2026,fp)$ for the two types of coarse-graining averages. With these definitions, the effective equations of motion for scales $>\u2113$ become

and

Here, we have introduced the “turbulent stress tensor”

“turbulent electromotive force (emf)”

“turbulent internal energy”

and “turbulent enthalpy transport”

with $h=u+p$ the enthalpy per volume.

Although Eqs. (55)–(58) are exactly equivalent to the coarse-grained ideal MHD equations (51)–(54), they appear “nonideal.” Integrating out the small-scale eddies has led to “renormalized” equations with new terms $T\u2113,$ $\epsilon \u2113,$ $\delta u\xaf\u2113,$ and $J\u2113h$, which would not appear in naive ideal MHD equations for the resolved fields $\rho \xaf\u2113,$ $v\u0303\u2113,$ $B\xaf\u2113,$ and $u\xaf\u2113$ at length-scales $>\u2113.$ The precise condition for the validity of the ideal equations (51)–(54) at length-scale $\u2113$ is that the standard dissipative transport terms are negligible compared with these new terms. For example, the ideal magnetic induction equation (57) is valid precisely when

One of the Onsager key insights was that contributions from coarse-graining such as $\epsilon \u2113$ can be expressed entirely in terms of “space-increments” (see Eyink and Aluie, 2006; Eyink and Drivas, 2018), such as

It follows that

In that case, the condition (60) becomes

which precisely defines the inertial-range of scales $\u2113$ for which the ideal induction equation is valid. This condition can be restated as $\u2113\u226b\u2113\eta $, where $\u2113\eta $ is the resistive dissipation length specified by $\u2113\eta \delta v(\u2113\eta )\u2243\eta .$ Similar conditions guarantee the validity of all of the other ideal equations as well (see Eyink and Drivas, 2018).

From the coarse-grained ideal MHD equations rewritten in this manner, one obtains a balance equation for “resolved mechanical-energy”

where the renormalization has introduced new apparently nonideal terms. In particular, the right side of this balance contains the “turbulent energy flux”

which represents the nonlinear cascade of mechanical energy to unresolved scales $<\u2113$. The energy lost from resolved mechanical modes reappears as total unresolved energy $u\xaf\u2113*=u\xaf\u2113+\delta u\xaf\u2113,$ where $u\xaf\u2113$ is the thermal kinetic energy of microscopic particles and $\delta u\xaf\u2113$ represents mechanical energy in unresolved turbulent eddies at scales $<\u2113.$ The balance equation for this total unresolved energy or “intrinsic internal energy” can be obtained by combining Eq. (62) with coarse-grained energy conservation to give

where $J\u2113u=J\u2113h+p\xaf\u2009\tau \xaf(\rho ,v)/\rho \xaf$ and $Q\u2113flux$ is the same term that appears in Eq. (62) but with the opposite sign.

If mechanical energy is dissipated in the limit $\mu ,\zeta ,\kappa ,\eta \u21920,$ then such dissipation will be apparent even to a “myopic” observer who sees only eddies $>\u2113$, and the dissipation rate must become independent of $\u2113$ for $\u2113\u21920.$ From the estimates of the form given by Eq. (61), it follows that $Q\u2113flux$ can be expressed entirely in terms of space-increments as discussed above, and therefore, $Q\u2113flux$ can be nonvanishing as $\u2113\u21920$ only if the solution fields develop a critical degree of singularity. For example, suppose that for all space-time points $(x,t)$

and similarly for $\delta v(r)$ with some exponent $hv$. Since

this term, which represents the cascade of magnetic energy, can remain nonvanishing as $\u2113\u21920$ only if

Similarly, the first term in $Q\u2113flux$ which represents the cascade of kinetic energy can remain nonvanishing only if $3hv\u22641$ and the third term, which generalizes to MHD the “baropycnal work” of Aluie (2013), requires at least one of the inequalities $h\rho +hv+hB\u22641$ or $h\rho +hv+hp\u22641$ in order to be nonvanishing. These conditions taken altogether imply that $hv<1$ as $\mu ,\zeta ,\kappa ,\eta \u21920$, and thus, the velocity field “cannot” remain Lipschitz in the limit.

It should be clear from the previous discussion that the balances of other quantities besides mechanical energy can suffer from dissipative anomalies. In particular, the induction equation (57) contains the turbulent electromotive force $\epsilon \u2113$ which not only drives the cascade of magnetic energy but also vitiates the freezing of the resolved magnetic field $B\xaf\u2113$ to the resolved velocity $v\xaf\u2113.$ If instead the magnetic field is referred to the Favre-averaged velocity $v\u0303\u2113,$ then the violation of flux-freezing is due to the combination $\epsilon \u2113\u2212(1/\rho \xaf\u2113)\tau \xaf\u2113(\rho ,v)\xd7B\xaf\u2113.$ Because generally $\epsilon \u2113\xb7B\xaf\u2113\u22600,$ it is not possible to express the turbulent emf as $\epsilon \u2113=\Delta v\u2113\xd7B\xaf\u2113$ for any choice of vector field $\Delta v\u2113,$ which would allow one to regard lines of $B\xaf\u2113$ as moving with a velocity $v\xaf\u2113*=v\xaf\u2113+\Delta v\u2113.$ Thus, it is not true, as is often stated that “because the ideal Ohm's law holds in the turbulent inertial-range, magnetic field-lines are frozen-in at those scales.” The error in this commonplace statement is that no precise operational meaning is ever given to the claim of “flux-freezing in the inertial-range of scales,” so that an experimentalist could in principle test it. As known from the exalted areas of physics such as relativity and quantum mechanics, the failure to formulate such claims in operationally meaningful terms can lead to paradoxes and inconsistencies. If one adopts the natural interpretation that the magnetic field $B\xaf\u2113$ of the eddies $>\u2113$ be frozen-in to the velocity fields $v\xaf\u2113$ or $v\u0303\u2113$ at those same scales, then the statements are certainly false.

We discuss this later in Sec. X D in relation to the concept of reconnection-mediated turbulence which implies that magnetic reconnection takes place only at some small scales where tearing gets important. It is obvious from the discussion above that magnetic reconnection takes place over the entire inertial range rather than some chosen small scales where dissipative processes take place. The reconnection at all scales of the turbulent cascade is necessary to explain the Richardson dispersion (see Sec. III B) and the numerical simulations demonstrating the gross violation of the classical flux freezing (see Sec. V B). In fact, this idea is implicitly at the very core of the LV99 theory.

A clear-minded skeptic might argue that such violations of flux-freezing in the turbulent inertial-range are mere artifacts of the coarse-graining and not “real.” In particular, coarse-graining could be applied even to a laminar ideal MHD flow and would lead to apparent violations of flux-freezing at length-scales $\u2113\u226bL\u2207,$ the “gradient-length” over which the MHD solution fields sensibly vary. If such violations are real, then, invoking renormalization-group invariance, they must be observed for “all” resolution scales. An experimentalist who measured the velocity and magnetic fields in a laminar ideal MHD flow would find that, as the length-scales resolved by measurements were refined more, standard flux-freezing would hold better. However, in a turbulent MHD flow, deviations from flux-freezing predictions would persist as the experimentalist improved measurement resolution down to smaller length-scales $\u2113$ within the inertial-range, even though ideal MHD holds there in the “coarse-grained sense.” Eventually, the experimentalist would resolve down to the resistive scale $\u2113\sigma $ and would find the same violations of flux-freezing predictions but now due to Ohmic electric fields. We note that the assumption that the turbulence is not damped up to the Ohmic diffusion scale is not a necessary one. In Sec. VIII A, we show that fast reconnection is possible in the situation when viscosity is larger than resistivity.

Considering the infinite-Reynolds limit $\mu ,\zeta ,\kappa ,\eta \u21920$, the ideal MHD equations would be valid at all scales $\u2113$ and the violations of flux-freezing would be due entirely to the turbulent emf $\epsilon \u2113.$ Eyink and Aluie (2006) have shown that magnetic flux conservation can be violated at an instant of time for an arbitrary small length scale $\u2113$ in the absence of any microscopic plasma nonideality only if at least one of the following necessary conditions is satisfied: (i) nonrectifiability of advected loops, (ii) unbounded velocity or magnetic fields, or (iii) existence of singular current sheets and vortex sheets that intersect in sets of large enough dimension. Both conditions (i) and (iii) can be expected to hold in MHD turbulence.

As we now review, this conclusion that flux-freezing can be violated in turbulent solutions of ideal MHD is corroborated by two entirely independent arguments. The first of these appeals to the phenomenon of spontaneous stochasticity of Lagrangian fluid particle trajectories for non-Lipschitz velocity fields. The second argument invokes the spontaneous stochasticity of the magnetic field-lines themselves and the concept of “slippage” of field-lines relative to the plasma fluid.

### C. Conventional and stochastic flux-freezing

The flux-freezing properties of smooth, laminar solutions of the ideal MHD equations are well known. One can write the induction equation for magnetic field $B$ as

where $Dt\u2261\u2202t+u.\u2207$. Assuming that the limiting solution for $\eta \u21920$ remains smooth, this equation may be combined with mass balance $Dt\rho +\rho (\u2207\xb7v)=0$ to yield

which succinctly expresses flux-freezing. The latter equation is solved explicitly by the Lundquist formula

where $X(a,t)$ is the position at time *t* of the Lagrangian fluid particle that was at point $a$ at time $t0$ so that

The mass-density ratio is represented by the determinant $det\u2009(\u2207aX(a,t))|X(a,t)=x=\rho (a,t0)/\rho (x,t).$

It is less widely known that analogous “stochastic flux-freezing” properties hold even if magnetic diffusivity is nonvanishing or $\eta >0.$ In that case, one must consider stochastic Lagrangian trajectories that satisfy

where $W\u0303(t)$ is a 3D Brownian motion that represents the effects of magnetic diffusivity. The exact solution of the magnetic induction equation is then represented by the “stochastic Lundquist formula”

where the overbar $(\xb7)\xaf$ represents an average over the Brownian motion $W\u0303(t)$ or, equivalently, a path-integral over the stochastic Lagrangian trajectories $X\u0303(t)$ themselves (see Eyink, 2009, 2011). If all solution fields remain smooth in the limit $\eta \u21920,$ then the ensemble of stochastic Lagrangian trajectories $X\u0303(a,t)$ collapses to the unique deterministic Lagrangian trajectory $X(a,t)$ and conventional flux-freezing is recovered in the limit.

This is not the case for turbulent solutions of ideal MHD obtained in the joint limit $\mu ,\zeta ,\kappa ,\eta \u21920.$ As we have observed in Sec. IV B, the velocity field $v(x,t)$ cannot remain uniformly Lipschitz if there is anomalous energy dissipation in that limit. In that case, there is no longer a unique Lagrangian trajectory for each fluid element in the infinite-Reynolds limit. Therefore, one may wonder with such an infinite number of limiting particle trajectories, which trajectory will be selected in the joint limit $\mu ,\zeta ,\kappa ,\eta \u21920$? The surprising answer is that Lagrangian trajectories, under such conditions, can remain random—the phenomenon of spontaneous stochasticity (Bernard *et al.*, 1998; Chaves *et al.*, 2003). This effect resembles “spontaneous symmetry-breaking” in condensed matter physics and quantum field theory, where below a critical temperature (such as the Curie temperature for a magnet), the mean order parameter remains nonvanishing even when the external symmetry-breaking field is taken to zero. Spontaneous stochasticity differs greatly from the usual randomness in turbulence theory associated with an ensemble of velocity fields because the Lagrangian fluid trajectories remain nonunique and stochastic for a “fixed” velocity field $v(x,t)$ and fixed particle position $a$ at the initial time $t0.$ The Laplacian determinism expected for classical dynamics breaks down because of the explosive separation of particles undergoing turbulent Richardson dispersion. This remarkable prediction is supported by direct empirical evidence for turbulent MHD flows (for details, see Eyink, 2011; Eyink *et al.*, 2013) (see also Sec. V B).

As a consequence of this surprising remnant or persistent randomness of Lagrangian trajectories, the stochastic flux-freezing relation given by Eq. (65) for nonideal MHD solutions does “not” reduce to conventional flux-freezing in the limit $\mu ,\zeta ,\kappa ,\eta \u21920.$ Instead, flux-freezing in such regimes remains valid in only a statistical sense associated with the spontaneous stochasticity of the fluid trajectories [Eyink, 2011; see also Jafari and Vishniac (2018a)]. Magnetic flux conservation in MHD turbulence at large kinetic and magnetic Reynolds numbers thus holds neither in the conventional sense nor is entirely broken. In ELV11, this stochastic flux-freezing relation was employed, together with GS95 predictions for turbulent Richardson dispersion, to rederive the predictions of LV99 for magnetic reconnection rates in the presence of MHD plasma turbulence. The detailed predictions depend upon the phenomenological theory of MHD turbulence which is adopted, but the qualitative features are independent of such assumptions and, in particular, the prediction that reconnection rates are independent of the precise values of microscopic resistivity or, similarly, of any effective “anomalous resistivity” arising from any other microscopic plasma processes.

### D. Field-line stochasticity and slippage

The original arguments of LV99 rested instead upon the random wandering of magnetic field-lines $\xi (s;x)$ with increasing arc length *s* as they move away from their initial position $x.$ Anticipating the concept of spontaneous stochasticity, LV99 realized that field-lines would disperse explosively so that separations between neighboring field-lines $\xi (s;x),$ $\xi (s;x\u2032)$ would become independent of the initial separation $|x\u2032\u2212x|$ with increasing *s*. This phenomenon allows field-lines from widely separated space regions to wander into the same thin current sheet with thickness governed by $\eta ,$ even in the limit as $\eta ,$ $\nu \u21920.$ In this way, LV99 argued that violations of flux-freezing by tiny amounts of resistivity can be enhanced to macroscopic scales. Boozer (2018) has emphasized that similar wandering due to the deterministic chaos of field-lines can lead to reconnection rates nearly independent of resistivity, but ubiquitous plasma turbulence can produce strictly fast reconnection with no dependence on resistivity whatsoever.

We may note, however, a formal inconsistency in the argument of LV99 in the limit of vanishing resistivity. On the one hand, they used the GS95 phenomenology of MHD turbulence to estimate the mean square dispersion of field-lines as proportional to $\epsilon (s/vA)3,$ where *ε* is the energy dissipation per mass and *v _{A}* is the Alfvén speed. On the other hand, GS95 theory predicts that the Ohmic electric field $R=(\eta /c)\u2207\xd7B$ has a magnitude $R\u223c(\eta /c)Brms/(\eta 3/\epsilon )1/4\u221d\eta 1/4\u21920$ as $\eta \u21920.$ Thus, Ohm's law $E+(v/c)\xd7B=R$ is predicted by GS95 phenomenology to become “pointwise ideal” in the limit $\eta ,$ $\nu \u21920,$ with not only the width of current sheets but also the magnitude of the nonideal Ohmic electric field vanishing in that limit. How can magnetic reconnection persist in the zero-resistivity limit when nonideal electric fields vanish uniformly? One obvious answer is that GS95 theory ignores small-scale intermittency which can produce current sheets with thickness $\u221d\eta $ and thus Ohmic fields that persist as $\eta \u21920$. However, a deeper solution to this puzzle lies in the recognition that magnetic reconnection is possible even at points where nonideal electric fields tend to zero!

A complete analysis of this problem was presented by Eyink (2015) who studied the rate at which magnetic-field lines “slip” through a nonideal plasma or, more accurately, the rate at which magnetic connections between plasma elements change due to nonideality. As argued already by Axford (1984), such a change in magnetic connections is the basal effect that underlies all other reconnection phenomena (diffusion of plasma through separatrices, topology change, conversion of magnetic energy to kinetic energy, etc.) This effect was quantified in Eyink (2015) by the “slip velocity” $\Delta w\u22a5(s;x)$, which is the relative velocity between the magnetic field line $\xi (s;x)$ and the fluid velocity $v$ at a distance *s* in the arc length away from the plasma element $x$ which that particular field-line threads. A simple equation was obtained in Eyink (2015) for the growth of slip-velocity $\Delta w\u22a5(s;x)$ in arc length *s* along the field-line

where $B\u0302(x,t)=B(x,t)/|B(x,t)|$ denotes the magnetic field direction and subscript “$\u22a5$” the vector component perpendicular to $B.$ Because this linear ODE can have solutions $\Delta w\u22a5(s;x)\u22600$ only if the righthand side is nonzero, the above expression indicates that flux-freezing holds if and only if $(\u2207\xd7R)\u22a5=0$. In fact, the relation $B\u0302\xd7(\u2207\xd7R)=0$ has long been known as the general condition for flux freezing (Newcomb, 1958). The quantity $\Sigma =\u2212c(\u2207\xd7R)\u22a5/|B|$ represents the rate of development of slip-velocity per unit length of field-line and was dubbed in Eyink (2015) the “slip-velocity source.” This expression has a topological interpretation as well, arising from the relation $\u2202tB\u0302=(\u2202tB)\u22a5/|B|$ (Jafari and Vishniac, 2019a, 2019b).

The key point in the resolution of the “vanishing Ohmic electric-field puzzle” is that, even if nonideal electric fields $R\u21920$, the corresponding slip-velocity source $\Sigma $ can remain nonvanishing or even diverge to infinity! For example, within the GS95 phenomenology of MHD turbulence, $\Sigma \u223c\eta /(\eta 3/\epsilon )1/2\u221d\eta \u22121/2\u2192\u221e$ as $\eta \u21920.$ Thus, the original argument of LV99 is self-consistent if GS95 phenomenology is employed throughout. More importantly, the above arguments show that the most intense current sheets due to small-scale intermittency are not required to produce flux-freezing breakdown and that field-lines even far from those strong sheets will slip through the background turbulence at very sizable rates. The slip velocities will certainly be the greatest at the most intense current sheets but, because of their rarity, more slippage may occur over time due to the weaker turbulence background. Eyink (2015) argued that such turbulent line-slippage accounts for the observed deviations from the Parker spiral model in the solar wind. This line-slippage is the same phenomenon which was termed “reconnection diffusion” introduced in Lazarian (2005) and numerically tested in Santos-Lima *et al.* (2010). This process was invoked to explain the removal of magnetic flux in star-formation (see more in Sec. IX B).

One may roughly translate the quantitative arguments of Eyink (2015), discussed above, into a more qualitative statement by saying that if the magnetic field is not perfectly frozen into the turbulent fluid, i.e., in the presence of stochastic flux freezing, the field should slip through the fluid. The differential equation governing this field-fluid slippage is given by Eq. (66) with the source term $(\u2207\xd7B)\u22a5/B$. Interestingly, it turns out that this term also acts as a source in the evolution equation of the unit vector $B\u0302$, which in turn implies that $(\u2207\xd7B)\u22a5/B$ has a topological interpretation. In fact, it is simple calculus to see that the derivative of the unit vector $B\u0302=B/|B|$ is given by

where $(.)\u22a5$ denotes the perpendicular component with respect to $B$. The evolution equation of $B\u0302$, therefore, is easily obtained using the induction equation

Jafari and Vishniac (2019b) used the relationship between the source term $(\u2207\xd7B)\u22a5/B$ and the field's topology, implied by Eq. (67) above, in order to define a measure for the spatial complexity, or self-entanglement, of the magnetic field $B$ as

where $B\u0302l=B\xafl/Bl$ and $B\u0302L=B\xafL/BL$ with $B\xafl$ and $B\xafL$ as the coarse-grained, or renormalized, components at arbitrary scales *l* and *L*. If we take *L* of order of the system size and $l\u226aL$ in the inertial range, $B\xafl(x,t)$ will be the average field of a fluid parcel of size $\u223cl$, while $B\xafL(x,t)$ can be thought of as the average field of a fluid parcel of scale $L\u226bl$. The former plays the role of a local field, whereas the latter acts as a global field, and therefore, $1\u2212B\u0302l\u22c5B\u0302L$ is in fact a local measure of the deviation of $B\xafl$ from $B\xafL$. The RMS value of this quantity measures the global spatial complexity of the field [for a mathematical treatment of magnetic topology and complexity, see Jafari and Vishniac (2019a)]. Magnetic reconnection, field-fluid slippage, and magnetic topology change can be studied statistically using the function *S*(*t*) and its time derivative $\u2202tS(t)$ (Jafari and Vishniac, 2019b; Jafari *et al.*, 2019b). It also turns out that *S*(*t*) is related to the Richardson diffusion of the magnetic field in turbulence; see Jafari *et al.* (2019a).

### E. Irrelevance of small scale physics in large scale turbulent reconnection

Although the previous arguments were presented in the context of resistive MHD, they are not restricted to plasmas with high collision rates but apply with equal force to inertial-range turbulence in collisionless plasmas, such as in the solar wind. The main difference is that in the “generalized Ohm's law” of a collisionless plasma

the nonideal electric field $R$ contains many additional terms from the Hall effect, electron pressure anisotropy, electron inertia, etc. While these microscopic nonidealities play a crucial role in reconnection of magnetic reversals near electron- and ion-scales, they are totally negligible and irrelevant to the reconnection of structures at scales $\u2113\u226b\rho i,$ the ion gyroradius. This has been demonstrated in detail, starting with either the generalized Ohm's law (Eyink, 2015) or with the full Vlasov-Maxwell-Landau kinetic equations (Eyink and Drivas, 2018) In both cases, it was shown by a Favre coarse-graining analysis that the generalized Ohm's law at scales $>\u2113$ takes the form

and that all of the coarse-grained nonideal terms on the right hand side are damped out by powers $(\delta i/\u2113)$ and $(\delta i/\u2113)2$ for $\u2113\u226b\delta i,$ the ion skin depth. The nonideal electric field $R\u0303\u2113$ is thus irrelevant in the renormalization-group sense, and an “ideal Ohm's law” holds in the coarse-grained sense

at inertial-range length-scales $\u2113\u226b\rho i.$

Exactly as in our discussion of the MHD model, the validity of an ideal Ohm's law in this generalized sense does not imply that magnetic flux is frozen-in at inertial-range scales. Indeed, the reconnection of magnetic structures many orders of magnitude larger than ion scales is observed in the turbulent solar wind (Gosling, 2012; Eyink, 2015), and such reconnection will be observed, obviously, whether ion-scale eddies are resolved or not! The principle of renormalization-group invariance implies once again that the rate of reconnection of such large-scale structures must be independent of the resolution scale $\u2113.$ Within the coarse-grained description, the ideal Ohm's law at scales $\u2113\u226b\rho i$ can be rewritten as

where two apparently “nonideal” terms now appear, one due to velocity fluctuations

and the other due to density fluctuations

Here, *n* and $v$ are the number density and bulk velocity of the ions. Of course, these two terms are actually the effects of the ideal Ohm's law, and together, they account for all of the reconnection and breakdown of flux-freezing at length-scales $\u2113\u226b\rho i.$ Simple analytical estimates show indeed with $\epsilon \u2113:=\epsilon \u2113v+\epsilon \u2113n$ that $|\epsilon \u2113|\u226b|R\u0303\u2113|$ and also $|(\u2207\xd7\epsilon \u2113)\u22a5|\u226b|(\u2207\xd7R\u0303\u2113)\u22a5|$ so that reconnection and line-slippage at scales $\u2113\u226b\rho i$ are due to ideal turbulence effects and not due to microscopic plasma nonideality.

Just to be completely clear, this conclusion does not mean that large-scale turbulent reconnection in a weakly collisional plasma such as the solar wind will be described by ideal MHD equations. Although an ideal Ohm's law holds, there are not enough collisions to make the ion distribution function close to a local Maxwellian and to make the ion pressure tensor an isotropic function of local density and temperature. Thus, to provide a detailed quantitative description of large-scale reconnection in the solar wind, one must resort either to the full Vlasov-Maxwell kinetic equations or to gyrofluid models that take advantage of the small ion gyroradius to develop closures for the ion pressure tensor. Ideal MHD should, however, be a suitable model to investigate some of the basic effects of turbulence at inertial-range scales. Since the energetically dominant wave component of solar wind turbulence consists of incompressible shear-Alfvén modes, quantitative measures of field-line stochasticity and turbulent Richardson dispersion should carry over with little change to this collisionless plasma environment.

### F. Insight into physics of the turbulent cascade

The explanation of the violation of the textbook concept of flux freezing provided above is related to the properties of turbulence first discovered by Onsager (1949). In his prophetic paper, he wrote that ideal equations indeed hold in a turbulent inertial range, but the ordinary formulation of the laws of motion in terms of differential equations becomes inadequate and must be replaced by a more general description. In other words, the “ideal equations” (Euler, MHD, etc.) hold in the inertial range in a “generalized sense,” but it does not hold in the standard, naive sense of partial differential equations. In particular, it does not hold in a sense that implies conclusions such as conservation of total kinetic energy for Euler or magnetic flux-freezing for ideal MHD.

This profound insight by Lars Onsager is not properly appreciated by the community. It is frequently assumed that the modes at scale $>L$ should satisfy the standard ideal equations, in the usual sense of partial differential equations, for about one eddy turnover time at that scale, $tL\u2009L/uL$, until small scales $l\nu $ are created where viscosity/resistivity (or microscopic plasma effects in a collisionless plasma) becomes non-negligible. This is not right way of thinking. In a turbulent flow, all scales in the range $L>l>l\nu $ are already present at the initial instant. Low-pass filtering the fields to observe the scales $>L$ does not mean that that scales $<L$ have been physically removed. Turbulent motions of the cascade are physically present and cannot be disregarded. For instance, in the solar wind at every instant of time, there is a broad-band spectrum with excitation at all scales.

If one wants to know what is happening at a given scale $l>l\nu $. low-pass filtering the fields at this scale is a purely passive operation of observing only the modes at length scales $>l$ and ignoring those at scales $<l$. However, the modes at scales less than *l* are still present and cannot be ignored. Therefore, the dynamics at the scale *l* is not directly affected by the viscosity/resistivity/etc. in the equations for scales larger than *l*, but those equations are not the standard ideal equations. In the above presentation dealing with magnetic reconnection, we wrote out these ideal equations in physical space as partial differential equations but had to add the subgrid terms that arise from the modes at scales less than *l*. These subgrid terms are exactly what one gets by coarse-graining the ideal PDEs, i.e., one may say that the ideal equations hold in the coarse-grained sense. This is equivalent (as pointed out by Lars Onsager) to saying that Fourier coefficients of the fields satisfy ordinary differential equations that follow from the ideal PDEs. However, none of these equivalent notions of generalized solution have the usual implications of the ideal PDE's, such as conservation of total kinetic energy or conservation of magnetic-flux.

The necessity of using “a more general description” is obvious from a trivial hydrodynamic example. If the large scales $>L$ in hydrodynamic turbulence “obeyed Euler equations” in the usual, naive sense, then the kinetic energy at scales $>L$ would be conserved in time. However, it is not true. This is the easiest way to see that obeying Euler equations does not mean in this context the usual, naive thing. It is exactly this experimental observation that large eddies decay away, at a rate independent of viscosity, which motivated Onsager mathematical analysis.

In other words, the turbulent cascade connects the scales, and because of the motions associated with this cascade, it is not possible to consider motions at large scales as ideal. The real world is described at appropriate macroscopic scales by nonideal fluid or kinetic PDEs with viscosity, resistivity, collisions, etc. If you resolve all motions of a turbulent flow down to those smallest scales of validity, you will see the solutions of the nonideal PDEs. All dissipation and reconnection arise from the nonideal terns within this fine-grained description. However, for many problems, including the problem of astrophysical magnetic reconnection, it is natural to consider only large scales $l\u226bl\nu $. At these scales, the direct effect of nonideal, e.g., resistivity and Hall, terms is negligible. However, “subgrid terms” should be accounted for. With these terms, the equations are valid in the generalized sense, and we explicitly accounted for these terms in our quantitative discussion above.

### G. Rigorous description of flux freezing violation: Summary

Turbulent reconnection predicts reconnection taking place at all scales in the turbulent media, and this entails flux freezing violation. This section provides a rigorous description of the magnetic field behavior in turbulent fluids and demonstrates clearly that for considering the reconnection at scale *l*, one should account for the turbulent motions at the scale less than *l*. This vividly demonstrates and quantifies the violation of the conventional flux freezing. The approach proves that the nonideal plasma effects can be safely disregarded if magnetic reconnection is studied at scales much larger than the corresponding plasma scales.

This section provides strong support for conclusion in LV99, resolves a number of problems related to the earlier treatment of magnetic field in turbulent fluids, and establishes solid mathematical foundations of the turbulent reconnection theory. It also clarifies the dynamics of the magnetic field in turbulent media and suggests practical ways for evaluating the violation of the magnetic flux freezing.

## V. TESTING OF TURBULENT RECONNECTION AND FLUX FREEZING VIOLATION

### A. Testing LV99 predictions with MHD codes

Numerical testings of reconnection in the presence of turbulence have a long history. Pioneering two-dimensional studies in a periodic box were done by Matthaeus and Lamkin (1985) and (1986). Although the turbulence was only injected initially and not driven, it persevered whole simulation time. The results indicated the strong effect of turbulence presence on the reconnection process. For instance, the rates of production of reconnected magnetic islands were insensitive to resistivity. Unfortunately, the numerical setup precluded the calculation of the long-term reconnection speed, and no studies were done on how the properties of turbulence itself affect reconnection. More recently, Watson *et al.* (2007) studied the effects of small-scale turbulence on two-dimensional reconnection; however, no significant effects of turbulence on reconnection were observed. Some possible explanations for such effects were provided by the authors. Of the most important is that the perturbations were injected in the flow-driven merging of the super-Alfvénic regime, resulting in a quick ejection of the fluctuations from the sheet and two-dimensional configuration limiting the degree of freedom of magnetic field lines.

Dynamics of magnetic fields is three-dimensional by nature. Any two-dimensional configuration will significantly restrict the magnetic field freedom, making impossible, for example, the conventional magnetic field line wandering. Moreover, two-dimensional and three-dimensional magnetized turbulence are very different in nature (see the discussion in Eyink *et al.*, 2011). Nevertheless, due to the computational challenges and in the absence of any quantitative model to be tested, simulations aimed at studying reconnection have been performed in two-dimensions for a long time. For instance, the two-dimensional study in Matthaeus and Lamkin (1986) stresses the importance of turbulence for modifying the character of magnetic reconnection and specifies heating and transport as the effect of particular significance and the formation of Petschek-type X-point in two-dimensional turbulence. Kim and Diamond (2001) showed that the transport of magnetic flux is not enhanced to the reconnection zone. At the same time, field wandering described in LV99 is the essential feature of three-dimensional reconnection. The dynamics of field lines is very different in two dimensions (see Eyink *et al.*, 2011). Therefore, 2D processes cannot provide us with the guide to what is going on in the real astrophysical situations.

Loureiro *et al.* (2009) studied two-dimensional turbulent reconnection performing high-resolution simulations in a periodic box for Lundquist numbers up to $5\xd7104$ concluding the existence of a critical threshold of the injection rate above which the reconnection rate is enhanced. Kulpa-Dybeł *et al.* (2010) performed numerical simulations of the LV99 model in two dimensions showing that the relationships between the reconnection speed and turbulence properties are different than those predicted by the LV99 model, and the reconnection is not fast in the presence of turbulence since it still depends on the resistivity. In any case, it is obvious that the difference in physics of magnetic turbulence in 2D and 3D makes two dimensional studies a subject of pure academic interest. The actual testing of turbulent reconnection requires one to use the numerical setup of the realistic 3D setup.

At the time of publication of the theoretical LV99 model of turbulent reconnection, it was not possible to verify its predictions immediately using numerical simulations. The model is three-dimensional intrinsically, requiring significant computational resources, and requires good resolution to resolve all scales from large scales where turbulence is injected into dissipation scales at which magnetic resistivity becomes important. Even though the relationships between the reconnection rate and turbulence properties could be tested using relatively low resolutions, any attempt to test the most important prediction of the rate being independent of the magnetic dissipation would be exposed to immediate critics. Knowing that numerical simulations would not reach astrophysical Lundquist numbers for decades, it was challenging to undertake such a task.

The first extensive testing of LV99 by means of numerical simulations was performed in Kowal *et al.* (2009), a decade after the publication of LV99. The tests were done using the MHD compressible code in a three-dimensional domain with open boundary conditions in the sub-Alfvénic turbulence regime. Note that for simplicity the periodicity was applied along the guide field direction, i.e., *Z*-direction, but this does not affect our results, as there is no outflow/inflow along *Z*-direction.

Additional tests in terms of different turbulence driving and exploring the parameter space have been performed in Kowal *et al.* (2012). These pioneering provided new numerical approaches to exploring 3D magnetic reconnection. In particular, the ways of measuring the 3D reconnection rate were introduced there.

The authors performed numerical models changing three basic parameters, the power *P _{inj}* and scale

*k*of driving and the resistivity coefficient

_{inj}*η*. Figure 8 shows the measured dependencies of the reconnection speed

*V*on the turbulence parameters

_{rec}*P*and

_{inj}*k*. Even though relatively low resolution was used, the authors confirmed the theoretical LV99 model predictions of the reconnection rate dependence on the turbulence parameters. Most importantly, however, the numerical studies of Kowal

_{inj}*et al.*(2009) have proven that the reconnection is indeed fast in the presence of turbulence since it does not depend on the resistivity

*η*, as shown in Fig. 9. Moreover, the same dependence was observed for models with weak and strong guide fields

*B*. In Kowal

_{z}*et al.*(2009), the authors also tested the role of the anomalous resistivity and found that there is no essential dependence of the reconnection rate on anomalous effects, as well. Apart from the confirmation of these important relations, Kowal

*et al.*(2012) reported a weak dependence of the reconnection rate on viscosity, $Vrec\u223c\nu \u22121/4$, later studied in more detail in Jafari

*et al.*(2018) in the context of large magnetic Prandtl numbers.

As expected in the LV99 model, the current sheet is broad with individual currents distributed widely within a three dimensional volume and the turbulence within the reconnection region is similar to the turbulence within a statistically homogeneous volume. The structure of the reconnection region was analyzed in more detail by Vishniac *et al.* (2012) based on the numerical work by Kowal *et al.* (2009). The results supported the LV99 assumptions about the reconnection region being broad, the magnetic shear is more or less coincident with the outflow zone, and the turbulence within it is broadly similar to turbulence in a homogeneous system.

### B. Demonstration of flux freezing violation

As we discussed earlier, the LV99 idea of turbulent reconnection is closely related to the gross violation of flux freezing in turbulent environments. As we discussed in Sec. III B, Richardson dispersion is a vivid manifestation of the flux freezing violation in turbulent fluids. Richardson dispersion in space was derived in LV99, and it was used for predicting the rate of turbulent reconnection. Figure 6 confirms the corresponding LV99 predictions. Indeed, it is clearly seen that over the inertial range, the separation between following magnetic field lines is growing as $s3/2$, where *s* is the distance calculated along the magnetic field line.

A direct testing of the temporal Richardson dispersion of magnetic field-lines was performed recently in Eyink *et al.* (2013). For this experiment, stochastic fluid trajectories were tracked backward in time from a fixed point in order to determine which field lines at earlier times would arrive to that point and be resistively “glued together.” For the flux-freezing situation, the one to one mapping is expected for different times, but the Richardson dispersion prescribes quite a different behavior. To test what is right, many time frames of an MHD simulation were stored so that equations for the trajectories could be integrated backward. The results of this study are presented in Fig. 10. It shows the trajectories of the arriving magnetic field-lines, which are clearly widely dispersed backward in time, more resembling a spreading plume of smoke than a single “frozen-in” line. The quantitative results of the time-dependent separation correspond well to the expectations for the temporal Richardson dispersion. This numerical experiment also demonstrates that magnetic reconnection happens at every scale of the cascade as we discussed in Sec. IV F.

### C. Self-driven turbulent reconnection

Although studies on reconnection with imposed driven turbulence have been done for a few decades already, the question of self-driven reconnection was approached by aims of numerical simulations within last several years only. The question was first raised in LV99 and further discussed in Lazarian and Vishniac (2009). It was realized that due to the existence of magnetohydrodynamic instabilities, turbulence could be generated directly by the inhomogeneous current sheet, leaving the imposed external turbulence obsolete. In the self-driven reconnection, those inhomogeneities of the current sheet should grow, resulting in variations of the current sheet thickness and the size and strength of the outflow from the local reconnection regions, which interact with the developed turbulence. While the reconnection feeds these interactions, the turbulence should grow. The developed turbulence, in return, influences the reconnection process, resulting in its enhanced rate.

The difficulty in studying self-driven reconnection comes from several factors. First of all, in order for small amplitude fluctuations to not simply propagate out of the box or dissipate, one has to use very small Ohmic resistivities $\eta \u226a10\u22123$, corresponding to large Lundquist numbers, $S\u226b103$. However, the thickness of current density scales as $\u223c\eta 1/2$, e.g., in the Sweet-Parker model, indicating the necessity of using very large resolutions in numerical models. For instance, for Lundquist number $S=VAL/\eta \u223c106$, assuming $VA\u223c1$ and $L\u223c1$, the current sheet thickness would be around $10\u22123L$, where *L* is the longitudinal size if the current sheet. Resolving such a small thickness numerically requires at least resolutions of 10^{4} cells per unit length *L*. Moreover, the goal is to reach a steady-state of self-driven reconnection, which cannot be achieved in a periodic box, and one has to use numerical domains with open boundary conditions, proper implementation of which is far from trivial. The steady-state also cannot be achieved using very short runs; therefore, one has to simulate a model for at least a few Alfvén times. Finally, one would like to study realistic three-dimensional configurations, which for the estimated resolutions would require enormous computational resources. There are only a few most powerful supercomputers which could deal with such a huge task nowadays.

The numerical simulations with turbulence generated through magnetic reconnection were presented in preprint by Beresnyak (2013) with the final extended version of the paper published in Beresnyak (2017). The studies were done in a three-dimensional periodic box using spectral incompressible MHD simulations up to resolution 1536^{3} with the highest Lunquist number of $S=6.4\xd7104$. The incompressible simulations used correspond to the limiting case of high beta plasma, i.e., plasma with sound velocity $\u226b$ the Alfvén one. While in 2D, the location of the current sheet was remained in the same position, in 3D, the location of the original current sheet position was largely forgotten as turbulence was developing.

The study reported magnetic reconnection that is very different from the magnetic reconnection studied earlier in 2D settings. First of all, the reconnection was observed to proceed at a steady rate, which is in contrast to the 2D simulations [see Loureiro *et al.* (2012)], exhibiting a significant time dependence. The flux ropes were observed in 3D, but, unlike magnetic islands in 2D, the flux ropes were turbulent inside. The number of these structures did not depend on the Lunquist number, contrary to the tearing reconnection described, e.g., in Uzdensky *et al.* (2010).

Clear turbulent structures were developed from the initial velocity noise with the kinetic energy of the order of $10\u22126$ of the magnetic one within several Alfvén times (see Fig. 2 in Beresnyak, 2017). The reported reconnection rate was measured to be $\u223c1.5%$ of the Alfvén speed (see Fig. 11). The reconnection rate was derived as $Vrec=d\Delta /dt$, where Δ is the measured of the turbulent layer thickness, and this result was justified using LV99 theory [a somewhat different derivation of the growth of Δ is obtained in Lazarian *et al.* (2015)]. The analysis of the data of Beresnyak (2017) concluded that the expression for the LV99 reconnection rate $Vrec=CLV\delta vL2/VA,x$ should have the constant *C _{LV}* very close to unity, with

*C*= 0.94 being claimed.

_{LV}Beresnyak (2017) also studied the properties of generated turbulence, reporting turbulence consistent with Goldreich-Sridhar (1995), i.e., kinetic energy power slope $\u223c\u22121.5\u2009to\u22121.7$ and scale-dependent anisotropy $\lambda \u2225/\lambda \u22a5\u223c\lambda \u22a5\u22121/3$, although the inertial range of compatible anisotropy was relatively short.

Later works include numerical simulations of stochastic reconnection done by Oishi *et al.* (2015) and Huang and Bhattacharjee (2016), both using compressible MHD approximation, in slightly different reconnection configurations, however. Oishi *et al.* (2015) used a numerical setup starting from the Sweet–Parker setup with a periodic box along the current sheet. This setup is somewhat unfortunate since the Sweet–Parker configuration requires open boundary conditions to allow inflow and outflow of the magnetic flux in order to maintain the steady-state. Due to periodic boundary conditions, the growth of the fluctuations at large scales may be attributed to the presence of the magnetic field setup which is not initially in the equilibrium. Nevertheless, they performed several long time simulations with Lundquist numbers up to $S\u223c3.2\xd7106$ with the conclusion that for sufficiently large Lundquist numbers $S\u2273105$, the thickness of current sheet $\delta SP$ becomes weakly dependent on *S* (see their Fig. 1). However, it should be noted that simulations with $S>105$ ($\eta <10\u22125$) where done with the maximum effective resolution 1024^{3}, much lower than the resolutions needed to resolve the thickness of current sheet estimated above; therefore, the conclusion might be attributed to the numerical effects. They also estimated the reconnection rate using the decay rate of the magnetic energy and concluding that 3D MHD stochastic reconnection can be fast without recourse to kinetic effects or external turbulence. They did not, however, provide any studied properties of reconnection-generated turbulence.

Huang and Bhattacharjee (2016) analyzed the statistics of velocity fluctuations generated by reconnection also in the presence of the initial velocity noise and used, as a starting point, a global Sweet–Parker configuration, but with nonperiodic boundaries along the current sheet and a guide field comparable in strength to the reconnecting component, which additionally prevented the interactions of the waves. Weak velocity perturbations were injected into such a domain in order to understand their effect on the reconnection rate. The reported turbulent power spectra of kinetic energy with slopes $\u22122.5\u2009to\u22122.3$ and observed nearly scale independent anisotropy of velocity fluctuations, both in contrast with the predictions of the Goldreich-Sridhar theory. As pointed out later by Kowal *et al.* (2017), the power spectra and structure functions calculated in Huang and Bhattacharjee (2016) from two-dimensional *XZ*-planes and later averaged cause strong bias, resulting in steeper slopes and much reduced anisotropy, when compared to the fully three-dimensional spectral and structure function analysis.

While the study in Huang and Bhattacharjee (2016) supports the LV99 expectation to the transfer of reconnection to the turbulent regime, some statements of this study contradict the findings in all other numerical studies by performed by different groups around the world, which also studied 3D reconnection. In particular, the statement of the similarity of 2D and 3D magnetic reconnection in Huang and Bhattacharjee (2016) is at odds with the findings by other authors who report that in 3D, the reconnection evolves differently from that in 2D (see Oishi *et al.*, 2015; Wang *et* *al*., 2015; Striani *et al.*, 2016; Beresnyak, 2017; Kowal *et al.*, 2017; Takamoto, 2018; Wang and Yokoyama, 2019; Kowal *et al.*, 2019). This statement in Huang and Bhattacharjee (2016) and their finding that the spectral slope and anisotropy of the obtained turbulence are different from expectations of MHD turbulence theory (see Sec. II B), which contradict the LV99 expectations. Below we provide the results of numerical testing of these claims in Huang and Bhattacharjee (2016).

Below we discuss the most recent studies of the reconnection-driven turbulence in Kowal *et al.* (2017). As the initial setup, they used a uniform box with a single discontinuity of the X-component of the magnetic field placed in the middle of the box along the XZ plane. Although they used periodicity along the X and Z directions, the box was open in the vertical direction, perpendicular to the current sheet, and significantly elongated in Y, with physical dimensions of $1.0\xd74.0\xd71.0$ in order to prevent any influence of the open boundary conditions. They let the turbulence grow from the initial noise of amplitude 0.01 applied to each velocity component for long runs up to $t=20tA$. Figure 12 shows three main cuts across the domain of the vorticity amplitude at the final time of one of the investigated models.

The choice of different vertical boundaries had some important consequences. Periodicity in all three directions with two separated current sheets imposed, as used in Beresnyak (2017) for instance, allows for the possibility of horizontal large scale interactions through the deformations of both current sheets. These interactions, especially at later times, can generate large scale motions in the perpendicular direction to the current sheet, providing additional reconnection energy input (see, e.g., Kowal *et al.*, 2011, for a similar setup with many current sheets, in which the current sheet deformations are well manifested). In this case, such large-scale interactions are not allowed; therefore, the only energy input comes from the small-scale reconnection events. They, contrary to other studies, investigated the dependence of the properties of developed turbulence on the *β*-plasma parameter. They reported Kolmogorov-like power spectra of generated velocity fluctuations (see Fig. 13), which also presented scale-dependent Goldreich–Sridhar anisotropy scaling at later times in large *β* models (shown in Fig. 14). They claimed that the turbulent statistics are similar to strong MHD turbulence, but they are significantly affected by the dynamics of reconnection induced flows, more clearly manifested in the low-*β* regime, where supersonic reconnection outflows were present. In the high-*β* regime, the Goldreich–Sridhar anisotropy scaling manifested at much earlier times due to the fluctuations generated by reconnection outflows propagating and interacting faster. They also estimated the reconnection rate to be twice as high as those estimated in Beresnyak (2017), which they justified by a more diffusive code.

Kowal *et al.* (2019) studied what exactly drives turbulence in their stochastic reconnection models. They considered two instabilities, tearing mode and Kelvin-Helmholtz instability. For each analyzed model, they determined the magnetic (in the case of tearing mode) and velocity (in the case of Kelvin-Helmholtz instability) shear location and performed a detailed analysis of the field profiles taking into account local geometry, i.e., the local coordinate system oriented with respect to the shear plane. Using such analysis, they could determine a number of quantities important for the development of instabilities. For instance, in order to study the tearing mode, they determined the local thickness of the current sheet, the length of the local sheet, the strength of the transverse, and the guide components of the magnetic field. For Kelvin-Helmholtz instability, they determined the shear amplitudes and the strength of the local magnetic field component which stabilized the instability or sonic and Alfvénic Mach numbers. These parameters allowed us also to estimate the growth rates for both instabilities.

The results indicated that tearing mode growth is significantly suppressed at all studied wave numbers. See, for example, Fig. 15, where the statistical distribution of the growth rates for the tearing mode (left panel) and Kelvin-Helmholtz instability (right panel) for three different perturbation wave numbers are shown. The strong suppression of the tearing mode can be explained by the fact that its growth rate actually decreases with the wave number *k*, while for Kelvin-Helmholtz instability, it increases. It could be also attributed to the presence of the transverse component of the magnetic field, easily generated by turbulence. Figure 16 shows the distribution of the transverse component relative strength $\xi =Bn/B$ vs $vA\delta $. It is well seen that a significant number of cells within current sheets are stable due to the strong transverse component.

It is evident that the study in Kowal *et al.* (2019) supports the expectation of the turbulent reconnection theory and disagrees with some of the conclusions in Huang and Bhattacharjee (2016). In particular, the spectrum and anisotropy that are reported are consistent with the expectations of the MHD turbulence theory (see Sec. II B). In addition, Kowal *et al.* (2019) testified that, in agreement with other earlier studies, the reconnection in 3D is different from that in 2D. In particular, the tearing modes that dominate the reconnection dynamics in 2D are replaced by perturbation induced by Kelvin-Helmholtz instability. In other words, the steady state self-driven reconnection is not driven by tearing/plasmoids. The issue why the results in Huang and Bhattacharjee (2016) are different requires further studies. One possibility is that the development of a true turbulent state requires more time than the current sheet was evolving in the simulations in Huang and Bhattacharjee (2016).

Studies of 3D self-driven turbulent magnetic reconnection with different initial conditions are important. In this respect, it is worth noting other studies of turbulence generating by stochastic reconnection are different than the Harris current sheet configurations. For example, turbulence developed through reconnection in a three-dimensional cylindrical plasma column was investigated in Striani *et al.* (2016). An example visualization is shown in Fig. 17. The authors report a fragmentation of the initial current sheet in small filaments, which leads to an enhanced reconnection rate, independent of the Lundquist number already at $S\u223c103$. Moreover, they observe at the first phase, before the initial perturbation develop turbulence, the classical Sweet-Parker dependence of the reconnection rate on the resistivity. Once the turbulence grows, the reconnection rate increases quickly and becomes independent of the Lundquist number (see Fig. 18). Apparently, this study also shows that 2D and 3D self-driven reconnections are different.

### D. Testing turbulent reconnection with the Hall-MHD code

One may ask a question whether the magnetic reconnection that we described above for the MHD case in both the externally driven turbulence and the case of self-driven turbulence is going to be different from the reconnection in actual plasmas. The LV99 theory gives a negative answer to this question and states that the microphysical plasma affects are not important at the scales much larger than the ion gyroradius or/and the ion skin depth $di=c/\omega pi$, where *ω _{pi}* is the ion plasma frequency. This conclusion corresponds to both the common knowledge that the MHD turbulence properties that should not depend on the physics at the dissipation scale and the accepted practice of representing magnetized flows on the scales much larger than

*d*using the MHD equations. However, the direct testing of the effects of plasma microphysics for magnetic reconnection on much larger scales is extremely challenging numerically. We address the approaches to this problem in the Sec. VI. In this section, we discuss simplified ways to perform this testing.

_{i}To test the LV99 prediction that plasma effects are not irrelevant for large scale turbulent reconnection, Kowal *et al.* (2009) performed the simulations of turbulent reconnection where the plasma effects were simulated through anomalous resistivity. No change in the reconnection was observed as the anomalous resistivity was varied.

For the case of reconnection with self-generated turbulence, Beresnyak (2018) performed Hall-MHD simulations with the resolution up to $2304\xd746082$. The simulations had different ratios of the ion-electron depth *d _{i}* to the grid size

*d*, namely, $di/d$ values were chosen to be 9.5 and 39. The latter ratio approaches the ratio of the ion to electron skin depths; the former is used to obtain follow the scaling for larger current sheet widths.

Figure 19 shows the dependence of the Hall-MHD reconnection rate *v _{r}* as a function of the width

*W*of the reconnection layer. Similar to the earlier MHD reconnection simulations in Beresnyak (2017), the reconnection rate in the Hall-MHD simulations was measured as the increase in the reconnection layer thickness [a quantitative study of this measure of reconnection for the simulations with periodic boundary conditions and its relation to the LV99 theory is provided by Lazarian

*et al.*(2015) in Sec. IV D there]

It can be shown that the differences between the setups with different $di/d$ values disappear when the width of the reconnection layer gets larger than $3.6di$.

The current layer in the simulations was getting fully turbulent with magnetic fluctuations exhibiting whistler turbulence (Sec. VIII C) below the *d _{i}* scale. The original rate of reconnection in the setup was larger than in the case of MHD self-excited reconnection. However, this rate was decreasing with the increase in the width of the reconnection layer. In other words, these findings support the notion that turbulent reconnection on the scales where one expects the MHD approximation to be valid (see a discussion of this in Eyink

*et al.*, 2011) does not depend on the underlying microphysical plasma processes. This corresponds to both the physical intuition and the accepted practice of using MHD for modeling the processes at scales much larger than the scales of the proton gyroradius.

We note that while the main achievement of the study by Beresnyak (2018) was to demonstrate that the effects of Hall-MHD get negligible for turbulent reconnection at scales $\u226bdi$, the study is also of importance for the reconnection at sufficiently small scales. Such a type of reconnection can take place, e.g., in the Earth's magnetosphere, and the study demonstrates how the turbulent plasma-scale reconnection gradually transfers into the reconnection that is determined entirely by the MHD turbulence. We, however, expect that the reconnection will depend on the plasma *β*, and the corresponding studies would be useful.

### E. MHD and Hall-MHD testing of turbulent reconnection: Summary

3D magnetic turbulent reconnection is notoriously difficult to study. The initial studies employed the externally driven turbulence and confirmed the analytical dependences for the reconnection rate predicted in LV99. Similarly, the violation of the conventional flux freezing predicted by the turbulent reconnection theory was also demonstrated in simulations with externally driven turbulence.

How reconnection proceeds in the situation when no turbulence is initially present is the focus of the ongoing research on the spontaneous 3D reconnection. This sort of reconnection is traditionally studied in 2D results, and such studies in the MHD limit established the tearing/plasmoid reconnection model. However, the 3D simulations clearly show that for sufficiently large Lundquist numbers, the reconnection layer in 3D gets turbulent, and this makes the reconnection in 2D and 3D is very different as it is established in many independent studies. 3D studies are computationally costly, but, unlike 2D ones, they are astrophysically relevant.

Quantitative measurements of the properties of the magnetic field in self-driven 3D reconnection are necessary. The first measurements of the turbulence spectrum in the 3D reconnection layer, turbulence anisotropy, and estimates of the growth rate of Kelvin-Helmholtz instability are very encouraging. These results support the predictions of turbulent reconnection theory. More studies of self-driven turbulent reconnection are necessary for exploring the phenomenon of flares induced by magnetic reconnection.

The work that shows that 3D self-driven reconnection naturally transfers to the turbulent regime testifies that the turbulent reconnection is a generic type of astrophysical large-scale reconnection. At the same time, tearing and various plasma effects can be important only in special circumstances, e.g., in the transient regime before the system gets turbulent. For studies of 3D reconnection, it is important to have high resolution. Otherwise, the numerical viscosity can suppress turbulence, completely distorting the real picture of astrophysical reconnection.

A recent study of the self-driven turbulent reconnection with the Hall term confirms the prediction of the turbulent reconnection theory that small-scale plasma effects do not change the rate of large scale reconnection in turbulent fluids. This is a confirmation of a practically important conclusion in LV99 that the large scale processes in turbulent fluids can be simulated using MHD codes.

## VI. TURBULENT RECONNECTION WITH KINETIC SIMULATIONS

Kinetic simulations offer another perspective at understanding the role of turbulence on reconnection. As energy cascades to smaller scales, the interactions between particles and fields could become increasingly important, including many collisionless processes. Fully kinetic simulations can track the motion of ions and electrons rigorously, and this self-consistent “closure” is different from the resistivity and viscosity treatment often invoked in fluid studies. In fact, a direct comparison (Makwana *et al.*, 2015) between the compressible MHD and kinetic simulations of the turbulent cascade process has demonstrated the validity of using particle-in-cell (PIC) codes to study the nonlinear processes in the continuum limit at large scales.

As described earlier in this review, it is instructive to quantify the field line diffusion behavior in kinetic simulations of reconnection. In addition to the often pre-existing turbulent background and the collisionless tearing driving scales, the 3D large-scale kinetic reconnection process also contains several characteristic scales associated with ions and electrons that naturally lead to kinetic instabilities on small scales (e.g., Bowers and Li, 2007; Daughton *et al.*, 2011; Guo *et al.*, 2015) and hydrodynamic instabilities at large scales (see Kowal *et al.*, 2019), which can produce fluctuations on various scales, and one might expect that the magnetic field diffusion at the microscopic plasma scales to be more complicated than its fluid counterpart. Below we describe some initial results based on one 3D PIC simulation of magnetic reconnection. It is very difficult to get fully turbulent fluid behavior for the outflow from the reconnection region. Therefore, in what follows the study is focused on exploring the separation of the field lines, which is related to the Richardson dispersion (see Sec. III B).

### A. Numerical method and initial setup

3D kinetic simulations were carried out using the VPIC code (Bowers *et al.*, 2008). VPIC is a first-principles electromagnetic charge-conserving code solving Maxwell's equations and the Vlasov equation in a fully relativistic manner. It has been carefully designed and optimized on recent peta-scale computers, which enables multidimensional simulations with more than 1 × 10^{12} particles and 5 × 10^{9} numerical cells. Using the newly implemented Trinity machine at LANL, we have successfully run simulations that contain 5.2 × 10^{12} particles with about 17 × 10^{9} cells.

The simulation presented here (Guo *et al.*, 2019) starts from a force-free current sheet with $B=B0tanh(z/\lambda )x\u0302+B0sech(z/\lambda )y\u0302$, where *B*_{0} is the strength of the reconnecting magnetic field and *λ* is the half-thickness of the current sheet. We assume a positron-electron pair plasma with $mi/me=1$. Initially, $\lambda =12de$, which implies that the growth rate for collisionless tearing is quite small. Here, $de=di=c/\omega pe=c/4\pi nie2/me$ is the electron/ion inertial length. The initial particle distributions are Maxwellian with uniform density *n*_{0} and temperature $Ti=Te=mec2$. The magnetization parameter is $\sigma =(\Omega ce/\omega pe)2=100$. The domain size is $Lx\xd7Lz\xd7Ly=1000de\xd7500de\xd7500de$. The grid sizes are $Nx\xd7Nz\xd7Ny=4096\xd72048\xd72048$. The 150 particles per cell per species were used. For electric and magnetic fields, periodic boundaries along the *x* and *y* directions and perfectly conducting boundaries along the *z*-direction were employed. For particles, we employ periodic boundaries along the *x* and *y*-directions and reflecting boundaries along the *z*-direction. The simulation lasts $\omega pet\u22481050$, which is about the time for an Alfvén wave to travel *L _{x}*. The time scale to traverse the current sheet thickness, however, $\omega pe\tau =2\lambda /c=24$, is much shorter.

Different from the earlier simulations, we inject an array of perturbations with different wavelengths with the total summed amplitude $(dB/B0)2\u223c0.1$ to initiate a background turbulence at the beginning of the simulation. [Injection of the initial perturbations follows the procedure described in Makwana *et al.* (2015).] These initially injected perturbations have wavelengths longer than the initial thickness of the current sheet.

### B. Analysis

Figure 20 depicts the overall global reconnection rate derived from this simulation, and Fig. 21 shows several snapshots of the total current density $|J|$, indicating the evolution of many flux ropes (in 3D) and the development of subsequent thin current sheets. The reconnection rate is evaluated using the method described in (Daughton *et al.*, 2014a), where the mixing of electrons originating from separate sides of the initial current layer is used as a proxy to rapidly identify the magnetic topology and track the evolution of magnetic flux. It was found that, while the distribution and dynamics of reconnection X-points are strongly modified by the injected fluctuations and self-generated fluctuations due to secondary instabilities, the peak rate is on the order of 0.1 times of the Alfvén speed, consistent with previous 2D and 3D kinetic studies. The reconnection changes by about a factor of 2 from its peak to its final quasisteady rate.

To quantify the magnetic field diffusion during the reconnection, the following procedure was adopted: using an output at a particular time snapshot, the two field lines were traced by starting from a pair position which is closely spaced (the typical initial separation is $s0=0.01di$). The separation between them $s(di)$ as a function of the field line path length $l(di)$ was calculated. Pairs at different initial distances (above and below) from the initial current sheet layer center are chosen randomly. Typically 1000 initial pairs to enhance the statistics were used. Figure 22 shows some sample magnetic field lines. There are 100 pairs within $0.01de$ of each other, whose starting locations are randomly chosen within a $2\xd72\xd72de$ box that is $2\lambda $ above the middle plane of the initial current sheet at $x=500,\u2009y=250,\u2009z=256de$. The time snapshot is at *t *=* *1280 or $\omega pet=126.3$, which is in the early stage of reconnection (before the peak, cf. Fig. 20). These field lines are integrated for up to $500de$ in total length, and we stop following them when they reach any simulation boundary. The initial turbulence is apparently causing a fraction of magnetic field lines to go cross the current sheet. This could have interesting implications for initiating reconnection.

Figure 23 gives the details of magnetic field separation behavior. Overall, we find that magnetic field lines indeed separate from each other at a rate faster than the regular diffusion process. The slopes from the power-law portion cluster around $1.24\u22121.35$. In addition, we have chosen two initial heights ($2\lambda $ and $4\lambda $) for tracing field lines, and their behavior is quite similar. Furthermore, such analyses were done using three different time snapshots capturing the at-peak and postpeak reconnection stages. Within the uncertainties of the fits, both the amplitudes and the exponents stay constant in time.

The exponent ranging between $1.24\u22121.35$ suggests that magnetic field lines exhibit superdiffusive behavior, consistent with the turbulent reconnection theory as discussed earlier. The relatively narrow range of this exponent also suggests that the initially injected turbulence might have strongly regulated the diffusion behavior, though the turbulence produced by 3D kinetic reconnection could have impacted the diffusive dynamics. More detailed analyses are needed to differentiate whether the injected and self-generated turbulence could lead to different field line diffusion behaviors and how they might interact with each other.

### C. Kinetic simulations of turbulent reconnection: Summary

By analyzing the magnetic field line behavior based on a large 3D kinetic simulation of reconnection that is subject to an initial turbulence injection, it was found that first, the overall global reconnection rate reaches 0.1 Aflvén speed, similar to previous kinetic simulations without the initial turbulence injection. This is perhaps not surprising, given that the kinetic physics is probably still dominating the dynamics despite the relatively large simulation box. Second, the presence of the initial turbulence, however, drastically shortens the stage through which the peak reconnection rate is achieved. This suggests that the turbulence can accelerate the “triggering” process. Third, magnetic field lines exhibit superdiffusive behavior throughout the reconnection process, generally consistent with the superdiffusion theory (LV99; Eyink *et al.,* 2011; Lazarian and Yan 2014) and previous 3D MHD reconnection studies (see Kowal *et al.*, 2017). More studies are needed to better understand the detailed dynamics and to delineate any possible different roles played by the initial background turbulence vs the self-generated turbulence via reconnection.

## VII. OBSERVATIONS OF TURBULENT RECONNECTION

Magnetic reconnection is notoriously difficult to study directly. In most cases, one can see the consequences of the reconnection, e.g., heat release or energetic particles accelerated, rather than the actual picture of reconnection. For turbulent large scale reconnection, the situation is really disadvantageous, as most of the *in situ* magnetic reconnection measurements have been done for the Earth's magnetosphere where the thickness of the current sheets is comparable to the ion Larmor radius. In such settings, one does not observe MHD-type turbulence, and therefore, magnetospheric measurements provide information about plasma reconnection dominating small scales rather than large scale reconnection relevant to most of the astrophysical settings. Nevertheless, there have been both observations and measurements in the relevant parameter regime and those support the turbulent reconnection model.

### A. Solar flares

Turbulent reconnection provides the most natural explanation for the variations of solar activity. Indeed, the reconnection speed that can be inferred from the observations of the magnetic activity in the solar atmosphere varies dramatically. This is difficult to account on the basis of any plasma reconnection process that prescribes a given reconnection rate. On the contrary, as we discussed earlier, the LV99 predictions suggest that the rate of magnetic reconnection varies with the level of turbulence, and this level of turbulence can be affected by the reconnection itself. The latter provides a natural explanation of solar flares (LV99, Lazarian and Vishniac, 2009).

One can estimate the reconnection rates for solar reconnection, assuming that turbulent velocity dispersion $Uobs,turb$ measured during the reconnection events is due to the strong MHD turbulence regime. This is reasonable as we expect strong bending of magnetic field lines in the reconnection region. Therefore

Similarly, the thickness of the reconnection layer should be defined as

For solar atmospheres, the expected ratio of Δ to the proton gyroradius *ρ _{i}* is 10

^{6}or even larger. Therefore, we expect to deal with the LV99 type of reconnection.

The expressions given by Eqs. (70) and (71) should be compared with observations. Using the corresponding data in Ciaravella and Raymond (2008), one can obtain the widths of the reconnection regions in the range from 0.08*L _{x}* up to 0.16

*L*and the observed Doppler velocities in the units of

_{x}*V*of the order of 0.1 (see more in Lazarian

_{A}*et al.*, 2015). It is easy to see that these values agree well with the predictions given by Eq. (71). In fact, the corresponding comparison was first done in Ciaravella and Raymond (2008) in order to observationally test the LV99 theory. However, the authors assumed that turbulence was isotropically driven, which was unrealistic for the solar flare case. While Ciaravella and Raymond (2008) provide, within the observational uncertainties, a satisfactory quantitative agreement between their observations and the LV99 theory prediction, it is demonstrated in Lazarian

*et al.*(2015) that this agreement can be further improved if the actual anisotropic driving is accounted for.

A more recent study exploring the growth of the thickness of the reconnection layer of the microflaring loop is presented in Chitta and Lazarian (2019). This study shows that the turbulent velocity measured using Doppler-broadened lines and the width of the reconnection region are in a good agreement with the predictions of the LV99 theory.

A very distinct prediction in the LV99 theory is the spread of magnetic reconnection from an active reconnection region to the adjacent reconnection regions. Indeed, the turbulence induced by the reconnection in one flaring region is expected to destabilize the adjacent regions that are undergoing relatively slow reconnection and induce flares of reconnection in them. It is evident that one should not expect any similar effects in the models of plasma reconnection. In Sych *et al.* (2009) and (2015), the authors observed the effect similar to one we just described. They explained quasiperiodic pulsations in observed flaring energy releases in an active region above the sunspot and proposed that the wave packets arising from the sunspots can trigger such pulsations via the LV99 mechanism.

Plasmoids associated with tearing reconnection have also been observed (Shibata and Tanuma, 2001) and on the basis of our experience with self-driven turbulent reconnection that we described in Sec. V C. We associate the plasmoid stage with the onset of turbulent reconnection.

### B. Reconnection in solar wind

Reconnection throughout most of the heliosphere is expected to be similar to that in the sun, and therefore, we do expect that the LV99 theory can be applicable to it. However, for the measurement interpretation, one should treat the small scale and large scale reconnection events separately. For example, there are now extensive observations of strong narrow current sheets in the solar wind (Gosling, 2012). The most intense current sheets observed in the solar wind have widths of a few tens of the proton inertial length *d _{i}* or proton gyroradius

*ρ*(whichever is larger). The reconnection associated with these small-scale current sheets of exhibits exhausts with widths of a few hundreds of ion inertial lengths (Gosling

_{i}*et al.*, 2007; Gosling and Szabo, 2008). Naturally, such small-scale reconnection events require collisionless physics for their description (Vasquez

*et al.*, 2007). Within the LV99 picture, these are small scale events representing the microphysics but not much affecting the large scale reconnection.

At the same time, there are very large-scale reconnection events in the solar wind, often associated with interplanetary coronal mass ejections and magnetic clouds or occasionally magnetic disconnection events at the heliospheric current sheet (Phan *et al.*, 2009; Gosling, 2012). These events have reconnection outflows with widths up to nearly 10^{5} of the ion *d _{i}* and exhibit a prolonged, quasistationary regime with reconnection lasting for several hours. Such large-scale reconnection is expected within the LV99 theory when very large flux-structures with oppositely directed components of the magnetic field interact with each other in the turbulent environment of the solar wind. The “current sheet” producing such large-scale reconnection in the LV99 theory contains itself many ion-scale, intense current sheets embedded in a diffuse turbulent background of weaker current.

The comparison between the reconnection events in High-Speed Solar Wind and those in MHD turbulence simulations was performed in Lalescu *et al.* (2015). In the latter simulations, the magnetic field undergo Richardson diffusion as proven in Eyink *et al.* (2013), and therefore, their turbulent reconnection corresponds to LV99 theory. Figure 24 shows that an inertial-range magnetic reconnection event in an MHD simulation is similar to those reconnection events observed in solar wind. The authors observed stochastic wandering of magnetic field lines in space, breakdown of conventional flux-freezing principle as a result Richardson dispersion of magnetic field lines, and also a broadened reconnection region containing several current sheets. The coarse-grain magnetic geometry in this study shown in Fig. 25 appears as a large scale reconnection event in the solar wind, while the detailed picture of magnetic field lines exhibits a much more stochastic turbulent nature of the event. Lalescu *et al.* (2015) concluded that the breakdown of conventional magnetic flux-freezing as an evidence of reconnection and also verified topology changes in both fine-grained and coarse-grained magnetic fields.

We believe that further studies of solar wind reconnection is a good way to test the LV99 model and the related theories of flux freezing violation and slippage of magnetic field lines that we discussed in the review. In view of this, we note that LV99 better is applicable to reconnection events more distant from the sun because the scales of the events increase. For instance, reconnecting flux structures in the inner heliosheath could have sizes up to ∼100 AU, much larger than the ion cyclotron radius $\u223c103$ km (see Lazarian and Opher, 2009).

### C. Heliospheric current sheet and Parker spiral

The violation of flux freezing required by the model of turbulent reconnection can also be directly observable. For instance, Eyink (2015) analyzed the data relevant to the region associated with the broadened heliospheric current sheet (HCS), noticed its turbulent nature, and argued that LV99 magnetic reconnection can explain the properties of the HCS. Naturally, more dedicated studies of the phenomenon are necessary.

Parker's spiral model (Parker, 1958) of the interplanetary magnetic field remains as one of the applications of the magnetic flux-freezing. The deviations from this model have been observed, however. Using yearly averaged magnetic field strengths calculated from the observations of Voyager 1 and 2, Burlaga *et al.* (2002) found a mean deviation of −2% between the observations and the Parker model. Yearly deviations were around $\xb120%$ with the maximum error (of about −40%). In addition to magnetic field strengths studied by Burlaga *et al.* (2002), an earlier work by Burlaga *et al.* (1982) considered also the magnetic field geometry which revealed notable deviations from the spiral model. This work found a different radial dependence for the radial and azimuthal magnetic field components from that predicted by the Parker model. These observations were based on the data collected by Voyager 1 and 2, while a comparison with earlier Pioneer 10 and 11 investigations during the less active solar activity period (1972 through 1976) confirmed the spiral model modulo variability on timescales shorter than a solar rotation period.

These early observations were recently confirmed by Khabarova and Obridko (2012). The extensive averaging used in the latter work eliminated the sizable fluctuations still observed in the daily averages of Burlaga *et al.* (1982); see their Fig. 1. Khabarova and Obridko (2012) interpreted their observations due to a quasicontinuous magnetic reconnection, occurring both at the heliospheric current sheet and at local current sheets inside the IMF sectors. They present extensive evidence that most nulls of radial and azimuthal magnetic field components (i.e., potential reconnection sites) are not associated with the heliospheric current sheet. This interpretation is consistent with the results of Eyink (2015) as magnetic slippage through the fluid as a result of pervasive turbulence in the near-ecliptic solar wind leads to a less tightly wound spiral and a stronger radial field strength than in the Parker model, as observed by Khabarova and Obridko (2012). This is intimately related to the reconnection diffusion proposed by Lazarian (2005) and first studied numerically by Santos-Lima *et al.* (2010) (see more in Sec. IX B).

In addition, Eyink (2015) analyzed the data relevant to the region associated with the broadened heliospheric current sheet (HCS), noticed its turbulent nature, and provided arguments on the applicability of the LV99 magnetic reconnection model to HCS. This seems to be a very promising direction of research to study turbulent reconnection in action using *in situ* spacecraft measurements.

### D. Observational testing turbulent reconnection: Summary

Quantitative studies of magnetic reconnection from observations and *in situ* spacecraft measurements are notoriously difficult. However, in the last 20 years, a number of predictions of the turbulent reconnection theory have been successfully tested. For instance, the prediction of the relationship between the level of turbulence and the width of the reconnection layer was confirmed with solar observations, and the deviations of the Parker spiral and heliospheric current sheet were explained with the turbulent reconnection theory. Very importantly, the predicted effect of inducing magnetic reconnection by Alfvénic perturbations coming from the neighboring reconnection event was observed.

The quantitative studies of the reconnection in solar wind showed that the reconnection events identified there have direct analogs in the reconnection events observed in numerical simulations of MHD turbulence. We believe that solar wind studies can provide the best testing ground for quantitatively comparing the prediction of the theory of turbulent reconnection with the *in situ* measurements of reconnection. At the same time, the magnetospheric reconnection events cannot serve the same purpose as they happen in different regimes that are dominated by plasma effects. The existing theory of turbulent reconnection cannot be directly applied for describing such reconnection events.

## VIII. TURBULENT RECONNECTION IN SPECIAL CONDITIONS

### A. Stochastic reconnection for large *Pr*_{m}

_{m}

Partially ionized, turbulent regions of the interstellar medium (ISM) are often associated with a high magnetic Prandtl number *Pr _{m}* (see, e.g., Xu and Lazarian, 2017b). This regime of large

*Pr*is also prevalent in the laboratory fusion plasmas, different phases of star formation, and in galaxies, protogalaxies, and clusters which often are filled with hot and rarefied plasmas (Jafari

_{m}*et al.*, 2018). Lazarian

*et al.*(2004) proposed a model for magnetic reconnection in partially ionized gases extending the arguments of stochastic reconnection (LV99) which was concerned with $Prm\u22431$. Another attempt, to study the magnetic reconnection in a large

*Pr*regime, was made by Lazarian

_{m}*et al.*(2015) who considered the role of viscosity in the diffusion of field lines caused by neutrals in a partially ionized gas. To reach a better agreement with numerical simulations (Kowal

*et al.*, 2009, 2012), in a more recent work, Jafari

*et al.*(2018) have considered the effects of viscosity on fluid motions along magnetic field lines which is basically a generalization of stochastic reconnection (LV99) for viscous fluids.

Stochastic reconnection was proposed in LV99 for high *β* plasmas with a magnetic Prandtl number of order unity, $Prm\u22431$. It predicts reconnection speeds of order the large scale turbulent eddy velocity. A modified version of the stochastic model for $Prm>1$ in high *β* plasmas, developed by Jafari *et al.* (2018), showed that both the width of the matter outflow layer and the ejection velocity are unaffected by viscosity. Nevertheless, if $Prm>1$, viscosity may suppress small scale reversals near and below the viscous damping scale. In that case, there will be a threshold for the suppression of large scale reconnection by viscosity when the Prandtl number is larger than the root of the Reynolds number, that is, $Prm>Re$. For $Prm>1$, this leads to the spectral index $\u223c\u22124/3$ for length scales between the viscous dissipation scale and eddies larger by roughly $Prm3/2$. Here, we briefly review this regime of a large magnetic Prandtl number studied in Jafari *et al.* (2018) (see also Jafari and Vishniac, 2018a).

The parallel wavelength at the viscous damping scale, $\lambda \u2225,d=k\u2225,d\u22121$, can be obtained noting that $k\u22a5,d\u22121\u223cl\u22a5(\nu /l\u22a5uT)3/4$. The latter vertical damping scale comes from the Kolmogorov scaling for the viscous damping length scale, $ld\u223c\u03f5\u22121/4\nu 3/4$, with $\u03f5\u2243uT3/l\u22a5\u2243uT2VA/l\u2225$, where we have also used $l\u22a5VA=l\u2225uT$. If we define the injection velocity as $uL=uTVA$, then we can write $\u03f5\u2243uT2VA/l\u2225=uL4/(l\u2225VA)$.

Analytical (see, e.g., Lazarian *et al.*, 2004; Jafari *et al.,* 2018) and numerical (Cho *et al.*, 2002b, 2003a; Schekochihin *et al.,* 2005) work indicate that a Prandtl number larger than unity will affect MHD turbulence in major ways. In MHD turbulence, there are two dissipation scales; the viscous damping scale *l _{d}* and the magnetic diffusion scale,

*l*. In fully ionized collisionless plasmas, both of these dissipation scales are very small. In partially ionized plasmas, they can be very different but always we have $ld\u226blm$ (Jafari

_{m}*et al.*, 2018; Lazarian

*et al.*, 2004). The kinetic energy is dissipated at a viscous damping scale, $ld\u223c\nu 3/4\u03f5\u22121/4$

where we have used the sub-Alfvénic energy transfer rate $\u03f5\u223cuL4/(l\u2225VA)\u2243uT2VA/l\u2225$ ($uL=uTVA$) with critical balance conditions (Goldreich and Sridhar, 1995; Lazarian and Vishniac, 1999; Jafari *et al.*, 2018) $uTl\u2225\u2243l\u22a5VA$. However, at this scale, magnetic energy does not dissipate as resistivity is smaller than viscosity by assumption. Instead, the magnetic damping scale is given by

The width of the reconnection layers, created by the viscosity damped turbulence, can be determined using the eddy size at the viscous damping scale *l _{d}*. Consequently, the width of such a reconnection layer would be roughly $k\u22a5,d\u22121=ld=\lambda \u22a5,d$. The corresponding parallel length at the damping scale, $\lambda \u2225,\u2009d$, can be obtained as

Below *l _{d}*, the energy cascade rate $\tau cas\u22121$ is scale independent, and thus, the constant magnetic energy transfer rate $bl2\tau cas\u22121$ leads to $bl\u223cconst.$

The above expression is the magnetic power spectrum by the assumption that the curvature of field lines is almost constant $k\u2225\u223cconst.$. This is consistent with a cascade driven by repeated shearing at the same large scale. Also, the component of the wavevector parallel to the perturbed field does not change much, and therefore, any change in *k* comes from the third direction (Jafari *et al.*, 2018).

To find the kinetic spectrum, we may use the filling factor $\varphi l$ as the fraction of the volume that contains strong field perturbations with a scale $k\u22121$ (Lazarian *et al.*, 2004). Assuming a constant filling factor, $\varphi l\u223cconst.$, the kinetic spectrum is given by $Ev(k)\u223ck\u22125$. Otherwise, we find $Ev(k)\u223ck\u22124$ (Jafari *et al.*, 2018; Jafari and Vishniac, 2018a).

To find the reconnection speed in a large magnetic Prandtl number regime, Jafari *et al.* (2018) used the following equation for the diffusion of field lines (see also Jafari and Vishniac, 2018a)

where $a=(l\u2225/l\u22a52),\u2009b=l1/3(VA/uT)2/3$, and $c=l1/2(uTVA)\u22121/4(VA/uT)3/4$, and finally, $dx=VAdt$. In order to compare to numerical data, we can use Eq. (76) with calibrated coefficients $a=a1l\u2225/l\u22a52,\u2009b=b1l\u22251/3(VA/uT)2/3$, and $c=c1l\u22251/2(uTVA)\u22121/4(VA/uT)3/4$. The numerical factors $a1=1.32,\u2009b1=0.36$, and $c1=0.41$ are constants introduced to match the solution to the available numerical data obtained by Kowal *et al.* (2012) in which the injection power $P=uT2VA/l\u2225$ and the Alfvén speed *V _{A}* are set to unity and the injection wavelength $kinj=l\u2225\u22121=8$. For these simulations, numerically, we have $a=a1/4$ and

*b*=

*b*

_{1}and

*c*=

*c*

_{1}. Figure 26 shows a comparison of the solution of Eq. (76) with the numerical data obtained by Kowal

*et al.*(2012). As mentioned before, we combine the solution $y=y(t)$ of Eq. (76) with the conservation of mass; $VR=Veject(\u27e8y2\u27e9)1/2/Lx$, where

*V*is the ejection speed of plasma out of the reconnection zone.

_{eject}For $Prm<1$, one expects $Veject\u223cVA$. However, for $Prm>1$, we may expect that the outflow is affected by viscosity if the viscous time scale is shorter than the outflow time scale, $\u27e8y2\u27e9/\nu <Lx/Veject$. Using $\u27e8y2\u27e9/l\u22a52\u223cLx/l\u2225,\u2009VAl\u22a5\u2243l\u2225uT$ and defining the Reynolds number as $Re=l\u22a5uT/\nu $, where *u _{T}* is the velocity characteristic of largest strongly turbulent eddies, one finds

This condition can hardly be satisfied; hence, viscosity is not expected to affect the outflow speed (Jafari *et al.*, 2018). However, it can still affect the reconnection rate by changing the outflow width. In fact, viscosity can affect the reconnection process either by hindering the outflow from the reconnection zone or by narrowing the reconnection channel width (Jafari and Vishniac, 2018a). If we assume that the perpendicular separation is completely inside the turbulence inertial cascade, then viscosity cannot not play any important role. That is if $y>ld$ even for *x* = *L _{x}*, the viscosity's effect will be negligible. In the dissipative range, the diffusion equation (76) has the solution

The condition $y\u2272ld$ then becomes

where we have used $Re=l\u22a5uT/\nu $. Consequently, the outflow width, $y=\u27e8y2\u27e91/2$, is unaffected unless the magnetic Prandtl number is exponentially larger than the Reynolds number.

For $Prm\u226b1$, one generally expects that viscosity will slow down the reconnection speed. Figure 26 compares the general solution of Eq. (76) with the numerical simulations performed by Kowal *et al.* (2009) and Kowal *et al.* (2012). Below the damping scale, Eq. (76) yields

which, for $x\u226b\lambda \u2225,d$, becomes Lyapunov exponential growth whose exponent is inversely proportional to the square root of the viscosity. Figure 27 illustrates this exponential growth of field line separation in a set of numerical simulations performed by Lazarian *et al.* (2004).

Another way that viscosity may affect reconnection is through changing the local small scale reversals. It is a very stringent condition for viscosity to affect the large scale outflow. This may imply that viscosity will almost never be important; thus, the question arises that if the small scale effect of viscosity can be significant.

Another effect of viscosity is associated with topologically distinct local reconnections at small scales. In stochastic reconnection, the global reconnection speed is the collective effect of numerous local reconnection events which occur simultaneously. This sets an upper limit on the global reconnection speed in Lazarian and Vishniac (1999)

where $vrec,eddy$ is the reconnection speed within an individual eddy and *N _{eddy}* is the number of independent eddies along a field line. Jafari

*et al.*(2018) obtained the following condition for viscosity to have an appreciable effect on local reversals

The implication is that for *Pr _{m}* = 1, the above condition is not an interesting limit. However, for

*Re*not very large and resistivity appreciably less than viscosity, it can be the controlling limit on reconnection (and not the outflow zone width). In particular, the dividing line between viscosity dominated reconnection and unimpeded reconnection is no longer overwhelmingly in favor of unimpeded reconnection. Therefore, the critical magnetic Prandtl number does not scale exponentially with the Reynolds number but close to its square root. The simple model outlined above, however, is based on the assumption of a large degree of intermittency below the Kolmogorov scale. This may not be very realistic. In addition, numerical simulations of the viscous regime have currently low resolutions, and their results are affected by a damping scale very close to the eddy scale. The small Reynolds numbers accessible by numerical simulations at present may be the reason why they hardly show any effect associated with viscosity (Jafari

*et al.*, 2018).

The magnetic structures are sheared below the viscous damping scale with a scale independent energy cascade rate, and we find a power spectrum for the magnetic field as $Eb(k)\u223ck\u22121$: Eq. (75). It follows that more power concentrates on the small scales than predicted by the Goldreich–Sridhar spectrum (Goldreich and Sridhar, 1995). With a viscosity of order, or slightly larger than the resistivity, the outflow width remains independent of the small scale physics, which is similar to the LV99 model. Even with viscosity appreciably larger than resistivity, the width of the outflow region is unaffected unless the magnetic Prandtl number is exponentially larger than the Reynolds number. The viscosity is even less important when the current sheet instabilities dominate over external turbulence. For significantly large Prandtl numbers, i.e., of order $Re$, the magnetic field perturbations cannot relax and a flat magnetic power spectrum extends below the viscous damping scale; see Fig. 28.

### B. Turbulent reconnection in relativistic fluid

The relativistic turbulent simulations in Takamoto *et al.* (2015) (henceforth TIL15) confirmed a general similarity between the relativistic reconnection and the nonrelativistic one. They demonstrated that the turbulent reconnection speed can be as large as $0.3c$, which thus enables highly efficient conversion of magnetic energy into kinetic motions and particle acceleration. The numerical results in TIL15 can be understood on the basis of the expression LV99 expression for turbulent reconnection. This expression can be written as

where for sub-Alfvénic turbulent, the energy cascading is given by the following expression (LV99):

However, for the relativistic case, one should take into account both the density change in the relativistic plasma and the modification of turbulence properties in the relativistic regime. The former can be obtained from the energy flux conservation requirement. With both analytical considerations and numerical simulations provided in TIL15, we have

where *β* was found in numerical simulations to be a function of plasma relativistic magnetization parameter $\sigma =B2/4\pi \rho hc2$, where $h=1+(4(p/\rho c2)$ is the specific enthalpy of the ideal relativistic gas, while *ρ* and *p* are the density and pressure, respectively.

The change in the outflow region Δ was shown to correspond to the original LV99 calculation but with injection *ϵ _{inj}* reduced compared to the value in the nonrelativistic case. This corresponds to the transfer of a larger fraction of energy to the fast modes which are subdominant in inducing magnetic field wandering, as indicated by simulations in Takamoto and Lazarian (2016, 2017). With these modifications, the corresponding expression of the turbulent reconnection can be written as

where $\alpha <1$ is the factor accounting for the decrease in the fraction of magnetic energy that induces magnetic field wandering.

Figure 29 illustrates the results of the numerical calculations in TIL15. This study shows that with the suggested modifications, the LV99 theory reflects the major features of turbulent relativistic reconnection.

It is evident from Eq. (86) that the theory of relativistic turbulent reconnection does require further development in order to decrease the uncertainties related to magnetic turbulence in relativistic fluids. For the time being, it is important for our further discussion that qualitatively relativistic turbulent reconnection is similar to its nonrelativistic counterpart. The existing numerical simulations in TIL15 provide us with guidance for studying the reconnection in GRBs. In particular, it is clear from the simulations that in high *σ* flows, the turbulent reconnection speed only slowly changes with the injection velocity and does not depend on the guide magnetic field, i.e., the common field component shared by the two reconnected fluxes. The injection scale of the turbulence induced through the kink instability is likely to be comparable to the scale of the magnetic field flux tubes, i.e., $l\u223cLx$. In this situation, by extrapolating the results in TIL15, one can claim that the reconnection rate is larger than $0.1cA$, where the relativistic Alfvén speed *c _{A}* is very close to the light speed.

### C. Whistler turbulence and reconnection

In a turbulent environment, magnetic field wandering also takes place at scales below the proton Larmor radius. For instance, the observations of the Earth's magnetosphere are frequently on scales less than the ion gyroradius. At such scales, electron MHD or the EMHD approximation can be relevant for describing turbulence and reconnection (see Mandt *et al.*, 1994).

The study of EMHD was reported more than 30 years ago (Kingsep *et al.*, 1987), and the turbulent spectrum $E(k)\u223ck\u22127/3$ was explored more than 20 years ago (Biskamp *et al.*, 1996; Biskamp *et al.*, 1999). A quantitative description of anisotropy in EMHD turbulence anisotropy was studied by Cho and Lazarian (2004, 2009).

We can derive the required scalings in analogy the approach that we discussed in Sec. II B for the derivation of MHD turbulence scalings. The EMHD approximation applies below the ion inertial length $di=c/\omega pi$, where *c* is the speed of light and *ω _{pi}* is the ion plasma frequency. At these scales, the ions constitute a smooth motionless background and electron flows carry all the current. (The ion inertial length is not the same as the ion Larmor radius but differs by a factor of $\beta \u22121/2$, where

*β*is the plasma

*β*. Here, we will neglect the difference.) As a result

where *n* is the electron density and *e* is the charge, respectively. In this limit, the magnetic field is dragged by the motions of the electrons.

The magnetic energy cascades to smaller scales at a rate

where *b _{l}* is the magnetic perturbation at the scale

*l*and

*v*is the electron velocity. Finally, $tcascade\u2248l/vl$ is the cascading time at the scale

_{l}*l*. Here, the scale

*l*refers to the scale perpendicular to the large scale field. We note that Whistler waves have comparable amounts of energy in components parallel to and perpendicular to the mean field and the parallel current is comparable to the perpendicular current. The parallel wavevector component, $k\Vert $, will in general be much smaller than the perpendicular component, $k\u22a5$, but this does not affect the estimates of the coherence time. The energy cascade argument depends only on the perpendicular component. Substituting

*t*in Eq. (88), one gets $bl\u221dl2/3$. The magnetic field wandering can be described as

_{cascade}where *B*_{0} is the mean magnetic field, *s* is the distance alone the field lines, and *y* is the rms perpendicular distance between two field lines. Substituting the scaling of *b _{l}* and using

*y*for

*l*, one gets the perpendicular deviation of magnetic field lines $y\u223cs3$. This is a larger exponent than the classical Richardson divergence of 3/2, which also applies to MHD turbulence.

Another way to understand this is to use the critical balance condition, which in this context means equating the cascade rate with the Whistler frequency

This implies that

and that the divergence of field lines follows the rule that two field lines separated by an eddy will, on the average, increase their separation by the width of the eddy over a distance comparable to the length of the same eddy.

The speed of reconnection on a scale *l* is

i.e., it actually increases on smaller scales. In terms of the strength of the energy cascade, it is

Due to the limited range of whistler turbulence, we do not expect this type of reconnection to play a big role for large scale astrophysical reconnection. In fact, the existing Hall-MHD simulations of reconnection, which show the spectrum of whistler turbulence, transfer for larger scale reconnection events to the turbulent reconnection in the MHD regime (see Sec. V D). However, we expect this type of reconnection to be present in the events that can be probed by *in situ* measurements in the Earth's magnetosphere. Thus, the numerical testing of Eq. (93) is of importance for understanding the magnetospheric physics.

### D. Special regimes: Summary

The theory of turbulent reconnection was formulated for the nonrelativistic reconnection in the MHD regime. The numerical studies of the relativistic regime show that the reconnection in the relativistic regime is rather similar to that in the nonrelativistic one. This arises from the fact that the properties of Alfvénic turbulence and therefore magnetic field wandering are similar in the two regimes. However, the effects of compressibility and the coupling of Alfvén and fast modes are more important for the relativistic reconnection.

The reconnection at a high Prandtl number demonstrates that the turbulent reconnection is a robust process and its rate is only marginally affected by the changes in the microscopic properties of the fluid. Indeed, the study shows that viscosity does not change the reconnection rate at large scales much larger than the viscous dissipation scale. We also note that the particular study is of practical interest, as high Prandtl number turbulence can represent some regimes of turbulence in partially ionized gas.

Finally, the very idea of turbulent reconnection theory that magnetic field wandering regulates the outflow region and determines the reconnection rate carries over to scales less than the ion gyroradius. The corresponding whistler (or kinetic Alfvén wave) turbulence can also induce turbulent reconnection. Such reconnection may be important when reconnection events at the plasma scales are considered. The numerical testing of the predictions that we provide in this section would be of interest, e.g., for magnetospheric reconnection research.

## IX. ASTROPHYSICAL IMPLICATIONS OF TURBULENT RECONNECTIONS

### A. Nonlinear turbulent dynamo and turbulent reconnection

Turbulent reconnection plays an important role for the theory of generating magnetic fields by turbulent motions, i.e., the theory of turbulent dynamo. The turbulent dynamo has been recognized as the most promising mechanism to account for the cosmic magnetism (Brandenburg and Subramanian, 2005). The dynamo process operating on the length scales smaller than the driving scale of turbulence, the so-called “small-scale” turbulent dynamo (SSTD), is especially important for the amplification and maintenance of the magnetic fields in the ISM and intracluster medium. Earlier studies on the SSTD have been restricted to the kinematic regime, where the magnetic field is dynamically insignificant and its effect on turbulent motions is negligible (Kazantsev, 1968; Kulsrud and Anderson, 1992). In the kinematic regime, the dynamo growth of magnetic fields can be analytically described by using the Kazantsev theory and considering the microscopic diffusion of magnetic fields, including the resistive diffusion and ambipolar diffusion (Schekochihin *et al.*, 2002; Kulsrud and Anderson, 1992; Xu and Lazarian, 2016; Xu *et al.*, 2019a). However, the kinematic dynamo gives rise to the magnetic fields with the correlation scale comparable to the dissipation scale of magnetic fluctuations. This cannot explain the observed magnetic fields with the magnetic energy concentrated at a large correlation scale in diverse astrophysical media (Leamon *et al.*, 1998; Vogt and Enßlin, 2005; Beck *et al.*, 2016). The nonlinear dynamo, as confirmed by numerical simulations, accounts for not only the amplification of magnetic field strength but also the increase in its correlation length and thus is astrophysically important.

Analytical studies on the nonlinear dynamo are challenging because of the strong back-reaction of magnetic fields on turbulent motions. Numerical simulations reveal that the nonlinear dynamo has a universal linear-in-time growth of magnetic energy and the growth rate is only a small fraction of the turbulent energy transfer rate (Ryu *et al.*, 2008; Cho *et al.*, 2009; Beresnyak *et al.*, 2009; Beresnyak, 2012). These numerical measurements can be used for testing the validity of analytical theories of nonlinear dynamo.

The efficiency of nonlinear dynamo is determined by the interaction between the dynamo stretching of magnetic fields by turbulent motions and the diffusion of magnetic fields. Earlier theoretical attempts to formulate the nonlinear dynamo focused on the microscopic diffusion of magnetic fields, which is applicable to the kinematic regime. For instance, for the nonlinear dynamo in a conducting fluid, Schekochihin *et al.* (2002) considered the resistive diffusion as the dominant diffusion mechanism and concluded that the resulting magnetic energy spectrum peaks at the resistive scale. In the case of a partially ionized medium, Kulsrud and Anderson (1992) adopted the ambipolar diffusion as the dominant diffusion mechanism and found that the magnetic fields amplified by the nonlinear dynamo have the spectral peak at the ambipolar diffusion scale. More recently, Schober *et al.* (2015) introduced the “effective magnetic diffusion” (Subramanian, 1999) to replace the microscopic diffusion used in the nonlinear dynamo. As a result, the magnetic energy spectrum peaks at the effective resistive scale. However, the theoretical predictions of these models on many characteristics including the dynamo efficiency, magnetic energy spectrum, and correlation length of magnetic fields are inconsistent with numerical findings.

As shown in dynamo simulations, MHD turbulence is developed as a consequence of nonlinear dynamo with a Kolmogorov spectrum of magnetic fluctuations (Brandenburg and Subramanian, 2005; Beresnyak, 2012). We discussed earlier that the violation of flux freezing predicted by the turbulent reconnection theory and supported by numerical simulations (see Secs. III B, IV A, and V B) in most cases exceeds orders of magnitude the magnetic diffusion that directly follows from nonideal effects, e.g., resistivity and plasma effects. Therefore, the process that was termed in Lazarian (2005) the reconnection diffusion (RD) (see Sec. IX B) dominates over microscopic diffusion processes in MHD turbulence. Naturally, one should include the RD to properly model the nonlinear dynamo, and this was done in Xu and Lazarian (2016). The latter paper analytically derived the evolution of magnetic energy $E$ for the nonlinear dynamo

where *ϵ* is the scale-independent turbulent energy transfer rate. The linear dependence on time comes from the constant energy transfer rate of both hydrodynamic turbulence and MHD turbulence. The small fraction 3/38 comes from the RD, which causes the inefficiency of nonlinear dynamo, and quantitatively agrees with numerical measurements (Cho *et al.*, 2009; Beresnyak, 2012). Note that in the above equation, 3/2 is used as the slope of the growing magnetic energy spectrum (Kazantsev, 1968). If a much steeper slope 4 is used (Eyink, 2010; Kraichnan and Nagarajan, 1967), then an even smaller dynamo growth rate is expected.

Due to the RD in MHD turbulence, the magnetic energy spectrum peaks at the driving scale of MHD turbulence, where the turbulent energy and magnetic energy are in equipartition. We note that the effect of magnetic reconnection on turbulent dynamo was also discussed in Kulsrud and Anderson (1992). They adopted Petschek's model for reconnection (Petschek, 1964), and thus, the peak scale of the magnetic energy spectrum is much less than that in Xu and Lazarian (2016), which cannot account for the large correlation scale of observed magnetic fields.

This peak scale increases with the growth of magnetic energy (Xu and Lazarian, 2016)

Here, *k*_{0} is the wavenumber where the spectrum peaks at the beginning of nonlinear dynamo at *t* = *t*_{0}. On smaller length scales, magnetic fields are dynamically important and result in the scale-dependent anisotropy of turbulence. On the other hand, the RD allows efficient diffusion and thus turbulent motions of magnetic fields. Therefore, instead of having a folded structure with magnetic field reversals only at the dissipation scale, the magnetic fields amplified by the nonlinear dynamo are turbulent magnetic fields with field reversals on all length scales within the inertial range of turbulence. Figure 30 illustrates a comparison between the analytical scalings of the magnetic energy spectrum obtained by Xu and Lazarian (2016) and numerical measurements in Brandenburg and Subramanian (2005) in cases with different magnetic Prandtl numbers *P _{m}*. We see that the nonlinear dynamo theory based on the RD well explains the behavior of nonlinear dynamo presented in numerical simulations.

The inefficiency of nonlinear dynamo has important astrophysical implications for the evolution and importance of magnetic fields in both early and present-day universe. In earlier studies on the turbulent dynamo during the primordial star formation, e.g., Schober *et al.* (2012), the RD was not taken into account. This leads to an efficient nonlinear dynamo and generation of dynamically significant magnetic fields over a short time scale. Xu and Lazarian (2016) applied the nonlinear dynamo theory including the RD and examined the dynamo evolution of magnetic fields in environments like the first stars. As shown in Fig. 31, among different evolutionary stages, the nonlinear dynamo has the longest time scale due to the effect of RD. Since the resulting dynamo time scale is longer than the system's free fall time, it is unlikely that the dynamo-amplified magnetic fields can regulate the gravitational collapse.