SCALING LAWS AND DIFFUSE LOCALITY OF BALANCED AND IMBALANCED MAGNETOHYDRODYNAMIC TURBULENCE

MHD turbulence has attracted attention of astronomers since mid-1960s. As most astrophysical media are ionized, plasmas are coupled to the magnetic fields (see, e.g., Biskamp 2003). A simple one-fluid description known as magnetohydrodynamics or MHD is broadly applicable to most astrophysical environments on macroscopic scales. On the other hand, turbulence has been observed in various circumstances and with a huge range of scales (see, e.g., Armstrong et al. 1995; Chepurnov & Lazarian 2010).

MHD turbulence, like hydrodynamics (Kolmogorov 1941), has a "standard" phenomenological model of energy cascade. This is the Goldreich–Sridhar model (Goldreich & Sridhar 1995, hereafter GS95), which uses a concept of critical balance, which maintains that turbulence will stay marginally strong down the cascade. The spectrum of GS95 is supposed to follow a −5/3 Kolmogorov scaling. However, shallower slopes have been reported by numerics, which motivated us to modify GS95 (see, e.g., Boldyrev 2005, 2006; Gogoberidze 2007).

The other problem of GS95 is that it is incomplete, as it does not treat the most general imbalanced, or cross-helical case. As turbulence is a stochastic phenomenon, an average zero cross-helicity does not preclude fluctuations of this quantity in the turbulent volume. Also, most of astrophysical turbulence is naturally imbalanced, due to the fact that it is generated by a strong localized source of perturbations, such as the Sun in the case of solar wind or the central engine in the case of AGN jets.

Several models of imbalanced turbulence appeared recently: Lithwick et al. (2007, hereafter LGS07), Beresnyak & Lazarian (2008), Chandran (2008), Perez & Boldyrev (2009, hereafter PB09), and Podesta & Bhattacharjee (2010). The full self-consistent analytical model for strong turbulence, however, does not yet exist. In this situation, observations and direct numerical simulations (DNS) of MHD turbulence will provide necessary feedback to theorists. We concentrated on two issues, namely, that (1) the energy power-law slopes of MHD turbulence cannot be measured directly from available numerical simulations, supporting an earlier claim in BL09b, and (2) the ratio of energy dissipation rates is a very robust quantity that can be used to differentiate among many imbalanced models. We believe that taken together they resolve confusion related to the subject.

In Section 2 we briefly describe numerical methods, in Section 3 we show and discuss dissipation rates as the most robust measures in numerical turbulence, in Section 4 we argue that it is impossible to measure asymptotic spectral slopes of turbulence directly from currently available DNS, in Section 5 we discuss the perspectives of using DNS to constrain models, and in Section 6 we present our conclusions.

We solved incompressible MHD or Navier–Stokes equations:

$$\begin{equation} \partial _t{\bf w^\pm }+\hat{S} ({\bf w^\mp }\cdot \nabla){\bf w^\pm }=-\nu _n(-\nabla ^2)^{n/2}{\bf w^\pm }+ {\bf f^\pm }, \end{equation} \tag{ 1 }$$

where $\hat{S}$ is a solenoidal projection and w^± (Elsasser variables) are w⁺ = v + b and w⁻ = v − b where we use velocity $\bf {v}$ and magnetic field in velocity units b = B/(4πρ)^1/2. The magnetic Prandtl number here is unity. The Navier–Stokes equation is a special case of Equations (1), where b ≡ 0 and, therefore, both equations are equivalent when w⁺ ≡ w⁻. The right-hand side of this equation includes a linear dissipation term which is called viscosity or diffusivity for n = 2 and hyper-viscosity or hyper-diffusivity for n>2 and the driving force f^±. The special case of f⁺ = f⁻ is a velocity driving. In most simulations from Table 1, we drove w⁺ and w⁻ randomly and independently. The driving consisted of many random components in Fourier space with wavevectors between 2 and 3.5, did not inject magnetic helicity, and is described in greater detail in BL09a. The amplitudes of f⁺ and f⁻ were different in simulations I1–6 (see Table 1), so in these simulations both energy and cross-helicity were injected. This driving supplies Elsasser energies, which physically means that we simulate small scales where independent Elsasser cascades are well established. This contrasts with, e.g., dynamo simulations where only kinetic energy is injected.

Table 1. MHD Simulations

Run	Resolution	f	Dissipation	+/−	(w+)2/(w−)2
B1	1024 · 30722	w±	−3.3 × 10−17k6	∼1	∼1
B2	768 · 20482	v	−3.1 × 10−16k6	∼1	∼1
B3	768 · 20482	w±	−3.1 × 10−16k6	∼1	∼1
B4	768 · 20482	w±	−6.7 × 10−5k6	∼1	∼1
I1	512 · 10242	w±	−1.9 × 10−4k2	1.187	1.35 ± 0.04
I2	7683	w±	−6.8 × 10−14k6	1.187	1.42 ± 0.04
I3	512 · 10242	w±	−1.9 × 10−4k2	1.412	1.88 ± 0.04
I4	7683	w±	−6.8 × 10−14k6	1.412	1.98 ± 0.03
I5	1024 · 15362	w±	−1.5 × 10−15k6	2	5.57 ± 0.08
I6	1024 · 15362	w±	−1.5 × 10−15k6	4.5	45.2 ± 1.5

Download table as:  ASCIITypeset image

We solved these equations with a pseudospectral code that was described in great detail in our earlier publications: BL09a and BL09b. Table 1 enumerates latest high-resolution runs, which were performed in the so-called reduced MHD approximation, where the w^± component parallel to the mean field (pseudo-Alfvén mode) is omitted and so are the parallel gradients in the nonlinear term ((δw^∓ · ∇_||)δw^±. The linear propagation term, i.e., ((δv_A · ∇_||)δw^±, however, is always finite. This is because we consider turbulence injected at such k_|| that v_Ak_|| stays finite; in other words, we use a computational box that is elongated in the x-direction by a factor proportional to v_A. In this situation, the actual value of the mean field B₀ or v_A drops out of the calculations. Physically this means that B₀ is "large enough" compared with perturbations. Under these assumptions one studies purely Alfvénic dynamics in a strong mean field, i.e., Alfvénic turbulence. We compared the properties of this pure Alfvénic turbulence to our earlier full MHD simulations (BL09a, BL09b) with a strong mean field (B₀ = 10, δB ∼ 1) and found that the spectra and anisotropies were very similar, in agreement with analytical argument above.

We started our high-resolution simulations with earlier lower-resolution runs that were evolved for a long time, typically hundreds of Alfvénic times, and reached the stationary state. The balanced runs were evolved, typically for 6–10 Alfvénic times, and the imbalanced runs were evolved for longer times, typically 10–40. The energy injection rates were kept constant in I1–6 and the fluctuating dissipation rate was within few percent of the former. The fluctuations in total energies (w⁺)² and (w⁻)² were the main source of uncertainty for I1–6 shown in Figure 1.

Zoom In Zoom Out Reset image size

One of the most robust quantities in numerical simulations of MHD turbulence is the energy cascading rate or the dissipation rate. In high Reynolds number, turbulence energy has to cascade through many steps before dissipating and the dissipation is negligible on the outer (large) scale. Therefore, the nonlinear energy cascading rate and the dissipation rate are used interchangeably.

In hydrodynamic turbulence, the dissipation rate and the spectrum of velocity are connected with the well-known Kolmogorov constant¹:

$$\begin{equation} E(k)=C_K \epsilon ^{2/3} k^{-5/3}. \end{equation} \tag{ 2 }$$

The important fact that strong hydrodynamic turbulence dissipates in one dynamic timescale l/v is reflected by C_K being close to unity (∼1.6). In MHD turbulence, however, there are two energy cascades (or "Elsasser cascades") and two dissipation rates, ⁺ and ⁻. The question of how these rates are related to velocity-like Elsasser amplitudes w⁺ and w⁻ is one of the central questions of the imbalanced MHD turbulence. Each model of the strong imbalanced turbulence advocates a different physical picture of cascading and provides a different relation between the ratio of energies (w⁺)²/(w⁻)² and the ratio of fluxes ⁺/⁻.

The Goldreich–Sridhar model (GS95) predicts that in the balanced case the cascading is strong and each wave is cascaded by the shear rate of the opposite wave, i.e.,

$$\begin{equation} \epsilon ^+=\frac{\big(w^+_l\big)^2w^-_l}{l},\ \ \ \ \epsilon ^-=\frac{\big(w^-_l\big)^2w^+_l}{l}. \end{equation} \tag{ 3 }$$

It is similar to the Kolmogorov cascade with w's replacing v. Although this model does not make predictions for the imbalanced case, one could hope that in the case of small imbalance these formulae will still work. In this case, we will obtain (w⁺)²/(w⁻)² = (⁺/⁻)². LGS07 argued that this relation will hold even for large imbalances.

For the purpose of this Letter, we mostly discuss the prediction of LGS07, (w⁺)²/(w⁻)² = (⁺/⁻)², and the prediction of PB09 that nonlinear timescales are equal for both waves, which effectively lead to² (w⁺)²/(w⁻)² = ⁺/⁻ (see the corresponding equation in PB09, $w^+/w^-=\sqrt{\epsilon ^+/\epsilon ^-}$). Note that the last prediction is also true for highly viscous flows (Re = Re_m ≪ 1). It could be rephrased that PB09 predict turbulent viscosity which is equal for both components.

Compared to spectral slopes, dissipation rates are robust quantities that require much smaller dynamical range and resolution to converge. Figure 1 shows the energy imbalance (w⁺)²/(w⁻)² versus the dissipation rate imbalance ⁺/⁻ for simulations I2, I4, I5, and I6. We also use two data points from our earlier simulations with large imbalances, A7 and A5 from BL09a. I1 and I3 are simulations with normal viscosity similar to I2 and I4. They show slightly less energy imbalances than I2 and I4 (Table 1).

We see that most data points are above the line which is the prediction of LGS07. In other words, one can deduce that numerics strongly suggest that

$$\begin{equation} \frac{(w^+)^2}{(w^-)^2}\ge \left(\frac{\epsilon ^+}{\epsilon ^-}\right)^2. \end{equation} \tag{ 4 }$$

Although there is a tentative correspondence between LGS07 and the data for small degrees of imbalance, the deviations for large imbalances are significant. Also, the numerics suggests that in the case of small imbalances the cascading smoothly transitions to the balanced case, i.e., the GS95 model, while in the case of strong imbalance it suggests that the strong component cascading rate is smaller than what is expected from strong cascading.

As to PB09 prediction, it is inconsistent with data for all degrees of imbalance including those with small imbalance and normal viscosity, i.e., I1 and I3.

Most of the attention of the theory has been directed toward the self-similar (or approximately self-similar) regime between dissipation and driving, which is colloquially known as "inertial range" (Kolmogorov 1941).

Although attempts to study this regime started long time ago, it is not until recently that simulations with resolution higher than 256³ have become commonplace. At this point, the interpretation of numerical simulations of MHD turbulence has been strongly affected by experience obtained from hydrodynamic simulations.

In hydrodynamics, the correspondence with the Kolmogorov slope has been fairly elusive³ for about a decade, even though the same simulations produced Kolmogorov constants that were close to what has been observed earlier in experiments. The fact that such correspondence has been found in simulations with relatively small, less than a thousand, Reynolds numbers demonstrated that hydrodynamic turbulence is fairly local and in order to reproduce asymptotic cascading only a few steps in log k space are necessary. The energy slopes, however, were affected by the bottleneck effect.

In MHD turbulence, the observed flat (power-law) energy spectra has been prematurely interpreted as "inertial range." However, as it turned out, the flat spectra of MHD turbulence is an indication of a lack of a good inertial range rather than its presence. Indeed, in simulations with hyperdiffusion we compared hydrodynamic and MHD slopes and found that while hydrodynamic energy spectra are highly distorted by the bottleneck effect, the MHD spectra stay very flat (see Figure 2).

Zoom In Zoom Out Reset image size

As it is not known a priori what is the contribution of the systematic error of the spectral slope measurement that comes from the bottleneck effect, it is therefore impossible to measure true asymptotic slopes directly. Also, it is incorrect to claim that the bottleneck effect is absent in simulations with second-order (natural) viscosity, as the existence of the effect was clearly demonstrated in numerical simulations (Kaneda et al. 2003).

Figure 2 shows a comparison between hydrodynamic and MHD energy slopes in 512³ simulations. As we see, the spectra show a variety of bottleneck effects, depending on the order of viscosity and type of simulation (MHD or hydro). Also, there are two types of spectrum, E_k and P_k (see BL09b), and while E_k is used in most numerical papers, it is P_k, a Fourier transform of the structure function, which is directly predicted from the Kolmogorov model. While in the asymptotic regime of exact power-law scaling, P_k and E_k have the same slope, in a realistic numerical simulation they differ quite a lot. From Figure 2, it is not immediately obvious that MHD slopes are shallower than hydro slopes. Most of the publications that made the aforementioned claim had performed only MHD simulations and compared the MHD slope with an asymptotic Kolmogorov slope, i.e., −5/3.

Our study (BL09b) reported that MHD turbulence is less local than hydrodynamic turbulence. This is clearly demonstrated by (1) the lack of visible bottleneck effect in MHD turbulence, while it was clearly present in hydro turbulence (Figure 2) and (2) the dependence of kinetic and magnetic spectra on driving.

An analytical bound for nonlocality can be obtained through the Hölger inequality and scalings of the turbulent fields (Aluie & Eiynk 2010). This bound is shown in Figure 3. From a practical viewpoint, however, this bound does not set a strict constraint on the "width" of the energy transfer window, T(k₀, k), which describes the energy transfer between wavevectors k₀ and k, as the maximally efficient transfer at k₀ could still be much lower than the estimate provided by Hölger inequality. We conclude that, from a practical standpoint, MHD turbulence can still be "diffuse local" i.e., less local than hydrodynamic turbulence despite this analytical bound.

Zoom In Zoom Out Reset image size

Although in this Letter we mostly relied on robust quantities, such as total energies and dissipation rates, we believe that numerical simulations have a wealth of data to be analyzed by theorists. One of the most important measures not mentioned in this Letter is the anisotropy of MHD turbulence. It had been considered in great detail in our earlier publication (BL09a). In particular, we refer the reader to the result of BL08 and BL09a that the anisotropy of the strong component is smaller than the anisotropy of the weak component. This fact is inconsistent with both the naive application of GS95 critical balance (which would have predicted the opposite) and the derivation in LGS07 that suggests that both waves have the same anisotropy. We believe it is incorrect to rely on only one particular measure obtained from simulations, such as the slope of a particular type of the spectrum of the total energy (see Figure 2 for the difference between E_k and P_k) to compare theories and simulations.

PB09 claim that the nonlinear timescales for both components are equal, i.e., there is a turbulent viscosity which is the same for both components, regardless of the degree of imbalance. This seems counterintuitive for transition to freely propagating Alfvénic waves (i.e., infinite imbalance). The formula in PB09, $w^+/w^-=\sqrt{\epsilon ^+/\epsilon ^-}$, suggests that the asymptotic (Re = Re_m ≫ 1) prediction for energy imbalance in this case will be the same as in the highly viscous case (Re = Re_m ≪ 1), i.e., (w⁺)²/(w⁻)² = ⁺/⁻. This is at odds with numerical evidence, which suggests (w⁺)²/(w⁻)² ⩾ (⁺/⁻)².

PB10 claimed that the disagreement between their model and numerics in BL09a is only exhibited in highly imbalanced simulations. This is incorrect. In fact, BL09a studied a range of imbalances γ = w⁺/w⁻ starting with 1 (balanced case) and also 2, 10, and 30. All of them, including one with small imbalance, showed significant inconsistencies with the PB09 model. Also, as we showed in this Letter, the asymptotic regime of very small imbalances shows the same inconsistency with PB09. Unfortunately, PB10 only dealt with the issue of the spectral slope, which is notoriously difficult to measure. We, however, believe that the total dissipation rates present a more robust measure and provide an acid test for any local theory of the imbalanced turbulence.

In BL09a, we discovered and described in detail an empirical fact that one cannot simulate large imbalances with normal (n = 2) dissipation. This empirical fact also suggests that the strong component, w⁺ in our notation, has much larger nonlinear timescale than the weak component w⁻. This is directly supported by the time evolution of w⁺ (Figure 4 of BL09a) and by the shorter inertial range of w⁺ (Figure 10 of BL09a). In a subsequent publication, PB10 found similar empirical evidence, in particular they claim that a resolution of 1024 is required to simulate imbalances of γ ∼ w⁺/w⁻ ∼ 2 and one needs to increase the resolution by a factor of γ^β where β is 2/3 or 3/4. However, this empirical fact directly contradicts to the claim of the PB09 model that both timescales are precisely the same, due to "dynamic alignment." As PB09 predict the same nonlinear timescales for both components, they must have the same dissipation cutoffs. This is contrary to what is observed in numerics.

Motivated by longer nonlinear timescales of the strong component, we used hyperdiffusion in this Letter as well as BL09a. The dissipation for n = 6 hyperdiffusion is a steep function of wavenumber (cutoff scales as ν^−1/(s+n+1), where s is the slope), this allowed us to ensure that dissipation is kept fairly close to the Nyquist frequency even for the strong component, thus allowing us to simulate large imbalances.

In this Letter, we mostly studied the regime of small imbalances in MHD turbulence. In this regime, the imbalanced MHD turbulence is similar to the balanced MHD turbulence in a sense that the cascading will be similar to one in the GS95 model (or LGS07), i.e., both waves will be cascaded strongly and the ratio of energies will be determined by (w⁺)²/(w⁻)² = (⁺/⁻)². Larger imbalances show that (w⁺)²/(w⁻)²>(⁺/⁻)², which suggests that the weak component does not have enough amplitude to provide strong cascading for opposing the strong component (BL08). Also we show that the PB09 model that claims the same nonlinear timescales for both components due to "dynamic alignment" strongly contradicts numerical evidence.

A.B. thanks IceCube project for support of his research and TeraGrid for computational resources. A.L. acknowledges the NSF grant AST-0808118 and support from the Center for Magnetic Self-Organization (LA-UR 10-05945).