SCALING PROPERTIES OF SMALL-SCALE FLUCTUATIONS IN MAGNETOHYDRODYNAMIC TURBULENCE

The energy distribution over scales in magnetic plasma turbulence is the important input ingredient in theories of the interstellar medium (Elmegreen & Scalo 2004; Brandenburg & Nordlund 2011), scintillation of galactic radio sources (Goldreich & Sridhar 2006; Coles et al. 2010), particle heating, and acceleration by magnetic plasma fluctuations in the solar wind (Chandran et al. 2010, 2011). At scales much larger than the plasma micro-scales (such as the ion gyroscale and the ion inertial length) many fundamental aspects of the plasma dynamics can be captured in the framework of magnetohydrodynamics (MHD; e.g., Biskamp 2003; Terry & Tangri 2009; Tobias et al. 2013), which can be effectively studied both analytically and numerically.

In spite of advances in present-day computer simulations, the Reynolds numbers of astrophysical flows (Re ∼ 10⁶–10¹⁶) exceed by many orders of magnitude the Reynolds numbers achieved numerically (Re ∼ 10⁴). In this situation, major interest is attracted to the scaling properties of MHD turbulence, which can be reliably established from numerical simulations. Such an approach motivates phenomenological models that can be extrapolated to astrophysically relevant scales. Examples include models of astrophysical dynamo action (e.g., Malyshkin & Boldyrev 2010; Brandenburg et al. 2012; Kolekar et al. 2012), models of magnetic reconnection at high Lundquist numbers in the solar wind and the solar corona (e.g., Longcope & Sudan 1994; Rappazzo et al. 2008; Ng et al. 2012; Zhdankin et al. 2013), and studies of turbulent mixing in the interstellar medium (e.g., Sur et al. 2014), etc.

The basic assumptions of universality and scale invariance that are common in studies of hydrodynamic turbulence are not as well-justified and, as a result, not well-understood in the MHD case, and they require careful investigation. In the hydrodynamic case, the Kolmogorov first self-similarity hypothesis implies that at scales much smaller than the driving scale, the energy spectrum of incompressible non-magnetized fluid turbulence has a universal form (Kolmogorov 1941; Obukhov 1941):

$$\begin{eqnarray} &&{E(k)}=C_k \epsilon ^{2/3}k^{-5/3} g_h(k\eta _h), \end{eqnarray} \tag{ 1 }$$

where is the mean energy dissipation rate,

$$\begin{eqnarray} &&\eta _h=\nu ^{3/4}\epsilon ^{-1/4} \end{eqnarray} \tag{ 2 }$$

is the Kolmogorov viscous scale, and ν is the fluid viscosity. The function g_h(x) is expected to be universal, that is, independent of the nature of the large-scale driving, and to satisfy g_h(0) = 1. With the driving applied at scale L, in the inertial interval of turbulence kL ≫ 1 ≫ kη_h the function g_h(x) → 1 as η_h → 0, thus leading to the well-known k^−5/3 Kolmogorov's inertial range spectrum. At the dissipation scales, kη_h ≳ 1, the form of the function g_h cannot be derived from scaling arguments. However, it has been constrained by detailed experimental and numerical measurements (e.g., Tsuji 2004; Donzis & Sreenivasan 2010) and phenomenologically modeled (e.g., Monin & Yaglom 1975).

In MHD turbulence, the Kolmogorov self-similarity relation (1) does not apply due to the presence of the Alfvén velocity $v_A=B_0/\sqrt{4\pi \rho }$ associated with the large-scale magnetic field, B₀. The large-scale magnetic field mediates the turbulent dynamics at small scales, therefore, the energy spectrum may essentially depend on the large scale (Iroshnikov 1963; Kraichnan 1965; Goldreich & Sridhar 1995; Ng & Bhattacharjee 1996; Galtier et al. 2000; Bhattacharjee & Ng 2001; Boldyrev 2005, 2006). In this case, the general form of the energy spectrum can be written as

$$\begin{eqnarray} &&E(k)=C^M_k \epsilon ^{2/3} k^{-5/3} g(k\eta, k\Lambda), \end{eqnarray} \tag{ 3 }$$

where η is the dissipation scale and Λ ∼ L is the scale parameter related to the large-scale organization of the flow—both can be different in different setups. The mediation of the small-scale interaction by the large-scale magnetic field implies that in the inertial interval, kη ≪ 1, one cannot require that g(0, kΛ) = 1. Rather, in order to establish the energy spectrum in this case one needs to study the nonlinear interaction of Alfvén wave packets in detail (Iroshnikov 1963; Kraichnan 1965; Goldreich & Sridhar 1995; Ng & Bhattacharjee 1996; Galtier et al. 2000; Bhattacharjee & Ng 2001; Boldyrev 2005, 2006).

We start with the case of balanced weak MHD turbulence, where the average energies in oppositely propagating Alfvén wave packets are the same. We assume a strong uniform background field, B₀ ≫ b_rms, and suppose that turbulence is isotropically excited at scale L, such that $v_{\rm rms}\sim b_{\rm rms}/\sqrt{4\pi \rho }$. The weakness of the interaction follows from the fact that the linear Alfvén frequency, v_A/L, is much larger than the frequency of nonlinear interaction, v_rms/L. It has been derived that the inertial–interval energy spectrum of balanced MHD turbulence scales with the field-perpendicular wavenumber as $k_{\perp }^{-2}$ (e.g., Galtier et al. 2000; Boldyrev & Perez 2009; Wang et al. 2011), which allows us to write the asymptotics of the spectral function in the form

$$\begin{eqnarray} &&g_w(0, k_\perp \Lambda)\sim (k_\perp \Lambda)^{-1/3}. \end{eqnarray} \tag{ 4 }$$

One can demonstrate that Λ ∼ L, and η_w = ν(v_A/Λ)^1/2. Thus, we observe that the inertial interval essentially depends on the outer scale.

A similar consideration applies in the regime of steadily driven balanced strong MHD turbulence. It has been found in numerical simulations and phenomenological models that the field-perpendicular energy spectrum in this case scales as $k_{\perp }^{-3/2}$ (Müller & Grappin 2005; Maron & Goldreich 2001; Haugen et al. 2004; Boldyrev 2005, 2006; Chen et al. 2011; Mason et al. 2008; Perez et al. 2012), which implies the following asymptotic form of the spectral function:

$$\begin{eqnarray} &&g_s(0, k_\perp \Lambda)\sim (k_\perp \Lambda)^{1/6}. \end{eqnarray} \tag{ 5 }$$

Here, Λ is related to the cross-helicity of the flow (for example, Λ ∼ L, if the magnetic and velocity fluctuations are driven at the outer scale in a non-correlated fashion), and

$$\begin{eqnarray} &&\eta _s = \nu ^{2/3} \Lambda ^{1/9} \epsilon ^{-2/9} \end{eqnarray} \tag{ 6 }$$

is the dissipation scale; see, e.g., Perez et al. (2012). In the imbalanced case, it has also been shown that the asymptotic form of the spectrum follows Equation (4) in the weak imbalanced case (Boldyrev & Perez 2009) and Equation (5) in the strong imbalanced case (Perez & Boldyrev 2009; Podesta & Bhattacharjee 2010).

Both examples demonstrate that the dependence on Λ in MHD turbulence is crucial for establishing the energy distribution in both the inertial and dissipation intervals. In both cases, the spectrum deviates from the Kolmogorov, k^−5/3, due to a reduction of the nonlinear interaction by a certain mechanism related to the large-scale magnetic field. In the case of weak turbulence, such a mechanism is the decorrelation of the triple-field products due to the short crossing time of counter-propagating Alfvén waves (e.g., Ng & Bhattacharjee 1996; Galtier et al. 2000; Bhattacharjee & Ng 2001). In the case of strong turbulence, a weakening of the nonlinear interaction is provided by the scale-dependent angular alignment between magnetic and velocity fluctuations, that is, progressive "Alfvénization" of the turbulence at small scales (e.g., Boldyrev 2005, 2006; Mason et al. 2008). As we argue in the next section, the dependence of the spectral function g_s(kη, kΛ) on the outer scale is crucial for the applicability of discrete numerical schemes for simulations of MHD turbulence.

In this section, we concentrate on strong MHD turbulence, and assume that the simulations are performed in a numerical scheme discretized at scale Δ, which can be the grid size of a finite-difference scheme, the inverse dealiasing cut-off of a pseudo-spectral scheme, etc. In the presence of a numerical cutoff, Δ, the general form of the function g_s is g_s(kη, kΛ, kΔ). The solution of the discrete scheme in general is different from the physical solution and it may have different scaling properties as it contains an additional dimensional parameter, Δ. However, it needs to converge to the physical solution as Δ → 0.

It should be recalled that in the hydrodynamic case, when the spectrum is independent of Λ, the g function can be written as g_h(kη_h, Δ/η_h). Therefore, as long as the numerical resolution is a fixed fraction of the dissipation scale, Δ/η_h = const, the form of the function g_h(kη_h) is universal (e.g., Gotoh 2002). In contrast, the presence of an additional scale, Λ, in the MHD case means that there are infinitely many g_s functions that provide the same inertial interval asymptotics, but have different behavior at the subrange scales. Consider, for example, a spectral function⁶

$$\begin{eqnarray} &&g_s =\left[\frac{1}{k\Lambda }+\frac{\Delta }{\sqrt{\eta _s\Lambda }}\right]^{-1/6}g_1(k\eta _s, k\Delta), \,\,\,\,\, \end{eqnarray} \tag{ 7 }$$

where g₁(0, 0) = 1. The expression (7) agrees with the inertial scaling of −3/2, while it steepens at small scales due to the finite discretization cutoff Δ. In this example, the scaling of the numerically measured energy spectrum changes at scale $k\sim \sqrt{\eta _s/\Lambda }(1/\Delta)$ from −3/2 to −5/3; however, this spectral steepening represents the property of the numerical scheme, not of the physical solution. In order to avoid the influence of numerical effects, one needs to ensure that the discretization cutoff is sufficiently small. For instance, in order to observe the inertial interval up to the dissipation scale k ∼ 1/η_s, one needs to require that $\Delta <\eta _s^{3/2}/\Lambda ^{1/2}$.

In this example, the discretization scale in the numerical simulations needs to decrease faster than the dissipation scale, for the simulations to be resolved. In the next section, we demonstrate that a similar situation is encountered in simulations of MHD turbulence.

The incompressible MHD equations can be written in the Elsässer form:

$$\begin{eqnarray} \left(\frac{\partial }{\partial t}\mp \boldsymbol {v}_A\cdot \nabla \right) \boldsymbol {z}^\pm +\left(\boldsymbol {z}^\mp \cdot \nabla \right)\boldsymbol {z}^\pm &=& -\nabla P +\nu \nabla ^2 \boldsymbol {z}^{\pm }+\boldsymbol {f}^\pm, \nonumber \\ \nabla \cdot {\boldsymbol {z}}^{\pm }&=&0, \end{eqnarray} \tag{ 8 }$$

where $\boldsymbol {z}^\pm =\boldsymbol {v}\pm \boldsymbol {b}$ are the Elsässer variables, and $\boldsymbol {v}$ and $\boldsymbol {b}$ are the fluctuating velocity and magnetic fields in units of the Alfvén velocity, ${\boldsymbol v}_A={\boldsymbol B}_{\boldsymbol {0/}}\sqrt{4\pi \rho _0}$. In these equations, P = (p/ρ₀ + b²/2), where p is the plasma pressure, ρ₀ is the background plasma density, and ${\boldsymbol \nu}$ is the fluid viscosity. For simplicity, the viscosity is equal to the magnetic diffusivity. The turbulence is driven at large scales by the forces $\boldsymbol {f}^\pm$. In the linear case, the plasma waves can be decomposed into shear Alfvén waves whose polarizations are perpendicular to both B₀ and the wave vector, $\boldsymbol {k}$, and pseudo-Alfvén waves whose polarizations are in the plane of B₀ and $\boldsymbol {k}$, and perpendicular to $\boldsymbol {k}$.

In the case of strong MHD turbulence, the pseudo-Alfvén modes are dynamically irrelevant for the turbulent cascade (e.g., Goldreich & Sridhar 1995). One can therefore filter out the pseudo-Alfvén modes by setting $\boldsymbol {z}^\pm _\Vert =0$, which reduces the equations to the Reduced MHD model:

$$\begin{eqnarray} &&\left(\frac{\partial }{\partial t}\mp \boldsymbol {v}_A\cdot \nabla _\Vert \right)\boldsymbol {z}^\pm _{\perp }+\left(\boldsymbol {z}^\mp _{\perp }\cdot \nabla _\perp \right)\boldsymbol {z}^\pm _{\perp } = -\nabla _\perp P+\nu \nabla ^2\boldsymbol {z}^\pm _{\perp } +\boldsymbol {f}_\perp ^\pm. \nonumber\\ \end{eqnarray} \tag{ 9 }$$

We note that in RMHD the fluctuating fields have only two vector components, but that each depends on all three spatial coordinates. Due to incompressibility, each field has only one degree of freedom, which can be expressed in terms of stream functions in the more standard form of the RMHD equations (Strauss 1976). The equivalence between RMHD and MHD in the strong turbulence regime has been shown in numerical simulations; for an extensive discussion see Mason et al. (2012). In light of this equivalence, we will refer to the numerical spectrum obtained from RMHD simulations as the MHD spectrum. We solve the RMHD Equations (9) in a periodic, rectangular domain with aspect ratio $L_{\perp }^2 \times L_\Vert$, where the subscripts denote the directions perpendicular and parallel to ${\boldsymbol {B_0}}$, respectively. We set L_⊥ = 2π, $L_\Vert /L_\perp =6$, and B₀ = 5 e_z. A fully dealiased three-dimensional pseudo-spectral algorithm is used on a grid with a resolution of $N_{\perp }^2\times N_\Vert$ mesh points.

Both Elsässer variables are driven by independent random forces ${\boldsymbol f^+}$ and ${\boldsymbol f^-}$ applied in Fourier space at wave numbers 2π/L_⊥ ⩽ k_{x, y} ⩽ 2(2π/L_⊥), $k_\Vert = 2\pi /L_\Vert$. The forces have no component along z and they are solenoidal in the xy plane. Their Fourier coefficients 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 time of the large-scale eddies. The variances $\sigma _{\pm }^2=\langle |\boldsymbol {f}^{\pm } |^2\rangle$ control the average rates of energy injection into the z⁺ and z⁻ fields. In this work, we discuss the regime of balanced MHD turbulence, i.e., σ⁺ ≈ σ⁻. The Reynolds number is defined as Re = v_rms(L/2π)/ν.

The spectra of MHD turbulence obtained in the simulations are shown in Figure 1. The solid line represents the energy spectrum for a 512³, Re = 2400 run, which is well-resolved. The dash–dotted line shows the same set up with decreased viscosity, which makes the dissipation scale approach the numerical cut-off. The scales at k > 15 are now significantly affected by the proximity to the dealiasing cut-off (see Perez et al. 2012). The proximity to the k-space cutoff is known to distort the spectral behavior at small scales in hydrodynamic simulations (e.g., Cichowlas et al. 2005; Frisch et al. 2008; Connaughton 2009; Grappin & Müller 2010); our MHD runs bear similarity with those results. Such a spectral distortion is a property of the numerical scheme and it should not be confused with the inertial interval behavior. Indeed, as one increases the number of grid points to 1024³ (thereby reducing Δ to Δ/2), while keeping all the physical parameters unchanged, the numerical distortion disappears, see the dashed line in Figure 1.

Zoom In Zoom Out Reset image size

We now note that in the resolved runs, the dissipation interval scales according to the formula (6), see Figure 2 and also Perez et al. (2012). However, in the unresolved runs, the scaling of the small-scale spectrum is different and is consistent with expression (2), as is seen in Figure 3. This is not surprising; as was discussed in the previous section, the scaling of the numerical solution close to the discretization cutoff may be different from the scaling of the physical solution if the former does not approximate the latter. Such a scaling of unresolved runs was previously observed in Beresnyak (2011, 2012, 2014), where it was incorrectly attributed to the scaling of the physical solution because the numerical convergence was not checked.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The presented results reveal an important property of the dissipation interval of MHD turbulence. As we have seen, two scales play a role here: the physical dissipation scale, η_s, given by Equation (6), and the scale, η_h, given by Equation (2). Suppose that the discretization scale, Δ, in numerical simulations is decreased proportionally to the scale, η_s, as the Reynolds number increases. Then, since the other scale, η_h, decreases faster, at some Reynolds number, we will have Δ ∼ η_h and the simulations will become unresolved. In this case, the numerical scheme will not approximate the physical solution at scales close to η_h. If starting from this point, Δ starts to decrease proportionally to η_h, the spectral distortion will be preserved, and the approximation will not improve (see Figure 3). However, if Δ continues to decrease slower than η_h, the approximation of the numerical scheme will continue to degrade. We therefore conclude that in order to resolve the dissipation interval of MHD turbulence, the discretization of the numerical scheme should satisfy Δ ≪ η_h, and it should decrease at least as fast as η_h to maintain the same level of approximation.

The physical explanation of the observed phenomenon is the following. As proposed in Boldyrev (2005, 2006), the −3/2 scaling of the inertial interval of MHD turbulence is related to the scale-dependent dynamic alignment between magnetic and velocity fluctuations, where the alignment angle decreases with the scale approximately as θ(r)∝r^1/4. Numerical simulations demonstrate that such alignment is present not only in the inertial interval, but also in the dissipation range as long as we do not approach the numerical cut-off too closely (Mason et al. 2008; Perez et al. 2012). This alignment leads to the scale-dependent weakening of the nonlinear interaction, changing the naive −5/3 scaling of the equation to the observed −3/2 scaling of the solution.

In the unresolved intervals, however, the correlation between magnetic and velocity fluctuations gets destroyed due to the proximity to the dealiasing cut-off. A similar behavior has been documented in numerical simulations of hydrodynamic turbulence, where an abrupt cutoff in the k space led to unphysical distortion of the energy spectrum (Cichowlas et al. 2005; Frisch et al. 2008; Connaughton 2009; Grappin & Müller 2010). In the extreme case, when the explicit dissipation was totally absent, the spectral cutoff would lead to an energy pile up and mode thermalization close to the cutoff scale. At such scales any spatial correlation of fluctuations existing in the turbulent flow would be lost.

It is reasonable to propose that a similar tendency exists in MHD turbulence when the physical dissipation is not strong enough in the vicinity of the numerical cutoff. This leads to the unphysical numerical distortion seen in Figure 3 and similar unphysical effects in the simulations of Beresnyak (2011, 2012, 2014). In this case, the correlation properties of the fluctuations can be affected as well. In particular, the angular alignment of magnetic and velocity fluctuations can weaken. This is indeed consistent with Figure 4, which illustrates the scaling of the alignment angle between the velocity and magnetic fluctuations measured at scale r (see the details of the measurement procedure in Perez et al. 2012). In the resolved run (dashed line) the alignment is well-preserved down to the scale r ∼ 0.005, corresponding to the scaling −3/2 observed in Figure 1. In contrast, in the unresolved runs (solid and dash–dotted lines), the alignment is already significantly spoiled at scale r ∼ 0.05, which is consistent with the transition to the −5/3 scaling at smaller scales (compare to Figure 3).

Zoom In Zoom Out Reset image size

In conclusion, we have demonstrated that the Kolmogorov first self-similarity hypothesis (1) does not apply to MHD turbulence. In contrast with the hydrodynamic case, the Fourier energy spectrum is not a universal function of the wavenumber normalized by the dissipation scale, kη_s. This happens due to the mediation of the small-scale turbulent energy cascade by the large-scale magnetic field, which introduces the large-scale dependence in both the inertial and the dissipation intervals. We have proposed that the lack of universality has implications for numerical simulations of MHD turbulence. In particular, when the explicit physical dissipation is not strong in the vicinity of the numerical cutoff, the measured spectrum gets distorted by the numerical effects. The solution of the discrete numerical scheme in the cutoff vicinity has a Reynolds number scaling that is different from the scaling of the physical solution. This imposes stringent conditions on the numerical resolution required for correct simulations of MHD turbulence. The cutoff scale should decrease faster than the dissipation scale as the Reynolds number increases in order for the numerical simulations to be adequately resolved.

This work was supported by NSF/DOE grant AGS-1003451, grant NNX11AJ37G from NASA's Heliophysics Theory Program, DoE grant DE-SC0003888, NSF grant NSF PHY11-25915, and the NSF Center for Magnetic Self-Organization in Laboratory and Astrophysical Plasmas. High Performance Computing resources were provided by the Texas Advanced Computing Center (TACC) at the University of Texas at Austin under the NSF-XSEDE Project TG-PHY110016.