# Scaling properties of three-dimensional magnetohydrodynamic turbulence

## Abstract

The scaling properties of three-dimensional magnetohydrodynamic turbulence are obtained from direct numerical simulations of decaying turbulence using 512<sup>3</sup> modes. The results indicate that the turbulence does not follow the Iroshnikov-Kraichnan phenomenology. The spectrum is consistent with $`k^{5/3}`$. In the case of hyperresistivity the structure functions exhibit a clear scaling range yielding absolute values of the scaling exponents $`\zeta _p`$, in particular $`\zeta _31`$, consistent with a recent analytical result. The scaling exponents agree with a modified She-Leveque model $`\zeta _p=p/9+1(1/3)^{p/3}`$, corresponding to Kolmogorov scaling but sheet-like geometry of the dissipative structures.

Magnetic turbulence is the natural state of a plasma in motion, especially in astrophysical systems. The convenient framework to describe such turbulence is magnetohydrodynamics (MHD). For high magnetic Reynolds number $`Rm=vl_0/\eta `$, where $`v`$ is a typical turbulent velocity, $`l_0`$ the integral scale and $`\eta `$ the magnetic diffusivity, there is a broad range of scales $`l`$ between $`l_0`$, and the dissipative scale length $`l_d`$, $`l_0ll_d`$, called the inertial range which exhibits characteristic self-similarity or scaling properties.

The concept of inertial-range scaling was introduced by Kolmogorov  for hydrodynamic turbulence, which is called the Kolmogorov (K41) phenomenology. Assuming homogeneity and isotropy of the turbulence as well as locality of the turbulent cascade process, he obtains $`ϵ(\delta v_l)^2/\tau _l=(\delta v_l)^3/l`$ yielding the scaling law $`\delta v_lϵ^{1/3}l^{1/3}`$. Here $`ϵ`$ is the energy dissipation rate and, to be specific, $`\delta v_l=[𝐯(𝐱+𝐥)𝐯(𝐱)]𝐥/l`$ is the longitudinal velocity increment. A direct consequence is the Kolmogorov energy spectrum $`E_kϵ^{2/3}k^{5/3}`$.

For MHD turbulence the Iroshnikov-Kraichnan (IK) phenomenology ,  takes into account the Alfvén effect, the coupling of small-scale velocity and magnetic fluctuations by the integral-scale field $`B_0`$. Hence the natural variables are the Elsässer fields $`𝐳^\pm =𝐯\pm 𝐁`$, which describe Alfvén waves. In the IK phenomenology the spectral transfer is reduced by the factor $`\tau _A/\tau _l`$, $`ϵ(\tau _A/\tau _l)(\delta z_l)^2/\tau _l`$, where $`\tau _A=l/v_A`$, $`v_A`$ = Alfvén velocity in the field $`B_0`$, $`\delta z_l\delta v_l\delta B_l(ϵv_A)^{1/4}l^{1/4}`$, and the IK energy spectrum becomes $`E_k(ϵv_A)^{1/2}k^{3/2}`$. This spectrum can also be written in the form $`E_kϵ^{2/3}k^{5/3}(kl_0)^{1/6}`$ with the integral scale $`l_0`$ defined by $`l_0=v_A^3/ϵ`$, which illustrates the nonlocal character of the energy cascade in MHD turbulence.

It is, however, well known that these qualitative scaling relations for $`\delta v_l`$ or $`\delta z_l`$ are not exactly valid in a statistical sense because of intermittency, which implies that the distribution of turbulent scales is not strictly self-similar. A quantitative measure is provided by the scaling exponents $`\zeta _p`$ of the structure functions, the moments of the field increments. For hydrodynamic turbulence She and Leveque  proposed a model leading to the expression $`\zeta _p^{\mathrm{SL}}=p/9+2[1(2/3)^{p/3}]`$, which fits the experimental results surprisingly well, reproducing in particular the exact result $`\zeta _3=1`$. This model has been modified for MHD incorporating the IK effect , , which yields $`\zeta _p^{\mathrm{IK}}=p/8+1(1/2)^{p/4}`$, in particular $`\zeta _4^{\mathrm{IK}}=1`$.

The IK phenomenology has been supported by direct numerical simulations of 2D MHD turbulence at moderate Reynolds numbers . However, recent developments in MHD turbulence theory cast some doubt on the general validity of the IK scaling. 2D simulations at considerably higher Reynolds numbers reveal an anomalous scaling behavior , , indicating that the results of Ref.  are not asymptotic. There have also been theoretical arguments in favor of a Kolmogorov scaling, e.g., , . Even more convincingly, exact relations have been derived for moments of certain triple products of $`\delta z_l`$ , which are shown to be proportional to $`l`$, i.e., $`\zeta _3=1`$, analogous to the well-known 4/5-relation in hydrodynamic turbulence, thus excluding the IK result $`\zeta _4=1`$. Scaling exponents for MHD turbulence have also been derived from observations in the solar wind . Here agreement with the IK exponents has been claimed , but in this comparison the observational results were normalized assuming $`\zeta _4=1`$. Actually the error bars seem to be too large to reach a definite conclusion.

To clarify the issue of scaling in 3D MHD turbulence direct numerical simulations are desirable with higher Reynolds numbers than studied previously, for instance in , , , , . In this Letter we present a numerical study of freely decaying turbulence with spatial resolution of $`512^3`$ modes. The scaling properties are analyzed by considering the time-averages of the normalized spectra and structure functions. We solve the incompressible MHD equations

$$_t𝐁\times (𝐯\times 𝐁)=\eta _\nu (1)^{\nu 1}^{2\nu }𝐁,$$
(1)   
$$_t𝐰\times (𝐯\times 𝐰)\times (𝐣\times 𝐁)=\mu _\nu (1)^{\nu 1}^{2\nu }𝐰,$$
(2)   
$`𝐰=\times 𝐯,𝐣=\times 𝐁,𝐯=𝐁=0,`$

by applying a pseudo-spectral method with spherical mode truncation as conveniently used in 3D turbulence simulations (instead of full dealiasing by the 2/3 rule used in most 2D simulations). The generalized magnetic Prandtl number $`\eta _\nu /\mu _\nu `$ has been set equal to unity. Initial conditions are

$$𝐁_𝐤=a\mathrm{e}^{k^2/k_0^2i\alpha _𝐤},𝐯_𝐤=b\mathrm{e}^{k^2/k_0^2i\beta _𝐤},$$
(3)   
which are characterized by random phases $`\alpha _𝐤`$, $`\beta _𝐤`$ and satisfy the conditions $`𝐤𝐁_𝐤=𝐤𝐯_𝐤=0`$, $`E=E^V+E^M=1`$ and $`E^V/E^M=1`$. Further restrictions on $`𝐁_𝐤`$ arise by requiring a specific value of the magnetic helicity $`H=d^3x𝐀𝐁`$. We believe that finite magnetic helicity is more typical than $`H0`$, since MHD turbulence usually occurs in rotating systems. The wavenumber $`k_0`$, the location of the maximum of the initial energy spectrum, is chosen $`k_0=4`$, which allows the inverse cascade of $`H_𝐤`$ to develop freely during the simulation time of 10 time units (about 7 eddy turnover times, defining the eddy turnover time as the time required to reach the maximum dissipation from the smooth initial state). Though this choice implies a certain loss of inertial range, the sacrifice is unavoidable in the presence of inverse cascade dynamics, since assuming $`k_01`$ would lead to magnetic condensation in the lowest-$`k`$ state, which would also affect the dynamics of higher-$`k`$ modes. Both normal diffusion $`\nu =1`$ and hyperdiffusion $`\nu =2`$ have been used, $`\nu =1`$ to discuss the spectral properties and $`\nu =2`$ to determine the scaling of the structure functions. All runs presented in this Letter have finite $`H`$, $`H/H_{\mathrm{max}}0.6`$, and negligible alignment. Table I lists the important parameters of the simulation runs, where the magnetic Taylor Reynolds number is $`Rm_\lambda =Rm^{1/2}`$. Since $`Rm_\lambda `$ is not stricly constant during turbulence decay but increases slowly $`Rm_\lambda t^{1/8}`$, we give the values taken at a specific time $`t=4`$.

We first discuss the spectral properties considering the angle-averaged energy spectrum $`E_k`$. Figure 1 shows the scatter plot of the normalized spectrum (the normalization is discussed below), compensated by $`k^{3/2}`$, taken from run 3 over the period $`t=410`$ of fully developed turbulence. The spectrum exhibits a clear scaling range of almost one decade with a spectral law, which is definitely steeper than the IK spectrum $`k^{3/2}`$, close to (in fact slightly steeper than) $`k^{5/3}`$ indicated by the dashed line. In order to form the time average in a system of decaying turbulence, the spectrum must be normalized to eliminate the time variation of the macroscopic quantities. In hydrodynamic turbulence the only such quantity is $`ϵ`$, which leads to the universal form of the Kolmogorov spectrum $`\widehat{E}(\widehat{k})=E_k/(ϵ\eta ^5)^{1/4}`$, $`\widehat{k}=kl_K`$, where $`l_K=(\eta ^3/ϵ)^{1/4}`$ is the Kolmogorov length. However, when normalized in this way the MHD energy spectrum is found to change during turbulence decay and even more strongly so when comparing runs of different $`Rm_\lambda `$. For finite magnetic helicity the spectrum may also depend on $`H`$, which introduces a second macroscopic length scale $`l_1=H/v_A^2`$ in addition to $`l_0`$, i.e, the spectrum may contain some function of $`l_0/l_1`$. To determine this function we propose the following argument. Since the Alfvén effect is clearly present in the simulations, kinetic and magnetic energy spectrum being nearly equal at small scales, while on the other hand the scaling $`E_kk^{5/3}`$ is observed, we modify the nonlinear transfer in the IK ansatz by a factor $`(l/l_1)^\delta `$,

$$\frac{\tau _A}{\tau _l}\left(\frac{l}{l_1}\right)^\delta \frac{\delta z_l^2}{\tau _l}=ϵ,$$
(4)   
and determine $`\delta `$ by requiring the observed scaling $`\delta z_ll^{1/3}`$, which gives $`\delta =1/3`$ and hence

$$E_kϵ^{2/3}(l_0/l_1)^{1/6}k^{5/3}.$$
(5)   
Also the dissipation scale length is slightly changed. Balancing nonlinear transfer and dissipation gives

$$l_K^H=l_K(l_0/l_1)^{1/8}.$$
(6)   
Using these relations we obtain a the normalized energy spectrum $`\widehat{E}(\widehat{k})=E_k/[(ϵ\eta ^5)^{1/4}(l_0/l_1)^{3/8}]=\widehat{k}^{5/3}F(\widehat{k})`$, $`\widehat{k}=kl_K^H`$. Normalized in this way the spectra at different times of run 3 coincide very well as seen in Fig. 1 and so do the time-averaged normalized energy spectra of runs 1-3 shown in Fig. 2, which vary only by the extent of the inertial range, apart from statistical oscillations. Relations (5), (6) are not valid for $`H0`$, where we expect the pure Kolmogorov normalization to be valid.

A more complete picture of the inertial-range distribution of turbulent structures is provided by the scaling exponents $`\zeta _p`$ of the structure functions, where the second order exponent $`\zeta _2`$ is related to the inertial-range spectral law $`k^{(1+\zeta _2)}`$. To be definite we consider the moments of the absolute value of the longitudinal increments $`\delta z_l^\pm `$ discussing only the runs 4 and 5 with the highest Reynolds numbers. The normalized structure functions $`\widehat{S}_p^\pm (\widehat{l})=|\delta z_l^\pm |^p/E^{p/2}`$, $`\widehat{l}=l/l_K^H`$, are averaged over time, exhibiting a similar weak scatter as for the spectrum in Fig. 1.

For normal diffusion $`\nu =1`$ no scaling range is visible. (Note that the structure function $`S_2`$ corresponds to the one-dimensional spectrum $`E_{k_x}`$, which has a shorter inertial range than the angle-averaged spectrum shown in Fig. 1.) For $`\nu =2`$, however, there is a scaling range $`30<\widehat{l}<200`$, as seen in Fig. 3, where the time-averaged curves $`\widehat{S}_p^+`$ are plotted for $`p=3,4`$. The inserts give the logarithmic derivatives, where the central quasi-constant parts determine the scaling coefficients, the dashed horizontal lines indicating the most probable values $`\zeta _3^+0.95`$, $`\zeta _4^+1.15`$. These results are consistent with the spectral law derived from Fig. 1 and are close to the analytical prediction $`\zeta _3=1`$. It is true that the analytical theory refers to third-order moments different from $`|\delta z|^p`$ discussed here, but the scaling coefficients should not depend thereof (the scaling range, however, does). One might object that the use of hyperdiffusion affects the inertial-range scaling, if the scaling range is not very broad. In fact, the energy spectrum law tends to be polluted by the bottleneck effect, which is particularly pronounced for hyperdiffusion (see e.g., ). Thus the energy spectrum in run 4 (not shown) is effectively flatter than $`k^{5/3}`$ expected from the value of $`\zeta _2`$. However, there is, to our knowledge, no argument for a similar effect in the structure functions.

Assuming the exact result $`\zeta _3=1`$ allows to obtain rather accurate values of $`\zeta _p`$ by using the property of ESS (extended self-similarity)  plotting $`S_p`$ as function of $`S_3`$. (It should be noted that ESS usually results in almost perfect scaling behavior, but the scaling coefficients thus derived vary in time, hence time averaging is required.) The results are shown in Fig. 4, which gives the ESS results of $`\zeta _p^+`$ for run 4 (diamonds) and the ESS values $`\xi _{3p}^+=\zeta _p^+/\zeta _3^+`$ from  for 2D MHD turbulence (triangles). (In 2D MHD the absolute values of $`\zeta _p`$ are found in  to decrease with $`Rm`$, while the relative values $`\xi _{3p}`$ appear to be independent of $`Rm`$.) The results indicate that in 3D MHD turbulence is less intermittent than 2D, but it is more intermittent than hydrodynamic turbulence, the continuous curve, which gives the She-Leveque result $`\zeta _p^{\mathrm{SL}}`$. As shown by Politano and Pouquet , the She-Leveque concept contains effectively three parameters: $`g`$ related to the scaling $`\delta z_ll^{1/g}`$, $`x`$ related to the energy transfer time at the smallest, the dissipative scales $`t_ll^x`$, and $`C`$, the codimension of the dissipative structures,

$$\zeta _p=\frac{p}{g}(1x)+C\left(1(1x/C)^{p/g}\right).$$
(7)   
Our results for the 3D MHD case suggest Kolmogorov scaling $`g=3`$, $`x=2/3`$, while different from hydrodynamic turbulence the dissipative structures are sheet-like, hence the codimension is $`C=1`$,

$$\zeta _p^{\mathrm{MHD}}=p/9+1(1/3)^{p/3}.$$
(8)   
This is the dashed curve in Fig. 4, which fits the numerical values very well.

In conclusion we have studied the spatial scaling properties of 3D MHD turbulence using direct numerical simulations with resolution of $`512^3`$ modes. The results indicate that the turbulence does not follow the Iroshnikov-Kraichnan (IK) phenomenology. The energy spectrum is consistent with a $`k^{5/3}`$ law. For hyperresistivity the structure functions exhibit a clear scaling range yielding absolute values of the scaling exponents $`\zeta _p`$, in particular $`\zeta _31`$, consistent with recent analytical predictions. The scaling exponents agree well with a modified She-Leveque model $`\zeta _p^{\mathrm{MHD}}=p/9+1(1/3)^{p/3}`$, corresponding to Kolmogorov scaling, but sheet-like geometry of the dissipative structures. The results are also consistent with observations of turbulence in the solar wind, which typically show a $`k^{1.7}`$ spectrum.

The authors would like to thank Andreas Zeiler for providing the basic version of the code, Antonio Celani for developing some of the diagnostics, and Reinhard Tisma for optimizing the code for the CRAY T3E.