{'fraction_non_alphanumeric': 0.036924325420977885, 'fraction_numerical': 0.011709126742602093, 'mean_word_length': 4.689973614775726, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Strong magnetohydrodynamic turbulence with cross helicity

Magnetohydrodynamics (MHD) provides the simplest description of magnetic plasma turbulence in a variety of astrophysical and laboratory systems. MHD turbulence with nonzero cross helicity is often called imbalanced, as it implies that the energies of Alfv\'en fluctuations propagating parallel and anti-parallel the background field are not equal. Recent analytical and numerical studies have revealed that at every scale, MHD turbulence consists of regions of positive and negative cross helicity, indicating that such turbulence is inherently locally imbalanced. In this paper, results from high resolution numerical simulations of steady-state incompressible MHD turbulence, with and without cross helicity are presented. It is argued that the inertial range scaling of the energy spectra (E^+ and E^-) of fluctuations moving in opposite directions is independent of the amount of cross-helicity. When cross helicity is nonzero, E^+ and E^- maintain the same scaling, but have differing amplitudes depending on the amount of cross-helicity.


I. INTRODUCTION
Magnetohydrodynamic (MHD) equations govern the dynamics of a plasma when the spatial scales of interest greatly exceed any intrinsic plasma scale. The MHD model is relevant in understanding waves, instabilities, and turbulence in a variety of astrophysical systems. For instance, magnetohydrodynamic (MHD) turbulence plays a key role in the theoretical modeling of the spectra of velocity and magnetic fluctuations in the Solar Wind [e.g., [1][2][3][4][5] as well as electron density fluctuations in the interstellar medium [e.g., 6,7].
Theoretical modeling of the spectra in strong MHD turbulence has developed on different fronts: phenomenology, statistical closures and more recently, highresolution numerical simulations. The first phenomenological model, a la Kolmogorov [8,9], of the energy spectrum was developed independently by Iroshnikov [10] and Kraichnan [11]. In these works they realized that the magnetic field of the large scale eddies acts in much the same way as a guide field that supports smaller scale Alfvénic fluctuations, and the turbulent energy transfer takes place as the result of many cumulative collisions between counter-propagating Alfvén wave packets traveling along the local magnetic field. One deficiency of this phenomenology is that it is based on the assumption of an isotropic spectral transfer, in contradiction with more recent results that reveal the anisotropic character of MHD turbulence [e.g., 12,13]. To account for anisotropy in the strong turbulence regime, Goldreich and Sridhar [14], hereafter referred as GS, proposed a new phenomenology in which they introduced the socalled critical balance condition, stating that the turbulence is strong when the time scale for wave propagation along the magnetic field matches the nonlinear interaction time. In the past decade, constantly increasing massively parallel computer clusters have become large enough to allow one to simulate the inertial range spectra of strong MHD turbulence, leading to an extensive volume of work addressing universal scaling laws in MHD turbulence [e.g., [15][16][17][18][19][20][21][22][23][24][25]. These simulations have motivated new phenomenological models and have resulted in a renewed interest in the fundamentals of MHD turbulence.
One aspect of MHD turbulence, which has recently drawn significant attention is the role that cross-helicity plays in the turbulent cascade [25][26][27][28][29][30]. Both total energy E = 1 2 (v 2 + b 2 )d 3 x and total cross-helicity H c = (v·b)d 3 x are conserved by nonlinear interactions; where v and b are velocity and magnetic fields. In MHD turbulence both energy and cross helicity cascade from large to small scales as a result of the nonlinear interaction of counter-propagating Alfvén fluctuations. If we denote E ± = 1 4 (v ± b) 2 d 3 x, which are the energies carried by Alfvén fluctuations propagating in opposite directions (the details are given below), the total energy and crosshelicity of the system are given by E = E + + E − and H c = E + − E − , respectively. Therefore, when the latter is nonzero, the energies of waves propagating along and against the guide field are not equal, and in that sense the turbulence is called imbalanced.
Independent numerical simulations by different groups have demonstrated that strong MHD turbulence is always locally imbalanced, in the sense that in the steady state, the turbulence develops correlated regions of positive and negative cross-helicity, irrespective of the overall cross-helicity of the system. Imbalance turbulence is also present in the solar wind, where velocity and magnetic fluctuations show high correlations of a preferred sign, that is, the normalized cross-helicity σ c = v · b /E = H c /E is close to unity. The preferred positive sign of σ c indicates that there is more energy in Alfvén waves propagating outwards from the Sun than propagating inwards [3].
Recently, several phenomenological models have addressed steady-state strong imbalanced MHD turbu-lence [23,[25][26][27]29], some with support from numerical simulations and observations. However, these models have lead to conflicting predictions. For instance, the theories by Lithwick et al [26] and Beresnyak and Lazarian [23] conclude that in imbalanced regions the Elsässer spectra have the same scalings The theory by Chandran [27] proposes that the spectra of E + (k ⊥ ) and E − (k ⊥ ) are different depending of the degree of imbalance. Finally, the analysis by Perez & Boldyrev [25] and Podesta & Bhattacharjee [29] finds that the spectra of E + (k ⊥ ) and E − (k ⊥ ) have different amplitudes but the same scalings In this work we present results from high resolution direct numerical simulations of strong and steadily driven MHD turbulence, aimed to elucidate the role of crosshelicity and resolve the controversies among different theories. Hereafter, strong turbulence is defined in the sense of Golreich and Sridhar [14]. This paper is organized as follows. Section II briefly describes the MHD equations and the relevance of using the Reduced MHD model for strong turbulence simulations. Section III describes in detail the proposed numerical strategy and discusses important numerical aspects that need to be considered when simulating strong imbalanced MHD turbulence. In section IV, the most relevant results from an extensive number of numerical simulations are presented. In section V a phenomenological model for the simulation results is proposed and section VI presents the conclusions.

II. MODEL EQUATIONS
In the presence of a guide field B 0 , the incompressible MHD equations describing the evolution of magnetic and velocity fluctuations, b(x, t) and v(x, t), can be written in terms of the so-called Elsässer variables [31], z ± = v ± b: 4πρ is the Alfvén velocity, ρ is the fluid density, P is the total pressure that is determined from the incompressibility condition, ∇ · z ± = 0, f ± are large-scale forcing, and ν is the viscosity which acts as an energy sink at small scales. In contrast with a uniform flow in hydrodynamic, which can be removed by a Galilean transform, the guide magnetic field, whether imposed externally or created by large-scale fluctuations, is essential in MHD turbulence. It cannot be removed by Galilean transform and it mediates nonlinear interactions at all smaller scales. The linear terms, (v A · ∇)z ± , describe advection of Alfvén wave packets along the guide field, while the nonlinear interaction terms, (z ∓ · ∇) z ± , are responsible for energy redistribution over scales. It can be observed from equations (1) that nonlinear interactions can only occur between z + and z − , and such in-teractions take place when the fields overlap or "collide" with each other. It is convenient to describe the resulting energy spectra in terms of the so-called Elsässer energies, defined as E ± = |z ± | 2 d 3 x/4.
In the limit of small amplitude of fluctuations, the incompressible MHD system (1) describes non-interacting linear Alfvén waves with dispersion relation ω ± (k) = ±k v A . The incompressibility condition requires these waves to be transverse, and they are typically decomposed into two polarizations, the shear-Alfvén wave (z ± S ) and the pseudo-Alfvén waves (z ± P ) given in Fourier space where e is the unit vector in the direction of the guide field. Strong MHD turbulence is dominated by fluctuations with k ⊥ ≫ k . Goldreich and Sridhar [14] argued that since for large k ⊥ the polarization of the pseudo-Alfvén fluctuations is almost parallel to the guide field, such fluctuations are coupled only to field-parallel gradients, which are small since k ≪ k ⊥ . Therefore, the pseudo-Alfvén modes do not play a dynamically essential role in the turbulent cascade. We can remove the pseudo-Alfvén modes by setting z ± = 0 in equations (1) to obtain where R denotes the Reynolds numbers (discussed below). In this system, the fluctuating fields have only two vector components,z ± = {z ± 1 ,z ± 2 , 0} (where we chose the z axis along the guide field B 0 ) but depend on all three spatial coordinates. It can be demonstrated that system (3) is equivalent to the Reduced MHD model, originally developed for tokamak plasmas by [32] and [33].

III. NUMERICAL APPROACH
The reduced MHD model (3) describing the nonlinear interactions of Shear-Alfvén modes is a valuable tool in theoretical and numerical studies of incompressible MHD turbulence. Based on this model, we now discuss the conditions that should be satisfied in order to correctly simulate strong MHD turbulence, both balanced and imbalanced.

A. Computational domain
In the presence of a strong guide field, MHD turbulence is inherently anisotropic. It is important to point out that such anisotropy will develop deep in the inertial range even if it is not present at the outer scale where the turbulence is driven. This can be understood from the Golreich and Sridhar [14] picture of strong turbulence. As the turbulence cascades to smaller scales, eddies become more and more shrunk in the field-perpendicular direction until fluctuations satisfy the critical balance condition k B 0 ∼ k ⊥ b λ , where k ∼ 1/l and k ⊥ ∼ 1/λ are the field-parallel and field-perpendicular wave vectors associated with an anisotropic eddy of parallel size l and perpendicular size λ, respectively, and b λ denoted rms magnetic fluctuations at scale λ. It is assumed that In a traditional cubic simulation box, when an anisotropic wave packet fits the field-parallel direction (z−dimension of the box) its field-perpendicular dimensions are much smaller than the x and y dimensions of the box. This is obviously not an optimal situation, since the field-perpendicular resolution required to resolve such a wave packet should be much higher than the field-parallel resolution. This decreases the effective field-perpendicular resolution of the simulations. For example, in the case a cubic 512 3 box with the guide field B 0 = 5b rms , one will have an equivalent fieldperpendicular resolution of only 512/5 ∼ 100.
For our simulations, we define the nonlinear interaction strength parameter so that the critical balance condition then implies χ ∼ 1.
Simulations of steadily-driven (dissipative) incompressible turbulence are generally based on the Fourier pseudo-spectral method. As with any spectral method, the solution to the differential equation is approximated by a truncated Fourier expansion, and the partial differential equation is converted to a set of ordinary differential equations in time for the N Fourier coefficients. Nonlinear terms in this representation become convolutions whose direct computation requires O(N 2 ) operations. In pseudo-spectral methods, the convolutions are computed in real (or configuration) space by means of a Fast Fourier Transform (FFT) that only requires O(N logN ) operations. For the numerical simulations presented in this work, a third order semi-implicit Runge-Kutta/Crank-Nicholson method was used for the time integration.
In order to allow for the inertial interval to develop, turbulence is driven at the lowest resolvable wave numbers, and the energy dissipates at large wave numbers determined by the Reynolds numbers. For simulations on a cubic periodic box of size L, the smallest wavenumbers along the field-parallel and field-perpendicular directions coincide, i.e., k ⊥ = k = 2π/L. Therefore, driving at the low k , k ⊥ results in an isotropic forcing and the nonlinear strength parameter at the forcing scale becomes which means that at least at the large scales, nonlinear interactions are weak. This would not be harmful if we had the resources to achieve arbitrarily high resolution, as the turbulence would proceed weakly until χ ∼ 1, and then it would become strong. However, simulations generally produce a rather limited inertial range, so that the parameter χ can hardly reach unity in such a setup.
As pointed out by [16], and applied in recent simulations [20,21,24], an effective way to avoid this is to use an anisotropic domain such that at the forcing scale the parameter χ is already of order unity, that is, the excited large-scale modes are already anisotropic and satisfy the critical balance condition. To achieve this we choose an elongated box L 2 ⊥ × L , so that the lowest field-perpendicular and field-parallel wave-numbers are k ⊥ = 2π/L ⊥ and k = 2π/L , respectively. In this case, forcing at the lowest k ⊥ , k leads to χ 0 = (L b λ )/(L ⊥ B 0 ), which is of order unity provided that In this way, the turbulence is excited in a strong regime and the cascade proceeds down to smaller scales preserving the critical balance condition.

B. Numerical Resolution
At first sight, it appears that elongating the box along the z direction to match the elongation of the eddies should not change the number of grid points required in this direction compared to the number of points in the x direction. Fortunately, the number of points in the z direction can be reduced. This follows from the fact that the turbulent spectrum declines quite slowly, as a powerlaw, in the k ⊥ direction, while it drops sharply in the k direction for k > k α ⊥ , where α is a some positive power not exceeding 1 [15,16,24,25,34] . This qualitatively different spectral behavior in k and k ⊥ directions allows one to reduce the numerical resolution by a factor of 2 to 4 in the parallel direction, see Table I. We checked that the restoration of the full resolution in the z direction does not change the results, while significantly increases the computing costs.

C. Periodic Boundary Conditions
The spectral method assumes periodic boundary conditions in all spatial directions. The periodic boundary conditions in the direction of the guide field (the z−direction) may raise the question of whether the magnetic field lines are periodic in such numerical simulations. If they were periodic then the Alfvén modes counter-propagating along a given magnetic field line would repeatedly interact only among themselves, and, therefore, they might not become sufficiently decorrelated between the consequent interactions. To answer this question we note that the numerical setup ensures periodicity of the fluctuations b, but not the magnetic field lines. Magnetic field lines are integral lines of a magnetic field, therefore, they are generally not periodic in the z−direction, since b(k = 0) is generally nonzero. Consequently, an eddy traveling along a magnetic field line interacts with many independent counterpropagating eddies. Reducing the parallel box size L may however lead to artificial overlap of an elongated eddy with itself. This may explain why reducing the parallel box size L below (5) somewhat spoils the spectrum at low wave numbers, as seen in forced simulations of Müller & Grappin [19]. Increasing L beyond (5)  In the inertial interval of turbulence, the guide field should be strong compared to the turbulent fluctuations, b λ ≪ B 0 . It is important to establish how small the fluctuations should be in order to exhibit the universal turbulent spectrum. The question is especially relevant for numerical simulations, as the large guide field implies small Alfvénic time and, therefore, large computing cost. This problem was numerically addressed in [20], where it was found that the transition to the universal regime occurs approximately at b λ /B 0 ∼ 1/3. For b λ /B 0 > 1/3 the energy spectrum is closer to −5/3, possibly indicating that the magnetic field does not qualitatively change the Kolmogorov dynamics for the scales attainable in the numerical simulations. For b λ /B 0 < 1/3, the guide field significantly affects the dynamics, and the spectrum changes to −3/2. The "sweet spot" often used in numerical simulations is b rms /B 0 ∼ 1/5; it has been checked that the smaller ratio, b λ /B 0 ∼ 1/10, does not lead to noticeably different results (see for example [19][20][21]).
It is also important to note that in numerical simulations, where inertial intervals have quite limited extent, the condition b rms /B 0 < 1/3 should be satisfied at the outer scale. Indeed, since b λ decreases with the scale quite slowly, say b λ ∝ λ 1/3 , if the condition b rms /B 0 < 1/3 is not satisfied at the outer scale, it can hardly be satisfied in the inertial interval.

E. Reynolds numbers
Probably the most significant limitation on numerical simulations is imposed by the consideration of imbalanced turbulence. Indeed, in the imbalanced case, γ = z + /z − > 1, the formally constructed Reynolds numbers Re ± λ = λz ∓ λ /ν corresponding to z + and z − fields at some scale λ would be essentially different (here z ± ≡ |z ± |). One can argue that the effective Re number is then the smaller of the two, see below. Therefore, the resolution requirements increase with the amount of imbalance in order to produce large inertial ranges. Assume that the number of grid points in a field-perpendicular direction scales with the Reynolds number as N ∼ Re β , where β = 2/3 or 3/4 depending on the spectral slope (3/2 or 5/3). Then increasing the imbalance γ by, say, a factor of 3 will require increasing the resolution by approximately a factor of 2. Noting that resolution of at least 1024 is required to simulate the imbalance γ ∼ 2 (see below), we conclude that significantly stronger im- balance is not achievable with present day computing power.

F. Random Forcing
Another question that arises in the strong turbulence regime concerns the type of forcing. In real systems, large-scale turbulence can be driven by instabilities exciting both velocity and magnetic harmonics, by external antennae exciting magnetic fields, etc.; this should not matter for the spectrum of small-scale fluctuations. Similarly, in numerical simulations one can force at large scales either velocity or magnetic field, or both of them, this does not change the result [21]. The goal is therefore not to mimic any real-life driving but rather to optimize the transition to the inertial interval.
In this section we discuss the important aspects to be considered when choosing a particular forcing. We assume a random forcef that has no component along z, it is solenoidal in the x − y plane and its Fourier coefficients outside the range are zero, where n z determines the width of the force spectrum in k , and L ⊥ = 2π. The Fourier coefficients inside that range are Gaussian random numbers with amplitude chosen so that the resulting rms velocity fluctuations are of order unity. The individual random values are refreshed independently at time intervals τ . The parameter n z controls the degree to which the critical balance condition is satisfied at the forcing scale. Note that we do not drive the k = 0 mode but allow it to be generated by nonlinear interactions.
In contrast with incompressible hydrodynamic system, an incompressible MHD system can support Alfvén waves. When the system is driven by a time dependent forcing, the most effectively driven modes are those resonating with the frequencies present in the forcing. Therefore, the spatial spectra of the large-scale velocity and magnetic fluctuations are not generally the same as the spatial spectrum of the force. Rather, they essentially depend on the both spatial and temporal spectra of the random forcing. Ultimately, it is the spectrum of the large-scale velocity and magnetic fluctuations, not the driving force, that should be controlled in numerical simulations. To illustrate this we performed a series of simulations in which we drive the modes given by equation (6) with n z = 8. We would expect this to excite all modes from k 0 = (2π/L z ) to 8k 0 . Figure 1 shows the field-parallel energy spectrum of fluctuations at k ⊥ = 3, for differing values of the force correlation time, showing that for short time correlations, the fluid response broadens.
This adds additional complexity as the way in which the turbulence is forced may affect the outcome of the simulation's results. In MHD turbulence there are regimes where the spectrum is expected to obey universal power laws, as well as transition regions where the turbulence changes character, for instance from weak to strong turbulence. Therefore, when simulating MHD turbulence, meaningless results could easily be obtained if the forcing is not chosen carefully. Not optimized forcing makes the transition from the forcing scale to the inertial range unnecessarily longer, as the turbulent fluctuations have to develop the anisotropy implied by critical balance. In the following we perform simulations with a short-time-correlated random forcing that drives the turbulence close to the critical balance condition at the large scales. A short-time correlated Gaussian random force has another important advantage, that it allows one to control the rates at which the energies of z + and z − modes are injected.

G. Inertial range and bottleneck
In spite of the significant growth in massively parallel supercomputers that has occurred over the last years, simulations of MHD turbulence still result in a rather limited inertial range. This generally raises the question of as to whether the simulations fully capture all scales from forcing to dissipation and whether the scaling exponents inferred from the energy spectrum are accurate enough to confirm or rule out the existing models of turbulence. In order to extend the inertial range, some groups have performed simulations with hyper-viscosity of different orders. However, using artificial viscosity might enhance possible bottlenecks (or bumps) that arise at high k due to an abrupt viscous suppression of the turbulent cascade in the dissipation range, e.g., [15,35]. This bottleneck can significantly affect measured scaling laws in simulations with inertial ranges of limited extent. We avoid the use of hyper-viscosity and perform direct numerical simulations of the RMHD equations (3).
In our results we identify the inertial range by performing a set of simulations from low to high Reynolds numbers (with increasing resolution). It is verified that the inertial range becomes larger as the Reynolds number increases. At larger k the inertial interval is followed by the dissipation range where the power-law scaling does not apply.
Spectral power laws in the simulations are detected by compensating the energy spectra with corresponding power laws. This method is an efficient way to compare simulations with existing theories and to observe any deviation from the expected scalings. As another advantage of compensating the energy spectrum, one notices that since the compensating factor increases with k ⊥ , it enhances any possible bottleneck region occurring at large k ⊥ . In the numerical results that we present in this paper, no evidence of a significant bottleneck effect is found, which is consistent with simulations by other groups, e.g., [16,19,20]. Note that for the resolution used in our paper, the corresponding HD simulations would already produce well developed and clearly observed bottleneck region [36].

H. Locality
The question of universality of MHD turbulence is closely related to the question of locality of nonlinear interactions. Locality loosely implies that only fluctuations of comparable length scales interact strongly with each other. More precisely, locality means that the properties of the inertial interval asymptotically become independent of the details of the driving and the dissipation as the Reynolds number increases. Locality is an essential property of hydrodynamic turbulence. In relation to MHD turbulence, it has been discussed in several works, e.g., [37][38][39][40]. Recently, Aluie and Eyink [41] showed both analytically and numerically that, similarly to hydrodynamic turbulence, scale-locality holds in MHD turbulence. This is consistent with the fact that numerical simulations of strong MHD turbulence performed with different forcing and dissipation mechanisms produce the same energy spectrum, e.g., [16,[19][20][21]24].

I. Integration time
Finally, we point out that numerical simulations of imbalanced MHD turbulence require longer integration time in order to accumulate good statistics. This may be related to the fact that the nonlinear interaction is reduced in imbalanced turbulence, see subsection III E. The numerical simulations presented in the next section indicate that for a modest imbalance of about γ ∼ 2 − 3, the relaxation time of the turbulence spectrum is about 20 formally estimated large-scale dynamical times (see below), while about 100 dynamical times are required to    Since the background magnetic field must be strong, we choose B 0 = 5 in the v rms units, (see the discussion in [20]). Time is normalized to the Alfvén transit time τ 0 = L /v A , where L is the field-parallel box size. This time is equivalent to the perpendicular transit time L ⊥ /v rms when v rms = 1. The Reynolds number is defined as Re = v rms (L ⊥ /2π)/ν and we have chosen the same value for the magnetic Reynolds number, Rm = b rms (L ⊥ /2π)/η, denoting both by R in (3). In each run, the average is performed over about 100 largescale-eddy turnover times. The results are presented in Fig. (3). Table I summarizes five representative simulations that incorporate all the aspects discussed in section III. Run A correspond to balanced turbulence, that is, σ c = 0. In runs labeled B, C and D, cross helicity is injected at the forcing scale in such a way that σ c reaches a steady state of ∼ 0.6. We use a short time-correlated forcing, which is on average 1/20 th of the Alfvén time of the excited modes, so that the energy injection rates for both z + and z − only depend on the variance of the imposed forcing, which is controlled in our simulations. In the imbalanced case, field-parallel box size is optimized to reach the critical balance at the large scales. Except for the Reynolds numbers, simulations B, C, and D have the exact same parameters including the energy injection rates, ǫ + and ǫ − . The field-perpendicular energy spectrum is obtained by averaging the angle-integrated Fourier spectrum, over field-perpendicular planes in all snapshots. Figure 3 shows the field-perpendicular energy spectra for each run. Our numerical setup offers significant advantages over full MHD simulations, and can be already seen in simulations of balanced turbulence, top frame in Fig. (3), run A. The energy spectra approach E ± (k ⊥ ) ∝ k −3/2 ⊥ , in good agreement with earlier numerical findings [e.g., 16,17,19,20], but requiring considerably less computational cost and producing slightly larger inertial intervals. Our most significant results are obtained for the imbalanced case. The bottom frame of Fig. 3 shows the spectra for three different Reynolds numbers (runs B, C, D). It is observed that the spectra E ± are pinned at the dissipation scales, which supports the phenomenological predictions by [27,[42][43][44]. We also find that the largescale parts of the both spectra are practically insensitive to the Reynolds numbers. These two important properties imply that as the Re numbers are further increased, the E ± spectra must become progressively more parallel in the inertial interval. This is indeed seen in our numerical simulations (B, C, D). Our numerical simulations indicate that both spectra approach the universal scaling of strong MHD turbulence E ± (k ⊥ ) ∝ k −3/2 ⊥ , while they have essentially different amplitudes and correspond to essentially different energy fluxes.

V. PHENOMENOLOGICAL MODELING
In this section we propose an explanation for the observed spectra. In the case of balanced MHD turbulence, our explanation of the k −3/2 ⊥ energy spectrum relies on the phenomenon of scale-dependent dynamic alignment, see [21,45]. (General details on dynamic alignment in MHD turbulence can be found in, e.g., [18,46]). Consider the eddy shown to the right of Fig. 4, obtained from simulations. In this eddy fluctuations are aligned within the small angle θ λ along x, while their directions and magnitudes change in an almost perpendicular direction, along y. In the case of strong balanced turbulence, the nonlinear interaction in such an eddy is then reduced by a factor θ λ for both z + and z − fields, and the corresponding nonlinear interaction time is estimated as The scaling of the fluctuating fields is then found from the requirement of constant energy fluxes: (z ± λ ) 2 /τ λ = const. One can argue [20,21,45] that the alignment angle decreases with scale as θ λ ∝ λ 1/4 , in which case the field-perpendicular energy spectrum is E(k ⊥ ) ∝ k −3/2 ⊥ . It turns out that the imbalanced simulations of previous section can also be explained within the framework of the scale-dependent dynamic alignment [25,30]. Let us assume that the alignment is present in the imbalanced case. Since the fields amplitudes z + and z − are now essentially different, the alignment angles are dif- ferent as well; we denote them θ + λ and θ − λ , see Fig. 5. These angles obey the important geometric constraint: as is clear from Fig. 5. The depletion of nonlinear interaction [47,48] is therefore different for z + and z − fields, which makes their nonlinear interaction times, τ ∓ λ ∼ 1/(z ± λ k ⊥ θ ± λ ), the same. The requirement of constant energy fluxes (z ± λ ) 2 /τ ± λ ∼ ǫ ± = const then ensures that z + λ /z − λ ∼ ǫ + /ǫ − , so both fields should have the same scaling, although different amplitudes. We also note that the Reynolds numbers that take into account the depletion of nonlinear interaction are the same for both fields, Re ± λ ∼ λz ∓ λ θ ∓ λ /ν, and they are reduced approximately by the factor γ = z + /z − with respect to the formal Reynolds number Re λ = λv λ /ν.
We note that the above relations of the kind z + λ /z − λ ∼ ǫ + /ǫ − hold locally for a particular domain, with given (positive or negative) alignment, while MHD turbulence consists of both positively and negatively aligned domains of various strengths. Since E + and ǫ + are concentrated mostly in positively aligned domains, while E − and ǫ − in negatively aligned ones, the quantities z ± λ and ǫ ± averaged over the global system should not a priori satisfy the same relations. The difference between the local and global quantities in MHD turbulence should be taken into account when one designs numerical tests.
Recently, Beresnyak & Lazarian [23,28] attempted to test the effects of the dynamic alignment in imbalanced turbulence by measuring the relations among global quantities z ± λ and ǫ ± . They did not observe the local scaling relations and concluded that the dynamic alignment proposed in [25,30] is absent. According to the explanation given above, such a conclusion is incorrect. The relations between the global quantities should depend on the details of the distribution of positive and negative domains in the turbulent system, see the discussion in [29]. For our present purposes we simply need the fact that if each domain has the same scaling of the fluctuating fields, the fields averaged over the whole turbulent system will have the same scaling as well.

VI. DISCUSSION.
We have presented a detailed numerical setup based on the Reduced MHD model (RMHD), and results from high resolution simulations of balanced and imbalanced MHD turbulence in steady state. We have used this numerical set up to address currently existing controversies regarding the spectra of imbalanced MHD turbulence.
The simulations are consistent with the theories and observations predicting same scaling for both Elsässer fields [23,25,26,29] and less consistent with the models predicting different scalings for E ± , [e.g., 27]. The measured scaling exponent in the simulations is close to the −3/2 supported by phenomenological models based on dynamic alignment [25,29]. The analysis of our simulation results may also explain somewhat puzzling numerical findings by Beresnyak and Lazarian [23,28], who report different spectra for the E ± Elsässer fields, and the intersection of the spectra rather than pinning at the dissipation scale. According to our results, the explanation might lie in the fact that in these simulations the imbalance was typically high, up to γ 2 = (z + ) 2 /(z − ) 2 ∼ 1000, and, therefore, the universal regime of imbalanced MHD turbulence was not reached, see our analysis in section III E.
The phenomenology of scale-dependent dynamic alignment can be applied to explain the observed spectra. In this phenomenology, the configuration space splits into eddies (domains) with highly aligned and anti-aligned magnetic and velocity fluctuations, where nonlinear interactions are reduced, as in left panel of Fig. 4. Even when the turbulence is balanced overall it still can be imbalanced locally, creating domains of positive and negative cross-helicity. In each of these regions the picture of imbalanced turbulence presented above applies. When averaged over all the regions, the spectra of balanced turbulence are reproduced.
This work was supported by the U.S. DOE Junior Faculty grant DE-FG02-07ER54932, by the DOE grant de-sc0001794, and by the NSF Center for Magnetic Self-Organization in Laboratory and Astrophysical Plasmas at the University of Wisconsin-Madison. High Performance Computing resources were provided by the Texas Advanced Computing Center (TACC) at the University of Texas at Austin under the NSF-Teragrid Project TG-PHY070027T.