On the energy spectrum of strong magnetohydrodynamic turbulence

The energy spectrum of magnetohydrodynamic turbulence attracts interest due to its fundamental importance and its relevance for interpreting astrophysical data. Here we present measurements of the energy spectra from a series of high-resolution direct numerical simulations of MHD turbulence with a strong guide field and for increasing Reynolds number. The presented simulations, with numerical resolutions up to 2048^3 mesh points and statistics accumulated over 30 to 150 eddy turnover times, constitute, to the best of our knowledge, the largest statistical sample of steady state MHD turbulence to date. We study both the balanced case, where the energies associated with Alfv\'en modes propagating in opposite directions along the guide field, E^+ and $E^-, are equal, and the imbalanced case where the energies are different. In the balanced case, we find that the energy spectrum converges to a power law with exponent -3/2 as the Reynolds number is increased, consistent with phenomenological models that include scale-dependent dynamic alignment. For the imbalanced case, with E^+>E^-, the simulations show that E^- ~ k_{\perp}^{-3/2} for all Reynolds numbers considered, while E^+ has a slightly steeper spectrum at small Re. As the Reynolds number increases, E^+ flattens. Since both E^+ and E^- are pinned at the dissipation scale and anchored at the driving scales, we postulate that at sufficiently high Re the spectra will become parallel in the inertial range and scale as E^+ ~ E^- ~ k_{\perp}^{-3/2}. Questions regarding the universality of the spectrum and the value of the"Kolmogorov constant"are discussed.


I. INTRODUCTION
Astrophysical plasmas are typically magnetized and turbulent, with turbulent fluctuations spanning a tremendous range of scales in which the energy spectrum follows a power law scaling [e.g., 1,2]. Incompressible magnetohydrodynamics (MHD) provides the simplest theoretical framework for studying magnetized plasma turbulence. The precise form of the MHD turbulence spectrum is crucial for a variety of processes in astrophysical systems with extended inertial intervals, such as plasma heating and wave-particle interactions, which are sensitive to small variations in the spatial scaling of the fluctuations [e.g., [3][4][5]. The incompressible MHD equations take the form where z ± = v ± b are the Elsässer variables, v is the fluctuating plasma velocity, b is the fluctuating magnetic field (in units of the Alfvén velocity), V A = B 0 / √ 4πρ 0 is the Alfvén velocity based upon the uniform background magnetic field B 0 , P = (p/ρ 0 + b 2 /2), p is the plasma pressure, ρ 0 is the background plasma density, ν is the fluid viscosity (which, for simplicity, we have taken to be equal to the magnetic diffusivity) and f ± represent forces that drive the turbulence at large scales. It can be shown that in the limit of small amplitude fluctuations, and in the absence of forcing and dissipation, the system describes non-interacting linear Alfvén waves with dispersion relation ω ± (k) = ±k V A . The incompressibility condition requires that these waves be transverse. Typically they are decomposed into shear Alfvén waves (with polarizations perpendicular to both B 0 and to the wavevector k) and pseudo-Alfvén waves (with polarizations in the plane of B 0 and k and perpendicular to k).
Nonlinear interactions (or collisions) between counterpropagating Alfvén wave packets distort the packets, splitting them into smaller ones until a scale is reached when their energy is converted into heat by dissipation [6]. The efficiency of the nonlinear interaction is controlled by the relative size of the linear and nonlinear terms in equation (1): The regime in which the linear terms dominate over the nonlinear terms is known as weak MHD turbulence, otherwise the turbulence is called strong. The Fourier energy spectrum of MHD turbulence can be derived analytically only in the limit of weak turbulence [e.g., [7][8][9][10][11][12][13][14]. However, it has been demonstrated both analytically and numerically that the energy cascade occurs predominantly in the plane perpendicular to the guiding magnetic field [8,10,15,16], which ensures that even if the turbulence is weak at large scales it encounters the strong regime as the cascade proceeds to smaller scales. Although weak turbulence may exist in some astrophysical systems [e.g., 13,[17][18][19], magnetic turbulence in nature is typically strong, for which an exact analytic treatment is not available. In this case, highresolution, well-optimized numerical simulations play a significant role in guiding our understanding of the turbulent dynamics. This provides the motivation for the present work. The ideal MHD system conserves the Elsässer energies E + = 1 4 (z + ) 2 d 3 x and E − = 1 4 (z − ) 2 d 3 x (equivalently, the total energy E = E v + E b = 1 2 (v 2 + b 2 )d 3 x = E + + E − and the cross-helicity H C = (v · b)d 3 x = E + − E− are conserved). The energies E + and E − cascade in a turbulent state toward small scales due to the nonlinear interactions of oppositely moving z + and z − Alfvén packets. MHD turbulence is called balanced when the energies carried by oppositely moving fluctuations E ± are equal, and it is called imbalanced when they are not the same. MHD turbulence in nature and in the laboratory is typically imbalanced. For instance, this is the case when the turbulence is generated by spatially localized sources, as is the case in the solar wind where more Alfvén waves propagate away from the Sun than towards it. The independent conservation of the two Elsässer energies (compared to only one conserved energy in hydrodynamics) has a profound consequence for the MHD dynamics [e.g., [20][21][22][23][24][25][26][27][28][29].
In this work we present the results of a series of direct numerical simulations of MHD and Reduced MHD for balanced and moderately imbalanced turbulence and investigate how the scalings of the Elsässer spectra behave as the Reynolds number is increased. We also present the first high-resolution direct comparison of simulations of MHD vs RMHD turbulence, demonstrating that the latter model completely captures the turbulence dynamics of strong MHD turbulence at roughly half the computational cost of a full MHD simulation. This paper is organized as follows. In section II we briefly describe the most recent phenomenological efforts to understand scaling laws in MHD turbulence, particularly in the imbalanced case. In section III we describe the numerical set up and the parameter regime for our simulations. In section IV we show measurements of the energy spectrum from a series of numerical simulations with varying Reynolds numbers. In section V we show measurements of scale-dependent dynamic alignment and establish its relation with the −3/2 scaling of the energy spectrum. In section VI we discuss the approach to the universal regime and the universality of Kolmogorov's constant in MHD. We show that dynamic alignment introduces a new robust scale-dependent quantity that enters the definition of the energy spectrum and uniquely sets the Kolmogorov constant. We propose that this new quantity is a consequence of cross-helicity conservation. Finally, in section VII we discuss our results.

II. MHD TURBULENCE PHENOMENOLOGY
For strong MHD turbulence, Goldreich & Sridhar [30] argued that the pseudo-Alfvén modes are dynamically irrelevant for the turbulent cascade (since strong MHD turbulence is dominated by fluctuations with k ⊥ ≫ k , the polarization of the pseudo-Alfvén fluctuations is almost parallel to the guide field and they are therefore coupled only to field-parallel gradients, which are small since k ≪ k ⊥ ). If one filters out the pseudo-Alfvén modes by setting z ± = 0, it can be shown that the re-sulting system is equivalent to the Reduced MHD model: We note that in RMHD the fluctuating fields have only two vector components, but that each depends on all three spatial coordinates. Moreover, because the z ± are assumed incompressible (∇ · z ± = 0), each field has only one degree of freedom which is more commonly expressed in terms of stream functions in the more standard form of the RMHD equations [31,32]. Conservation of both the Elsässer energies means that once an imbalance has been created it cannot be destroyed by the MHD dynamics. It is also well known that decaying MHD turbulence, affected only by the dissipation, becomes increasingly more imbalanced with time [e.g., 20,26,27]. Several analytic and numerical studies have shown that imbalance is also an inherent property of driven MHD turbulence even if the turbulence is forced without introducing a net imbalance at the largest scales -the turbulent domain spontaneously fragments into local imbalanced domains where the cross helicity is either positive or negative [21-26, 28, 29].
In imbalanced domains, the directions of the magnetic and velocity fluctuations are not independent, rather, they are either aligned or counter-aligned to a certain degree [58]. The organization of such a domain is the following: the directions of both the magnetic and velocity fluctuations vary within a small angle (comparable to the alignment angle) throughout the domain, while their amplitudes change predominantly in the direction normal to their polarizations. Such positively and negatively aligned domains appear to be the building blocks of MHD turbulence, whether it is balanced overall or not.
The origin of such domains can be qualitatively understood from the conservation of energy and cross-helicity in an ideal MHD system. When small dissipation is present and the system is unforced, it can be argued that energy decays faster than cross-helicity. This selective decay would eventually lead to Alfvénization of the flow, that is, to progressively stronger alignment (or counteralignment, depending on the initial state) between the directions of the magnetic and velocity fluctuations, e.g., [20,21,33,34]. In a perfectly aligned (counter-aligned) state either z + or z − is identically zero, and the nonlinear interaction vanishes. In a driven state, characterized by strong nonlinear interaction and a constant energy flux over scales, the alignment cannot be perfect. Rather, it turns out that alignment depends on the scale, the smaller the scales the better the alignment. Below we will demonstrate this phenomenon in numerical simulations. From a more qualitative point of view, one can argue that whenever a partly aligned domain appears, nonlinear interaction inside such a domain gets reduced, and its evolution time increases compared to non-aligned domains. Therefore aligned domains persist longer, which explains the tendency of a turbulent flow to exhibit such self-organization. These aligned domains are the domains where the essential energy of the turbulence is contained, and they are typically well seen in numerical simulations. Solar wind observations also show that globally balanced turbulence is made up of locally imbalanced patches at all scales [35][36][37].
In the aligned or imbalanced domains, the Elsässer energies are unequal, and one can ask whether their spectra have to be the same. This raises questions of whether MHD turbulence is universal and scale-invariant. Indeed, if imbalanced domains have different spectra that depend on the degree of imbalance, their superposition may not have a universal scaling.
Phenomenological treatment of strong imbalanced MHD turbulence is complicated by the fact that one can formally construct two time scales for the nonlinear energy transfer: The times of nonlinear deformation of the z ± packets at some spatial scale λ are τ ± ∼ λ/z ∓ λ , which can be significantly different in the case of strong imbalance, [e.g., 20,27]. In recent years, several phenomenological models attempting to accommodate this difference have been proposed. However, the theories have generated conflicting predictions because they use different assumptions regarding the physics of the nonlinear energy cascade. For example, the theory by Lithwick et al. [38] concludes that in the imbalanced regions the Elsässer spectra have the scalings ; the same spectra were also suggested by Beresnyak and Lazarian [39]. The theory by Chandran [40] proposes that the spectra of E + (k ⊥ ) and E − (k ⊥ ) are different depending of the degree of imbalance, while the theories by Perez and Boldyrev [28] and Podesta and Bhattacharjee [41] find that the spectra of E + (k ⊥ ) and E − (k ⊥ ) have different amplitudes but the same scalings One would expect that numerical simulations could clarify the picture. However, the first numerical simulations of strongly imbalanced MHD turbulence [e.g., 28,39,42] also produced conflicting results regarding which power law E ± should follow. The conflicting numerical findings apparently reflect the fact that imbalanced MHD simulations require significantly more computational effort compared to the balanced cases [43]. This happens since in the imbalanced domains the nonlinear interaction is depleted and the Reynolds and magnetic Reynolds numbers are reduced. This can be formally seen from the fact that, in a strongly imbalanced domain with z + ≫ z − , the z + field is advected by a lowamplitude z − field, and therefore z + becomes directly affected by the dissipation at smaller wave-vectors (compared with the balanced case), which reduces its inertial interval. Now, z − is advected by a strong z + , but z + is significantly affected by the dissipation, so the inertial interval of z − becomes spoiled as well.
In order to produce large inertial intervals simultaneously for both Elsässer fields when strongly imbalanced domains are present in the flow, one therefore needs to have a significantly higher Reynolds number as compared to the balanced case. However, as one increases the Reynolds number, one needs to increase the numerical resolution in order to appropriately resolve the small scales and to make sure the numerical run is stable. Therefore, the larger the imbalance, the larger the numerical resolution required to describe correctly the Elsässer spectra. Fortunately, it has been argued that Reduced MHD can be used to investigate the universal properties of MHD turbulence, which offers the advantage that an RMHD simulation can be achieved at half the cost of an MHD simulation.

III. NUMERICAL SETUP
We solve the MHD equations (1) and their RMHD counterpart (2) in a periodic, rectangular domain with aspect ratio L 2 ⊥ × L , where the subscripts denote the directions perpendicular and parallel to B 0 , respectively. We set L ⊥ = 2π, L /L ⊥ = 6 or 10 and B 0 = 5e z . A fully dealiased 3D pseudo-spectral algorithm is used to perform the spatial discretization on a grid with a resolution of N 2 ⊥ × N mesh points. We note that the domain is elongated in the direction of the guide field in order to accommodate the elongated wave-packets and to enable us to drive the turbulence in the strong regime while maintaining an inertial range that is as extended as possible (see [44]). This is a physical requirement that should be satisfied no matter what model system, full MHD or reduced MHD, is used for simulations.
In the case of reduced MHD though, when the z ± components are explicitly removed, the resulting system (2) is invariant with respect to simultaneous rescaling of the background field B 0 and the field-parallel spatial dimension of the system, if one neglects the dissipation terms. Therefore, for any strength of the background field B 0 ≫ 1, one can rescale the field to B 0 = 1 and the field-parallel box size to L = L ⊥ , that is, conduct the simulations in a cubic box. We should note however that the dissipation terms in (2) are not invariant and they should be changed accordingly under such rescaling.
To save on computational cost we have reduced the field-parallel numerical resolution for some simulations, i.e., the numerical grid is anisotropic with L /N > L ⊥ /N ⊥ . This is appropriate since the energy cascade proceeds much faster in the field-perpendicular direction, and the energy spectra decline relatively slowly in the field-perpendicular direction and relatively fast in the field-parallel direction. Energies at large k are therefore reduced and a lower field-parallel resolution is not expected to alter the behavior of the spectra in the inertial interval. An isotropic resolution with the value imposed by the field-perpendicular dynamics would therefore be wasteful.
We should however caution that a reduced resolution (or, equivalently, unreasonably high Reynolds number for a given resolution) may contaminate the dissipative physics, even if the inertial interval is unaffected. For example, if the precise scaling behavior in the dissipation interval is of interest, as is the case for extended scaling laws such as the dynamic alignment angle, somewhat smaller Reynolds numbers may need to be chosen. As a general rule, whether the numerical simulations are conducted to investigate the inertial or the dissipation interval, a resolution study must be performed in order to establish the optimal Reynolds number for a given task. In particular, it has to be verified that increasing the numerical resolution while keeping the physical parameters such as Reynolds number, forcing mechanism, etc. unchanged does not affect the studied spectra, e.g., [45]. This point will be illustrated below in the balanced case.
The turbulence is driven at the largest scales by colliding Alfvén modes [59]. We drive both Elsässer populations by applying statistically independent random forces f + and f − in Fourier space at wave-numbers 2π/L ⊥ ≤ k ⊥ ≤ 2(2π/L ⊥ ), k = 2π/L . The forces have no component along z and are solenoidal in the xy-plane. All of the Fourier coefficients outside the above range of wavenumbers are zero and inside that range are Gaussian random numbers with amplitudes chosen so that v rms ∼ 1. The individual random values are refreshed independently on average approximately 10 times per turnover of the large-scale eddies. The variances σ 2 ± = |f ± | 2 control the average rates of energy injection into the z + and z − fields. We take σ + > σ − and in the statistically steady state we measure the degree of imbalance through the pa- Thus h = 0 corresponds to balanced turbulence and h = 1 defines maximally imbalanced turbulence. Time is normalized to the large scale eddy turnover time τ 0 = L ⊥ /(2πv rms ). The field-perpendicular Reynolds number is defined as Re ⊥ = v rms (L ⊥ /2π)/ν ≈ 1/ν. In order to accommodate the reduced field-parallel resolution we have also modified the diffusion operator in equations (1) and (2), i.e., we have replaced ν∇ 2 with ν(∂ xx + ∂ yy ) + ν ∂ zz .
The system is evolved until a stationary state is reached, which is confirmed by observing the time evolution of the total energy of the fluctuations. The data are then sampled in intervals of the order of the eddy turnover time. All results presented correspond to averages over 30-150 samples for each run. As shown in Table III, we conduct a number of MHD and RMHD simulations in the balanced and imbalanced regime in order to investigate the scaling of the energy spectra as the field-perpendicular Reynolds number increases.

IV. MEASUREMENTS OF THE ENERGY SPECTRUM
The field-perpendicular energy spectrum is obtained by averaging the angle-integrated Fourier spectrum,  over field-perpendicular planes in all samples. Identifying the inertial range in numerical simulations with limited resolution is generally difficult, due to the relatively modest separation between the forcing and dissipation scales that current super-computers can afford. For instance, a measurement of the turbulence spectrum for a single Reynolds number is not enough to ensure that the simulated turbulence has converged to the asymptotic universal scaling. Instead, one carries out a set of numerical simulations with increasing resolution and Reynolds number. The spectra are then compensated by the different phenomenological predictions and the preferred model is distinguished by the best fit. In Figures   1 to 4 the inertial range is identified by the flat regions of the spectra compensated by k 3/2 , which extend further to the right with increasing Reynolds number (and resolution). Figures 1 and 2 show the total field-perpendicular energy spectrum E(k ⊥ ) in the balanced regime for the RMHD and MHD cases, respectively. The RMHD and MHD spectra are remarkably similar, confirming that the pseudo-Alfvén modes are dynamically insignificant and that the RMHD approximation is valid. In both cases the total energy spectrum remains of the form as the Reynolds number increases, with the inertial range starting at k ≈ 4 and extending up to k 30 in the highest Reynolds number case. In neither RMHD or MHD is there any evidence of a build up of energy close to the dissipative wave-numbers-often referred to as a bottleneck effect-with both spectra falling off smoothly in the dissipative range. Figures 3 and 4 show the field-perpendicular Elsässer spectra in the imbalanced regime for the RMHD and MHD cases, respectively. Again the behavior of both spectra in the RMHD and MHD regimes are very similar. In both cases, it is seen that while E − keeps the scaling as the Reynolds number increases, the scaling of E + (k ⊥ ) is more difficult to pin down. Indeed, both the RMHD and MHD results for Re = 2200 yield a steeper spectrum for E + (k ⊥ ), with an exponent possibly nearer to −5/3 than −3/2. However, we believe that there is no real significance to the value of −5/3 here, the exponent is simply steeper than −3/2. Indeed, in both cases, as the Reynolds number is increased E + (k ⊥ ) appears to flatten, which means that E + (k ⊥ ) has not fully established the universal scaling behavior yet. Since E + (k ⊥ ) and E − (k ⊥ ) are pinned (i.e., converge to each other) at the dissipation scales and are anchored (i.e., independent of the Reynolds number) at the driving scales, we postulate that at sufficiently high Re (where the inertial range is extensive) the spectra will become parallel in the inertial range and attain the scal- . Numerical tests of this prediction must await a significant increase in computational power.

V. MEASUREMENTS OF DYNAMIC ALIGNMENT
An important test that can be performed in the presented simulations concerns the so-called dynamic alignment angle. This angle is defined by the following ratio of the two specially constructed structure functions [25]: where δv ⊥ (l) and δb ⊥ (l) are the field-perpendicular velocity and magnetic field increments, respectively, corresponding to the field-perpendicular scale separation l.
(We note that in definition (4) we have assumed that the angle is small, and hence no distinction between θ(l) and sin θ(l) is made. Hereafter, by θ(l) we will always understand the quantity (4)).
As proposed in [23,24] the alignment angle θ(l) has a nontrivial scaling with l, which may explain the observed −3/2 scaling exponent of the energy spectrum. As discovered in [45], the scale dependent dynamic alignment exists not only in the inertial interval, but it also extends into the dissipation range and is limited only by the grid size of the numerical scheme. We will demonstrate that the alignment angle scaling provides a sensitive test probing the turbulent cascade deep in the dissipation interval.
In particular we will see that if the simulated dissipation range is under-resolved (e.g., as a result of the use of too large a Reynolds number or strongly anisotropic resolution), the dynamic alignment can be easily spoiled at the dissipation scales even if it is present in the inertial interval.
The measurements of the alignment angle are presented in Figure 5. The first panel shows three simulations (RB2d,c,b in Table III) performed at the same numerical resolution of 1024 3 but with different Reynolds numbers Re = 1800, 3200, 6000. Plots for Re = 1800, 3200 show a remarkable property of the alignment scaling: It extends deep down into the dissipation region, practically up to the scale of the numerical discretization, independently of the Reynolds number (see also [45]). However, this behavior is spoiled if the Reynolds number is pushed to very high values, at which the dissipation interval becomes under-resolved. In this case, the scaling starts to degrade at large wave-numbers, as is seen in the case Re = 6000.
The alignment scaling is however restored back to its original value if the numerical resolution is increased to 2048 3 , so that the dissipation scales become well resolved again. This is seen from comparison of the plot for RB2b (1024 3 , Re = 6000) in the first panel of Figure 5 with the plot for RB3d (2048 3 , Re = 5700) in the second panel. Further increase of the Reynolds number in the second panel of this figure demonstrates that the alignment scaling is stable up to Re = 9000 (RB3c, 2048 3 ), however, it starts to degrade at large wave-numbers for higher Reynolds numbers Re = 15000 (RB3b, 2048 3 ), in complete analogy with the behavior depicted in the upper panel of Fig. 5 at smaller resolution. The alignment angle is spoiled even more in the run RB3a (2048 2 × 512, Re = 15000) in the same figure where we simultaneously decrease the field-parallel numerical resolution, making the dissipation interval even more under-resolved. Note however, that in both the first and the second panels of Fig. 5, the heaviest distortion of the alignment behavior occurs in the dissipation region, while the inertial interval (approximately contained between the two vertical lines) is relatively unaffected. This may explain why an under-resolved dissipation interval is not manifest in the scaling of energy spectra, as seen in Figure 6. Fig. 7 shows three well resolved simulations with numerical resolutions increasing from 512 3 to 1024 3 to The dynamic alignment scaling extends well into the dissipation range, up to scales close to the grid cell (roughly l ∼ 3 grid cells). When the Reynolds number is pushed to very high values (so that the dissipation interval becomes under-resolved) or the numerical resolution in the field-parallel direction is reduced, the alignment-angle scaling degrades at small scales. The vertical lines show the approximate boundaries of the inertial interval (cf. Fig. (8)). The straight dotted line has a slope of 1/4. 2048 3 . We observe that the scaling interval of the alignment angle becomes progressively longer and its scaling index stays close to the predicted value 1/4 [23] with little or practically no dependence on the Reynolds number [60]. This means that we observe a truly universal scaling behavior of the dynamic alignment. The lower panel of Fig. 7 shows the same curves where the spatial scale is normalized by the dissipation length. We observe that the flattened parts of the curves at small scales do FIG. 6: Energy spectra for runs RB3a (solid) and RB3d (dash). Simulation RB3a on 2048 2 × 512 has an unresolved dissipation at the expense of longer inertial interval. Simulation RB3d is performed at lower Re to capture alignment in the dissipation region, with a shorter inertial range.
not overlap under such rescaling, which supports our observation mentioned above that the extent of the scaling interval is not defined solely by the dissipation scale, but rather depends on the numerical discretization step.

VI. ENERGY SPECTRUM: KOLMOGOROV CONSTANT AND DISSIPATION SCALE
For a more complete study of the energy spectrum, one can also evaluate the amplitude of the spectrum and the dissipation scale for each simulation and verify that they agree with a given phenomenology. Since our spectral scaling conforms to the phenomenology of Boldyrev [23,24], we now study in more detail the scaling associated with this model. First, we need to derive the expression for the energy spectrum, which is done in the following way [23]. The time of nonlinear interaction at field-perpendicular scale λ in this model is τ ∼ λ/(v λ θ λ ), where v λ denotes the typical (rms) velocity fluctuations, θ λ = θ 0 (λ/L ⊥ ) 1/4 is the scale-dependent alignment angle between magnetic and velocity fluctuations, which was studied in the previous section, and θ 0 is the typical alignment angle at the outer scale (forcing scale) L ⊥ . The rate of energy cascade is then evaluated as One however notices that the amplitude of the energy spectrum is not uniquely defined in this equation, since the outer-scale quantities θ 0 and L ⊥ essentially depend on the forcing routine. This is understood from the following example. Assume that the large-scale force drives only unidirectional Alfvén waves z + , for which v is perfectly aligned with b and θ 0 = 0. Then the wave energy will grow without bound, since the nonlinear interaction leading to the energy cascade and eventual dissipation at small scales is absent.
Even when a particular forcing routine is specified, the definitions of the values of θ 0 and L ⊥ are still subjective  (4) in balanced RMHD. The frames show numerical simulations with increasing Reynolds number and numerical resolution with properly resolved dissipation ranges (runs RB1c (dashdotted), RB2c (dashed), and RB3d (solid)). In the lower plot, the scale l is rescaled by the dissipation length (see section VI for precise definitions). It can be seen from here that the region of scale dependent dynamic alignment increases as smaller scales are made available by increased numerical resolution. The extent of the alignment region is not limited by the dissipation scale, but rather depends on the grid size of the numerical scheme. The straight dotted line has the slope 1/4. since they essentially rely on the outer-scale properties of turbulence rather than on the measurements of the inertial interval. We now propose that this problem can be remedied in an efficient way. For that we notice that there exists a well-defined quantity that is remarkably stable (scale-independent) in the inertial interval: where θ(l) is defined in (4), see the discussion in the preceding section. In this definition one can use any scale l from the inertial interval or dissipation interval if the numerical simulations are well resolved. A somewhat simpler rule can be used in numerical (or observational) studies, where one does not have to know a priori what scales correspond to the inertial interval and does not have the luxury of having the plot in Fig. 7 available. In this case l in formula (5) can be chosen to be the Taylor micro-scale based on either the magnetic or the velocity fluctuations, l = v rms /|∇×v| rms or l = b rms /|∇×b| rms , assuming the magnetic Prandtl number is of order one. We therefore propose the following normalization of the energy spectrum: where Λ is defined by (5). The scale Λ that is defined solely through the inertial-interval quantities, incorporates the essential information about the cross-helical structure of MHD turbulence. It is not uniquely defined by the outher scale of the turbulence, rather it also depends on the large-scale driving mechanism. Therefore, the inertial-interval energy spectrum is defined by the two quantities ǫ and Λ, characterizing the energy cascade rate and the level of cross-helical organization of the flow. The presence of the two quantities characterizing the spectrum of MHD turbulence (as oppose to only one quantity in hydrodynamic turbulence) is the manifestation of the two conserved quantities cascading toward small scales in MHD turbulence: energy and cross-helicity. We expect that the constant C k in (6) may be "universal," that is, largely independent of the character of the driving, analogous to the Kolmogorov constant in hydrodynamical turbulence. This constant can be measured in our simulations in the following way. First, we specify l that we use to measure the alignment scale Λ in (5). According to our plots in Figs. 5 and 7, we may choose l = 0.07L ⊥ , say, as a scale belonging to the inertial interval and not yet affected by the numerical resolution effects. Then, for simulations RB1a, RB2a, RB3a we find Λ = 1.34L ⊥ , 1.41L ⊥ , 1.48L ⊥ , respectively.
The dissipation rate can be evaluated based on the energy spectrum (6) as follows: Our numerical results confirm that the integral of ν k 2 leads to a negligible correction to the dissipation rate, and therefore it can be omitted, and we can use the field-perpendicular spectrum E(k ⊥ ) = E(k , k ⊥ )dk . Then, for simulations RB1a, RB2a, RB3a we find: ǫ = 0.15, 0.15, 0.16. The dissipation scale can be found (or defined) based on the energy spectrum. Omitting the dimensionless constants, we then accept, by definition, We can demonstrate that our simulations agree with this scaling by plotting the energy spectra in the balanced case (RB1a,2a,3a) versus the wave-vector normalized with the dissipation scale (8), where we measure  [61]. Note that with the wave vector normalized with the single parameter η, the whole spectra collapse onto each other, thus providing additional evidence that the universal functional behavior of the spectrum is obtained in our simulations. The lower plot in figure 8 shows that the length of the inertial range increases as l 0 /l d ∼ Re 2/3 , also in good agreement with the estimate for the dissipation scale (8). The "Kolmogorov constant" C k can be evaluated from the upper plot as C k ≈ 2.

VII. DISCUSSION
We have presented results from state-of-the-art direct numerical simulations of balanced and imbalanced driven MHD turbulence. The simulations are achieved at the extremely large numerical resolution of 2048 3 and the longest running time, with many runs span-ning more than a hundred eddy turnover times in the steady state. The simulations were performed using two pseudo-spectral codes, one solving the MHD equations and the other solving the RMHD equations. In theories and simulations of MHD turbulence, it has long been argued that RMHD provides a correct and accurate framework for investigating the universal properties of MHD turbulence both in the weak and strong turbulence regimes. We have presented a direct comparison of high-resolution numerical simulations of MHD vs RMHD turbulence using two independently developed pseudospectral codes with identical parameters. It is shown that in the strong turbulence regime, in both the balanced and imbalanced state, the energy spectrum of the Elsässer variables in MHD and RMHD are in remarkable agreement (for details of a lower resolution comparison, including the individual velocity and magnetic spectra and the alignment angle, see [46]). These results are of essential value for MHD turbulence research, as simulating MHD turbulence can be accomplished using RMHD codes that generally incur a smaller computational cost.
In the balanced case, the simulated energy spectra of E + and E − show a clearly identifiable inertial range, consistent with a slope of k −3/2 ⊥ for both E + and E − . It is observed from Figures 1 and 2 that the compensated energy spectra show a flat region that extends as the Reynolds number is increased. This is consistent with previous, lower resolution simulations of strongly magnetized MHD turbulence, e.g., [16,[47][48][49][50][51][52]. In the imbalanced case, the interpretation of the numerical results is not as straightforward. Figures 3 and 4 show that the energy spectrum of E − remains reasonably close to k −3/2 ⊥ , only slightly changing its overall amplitude for small Reynolds numbers. As for the E + spectrum, the compensated spectrum shows a slope slightly steeper than −3/2 which however flattens as the Reynolds number increases. Another observation from the large Reynolds number imbalanced numerical simulations is that the spectra of E + and E − are "anchored" at large scales and "pinned" at the dissipation scale. From these results we propose that the energy spectra of E + becomes asymptotically closer to k −3/2 ⊥ as the Reynolds number is increased. Much higher resolutions, exceeding the capabilities of today's supercomputers, are required to conclusively demonstrate this conjecture.
Finally, during the refereeing process, our attention was drawn to recent publications by the group of Beresnyak & Lazarian [53][54][55], in which the authors address issues similar to the ones contained in this paper. Most of the conclusions of those papers appear to be at odds with ours (and with similar results or other groups, e.g., [16,[47][48][49][50][51]). We however note that the actual numerical results presented in [53][54][55] agree with ours in the range of scales that we study, while they differ from ours at very large wavenumbers, e.g., k 50 in the runs with highest resolution. Beresnyak & Lazarian suggest that the true inertial interval exists only at these large wavenumbers where they perform their measurements of the scaling relations. The formal cause of the disagreement of our conclusions with those by Beresnyak & Lazarian is thus the numerical measurements being performed in essentially different regions of the phase space.
The question however remains as to what causes the results of the numerical simulations by Beresnyak & Lazarian to disagree with ours at small, sub-inertial scales. According to our analysis, the answer is the following: the k-space intervals on which references [53][54][55] base their conclusions are significantly affected by numerical effects due to the numerical setup they use. It is not appropriate in their simulations to use those intervals for addressing either the inertial or the dissipation regimes. We however note that the dissipation-range dynamics and the behavior of the numerical solution of the MHD equations close to the numerical cutoff is an interesting and not well studied question. It is therefore worth addressing the differences between our simulations and those by [53][54][55] in more detail. Since such analysis is not the main objective of the present work, we have presented the corresponding discussion in the Appendix. In this Appendix we comment on the numerical reconstruction of solution of the MHD equations at small scales, that is, scales within the dissipation range and close to the numerical cutoff in k-space (the dealiasing cutoff in a pseudo-spectral code). Recent publications by Beresnyak & Lazarian [53][54][55] found that the energy spectrum in this region (roughly corresponding to k 50 in their highest-resolution runs) can have a peculiar structure that is inconsistent with the structure and the scaling found in our numerical simulations. Much confusion was created by the suggestion by [54,55] that this high-k region is, in fact, the true inertial interval of MHD turbulence, while the region that is studied in our works (corresponding to 4 k 30 in highest-resolution runs) is a "non-converged" forcing-dominated region.
It is therefore useful to address the small-scale numerical solution of the MHD equations in more detail and to relate our findings and conclusions to those presented in the recent works by Beresnyak and Lazarian [53][54][55]. These references claim that the energy spectral index of MHD turbulence is −5/3 and that there is no conclusive evidence for dynamic alignment in the numerical results. In discussing what could lead to such (in our opinion, erroneous) conclusions it is useful to distinguish two factors. One is related to differences that arise because the simulations by Beresnyak & Lazarian that allegedly are identical to ours, in fact are not identical at all because of differences in the details of the numerical setup. The other is related to the methods that are used to analyze the results and, ultimately, support one claim or another. Both play a role in the origin of the disagreement.
First, we concentrate on issues that result from the different setup of the numerical simulations. In our previous publications (e.g., [43,46]) we have discussed at length those aspects of the simulation design that are essential for accurately capturing the physics of the strong turbulent cascade. It is not necessary to repeat those discussions here, however, it is important to point out that many of the simulations of Beresnyak and Lazarian [53][54][55] differ from ours through their choice of numerical hyper-dissipation, significantly smaller viscosities for a given numerical resolution, and a considerably smaller statistical ensemble from which averages are computed. Each of these factors is potentially detrimental for the observation of the correct scaling behavior. For example, the measurements of the alignment angle that are shown in Figure 3 of reference [54] and Figure 2 of reference [55] lead Beresnyak to conclude that dynamic alignment is not present in MHD turbulence as the alignment angle saturates, i.e. flattens as a function of l at small l, when the Reynolds number increases. However, those plots exhibit a behavior that is qualitatively similar to that displayed in our Figure 5, where insufficient numerical resolution is demonstrated to affect the alignment angle at small scales. It therefore reasonable to conclude that the observed flattening of the alignment angle in the simulations of references [54,55] is an artifact of unresolved dissipation scales and, possibly, part of the inertial-range scales, rather than a physical effect.
The influence of hyper-dissipation may be similarly assessed from comparing the energy spectra obtained in our work with the energy spectra obtained in, say, reference [54]. Our spectra in Figures 1, 2, 6 & 8 exhibit an extended interval with the scaling k −3/2 , identified as the inertial interval, followed by a steep decline, identified as the dissipation range. The spectra in Figure 2 of reference [54] also show an extended interval with the scaling k −3/2 (interpreted in reference [54] as a "nonconverged" range) followed at large wavenumbers by a very short steepening (interpreted in [54] as the "inertial interval") and then flattening and ultimate cut-off. In our view, such an interpretation is incorrect; the spectral behavior observed in reference [54] close to the dissipation region is not a property of the inertial interval, but rather is evidence of the so-called bottleneck effect that is expected when numerical hyper-dissipation is present. Indeed, as discussed in references [56,57], an energy spectrum abruptly terminated in k-space by hyper-dissipation or by other Galerkin-type truncation mechanisms, exhibits an inertial interval followed by a pseudo-dissipation region (steepening of the spectrum), then by a partly thermalized region (a rise in the spectrum), and then by a far dissipation range (ultimate cutoff). The measurements presented in References [54,55] are consistent with such spectral behavior, which motivates a natural explanation of their results as an inertial interval with the −3/2 scaling, modified by a substantial bottleneck effect close to the dissipation scales. Moreover, a thermalization brought about by sharp termination of the spectrum in the k-space tends to decorrelate small-scale fluctuations, which otherwise would remain strongly aligned throughout the dissipation interval, cf. our Figure 7. This is also consistent with the significant loss of dynamic alignment at small scales that is observed in references [54,55].
A more detailed comparison of our results can be made with those MHD simulations by Beresnyak [55] that employ a physical Laplacian dissipation (simulations R8 & R9 in [55]). By evaluating the Reynolds numbers for those calculations in the same way that it is done in our work, Re = v rms L/(2πν) with v rms ≈ 1, we find that simulation R8, with a resolution 768 3 mesh points, is performed at the Reynolds number Re ≈ 8000, while calculation R9 (resolution 1536 3 ) is performed at Re ≈ 20000. According to our results in Figure 5, in the simulations with a resolution of 1024 3 mesh points the dissipation interval is under-resolved already at Re ≈ 6000, while in the 2048 3 simulations the dissipation interval is underresolved at Re ≈ 15000. Thus, the runs R8 & R9 of reference [55] that are most similar to ours have lower numerical resolutions while higher Reynolds numbers. Therefore, they have essentially unresolved dissipation intervals and, possibly, parts of the inertial intervals.
The lack of resolution at the bottom of the inertial intervals in the simulations R8 & R9 can also be seen from the alignment-angle curves shown in Figure 2 of [55]. Under the rescaling applied in that figure, the curves should approach each other in the inertial interval, as they do in our Fig. 7, lower panel. In contrast, one can see only a short region of in Figure 2 of [55] (runs R8 & R9) where the curves approach each other, approximately within the range 20 l/η 40. Apparently, this is the only piece of the inertial interval that is resolved, and in this interval the scaling exponent of the angle indeed approaches 1/4, see Figure 3 in reference [55], as expected according to our results.
We now turn to the second factor that contributes to the differing conclusions drawn by the Beresnyak & Lazarian group, namely the method of analysis. We recall that the objective is to determine the scaling behaviour within the inertial range. Concerning the energy spectrum, we assess whether the numerical data preferen- by compensating the numerical data by k 3/2 and by k 5/3 in turn. For the correct model, the inertial range then corresponds to the range of scales over which the compensated spectrum is flat. We always find that k −3/2 ⊥ provides the better fit, with the inertial range starting at k ⊥ ≈ 4 and extending up to k ⊥ 30 at highest resolution. As the Reynolds number increases, numerical convergence is demonstrated by the fact that this region maintains its amplitude and scaling and increases in extent to larger wavenumbers, see, e.g. our Figure 1.
In contrast, Beresnyak [55] uses an indirect method to select the preferred spectral exponent. He uses the two phenomenological models that describe the inertial range characteristics to predict the dissipation scales (η), plots the compensated spectrum as a function of the dimensionless wavenumber kη, and identifies the preferred model as that which displays the better convergence properties at large wavenumbers kη as the Reynolds number increases. Figure 1 of Beresnyak [55] lead him to conclude that it is the −5/3 model that displays the better convergence properties at large k.
It can be shown, however, that the convergence at small scales observed in [54,55] is a simple artifact of the numerical setup adopted in [54,55], rather than a physical effect. To explain this, we note that any discrete numerical scheme solves only the corresponding discrete algebraic equations. If the numerical setup is done correctly, the numerical solution approximates the physical one independently of the discretization step. If, however, a special numerical setup is adopted where η is rigidly tied to the grid size such that ηN is kept the same in all runs (as is done in [54,55]), then the numerical solution plotted as a function of kη is always affected by the discretization in the same way, thus consistently reproducing the same small-scale numerical effects that are present in the setup. The convergence at large kη is then the convergence among solutions of the given numerical scheme, which should not be confused with the convergence to the physical solution.
To illustrate this effect in our simulations we replot the spectra presented in Fig. 8 choosing the Kolmogorov normalization scale η K41 = ν 3/4 ǫ −1/4 . Due to a particular choice of viscosities in our runs depicted in Fig. 8, in this case η K41 happens to double every time the resolution decreases by a factor of 2, thus ensuring that η K41 N ⊥ = const, see Fig. 9. It is therefore not surprising that all the curves converge in the vicinity of the numerical dealiasing cutoff corresponding to kη K41 ≈ 0.8, while they do not converge in the inertial interval and in the most of the dissipation interval. A similar, by design, FIG. 9: The spectra presented in Fig. 8 rescaled with a new parameter ηK41 = ν 3/4 ǫ −1/4 . For this particular choise, ηK41 is proportional to the step of numerical discretization, that is, N ⊥ ηK41 = const. According to our estimate of the onset of the dissipation region in Fig. 8, the corresponding region in the present plot starts at k ⊥ ηK41 ≈ 0.04, while the numerical dialiasing cutoff that is always imposed at k ⊥ = N ⊥ /3, is seen at k ⊥ ηK41 ≈ 0.8. We observe that the convergence is present in the visinity of the dealiasing cutoff, k ⊥ ηK41 0.3, while it is absent in the inertial interval and in most of the dissipation interval.
FIG. 10: Comparison of the spectra obtained in numerical runs RB2c and RB3a, rescaled by the Kolmogorov dissipation scale ηK41 = ν 3/4 ǫ −1/4 (upper panel) and by the dissipation scale obtained in the dynamic alignment theory ηDA = ν 2/3 Λ 1/9 ǫ −2/9 , see (8) (lower panel). In these runs, the numerical discretization is not related to the dissipation scale as N ⊥ η = const. With no spurious numerical convergence imposed by the numerical setup, the −5/3 model does not fit the data, while the dynamic alignment model shows a good agreement in the inertial interval and in the dissipation range, up to the scales where the numerical effects eventually become significant.
Such convergence at very small scales is a spurious numerical effect, which does not reflect the convergence of the physical solutions, and which cannot give preference to any phenomenological model. When the viscosities in different runs do not conform to the special condition ηN ⊥ = const, the spurious convergence disappears, and the −5/3 model does not fit the data, while the −3/2 model still provides a good fit in the inertial and dissipation intervals, see Fig. 10.
We therefore conclude that the numerical simulations by Beresnyak & Lazarian group [53][54][55] are likely signifi-cantly affected by numerical effects at small scales where their measurements are performed. This is notwithstanding the statements made in [53][54][55] that the simulations are resolved in those works. These statements, in our opinion, are not supported by the factual numerical data presented in these papers. Until the effects of hyperdissipation are better understood and numerical convergence is demonstrated in the settings of [53][54][55], it is hard to assess fully the degree to which the numerical findings of [53][54][55] can be compared with our results, or with similar results of other groups [e.g., 16,[47][48][49]51].