Inverse cascade in decaying 3D magnetohydrodynamic turbulence

We perform direct numerical simulations of three-dimensional freely decaying magnetohydrodynamic (MHD) turbulence. For helical magnetic fields an inverse cascade effect is observed in which power is transfered from smaller scales to larger scales. The magnetic field reaches a scaling regime with self-similar evolution, and power law behavior at high wavenumbers. We also find power law decay in the magnetic and kinematic energies, and power law growth in the characteristic length scale of the magnetic field.


I. INTRODUCTION
Within cosmology, astrophysics or geophysics one often needs to deal with electrically conducting plasmas at high kinematic and magnetic Reynolds numbers where magnetic fields are dynamically important.Indeed, much of the turbulence in the interstellar medium is magnetohydrodynamic in nature.
Hydromagnetic turbulence has been explored extensively in connection with the generation of large scale magnetic fields in astrophysical bodies such as planets, stars, accretion discs and galaxies through dynamo theories.Non-driven, freely decaying turbulence may also be of interest in connection with both the physics of the interstellar medium and cosmology.Our interest was inspired by the cosmology of primordial magnetic fields, which is sometimes considered as a possible source for providing the seed for the galactic dynamo [1].
There have been various related works on decaying MHD turbulence, by authors interested in different contexts [2][3][4][5][6][7].Most directly comparable to our work, Biskamp and Müller [6] studied the energy decay in incompressible 3D magnetohydrodynamic turbulence in numerical simulations at relatively high Reynolds number, and in a companion letter [7] studied the scaling properties of the energy power spectrum.
We are here especially interested in the inverse cascade of magnetic helicity whereby magnetic energy is transferred from small to large scale fluctuations.This is important for a primordial magnetic field to reach a large enough scale with sufficient amplitude to be relevant for seeding the galactic dynamo [8].
It should be noted that due to the conformal invariance of MHD in the radiation era the MHD equations in an expanding universe can be converted into the relativistic MHD equations in flat spacetime by an appropriate scaling of the variables and by using conformal time [9].The equations of [9] differ slightly from the ordinary non-relativistic MHD equations.However, in order to facilitate comparison with earlier work, we use the non-relativistic equations.
We perform 3D simulations both with and without magnetic helicity, starting from statistically homogeneous and isotropic random initial conditions, with power spectra suggested by cosmological applications.We find a strong inverse cascade in the helical case, with equivocal evidence for a weak inverse cascade when only helicity fluctuations are present.In the helical case we also find a self-similar power spectrum with an approximately k −2.5 behavior at high k.We present energy decay laws which are comparable to those found in the incompressible case by Biskamp and Müller [6], and in the compressible case by Mac Low et al. [5].

II. THE MODEL
We consider the equations for an isothermal compressible gas with a magnetic field, which is governed by the momentum equation, the continuity equation, and the induction equation, written here in the form where B = ∇ × A is the magnetic field in terms of the magnetic vector potential A, u is the velocity, J is the current density, ρ is the density, µ is the dynamical viscosity, and η is the magnetic diffusivity.
The code for solving these equations [10] uses a variable third order Runge-Kutta timestep and sixth order explicit centered derivatives in space.All our runs are performed on a 120 3 grid, and we use periodic boundary conditions, which means that the average plasma density ρ 0 = ρ 0 is conserved during runs.Here ρ 0 is the value of the initially uniform density, and the brackets denote volume average.
We adopt nondimensional quantities by measuring u in units of c, where c is the speed of light, k in units of k 1 , where k 1 is the smallest wavenumber in the box, which has a size of L BOX = 2π, density in units of ρ 0 = 1, and B is measured in units of √ µ 0 ρ 0 c, where µ 0 is the magnetic permeability.This is equivalent to putting c = k 1 = ρ 0 = µ 0 = 1.In the following we will refer to the mean kinematic viscosity ν which we define as ν ≡ µ/ρ 0 .The sound speed c s takes the value c s = 1/ √ 3, as appropriate for a relativistic fluid.With c = 1, the unit of time is such that the light crossing time of the box is 2π.
Our equations are similar to those for the relativistic gas in the early universe using scaled variables and conformal time for non-relativistic bulk velocities [9].We expect our results to change little using the true relativistic equations, as our advection velocity is at most only mildly relativistic, and this only at the beginning of the simulation.

III. ON THE ROLE OF THE INVERSE CASCADE
The magnetic helicity H M is given by and characterizes the linkage between magnetic field lines.H M is conserved in the absence of ohmic dissipation, although it is still possible to have local, small scale helicity fluctuations.Helicity plays an important role in dynamo theory [11,12], where turbulence is driven.In many astrophysical and cosmological situations the magnetic Reynolds number Re M is very large.We define the magnetic Reynolds number as Re M = Lv/η, where L and v are the typical length scale and velocity of the system under consideration and η is the resistivity.The magnetic Reynolds number is a measure of the relative importance of flux freezing versus resistive diffusion.In a cosmological context this number can be extraordinarily large: causality imposes the weak limit L ≤ ct and relativity demands v < c.With conductivities relevant to the era when the electroweak phase transition took place [13], one can in principle obtain a magnetic Reynolds number of about 10 16 .This is often taken to mean that the magnetic field is frozen into the plasma, and the scale length of the field increases only with the expansion of the Universe.
However, this simple picture does not necessarily give a full description of the dynamics because the MHD equations, especially at high Reynolds numbers where nonlinear terms are important, exhibit turbulent behavior, which can lead to a redistribution of magnetic energy over different length scales [9].Energy in a turbulent magnetic field can undergo an inverse cascade and be transferred from high frequency modes to low frequency modes, increasing the overall comoving correlation length [11].This process is due to the nonlinear terms giving rise to interactions between many different length scales.
We will take the initial primordial power spectrum as given and address the question of how such a primordial spectrum evolves as a consequence of the nonlinear equations of motion.

IV. INITIAL CONDITIONS
Since one of the aims of the present work is to investigate the role of magnetic helicity in the inverse cascade we describe how the initial conditions for our simulations were set up.We chose our initial condition by setting up magnetic fluctuations with an initial power spectrum Fourier space (and averaged over shells of constant k = |k|), for low values of the wavenumber k, using a exponential cutoff k c .
[The shell-averaged power spectrum, P M (k), is not to be confused with the shell-integrated energy spectrum, The magnetic field fluctuations are drawn from a Gaussian random field distribution fully determined by its power spectrum in Fourier space according to the following procedure.For each grid point we use the corresponding wavenumber to select an amplitude from a Gaussian distribution centered on zero and with the width where k = |k|.We then transform the field back into real space to obtain the field at each grid point.This is done independently for each field component.
There is a requirement in cosmology that n ≥ 2, which is set by causality demanding that the correlation function of the magnetic field vanishes at large distances, and the fact that the magnetic field is divergence-free [14].In the simulations presented we chose the slope of the power spectrum to be n = 2.We also chose k c = 30, unless specified otherwise, which gives a power spectrum peaked at a relatively large value of k.Biskamp and Müller [6,7] started with a spectrum peaked at k c = 4, which may account for the different slope in the late-time power spectrum which we observe (see Section V A).
Our velocity power spectrum was chosen in a similar way, but with n = 0 corresponding to white noise at large scales (there is no requirement for incompressibility in the early Universe).The initial magnetic energy was taken equal to the kinetic energy, and had the value 5 × 10 −3 in all runs, as the primordial field is thought likely to be weak.
In order to introduce a non-zero average magnetic helicity into the system it is useful to represent the vector potential in terms of its projection onto an orthogonal basis formed by ê+ , ê− and k.The two basis vectors ê+ and ê− can be chosen to be the unit vectors for circular polarization, right-handed and left-handed respectively.That is ê± = ê1 ± iê 2 where ê1 and ê2 are unit vectors orthogonal to each other and to k.They are given by ê1 Note that since where s = ±1, this corresponds to an expansion of the magnetic vector potential into helical modes.
Using these basis vectors it is easily seen that the magnetic energy spectrum is where the amplitude of the magnetic field is given by and the expression for the magnetic helicity spectrum where The function H M (k) is a sensitive measure of the correlation between the vector potential and the magnetic field.H M (k) may, of course, be positive in one part of Fourier space and negative in another part.It is, however, bounded in magnitude by the inequality A field which saturates the above inequality is maximally helical.
The amplitudes A ± k can be chosen independently, provided A * ± −k = A ± k , which is just the condition that the vector potential be real.Therefore it is possible to adjust the amplitudes |A + k | and |A − k | freely and in so doing obtaining a magnetic field with arbitrary magnetic helicity.With our method we are able to put statistically random but maximally helical fields in our initial conditions.In our runs with initial helicity we take H M = H max .
Because we evolve our dynamical fields on a discrete lattice we have to be careful when using derivative operations in Fourier space.In general, the wave vector, which is an eigenvalue of the derivative operator, needs to be replaced by some function k eff (k), which is an eigenvalue of the discrete derivative operator on the lattice.In our case we have, for the sixth order explicit centered derivative In order to be consistent with the scheme used in the simulation, we use k eff (k) when calculating the initial condition in Fourier space.

V. RESULTS
In all runs the mean kinematic viscosity ν and the resistivity η were chosen to be equal with values between ν = η = 5 × 10 −4 − 5 × 10 −5 .In our simulations we typically obtain Reynolds numbers of the order of 100 − 200.The Reynolds numbers in our simulations are evaluated using the magnetic Taylor microscale which we calculate here as the ratio of the rms magnetic field and the rms current density, L T = 2πB rms /J rms .The 2π factor is here included so that L T represents the typical wave length (and not the inverse wavenumber) of structures in the current density.

A. Spectral evolution
The inverse magnetic cascade for decaying MHD turbulence is best visualized in terms of magnetic energy spectra E M (k) because information on nonlinear interaction between different scales is contained in E M (k).In Fig. 1 we show a run with initial magnetic helicity.In Fig. 1 we see evidence for a dual energy transfer both toward higher and lower wavenumbers.The inverse cascade is characterized by the transfer of energy from small scale structures in the magnetic field to larger ones.In Fig. 1 this behavior is clearly seen as indicated by the rise in the energy spectrum at small wavenumbers.Some energy is also being transported to smaller scales where the spectrum is decaying due to diffusive effects.We also note that at wavenumbers above the peak k p (t) the spectrum develops a power law shape.This power law has approximately a k −2.5 slope.This differs from the approximately k −5/3 law found by Müller and Biskamp [7].We suggest that this is due to finite size effects, which affect the spectrum if the initial scale separation between k p and the smallest wavenumber in the box (k = 1) is insufficient, and if the flow is strongly helical so that its spectrum is governed by inverse cascading.In order to check this we have performed a run with larger initial length scale, k c = 5.In this case the magnetic energy spectrum develop into an approximate k −5/3 law at late times.However, this occurs only after the peak of the spectrum has left the simulation box, i.e. after finite size effects have begun to play a role.
To check if the magnetic field evolution is self-similar one can make the following ansatz for the energy spectrum Here ξ is the characteristic length scale of the magnetic field, taken to be the magnetic Taylor microscale defined above, and q is a parameter whose value is some real number.We call g M (kξ) the magnetic scaling function.FIG. 2. The magnetic scaling function gM(kξ) described in the text, equation ( 13), versus kξ.The straight lines indicate the power laws ∝ (kξ) 4.0 and ∝ (kξ) −2.5 respectively.
In Fig. 2 we have plotted ξ(t) q E M (k, t) versus the scaled variable kξ(t).The value of the parameter q in this run is q = 0.7.It is seen that for each different value of time t, the data collapses onto a single curve given by the scaling function g M (kξ), demonstrating the self-similarity of the magnetic field evolution.
We also performed runs in which the magnetic helicity was zero, in the statistical sense.Magnetic helicity was present due to fluctuations, but was of very small amplitude.In these runs no significant inverse cascade was observed.Fig. 3 shows the energy spectrum for such a run with only small magnetic helicity fluctuations present in the initial conditions.It is seen that only a weak inverse cascade is present at the lowest wavenumbers, much smaller than in the helical case.However, that it seem to be present at all is interesting as the effect could become more pronounced for higher Reynolds numbers.It is possible that this effect is due to the magnetic helicity fluctuations even though they were small.One simulation was performed with identically zero initial magnetic helicity fluctuations.In this case random fluctuations develop rapidly and no differences between the two cases, were observed.

B. Energy decay
In Fig. 4 we show the time evolution of the magnetic energy E M (t) and the kinetic energy E K (t) for a run with initial helicity and a k 4 initial energy spectrum slope.FIG. 4. Time evolution of the magnetic energy EM(t) and the kinetic energy EK(t) in the case where there is initial magnetic helicity.ν = η = 5 × 10 −5 .The straight lines indicate the power laws ∝ t −0.7 and ∝ t −1.1 respectively.
It is seen that the asymptotic decay rate for E M (t) is approximately t −0.7 .The Reynolds number for this run was around Re ∼ 200 at late times.In another run with Re ∼ 100 the decay rate was seen to be t −0.8 , so there seems to be a dependence of the decay rate of the magnetic field on the Reynolds number and perhaps the resulting slope of the spectrum.
The kinetic energy also decays with a power law behavior at late times.In the case of runs with initial helicity the kinetic energy E K (t) decays with a different, faster rate than E M (t).The asymptotic decay rate is close to t −1.1 .In runs without initial helicity the decay rates of E M (t) and E K (t) are approximately the same, close to t −1.1 .
In our runs with E K = E M initially, the kinetic energy spectrum shows no evidence of an inverse cascade at any scale.However, when the initial velocity distribution is zero the kinetic spectrum grows on all scales initially and in the low wavenumber region the energy continues to grow even after the high wavenumber modes start to decay.

C. Coherence length evolution
During the course of the simulations the initially small scale structures gain in size.A convenient length scale is the magnetic Taylor microscale L T , which was defined above.This length scale is mostly characteristic of the small scales but even they grow during the course of the simulations.
In Fig. 5 we show the evolution of L T for a run with initial helicity.The asymptotic behavior of the length scale is seen to grow approximately as L T ∼ t 0.5 .
In runs with non-helical initial conditions the growth of the magnetic Taylor microscale is slower than in the case of helical initial conditions.In this case the magnetic Taylor microscale grows approximately as L T ∼ t 0.4 .
The discussion so far has mainly been concerned with the evolution of causally generated magnetic fields using an initial k 4 slope in the magnetic energy spectrum.Now we briefly comment on the other cases we have looked at.For a white noise initial spectrum E M (k) ∼ k 2 , the evolution is qualitatively and quantitatively similar to the causal case.For helical fields we observe an inverse cascade, while for non-helical fields a much smaller inverse cascade is present only for the lowest modes.

VI. DISCUSSION
Our simulations show the decay rate of magnetic energy for compressible turbulence being sensitive to the initial helicity of the magnetic field configuration.A similar result was found in [6] in the case of incompressible turbulence.The fact that magnetic helicity is conserved (except for resistive changes), and the magnetic energy decays slower for helical fields, is connected with the observed inverse cascade in which magnetic energy is transported toward larger scales because of nonlinear dynamics.
The decay of kinetic energy does not seem to depend on the initial helicity and its decay rate E K (t) ∼ t −1.1 is consistent with the earlier work of [5,6].Note that in the helical case we observe the kinetic energy decaying more rapidly than the magnetic one; this behavior was also found in [6].
While these results are not directly applicable to the evolution of primordial magnetic fields in the early universe, they do suggest that nonlinear magnetohydrodynamical effects may play an important role in this case.
In any case, it is interesting to compare our results with the work of other authors interested in the decay properties of cosmological magnetic fields [15][16][17][18].Ideal MHD has a scale invariance which leads to the scaling law [15,17] where ψ is an unknown function, related to g M .Assuming it is peaked somewhere, and h < 0, the characteristic scale of the field goes as L(t) ∼ t 1/(1−h) .It is also often assumed that ψ(0) exists and is non-zero: thus h is determined by the initial power spectrum.Hence for a magnetic power spectrum of index n, h = −(n+3)/2 and This law can also be recovered by assuming that the characteristic time scale for the decay of turbulence on a scale l is the eddy turnover time τ l = l/v l , where v l ∼ l −(n+3) is the velocity averaged on a scale l [16].If the characteristic scale of the field is that scale which is just decaying, then τ L ∼ t, and we again find Eq.( 15).One should note that these arguments ignore helicity conservation.We recall that our non-helical runs had n = 2 for the magnetic power spectrum and n = 0 for the velocity power spectrum.The observed growth law for the magnetic Taylor microscale, t 0.4 , is not consistent with the predicted power law for n = 2, although it does square with the growth law for n = 0, and it is possible that the growth in the magnetic field length scale is being controlled by the velocity field.Simulations at higher Reynolds numbers seem required to resolve this issue.
One would expect on integrating the helicity power spectrum that H M ∼ L I E M , where L I is the integral scale.We would expect that L I ∼ L T , and hence, if magnetic helicity is conserved, However, magnetic helicity is not conserved exactly: we observe a decrease in H M by a factor of about 2 in a run with viscosity ν = 5 × 10 −5 .Indeed, with L T ∼ t 0.5 we find a somewhat steeper relation: E M ∼ t −0.7 ∼ L −1.4 T .Finally, it is interesting to note that Son's numerical simulations of decaying turbulence [16], performed in the Eddy-Damped Quasi-Normal Markovian (EDQNM) approximation, show some evidence of a power law developing at high k, the slope being close to k −2.5 , although there was no net helicity present, and no inverse cascade.Furthermore, Field and Carroll [18], again using the EDQNM approximation, found that there were selfsimilar solutions with E M ∼ t −2/3 ∼ L −1 T .

VII. CONCLUSIONS
We have shown that for an isothermal and compressible magnetized turbulent fluid, when undergoing a process of free decay, a substantial inverse cascade is present for helical magnetic field configurations, which transfer energy from smaller scale magnetic fluctuation to larger scale ones.For non-helical magnetic fields only a weak inverse cascade was observed on the largest scales.
The energy spectrum of the magnetic field shows evidence for a self-similar evolution with a development of a power law of roughly k −2.5 beyond the peak.Decay laws for both the kinematic and magnetic energy were found.The kinetic energy decay was approximately t −1.1 for both helical and non-helical magnetic fields.The decay of the magnetic field energy was found to be strongly dependent on the the initial helicity, decaying roughly as t −0.7 and t −1.1 for helical and non-helical initial conditions respectively.For the helical case, the magnetic energy decay rate showed a dependence on the Reynolds number, with a slower decay rate for larger Reynolds numbers.
We also observed power law behavior in the characteristic length scale of the magnetic field, defined as the Taylor microscale L T .In the helical case L T ∼ t 0.5 , whereas for non-helical fields the growth was somewhat slower, L T ∼ t 0.4 , and we ascribe the faster growth rate to the presence of the inverse cascade in the helical case.