Cross-helically forced and decaying hydromagnetic turbulence

TL;DR

This paper investigates how strong cross helicity influences the evolution of kinetic and magnetic energy spectra in magnetohydrodynamic turbulence, revealing inverse energy transfer and exponential decay phenomena.

Contribution

It provides new insights into the spectral evolution of MHD turbulence under strong cross helicity, including inverse transfer and decay suppression mechanisms.

Findings

Weak inverse transfer of kinetic energy in forced turbulence.

Exponential decay of energy over a range of wavenumbers.

Dominance of kinetic over magnetic energy at small scales.

Abstract

We study the evolution of kinetic and magnetic energy spectra in magnetohydrodynamic flows in the presence of strong cross helicity. For forced turbulence, we find weak inverse transfer of kinetic energy toward the smallest wavenumber. This is plausibly explained by the finiteness of scale separation between the injection wavenumber and the smallest wavenumber of the domain, which here is a factor of 15. In the decaying case, there is a slight increase at the smallest wavenumber, which is probably explained by the dominance of kinetic energy over magnetic energy at the smallest wavenumbers. Within a range of wavenumbers covering almost an order of magnitude the decay is purely exponential, which is argued to be a consequence of a suppression of nonlinearity due to the presence of strong cross helicity.

Equations

$$\frac{\partial\mbox{\boldmath$A$}{}}{\partial t}=\mbox{\boldmath$U$}{}\times\bm{B}-\eta\mu_{0}\mbox{\boldmath$J$}{}+\mbox{\boldmath$f$}{}_{\rm M},$$

$$\frac{{\rm D}{}\mbox{\boldmath$U$}{}}{{\rm D}{}t}=-c_{\rm s}^{2}\mbox{\boldmath$\nabla$}{}\ln\rho+\frac{1}{\rho}\mbox{\boldmath$J$}{}\times\bm{B}+\frac{1}{\rho}\mbox{\boldmath$\nabla$}{}{\bm{\cdot}}(2\rho\nu\mbox{\boldmath${\sf S}$}{})+\mbox{\boldmath$f$}{}_{\rm K},$$

$$\frac{{\rm D}{}\ln\rho}{{\rm D}{}t}=-\mbox{\boldmath$\nabla$}{}{\bm{\cdot}}\mbox{\boldmath$U$}{}.$$

$$\displaystyle\mbox{\boldmath$f$}{}_{\rm M}\!\!$$

$$\displaystyle\mbox{\boldmath$f$}{}_{\rm K}\!\!$$

$$N_{\mathrm{M}}$$

$$N_{\mathrm{K}}$$

$$\mbox{\boldmath$f$}{}_{\bm{k}(t)}=\frac{\bm{k}(t)\times\mbox{\boldmath$e$}{}(t)}{\sqrt{\bm{k}(t)^{2}-\left[\bm{k}(t){\bm{\cdot}}\mbox{\boldmath$e$}{}(t)\right]^{2}}},$$

$$\text{mbox}{\mathrm{Pr}}_M=\nu /\eta ,\text{mbox}{\mathrm{Re}}_M=u_{\mathrm{rms}}/\eta k_{\mathrm{f}},\text{mbox}\mathrm{Lu}=v_{\mathrm{A}}/\eta k_{\mathrm{f}},$$

$$\sigma_{c}=\frac{2\langle\mbox{\boldmath$u$}{}\cdot\bm{B}\rangle}{q\langle\mbox{\boldmath$u$}{}^{2}\rangle+\langle\bm{B}^{2}\rangle/q}$$

$$k_{\mathrm{d}}(t)=({ϵ}_M/{\eta }^3{)}^{1/4},$$

$$k_{\mathrm{I}}^{-1}(t)=\int k^{-1}{E}_{\mathrm{M}}\mathrm{d}k/\int E_{\mathrm{M}}\mathrm{d}k.$$

Full text

††articletype: Original Article

\fundingInfo

University of Colorado. NSF Astronomy a Astrophysics Grants Program, 1615100.

1 Introduction

Conservation laws fundamentally affect the cascade properties of hydrodynamic and magnetohydrodynamic (MHD) turbulence. This can be seen both in forced and decaying turbulence, but the effects are often more dramatic in the decaying case. Decaying MHD turbulence is strongly affected by the presence of non-vanishing magnetic helicity [see, e.g., Tevzadze et al., 2012, Kahniashvili et al., 2013]. The study of decaying MHD turbulence has recently received increased attention in connection with the study of the decay of primordial magnetic fields during the early Universe, because the possibility of an inverse cascade leads to progressively larger length scales of the turbulent magnetic field and can reach kiloparsec scales at the present time [Brandenburg et al., 1996], if it was of sub-horizon scales at the time of the electroweak phase transition some $10^{-11}\,{\rm s}$ after the big bang; see Durrer & Neronov [2013] and Subramanian [2016] for recent reviews. The possibility of an inverse cascade is believed to be related to the conservation of magnetic helicity [Frisch et al., 1975]. However, since the work of Woltjer [1958], we know that there is another important conserved quantity, the cross helicity. The question arises whether finite cross helicity can have similar effects. In particular we wish to know whether cross helicity causes a slow-down of the decay of MHD turbulence and whether there is accelerated growth of the correlation length when cross helicity becomes important, as was found by Christensson et al. [2001], Banerjee & Jedamzik [2004], and Tevzadze et al. [2012], for example. In those cases, it was found that the fractional magnetic helicity, i.e., the magnetic helicity normalized by the magnetic energy and the correlation length, grows proportional to $t^{1/2}$, because magnetic helicity stays constant, the energy decreases like $t^{-1}$, and the correlation length increases like $t^{1/2}$. A qualitatively similar behavior has been found in the shell model work of Frick & Stepanov [2010], although in their case the growth of fractional helicity may have been a consequence of a true instability, whose origin is not yet understood, however.

As is well known, in homogeneous incompressible MHD turbulence the normalized cross helicity tends to grow in magnitude, which is associated with weakening of the nonlinearities in the MHD equations [e.g., Dobrowolny et al., 1980a, b, Grappin et al., 1982, 1983, Matthaeus et al., 1983, Pouquet et al., 1986, Matthaeus et al., 2008, Servidio et al., 2008]. The related question of a slow-down of the growth of cross helicity has already been partly addressed by Sur & Brandenburg [2009], who showed that for the non-helical Archontis flow [Archontis, 2000, Archontis et al., 2003], which is driven by a forcing that is proportional to $(\sin kz,\sin kx,\sin ky)$, the mean magnetic field grows at first exponentially and then develops a slow saturation behavior. However, this slow-down does not appear to depend on the microphysical values of kinematic viscosity or magnetic diffusivity.

For completeness, we note that cross helicity can also be produced in strongly stratified turbulence in the presence of magnetic fields parallel to the direction of gravity [Rüdiger et al., 2011], which also leads to the growth of magnetic fields at large length scales, which is suggestive of inverse transfer behavior [Brandenburg et al., 2014]. In that case, however, cross helicity develops more rapidly and it also disappears quickly, if the external magnetic field is removed. This case is inhomogeneous owing to the presence of gravitational stratification and will not be considered in the present paper.

2 Simulation setup

Our main objective is to compute the decay of MHD turbulence in the presence of cross helicity. As initial conditions we take the result of an earlier simulation where both velocity $u$ and magnetic field $\bm{B}$ were driven by the same forcing function, which ensures that finite cross helicity \mbox{\boldmathu}{}\cdot\bm{B} is injected into the system. We chose the injection wavenumber to be large enough so that there is a chance to see inverse transfer from the forcing wavenumber $k_{\rm f}$ to the smallest wavenumber of the domain, $k_{1}$. The ratio $k_{\rm f}/k_{1}$ is referred to as the scale separation ratio, and we choose $k_{\rm f}/k_{1}=15$ in all our simulations.

Our setup for the driven case is identical to that of Brandenburg & Rädler [2013], except that they included also rotation and used a smaller scale separation ratio of $k_{\rm f}/k_{1}=5$. As in earlier work [Kahniashvili et al., 2013, Tevzadze et al., 2012], we use an isothermal equation of state so that the pressure $p$ and the mass density $\rho$ are proportional to each other, i.e., $p=\rho c_{\rm s}^{2}$, where $c_{\rm s}={\rm const}{}$ is the isothermal sound speed. The magnetic field $\bm{B}$, the fluid velocity $U$, and the mass density $\rho$ obey

[TABLE]

[TABLE]

[TABLE]

Here, $A$ is the magnetic vector potential, so \mbox{\boldmath\nabla}{}\times\mbox{\boldmathA}{}=\bm{B} is the magnetic field, and \mbox{\boldmathJ}{}=\mbox{\boldmath\nabla}{}\times\bm{B}/\mu_{0} is the current density, where $\mu_{0}$ is the vacuum permeability, $\eta$ is the magnetic diffusivity, {\rm D}/{\rm D}{}t=\partial/\partial t+\mbox{\boldmathU}{}{\bm{\cdot}}\mbox{\boldmath\nabla}{} is the advective time derivative, {\sf S}_{ij}={\textstyle{1\over 2}}(U_{i,j}+U_{j,i})-{\textstyle{1\over 3}}\delta_{ij}\mbox{\boldmath\nabla}{}{\bm{\cdot}}\mbox{\boldmathU}{} are the components of the trace-less rate of strain tensor, $\nu$ the kinematic viscosity, and \mbox{\boldmathf}{}_{\rm M} and \mbox{\boldmathf}{}_{\rm K} define the magnetic and kinetic forcings, respectively. These will be specified below. We consider small Mach numbers, so compressibility effects are negligible. No DC magnetic field is imposed.

Equations (1)–(3) are solved numerically in a cubic domain of side length $L$ using periodic boundary conditions. Thus, $k_{1}=2\pi/L$ is the smallest possible (non-zero) wavenumber. During the first part, before studying the decay, we drive the system in a cross-helical fashion such that \mbox{\boldmathf}{}_{\rm K}=\mbox{\boldmath\nabla}{}\times\mbox{\boldmathf}{}_{\rm M} (ignoring a correction factor for units), with random functions that are $\delta$-correlated in time.

We approximate a forcing that is $\delta$-correlated in time by adding after each time step of size $\delta t$ the contributions \delta t\mbox{\boldmathf}{}_{\rm M} and \delta t\mbox{\boldmathf}{}_{\rm K} to $A$ and $U$, respectively, and change \mbox{\boldmathf}{}_{\rm M} and \mbox{\boldmathf}{}_{\rm K} randomly from one step to the next [Brandenburg, 2001]. Thus, we put

[TABLE]

where $N_{\rm M}$ and $N_{\rm K}$ are given by

[TABLE]

and ${\cal N}_{\rm M}$ and ${\cal N}_{\rm K}$ are dimensionless amplitudes, $\rho_{0}$ is the initial mass density, considered to be uniform, $k_{\rm f}$ the average forcing wavenumber, $\delta t$ the size of the time step, and \mbox{\boldmathf}{}_{\bm{k}} is given by

[TABLE]

where \mbox{\boldmathe}{}(t) is a unit vector which is in the same sense random as $\bm{k}(t)$ but not parallel to it [Haugen et al., 2004]. The wavevector $\bm{k}$ and the phase $\phi$ are random functions of time, i.e., $\bm{k}=\bm{k}(t)$ and $\phi=\phi(t)$, such that their values within a given time step are constant. Note that \mbox{\boldmath\nabla}{}{\bm{\cdot}}\mbox{\boldmathf}{}_{\rm M}=\mbox{\boldmath\nabla}{}{\bm{\cdot}}\mbox{\boldmathf}{}_{\rm K}=0. The wavevectors $\bm{k}$ are chosen such that their moduli $k=|\bm{k}|$ lie in a band of width $\delta k$ around a mean forcing wavenumber $k_{\rm f}$, that is, $k_{\rm f}-\delta k/2\leq k\leq k_{\rm f}+\delta k/2$, and we choose $\delta k=k_{1}$. In the limit of small time steps, which we approach in our calculations, the forcing may be considered as $\delta$-correlated.

When a statistically steady state is reached, we set $N_{\rm M}=N_{\rm K}=0$ and study in that way decaying turbulence. We describe the statistically stationary state of our simulations using the magnetic Prandtl number $\mbox{\rm Pr}_{M}$, the magnetic Reynolds number $\mbox{\rm Re}_{M}$, and the Lundquist number Lu,

[TABLE]

with $u_{\rm rms}$ and $v_{\rm A}=b_{\rm rms}/\sqrt{\mu_{0}\rho_{0}}$ being defined using averages over the full computational volume. While $\mbox{\rm Pr}_{M}$ is an input parameter, $\mbox{\rm Re}_{M}$ and Lu are used as diagnostics. We analyze our results in terms of kinetic and magnetic energy spectra, $E_{\rm K}(k,t)$ and $E_{\rm M}(k,t)$ that are normalized such that \int E_{\rm K}(k,t)\,{\rm d}{}k=\langle\mbox{\boldmathU}{}^{2}/2\rangle\equiv{\cal E}_{\rm K} and $\int E_{\rm M}(k,t)\,{\rm d}{}k=\langle\bm{B}^{2}/2\mu_{0}\rangle\equiv{\cal E}_{\rm M}$ are the mean kinetic and magnetic energy densities, respectively.

For our numerical simulations we use the Pencil Code111https://github.com/pencil-code, DOI:10.5281/zenodo.2315093, which is a high-order public domain code for solving partial differential equations, including the hydromagnetic equations given above. It uses sixth order finite differences in space and the third order 2N-RK3 low storage Runge–Kutta time stepping scheme of Williamson [1980].

3 Results

3.1 Driven turbulence

We begin by presenting results of simulations with a resolution of $512^{3}$ meshpoints. In Figure 1 we show kinetic and magnetic energy spectra, in both uncompensated and compensated forms. Note that at the smallest wavenumbers ($k/k_{\rm f}=1/15$) the spectral kinetic energy exceeds the magnetic spectral energy and its spectrum is approximately flat. This is reminiscent of simulations that exhibit an inverse cascade (or at least inverse transfer) of kinetic energy owing to what is known as the anisotropic kinetic alpha (AKA) effect of Frisch et al. [1987]; see also Brandenburg & von Rekowski [2001] for simulations at larger Reynolds numbers. Here, however, no AKA effect is expected to occur. Thus, the slight uprise of $E_{\rm K}(k,t)$ at low $k$ is a new result that occurs now with the addition of cross helicity. It is therefore tempting to associate it with the conservation of cross helicity.

3.2 Decaying turbulence

We use a snapshot from the statistically steady phase of the forced simulation with the forcing being turned off. During the subsequent decay, the normalized cross helicity,

[TABLE]

tracks the evolution of $v_{\rm A}/u_{\rm rms}$, both of which attain peak values reaching 0.996 shortly after turning off the forcing and then decay to about 0.8 after some 500 sound travel times. In Equation (10), we have introduced $q=\sqrt{\mu_{0}\rho}$ for dimensional reasons [Zhang & Brandenburg, 2018].

In Figure 2 we show magnetic and kinetic energy spectra in logarithmically spaced time intervals during the decay phase at $t/c_{\rm s}k_{1}=5$, 10, 20, 50, 100, 200, and 500 sound travel times, corresponding to $t/v_{\rm A0}k_{\rm f}=1.23$, …, 123 Alfvén times, where $v_{\rm A0}$ is the initial Alfvén speed. Note that the spectral kinetic and magnetic energies decrease at large wavenumbers, while at small wavenumbers the kinetic energy decreases and the magnetic energy increases so that the two quantities approach each other at large scales (although they are still far from being equal). We recall that MHD absolute equilibrium studies predict that the Alfvén ratio $r_{\rm A}={\cal E}_{\rm K}/{\cal E}_{\rm M}\leq 1$ [Stribling & Matthaeus, 1990], although the relationship of these results to the particular initially forced situation considered here, and decaying situations in general [Stribling & Matthaeus, 1991], warrants further investigation. In our case, instead, $r_{\rm A}$ increases from unity to about 1.5 during the course of the simulation.

The magnetic Prandtl number is unity and the Reynolds number based on the smallest wavenumber in the domain decreases from about 5000 to about 200 after about 500 sound travel times. The decay is neither algebraic nor exponential; see Figure 3. This is explained by the fact that the decay of magnetic and kinetic energies is actually composed of a continuous sequence of uncoupled modes, each decaying exponentially with their own resistive decay rate, $\mu k^{2}$, where $\mu=\nu+\eta$ is the sum of kinematic viscosity and magnetic diffusivity.

In Figure 4 we plot the time dependence of the magnetic dissipation wavenumber,

[TABLE]

where $\epsilon_{\it M}=2\eta\int E_{\rm M}k^{2}\,{\rm d}{}k$ is the instantaneous mean energy dissipation rate. We also plot the evolution of the integral wavenumber,

[TABLE]

Both quantities are evaluated from the actual spectra.

In Figure 5 we show the decay of magnetic and kinetic energies for a few selected wavenumbers. The instantaneous decay rates $\lambda(t)$ are shown in the lower panel. We determine $\lambda(t,k)$ as an average of the instantaneous decay rates for a suitable time interval when the decay is indeed exponential. The resulting decay rates $\overline{\lambda}$ are shown in Figure 6 and compared with the visco-resistive decay rate.

Note that, within a certain wavenumber interval, the decay does indeed follow a mode-by-mode exponential decay, so there is no coupling between different wavenumbers, except for $k/k_{1}>20$, where the decay is slower than what is expected based on a purely visco-resistive decay. This must be a consequence of a depletion of the \mbox{\boldmathu}{}\times\bm{B} nonlinearity when \mbox{\boldmathu}{}\cdot\bm{B} is maximized.

4 Conclusions

Our work has shown that in forced MHD turbulence with strong cross helicity, kinetic energy displays a slight uprise at the smallest wavenumbers. This is reminiscent of the inverse transfer found in turbulence with a non-Galilean invariant forcing, giving rise to an AKA effect. Here, however, the forcing is Galilean invariant because our forcing function is $\delta$ correlated in time, so it has no memory of the previous forcing time step and therefore no proper motion of the forcing field can be defined. We are therefore tempted to associate this small uprise of spectral kinetic energy with the presence of cross helicity.

In the decaying case, we find that the kinetic energy decays at all wavenumbers, i.e., there is no evidence for inverse transfer or inverse cascade behavior. This raises doubts about our tentative conclusion regarding the forced case, where the small uprise might also be just the result of the finiteness of the domain and would disappear at larger scale separation or for larger domains. We have not studied this here, but refer instead to a similar observation by Yokoi & Brandenburg [2016], who simulated rotating forced hydrodynamic turbulence in the presence of a profile of kinetic helicity and found inverse transfer, but only when the domain was not too large.

Unlike $E_{\rm K}(k,t)$, which shows a decay at all $k$, $E_{\rm M}(k,t)$ does actually show a slight increase at the smallest wavenumbers. This is likely a consequence of the presence of finite kinetic energy at large-scale scales and is reminiscent of the inverse transfer found in nonhelical MHD turbulence found in the magnetically dominated case [Brandenburg et al., 2015], where the magnetic energy has a $k^{4}$ subinertial range, which is steeper than the $k^{2}$ subinertial range found for kinetic energy. Again, then, the kinetic energy exceeds the magnetic energy at the smallest wavenumbers.

In all the cases presented here, the relation with ‘ordinary’ (i.e., low cross helicity) turbulence is unclear, because the presence of strong cross helicity implies a suppression of nonlinearity [e.g., Dobrowolny et al., 1980a]. As a result, the decay law is, within a certain wavenumber interval, purely exponential and not like a power law, as in ordinary turbulence. As is well-known, when the nonlinear terms are negligible each Fourier mode will undergo exponential decay at the appropriate wavenumber-dependent decay rate associated with the linear dissipation terms.

When we started our work, we regarded an initial integral wavenumber of $15\,k_{1}$ as large enough. This may need to be reconsidered in future work, because in related studies of decaying hydromagnetic turbulence, an initial scale separation of 60 has meanwhile been found to be more adequate [Brandenburg et al., 2015]. Also, a resolution of $512^{3}$ mesh points may not be enough for studying the possibility of inverse transfer. Most importantly, perhaps, it would be useful to work with more controlled initial conditions that have well determined sub-inertial and inertial range spectra. Such studies have been done in the presence of magnetic helicity [Brandenburg & Kahniashvili, 2017, Brandenburg et al., 2017], but not yet in the presence of cross helicity. Also, our way of initializing cross helicity by setting \mbox{\boldmathA}{}\propto\mbox{\boldmathU}{} may be rather special. Another possibility is to drive cross helicity through the application of gravity and a magnetic field aligned with it, as done in the earlier work of Rüdiger et al. [2011] and Brandenburg et al. [2014]. This would lead to a stratified system in which the occurrence of inverse transfer can directly be associated with the formation of spots that have been found to form as a generic result of the negative effective magnetic pressure instability for vertical magnetic fields see Fig. 1 of Brandenburg et al. [2013] and Brandenburg et al. [2016] for a review.

Acknowledgements

This research was supported in part by the Astronomy and Astrophysics Grants Program of the National Science Foundation (grant 1615100), and the University of Colorado through its support of the George Ellery Hale visiting faculty appointment. We acknowledge the allocation of computing resources provided by the Swedish National Allocations Committee at the Center for Parallel Computers at the Royal Institute of Technology in Stockholm.