# Decay laws for three-dimensional magnetohydrodynamic turbulence

## Abstract

Decay laws for three-dimensional incompressible magnetohydrodynamic turbulence are obtained from high-resolution numerical simulations using up to $`512^3`$ modes. For the typical case of finite magnetic helicity $`H`$ the energy decay is found to be governed by the conservation of $`H`$ and the decay of the energy ratio $`\mathrm{\Gamma }=E^V/E^M`$. One finds the relation $`(E^{5/2}/ϵH)\mathrm{\Gamma }^{1/2}/(1+\mathrm{\Gamma })^{3/2}=const`$, $`ϵ=dE/dt`$. Use of the observation that $`\mathrm{\Gamma }(t)E(t)`$ results in the asymptotic law $`Et^{0.5}`$ in good agreement with the numerical behavior. For the special case $`H=0`$ the energy decreases more rapidly $`Et^1`$, where the transition to the finite-$`H`$ behavior occurs at relatively small values.

Many plasmas, especially in astrophysics, are characterized by turbulent magnetic fields, the best-known and most readily observable example being the solar wind. The convenient framework to describe such turbulence is magnetohydrodynamics (MHD). Here one ignores the actual complicated dissipation processes, which occur on the smallest scales and would usually require a kinetic treatment, assuming that the main turbulent scales are essentially independent thereof. Instead dissipation is modeled by simple diffusion terms. If, moreover, interest is focussed on the intrinsic turbulence dynamics, one can also ignore the largest scales in the system, which depend on the specific way of turbulence generation, restricting consideration to a small open homogeneous domain of the globally inhomogeneous turbulence.

Homogeneous MHD turbulence has become a paradigm in fundamental turbulence research, which has been receiving considerable attention. It is well known that 2D and 3D MHD turbulence have many features in common concerning, in particular, the cascade properties. In both cases there are three quadratic ideal invariants: the energy $`E=\frac{1}{2}(v^2+B^2)𝑑V`$, the cross helicity $`K=𝐯𝐁𝑑V`$, and a purely magnetic quantity, the magnetic helicity $`H=𝐀𝐁𝑑V`$ in 3D and the mean-square magnetic potential $`H^\psi =\psi ^2𝑑V`$ in 2D, which both exhibit an inverse cascade. Many theoretical predictions do not distinguish between 2D and 3D, concerning, e.g., the tendency toward velocity and magnetic field alignment or the spectral properties. Thus it is not surprising, that numerical studies of MHD turbulence have been mostly concentrated on two-dimensional simulations, where high Reynolds numbers can be reached much more readily, see e.g., . While 2D simulations are now being performed with up to $`N^2=4096^2`$ modes (or, more accurately, collocation points) , studies of 3D MHD homogeneous turbulence have to date been restricted to relatively low Reynolds numbers using typically $`N^3=64^3`$ modes, e.g., , , precluding an inertial range scaling behavior. Also in Ref. , where a somewhat higher Reynolds number could be reached by using $`180^3`$ modes, attention was focussed primarily on the process of turbulence generation from smooth initial conditions and the properties of the prominent spatial structures, current and vorticity sheets.

In this Letter we present results of a numerical study of freely decaying 3D MHD turbulence with spatial resolution up to $`512^3`$ modes. We discuss the decay laws of the integral quantities, in particular the energy $`E`$ and the ratio of kinetic and magnetic energies $`\mathrm{\Gamma }=E^V/E^M`$, and their dependence on the quasi-constant value of $`H`$. The energy decay is found to follow a simple law, which is determined by $`\mathrm{\Gamma }(t)`$ and $`H`$. While most previous studies have been restricted to the case of negligible magnetic helicity $`H0`$, we focus attention on the properties of the turbulence for finite $`H`$, which is more typical for naturally existing MHD turbulence occuring mostly in rotating systems. We find that for finite $`H`$ the energy decays significantly more slowly than for $`H0`$. This behavior is primarily caused by the rapid decrease of the energy ratio $`\mathrm{\Gamma }`$, which has the same decay time as the energy.

The 3D incompressible MHD equations, written in the usual units,

$$_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 𝐁,`$

are solved in a cubic box of size $`2\pi `$ with periodic boundary conditions. The numerical method is a pseudo-spectral scheme with spherical mode truncation as conveniently used in 3D turbulence simulations (instead of full dealiasing by the 2/3 rule chosen in most 2D simulations). 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`$ as well as $`E=1`$ and $`\mathrm{\Gamma }=1`$. Further restrictions on $`𝐁_𝐤`$ and $`𝐯_𝐤`$ arise by requiring specific values of $`H`$ and $`K`$, respectively. The wavenumber $`k_0`$, the location of the maximum of the initial energy spectrum, is chosen as $`k_0=4`$, which allows the inverse cascade of $`H_𝐤`$ to develop freely during the simulation time of 10-20 eddy turnover times. This implies a certain loss of inertial range, i.e., a reduction in Reynolds number, but the sacrifice is unavoidable in the presence of inverse cascade dynamics. Choosing $`k_01`$ would lead to magnetic condensation in the lowest-$`k`$ state, which would affect the entire turbulence dynamics. We have used both normal diffusion $`\nu =1`$ and hyperdiffusion $`\nu =2`$. Apart from the fact that inertial ranges are wider and $`H`$ is much better conserved for $`\nu =2`$ than for $`\nu =1`$, no essential differences are found between the two cases. The generalized magnetic Prandtl number $`\eta _\nu /\mu _\nu `$ has been set equal to unity. Table I lists the most important parameters of the simulation runs.

The energy decay law is a characteristic property of a turbulent system. In hydrodynamic turbulence the decay rate depends on the energy spectrum at small $`k`$. Assuming time invariance of the Loitsianskii integral $`=_0^{\mathrm{}}𝑑ll^4v_l(x+l)v_l(x)`$ the energy has been predicted to follow the similarity law $`Et^{10/7}`$ . The invariance of $``$ has, however, been questioned, see e.g., . Both closure theory  and low-Reynolds number simulations  yield a significantly slower decrease, $`Et^1`$. Experimental measurements of the energy decay law $`t^n`$ are rather difficult and do not give a uniform picture, $`n`$ ranging between 1.3  and 2 .

The invariance of the Loitsianskii integral has recently also been postulated for MHD turbulence , where $`_{\mathrm{MHD}}`$ is defined in analogy to $``$ in terms of the longitudinal correlation function $`z_l^\pm (x+l)z_l^\pm (x)`$ of the Elsaesser fields $`𝐳^\pm =𝐯\pm 𝐁`$. Since $`z^2E`$, this assumption gives $`_{\mathrm{MHD}}L^5E=const`$, where $`L`$ is the macroscopic scale length of the turbulence. In addition the expression for the energy transfer $`dE/dt=ϵz^4/LB_0`$ was used, which formally accounts for the Alfvén effect ,. These relations give $`(dE/dt)B_0/E^{11/5}=const`$ and hence $`Et^{5/6}`$, treating $`B_0`$ as constant. One may, however, argue that the Alfvén effect is only important on small scales $`lL`$, while on the scale $`L`$ of the energy-containing eddies $`B_0`$ is not constant but $`B_0E^{1/2}`$ (except for the case that $`B_0`$ is an external field, which would, however, make the turbulence strongly anisotropic), hence $`ϵE^{3/2}/L`$, which would give the same result $`n=10/7`$ as predicted for hydrodynamic turbulence. Low-resolution numerical simulations  indicate $`n1`$, which is also found in recent simulations of compressible MHD turbulence .

For finite magnetic helicity $`H`$ provides a constant during energy decay, which for high Reynolds number is more robust than the questionable invariance of the Loitsianskii integral. It is true that in contrast to the 2D case, where $`E^M`$ and $`H^\psi `$ are tightly coupled, such that $`E^M0`$ implies $`H^\psi 0`$, in 3D a state with $`H=0`$ and finite magnetic energy is possible. But this is only a special and not typical case, since in nature magnetic turbulence usually occurs in rotating systems, which give rise to finite magnetic helicity.

If the process of turbulence decay is self-similar, which also implies that the energy ratio $`\mathrm{\Gamma }`$ remains constant, the energy decay law follows from a simple argument . With the scale length $`L=E^{3/2}/ϵ`$, the dominant scale of the energy-containing eddies, we have

$$HE^MLEL,$$
(4)   
since owing to the assumed self-similarity $`E^ME^VE`$. Inserting $`L`$ gives

$$\frac{dE}{dt}=ϵ\frac{E^{5/2}}{H},$$
(5)   
which has the similarity solution $`Et^{2/3}`$. In Fig. 1 the ratio $`E^{5/2}/(ϵH)`$ is plotted for the runs from Table I with $`H0`$ and small initial correlation $`\rho _0`$. The figure shows that this quantity is not constant, but increases in time. Moreover, there is a significant scatter of the different curves. Integration yields a slower asymptotic energy decay than predicted $`n0.50.55`$. (The log-log representation of $`E(t)`$, often given in the literature to make a power law behavior visible, is misleading, since the major part of such a curve refers to the transition period of turbulence generation. The solution $`(tt_{})^n`$ approaches the power law $`t^n`$ only asymptotically for $`tt_{}`$, where $`t_{}`$ is not accurately known. We therefore prefer to plot the decay law in the primary differential form.)

We can attribute this discrepancy to the fact that the turbulence does not decay in a fully self-similar way. Indeed the energy ratio $`\mathrm{\Gamma }`$ is found to decrease rapidly, in contrast to the 2D case, where $`\mathrm{\Gamma }`$ decays much more slowly, typically logarithmically , . (The ratio of viscous and resistive dissipation $`ϵ^\mu /ϵ^\eta `$, however, remains constant just as in the 2D case , which simply reflects the basic property, that dissipation takes place in current sheets and that these are also vorticity sheets, i.e., the location of viscous dissipation.) Let us incorporate the dynamic change of $`\mathrm{\Gamma }`$ in the theory of the energy decay. Assuming that the most important nonlinearities arise from the $`𝐯`$ contributions in the MHD equations, Eq. (5) is replaced by

$$ϵ(E^V)^{1/2}\frac{E}{L}=\frac{\mathrm{\Gamma }^{1/2}}{(1+\mathrm{\Gamma })^{3/2}}\frac{E^{5/2}}{H},$$
(6)   
using the relation (4). Figure 2 shows that $`(E^{5/2}/ϵH)\mathrm{\Gamma }^{1/2}/(1+\mathrm{\Gamma })^{3/2}`$ is indeed nearly constant for $`t>2`$, when turbulence is fully developed, and the scatter in Fig. 1 is strongly reduced. Hence relation (6) is generally valid for finite magnetic helicity. It is also independent of the magnitude of the dissipation coefficients and character of the dissipation ($`\nu =1`$ or 2), as long as $`H`$ is well conserved.

Also the time evolution of the energy ratio $`\mathrm{\Gamma }`$ exhibits a uniform behavior which is demonstrated in Fig. 3. The slight shift of the uppermost curve corresponding to the smallest value of $`H`$ (run 4), is due to the smaller drop of $`\mathrm{\Gamma }`$ during the very first phase of turbulence generation $`t<0.5`$ not included in the figure. Moreover, we find that $`\mathrm{\Gamma }(t)`$ is proportional to $`E(t)`$, $`\mathrm{\Gamma }cE/H`$, $`c=0.10.15`$, as seen in Fig. 4, where $`\mathrm{\Gamma }/(E/H)`$ is plotted. Inserting this result in Eq. (6) we obtain the differential equation for $`E`$, which in the asymptotic limit $`\mathrm{\Gamma }1`$ becomes

$$\frac{dE}{dt}0.5\frac{E^3}{H^{3/2}}$$
(7)   
with the similarity solution $`Et^{0.5}`$. For finite $`\mathrm{\Gamma }`$ the theory predicts a somewhat steeper decay flattening asymptotically to $`t^{0.5}`$ as $`\mathrm{\Gamma }`$ becomes small, which is exactly the behavior of $`E(t)`$ observed in the simulations. (Note, that if $`E(t)`$ is plotted on the traditional log-log scale, which overrates the transition period $`t1`$, a steeper decay would be suggested.) The relation $`\mathrm{\Gamma }E`$ now gives also the similarity law for the kinetic energy $`E^Vt^1`$.

This theory does not apply to the special case $`H=0`$. Here we find indeed a different decay law, $`Et^1`$ from run 3, which is consistent with previous simulations at lower Reynolds numbers  and with the prediction in Ref. . The transition to the slower decay for finite $`H`$ occurs at relatively small values, 0.1–0.2 of the maximum possible value.

We have also studied the effect of an initial velocity and magnetic field alignment $`\rho _0=K/E`$. For small $`\rho _0<0.1`$ the alignment, after increasing initially, tends to saturate at some small value, which is due to the fact that $`K`$ is less well conserved than $`H`$. For higher $`\rho _0>0.3`$ (runs 9 and 10 in Table I) the alignment becomes very strong, which as expected slows down the energy decay drastically.

In conclusion we have presented a new phenomenology of the energy decay in 3D incompressible MHD turbulence, which agrees very well with direct numerical simulations at relatively high Reynolds numbers. We consider in particular the case of finite magnetic helicity $`H`$, which is typical for naturally occuring magnetic turbulence. The energy decay is governed by the conservation of $`H`$ and the time evolution of the energy ratio $`\mathrm{\Gamma }=E^V/E^M`$. We find that the relation $`(E^{5/2}/ϵH)\mathrm{\Gamma }^{1/2}/(1+\mathrm{\Gamma })^{3/2}const`$ is satisfied for most $`H`$-values and is independent of the magnitude of the dissipation coefficients and the order of the diffusion operator, provided the Reynolds number is sufficiently high such that $`H`$ is well conserved. The kinetic energy is found to decrease more rapidly than the magnetic one, in contrast to the behavior in 2D, in particular we find $`\mathrm{\Gamma }E`$. This proportionality leads to a simple energy decay law, $`dE/dtE^3`$, or $`Et^{0.5}`$. We also obtain the similarity law for the kinetic energy $`E^Vt^1`$. For the special case $`H=0`$ the energy decays more rapidly, $`Et^1`$, which agrees with previous simulations at lower Reynolds numbers. The transition to the finite-$`H`$ behavior occurs at relatively small values of $`H`$.

Results concerning the spatial scaling properties of 3D MHD turbulence will be published in a subsequent paper.

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 computer.