# Generation of chiral asymmetry via helical magnetic fields

Jennifer Schober<EMAIL_ADDRESS>Laboratoire d’Astrophysique, EPFL,
CH-1290 Sauverny, Switzerland  Tomohiro Fujita
<EMAIL_ADDRESS>Department of Physics, Kyoto University,
Kyoto, 606-8502, Japan  Ruth Durrer<EMAIL_ADDRESS>Départment de
Physique Théorique and Center for Astroparticle Physics, Université de Genève,
Quai E. Ansermet 24, CH-1211 Genève 4, Switzerland

###### Abstract

It is well known that helical magnetic fields undergo a so-called inverse
cascade by which their correlation length grows due to the conservation of
magnetic helicity in classical ideal magnetohydrodynamics (MHD). At high
energies above approximately $10$ MeV, however, classical MHD is necessarily
extended to chiral MHD and then the conserved quantity is
$\langle\mathcal{H}\rangle+2\langle\mu_{5}\rangle/\lambda$ with
$\langle\mathcal{H}\rangle$ being the mean magnetic helicity and
$\langle\mu_{5}\rangle$ being the mean chiral chemical potential of charged
fermions. Here, $\lambda$ is a (phenomenological) chiral feedback parameter.
In this paper, we study the evolution of the chiral MHD system with the
initial condition of nonzero $\langle\mathcal{H}\rangle$ and vanishing
$\mu_{5}$. We present analytic derivations for the time evolution of
$\langle\mathcal{H}\rangle$ and $\langle\mu_{5}\rangle$ that we compare to a
series of laminar and turbulent three-dimensional direct numerical
simulations. We find that the late-time evolution of
$\langle\mathcal{H}\rangle$ depends on the magnetic and kinetic Reynolds
numbers ${\rm Re}_{{}_{\mathrm{M}}}$ and ${\rm Re}_{{}_{\mathrm{K}}}$. For a
high ${\rm Re}_{{}_{\mathrm{M}}}$ and ${\rm Re}_{{}_{\mathrm{K}}}$ where
turbulence occurs, $\langle\mathcal{H}\rangle$ eventually evolves in the same
way as in classical ideal MHD where the inverse correlation length of the
helical magnetic field scales with time $t$ as $k_{\mathrm{p}}\propto
t^{-2/3}$. For a low Reynolds numbers where the velocity field is negligible,
the scaling is changed to $k_{\mathrm{p}}\propto
t^{-1/2}\mathrm{ln}\left(t/t_{\mathrm{log}}\right)$. After being rapidly
generated, $\langle\mu_{5}\rangle$ always decays together with
$k_{\mathrm{p}}$, i.e. $\langle\mu_{5}\rangle\approx k_{\mathrm{p}}$, with a
time evolution that depends on whether the system is in the limit of low or
high Reynolds numbers.

Magnetohydrodynamics – turbulence – relativistic processes – magnetic fields –
early universe

##  I Introduction

Natural systems can be modeled as fluids when their macroscopic spatial
extension is much larger than the typical mean free path of particle
collisions. This often applies in geophysics and astrophysics including
planets, stars, the interstellar medium, and galaxies. Hydrodynamics is even
applicable in cosmology, in particular, when modeling the plasma of the early
Universe and the early stages of cosmic structure formation. Often, natural
fluids are highly turbulent which is quantified by large hydrodynamic Reynolds
numbers ${\rm Re}_{{}_{\mathrm{K}}}$. The latter measures the ratio of
advection and diffusion effects and is defined as ${\rm
Re}_{{}_{\mathrm{K}}}\equiv u_{\mathrm{rms}}/(k_{\mathrm{f}}\nu)$ where
$u_{\mathrm{rms}}$ is the rms velocity, $k_{\mathrm{f}}$ is the forcing wave
number, and $\nu$ is the viscosity. Indeed, ${\rm Re}_{{}_{\mathrm{K}}}\gg 1$
in many astrophysical fluids resulting from an efficient turbulent driving
that has its origin, for instance in galaxies, in supernova explosions (Mac
Low and Klessen, 2004) and/or accretion flows (Klessen and Hennebelle, 2010).
Moreover, the local Universe appears to be permeated with magnetic fields on
all length scales accessible to observations: They span from planets Stevenson
(2003) and stars Donati and Landstreet (2009), including our Sun Solanki et
al. (2006), to the interstellar medium Crutcher (2012), galaxies Beck (2012),
up to galaxy clusters Govoni and Feretti (2004) and possibly cosmic voids
Neronov and Vovk (2010). The most established mechanisms to amplify and
maintain magnetic fields in stars and galaxies are magnetohydrodynamical (MHD)
dynamos that convert kinetic energy into magnetic energy (Brandenburg and
Subramanian, 2005). In general, MHD turbulence describes the dynamics of many
astrophysical and cosmological flows.

Without an energy input, turbulence and magnetic fields decay freely which can
be of interest in various astrophysical applications. One such example of non-
driven MHD turbulence is the evolution of magnetic fields in the very early
Universe before the epoch of recombination. Scenarios of primordial
magnetogenesis include specific models for inflation (e.g. Turner and Widrow,
1988) and the cosmological phase transitions (e.g. Sigl et al., 1997), but the
subsequent evolution of primordial magnetic fields is governed by the laws of
decaying MHD turbulence 111It is worth noting that in a flat expanding
Universe, the equations of relativistic MHD simplify to the nonrelativistic
MHD equations if all quantities are replaced by co-moving quantities and
instead of physical time conformal time is used Brandenburg et al. (1996).,
(see e.g. Durrer and Neronov, 2013, for a review). Whether primordial magnetic
fields can survive until they might e.g. serve as seed fields for galactic
dynamos or remain as significant relics in present-day cosmic voids, depends
on the prospect of transferring magnetic energy from small to large spatial
scales. Such a transfer of magnetic energy is known as inverse cascade and is
well studied within MHD turbulence (Pouquet et al., 1976; Christensson et al.,
2001).

The scaling laws of the inverse cascade depend crucially on the magnetic
helicity $\int_{V}\mathbf{A}\cdot\mathbf{B}~{}\mathrm{d}V$ which is a
topological property of the magnetic field. Here, the volume integral is taken
over the product of the vector potential $\mathbf{A}$ and the magnetic field
$\mathbf{B}=\nabla\times\mathbf{A}$. Since
$\int_{V}\mathbf{A}\cdot\mathbf{B}~{}\mathrm{d}V$ is a conserved quantity in
ideal MHD, a decay of the magnetic field necessarily results into an increase
of its correlation length $\ell$. Using a phenomenological approach, similar
to the Kolmogorov theory of nonmagnetized turbulence, magnetic helicity
conservation implies scalings with time $t$ of $\ell\propto t^{2/3}$ and
$B\propto t^{-1/3}$ (Biskamp and Müller, 1999; Müller and Biskamp, 2000). It
has been demonstrated in three-dimensional numerical simulations that inverse
transfer of energy occurs even in nonhelical decaying MHD turbulence, however
with less efficiency (Brandenburg et al., 2015; Reppin and Banerjee, 2017).
The scaling laws of decaying MHD turbulence, estimated from phenomenological
arguments or extracted from numerical simulations, have been applied, for
instance, to the evolution of primordial magnetic fields from their generation
until recombination (Kahniashvili et al., 2013; Brandenburg et al., 2017a).

The assumption of conserved magnetic helicity, however, breaks down when
charged fermions can be considered as being massless. This is the case at high
energies that are reached in the very early Universe. In fact, above
approximately $10$ MeV (Boyarsky et al., 2012), MHD necessarily needs to be
extended to chiral MHD where only the sum of magnetic helicity and fermionic
chirality is conserved. Fermionic chirality is the divergence of the chiral
current $j_{5}$ which is not conserved due to the chiral anomaly, a pure
quantum effect with macroscopic consequences. It can be quantified by the
chiral chemical potential $\mu_{5}\equiv\mu_{\mathrm{L}}-\mu_{\mathrm{R}}$
with $\mu_{\mathrm{R}}$ and $\mu_{\mathrm{L}}$ are the chemical potentials of
right- and left-handed fermions, respectively. A non-zero $\mu_{5}$ in the
presence of a magnetic field, leads to the chiral magnetic effect (CME)
Vilenkin (1980) which is a macroscopic quantum effect within the standard
model of particle physics and implies an additional electric current along a
magnetic field (Son and Surówka, 2009). The CME leads to a magnetic field
instability (Joyce and Shaposhnikov, 1997) which is the subject of many
studies (Fröhlich and Pedrini, 2000, 2002; Semikoz and Sokoloff, 2004; Semikoz
et al., 2009; Boyarsky et al., 2012; Semikoz et al., 2012; Tashiro et al.,
2012; Dvornikov and Semikoz, 2012, 2013, 2014, 2017; Manuel and Torres-Rincon,
2015; Gorbar et al., 2016; Pavlović et al., 2016; Pavlović et al., 2017;
Rogachevskii et al., 2017; Figueroa et al., 2019). Recently, the nonlinear
dynamics of a chiral plasma was studied also in direct numerical simulations
(DNS) with a focus on chiral dynamos (Brandenburg et al., 2017b; Schober et
al., 2018; Masada et al., 2018). The energies necessary for chiral effects are
reached in the early Universe and in protoneutron stars, but also in heavy ion
collisions (Kharzeev et al., 2008). Furthermore, chiral MHD is relevant for
modeling the dynamics of electronic solid-state materials like Weyl semimetals
(Galitski et al., 2018).

The extension of MHD to chiral MHD raises the following questions: How is the
inverse cascade affected by the new degree of freedom, the chiral chemcial
potential? How much chiral asymmetry can be generated from an initial helical
magnetic field? This scenario has been explored by Hirono et al. (2015) who
considered a plasma composed of charged fermions with initial magnetic
helicity and vanishing chiral asymmetry. Within their assumption of negligible
velocity fields, they have identified a three-stage evolution: First, the
magnetic helicity is transferred to fermionic chirality due to the
conservation law of chiral MHD. Second, the total helicity is dominated by
fermionic chirality which eventually leads to a CME-assisted inverse cascade
of magnetic helicity. And third, at late times, Hirono et al. (2015) report a
self-similar evolution of $\mu_{5}$ and the peak of the magnetic energy
spectrum proportional to $t^{-1/2}$. This self-similar evolution during the
decay of a large $\mu_{5}$ has been observed in lattice simulations Figueroa
et al. (2019); Mace et al. (2019). A remaining open question is, however, how
such an evolution of chiral MHD is modified in presence of turbulence, where
there can be a strong coupling between the magnetic field and the velocity
field. Understanding decaying chiral helical MHD turbulence and its
differences to the classical MHD scenario is the goal of the present study.

To this end, we investigate the evolution of the magnetic field $\bm{B}$, the
velocity field $\bm{U}$ and the chiral chemical potential $\mu_{5}$ in both
the laminar and the turbulent regime. Our initial conditions are a vanishing
chiral asymmetry and velocity field, $\bm{U}=\mu_{5}=0$ and a maximally
helical magnetic field. These initial conditions are realized, for instance,
in various inflationary magnetogenesis models where a pseudo-scalar field
generates magnetic fields in a parity violating manner (Field and Carroll,
2000; Anber and Sorbo, 2006; Durrer et al., 2011; Caprini and Sorbo, 2014;
Fujita et al., 2015; Adshead et al., 2016; Fujita and Durrer, 2019), but does
not introduce a chemical potential for fermions. Note that if the fermion
masses are negligible, not only the helical magnetic field but also the chiral
asymmetry can be generated during inflation such that the net helicity plus
chirality is conserved precisely due to the chiral anomaly Domcke and Mukaida
(2018). This alternative initial condition is beyond the scope of the present
study but a target of our future work.

The paper is structured as follows. In Section II we briefly review the
inverse cascade in classical MHD and the system of equations in chiral MHD and
introduce our numerical methods. In Section III we present an analytical
derivation of the self-similar inverse cascade in chiral MHD with a vanishing
velocity field and confirm the validity of our analytical results with three-
dimensional numerical simulations. The transition from a system with vanishing
velocity field to a regime where turbulence is driven efficiently via the
Lorentz force exerted by the helical magnetic field, is presented in Section
IV. For the limit of large Reynolds numbers, we use a phenomenological
approach to find solutions for the evolution of $\mu_{5}$. The analytical
solutions are compared to results from turbulent DNS. We draw our conclusions
in Section V.

##  II Theoretical background and methods

###  II.1 Review of the classical inverse cascade

The dynamics of magnetized fluids in the one-fluid magnetohydrodynamical limit
is described by the following set of equations:

$\displaystyle\frac{\partial\bm{B}}{\partial t}$ $\displaystyle=$
$\displaystyle\nabla\times\left[{\bm{U}}\times{\bm{B}}-\eta\,\left(\bm{\nabla}\times{\bm{B}}\right)\right],$
(1) $\displaystyle\rho\frac{D\bm{U}}{Dt}$ $\displaystyle=$
$\displaystyle(\bm{\nabla}\times{\bm{B}})\times\bm{B}-\bm{\nabla}p+\bm{\nabla}{\bm{\cdot}}(2\nu\rho\mbox{\boldmath${\sf
S}$}{}),$ (2) $\displaystyle\frac{D\rho}{Dt}$ $\displaystyle=$
$\displaystyle-\rho\,\bm{\nabla}\cdot\bm{U}.$ (3)

Here, $\bm{B}$ is the magnetic field, $t$ is time, $\bm{U}$ is the velocity
field, and $\rho$ is the mass density. Furthermore, $p$ is the hydrodynamic
pressure, ${\sf
S}_{ij}=1/2(U_{i,j}+U_{j,i})-1/3~{}\delta_{ij}{\bm{\nabla}}{\bm{\cdot}}\bm{U}$
are the components of the tracefree strain tensor ${\sf S}$, where commas
denote partial spatial derivatives, and $D/Dt=\partial/\partial
t+\bm{U}\cdot\bm{\nabla}$ is the advective derivative. The Ohmic resistivity
is denoted by $\eta$ and $\nu$ is the viscosity. The set of equations is
closed by an isothermal equation of state, meaning that the pressure is
related to the density via $p=c_{\rm s}^{2}\rho$, where $c_{\rm s}$ is the
sound speed.

An important role for the evolution of magnetic fields is played by magnetic
helicity which is defined as $\int_{V}\mathcal{H}~{}\mathrm{d}V$ with
$\mathcal{H}\equiv\bm{A}\cdot\bm{B}$. The integral is taken over a periodic
volume $V$ or over an unbounded volume with the fields falling off
sufficiently rapidly at spatial infinity so that a boundary term can be
neglected. In these cases, magnetic helicity is gauge invariant. Its evolution
equation can be derived by multiplying Faraday’s law with its uncurled version
for the vector potential and yields

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\int_{V}\mathcal{H}~{}\mathrm{d}V=-2\eta\int_{V}\bm{J}\cdot\bm{B}~{}\mathrm{d}V.$
(4)

A remarkable consequence of magnetic helicity conservation at $\eta\rightarrow
0$ (faster than the current helicity $\int_{V}\bm{J}\cdot\bm{B}~{}\mathrm{d}V$
may possibly diverge) is the inverse cascade of energy for a fully helical
magnetic field (Frisch et al., 1975; Pouquet et al., 1976).

The highly nonlinear evolution of helical decaying MHD turbulence has been
studied intensely with DNS. For incompressible 3D magnetohydrodynamic
turbulence at relatively high ${\rm Re}_{{}_{\mathrm{M}}}$, the energy decay
(Biskamp and Müller, 1999) as well as scaling relations of the energy power
spectrum have been analyzed (Müller and Biskamp, 2000). The role of magnetic
helicity in the inverse cascade was investigated by Christensson et al.
(2001). With their DNS, Christensson et al. (2001) found evidence for a self-
similar evolution of magnetic energy spectrum with a development of a power
law of roughly $k^{-2.5}$ beyond the peak and analyzed decay laws for both the
kinematic and magnetic energy. The scaling relations of a helical magnetic
field can, in the limit of high Reynolds numbers, be derived by using a
Kolmogorov-type phenomenological approach (see e.g. Biskamp, 2003). In
particular, the magnetic energy evolves as $\langle\bm{B}^{2}\rangle/2\propto
t^{-2/3}$ and the correlation length of the magnetic field as $\xi\propto
t^{2/3}$ which has been confirmed by DNS (e.g. Kahniashvili et al., 2013). An
inverse transfer of magnetic energy has also been found for nonhelical
magnetic fields, however, it is less efficient than in the fully helical case
(Kahniashvili et al., 2013; Zrake, 2014; Berera and Linkmann, 2014;
Brandenburg et al., 2015; Reppin and Banerjee, 2017).

Table 1: Summary of all runs presented in this paper. The reference runs for laminar (R1 and R1mhd) and turbulent (R8 and R8mhd) simulations which are presented in detail in Figures 1 and 10, respectively, are highlighted by bold font. The amplitude of the initial magnetic power spectrum is exactly the same for all runs R1–R7. Runs R7 and R7mhd have, however, a larger initial rms magnetic field strength $B_{0}$ which is due to the higher resolution, meaning the larger number of modes available in R7 (the exponentially suppressed tail extends to higher wave numbers). For runs R8, R8b, and R8mhd, a larger initial amplitude has been set.  |  |  | Input parameters: |  |  |  |  | Measured parameters: |  |
---|---|---|---|---|---|---|---|---|---|---
Name | MHD | resolution | $\dfrac{B_{\mathrm{rms},0}}{10^{-2}}$ | $k_{\mathrm{p},0}$ | $\mu_{5,0}$ | $\dfrac{B_{\mathrm{rms},0}}{\eta}$ | $\dfrac{\lambda B_{\mathrm{rms},0}^{2}}{k_{\mathrm{p},0}}$ | ${\rm Re}_{{}_{\mathrm{M}}}^{{\mathrm{min}(\mathcal{H})}}$ | ${\rm Re}_{{}_{\mathrm{M}}}^{{k_{\mathrm{p}}=1}}$ | ${\rm Re}_{{}_{\mathrm{M}}}^{\mathrm{max}}$
R1a | chiral | $320^{3}$ | $1.153$ | $85$ | 0 | $11.53$ | $1.662$ | $7.2\times 10^{-6}$ | $1.9\times 10^{-5}$ | $9.6\times 10^{-4}$
R1 | chiral | $\bf 320^{3}$ | $\bf 1.153$ | $\bf 85$ | $\bf 0$ | $\bf 11.53$ | $\bf 16.618$ | $\bf 2.7\times 10^{-5}$ | $\bf 1.4\times 10^{-1}$ | $\bf 1.6\times 10^{-1}$
R1b | chiral | $320^{3}$ | $1.153$ | $85$ | 0 | $11.53$ | $166.176$ | $1.6\times 10^{-3}$ | $2.3\times 10^{-1}$ | $2.3\times 10^{-1}$
R1mhd | classic | $\bf 320^{3}$ | $\bf 1.153$ | $\bf 85$ | $\bf-$ | $\bf 11.53$ | $\bf-$ | $\bf 2.4\times 10^{-6}$ | $\bf 2.9\times 10^{-5}$ | $\bf 1.0\times 10^{-3}$
R2 | chiral | $320^{3}$ | $1.153$ | $85$ | 0 | $23.06$ | $16.618$ | $1.2\times 10^{-4}$ | $-$ | $2.3\times 10^{-1}$
R3 | chiral | $320^{3}$ | $1.153$ | $85$ | 0 | $115.3$ | $16.618$ | $2.6\times 10^{-3}$ | $3.4\times 10^{0}$ | $3.6\times 10^{0}$
R4 | chiral | $320^{3}$ | $1.153$ | $85$ | 0 | $230.6$ | $16.618$ | $9.8\times 10^{-3}$ | $7.1\times 10^{0}$ | $7.1\times 10^{0}$
R5 | chiral | $320^{3}$ | $1.153$ | $85$ | 0 | $576.5$ | $16.618$ | $8.1\times 10^{-2}$ | $-$ | $6.6\times 10^{0}$
R6 | chiral | $320^{3}$ | $1.153$ | $85$ | 0 | $1153.0$ | $16.618$ | $1.6\times 10^{0}$ | $-$ | $9.9\times 10^{0}$
R7 | chiral | $512^{3}$ | $1.400$ | $85$ | 0 | $2800.0$ | $24.5$ | $1.4\times 10^{1}$ | $-$ | $2.5\times 10^{1}$
R7mhd | classic | $512^{3}$ | $1.400$ | $85$ | $-$ | $2800.0$ | $-$ | $8.3\times 10^{0}$ | $-$ | $1.2\times 10^{1}$
R8 | chiral | $\bf 512^{3}$ | $\bf 4.667$ | $\bf 85$ | $\bf 0$ | $\bf 9333.6$ | $\bf 24.501$ | $\bf 1.1\times 10^{2}$ | $\bf 3.2\times 10^{2}$ | $\bf 3.2\times 10^{2}$
R8b | chiral | $512^{3}$ | $4.667$ | $85$ | $0$ | $9333.6$ | $2450.140$ | $5.5\times 10^{1}$ | $-$ | $7.8\times 10^{2}$
R8mhd | classic | $\bf 512^{3}$ | $\bf 4.667$ | $\bf 85$ | $\bf-$ | $\bf 9333.6$ | $\bf-$ | $\bf 2.6\times 10^{2}$ | $\bf 2.7\times 10^{2}$ | $\bf 2.7\times 10^{2}$

###  II.2 Chiral MHD equations

At high energies, MHD necessarily needs to be generalized to chiral MHD in
which the chiral asymmetry appears as a new degree of freedom. Here, an
asymmetry between the number densities of left- and right-handed fermions
gives rise to the CME that results in an electric current proportional to the
magnetic field and a chiral chemical potential

$\displaystyle\mu_{5}^{\mathrm{phys}}\equiv\mu_{{}_{\rm L}}-\mu_{{}_{\rm R}},$
(5)

where $\mu_{{}_{\rm L}}$ and $\mu_{{}_{\rm R}}$ are the chemical potentials of
left- and right-handed fermions, respectively. In the following, we will
replace $\mu_{5}^{\mathrm{phys}}$ by a normalized version,
$\mu_{5}\equiv(4\alpha_{\rm em}/\hbar c)\mu_{5}^{\mathrm{phys}}$, that has the
same units as a wave number.

Chiral asymmetry is coupled to magnetic helicity and significantly modifies
the phenomenology of the plasma. The set of chiral MHD equations is given by
(Rogachevskii et al., 2017):

$\displaystyle\frac{\partial\bm{B}}{\partial t}$ $\displaystyle=$
$\displaystyle\bm{\nabla}\times\left[{\bm{U}}\times{\bm{B}}-\eta\,\left(\bm{\nabla}\times{\bm{B}}-\mu_{5}{\bm{B}}\right)\right],$
(6) $\displaystyle\rho\frac{D\bm{U}}{Dt}$ $\displaystyle=$
$\displaystyle(\bm{\nabla}\times{\bm{B}})\times\bm{B}-\bm{\nabla}p+\bm{\nabla}{\bm{\cdot}}(2\nu\rho\mbox{\boldmath${\sf
S}$}{}),$ (7) $\displaystyle\frac{D\rho}{Dt}$ $\displaystyle=$
$\displaystyle-\rho\,\bm{\nabla}\cdot\bm{U},$ (8)
$\displaystyle\frac{D\mu_{5}}{Dt}$ $\displaystyle=$
$\displaystyle\mathscr{D}_{5}\,\Delta\mu_{5}+\lambda\,\eta\,\left[{\bm{B}}{\bm{\cdot}}(\bm{\nabla}\times{\bm{B}})-\mu_{5}{\bm{B}}^{2}\right].$
(9)

Here, the chiral vortical effect, the chiral separation effect, and chirality
flipping are neglected. The latter is well justified at sufficiently high
temperatures while the former might not be a very good approximation in the
case of high Reynolds number where large vortical velocities can be generated.
The new equation (9) includes a diffusion term with the diffusion constant
$\mathscr{D}_{5}$ that is only relevant when $\mu_{5}$ is strongly varying in
space. In this study, the evolution of $\mu_{5}$ will be mostly affected by
the electromagnetic field via the second term on the right-hand side of
Equation (9). The strength of the feedback is controlled by the chiral
feedback parameter

$\displaystyle\lambda=3\hbar c\left(\frac{8\alpha_{\rm em}}{k_{\rm
B}T}\right)^{2},$ (10)

which is valid for $k_{\rm B}T\gg\max(|\mu_{L}|,|\mu_{R}|$) (Boyarsky et al.,
2015). In the following, $\lambda$ will be considered constant, yet one should
keep in mind that it scales with temperature $T^{-2}$ (see, e.g. Boyarsky et
al., 2012).

The system of Equations (6)–(9) implies a conservation law:

$\displaystyle\frac{\partial}{\partial
t}\left(\frac{\lambda}{2}\mathcal{H}+\mu_{5}\right)$
$\displaystyle+\bm{\nabla}{\bm{\cdot}}\left(\frac{\lambda}{2}\left({\bm{\bm{E}}}\times{\bm{A}}+\bm{B}\,\Phi\right)-\mathscr{D}_{5}\bm{\nabla}\mu_{5}\right)=0,$
(11)

where
${\bm{\bm{E}}}=-c^{-1}\,[{\bm{\bm{U}}}{\bm{\times}}{\bm{B}}+\eta\,(\mu_{5}{\bm{B}}-{\bm{\nabla}}{\bm{\times}}{\bm{B}})]$
is the electric field and $\Phi$ is the electrostatic potential. Hence, the
total chirality $\langle\mathcal{H}\rangle+2\langle\mu_{5}\rangle/\lambda$ is
a conserved quantity, where $\langle\mu_{5}\rangle$ is the mean value of the
chiral chemical potential and $\langle\mathcal{H}\rangle\equiv
V^{-1}\int\bm{A}\cdot\bm{B}\,dV$ is the mean magnetic helicity density in the
volume $V$.

###  II.3 Numerical methods

To go beyond the limitations of analytical calculations, we use the Pencil
Code 222http://pencil-code.nordita.org/ for solving Equations (1)–(3) for
classical MHD and Equations (6)–(9) for chiral MHD, respectively. The system
of equations is solved in a three-dimensional periodic domain of size
$L^{3}=(2\pi)^{3}$ via a third-order accurate time-stepping method of
Williamson (1980) and sixth-order explicit finite differences in space
(Brandenburg and Dobler, 2002; Brandenburg, 2003). The time step is specified
as the Courant time step, however, for our reference runs we initially use
very small manually set time steps to resolve the very early time evolution in
more detail. After the initial phase, the automatic Courant time step is used
in the reference runs. The resolution is varied between $320^{3}$ and
$512^{3}$. The smallest wave number covered in the numerical domain is
$k_{1}=2\pi/L=1$ which we use as normalization of length scales. All
velocities are normalized to the sound speed $c_{\rm s}=1$ and the mean fluid
density to $\overline{\rho}=1$. Time is normalized by the diffusion time
$t_{\eta}=\eta^{-1}$, where $\eta$ is the Ohmic resistivity.

In this study, all runs are initialized with vanishing chiral chemical
potential $\mu_{5}$ and a strong helical random magnetic field 333For the
setting the initial condition of the helical random magnetic fields, we use
the routine power_randomphase_hel in the Pencil Code.. In practice, the
magnetic field is set up via the vector potential $\bm{A}(\boldsymbol{x})$
which is constructed from a random and $\delta$-correlated three-dimensional
vector field in real space. The magnetic field is calculated from the Fourier
transform of $\bm{A}(\boldsymbol{x})$ via
$\bm{B}(\boldsymbol{k})=i\boldsymbol{k}\times\bm{A}(\boldsymbol{k})$. Then the
magnetic field is scaled by functions of $k$ to adjust the shape of the
magnetic energy spectrum $E_{\rm M}(k)=2\pi B^{2}(k)k^{2}$ for which we use a
slope proportional to $k^{4}$ between $1$ and the initial wave number of the
energy-carrying eddies $k_{\mathrm{p},0}=85$, i.e. the initial peak of the
magnetic energy spectrum. For $k>k_{\mathrm{p},0}$ in all runs, the spectrum
is suppressed. Multiplication by the operator
$P_{ij}(\boldsymbol{k})-i\sigma_{\mathrm{M}}\hat{k}_{l}$, where
$P_{ij}=\delta_{ij}-\hat{k_{i}}\hat{k_{j}}$ is the projection operator,
ensures a fully helical magnetic field for $\sigma_{\mathrm{M}}=\pm 1$.
Finally, the energy spectrum is normalized such that
$\langle\bm{B}^{2}\rangle/2=\bm{B}_{\mathrm{rms}}^{2}/2=\int_{1}^{k_{\mathrm{max}}}E_{\rm
M}(k)\,\mathrm{d}k$, where the integration is performed over the entire
numerical domain, i.e. from $k=1$ up to the maximally resolved wave number
$k_{\mathrm{max}}$. We note that these initial conditions of the magnetic
field are similar to the ones used in previous studies of decaying MHD
turbulence, like in Brandenburg et al. (2017a).

No external forcing is applied to drive turbulence in our simulations, i.e.,
the velocity field is purely driven via the Lorentz force that is exerted on
the flow through the magnetic field. The transition to a turbulent plasma
occurs when the magnetic and kinetic Reynolds numbers, ${\rm
Re}_{{}_{\mathrm{M}}}=u_{\mathrm{rms}}/(k_{\mathrm{f}}\eta)$ and ${\rm
Re}_{{}_{\mathrm{K}}}=u_{\mathrm{rms}}/(k_{\mathrm{f}}\nu)$, respectively,
become much larger than unity. Here, $u_{\mathrm{rms}}$ is the rms velocity
and $k_{\mathrm{f}}$ is the wave number on which kinetic energy is injected in
the system. For magnetically driven turbulence, $k_{\mathrm{f}}$ corresponds
to the inverse correlation length of the magnetic field and we will use
$k_{\mathrm{f}}=k_{\mathrm{p}}(t)$. Viscosity $\nu$ and Ohmic resistivity
$\eta$ are implemented explicitly in the code. To explore systems with
different Reynolds numbers, the values of $\nu$ and $\eta$ are systematically
changed, while their ratio, i.e. the magnetic Prandtl number ${\rm
Pr}_{{}_{\mathrm{M}}}$, is set to unity for all of the simulations. We note,
that the choice of ${\rm Pr}_{{}_{\mathrm{M}}}=1$ does not reflect the
situation in most astrophysical applications. However, ${\rm
Pr}_{{}_{\mathrm{M}}}\ll 1$ or ${\rm Pr}_{{}_{\mathrm{M}}}\gg 1$ are
notoriously difficult to treat in DNS, since that requires a large separation
of scales. In what follows, we therefore only mention the magnetic Reynolds
number which in our settings equals the kinetic one.

An overview of the input parameters and characteristic numbers of all runs
discussed in this work is presented in Table 1.

##  III Inverse cascade in chiral MHD with a vanishing velocity field

In this section, we discuss the evolution of a decaying helical magnetic field
for simulations where the velocity field can be neglected throughout the
entire simulation time. We note, however, that Equation (7) is nevertheless
included in the DNS.

Figure 1: Comparing a run of a decaying helical magnetic field in classical
MHD (Left, run R1mhd) with chiral MHD (Right, run R1). The velocity field in
these simulations is negligible. Top: Time evolution of the mean magnetic
helicity $\langle\mathcal{H}\rangle$ and the wave number on which the magnetic
energy spectrum has its maximum $k_{\mathrm{p}}$. The time $t$ is normalized
to the resistive time $t_{\eta}$. For chiral MHD, we also show the mean chiral
chemical potential $\langle\mu_{5}\rangle$ over the chiral feedback parameter
$\lambda$ as well as the conserved quantity
$\langle\mathcal{H}\rangle+2\langle\mu_{5}\rangle/\lambda$. The solid vertical
lines indicate the time when $k_{\mathrm{p}}$ has reached the minimum wave
number possible in the numerical domain, $k=1$, and the dotted vertical lines
indicates the time at which the minimum of magnetic helicity is reached.
Middle: $k_{\mathrm{p}}$ and $\langle\mu_{5}\rangle$ normalized by their
theoretically expected scaling in the self-similar evolution phase (see
Section III.1). Bottom: Evolution of the magnetic energy spectrum. The thick
dotted black lines show the initial spectra and the thick solid black lines
show the final spectra of the simulations.

###  III.1 Analytical solutions in the limit of dynamically insignificant
velocity fields

Provided that helical magnetic fields are nonzero but the chiral chemical
potential is vanishing $\mu_{5}(t_{0})\equiv\mu_{5,0}=0$ at the initial time,
the second term in the right-hand side in Equation (9) sources $\mu_{5}$.
$|\mu_{5}|$ may grow until it reaches the maximum value at which the second
term and the third term cancel each other,

$|\mu_{5}(t)|\simeq\left|\frac{\bm{B}\cdot(\nabla\times\bm{B})}{B^{2}}\right|\simeq
k_{p}(t),$ (12)

where we ignore the diffusion term $\mathscr{D}_{5}\Delta\mu_{5}$ in Equation
(9). Due to the conservation law
$\langle\mathcal{H}\rangle+2\langle\mu_{5}\rangle/\lambda=\mathrm{const.}$,
however, $\mu_{5}$ cannot exceed the initial value of the conserved quantity,

$|\mu_{5}(t)|\lesssim\lambda|\mathcal{H}(t_{0})|\simeq\lambda
B^{2}_{0}/k_{\mathrm{p},0},$ (13)

where $B_{0}$ is the strength of the initial magnetic field with maximum
energy on the wave number $k_{\mathrm{p},0}$. Therefore, the chiral chemical
potential generated from magnetic helicity is given by

$\displaystyle|\mu_{5}(t)|\simeq\mathrm{min}\left(k_{\mathrm{p}}(t),\lambda
B_{0}^{2}/k_{\mathrm{p},0}\right).$ (14)

Since $k_{\mathrm{p}}(t)$ decreases in time due to the diffusion of magnetic
energy on large $k$ (small scales), even if $\mu_{5}$ reaches $\lambda
B_{0}^{2}/k_{\mathrm{p},0}$ at early times, it starts decaying when
$k_{\mathrm{p}}(t)$ becomes smaller than $\mu_{5}$ because of the third term
in Equation (9). Thus $\mu_{5}(t)$ eventually follows $k_{\mathrm{p}}(t)$.
Depending on the initial conditions, two different scenarios can be identified
for decaying chiral MHD, as long as the velocity field can be neglected.

If $k_{\mathrm{p},0}>\lambda B_{0}^{2}/k_{\mathrm{p},0}$ the system evolves in
three phases:
(i) Production of $\langle\mu_{5}\rangle$ at the expense of magnetic helicity
$\langle\mathcal{H}\rangle$, with an efficiency depending on $\lambda$ up to
$\lambda B_{0}^{2}/k_{\mathrm{p},0}$.
(ii) Once $k_{\mathrm{p}}(t)$ has decayed to
$\langle\mu_{5}\rangle\simeq\lambda B_{0}^{2}/k_{\mathrm{p},0}$, a chiral
dynamo, the so-called $v_{\mu}^{2}$ dynamo as discussed in Rogachevskii et al.
(2017), leads to an exponentially fast restoration of
$\langle\mathcal{H}\rangle$.
(iii) A self-similar inverse cascade sets in with
$k_{\mathrm{p}}\approx\langle\mu_{5}\rangle$, during which
$\langle\mu_{5}\rangle$ is converted into $\langle\mathcal{H}\rangle$
according to the conservation law.
In the other case with $k_{\mathrm{p},0}<\lambda B_{0}^{2}/k_{\mathrm{p},0}$,
the growth of $\langle\mu_{5}\rangle$ stops when it becomes comparable to
$k_{\mathrm{p}}(t)$ during the phase (i) and the system immediately enters the
phase (iii) by skipping (ii). In this case, the magnetic helicity always
dominates the conserved quantity and the evolution of magnetic fields is not
significantly altered by $\langle\mu_{5}\rangle$.

As we show in Appendix A that during this self-similar evolution in the phase
(iii), $\mu_{5}$ and $k_{\mathrm{p}}$ evolve as

$\displaystyle|\mu_{5}(t)|\approx
k_{\mathrm{p}}(t)\approx\left[\frac{3+n}{4\eta
t}~{}\ln\left(\frac{t}{t_{\mathrm{log}}}\right)\right]^{1/2},$ (15)

where $n$ denotes the slope of the initial magnetic helicity spectrum
$\langle\mathcal{H}(t_{0})\rangle_{k}\propto k^{n}$. Our DNS are initiated
with $n=4$. We found that the logarithmic correction time is roughly given by
$t_{\mathrm{log}}=(2k_{\mathrm{p},0}^{2})^{-1}$ in our simulations, which is
written as $t_{\mathrm{log}}=(2k_{p,0}^{2}/k_{1}^{2})^{-1}t_{\eta}$ when using
explicit units. This result up to the logarithmic correction term has also
been found in Hirono et al. (2015).

###  III.2 Comparison of classical MHD with a three-phase chiral MHD scenario
in DNS

The reference run for a three phase scenario of a decaying magnetic field in
chiral MHD is R1. In Figure 1, R1 is compared to a classical MHD analog
(R1mhd). The parameters and initial conditions in both runs are the same but
in R1 the plasma evolves according to Equations (6)–(9) and in R1mhd according
to Equations (1)–(3). R1 is presented in the right panels of Figure 1 and
R1mhd is presented in the left panels.

The time evolution of the mean magnetic helicity, $\langle\mathcal{H}\rangle$,
and the wave number on which the magnetic energy spectrum has its maximum,
$k_{\mathrm{p}}$, are presented in the top row of Figure 1. In the classical
nonideal MHD case the magnetic helicity decreases by approximately eight
orders of magnitude during one resistive time, i.e. until $t\approx 1$.
Resistivity acts on small spatial scales, e.g. large wave numbers $k$. This
leads to a decrease of the magnetic energy on large $k$ and therefore a move
of the peak scale of the magnetic energy spectrum, $k_{\mathrm{p}}$, to
smaller $k$. Note, that $k_{\mathrm{p}}$ has discrete values only, leading to
steps in its time evolution that become more evident at late times when
$k_{\mathrm{p}}$ approaches $1$. Since the velocity field is negligible during
the entire run, there is no inverse transfer of magnetic energy, as can be
seen in the evolution of the magnetic energy spectrum, see Figure 1c. The time
evolution of $k_{\mathrm{p}}$, normalized to the theoretically predicted value
for chiral MHD given in Equation (15), is presented in the middle row of
Figure 1. For R1mhd, Equation (15) is not valid and therefore the orange
dashed line in Figure 1b moves away from $1$ with increasing time.

The time evolution of $\langle\mathcal{H}\rangle$ in R1 is significantly
different from the one in classical MHD, as can be see in Figure 1d. First,
$\langle\mathcal{H}\rangle$ decreases by roughly two orders of magnitude. At
the same time, a mean chiral chemical potential $\langle\mu_{5}\rangle$ is
generated, such that the sum
$\langle\mathcal{H}\rangle+2\langle\mu_{5}\rangle/\lambda$ is conserved during
the entire run. The value of $k_{\mathrm{p}}$ decreases in time in R1, but not
as quickly as in R1mhd. The three phases described in Section III.1 can be
clearly distinguished in Figure 1d: Phase (i) during which
$\langle\mathcal{H}\rangle$ decreases and which ends at $t\approx 4\times
10^{-2}$ is followed by a phase of dynamo amplification, phase (ii). For
$t\gtrsim 4\times 10^{-1}$, $\langle\mu_{5}\rangle$ and $k_{\mathrm{p}}$
evolve in a self-similar way, what was defined as phase (iii). During this
phase, the evolution of $\langle\mu_{5}\rangle$ and $k_{\mathrm{p}}$ is
reasonably well described by Equation (15), as can be seen in Figure 1e. For
comparison with the scaling of $k_{\mathrm{p}}\propto t^{-1/2}$, we have added
the black dashed line in Figure 1e from which the simulation data clearly
deviates in phase (iii). The time evolution of the magnetic energy spectrum in
run R1 (right bottom panel) is very different from the one in R1mhd (Figure
1c). The main difference occurs at late times, where in R1mhd, the magnetic
energy first grows on $k\approx 5$ and then moves to smaller wave numbers in
an CME-assisted inverse cascade.

###  III.3 Dependence on the chiral feedback parameter $\lambda$

Figure 2: The mean chiral chemical potential, $\langle\mu_{5}\rangle$, as a
function of time for runs R1a, R1, and R1b which differ only in their value of
$\lambda$ (see Table 1). The solid vertical lines indicate the time when
$k_{\mathrm{p}}$ has reached the minimum wave number possible in the numerical
domain, $k=1$, and the dotted vertical lines indicates the time at which the
minimum of magnetic helicity is reached. The horizontal gray dotted line
indicates the threshold for $\mu_{5}$ above which a dynamo instability occurs
in the numerical box.  Figure 3: The mean magnetic helicity,
$\langle\mathcal{H}\rangle$, as a function of time for runs R1mhd, R1a, R1,
and R1b (see Table 1). The horizontal gray dotted line indicates the initial
value of $\langle\mathcal{H}\rangle$ which is equal for all runs presented in
this figure. The vertical lines indicate the same characteristic times as in
Figure 2.  Figure 4: The ratio of the peak scale of the magnetic energy
spectrum $k_{\mathrm{p}}$ measured from the simulation data over the
theoretical prediction (see Equation 15) as a function of time for runs R1,
R1a, R1b, and R1mhd (see Table 1). The vertical lines indicate the same
characteristic times as in Figure 2.

The generation of a chiral asymmetry and the subsequent evolution of the
plasma depends strongly on the chiral feedback parameter $\lambda$ or, in
dimensionless units, $\lambda B_{0}^{2}/k_{\mathrm{p},0}$ as compared to the
value of $k_{\mathrm{p},0}$. As discussed in Section III.1, for $\lambda
B_{0}^{2}/k_{\mathrm{p},0}\ll k_{\mathrm{p},0}$, we expect a three-phase
evolution, while for $\lambda B_{0}^{2}/k_{\mathrm{p},0}\gg k_{\mathrm{p},0}$,
our models suggest a 2-phase scenario. In this section, we present a
comparison of R1 with a run with a smaller value of $\lambda$, R1a, and a run
with a larger value of $\lambda$, R1b.

The time evolution of $\mu_{5}$ in runs R1a, R1, and R1b is presented in
Figure 2. The largest maximum value of $\mu_{5}\approx k_{\mathrm{p},0}=85$ is
reached in run R1b, for which $k_{\mathrm{p},0}=85<\lambda
B_{0}^{2}/k_{\mathrm{p},0}\approx 166.2$. This run, almost instantaneously,
enters the self-similar phase where both $\mu_{5}$ and $k_{\mathrm{p}}$ decay
proportional to $t^{-1/2}\left(\mathrm{ln}(t/t_{\mathrm{log}})\right)^{1/2}$.
This scaling is indicated as a grey line in Figure 2. During the entire run
time of R1b, magnetic helicity is sourced by $\mu_{5}$ and, therefore, a CME-
assisted inverse cascade of magnetic energy occurs.

The situation is different in the reference run R1, where
$k_{\mathrm{p},0}=85>\lambda B_{0}^{2}/k_{\mathrm{p},0}\approx 16.6$. In R1, a
value of $\mu_{5}\approx 8$ is generated quickly and stays constant up to
$t\approx 0.3$. At that time the value of $\langle\mathcal{H}\rangle$ is up to
two orders of magnitude less than its initial value; see the time evolution of
$\langle\mathcal{H}\rangle$ in Figure 3. Via a chiral dynamo
$\langle\mathcal{H}\rangle$ grows exponentially in time around $t\approx 0.2$.
Subsequently, R1 enters the self-similar evolution phase with an CME-assisted
inverse cascade.

Run R1a has the lowest value of $\lambda$ where $k_{\mathrm{p},0}=85>\lambda
B_{0}^{2}/k_{\mathrm{p},0}\approx 1.7$. For these parameters, a maximum value
of $\mu_{5}\approx 0.8$ is generated. For this run, unlike in our reference
runs, e.g. R1, we have applied an automatic time stepping method during the
entire simulation time. As a result, the value of $\mu_{5}$ increases from $0$
to $\approx 0.8$ within the first timestep, leaving the $\mu_{5}$ generation
phase unresolved. With the maximum $\mu_{5}$ being less than $1$, there can be
no chiral dynamo instability within the simulation domain. Therefore, after
its initial generation, $\mu_{5}$ remains approximately constant throughout
the entire simulation time. Simultaneously, $\langle\mathcal{H}\rangle$
decreases in time as can be seen in Figure 3. For comparison, also the MHD run
R1mhd is presented Figure 3. Here, the value of $\langle\mathcal{H}\rangle$
decreases constantly at a rate that is only slightly larger than for R1a. At
$t\approx 0.45$, the peak of the magnetic energy spectrum in R1a reaches the
minimum wave number in the simulation domain, $k_{\mathrm{p}}=1$. A dynamo
instability for $\mu_{5}=0.8$ would occur at $k_{5}=\mu_{5}/2=0.4$. Since this
is outside of the numerical domain, a dynamo and a subsequent CME-assisted
inverse cascade is not seen in R1a. We stress, however, that this is purely
caused by the finiteness of the numerical domain. For infinite systems, a
three-phase scenario is expected for all high-energy plasmas with parameters
such that $k_{\mathrm{p},0}>\lambda B_{0}^{2}/k_{\mathrm{p},0}$.

The time evolution of the peak scale of the magnetic energy spectrum
$k_{\mathrm{p}}$ in R1, R1a, R1b, and R1mhd is presented in Figure 4. Here,
$k_{\mathrm{p}}$ is normalized by the analytical solution in the self-similar
phase, Equation (15). For R1b,
$k_{\mathrm{p}}/\left[7/(4t)\,\ln(t/t_{\mathrm{log}})\right]^{1/2}$ has a
constant value of $\approx 0.9$ for $t\gtrsim 3\times 10^{-3}$ which is
equivalent to a few timesteps of the simulation. Hence, Equation (15)
describes the evolution during the CME-assisted inverse cascade well. We show
the direct comparison with the scaling $k_{\mathrm{p}}\propto t^{-1/2}$ as the
black dashed line in Figure 4. The deviation from the $k_{\mathrm{p}}\propto
t^{-1/2}$ scaling is clearly visible in our simulations once they have entered
phase (iii) in which $k_{\mathrm{p}}$ and $\langle\mu_{5}\rangle$ evolve self-
similarly. This is strong evidence for the need of a logarithmic correction
that emerges naturally in our analytic derivation that is given in the
appendix.

For the classical MHD simulation, R1mhd,
$k_{\mathrm{p}}/(7/(4t)\,\ln(t/t_{\mathrm{log}}))^{1/2}$ is, at maximum, $0.5$
for $t\approx 3\times 10^{-3}$ and later decreases as $k_{\mathrm{p}}\propto
t^{-1/2}$. In R1a, which has the lowest chiral feedback parameter,
$k_{\mathrm{p}}/(7/(4t)\,\ln(t/t_{\mathrm{log}}))^{1/2}$ evolves very similar
to the MHD case, R1mhd. Initially, also
$k_{\mathrm{p}}/(7/(4t)\,\ln(t/t_{\mathrm{log}}))^{1/2}$ in R1 evolves similar
as in R1mhd. But at $t\approx 0.2$, the transition to phase (iii) occurs and
$k_{\mathrm{p}}/(7/(4t)\,\ln(t/t_{\mathrm{log}}))^{1/2}$ in R1 evolves similar
as in R1b.

##  IV Inverse cascade in chiral MHD with turbulence

In this section we explore the transition from laminar to turbulent flows. In
particular, we are interested in how turbulence modifies the three-phase
scenario of a decaying helical magnetic field in chiral MHD that was
established in Section III. Therefore we run a series of simulations where the
viscosity and Ohmic resistivity are systematically decreased. In the limit of
large Reynolds numbers, analytical estimates can be compared to the results
from DNS.

###  IV.1 Reynolds numbers in DNS of decaying (chiral) MHD turbulence

During decaying (chiral) MHD, the magnetic Reynolds number ${\rm
Re}_{{}_{\mathrm{M}}}$ is a function of time because (i) the decaying magnetic
field drives a velocity field which changes in time and (ii) the
characteristic wave number on which magnetic forcing occurs corresponds to the
correlation length of the magnetic field. The latter increases in time due to
the inverse cascade which occurs when the magnetic field is helical. In the
following, we approximate the correlation length of the magnetic field by the
scale at which the magnetic energy spectrum reaches its maximum,
$k_{\mathrm{p}}$, and define the time-dependent magnetic Reynolds number as

$\displaystyle{\rm
Re}_{{}_{\mathrm{M}}}(t)=\frac{u_{\mathrm{rms}}(t)}{k_{\mathrm{p}}(t)\eta}.$
(16)

The time evolution of ${\rm Re}_{{}_{\mathrm{M}}}$ in the majority of
simulations (all except R1a, R1b, and R8b) from this study is shown in Figure
5. Especially for the DNS with high diffusion, ${\rm Re}_{{}_{\mathrm{M}}}$
changes significantly during the simulation time. In our reference run for
chiral MHD with a vanishing velocity field, R1, ${\rm Re}_{{}_{\mathrm{M}}}$
decreases from a value of ${\rm Re}_{{}_{\mathrm{M}}}\approx 10^{-3}$ at the
beginning to ${\rm Re}_{{}_{\mathrm{M}}}\approx 10^{-5}$ at $t\approx 0.1$,
and then increases again, reaching ${\rm Re}_{{}_{\mathrm{M}}}\approx 10^{-1}$
at the final time of the simulation $t\approx 4$. The time dependence in runs
where ${\rm Re}_{{}_{\mathrm{M}}}$ is larger than unity in the beginning are
less dramatic. In the most turbulent run, R8, the magnetic Reynolds number
increases only by a factor of approximately $10$.

Figure 5: Magnetic Reynolds number ${\rm Re}_{{}_{\mathrm{M}}}$ as a function
of time. The colors refer to different simulations of classical and chiral MHD
simulations as given in the legend; see Table 1 for details. In the same
colors as the lines showing ${\rm Re}_{{}_{\mathrm{M}}}$, we indicate as
vertical dotted lines the time at which the minimal magnetic helicity is
reached and as vertical solid lines the time when the peak of the magnetic
energy spectrum $k_{\mathrm{p}}$ has reached the minimal wave number possible
in our numerical domain, $k_{p}=1$, respectively. Arrows on the vertical lines
indicate that these are characteristic times for classical MHD runs.

To distinguish the level of turbulence in different simulations, we may use
the maximum Reynolds number, respectively, which is given as

$\displaystyle{\rm
Re}_{{}_{\mathrm{M}}}^{\mathrm{max}}=\max\left(\frac{u_{\mathrm{rms}}(t)}{k_{\mathrm{p}}(t)\eta}\right).$
(17)

However, ${\rm Re}_{{}_{\mathrm{M}}}^{\mathrm{max}}$ is not a useful
characteristic of a simulation because it depends very much on the time at
which the simulation is stopped. A more consistent way of comparing different
simulations is by using Reynolds numbers that are defined at characteristic
times during the evolution. In the following, we will use the value of the
${\rm Re}_{{}_{\mathrm{M}}}$ at the time $t_{k_{\mathrm{p}}=1}$ at which the
peak of the magnetic energy spectrum reaches the minimum wave number within
the numerical domain:

$\displaystyle{\rm Re}_{{}_{\mathrm{M}}}^{{k_{\mathrm{p}}=1}}={\rm
Re}_{{}_{\mathrm{M}}}(t=t_{k_{\mathrm{p}}=1}).$ (18)

Additionally, we will consider the Reynolds number at the time
$t_{\mathrm{min}(\mathcal{H})}$ when the magnetic helicity reaches its
minimum:

$\displaystyle{\rm Re}_{{}_{\mathrm{M}}}^{{\mathrm{min}(\mathcal{H})}}={\rm
Re}_{{}_{\mathrm{M}}}(t=t_{\mathrm{min}(\mathcal{H})}).$ (19)

The values of ${\rm Re}_{{}_{\mathrm{M}}}^{\mathrm{max}}$, ${\rm
Re}_{{}_{\mathrm{M}}}^{{k_{\mathrm{p}}=1}}$, and ${\rm
Re}_{{}_{\mathrm{M}}}^{{\mathrm{min}(\mathcal{H})}}$ for all DNS presented in
this work are listed in the last three columns of Table 1.

Figure 6: Ratio of magnetic over kinetic energy as a function of time for the
same simulations as presented in Figure 5. See the caption of Figure 5 for a
description of the thin vertical lines.  Figure 7: The peak scale of the
magnetic energy spectrum as a function of time for the same simulations as
presented in Figure 5. See the caption of Figure 5 for a description of the
thin vertical lines.

###  IV.2 Transition from low to high ${\rm Re}_{{}_{\mathrm{M}}}$ in DNS of
decaying (chiral) MHD turbulence

Figure 8: Comparing runs R1, R2, R3, R4, R5, and R6, presented by blue data
points. In orange color, run R8 is shown which has the largest ${\rm
Re}_{{}_{\mathrm{M}}}$. Note, however, that R8 has a different value of
$\lambda B_{0}^{2}/(k_{\mathrm{p},0})$ in comparison to R1–R6; see Table 1.
Key properties of the runs are plotted as a function of the magnetic Reynolds
number ${\rm Re}_{{}_{\mathrm{M}}}$ for which three characteristic values are
presented: ${\rm Re}_{{}_{\mathrm{M}}}^{{k_{\mathrm{p}}=1}}$ is the magnetic
Reynolds number at the time when the peak of the magnetic energy spectrum
reaches the minimum value in the numerical domain, ${\rm
Re}_{{}_{\mathrm{M}}}^{{\mathrm{min}(\mathcal{H})}}$ is the magnetic Reynolds
number at the time when the minimum of $\langle\mathcal{H}\rangle$ is reached,
and ${\rm Re}_{{}_{\mathrm{M}}}^{\mathrm{max}}$ is the maximum Reynolds number
during the entire simulation time. The latter depends strongly on time when
the simulation is stopped and is not suitable for a comparison between
different runs.
a) Maximum value of $2\langle\mu_{5}\rangle/\lambda$ over the conserved
quantity $\langle\mathcal{H}\rangle+2\langle\mu_{5}\rangle/\lambda$.
b) Maximum value of $\langle\mu_{5}\rangle$ generated in the simulation over
the theoretically predicted value in the kinematic limit
$\mathrm{min}(k_{\mathrm{p},0},\lambda B_{0}^{2}/k_{\mathrm{p},0})$.  Figure
9: Using the same conventions as for Figure 8 we show:
a) Difference between the maximum and the minimum value of the logarithm of
the magnetic helicity $\langle\mathcal{H}\rangle$.
b) Time until the minimum of $\langle\mathcal{H}\rangle$ is reached
$t_{\mathrm{min}(\bf{A}\cdot\bf{B})}$.
c) Maximum of the ratio of magnetic over kinetic energy for
$t>t_{\mathrm{min}(\bf{A}\cdot\bf{B})}$.

Starting from our reference run of chiral MHD with negligible velocity field,
R1, we systematically decrease the values of $\eta$ and $\nu$ in runs R2–R8 in
order to explore the transition to the turbulent regime where the velocity
field is expected to impact significantly the evolution of the magnetic field.
As a characteristic parameter for the degree of nonlinearity in the Equations
(6)–(9), we list the ratio of the initial magnetic field strength $B_{0}$ over
$\eta$ for all DNS in the 7th column of Table 1. With $B_{0}$ being the Alfvén
velocity within our unit system, $B_{0}/\eta$ can be considered as the initial
Alfvénic Reynolds number. The time evolution of the ${\rm
Re}_{{}_{\mathrm{M}}}$ for runs R2–R8 is presented in Figure 5. The value of
${\rm Re}_{{}_{\mathrm{M}}}$ is larger than $1$ in the entire simulation time
of R6, R7, and R8. In the latter run, the simulation reaches ${\rm
Re}_{{}_{\mathrm{M}}}\approx 270$.

When increasing the Reynolds number in the simulations, we observe two trends.
First, the maximal ratio of magnetic over kinetic energy density decreases
when $\eta$ and $\nu$ are decreased 444This excludes, of course, the very
first time step of the simulations where the kinetic energy is zero., see
Figure 6. Second, the scaling of the peak scale of the magnetic energy with
time changes from
$k_{\mathrm{p}}\propto(7/(4t)\,\mathrm{log}(t/t_{\mathrm{log}}))^{1/2}$ at
late times to $k_{\mathrm{p}}\propto t^{-2/3}$; see Figure 7.

The changes in the plasma evolution at different Reynolds numbers is also
clearly visible in Figure 8. Different measured characteristics of the
simulations are presented here as a function of ${\rm Re}_{{}_{\mathrm{M}}}$.
According to our discussion in Section IV.1, for each simulation these
parameters are plotted as a function of ${\rm
Re}_{{}_{\mathrm{M}}}^{\mathrm{max}}$ (open dots), ${\rm
Re}_{{}_{\mathrm{M}}}^{{k_{\mathrm{p}}=1}}$ (filled dots), and ${\rm
Re}_{{}_{\mathrm{M}}}^{{\mathrm{min}(\mathcal{H})}}$ (stars). In Figure 8a,
the maximum value of $2\langle\mu_{5}\rangle/\lambda$ over the conserved total
chirality $\langle\mathcal{H}\rangle+2\langle\mu_{5}\rangle/\lambda$ is
presented. While $2\langle\mu_{5}\rangle/\lambda$ is almost 100 percent of the
total chirality at one time of the plasma evolution for ${\rm
Re}_{{}_{\mathrm{M}}}<1$, the maximum of the ratio
$\left(2\langle\mu_{5}\rangle/\lambda\right)/\left(\langle\mathcal{H}\rangle+2\langle\mu_{5}\rangle/\lambda\right)$
drops to about $1/2$ once ${\rm Re}_{{}_{\mathrm{M}}}$ becomes significantly
larger than $1$. The mean magnetic helicity, on the other hand, always
dominates the total chirality at one point in time for all simulations
regardless of their degree of turbulence as expected from the choice of your
initial conditions. In simulations with ${\rm Re}_{{}_{\mathrm{M}}}>1$, the
maximum value of $2\langle\mu_{5}\rangle/\lambda$ is never reached, as is
shown in Figure 8b. The difference between the logarithm of the maximum and
the minimum of $\langle\mathcal{H}\rangle$, see Figure 9a, and also the time
needed to reach the minimum of $\langle\mathcal{H}\rangle$ drops for ${\rm
Re}_{{}_{\mathrm{M}}}>1$, see Figure 9b. Overall, we observe a decrease of the
maximal ratio of magnetic over kinetic energy in our simulations with
increasing magnetic Reynolds number. This ratio is plotted in Figure 9c and
decreases continuously with decreasing Ohmic resistivity and not suddenly at
the transition ${\rm Re}_{{}_{\mathrm{M}}}\approx 1$.

###  IV.3 Analytical estimates for the limit of large ${\rm
Re}_{{}_{\mathrm{M}}}$

Let us estimate analytically the time evolution of $\mu_{5}$ for high magnetic
Reynolds numbers. A high ${\rm Re}_{{}_{\mathrm{M}}}$ implies that the first
term is more important than the second term in the right-hand side of Equation
(6). In the case of the vanishing chiral chemical potential at the initial
time, $\mu_{5,0}=0$, as we have seen in the previous section, the third term
never overwhelms the second term in Equation (6), because $\mu_{5}(t)$ would
decay for $k_{\mathrm{p}}\ll\mu_{5}$. Therefore, the evolution of the magnetic
field is governed by the interaction to the fluid velocity $\bm{U}$ in the
same way as the (nonchiral) classical ideal MHD, and the magnetic fields
undergo the classical inverse cascade.

Figure 10: Comparing runs with high Reynolds numbers: Classical MHD (left
panels, run R8mhd) and chiral MHD (right panels, run R8). Top: Time evolution
of the mean magnetic helicity $\langle\mathcal{H}\rangle$ and the wave number
at which the magnetic energy spectrum has its maximum $k_{\mathrm{p}}$. For
chiral MHD, we also show $\langle\mu_{5}\rangle/\lambda$ as well as the
conserved quantity $\langle\mathcal{H}\rangle+2\langle\mu_{5}\rangle/\lambda$.
The solid vertical lines indicate the time when $k_{\mathrm{p}}$ has reached
its minimum, $k=1$, and the dotted vertical lines indicate the time at which
$\langle\mathcal{H}\rangle$ is minimal. 2nd row: $k_{\mathrm{p}}$ and
$\langle{\mu_{5}}\rangle$ normalized by the theoretically expected scaling
during self-similar evolution (see Section III.1). 3rd row: Evolution of the
magnetic energy spectrum. The thick dotted black lines show the initial
spectra and the thick solid black lines show the final spectra of the
simulations. Bottom: The same but for the kinetic energy spectrum.

As a simple model of the classical inverse cascade (see, e.g. Durrer and
Neronov, 2013), we consider the following behaviors,

$\displaystyle k_{\mathrm{p}}(t)$
$\displaystyle=\left\\{\begin{array}[]{cc}k_{\mathrm{p},0}&(t\leq
t_{\mathrm{I}})\\\
k_{\mathrm{p},0}(t/t_{\mathrm{I}})^{-2/3}&(t_{\mathrm{I}}\leq t)\\\
\end{array}\right.,$ (22) $\displaystyle B(t)$
$\displaystyle=\left\\{\begin{array}[]{cc}B_{0}&(t\leq t_{\mathrm{I}})\\\
B_{0}(t/t_{\mathrm{I}})^{-1/3}&(t_{\mathrm{I}}\leq t)\\\ \end{array}\right.,$
(25)

where $t_{\mathrm{I}}$ denotes the onset time of the classical inverse
cascade. Inserting this into Equation (9) and solving for $\mu_{5}$, we obtain

$\displaystyle\mu_{5}(t\leq t_{\mathrm{I}})$
$\displaystyle=k_{\mathrm{p},0}\left[1-e^{-\xi t/t_{\mathrm{I}}}\right],$ (26)
$\displaystyle\mu_{5}(t\geq t_{\mathrm{I}})$
$\displaystyle=k_{\mathrm{p},0}\Big{[}-3\xi(t/t_{\mathrm{I}})^{-1/3}$
$\displaystyle+e^{-3\xi(t/t_{\mathrm{I}})^{1/3}}\Big{\\{}\left(e^{\xi}(1+3\xi)-1\right)e^{2\xi}$
$\displaystyle+9\xi^{2}\left({\rm Ei}(3\xi(t/t_{\mathrm{I}})^{1/3})-{\rm
Ei}(3\xi)\right)\Big{\\}}\Big{]},$ (27)

where $\xi\equiv\lambda\eta B_{0}^{2}t_{\mathrm{I}}$ and ${\rm Ei}(x)$ is the
exponential integral function. It is interesting to consider the late time
limit of this solution,

$\mu_{5}(t\gg\xi^{-3}t_{\mathrm{I}})\simeq
k_{\mathrm{p},0}(t/t_{\mathrm{I}})^{-2/3},$ (28)

which exactly matches $k_{\mathrm{p}}(t)$ in Equation (22). Although
$\mu_{5}(t)\approx k_{\mathrm{p}}(t)$ holds at late times irrespective of
${\rm Re}_{{}_{\mathrm{M}}}$, their time evolution Equations (15) and (28) are
quite different. Note that since the sudden change of the behaviors in
Equations (22) and (25) at $t=t_{\mathrm{I}}$ are crude approximations, we
anticipate a slight deviation between the analytic estimates of $\mu_{5}(t)$
and the DNS results there.

###  IV.4 Simulations of chiral helical MHD turbulence

In Figure 10, the run with lowest diffusion, hence highest ${\rm
Re}_{{}_{\mathrm{M}}}$, (R8, left panels) is compared to a classical MHD
analogue (R8mhd, right panels). The analysis is exactly the same as in Figure
1, except for the addition of the kinetic energy spectra in the last row of
Figure 10.

Due to the small but finite value of the resistivity in R8mhd,
$\langle\mathcal{H}\rangle$ decays by roughly a factor of three over the
entire simulation time. Nevertheless, the magnetic helicity, in combination
with turbulence leads to an efficient inverse cascade in energy which can be
seen in the evolution of the magnetic energy spectrum in Figure 10c. The
scaling of $k_{\mathrm{p}}$ proportional to $t^{-2/3}$ as expected for the
turbulent inverse cascade of in helical MHD turbulence, sets in at time
$t\approx 10^{-5}$. This coincides roughly with the minimum of the eddy turn
over time in the simulation, hence we will use

$\displaystyle t_{\mathrm{I}}\approx
t_{\mathrm{eddy}}^{\mathrm{min}}\equiv\min\left(\frac{1}{k_{\mathrm{p}}(t)u_{\mathrm{rms}}(t)}\right).$
(29)

The time evolution of $k_{\mathrm{p}}$ normalized by
$k_{\mathrm{p,0}}(t/t_{\mathrm{eddy}}^{\mathrm{min}})^{-2/3}$ is presented in
Figure 10b The scaling with $t^{-2/3}$ is observed in our DNS for times later
than approximately $t_{\mathrm{eddy}}^{\mathrm{min}}$. The value of
$k_{\mathrm{p}}$ reaches the minimum value of the box after $t\lesssim 0.2$.

The time evolution and energy spectra of the chiral MHD run with highest
magnetic Reynolds number, R8, that are presented in the right panels of Figure
10 are very similar to the ones in the classical MHD run R8mhd. Up $t\approx
10^{-3}$, $k_{\mathrm{p}}$, $\langle\mathcal{H}\rangle$, and
$t_{\mathrm{eddy}}$ evolve identically in turbulent MHD and turbulent chiral
MHD. However, in R8 a $\langle\mu_{5}\rangle$ is generated and restores a
small amount of $\langle\mathcal{H}\rangle$. The energy spectra in R8mhd and
R8 are indistinguishable 555Note, that small gaps in time seem to occur in the
energy spectra. This is an artefact from the code that calculates energy
spectra at a manually fixed time interval. The simulations were restarted
several time and, since they cover many orders of magnitude in time, the
frequency of writing spectra is reduced at every restart to save computing
power. The gap appears after a restart when the sampling time step was reduced
slightly too much to homogeneously fill the logarithmically spaced array of
example spectra.. As expected from Equation (26), $\langle\mu_{5}\rangle$
increases linearly in time in the beginning. However the scaling proportional
to $t^{-2/3}$ as expected for late times according to Equation (28) is not
observed in R8. This is caused by the fact that the peak of the magnetic
energy spectrum has moved to the minimum wave number, $k_{\mathrm{p}}=1$,
before the scaling of $\langle\mu_{5}\rangle$ could converge to the one of
$k_{\mathrm{p}}$.

To test the late time scaling of $\langle\mu_{5}\rangle\propto t^{-2/3}$ we
have repeated run R8 with a larger value of $\lambda$. For larger $\lambda$,
the condition $t\gg\xi^{-3}t_{\mathrm{I}}$ is fulfilled while the inverse
cascade still proceeds within the numerical domain. Run R8b has a value of
$\lambda$ that is $10^{2}$ times larger than the one in R8. We compare these
two runs with the MHD analog, R8mhd, in Figure 11. The time evolution of
$k_{\mathrm{p}}$ in all three runs is almost identical, reaching a scaling of
$k_{\mathrm{p}}\propto t^{-2/3}$ at $t\gtrsim 10^{-5}$. Except for the time
around the onset of the inverse cascade $t_{\mathrm{I}}$, the time evolution
of $\langle\mu_{5}\rangle$ measured in DNS (solid lines) agrees very well with
the theoretically predicted curves (dotted lines). However, the value in DNS
is approximately larger by a factor of $1.5$ compared to the result from
Equation (26) at early times. This behavior might be corrected when including
the exact shape of the initial magnetic energy spectrum. The evolution of
$\langle\mu_{5}\rangle$ in R8b is very well described by Equation (27) for
$t\gtrsim 10^{-4}$.

Figure 11: The mean chiral chemical potential, $\langle\mu_{5}\rangle$, as a
function of time for runs R8 and R8b which differ only in their value of
$\lambda$ (see Table 1). Additionally, the evolution of $k_{\mathrm{p}}$ is
shown as dashed lines. For direct comparison also $k_{\mathrm{p}}$ from the
corresponding classical MHD run, R8mhd, is presented. The solid vertical lines
indicate the time when $k_{\mathrm{p}}$ has reached the minimum wave number
possible in the numerical domain, $k=1$, and the dotted vertical lines
indicates the time at which the minimum of magnetic helicity is reached. The
horizontal gray dotted line indicates the threshold for $\mu_{5}$ above which
a dynamo instability occurs in the numerical box.

##  V Conclusion

In this study we have explored the evolution of a decaying fully helical
magnetic field in a high-energy plasma in which the chiral magnetic effect can
occur. The chiral magnetic effect is a macroscopic quantum effect that
describes the emergence of an electric current in the presence of a chiral
asymmetry, e.g. a difference between the chemical potential of left- and
right-handed fermions, $\mu_{5}$, and a magnetic field. The dynamics of such a
plasma is determined by the system of equations of chiral MHD. We have
investigated how a magnetic field decays in chiral MHD and how a chiral
asymmetry develops, from vanishing initial $\mu_{5}$.

When the velocity field plays no dynamically important role, we have
identified a three phase evolution: (i) a mean chiral chemical potential
$\langle\mu_{5}\rangle$ is produced at the expense of the mean magnetic
helicity $\langle\mathcal{H}\rangle$, (ii) once $\langle\mu_{5}\rangle$
exceeds the inverse correlation length of the helical magnetic field
$k_{\mathrm{p}}$ a chiral dynamo reestablishes $\langle\mathcal{H}\rangle$,
and (iii) a self-similar inverse cascade where $|\mu_{5}(t)|\approx
k_{\mathrm{p}}(t)\approx\left[(3+n)/(4\eta
t)~{}\ln\left(t/t_{\mathrm{log}}\right)\right]^{1/2}$, where $n$ is the slope
of the initial magnetic helicity spectrum. A similar scenario has been
reported by Hirono et al. (2015). However, they did not include the
logarithmic correction in phase (iii) that agrees well with our DNS.

Our simulations, performed with the Pencil Code that has explicit viscosity
and Ohmic resistivity, allow us to systematically explore a decaying magnetic
field in chiral MHD and, in particular, to probe the transition from low to
high magnetic Reynolds numbers. We find that the magnetic energy spectrum
evolves more and more similar in chiral MHD and classical MHD when ${\rm
Re}_{{}_{\mathrm{M}}}$ is increased by lowering the dissipation in DNS. For
${\rm Re}_{{}_{\mathrm{M}}}\gg 1$, we observe an inverse cascade of the
magnetic field, where the peak of the magnetic energy spectrum
$k_{\mathrm{p}}\propto t^{-2/3}$, regardless of the existence of an additional
degree of freedom in from of a $\mu_{5}$. In the simulations of chiral MHD, a
$\mu_{5}$ is initially generated linearly in time. For
$t\gg\xi^{-3}t_{\mathrm{I}}$ with $\xi\equiv\lambda\eta
B_{0}^{2}t_{\mathrm{I}}$ we find a scaling of $\mu_{5}\simeq
k_{\mathrm{p},0}(t/t_{\mathrm{I}})^{-2/3}$, i.e., like for low Reynolds
number, the evolution of $\mu_{5}$ follows the one of $k_{\mathrm{p}}$,
$k_{\mathrm{p}}\simeq\mu_{5}$.

We have shown that the chiral magnetic effect not only modifies the inverse
cascade of magnetic fields but it also leads to the generation of a chiral
anomaly. Such an anomaly may manifest itself in the neutrino density in the
late Universe and it might be relevant, e.g., for the the number of effective
degrees of freedom, $N_{\rm eff}$, at recombination. A study of this
possibility is referred to future work. Our results may also have other
important implications for the evolution of relativistic plasmas in the early
Universe and protoneutron stars.

###### Acknowledgements.

The authors would like to acknowledge the Mainz Institute for Theoretical
Physics (MITP) of the DFG Cluster of Excellence PRISMA+ (Project ID 39083149),
for enabling us to complete a significant portion of this work. JS
acknowledges the funding from the Swiss National Science Foundation under
Grant No. 185863, European Union’s Horizon 2020 research and innovation
program under the Marie Skłodowska-Curie grant No. 665667, and the support by
the National Science Foundation under Grant No. NSF PHY-1748958. The work of
TF was supported by JSPS KAKENHI No. 17J09103 and No. 18K13537. RD is
supported with the Swiss National Science Foundation under Grant No. 200020
182044\. The simulations presented in this work were performed on resources at
Chalmers Centre for Computational Science and Engineering (C3SE) provided by
the Swedish National Infrastructure for Computing (SNIC) as well as on the
Baobab cluster at the University of Geneva.

##  Appendix A Derivation of the self-similar evolution in the limit of a
vanishing velocity field

Here we shall derive the self-similar solution of $\mu_{5}(t)$ and
$k_{\mathrm{p}}(t)$ in the phase (iii), Equation (15). Ignoring the velocity
field in Equation (6), the equation of motion for magnetic field reads

$\partial_{t}\bm{B}=\eta\bm{\nabla}\times\left(-\bm{\nabla}\times\bm{B}+\mu_{5}\bm{B}\right).$
(30)

The Fourier transformation of the magnetic field is written as

$\bm{B}(t,\bm{x})=\sum_{\lambda=\pm}\int\frac{{\rm
d}^{3}k}{(2\pi)^{3}}e^{i\bm{k\cdot
x}}\bm{e}_{\lambda}(\hat{\bm{k}})B_{\lambda}(t,k),$ (31)

where $\lambda=\pm$ is the label of the circular polarization and
$\bm{e}_{\lambda}(\hat{\bm{k}})$ is the circular polarization vector that
satisfies $i\bm{k}\times\bm{e}_{\pm}(\hat{\bm{k}})=\pm
k\bm{e}_{\pm}(\hat{\bm{k}})$. Using the relation between the magnetic helicity
and $B_{\pm}$ in Fourier space,
$\mathcal{H}_{k}\equiv\left(B_{+}^{2}-B_{-}^{2}\right)/k$, one can recast
Equation (30) into the equation of motion for the magnetic helicity,

$\partial_{t}\mathcal{H}_{k}+2\eta
k^{2}\mathcal{H}_{k}-4\eta\mu_{5}\,\rho_{B}(k)=0,$ (32)

where the spatial fluctuations of $\mu_{5}$ are neglected (i.e.
$\mu_{5}\rightarrow\langle\mu_{5}\rangle$) and
$\rho_{B}(k)\equiv(B_{+}^{2}+B_{-}^{2})/2$ is the magnetic energy density in
Fourier space. When the magnetic field is maximally helical, $B_{\pm}^{2}\gg
B_{\mp}^{2}$, $\rho_{B}(k)\simeq\sigma_{\mathcal{H}}\,k\mathcal{H}_{k}/2$,
where $\sigma_{\mathcal{H}}\equiv{\rm sign}[\mathcal{H}_{k}]$ ensures
$\rho_{B}$ is always positive. In this case, the above equation is reduced
into

$\partial_{t}\mathcal{H}_{k}+2k\eta(k-|\mu_{5}|)\mathcal{H}_{k}=0,\quad{\rm(maximal\,helical)}$
(33)

where we used $\sigma_{\mathcal{H}}\mu_{5}=|\mu_{5}|$ because the magnetic
helicity and $\mu_{5}$ have the same sign, if $\mu_{5}$ is produced by the
magnetic helicity as we assume. This equation has the formal analytic solution

$\mathcal{H}_{k}(t)=\mathcal{H}_{k}(t_{0})\exp\left[2k\eta\left(-k(t-t_{0})+\theta(t)\right)\right],$
(34)

with
$\theta(t)\equiv\int^{t}_{t_{0}}\mathrm{d}t^{\prime}|\mu_{5}(t^{\prime})|.$ In
the phase (iii) the magnetic helicity dominates the conserved quantity and
hence $\mathcal{H}(\gg 2\mu_{5}/\lambda)$ is independently conserved

$\int\mathrm{d}k\,k^{2}\mathcal{H}_{k}(t_{0})\exp\left[2k\eta\left(-k(t-t_{0})+\theta(t)\right)\right]={\rm
const}.$ (35)

Taking the time derivative of this equation and dropping time dependent but
nonvanishing factors, one finds

$\int^{k_{*}}_{0}\mathrm{d}k\,k^{3+n}\left[|\mu_{5}(t)|-k\right]e^{-2\eta(t-t_{0})\left(k-\frac{\theta(t)}{2(t-t_{0})}\right)^{2}}=0.$
(36)

Here we assume a power-law helicity slope with an UV-cutoff at $k_{*}$

$\mathcal{H}_{k}(t_{0})=\mathcal{H}_{k_{*}}(t_{0})\left(\frac{k}{k_{*}}\right)^{n}\Theta(k_{*}-k),$
(37)

where $\Theta(x)$ is the Heaviside function. Changing the dummy variable from
$k$ into $p\equiv\sqrt{2\eta(t-t_{0})}k$, we obtain

$|\mu_{5}(t)|=\frac{1}{\sqrt{2\eta(t-t_{0})}}\,\frac{\int^{p_{*}}_{0}\mathrm{d}p\,p^{4+n}e^{-\left(p-\sqrt{\frac{\eta}{2(t-t_{0})}}\,\theta(t)\right)^{2}}}{\int^{p_{*}}_{0}\mathrm{d}p\,p^{3+n}e^{-\left(p-\sqrt{\frac{\eta}{2(t-t_{0})}}\,\theta(t)\right)^{2}}},$
(38)

with $p_{*}\equiv\sqrt{2\eta(t-t_{0})}k_{*}$.

To simplify this expression, we make an additional approximation. For this, we
restrict ourselves into a late time regime, $t\gg t_{0}$. Then the upper limit
of the integrals $p_{*}\propto t^{1/2}$ can be approximated by $\infty$, and
$\sqrt{\eta/2t}\,\theta$ is also considered as large, as we will confirm a
posteriori. The integrals are computed for
$X\equiv\sqrt{\eta/(2(t-t_{0}))}\,\theta(t)\gg 1$ as

$\frac{\int^{\infty}_{0}\mathrm{d}p\,p^{4+n}\exp\left[-\left(p-X\right)^{2}\right]}{\int^{\infty}_{0}\mathrm{d}p\,p^{3+n}\exp\left[-\left(p-X\right)^{2}\right]}=X+\frac{3+n}{X}+\mathcal{O}(X^{-3}).$
(39)

Therefore Equation (38) is simplified to

$\theta^{\prime}(t)=|\mu_{5}(t)|\simeq\frac{\theta(t)}{2t}+\frac{3+n}{2\eta\theta(t)}.$
(40)

Note that in order for $\mathcal{H}$ not to develop an infrared singularity we
must require $n>-3$ so that both terms in (40) are always positive. The
solution of this differential equation is

$\displaystyle\theta(t)$
$\displaystyle\simeq\sqrt{\frac{t}{\eta}}\sqrt{\mathcal{C}+(3+n)\ln\left(\frac{t}{t_{\mathcal{C}}}\right)},$
(41) $\displaystyle\mu_{5}(t)$ $\displaystyle\simeq\frac{1}{2\sqrt{\eta
t}}\frac{3+n+\mathcal{C}+(3+n)\ln(t/t_{\mathcal{C}})}{\sqrt{\mathcal{C}+(3+n)\ln(t/t_{\mathcal{C}})}},$
(42)

where $\mathcal{C}$ is an integration constant and $t_{\mathcal{C}}$ is
degenerate with $\mathcal{C}$. The approximation used above,
$X=\sqrt{\eta/2t}\,\theta\gg 1$, is valid for a sufficiently late time,

$\sqrt{\frac{\eta}{t}}\theta(t)\simeq\sqrt{\mathcal{C}+(3+n)\ln\left(\frac{t}{t_{\mathcal{C}}}\right)}\gg
1.$ (43)

This also allows us to further simplify $\mu_{5}$ as

$\mu_{5}\simeq\left[\frac{3+n}{4\eta t}\ln\left(\frac{t}{t_{\rm
log}}\right)\right]^{\frac{1}{2}},$ (44)

where the integration constant is rewritten as
$\mathcal{C}=(3+n)\ln(t_{\mathcal{C}}/t_{\rm log})$, assuming $n\neq-3$. Note
that this logarithmic correction which slightly slows down the decay of
$\mu_{5}$ becomes more significant as $n$ increases. It implies that also the
inverse cascade (i.e. the transportation to larger scales) of the peak scale
$k_{p}$ is slowed-down. This is because it takes more time for a large-scale
helicity modes to grow large enough to ensure the conservation law, when the
initial helicity has a bluer spectrum, i.e. more power on smaller scales. If
the initial helicity is scale invariant, $n=-3$, the logarithmic correction
vanishes.

## References

  * Mac Low and Klessen (2004) M.-M. Mac Low and R. S. Klessen, Rev. Mod. Phys. 76, 125 (2004), URL https://link.aps.org/doi/10.1103/RevModPhys.76.125.
  * Klessen and Hennebelle (2010) R. S. Klessen and P. Hennebelle, A & A 520, A17 (2010), eprint 0912.0288.
  * Stevenson (2003) D. J. Stevenson, Earth and Planet. Sci. Let. 208, 1 (2003).
  * Donati and Landstreet (2009) J.-F. Donati and J. D. Landstreet, Ann. Rev. Astron. & Astrophys. 47, 333 (2009).
  * Solanki et al. (2006) S. K. Solanki, B. Inhester, and M. Schüssler, Reports on Progress in Physics 69, 563 (2006), eprint 1008.0771.
  * Crutcher (2012) R. M. Crutcher, Ann. Rev. Astron. & Astrophys. 50, 29 (2012).
  * Beck (2012) R. Beck, Space Sci. Rev. 166, 215 (2012).
  * Govoni and Feretti (2004) F. Govoni and L. Feretti, International Journal of Modern Physics D 13, 1549 (2004).
  * Neronov and Vovk (2010) A. Neronov and I. Vovk, Science 328, 73 (2010), eprint 1006.3504.
  * Brandenburg and Subramanian (2005) A. Brandenburg and K. Subramanian, Phys. Rept. 417, 1 (2005).
  * Turner and Widrow (1988) M. S. Turner and L. M. Widrow, Phys. Rev. D 37, 2743 (1988).
  * Sigl et al. (1997) G. Sigl, A. V. Olinto, and K. Jedamzik, Phys. Rev. D 55, 4582 (1997).
  * Durrer and Neronov (2013) R. Durrer and A. Neronov, Astron. Astrophys. Rev. 21, 62 (2013), eprint 1303.7121.
  * Pouquet et al. (1976) A. Pouquet, U. Frisch, and J. Léorat, Journal of Fluid Mechanics 77, 321–354 (1976).
  * Christensson et al. (2001) M. Christensson, M. Hindmarsh, and A. Brandenburg, Phys. Rev. E 64 (2001).
  * Biskamp and Müller (1999) D. Biskamp and W.-C. Müller, Phys. Rev. Lett. 83 (1999).
  * Müller and Biskamp (2000) W.-C. Müller and D. Biskamp, Phys. Rev. Lett. 84, 475 (2000), URL https://link.aps.org/doi/10.1103/PhysRevLett.84.475.
  * Brandenburg et al. (2015) A. Brandenburg, T. Kahniashvili, and A. G. Tevzadze, Phys. Rev. Lett. 114, 075001 (2015), eprint 1404.2238.
  * Reppin and Banerjee (2017) J. Reppin and R. Banerjee, Phys. Rev. E 96, 053105 (2017), eprint 1708.07717.
  * Kahniashvili et al. (2013) T. Kahniashvili, A. G. Tevzadze, A. Brandenburg, and A. Neronov, Phys. Rev. D 87, 083007 (2013).
  * Brandenburg et al. (2017a) A. Brandenburg, T. Kahniashvili, S. Mandal, A. R. Pol, A. G. Tevzadze, and T. Vachaspati, Phys. Rev. D 96, 123528 (2017a), eprint 1711.03804.
  * Boyarsky et al. (2012) A. Boyarsky, J. Fröhlich, and O. Ruchayskiy, Phys. Rev. Lett. 108 (2012).
  * Vilenkin (1980) A. Vilenkin, Phys. Rev. D 22, 3080 (1980).
  * Son and Surówka (2009) D. T. Son and P. Surówka, Physical Review Letters 103, 191601 (2009), eprint 0906.5044.
  * Joyce and Shaposhnikov (1997) M. Joyce and M. E. Shaposhnikov, Phys. Rev. Lett. 79, 1193 (1997).
  * Fröhlich and Pedrini (2000) J. Fröhlich and B. Pedrini, in _Mathematical Physics 2000_ , edited by A. S. Fokas, A. Grigoryan, T. Kibble, and B. Zegarlinski (World Scientific Publishing Company, 2000), International Conference on Mathematical Physics 2000, Imperial college (London).
  * Fröhlich and Pedrini (2002) J. Fröhlich and B. Pedrini, in _Statistical Field Theory_ , edited by A. Cappelli and G. Mussardo (Kluwer, 2002).
  * Semikoz and Sokoloff (2004) V. B. Semikoz and D. D. Sokoloff, A & A 433, L53 (2004).
  * Semikoz et al. (2009) V. B. Semikoz, D. D. Sokoloff, and J. W. F. Valle, Phys. Rev. D80, 083510 (2009).
  * Boyarsky et al. (2012) A. Boyarsky, O. Ruchayskiy, and M. Shaposhnikov, Phys. Rev. Lett. 109, 111602 (2012).
  * Semikoz et al. (2012) V. Semikoz, D. Sokoloff, and J. Valle, JCAP 1206, 008 (2012).
  * Tashiro et al. (2012) H. Tashiro, T. Vachaspati, and A. Vilenkin, Phys. Rev. D86, 105033 (2012).
  * Dvornikov and Semikoz (2012) M. Dvornikov and V. B. Semikoz, JCAP 1202, 040 (2012).
  * Dvornikov and Semikoz (2013) M. Dvornikov and V. B. Semikoz, Phys. Rev. D87, 025023 (2013).
  * Dvornikov and Semikoz (2014) M. Dvornikov and V. B. Semikoz, JCAP 1405, 002 (2014).
  * Dvornikov and Semikoz (2017) M. Dvornikov and V. B. Semikoz, Phys. Rev. D95, 043538 (2017).
  * Manuel and Torres-Rincon (2015) C. Manuel and J. M. Torres-Rincon, Phys. Rev. D92, 074018 (2015).
  * Gorbar et al. (2016) E. V. Gorbar, I. Rudenok, I. A. Shovkovy, and S. Vilchinskii, Phys. Rev. D94, 103528 (2016).
  * Pavlović et al. (2016) P. Pavlović, N. Leite, and G. Sigl, JCAP 1606, 044 (2016).
  * Pavlović et al. (2017) P. Pavlović, N. Leite, and G. Sigl, Phys. Rev. D 96, 023504 (2017), eprint 1612.07382.
  * Rogachevskii et al. (2017) I. Rogachevskii, O. Ruchayskiy, A. Boyarsky, J. Fröhlich, N. Kleeorin, A. Brandenburg, and J. Schober, Astrophys. J. 846, 153 (2017).
  * Figueroa et al. (2019) D. G. Figueroa, A. Florio, and M. Shaposhnikov, Journal of High Energy Physics 2019, 142 (2019), eprint 1904.11892.
  * Brandenburg et al. (2017b) A. Brandenburg, J. Schober, I. Rogachevskii, T. Kahniashvili, A. Boyarsky, J. Fröhlich, O. Ruchayskiy, and N. Kleeorin, ApJL 845, L21 (2017b), eprint 1707.03385.
  * Schober et al. (2018) J. Schober, I. Rogachevskii, A. Brandenburg, A. Boyarsky, J. Fröhlich, O. Ruchayskiy, and N. Kleeorin, Astrophys. J. 858, 124 (2018), eprint 1711.09733.
  * Masada et al. (2018) Y. Masada, K. Kotake, T. Takiwaki, and N. Yamamoto, ArXiv e-prints (2018), eprint 1805.10419.
  * Kharzeev et al. (2008) D. E. Kharzeev, L. D. McLerran, and H. J. Warringa, Nucl. Phys. A803, 227 (2008).
  * Galitski et al. (2018) V. Galitski, M. Kargarian, and S. Syzranov, Phys. Rev. Lett. 121, 176603 (2018), eprint 1804.09339.
  * Hirono et al. (2015) Y. Hirono, D. E. Kharzeev, and Y. Yin, Phys. Rev. D 92, 125031 (2015), eprint 1509.07790.
  * Mace et al. (2019) M. Mace, N. Mueller, S. Schlichting, and S. Sharma, arXiv e-prints arXiv:1910.01654 (2019), eprint 1910.01654.
  * Field and Carroll (2000) G. B. Field and S. M. Carroll, Phys. Rev. D 62, 103008 (2000), eprint astro-ph/9811206.
  * Anber and Sorbo (2006) M. M. Anber and L. Sorbo, JCAP 2006, 018 (2006), eprint astro-ph/0606534.
  * Durrer et al. (2011) R. Durrer, L. Hollenstein, and R. K. Jain, JCAP 1103, 037 (2011), eprint 1005.5322.
  * Caprini and Sorbo (2014) C. Caprini and L. Sorbo, JCAP 1410, 056 (2014), eprint 1407.2809.
  * Fujita et al. (2015) T. Fujita, R. Namba, Y. Tada, N. Takeda, and H. Tashiro, JCAP 1505, 054 (2015), eprint 1503.05802.
  * Adshead et al. (2016) P. Adshead, J. T. Giblin, T. R. Scully, and E. I. Sfakianakis, JCAP 1610, 039 (2016), eprint 1606.08474.
  * Fujita and Durrer (2019) T. Fujita and R. Durrer, JCAP 1909, 008 (2019), eprint 1904.11428.
  * Domcke and Mukaida (2018) V. Domcke and K. Mukaida, JCAP 1811, 020 (2018), eprint 1806.08769.
  * Frisch et al. (1975) U. Frisch, A. Pouquet, J. LÉOrat, and A. Mazure, Journal of Fluid Mechanics 68, 769–778 (1975).
  * Biskamp (2003) D. Biskamp, _Magnetohydrodynamic Turbulence_ (Cambridge University Press, 2003).
  * Zrake (2014) J. Zrake, The Astrophysical Journal 794, L26 (2014), URL https://doi.org/10.1088%2F2041-8205%2F794%2F2%2Fl26.
  * Berera and Linkmann (2014) A. Berera and M. Linkmann, Phys. Rev. E 90, 041003 (2014), eprint 1405.6756.
  * Boyarsky et al. (2015) A. Boyarsky, J. Fröhlich, and O. Ruchayskiy, Phys. Rev. D 92, 043004 (2015).
  * Williamson (1980) J. H. Williamson, J. Comp. Phys. 35, 48 (1980).
  * Brandenburg and Dobler (2002) A. Brandenburg and W. Dobler, Comp. Phys. Comm. 147, 471 (2002), eprint astro-ph/0111569.
  * Brandenburg (2003) A. Brandenburg, in _Advances in Nonlinear Dynamics_ , edited by A. Ferriz-Mas and M. Núñez (CRC Press, 2003), pp. 269–344, ISBN 9780203493137, ISSN 2155-4692, URL http://dx.doi.org/10.1201/9780203493137.ch9.
  * Brandenburg et al. (1996) A. Brandenburg, K. Enqvist, and P. Olesen, Phys. Rev. D54, 1291 (1996), eprint astro-ph/9602031.