Dynamic alignment and plasmoid formation in relativistic magnetohydrodynamic turbulence

We present high resolution 2D and 3D simulations of magnetized decaying turbulence in relativistic resistive magneto-hydrodynamics. The simulations show dynamic formation of large scale intermittent long-lived current sheets being disrupted by the tearing instability into plasmoid chains. These current sheets are locations of enhanced magnetic field dissipation and heating of the plasma. We find magnetic energy spectra $\propto k^{-3/2}$, together with strongly pronounced dynamic alignment of Elsasser fields and of velocity and magnetic fields, for strong guide-field turbulence, whereas we retrieve spectra $\propto k^{-5/3}$ for the case of a weak guide-field.


INTRODUCTION
Turbulence provides a route for the energy cascade and dissipation in a wide range of astrophysical plasmas. This is relevant for astrophysical systems like black hole accretion disk-jet systems (e.g., Ripperda et al. 2020Ripperda et al. , 2021Mahlmann et al. 2020), magnetar magnetospheres (Beloborodov 2020) and pulsar wind nebulae (e.g., Lyubarsky 1992;Begelman 1998). These astrophysical systems are typically relativistic, meaning that the magnetization σ = B 2 /(4πωρc 2 ) ≥ 1, where B is the magnetic field strength, ρ is the plasma density, and ω is the relativistic enthalpy density, indicating that the magnetic energy density is larger than the plasma energy density. This results in an Alfvén speed v A = σ/(σ + 1)c that is close to the speed of light c.
Most turbulence studies have been in the realm of nonrelativistic magnetohydrodynamics (MHD) when the Alfvén speed, v A , is much lower than the speed of light, c. Iroshnikov (1963); Kraichnan (1965) showed that the energy cascade from large to small scales is caused by the mutual shear of counter-propagating Alfvén waves.
Corresponding author: Alexander Chernoglazov alexander.chernoglazov@gmail.com * Joint Princeton/Flatiron Postdoctoral Fellow Thirty years later Goldreich & Sridhar (1995, 1997 suggested that turbulent systems are in the critical balance regime meaning that an eddy is significantly deformed during one Alfvén-crossing time. This also means that the turbulent eddies are elongated along the background magnetic field. The first steps towards a theory of relativistic turbulence were taken recently by Chandran et al. (2018), and they demonstrated that the relativistic picture is very similar to the Newtonian limit (more details are presented in Section 2). Boldyrev (2005Boldyrev ( , 2006 suggested that turbulent eddies are anisotropic in all three directions: they are elongated along the guide magnetic field and have two different sizes in the guide field-perpendicular plane. The ratio of these two sizes is called dynamic alignment angle. These eddies are progressively more elongated at smaller scales. Recent theories (e.g., Boldyrev & Loureiro 2017;Mallet et al. 2017) proposed that the elongated eddies at small enough scale become unstable to the tearing instability, causing a steepening of the turbulence spectrum.
In their recent paper, Dong et al. 2018 demonstrated the formation of reconnecting current sheets in twodimensional (2D) decaying non-relativistic turbulence. They also demonstrated the formation of a turbulence spectrum and dynamic alignment in agreement with Boldyrev's theory. It is however as of yet unclear whether these findings persist in the case of realistic three-dimensional (3D) turbulence. In 3D MHD, despite prominent current sheet formation (Zhdankin et al. 2013), it remains unclear whether reconnection can occur in the fast regime, when the dissipation efficiency is independent of resistivity. This regime is associated with the formation of plasmoid chains, resulting in a universal reconnection rate of order 0.01 (Bhattacharjee et al. 2009;Uzdensky et al. 2010).
Plasmoid-mediated reconnection in relativistic plasmas can accelerate particles to non-thermal energies (e.g., Sironi & Spitkovsky 2014;Guo et al. 2014;Werner et al. 2015), responsible for the high-energy emission in many environments of compact objects (e.g., Cerutti et al. 2015;Beloborodov 2017). Recent studies of relativistic turbulence in collisionless plasmas have shown efficient particle acceleration (Zhdankin et al. 2017;Comisso & Sironi 2018) and the formation of reconnecting current sheets, which are important for the process of initial particle acceleration (Comisso & Sironi 2018) both in 2D and 3D. The high-energy power-law tail of the distribution function has been shown to get steeper quickly for smaller ratios of the turbulent component of the field to the guide field at the outer scale, δB/B 0 . This observation further motivates the exploration of current sheet properties at moderate δB ∼ B 0 , when particle acceleration is efficient.
The highly magnetized relativistic (v A ≈ c) MHD limit has been largely unexplored, and it is unclear whether dynamic alignment forms in this regime, and whether it plays an important role for the current sheet formation for situations where δB ∼ B 0 . Neither any presence of dynamic alignment nor plasmoid unstable current sheets were shown in the first relativistic ideal MHD simulations by Zrake & MacFadyen (2012).
In this Letter we present numerical relativistic resistive MHD simulations of decaying turbulence in highlymagnetized plasma both in 2D and 3D. We demonstrate that dynamic alignment forms both in 2D and 3D. We show intermittent long-lived current sheets form naturally in the turbulence and become plasmoid-unstable.

THEORETICAL OVERVIEW
The study of non-relativistic turbulence is usually done with a reduced MHD approach. This method employs a few assumptions: a uniform, strong, in comparison to the perturbation δB, guide field B 0 and incompressibility of the flow (c s → ∞, where c s is the sound speed). Under these assumptions, the only waves of interest are perpendicularly-polarized Alfvén waves, propagating along the guide field. The reduced form of MHD equations in this limit reads (Elsasser 1950): (1) where P is the total pressure, and z ± = δv ± δB/ √ 4πρ 0 are the Elsasser fields, representing counter-propagating Alfvén waves.
The ideal relativistic MHD equations consist of mass and stress-energy conservation laws, and the induction equation for the magnetic field evolution: Here µ, ν are 4-dimensional space-time indices, such that u µ = (Γ, Γv) is the four-velocity vector, Γ is Lorentz factor, and T µν is the stress-energy tensor with E = ρωc 2 + b 2 and ω = 1 + (γ/(γ − 1))P/ρ is the relativistic enthalpy, η µν = diag{−1, 1, 1, 1}, the flatspacetime Minkowski metric. b µ is the magnetic field four-vector and b 2 = b µ b µ . Introducing the relativistic Elsasser fields and modified pressure term, Π = (2P +b 2 )/(2E), one can rewrite the relativistic MHD equations in the Elsassertype form: (Chandran et al. 2018;TenBarge et al. 2021) In contrast to the non-relativistic case, one cannot formally introduce an incompressible limit in relativistic MHD wherein there is a maximum speed of propagation, c. The finite speed of light prevents easy elimination of the fast magnetosonic modes (Takamoto & Lazarian 2017). However, it is still possible to order them out in the highly anisotropic limit, k ⊥ k || , which implies ω F ∼ kc ω A ∼ k || c, where ω F and ω A are the frequencies of the fast magnetosonic and Alfvén modes, correspondingly (TenBarge et al. 2021). Here, we are interested in the highly-magnetized limit, where the hot magnetization parameter, σ, is large: This implies b 2 P , or Π = 1/2, and the elimination of the slow magnetosonic mode. The second term of the relativistic Elsasser variable z µ ± is a unit vector in the direction of the four-magnetic field vector as E ≈ √ b 2 in this limit. TenBarge et al. (2021) discusses the Elsasser-type equations of the relativistic MHD in the highly anisotropic limit, k ⊥ k , constructed in the average fluid rest frame u i = 0 (Chandran et al. 2018). In three-vector form the result is particularly similar to the reduced non-relativistic MHD equations: Due to the close resemblance between the Newtonian and relativistic set of reduced MHD equations, once the anisotropic cascade reaches sufficiently small scales, where k ⊥ k is satisfied, one can expect that relativistic MHD turbulence is statistically similar to its Newtonian counterpart. In the case of interest, σ 1, the last term of (8) is negligible, and the Alfvén speed is close to the speed of light, v A ≈ c, such that the equations are particularly simple: The applicability of equation (9) is limited to regimes where δP δB 2 . However, when plasma is heated to relativistic temperatures, e.g., in reconnection layers, this assumption is not justified, at least locally. Since there is no a formally incompressible limit in relativistic systems, and our interest in systems with δB ∼ B 0 , which is particularly challenging to explore analytically, we turn to numerical simulations to confirm these expectations.
An important feature of the highly magnetized MHD turbulence is the overall dominance of the magnetic and electric field fluctuations, δE B , over the kinetic energy, δE kin . This can be seen from the following relations for a single Alfvén wave in the relativistic regime (σ 1, v A ∼ c): where the inequality is used, Current sheets are important dissipative structures in magnetized turbulence, and it is useful to compare their behavior in non-relativistic and relativistic regimes. In a near-stationary current sheet, the reconnection rate is the ratio of the inflow velocity to the outflow velocity v in /v out (Parker 1957;Sweet 1958). If the plasma can be assumed incompressible, it then follows that v in /v out = δ/L where δ is the thickness and L is the length of the current sheet. The thickness of a Sweet-Parker current sheet is determined from continuity of the resistive and ideal electric fields as δ = η/v in , where η is the resistivity. The outflow speed can be approximated as the Alfvén speed v out ∼ v A , which in a relativistic plasma is v A ∼ c. This results in a relativistic reconnection rate of v in /v out ∼ v in /c ∼ η/(Lv in ), i.e., a result identical to the non-relativistic case (Lyubarsky 2005). The reconnection rate in a Sweet-Parker sheet then scales as ∼ S −1/2 , where S = Lv A /η is the Lundquist number. For large Lundquist numbers, typical in astrophysical sources, reconnection is mediated by the plasmoid instability, which in non-relativistic settings gets triggered at S crit ≥ 10 4 , leading to a saturation of the reconnection rate at ≈ 0.01 (Loureiro et al. 2007;Bhattacharjee et al. 2009). It was shown semi-analtyically and numerically, by solving the full set of resistive relativistic MHD equations, that this result holds for highly magnetized relativistic plasmas (Del Zanna et al. 2016;. Since magnetic field fluctuations dominate over kinetic energy, resistive dissipation dominates over viscous dissipation in the inertial range of turbulence. We show that resistive dissipation in highly-magnetized MHD plasmas is also dominant in the exhausts of reconnection layers. One can see this by comparing the rate of resistive energy dissipation ∼ ηB 2 /4πδ 2 exhaust , to the rate of viscous dissipation ∼ 2νE kin /δ 2 exhaust , where ν is the viscosity, E kin is the kinetic energy in the exhaust region, B 2 /8π is magnetic energy and δ exhaust is the typical width scale in the exhaust region. The ratio of resistive to viscous dissipation rates is then η/ν ·(B 2 /E kin ) ∼ (B 2 /E kin ) 1 . In a non-relativistic plasma, this ratio is always of order O(1) since E kin ∼ ρv 2 A and v 2 A ∼ B 2 /ρ. However, in highly magnetized relativistic plasma, E kin ∼ ργ exhaust c 2 ∼ ρ √ σc 2 (Lyubarsky 2005) while B 2 ∼ ρσc 2 , and hence the ratio between resistive and viscous dissipation rates 1 Our argument is applicable for the case of a scalar viscosity ν, i.e., independent of the gas pressure which is high in the exhaust region. Whether this assumption is accurate for the effective viscosity of a collisionless relativistically hot plasma is an interesting question for further investigation.
is proportional to √ σ 1 such that resistive dissipation dominates.

NUMERICAL METHOD AND SETUP
We solve the set of special relativistic resistive MHD (SRRMHD) equations with the Black Hole Accretion Code (BHAC, Porth et al. 2017;Olivares et al. 2019) and an Implicit-Explicit (IMEX) time stepping scheme to evolve the stiff resistive Ohm's law . We employ a constant and uniform resistivity η, which provides the simplest prescription to allow resolved magnetic reconnection.
The SRRMHD equations are numerically evolved in a periodic domain of size L 2 in 2D and L 2 × L z in 3D. We initialize an out-of-plane (in theẑ-direction) guide magnetic field and an in-plane (x − y) magnetic field perturbation δB ⊥ : , k m = 2πm/L, and φ mn , ϕ mn are random phases. We set N = 8 initial waves in each direction for 2D runs (64 initial modes) and N = 4 for 3D runs, in order to allow for a larger inertial range in 3D simulations. The outer (or energy containing) scale is then l ⊥ = L 0 /8 for 2D simulations and l ⊥ = L 0 /4 for 3D. The turbulence at smaller scales forms self-consistently via energy cascading. In 3D we modulate δB ⊥ with two modes ∝ sin(k l z +ψ mnl ), where ψ mnl is also a random phase. The normalization coefficient is . We initialize the plasma at rest, with velocity field v = 0, and with a uniform gas pressure p 0 and rest mass density ρ 0 . We set an adiabatic index γ = 4/3, assuming an ideal relativistic gas. Similar initial conditions for the magnetic field were employed in relativistic particle-in-cell (PIC) (Comisso & Sironi 2018;Nättilä & Beloborodov 2020) turbulence simulations. For all simulations we set L = 1.
In order to characterize the strength of both the guide and the in-plane magnetic field we introduce two magnetization parameters A summary of the performed runs is given in Table  1. We employ an elongated box with L = 1, L z = 3 for run 3D[d] to enforce the critical balance condition δB ⊥ /L ≈ B 0 /L z at the outer scale. We set the resistivity to either η = 10 −5 , 10 −6 in the 2D setup, which corresponds to Lundquist numbers S ≈ 10 4 , 10 5 for the largest current sheets of length L cs ≈ 0.1, and v A /c ≈ 1. The simulation with η = 10 −6 is well-above the critical Lundquist number limit S crit , while the simulation with η = 10 −5 is approximately at the limit S ≈ S crit . Potentially, current sheets can become plasmoid unstable at a smaller critical Lundquist number in a turbulent flow (Loureiro et al. 2009). We explore whether this effect is significant in 2D relativistic turbulence with our η = 10 −5 simulation.
In order to ensure that the resistive length scales are resolved, and results are converged with numerical resolution, we develop a novel adaptive mesh refinement (AMR) strategy (see Appendix A). We benchmark our 2D results with a short simulation (until t = 0.5L/c) on a uniform grid, with a resolution of 65536 2 . In 3D it is impossible to fully converge due to numerical limitations, and instead we employ the highest feasible resolution that allows to capture the development of the plasmoid instability in the longest current sheets. One high-resolution run is performed with 3200 3 grid points to probe the formation of plasmoid chains. We additionally present a study with different values of the magnetization parameter at a resolution of 2048 3 grid points.
The SRRMHD algorithm relies on viscosity ν at the grid level, such that the magnetic Prandtl number Pr m = ν/η 1 for 2D simulations, and Pr m 1 for 3D simulations with a marginally resolved resistive scale, assuming resistive and viscous scales are similar and governed by the finite grid resolution. This choice is further motivated by the fact that viscous effects are subdominant for highly-magnetized plasmas, as we have demonstrated in Section 2.

RESULTS
We present the results of simulations which we run until t = 2L/c, such that the turbulence is fully developed and settles to a quasi-steady state. For the case of decaying turbulence considered in this Let- Table 1. Summary of simulation parameters.

Sim
Res 1 1 2 × 3 Uni ter, we define a quasi-steady state when the spectral slope is constant for at least one outer scale eddy turnover time, ∼ L/c, while the total energy Here, B k is the amplitude of the Fourier mode of the magnetic field with wavenumber k. It takes ∆t ≈ 0.3L/c for the energy to cascade from the initial low k modes to the resistive scale. At ∆t ≈ 1 − 2L/c the spectrum flattens until it reaches a quasi-steady state. This behavior of the power spectrum is illustrated in the attached video 2 . In order to test convergence of the simulation, we compare spectra of magnetic energy, E(k)dk = k∈dk B k · B * k /8π, for different resolutions: if the onset of the inertial range cutoff does not change with increasing resolution (the vertical lines in Figure 5c), i.e., if the cutoff is determined by the resolved resistive scale, the simulation is considered converged. More details about the AMR strategy and convergence tests are presented in the Appendix A.
In Figure 1a we present the distribution of the outof-plane electric current density j z ∼ (∇ × B) z for a 2D simulation with an effective resolution of 65536 2 grid points, 2D[a], and 3 AMR levels, for a resistivity η = 10 −6 . Here, we set magnetizations σ 0 = δσ = 5, equivalent to a total magnetization σ = 10. Very long current sheets emerge at the interfaces of large merging eddies. The length of a current sheet is mainly defined by the size of the largest eddies present in the system. Estimating the length of these current sheets as L sheet /L ≈ 0.1 and accounting for the relativistic Alfvén speed v A /c = σ/(σ + 1) ≈ 1, we find the Lundquist number to be S ≈ 10 5 S crit = 10 4 . These current sheets are plasmoid-unstable and break up into current sheets of smaller length scales, such that their Lundquist numbers S local ≈ S crit . This results in a maximum number of ∼ 10 plasmoids, which is consistent with results shown in Figure 1a. We also perform simulations for σ 0 = δσ = 1 and σ 0 = 1, δσ = 5 and resistivity η = 10 −6 . Plasmoid-unstable current sheets form ubiquitously for all of these settings. For resistivity η = 10 −5 we observe only few plasmoids (for the longest current sheets) in the whole domain indicating that the critical Lundquist number S crit ≈ 10 4 holds for the plasmoid instability in current sheets in a 2D turbulent flow. By varying numerical resolution, we find that the onset of the plasmoid instability occurs at lower resolution for the cases with weaker guide field, σ 0 ≤ δσ, motivating our choice to perform our highest resolution 3D simulation for σ 0 /δσ = 1/5 (run 3D[a]).
2D and 3D weak guide field (3D[a]) simulations show pronounced reconnection-mediated mergers of smaller eddies. This process has also been recently observed in simulations of merging non-helical flux tubes (Zhou et al. 2020). Our very long 2D simulations with a (smaller) resolution of 3200 2 demonstrate that the terminal state of the turbulence has two large eddies of opposite magnetic helicity A · Bdx remaining in the simulation box. 3D simulations show similar behavior.
In order to identify current sheets, we choose a threshold in the current density, ξ, and consider a point x to be in the current sheet if j z (x) > ξj rms , where j rms is the root-mean-square of the electric current j z in the domain (Zhdankin et al. 2013). The long current sheets have an intermittent nature and occupy about 0.2-0.5% of the domain in 2D, yet they are responsible for 20 − 25% of the magnetic field dissipation, ∝ ηj 2 , for η = 10 −6 . Our results are insensitive to the exact value of ξ, as long as ξ 5. For a larger value of resistivity, η = 10 −5 , only few plasmoids form in the whole simulation box (see Figure 1b). Comparing to the η = 10 −6 case, the current sheets for η = 10 −5 are thicker and, hence, have a lower current density amplitude. In this case we find only 10% of the dissipation to happen in the localized current sheets.
The anisotropic properties of the turbulence can be quantified by measuring the dynamic alignment angle of eddies in the plane perpendicular to the guide field (Boldyrev 2005(Boldyrev , 2006.
We employ a Monte-Carlo method to compute dynamic alignment angle as a function of a point-separating vector.
Following the method proposed by Mason et al. (2006); Perez et al. (2012), we compute two structure functions S 1 where δv ⊥ (l) and δB ⊥ (l) are increments of the velocity and the magnetic field perpendicular to the local guide field at scale l. The alignment angle is defined as where l and ξ are sizes of the eddy in the guide fieldperpendicular plane. Note that both S 1 1 and S 1 2 are first order structure functions S n (l) = |f (r + l) − f (r)| n , and one can define the alignment angle for any n as θ ≈ (S n 1 (l)/S n 2 (l)) 1/n . It turns out that the slope of the alignment angle function is dependent on the order of the structure function n (Mallet et al. 2016) revealing the intermittent nature of dynamic alignment (Frisch 1995). We present more details on the intermittency of dynamic alignment in Appendix C. The slope of the dynamic alignment angle as a function of the size of eddies is tightly connected to the power-law of the magnetic energy spectrum, P B (k), which is predicted to be k −3/2 ⊥ for turbulence that is anisotropic in the plane perpendicular to the guide field (Boldyrev 2006). By introducing a non-linear time, τ c (Boldyrev 2005), we can relate the power spectrum with the alignment angle: where ε is the energy cascading rate, δB l is the increment of the magnetic field at a scale l in the plane perpendicular to the guide field. For sin θ(l) ∼ l 1/4 ∼ k −1/4 ⊥ , it reproduces Boldyrev's spectrum E(k ⊥ ) ∼ k −3/2 ⊥ . In the case of no alignment being present, θ(l) ∼ const, it reduces to the Goldreich & Sridhar (1995) In Figure 1c we show the magnetic power spectrum in the steady state and multiply the result by k 3/2 ⊥ , to make the difference between the power law indices −3/2 and −5/3 more pronounced. Figure 1c clearly demonstrates that the spectrum is closer to k −3/2 ⊥ in the inertial range. We define the steady state of decaying turbulence when the spectral slope is constant in time, at t ≥ 0.5L/c, allowing us to compare our results with steady state theory. Note that the total energy B (t) = (B 2 − B| 2 k=0 )/8πdV decreases in time while the normalized spectrumẼ(k)dk = k∈dk B k ·B * k /8π b , which we present in all spectrum plots, is constant in time. It is worth mentioning that in non-relativistic reduced MHD simulations one typically analyzes the spectrum of the total kinetic and magnetic energy (Perez et al. 2012). In our highly magnetized relativistic simulations however, , and we confirmed that the contribution of the kinetic energy is negligible for both 2D and 3D simulations (see spectra in Figures 1c and  2c). To preserve the kinetic to magnetic field energy ratio, we also normalize the kinetic energy by b (t). In agreement with the k −3/2 ⊥ power spectrum, the v − B dynamic alignment (Figure 1d) demonstrates a perfect match with Boldyrev's prediction, θ(l) ∼ l 1/4 , at the intermediate scales, l res l l max , where l max = L 0 /8 is defined by the number of modes in the initial conditions, and l res is defined by the resistive scale. Chandran et al. (2015) proposed that mutual shear of counter-propagating Elsasser fields δz ± is responsible for the dynamic alignment. They predict that these two fields create a progressively decreasing alignment angle, while the slope becomes flatter. To test this hypothesis, we measure the alignment angle between two Elsasser fields: Straightforward application of the non-relativistic Elsasser field expression, δz giving that their ratio θ z + ,z − ∼ δv/δB 1 in highly magnetized plasma. However, one should use the relativistic formulation of Elsasser fields (5) in this regime, where u and b/ √ E can be comparable. The dynamic alignment angle between the relativistic Elsasser fields is flatter than l 0.25 at t = 2, for η = 10 −6 ( Figure 1d). The average slope is close to the l 0.1 result, as predicted by Chandran et al. (2015), although it displays an unexpected break at intermediate scales.
The smallest averaged dynamic alignment angle, θ v−B , in the simulation with η = 10 −6 is 0.175, and it is approximately constant for small scales. Deviations from Boldyrev's scaling l 0.25 are visible at scales where resistive effects become important. Note that this is also where the inertial range of the spectra ends. The plasmoid-unstable current sheets we observe in the simulation possess much smaller alignment angles θ ≈ 0.01, in accordance with Loureiro et al. (2007). The presence of such current sheets with alignment angles of an order of magnitude smaller than the minimal averaged alignment angle that we find, implies the intermittent nature of these sheets (Dong et al. 2018). Formation of intermittent plasmoid-unstable current sheets can be responsible for a steepening of the spectrum at the end of the inertial range, which we observe in the range k ⊥ ≈ 300 − 1200 at t = 1 in Figure 1c. However, we assume that the scale separation in our simulations is not enough to robustly confirm the k −11/5 ⊥ prediction by Boldyrev & Loureiro (2017) and Mallet et al. (2017) for the non-relativistic reconnection-mediated regime. We also do not observe the increase of the alignment angle at small scales l corresponding to wave-vectors k ⊥ in the steepening range, as predicted in Boldyrev & Loureiro (2017).
Since the onset of the plasmoid instability occurs at lower resolution in 2D simulations if a weaker guide field is assumed, we run a 3D simulation (3D[a]) with σ 0 = 1, δσ = 5, and highest resolution of 3200 3 grid points. For 2D simulations we confirm that full plasmoid chains form for smaller values of δB ⊥ /B 0 as well, but higher resolutions are required to resolve the instability. We refer to the case with initial δB ⊥ /B 0 = √ 5/1 (run 3D[a]) as a weak guide field, and to the case with initial δB ⊥ /B 0 = 1/3 (run 3D[d]) as a strong guide field. We note that by t = 1 − 2, when we analyze the simulations, the turbulent component of the field decayed to For the strong guide field case, the energy cascade is developing mainly in k ⊥ ⊥ẑ, and the full 3D analysis can be reduced to a 2D analysis in a set of x − y planes (e.g., Perez et al. 2012). For simplicity, in the case of the weak guide field we also compute the spectrum for wavevectors k ⊥ perpendicular to B 0 using the same method (a more accurate calculation would use structure functions which take into account a locally varying guide field, Cho & Vishniac 2000). In order to provide a statistically significant result, we average the 2D spectrum and dynamic alignment angle in the set of x − y planes taken at various z. We confirm that the spectrum and the alignment angles are independent of the choice of the sampling planes if N planes N z /3, where N z is the number of grid points in the direction along z.
We consider the turbulence at t = 2 to be in a steady state, i.e., the dynamic alignment is fully formed (see Figure 2). We confirm the steady state shape of the dynamic alignment angle function beyond t = 2 with longer simulations at a lower numerical resolution, 2048 3 (runs 3D[b], 3D[c]). The slope of the v − B alignment angle is close to the predicted l 0.25 for the smaller eddies and is less pronounced for eddies of the system size scale (see Figure 2d). In simulation 3D[a] with the initially weaker guide field δB ⊥ /B 0 = √ 5/1, at t ≈ 1, the alignment angle curve is significantly shallower, consistent with the steady state in driven non-relativistic turbulence at δB ⊥ /B ∼ 1 (Mason et al. 2006). At this time, the strength of turbulent fluctuations is similar to the strength of the guide field, |δB ⊥ | ≈ |B z | . Further dissipation of the magnetic energy leads to |δB ⊥ | ≈ 0.7 |B z | at t = 2, and a steeper alignment angle curve. The spectrum of the turbulence develops simultaneously with the dynamic alignment. The slope of the z + − z − dynamic alignment angle, θ z + ,z − (l), is comparable to θ(l) for the strong guide field (run 3D[d], t = 2, |δB ⊥ | / |B z | ≈ 0.2). For the weak guide field (3D[a], t = 2, |δB ⊥ | / |B z | ≈ 0.7), the z + − z − alignment is very weakly pronounced. At the same time, the slope of the energy spectrum of 3D[a] is closer to −5/3 as predicted by Goldreich-Sridhar theory with no dynamic alignment. It could be considered as an indication that the dynamic alignment of Elsasser fields δz + , δz − rather than the one of v, B reduces the non-linearity.
3D simulations show less pronounced boundaries of large-scale eddies, but the intermittent large current sheets are still present in the system with the weak guide field. Figure 2a and the linked video 3 demonstrate the distribution of the electric current j z in the planes perpendicular to the guide field, B z , at t = 1, and Figure 2b and the accompanying video 4 show the same at t = 2. Similarly to the 2D results, intense current sheets occupy up to 4−5% of the total volume of the domain 5 and lead to 20% of the total dissipation of the magnetic energy. Intermittent long current sheets are clearly plasmoidunstable as shown by the insets in Figure 2. A few initial eddies are still clearly seen at t = 1, but many long intermittent current sheets are unaffected by the choice of the initial conditions. At t = 2 no visible features are associated with the initial conditions (see Figure 2b).
Overall, the structure of the electric current in the 3D weak guide field simulation looks similar to the one in 2D (Figure 1a): there is a number of well-pronounced long, plasmoid unstable current sheets formed at the outer scale. Their formation is likely associated with the mergers and subsequent reconnection of large coher- ent structures (Hosking & Schekochihin 2020). Unlike in the weak guide field regime, the strong guide field simulation 3D[d] shows the statistical properties of "aligned" critically-balanced turbulence: the k −3/2 ⊥ spectrum and a pronounced dynamic alignment (dashed lines in Figure 2c and d). The spatial distribution of the electric current is more uniform in this case (see Figure 3). The absence of very long current sheets is consistent with the observation of a very few plasmoids in the simulation (see insets of Figure 3). A possible explanation can be found in the small ratio of the length, L sheet ∼ 0.05, for the sheets shown in the insets of Figure 3, to the width of these sheets, which at our resolution is still limited by the numerical diffusion. We anticipate that the plasmoid instability can be more reliably captured at much higher spatial resolution: for the typical length, L sheet ∼ 0.05, and the width-to-length ratio θ ≈ 0.01, one requires (N δ /(L sheet θ)) 3 ≈ 10000 3 grid points, where N δ ≈ 5 cells is the minimally desired resolution per width of the plasmoid-unstable current sheet.
The structure of a representative current sheet for the weak guide field simulation 3D[a] is presented in Figure  4. The volume render represents the current density amplitude, and solid black lines show selected magnetic field lines. The lower threshold for the volume rendering is chosen to be ≈ 2j rms , in order to remove the upstream regions without significant current. The initial (seed) points for the integration of magnetic field lines are set inside two randomly chosen plasmoids. Wrapped, helical magnetic field is responsible for the large current density inside the plasmoids. The helical structure allows longer plasmoids (or, flux tubes) to be kink-unstable if their length is large enough to exceed the Kruskal-Shafranov stability limit. This instability likely limits the life time of plasmoids in current sheets and their axial extension. A zoom into the 3D structure of a plasmoid is shown in Figure 4b.
Acceleration of the flow from the X-point of a Sweet-Parker current sheet up to Alfvénic speed creates a velocity shear, which may become unstable to the Kelvin-Helmholtz instability (KHI). The analytical nonrelativistic stability criterion ∆u v A (Loureiro et al. 2013) suggests that the strong upstream magnetic field can lead to the stabilization of the KHI for the velocity shear ∆u and the Alfvén speed v A determined by the upstream magnetic field strength. A similar criterion was derived for a simplified model with B||v in a fully relativistic case (Osmanov et al. 2008). Thus, we expect current sheets in highly-magnetized turbulence to be stabilized by the strong upstream magnetic field. To confirm this prediction, we conduct localized numerical experiments with conditions inferred from turbulence simulations (see Appendix B for the description of the setups) that confirm that long plasmoid-unstable current sheets are Kelvin-Helmholtz-stable both in 2D and 3D simulations.

CONCLUSIONS
In this Letter we present the first 2D and 3D numerical SRRMHD simulations of highly magnetized decaying turbulence. We calculate statistical properties of the turbulence, by analyzing a quasi steady state at two Alfvén-crossing times of the simulation box. We show that the spectrum of magnetic energy in both cases is close to Boldyrev's spectrum, k −3/2 ⊥ , and the v − B dynamic alignment angle follows an l 1/4 dependence. Despite the dynamic alignment angle of v and B fields in 2D is perfectly following Boldyrev's prediction, its formation cannot be explained by the uncertainty principle originally employed by Boldyrev (2006). On the other hand, intermittent structures are vastly present in the simulations, favoring the theory of mutual shearing of Elsasser fields by Chandran et al. (2015): an in-depth analysis of this approach is presented in Appendix C. We demonstrate that long-lived intermittent current sheets form dynamically throughout the evolution. These sheets are plasmoid unstable and KH-stable. They occupy a very small fraction of the numerical do- main but provide a significant fraction of the total magnetic field dissipation.
In our simulations we only employ explicit resistivity while viscosity is dictated by the finite grid resolution. We expect that the magnetic energy dominates the kinetic energy at all scales, and dissipation is governed by resistivity. It will be useful to perform simulations with explicit viscosity and fixed magnetic Prandtl number Pr m in the future studies, and to consider the transrelativistic regime, σ ∼ 1. These studies can be applied to turbulence in the accretion disk-jet boundary with moderate magnetization (Ripperda et al. 2020).
In order to study the properties of intermittent current sheets in a statistical steady state, it is important to study driven turbulence in highly magnetized plasmas σ 1. High magnetization leads to efficient heating of the plasma due to the dissipation of magnetic energy and a significant drop of σ. To mediate the effect of runaway heating, radiative cooling of the plasma should be incorporated in the simulations (Zhdankin et al. 2021).
The limitation of computational resources does not allow to reach numerical resolutions significantly above 10000 3 in the nearby future. This is too low to reach alignment angles substantially below θ ∼ 0.1 at the smallest scales. On the other hand, the intriguing similarity of statistical properties of 2D and 3D turbulence in our simulations makes it interesting to perform even higher resolution simulations of 2D turbu-lence. The most significant milestone will be a resolution of ∼ (10 8 ) 2 which allows progressively elongated eddies to reach an alignment angle θ ∼ 0.01 corresponding to the plasmoid instability of these eddies. The steepening of the turbulence spectrum due to the onset of the plasmoid instability in intermittent current sheets (or due to the linear tearing instability in elongated eddies, Boldyrev & Loureiro 2017) can be measured reliably at resolutions of ∼ (10 6 ) 2 , realistically attainable in the nearby future, in particular with the AMR criterion we propose here.

ACKNOWLEDGEMENTS
We acknowledge useful discussions with Lev Arzamasskiy, Amitava Bhattacharjee, Benjamin Chandran, Luca Comisso, Mikhail Medvedev, Joonas Nättilä, Jason TenBarge, James Stone, and help in navigating through 3D-visualization by Hayk Hakobyan. The authors acknowledge insightful comments by the anonymous referee which helped to significantly improve the manuscript. The computational resources and services used in this work were provided by facilities supported by the Scientific Computing Core at the Flatiron Institute, a division of the Simons Foundation; and by the VSC (Flemish Supercomputer Center), funded by the Research Foundation Flanders (FWO) and the Flemish Government -department EWI. A.C. gratefully acknowledges support and hospitality from the Simons For the fully resolved and converged 2D simulations we present the adaptive mesh refinement (AMR) criterion we designed to accelerate the simulations and to simultaneously capture the main properties of the turbulence and dissipative structures. The main principle is that the largest eddies are resolved by many cells at low resolution. To capture the physics at smallest scales, one needs to refine the resolution in the smallest eddies, capturing both velocity and magnetic field gradients. We define the characteristic sizes of the eddies as The refinement routine is called if the size of any of the two quantities is less than a threshold value: l v,B < α∆x at the point. Coefficient α is chosen to be such that the threshold scale is larger than the numerical resistive scale. In the simulations we use α = 8, which is larger than the numerical resistive scale in simulations 2D[a] and 2D [c], and ∆x is the grid spacing at a given resolution. Coarsening of the grid in the numerical domain is only allowed if both quantities at a given grid point are larger than the threshold. Since the electric current density is roughly given by the gradient of the magnetic field, j ∼ ∇ × B, regions of the large electric current density (indicating current sheets) are automatically refined. Since the inverse cascade is very pronounced in 2D simulations, AMR shows very high efficiency at early times, when the spectrum is being formed, and at later times, when small eddies merge in larger ones (see Figure 5a for η = 10 −6 ). Since the resistive scale is much larger for η = 10 −5 , the coverage by the highest resolution level does not exceed 15% in this case.
The threshold value is tested for a resolution of 32768 2 grid points by comparing spectra of the magnetic field energy for uniform grid and AMR-enabled runs (where the effective resolution for the AMR runs indicates the total resolution if the whole domain were refined to the highest AMR level) at the same moment in time (see Figure 5b). Interestingly, the most accurate spectra are produced by simulations where the refinement algorithm is called only every 50 − 100 time-steps, most likely due to less numerical noise being generated during the refining and coarsening of the grid and re-interpolation. The frequency of the refinement calls is defined to ensure that the finest structures are always located inside the refined grid block during their motion in the domain. For the bulk velocity of the fluid u ≈ 0.1c, and CFL number 0.4, an element of the fluid travels about 20 cells between two calls of the refinement, while the minimum size of a refined grid is 32 2 cells.
This AMR strategy does not seem to be effective in 3D simulations due to the overall low grid resolution, compared to the extreme resolutions employed in 2D: AMR automatically chooses the resolution needed to resolve all the features in the block. Since the size of all features in the flow is rather defined by the numerical resolution than by an explicit resistivity, AMR refines the whole domain up to the highest available resolution. It is impossible to find a reasonable threshold α for the 2D counterpart of the highest resolution 3D run with 3200 3 grid points: any chosen α either truncates the inertial range of the spectrum or refines the whole domain shortly after the start of the simulation. Figure 5c demonstrates that a base resolution of 32000 2 grid points is needed to fully resolve the resistive scale for η = 10 −6 and keep the inertial range of the turbulence unaffected by the resolution. In order to demonstrate this, we compare spectra for resolutions with 16384 2 , 32768 2 , 65536 2 points.

B. KELVIN-HELMHOLTZ STABILITY OF CURRENT SHEETS
In order to study the stability of the magnetized shear flow in our plasmoid-unstable current sheets in 2D, we calculate the value of the in-plane reconnecting magnetic field component B || and the out-of-plane component B z as well as velocity field components v || , v z for each of the three slices shown in Figure 6a by green lines. For all these slices we find |B ⊥ | |B || | and |v ⊥ | |v || |, where || represents the direction parallel to the current sheet at a given point in the slice (the arrows in Figure 6a indicate parallel and perpendicular directions, and the z-direction is out-of-plane). We show the typical behavior of these parameters in Figure 6b, which implies that the flow satisfies the non-relativistic stability criterion |δv| < |δB|.
For each slice across the current sheet, we run a local simulation of the shear flow, with one flow having parameters (ρ, B || , B z , v || , v z ) given by the upstream of the current sheet, and its counter-flow having parameters from the interior of the current sheet, particularly, B || = 0. Zero parallel magnetic field in one of the two interacting flows prohibits reconnection at their interface in these experiments, but allows to study the KHI. Plasma pressure is adjusted to maintain the force balance across the interface of the flows.
We run simulations with a resolution of 4096 2 grid points for 20 light-crossing times along the sheet, which exceeds the life time of intermittent current sheets in the full 2D turbulence simulation. For the intermediate slice (shown by the solid line in Figure 6a), we run an AMR-enabled simulation with an effective resolution of 32768 2 grid points, which resolves all the scales up to the resistive scale for a resistivity η = 10 −6 (see Appendix A for the resolution study). Figure 7. Intermittency of the turbulence: The two top rows show the conditional probability distribution function (PDF) of the dynamic alignment angle for a fixed amplitude of the Elsasser field δz + /δz + l and at a fixed scale l (the top row, l = 0.1, corresponds to the low values of k ⊥ , approximately at the beginning of the intertial range; the middle row, l = 0.01,-corresponds to the high-k ⊥ end of the inertial range), measured at t = 2. Here, δz + l = exp ln δz + |l is the geometrical mean of δz + at a fixed scale l. The bottom row shows the PDF of δz + /δz + l for 30 logarithmically distributed scales between l = 0.9L (red line) to l = 0.01L (blue line).
This finest grid covers the whole interface of the flows at any moment of the simulation. In all of these experiments we do not observe any instability growth. This implies that the in-plane magnetic field in the upstream of the current sheet is capable of preserving KH-stability in both the upstream and downstream of the current sheet.
In order to explore the KH-stability of the 3D current sheet, we select slices between the two plasmoids (red-blue line in Figure 4), and perform test simulations in a 3D setup with a geometry similar to the one described above for 2D simulations. This slice is also normal to the surface of the current sheet at the point of their intersection. We test the sheet's stability by running a 3D simulation with a resolution of 1024 3 grid points for 20 light-crossing times along the shear interface, with parameters corresponding to the slice shown in Figure 4c. Since the lack of scale invariance is the most prominent sign of intermittency, we focus on the scale dependence of the Elsasser field increments, δz + . The results presented in this Appendix are similar for δz − as we expect the turbulence to be balanced, δz + ∼ δz − . As shown by Mallet et al. (2015), the scale invariance can be characterized by the similarity of the conditional probability distribution functions (PDFs) P(δz + |l) of δz + at different scales l. We computed these PDFs for 2D (run 2D[a]) and 3D (runs 3D[a] and 3D[d]) simulations. To measure the PDF, we use a set of 30 logarithmically spaced scales {l i } from 0.01L to 0.9L. The smallest scale l 1 , corresponding to k ⊥ ≈ 100, lies deeply in the inertial range of the energy spectrum, while the largest scale of the considered set, l 30 , is located at the outer scale of the turbulence. As the bottom row of Figure 7 shows, the smaller-scale eddies (darker lines) have higher probability to reach large normalized amplitudes of δz + . The flattening of the PDFs at smaller scales can be attributed to the sheet-like structures emerging at these scales. For 2D and weak guide field 3D simulations, the presence of long current sheets can also explain a flatter tail of the PDF at the larger scales, while the PDF for the strong guide field 3D case has an abrupt cutoff at the high δz + for the same eddy sizes. We normalized δz + (l) by a geometrical mean of δz + = exp ln δz + |l at a given scale l, as it is less sensitive to outliers than an arithmetical mean.
The intermittent, scale-dependent, nature of the dynamic alignment can also be shown by measuring the PDF of the dynamic alignment angle at given scales, as considered by Dong et al. (2018). We are, however, also interested in testing the assumption of Chandran et al. (2015) that large δz + rotates δz − into alignment, while balanced collisions δz + ∼ δz − ∼ δz ± are not aligned. This anti-correlation of the alignment angle with the amplitude of δz + contradicts the intuitive explanation of the dynamic alignment by an uncertainty principle. To test this, we measure the conditional PDF of the dynamic alignment angle P(θ|l, δz + ) for a given scale l and the amplitude of the Elsasser field δz + /δz + . The middle and top rows of Figure 7 show that the the prediction is matched perfectly for strong guide field 3D turbulence: the larger δz + , the more aligned δz + and δz − are. For 2D and weak guide field 3D turbulence there is a deviation from this prediction at the outer scale: while the statement holds for intermediate amplitudes of δz + , at high amplitudes eddies become uncorrelated again. The most powerful increments δz + are associated to current sheets and plasmoids, and one can expect that circular plasmoids have an alignment angle (the ratio of two length scales of the eddy) θ ∼ 1 that can explain decorrelation of the alignment angle at high δz + /δz + .