NORDITA-2023-014, RESCEU-5/23, KEK-TH-2504, MS-TP-23-13

Chiral magnetohydrodynamics with zero total chirality

Axel Brandenburg,^{1, 2, 3, 4} Kohei Kamada,^{5} Kyohei Mukaida,^{6, 7} Kai Schmitz,^{8} and Jennifer Schober^{9}

1Nordita, KTH Royal Institute of Technology and Stockholm University, 10691 Stockholm, Sweden

2The Oskar Klein Centre, Department of Astronomy, Stockholm University, AlbaNova, SE-10691 Stockholm, Sweden 3School of Natural Sciences and Medicine, Ilia State University, 0194 Tbilisi, Georgia 4McWilliams Center for Cosmology and Department of Physics, Carnegie Mellon University, Pittsburgh, Pennsylvania 15213, USA 5Research Center for the Early Universe (RESCEU), Graduate School of Science, The University of Tokyo, Hongo 7-3-1, Bunkyo-ku, Tokyo 113-0033, Japan 6KEK Theory Center, Tsukuba 305-0801, Japan 7Graduate University for Advanced Studies (Sokendai), Tsukuba 305-0801, Japan 8Institute for Theoretical Physics, University of Münster, 48149 Münster, Germany 9Institute of Physics, Laboratory of Astrophysics, École Polytechnique Fédérale de Lausanne (EPFL), 1290 Sauverny, Switzerland (Dated: July 19, 2023)

We study the evolution of magnetic fields coupled with chiral fermion asymmetry in the framework of chiral magnetohydrodynamics with zero initial total chirality. The initial magnetic field has a turbulent spectrum peaking at a certain characteristic scale and is fully helical with positive helicity. The initial chiral chemical potential is spatially uniform and negative. We consider two opposite cases where the ratio of the length scale of the chiral plasma instability (CPI) to the characteristic scale of the turbulence is smaller and larger than unity. These initial conditions might be realized in cosmological models such as certain types of axion inflation. The magnetic field and chiral chemical potential evolve with inverse cascading in such a way that the magnetic helicity and chirality cancel each other at all times. The CPI time scale is found to determine mainly the time when the magnetic helicity spectrum attains negative values at high wave numbers. The turnover time of the energycarrying eddies, on the other hand, determines the time when the peak of the spectrum starts to shift to smaller wave numbers via an inverse cascade. The onset of helicity decay is determined by the time when the chiral magnetic effect becomes efficient at the peak of the initial magnetic energy spectrum, provided the CPI does not grow much. When spin flipping is important, the chiral chemical potential vanishes and the magnetic helicity becomes constant, which leads to a faster increase of the correlation length, as expected from magnetic helicity conservation. This also happens when the initial total chirality is imbalanced. Our findings have important implications for baryogenesis after axion inflation.

I. INTRODUCTION

Relativistic plasmas are described by the evolution equations of chiral magnetohydrodynamics (MHD) [1– 9]. Chirality enters in two distinct ways: first, through a nonvanishing chiral chemical potential, μ̃_{5}, and second, through nonvanishing magnetic helicity density, A · B, where B = ∇ × A is the magnetic field expressed in terms of the vector potential A.

It has been known for some time that fermion chirality can be transferred into magnetic helicity and vice versa through the chiral anomaly [10, 11]. The transfer of fermion chirality to magnetic helicity occurs through an instability [12] known as the chiral plasma instability (CPI) [13]. This instability is the fastest at a specific wave number, whose value depends on the chiral chemical potential. The transfer from magnetic helicity to chiral chemical potential does not involve any instability, but occurs just through a nonvanishing nonlinear source term in the evolution equation for the chiral chemical potential [3, 14, 15]. These differences in the evolutions of the chiral chemical potential and magnetic field can lead to nontrivial dynamics, which has triggered a lot of research [16–18]. Since fermion chirality is tightly related to the baryon and lepton asymmetries at high temperature in the early Universe, their co-evolution with magnetic helicity in the context of cosmology has also extensively studied [19–27].

Previous investigations mostly assumed an initial imbalance between fermion chirality and magnetic helicity. This leads to a conversion of fermion chirality to a maximally helical magnetic field [3]. Also just spatial fluctuations can lead to magnetic field production [28, 29]. In many investigations, however, the initial fermion chirality is nonvanishing while initial magnetic helicity is zero or vice versa. Such chiral asymmetry, which can trigger the CPI, could be generated [30–32] in GUT baryogenesis in the early Universe [33–37] or weak interactions in compact stars [38–42] (see also Ref. [43] and references therein). However, numerical studies on other interesting initial conditions are still lacking, where fermion chirality is exactly opposite to magnetic helicity. Such an initial condition is expected if the chiral symmetry in the fermion sector is only broken through the topological density, ∂_{μ}J_{5}^{μ} = −e^{2}F_{μν} F̃^{μν}/(8π^{2}ℏ^{2}c), or the chiral anomaly [10, 11], with J_{5}^{μ} being the chiral current and e^{2}F_{μν} F̃^{μν}/(8π^{2}ℏ^{2}c) being the topological density. Since the topological density can be written as a total derivative of the magnetic helicity density, the sum of chiral asymmetry and magnetic helicity vanishes when they are generated [44].

Configurations with vanishing total chirality are interesting not only in the context of chiral MHD, but also in particle physics and cosmology. At a high enough temperature realized in the early Universe, the electron Yukawa interaction becomes inefficient for T ≳ 10^{5} GeV [45, 46]. There we find the conservation of the total chirality because of ∂_{μ}J_{e_{R}}^{μ} = −g_{Y}^{2} Y_{μν} Ỹ^{μν}/(16π^{2}ℏ^{2}c) with J_{e_{R}}^{μ} being the right-handed electron current and Y^{μν} being the field strength of the hypercharge gauge field with gauge coupling g_{Y}. For instance, in a certain class of axion inflation, configurations with zero net chirality are generated during inflation [44], which can be the origin of the observed baryon asymmetry of the Universe [47– 49] and it could explain the proposed intergalactic magnetic field; see, however, Ref. [26] for the baryon overproduction problem and Ref. [50] for the too large baryon isocurvature problem. The main purpose of this paper is to perform a full numerical chiral MHD simulation under the initial condition of vanishing total chirality and provide a better understanding of the nonlinear dynamics in this case.

Before we begin our investigations, it is useful to recall the main findings of earlier work where the total chirality was mostly different from zero. Following the work of Ref. [14], who studied a system consisting of the gauge field and the chiral chemical potential, but without fluid velocity fields, and with the initial condition ⟨A · B⟩ ≠ 0, μ̃_{5} = 0, three stages can be identified: (i) exponential decline of the magnetic helicity together with an increase of μ̃_{5}, followed by (ii) a continued decrease of the typical peak wave number k_{p}, while μ̃_{5} stays at its maximum value with ⟨A · B⟩ being essentially zero, and (iii) a phase when all the fermion chirality μ̃_{5} gets transferred back to magnetic helicity. As expected, owing to magnetic helicity conservation, and because the magnetic field from the CPI is maximally helical, the magnetic energy density ⟨B^{2}⟩/2 decays at late times such that ⟨B^{2}⟩ξ_{M} ≈ const, where ξ_{M} ≡ k_{p}^{−1} is the magnetic correlation length. In other words, both ⟨B^{2}⟩ and k_{p} decay in the same fashion, but, unlike the expected t^{−2/3} scaling found previously for helical turbulence [51–54], they find a t^{−1/2} scaling both for ⟨B^{2}⟩ and k_{p}. For sufficiently strong initial magnetic fields, the magnetic Reynolds number can be much larger than unity and the eddy turnover scale much longer than the estimated inverse peak momentum scale, if equipartition between the magnetic fields and fluid velocity fields is established. This suggests that the effect of the fluid velocity cannot be negligible in general.

The earlier analytic study of Ref. [14] was revisited using direct numerical simulations of chiral MHD [15]. At large magnetic Reynolds numbers, the authors found clear evidence for a t^{−2/3} scaling of both ⟨B^{2}⟩ and k_{p} at late times. They also found that the initial evolution is not exponential, as suggested in Ref. [14], but linear in time. However, they only considered the case where the initial fermion chirality was zero. When it is finite and balancing exactly the magnetic helicity, the magnetic field decays in a way similar to the case of a strong, nonhelical field [55], where the decay is governed by the conservation of the Hosking integral [56–58]. This integral describes the strength of magnetic helicity fluctuations on different length scales and has the dimensions of cm^{9} s^{−4}, which implies the scalings ξ_{M} ∝ t^{4/9} and ⟨B^{2}⟩ ∝ t^{−10/9} [56]. The general validity of the Hosking integral was further demonstrated by applying a corresponding analysis to the decay of a nonhelical magnetic field in neutron star crusts [59], where the magnetic field evolution is covered by the Hall effect [60].

Our goal here is to bridge the gap between the two extremes, where the initial chirality is either only in the fermions or only in the magnetic field, and to consider the intermediate case where fermion chirality and magnetic helicity balance to zero, extending the study of the present authors [55]. This is another case where the decay of ⟨B^{2}⟩ and k_{p} are described by a correspondingly adapted Hosking integral of the total chirality. In the following, we therefore refer to the Hosking integral with the chiral chemical potential included as the “adapted” Hosking integral; see Ref. [55] for detail.

As mentioned, our findings on the evolution of the system with vanishing total chirality has a significant impact on the present baryon asymmetry of the Universe. Another goal of the present paper is then to clarify how the non-trivial co-evolution of the magnetic field and fermion chirality affect the model space of axion inflation consistent with the present Universe, which has not been explored before.

We begin, by presenting the basic equations and the mathematical setup of our simulations in Sect. II. We then discuss the parameter dependence of characteristic time scales, consider also the effect of spin flipping, and finally cases where the perfectly vanishing chirality balance is relaxed in Sect. III. Applications to the early Universe are discussed in Sect. IV. Conclusions are presented in Sect. V.

II. CHIRAL MAGNETOHYDRODYNAMICS

A. Chiral magnetic effect

Using Lorentz-Heaviside units, the Ampère-Maxwell equation for the QED-like model in the MHD limit (omitting the displacement current) reads

∇ × B = \frac{1}{c}J.  (1)

The electric current J is the sum of the Ohmic current and the chiral magnetic effect (CME) [61–63],

J = \frac{σ}{c} (cE + u × B) + \frac{e^{2}}{2π^{2}ℏ^{2}c}μ̃_{5}B,  (2)

where we consider the case with μ̃_{5} ≪ (e^{2}/ℏc)k_{B}T. By rewriting cE = −∂A/∂t in the Weyl gauge, e^{2}/4πℏc ≡ α, Eq. (2) is rewritten as

\frac{∂A}{∂t} = \frac{c^{2}}{σ} (μ_{5}B − ∇ × B) + u × B,  (3)

where we defined [6]

μ_{5} ≡ \frac{2α}{πℏc}μ̃_{5}.  (4)

This expression agrees with Eq. (32) of Ref. [6], except for a factor of 2 resulting from our definition μ̃_{5} = (μ̃_{R} − μ̃_{L})/2 in terms of the chemical potentials for right- and left-handed fermions [43]. The additional 4π factor in the numerator of the expression in Ref. [6] is a consequence of their use of cgs units.

B. Model description and basic equations

We perform simulations in a cubic domain of size L^{3} with side lengths L and triply-periodic boundary conditions. The mass in the domain is therefore constant, so the mean density \overline{ρ} is always equal to its initial value ρ_{0} and put to unity in all cases. The lowest wave number in the domain is k_{1} = 2π/L. Using N^{3} mesh points, the largest wave number in the simulations is the Nyquist wave number k_{Ny} = k_{1}N/2.

In the following, we set c = 1, so J = ∇ × B. To include the effects of the cosmic expansion with scale factor a(t) ∝ t^{1/2} in the radiation-dominated era, which we assume to be a spatially flat Friedmann Universe, we use correspondingly scaled quantities and conformal time, η(t) = ∫ dt/a(t), in which the evolution equations of MHD are the same as in the absence of expansion [64]. In order to obtain the physical quantities, we can simply normalize the corresponding comoving quantities with the appropriate powers of the scale factor a. Furthermore, using λ = 3ℏ(2α/πk_{B}T)^{2} and including spin flipping and spatial diffusion, our chiral anomaly equation is

\frac{∂μ_{5}}{∂η} +∇·(μ_{5}u) = \frac{λ}{σ} (J − μ_{5}B)·B+D_{5}∇^{2}μ_{5}−Γμ_{5}, (5)

where D_{5} is an empirical diffusion coefficients for the chiral chemical potential. Here we used the relationship between the chiral chemical potential and the number density,

n_{5} ≡ n_{R} − n_{L} = 2 × \frac{μ̃_{5}}{6ℏ^{3}} (k_{B}T)^{2} = \frac{πμ_{5}}{6αℏ^{2}} (k_{B}T)^{2}, (6)

and used J_{5}^{μ} = (n_{5}, n_{5}u − D_{5}∇n_{5}) for the chiral 4current.

Owing to the chiral anomaly [10, 11], the total chirality is conserved in the absence of spin flipping interaction [3, 6]. It is then convenient to introduce the mean magnetic chirality equivalent as

⟨μ_{M}⟩ ≡ \frac{1}{2}λ⟨A · B⟩,  (7)

so that the conservation law derived from Eqs. (3) and (5) can be stated in the form

μ_{tot} = ⟨μ_{5}⟩ + ⟨μ_{M}⟩ = const.  (8)

We complement Eqs. (3) and (5) by the momentum and continuity equations [6, 7, 65]

\frac{Du}{Dη} = \frac{2}{ρ}∇ · (ρνS) − \frac{1}{4}∇ ln ρ + \frac{u}{3} (∇ · u + u · ∇ ln ρ)
− \frac{u}{ρ} [u · (J × B) + ηJ^{2}] + \frac{3}{4ρ}J × B,  (9)
\frac{∂ ln ρ}{∂η} = −\frac{4}{3} [∇ + (∇ ln ρ)] · u + \frac{1}{ρ} [u · (J × B) + ηJ^{2}],

where D/Dη ≡ ∂/∂η + u · ∇ is the advective derivative, S_{ij} = (∂_{i}u_{j} + ∂_{j}u_{i})/2 − δ_{ij}∇ · u/3 are the components of the rate-of-strain tensor, ν is the viscosity, and p is the pressure, which is assumed to be proportional to the density, i.e., p = ρc_{s}^{2}, with c_{s} = 1/\sqrt{3} being the sound speed for the ultrarelativistic fluid.

For all our simulations, we use the Pencil Code [66], where the relevant equations are readily implemented. We use N^{3} = 1024^{3} mesh points for most of the runs, and N^{3} = 2048^{3} mesh points for one particular run. In a small number of cases, we have included the slope-limited diffusion (SLD) scheme of Ref. [67, 68]. In those cases, SLD acts in addition to the ordinary viscous and diffusive processes stated in the equations above, but prevents the code from crashing during an early more violent phase when the mesh resolution is insufficient to dissipate the energy at high wave numbers. At later times, however, this additional numerical device has little effect. Below, we demonstrate in one case that the solutions with and without SLD yield the same result.

C. Diagnostic quantities

We introduce two characteristic times in our simulations, which are the time scale of the CPI and the magnetic diffusion time,

η_{CPI} = σμ_{50}^{−2} and η_{diff} = σk_{0}^{−2},  (10)

respectively. The ratio (η_{diff}/η_{CPI})^{1/2} = |μ_{50}|/k_{0} characterizes the degree of scale separation between the scales of magnetic helicity and fermion chirality. We also define the turnover time of the energy-carrying eddies, which would determine the onset of turbulent inverse cascading,

η_{turb} = (u_{rms}^{max}k_{0})^{−1} ,  (11)

where u_{rms}^{max} is the maximum value (in time) of the rms velocity.

Next, we introduce several parameters with a dimension of velocity. The nature of the CPI is characterized by the following parameters [7]

v_{λ} = |μ_{50}|/(\overline{ρ}λ)^{1/2} and v_{μ} = |μ_{50}|/σ. (12)

The former represents the ratio of the length scale of the magnetic field at saturation of the CPI to the CPI time scale, while the latter represents the ratio of the length scale of the initial instability to the CPI time scale. The ratio v_{λ}/v_{μ} = σ/(\overline{ρ}λ)^{1/2} characterizes the length of the k^{−2} spectrum that develops if the CPI operates without a strong pre-existing field [7]. In the unbalanced case, u_{rms}^{max} is approximated by v_{λ}. In the present case, however, it does not seem to play any role. Instead, to compute u_{rms}^{max}, we approximate the velocity field by the initial magnetic field such that B_{rms}^{2} ≃ ρ̄u_{rms}^{2}. Using Eqs. (7) and (8), we estimate

B_{rms}^{(0)} ≈ (k_{0}|⟨A · B⟩|)^{1/2} ≈ (\frac{2k_{0}|μ_{50}|}{λ})^{1/2} ,  (13)

which thus defines a new quantity ṽ_{λ} as

ṽ_{λ} ≡ (\frac{2k_{0}|μ_{50}|}{ρ̄λ})^{1/2} (≈ \frac{B_{rms}^{(0)}}{ρ̄^{1/2}} ) .  (14)

A predictive estimate for the turnover time of the energycarrying eddies is thus

η_{λ} = (ṽ_{λ}k_{0})^{−1} = (\frac{ρ̄λ}{2k_{0}^{3}|μ_{50}|})^{1/2} ,  (15)

which is later used to predict the time when the inverse cascade sets in.

In this work, an important diagnostics is the magnetic energy spectrum, E_{M}(k). It is normalized such that ∫ E_{M}(k) dk = ⟨B^{2}⟩/2 ≡ E_{M} where E_{M} is the magnetic energy density^{1}. The kinetic energy spectrum E_{K}(k) is defined similarly, i.e., ∫ E_{K}(k) dk = ⟨ρu^{2}⟩/2 ≡ E_{K}. We also define the magnetic helicity spectrum H_{M}(k), which is normalized such that ∫ H_{M}(k) dk = ⟨A · B⟩. In our simulations, k|H_{M}(k)|/2 approaches E_{M}(k) near the maximum. In fact, the spectra H_{M}(k) and E_{M}(k) satisfy the realizability condition [69],

k|H_{M}(k)|/2 ≤ E_{M}(k).  (16)

1 In terms of the mode function in the polarization basis, A(x, t) ≡ ∫ d^{3}k/(2π)^{3/2} ∑_{λ=±} A_{λ}(k, t)e^{λ}(k)e^{ikx}, E_{M} is given as E_{M}(k) = ∑_{s=±}(k^{4}/4π^{2})|A_{s}(k)|^{2}. We also have H_{M}^{s}(k) = (k^{3}/2π^{2})|A_{s}(k)|^{2} and H_{M}(k) = (k^{3}/2π^{2}) ∑_{s=±} s|A_{s}(k)|^{2}.

When this inequality is saturated for specific wave numbers, we say that the magnetic field is locally fully helical. After some time, the magnetic helicity spectrum is characterized by two subranges, one with positive and one with negative values of H_{M}(k), which are separated by the wave number k_{±}(η), where the sign changes. In addition to the evolution of k_{±}(η), we characterize the spectrum and its evolution by the numbers k_{I}(η) and k_{II}(η), which are the wave numbers of the first positive and second negative peak of H_{M}(k).^{2} The intermediate wave number k_{±}(η) is numerically often better determined than k_{II}(η), especially at early times.

The wave number of the first peak of the spectrum is close to the initial inverse correlation length,

ξ_{M} = E_{M}^{−1} ∫ k^{−1}E_{M}(k) dk.  (17)

In fully helical turbulence, the value of ξ_{M}(η) tends to increase with time in a power law fashion, ξ_{M} ∝ η^{q}, where q = 4/9 in our cases of balanced chirality [55]; see also Sec. II E. Note that in our setup the positive helicity modes always dominate the energy density of the magnetic field, and hence approximately we have ξ_{M} ≃ k_{I}^{−1}. It is convenient to introduce the mean magnetic chirality for the positive helicity modes for k < k_{±} and the negative ones for k > k_{±} as

⟨μ_{M}^{+}⟩ = \frac{λ}{2} ∫_{0}^{k_{±}} H_{M}(k) dk,  (18)
⟨μ_{M}^{−}⟩ = −\frac{λ}{2} ∫_{k_{±}}^{∞} H_{M}(k) dk.  (19)

The conservation law takes then the form

⟨μ_{5}⟩ + ⟨μ_{M}^{+}⟩ − ⟨μ_{M}^{−}⟩ = μ_{tot},  (20)

where μ_{tot} = μ_{50} + μ_{M0} = μ_{50} + μ_{M0}^{+} − μ_{M0}^{−} is given by the initial values.

When we study the effect of spin flipping, we invoke a nonvanishing flipping rate with

Γ = \{ \begin{matrix}Γ_{f0} & for η_{flip} ≤ η ≤ η_{off} \\ 0 & otherwise, \end{matrix}  (21)

where η_{flip} denotes the time when spin flipping is turned on, and in a few cases we allow for a finite value of η_{off}, which denotes the time when spin flipping is later turned off again.

D. Initial conditions

In our numerical experiments, the initial magnetic field is fully helical with positive magnetic helicity and random

2 In the present work, k_{I} ≈ k_{p} ≡ ξ_{M}^{−1}, but the latter is based on the magnetic energy spectrum while the former is based on the magnetic helicity spectrum.

FIG. 1. Visualizations of A · B λ/2 (upper row) and μ_{5} (lower row) on the periphery of the computational domain for Run O at η ≈ 4600 (left) 15,000 (middle), and 46,000 (right).

phases. The initial magnetic energy spectrum is a broken power law

E_{M}(k, η_{0}) ∝ \begin{cases}k^{4} & for k < k_{0}, \\ k^{−5/3} & for k > k_{0}, \end{cases}  (22)

where the initial peak is identified as k_{0} = k_{I}(η_{0}). The IR spectrum is motivated by causality constraints [70], while the UV spectrum is taken as a Kolmogorov-type spectrum. The strength of the magnetic field is adjusted such that the initial magnetic chirality obeys μ_{M0} = −μ_{50} such that μ_{tot} = 0. The chiral chemical potential is initially assumed to be uniformly distributed in space. Its initial value is always negative, i.e., μ_{50} < 0. However, even for an initially uniform chiral chemical potential, there is a specific length scale associated with the value of μ_{5} through the wave number of the most unstable mode of the CPI, k = |μ_{50}|/2. The initial velocity is assumed vanishing in all cases.

E. Theoretical predictions

As was recently shown in Ref. [55], the present case of zero total chirality, where the magnetic helicity is canceled by fermion chirality, is remarkably similar to the case of ordinary MHD without chemical potential and zero magnetic helicity. In both cases, as already alluded to in the introduction, one can define a correlation integral of the total chirality, which is a quantity with dimensions cm^{9} s^{−4} and is dubbed the adapted Hosking integral. The evolution of the system can be explained by the conservation of this quantity. With a self-similar evolution of the magnetic spectrum being assumed, this yields the scalings ξ_{M} ∝ η^{4/9} and ⟨B^{2}⟩ ∝ η^{−10/9} for the typical length scale and the magnetic energy density, respectively [56]. Note that the conservation of the adapted Hosking integral suggests

ξ_{M}^{5}⟨B^{2}⟩^{2} = const, or k_{I}^{−3}E_{M}(k_{I})^{2} = const, (23)

if the magnetic energy density is dominated by the positive helicity mode, which is peaked at k = k_{I}. For the magnetic field with an IR spectrum ∝ k^{4}, as motivated from the causality constraints, the evolution of the magnetic field exhibits inverse cascading.

The big difference between ordinary MHD without helicity on the one hand and chiral MHD with helicity balanced by fermion chirality on the other hand is that in the latter, both the magnetic helicity and the fermion chirality are decaying, which we shall call anomalous chirality cancellation (ACC). In the former, by contrast, the Hosking integral based just on the ordinary magnetic helicity density is conserved. In the latter, contrary to the naively expected exponential decay of fermion chirality due to the CPI in chiral MHD, we actually have a much slower power-law decay proportional to η^{−2/3}, since the magnetic helicity is roughly estimated by ξ_{M}⟨B^{2}⟩, and likewise for |⟨μ_{5}⟩| [55]. Here we have considered the case where the real space realizability condition of magnetic helicity [53], |H_{M}| ≤ 2E_{M}ξ_{M}, is nearly saturated at ACC onset. Once this power law decay of the chirality starts, the CPI rate, ⟨μ_{5}⟩^{2}/σ, decays faster than η^{−1}, which suggests that the CPI does not grow anymore. Hence, the magnetic energy is always dominated by helicity modes of the same sign as the initial ones, which, in our case, are positive helicity modes.

The adapted Hosking integral makes sense only when the communication between the helicity and chirality through the CME becomes effective at the characteristic scale. Therefore we expect that the scaling evolution discussed above starts at the time scale of the CME at the peak scale. With the evolution equation for the magnetic field, equivalent to Eq. (3),

\frac{∂B}{∂η} = \frac{1}{σ} [∇^{2}B + ∇ × (μ_{5}B)] + ∇ × (u × B) (24)

(where second term in the right-hand side represents the CME), we estimate η_{ACC} as the time when the following condition is satisfied:

η_{ACC} ≃ \frac{σ}{μ_{5}(η_{ACC})k_{I}(η_{ACC})}.  (25)

Note that from Eq. (24) we can also confirm that the magnetic field has an instability (the CPI) for one of the two circular polarization modes with k = |μ_{50}|/2 being the most unstable mode. The instability rate is roughly given by μ_{50}^{2}/σ, which determines η_{CPI}.

The evolution of the system is classified into two cases, determined by the comparison between η_{λ} and η_{ACC} estimated by the initial conditions of k_{I} and μ_{5}. For relativistic plasmas with ρ̄ ≃ (π^{2}g_{∗}/30)T^{4}, where g_{∗} is number of the relativistic degrees of freedom and σ ≃ T/(α log α^{−1}) [71, 72], we have η_{ACC} < η_{λ} for k_{0} ≪ |μ_{50}| [more precisely, k_{0} ≪ (ρ̄λ/4σ^{2})|μ_{50}|, which is independent of temperature], and vice versa. For k_{0} ≪ |μ_{50}|, we have the following estimates for the evolution of the system,

1. The system is frozen when η < η_{CPI}.

2. The CPI starts to grow at η ≃ η_{CPI}.

3. If the CPI does not sufficiently amplify the negative helicity modes such that k_{I} is unchanged, the chiral chemical potential starts to decay at η = η_{ACC}(> η_{CPI}) with

η_{ACC} ≃ \frac{σ}{|μ_{50}|k_{0}}  (26)

in a mild way.

4. When η ≃ η_{λ}(> η_{ACC}), the system starts to evolve according to the scaling law found in Ref. [55],

k_{I} ∝ η^{−4/9}, E_{M} ∝ η^{−10/9}, and  (27)
⟨μ_{5}⟩ = −⟨μ_{M}^{+}⟩ + ⟨μ_{M}^{−}⟩ ∝ η^{−2/3}.  (28)

In the case of mild hierarchy^{3}, |μ_{50}| ≳ k_{0}, the case we mainly study in the next section, the assumption of inefficient CPI is guaranteed.

For k_{0} ≫ |μ_{50}|, on the other hand, we expect the following evolution of the system.

1. The system is frozen at η < η_{λ}.

2. The magnetic field evolves according to the inverse cascade at η ≃ η_{λ} in a similar way as the usual inverse cascade for a nonchiral helical magnetic field,

k_{I} ∝ η^{−2/3}, E_{M} ∝ η^{−2/3}, and  (29)
⟨μ_{5}⟩ = −⟨μ_{M}^{+}⟩ + ⟨μ_{M}^{−}⟩ = const,  (30)

since the CME is not effective at k ≃ k_{I} so that the magnetic helicity and chirality are individually conserved. Since k_{I} decays, eventually it becomes smaller than |μ_{5}|, and the system enters the regime similar to the previous case.

3. If the CPI does not sufficiently amplify the negative helicity modes, the CME becomes effective at η ≃ η_{ACC} (> η_{λ}), which is now evaluated as

η_{ACC} ≡ \frac{σ^{3}}{|μ_{50}|^{3}k_{0}^{3}} η_{λ}^{−2} ≃ \frac{2σ^{3}}{ρ̄λμ_{50}^{2}} .  (31)

Here, we have used Eq. (25) and k_{I}(η) = k_{0}(η/η_{λ})^{−2/3}, as well as Eq. (14). When η > η_{ACC} we have the inverse cascade with the conservation of the adapted Hosking integral,

k_{I} ∝ η^{−4/9}, E_{M} ∝ η^{−10/9}, and  (32)
⟨μ_{5}⟩ = −⟨μ_{M}^{+}⟩ + ⟨μ_{M}^{−}⟩ ∝ η^{−2/3}.  (33)

Note that in this case, the CPI would not grow much due to the earlier onset of the chirality decay.

The assumption of inefficient CPI is guaranteed if the mild hierarchy |μ_{5}| ≳ k_{I} at η_{CPI} or η_{ACC} still holds. However, we have |μ_{5}|/k_{I}|_{η_{CPI}} = (2σ^{2}/ρ̄λ)^{1/3}, which is O(10^{2}) for the plasma of the Standard Model (SM) of particle physics with α ≃ 10^{−2} and g_{∗} ≃ 10^{2}. In this case, there is a relatively large hierarchy between the chiral chemical potential and the peak scale of the magnetic field. We then might expect an earlier onset of the chirality decay triggered by the CPI. Namely, we have another possibility of the scaling law as follows,

3 With “mild hierarchy”, we have in mind a modest scale separation |μ_{50}|/k_{0} = O(10).

3.′ After some epoch of the onset of CPI,

η ≃ η_{5dec} = c_{A} η_{CPI}, with c_{A} = O(10),  (34)

the system enters the regime of the ACC. If the conservation of the adapted Hosking integral, with the rapid communication between the chirality and helicity through the CME, governs the evolution of the system, we would still have the scaling laws Eqs. (32) and (33), but the evolution nevertheless would have relatively large uncertainty, because of the exponential instability from the CPI. For later purposes, we keep the ambiguity in the scaling evolution and introduce a scaling index q_{5} such that

⟨μ_{5}⟩ = −⟨μ_{M}^{+}⟩ + ⟨μ_{M}^{−}⟩ ∝ η^{−q_{5}},  (35)

which will be used in Sec. IV.

In Table I, we summarize the characteristic time scales relevant for the evolution of the system.

Some of the features described above will be confirmed by direct numerical simulations in Sect. III. They can have important consequences for baryon production, as will be discussed at the end of the paper.

III. RESULTS

In this section, we show the results of the direct numerical simulation. We first study the case with k_{0} ≪ |μ_{50}| until Sec. III G. In Sec. III H we study the case with k_{0} ≫ |μ_{50}|. Some of our observations will turn out to be consistent with the theoretical prediction discussed in Sec. II E. We will also see some other features that have not been addressed there.

A. Visualization of magnetic and fermion chiralities

We begin by discussing the simulation of Ref. [55] with k_{0} ≪ |μ_{50}|, which we refer to as Run O. In Fig. 1, we present visual impressions of magnetic and fermion chiralities in Run O at different times. We see that the turbulent structures gradually grow in size and the extreme values away from zero decrease as time goes on. Furthermore, μ_{5} and A · B λ/2 have predominantly opposite signs, as expected. Locally, however, there is no correspondence between the two fields. This is because the vanishing total chirality is only a statistical property.

B. Evolution of characteristic scales

As discussed in Ref. [55], it is important to allow for sufficient scale separation between the smallest available wave number k_{1} ≡ 2π/L and the initial wave number of the peak, k_{0}. It is also important that there is enough separation between k_{0} and the initial wave number of

FIG. 2. Magnetic energy (red lines) and kinetic energy (blue lines) spectra for Run O at times η = 460, 4600, and 46,000. The dotted lines denote the earliest outputted time η = 0.3. In the upper part, the two-sided arrows indicate the requirements for scale separation at small and large k to obtain the k_{p} ∝ η^{−4/9} decay and to resolve |μ_{50}|, respectively.

the CPI, |μ_{50}|/2, to confirm distinct features of the evolution of the system. Both, k_{0} and |μ_{50}|/2, in turn, must be much smaller than the largest available wave number k_{Ny} = k_{1}N/2. Sufficient scale separation between k_{1} and k_{0} is particularly important for obtaining the theoretically expected increase of ξ_{M} ∝ η^{4/9} along with the decay of E_{M} ∝ η^{−10/9}, based on the conservation of the Hosking integral adapted to the total chirality. Indeed, in Run O, an optimized balance between the two scale separation requirements has been achieved.

With the start of the simulation, the helical random magnetic field, which is present initially, drives turbulent motions through the Lorentz force. Those motions are in approximate equipartition with the magnetic field at high wave numbers; see Fig. 2, where we compare kinetic and magnetic energy spectra at different times. In this figure, we also mark the two scale separation ratios. This observation supports the estimate of u_{rms}^{max} to ṽ_{λ}; see Eq. (14).

As already discussed in Ref. [55], even though there is vanishing net chirality, ⟨μ_{M}⟩ + ⟨μ_{5}⟩ = 0, there is still some degree of inverse cascading, just like in nonhelical magnetically dominated turbulence [54, 73]. We see this clearly in Fig. 3, where the position of the magnetic peak, k_{I}(η), gradually moves to smaller values. At the same time, the height of the peak decreases, following an approximate power law ∝ k^{β}, with β = 3/2; see Fig. 3. This can be explained by the conservation of the Hosking integral [56, 58]; see also Eq. (23). The exponent β = 3/2 is characteristic of the fact that the net chirality vanishes, even though near the peak itself the field is locally fully helical, as we see from the proximity of k|H_{M}(k)|/2 and E_{M}(k); see Eq. (16).

The newly injected magnetic helicity from the CPI

TABLE I. Relevant time scales defined in this paper.

Time scale Expression  Equation Explanation
η_{CPI}  σμ_{50}^{−2}  Eq. (10) time scale of the CPI
η_{diff}  σk_{0}^{−2}  Eq. (10) magnetic diffusion time
η_{turb}  (u_{rms}^{max}k_{0})^{−1}  Eq. (11) turnover time of the energy-carrying eddies
η_{λ}  (ṽ_{λ}k_{0})^{−1} = [ρ̄λ/(2k_{0}^{3}|μ_{50}|)]^{1/2} Eq. (15) predicted turnover time of the energy-carrying eddies
η_{ACC}  σ/[μ_{5}(η_{ACC})k_{I}(η_{ACC})]  Eq. (25) onset time of the ACC
η_{5dec}  c_{A} η_{CPI}  Eq. (34) onset time of the ACC triggered by the CPI

FIG. 3. Magnetic energy (solid lines) and normalized helicity spectra kH_{M}(k)/2 (dotted lines with red and blue symbols for positive and negative helicity spectra, respectively) for Run O at times η = 0.3, 460, 1500, 4600, 15,000, and 46,000. The peaks k_{I} (peaks of the red curves) and k_{II} (peaks of the blue curves) evolve underneath the envelopes ∝ k^{3/2} and ∝ k^{8/3}, respectively.

leads to a growth of the magnetic field at large wave numbers. It manifests itself mostly through the build-up of negative magnetic helicity at high wave numbers. At some point, we also see a gradual propagation of the secondary peak k_{II} toward smaller k, which has not been addressed in Sec. II E. It lies underneath an envelope with an approximate k^{8/3} slope; see Fig. 3. At present, the exponent 8/3 is just empirical and there is no theory for it. It should be noted, however, that in other cases with a shorter inertial range, we have found larger exponents. Thus, the exponent could also be smaller when the inertial range is larger, i.e., when there is more scale separation and ⟨μ_{5}⟩ξ_{M} is larger.

Another characteristic wave number is k_{±}, where the sign of the spectral magnetic helicity changes. It is used in the definitions of ⟨μ_{M}^{+}⟩ and ⟨μ_{M}^{−}⟩ in Eqs. (18) and (19).

In Fig. 4, we plot the evolution of the characteristic wave numbers k_{I}, k_{±}, and k_{II}. We clearly see the k_{I} ∝ η^{−4/9} decay predicted by the conservation of the Hosking integral adapted to the total chirality [55]. It emerges after a time η_{I}, which is expected to be close to η_{λ} (and

FIG. 4. Comparison of k_{I} (red), k_{±} (orange), and k_{II} (blue) for Run O. The dashed-dotted line indicates η^{−4/9}. The green dashed line shows ⟨μ_{5}⟩ and the green dotted line shows the rms value μ_{5}^{rms}.

FIG. 5. Similarly to Fig. 4, but for Run I. The red dasheddotted line indicates here the η^{−1/3} scaling, which describes the k_{I} scaling better than the η^{−4/9} scaling indicated by the red dotted line. The orange and red dots indicate the crossings of the extrapolated tangents on which the times η_{±}^{(i)}, η_{I}, and η_{±}^{(ii)} are based.

TABLE II. Summary of the runs discussed in this paper. Except for Run P, where λ = 500, we have in all other cases λ = 2 × 10^{4}. Runs A and B below the last horizontal line have nonvanishing net chirality and are discussed at the end of the paper. The asterisk on the value of k_{1} = 0.01 Run J” indicates that the resolution is N^{3} = 2048^{3}, so the Nyquist wave number is here the same as for Run J’ with k_{1} = 0.02.

Run k_{0} −μ_{50} σ^{−1} SLD η_{CPI}  η_{diff} η_{μ_{M}^{+}} q_{I} η_{flip} η_{±}^{(i)} η_{I} η_{±}^{(ii)} η_{μ_{M}^{−}} |⟨μ_{M}^{−}⟩|^{max} k_{±}/k_{I} k_{II}/k_{I} v_{λ} u_{rms}^{max} k_{1}
VI 1 160 5 × 10^{−4} no 0.08 2 × 10^{3} 3 1/3 — 0.3 45 14 7  43  1.4 2.0 1.13 0.095 0.2
V 1 80 5 × 10^{−4} no 0.3 2 × 10^{3} 4 1/3 — 0.5 55 40 4  13  1.5 1.8 0.57 0.050 0.2
IV 1 50 5 × 10^{−4} no 0.8 2 × 10^{3} 6 1/3 — 1.6 50 55 50  5.1  1.5 2.0 0.35 0.050 0.1
III 1 30 5 × 10^{−4} no 2.2 2 × 10^{3} 14 1/3 — 1.7 80 108 150  1.5  1.8 2.3 0.21 0.018 0.2
II+ 1 20 5 × 10^{−4} no  5  2 × 10^{3} 18 1/3 — 6 75 160 200  0.46  1.8 3.1 0.141 0.031 0.1
II 1 20 2 × 10^{−4} no 12.5 5 × 10^{3} 60 1/3 — 9 75 120 200 0.009  4.4 6.7 0.141 0.032 0.1
II− 1 20 1 × 10^{−4} no  25 1 × 10^{4} 160 1/3 — 14 75 140 200 0.003  6.7 10 0.141 0.031 0.1
I  1 10 2 × 10^{−4} no  50 5 × 10^{3} 125 1/3 — 30 80 160 250 0.009  6.7 9.6 0.071 0.030 0.05
O 1 10 2 × 10^{−4} no  50 5 × 10^{3} 125 4/9 — 20 70 120 300 0.008  7.3 9.5 0.071 0.0123 0.02
O’ 1 10 2 × 10^{−4} yes 50 5 × 10^{3} 125 4/9 — 20 110 120 400 0.015  6.4 8.4 0.071 0.0103 0.02
L 1 10 2 × 10^{−4} yes 50 5 × 10^{3} 125 4/9 — 180 400 500 1000 0.027 6.7 8.7 0.071 0.0079 0.01
M 1 7 2 × 10^{−4} yes 102 5 × 10^{3} 165 1/3 — 260 220 800 800 0.006  5.8 7.2 0.049 0.0065 0.01
N 1 5 2 × 10^{−4} yes 200 5 × 10^{3} 235 1/3 — 350 200 800 1000 0.0015 6.3 7.8 0.035 0.0055 0.01
N’ 1 5 2 × 10^{−4} yes 200 5 × 10^{3} 200 4/9 — — 800 3000 15000 0.00004 5.5 2.4 0.035 0.0035 0.005
F  1 5 2 × 10^{−4} yes 200 5 × 10^{3} — 4/9 100 — 250 9000 —  —  30  30 0.035 0.0055 0.01
J  1 5 5 × 10^{−4} no  80 2 × 10^{3} 71 4/9 — 230 300 500 700 0.0003 6.1 7.6 0.035 0.0068 0.01
J” 1 5 5 × 10^{−4} no  80 2 × 10^{3} 71 4/9 — 90 300 500 460 0.0006 6.1 7.6 0.035 0.0071 0.01*
J’ 1 5 5 × 10^{−4} no  80 2 × 10^{3} 76 1/3 — 95 120 500 460 0.0005 5.8 7.7 0.035 0.0070 0.02
P  1 0.1 2 × 10^{−4} no 5 × 10^{5} 5 × 10^{3} 10^{4} 3/5 — – 160 — — 3 × 10^{−9} —  — 0.001 0.0070 0.02
G 0.5 10 2 × 10^{−4} no  50 2 × 10^{4} 200 1/3 — 30 360 360 300 0.044  4.6 6.6 0.071 0.0188 0.05
H 0.2 10 2 × 10^{−4} no  50 1.2 × 10^{5} 375 1/3 — 75 2000 2000 3000 0.42  2.0 2.7 0.071 0.0074 0.02
A 1 10 2 × 10^{−4} no  50 5 × 10^{3} — 4/9 — 8 110 210 —  —  7.3 9.5 0.071 0.0109 0.02
B  1 10 2 × 10^{−4} no  50 5 × 10^{3} — 4/9 — 20 90 120 250  —  7.8 8.9 0.071 0.0137 0.02

also η_{turb}); see Eq. (15). In Run O we find η_{I} ≈ 100.

The evolution of k_{±} and k_{II} can be seen more clearly when the Nyquist wave number is larger. We therefore discuss in Fig. 5 another run, also with N = 1024^{3} mesh points, but now with k_{1} = 0.05 (instead of 0.02), so k_{Ny} = 25.6, which is now a little over five times larger than |μ_{50}|/2 = 5. In Table II, this run is referred to as Run I, which differs from the previously discussed Run O mainly in the value of k_{1}. It also has a shallower scaling of the correlation length, ξ_{M} ∝ k_{I}^{−1} ∝ η^{1/3}, which seems to be an artifact caused by insufficient scale separation, i.e., the value of k_{1} is not sufficiently small. Empirically, we find that if k_{0}/k_{1} ≫ 20, there is an inverse cascade with ξ_{M} ∝ k_{I}^{−1} ∝ η^{4/9}. The parameters η_{I}, η_{±}^{(i)}, and η_{±}^{(ii)}, listed in Table II, are discussed below. We also give here the values of v_{μ} and v_{λ}/v_{μ}, as well as the number of mesh points and the length of the run, η_{max}, which can be useful in assessing the reliability of the numerical results. Run O’ is similar to Run O, except that here, SLD has been added. The two runs are virtually indistinguishable.

The evolution of the peaks of the spectrum can be summarized as follows. (i) After the start of the run, the CPI induces a growth of the negative helicity modes at the secondary peak k_{II}, which stays constant until η = η_{±}^{(i)}, and then starts to decrease with time in a power law fashion, k_{II} ∝ η^{−q_{IIa}} with q_{IIa} ≈ 1 in all cases. (ii) The original large-scale spectrum is unchanged until some time η = η_{I} and starts to decrease via an inverse cascade with k_{I}(η) ∝ η^{−q_{I}}, where q_{I} is expected to be equal to the exponent q = 4/9 found in Ref. [55]. (iii) At time η = η_{±}^{(ii)}, the decay of the secondary peak becomes slower with a smaller index, k_{II} ∝ η^{−q_{IIb}}, with q_{IIb} < q_{IIa} ≈ 1. Those parameters are summarized in Table II.

The plot of characteristic wave numbers k_{I}, k_{±}, and k_{II} in Fig. 5 shows three distinct times η_{±}^{(i)} < ∼ η_{I} < ∼ η_{±}^{(ii)}, where k_{±} begins to decrease first rapidly, at η = η_{±}^{(i)}, and later, at η = η_{±}^{(ii)}, more slowly, approximately like η^{−4/9}, i.e., q_{IIb} ≈ q = 4/9. The decay of k_{II} closely follows that of k_{±}. The decay of k_{I}, on the other hand, does not show the rapid decay phase that we see in k_{±} and k_{II}, but turns directly into the approximate η^{−4/9} decay at η = η_{±}^{(i)}.

FIG. 6. Dependence of η_{±}^{(i)}, η_{I}, and η_{±}^{(ii)} on η_{CPI}. η_{±}^{(i)} shows an approximate η_{CPI}^{2/3} dependence along two branches that are separated by a factor of about 6. η_{I} and η_{±}^{(ii)} are essentially independent of η_{CPI}. The inset shows that η_{±}^{(i)} scales inverse quadratically with |μ_{50}|.

C. Onset of inverse cascading

It is of interest to vary the separation between |μ_{50}|/2 and k_{0} to see the dependence of the relevant characteristic times on these wavenumbers. We have performed simulations for different values and consider runs where we change k_{I} and keep k_{50} fixed, and others where we change k_{50} and keep k_{0} fixed. It both cases, of course, since we want to satisfy ⟨μ_{5}⟩ + ⟨μ_{M}⟩ = const, we need to adjust the amplitude of the initial magnetic field correspondingly. The results are summarized in Table II and plotted in Figs. 6 and 7.

One may presume that η_{±}^{(i)} is roughly estimated by η_{CPI} since the grow of negative helicity modes becomes effective at that time. We see, however, that, while η_{±}^{(i)} decreases quadratically with increasing |μ_{50}|, the dependence on η_{CPI} = σμ_{50}^{−2} is shallower than linear and follows approximately an η_{CPI}^{2/3} scaling; see Fig. 6. Thus, k_{II} starts to decline more rapidly when |μ_{50}| is large, although it is unclear why this exponent is here ≈ 2/3. On the other hand, we see that the five data points with k_{1} = 0.01 (Runs L, M, N, J, and J” with smaller |μ_{50}|) lie on another η_{CPI}^{2/3} line that is shifted upward by a factor of about 6 relative to the runs with larger k_{1}. The reason for this is that for large values of η_{CPI}, it became necessary to decrease the value of k_{1}. This decreased the Nyquist wave number since N remained unchanged, which can cause artifacts in the values of k_{±}. Small values of k_{1} also facilitates the η^{4/9} scaling of ξ_{M} and related length scales; see the comparison between Runs N and N’ in Table II. This shows that η_{±}^{(i)} is currently very sensitive to these restrictions which will be alleviated in future with larger computational power. Nevertheless, there is

FIG. 7. Dependence of η_{I} ≈ η_{±}^{(ii)} and η_{±}^{(i)} on (a) η_{turb} and (b) η_{λ}, as well as (c) the dependence of η_{turb} on η_{λ}.

clearly a trend for an uprise in the dependence of η_{±}^{(i)} on η_{CPI} for large values.

Next, we examine the dependence of η_{I} and η_{±}^{(ii)} on k_{0} and μ_{50}. Figure 6 shows that the time η_{I} of the onset of the decline of k_{I} does not strongly depend on the value of μ_{50}. Likewise, the time η_{±}^{(ii)} when the decay of k_{II} slows down, does not strongly depend on μ_{50}. Again, however, there is an upward shift of data points for the four runs, for which k_{1} = 0.01. As discussed in Sec. II E, we expect that η_{I} is close to η_{turb} and η_{λ}. The upper two panels of Fig. 7 show the dependence of η_{±}^{(i)}, η_{I}, and η_{±}^{(ii)} on η_{turb} and η_{λ}, respectively. From these plots, we estimate that

η_{I} ≈ 1.4 η_{turb} ≈ 2.2 η_{λ}.  (36)

In the lowest panel of Fig. 7, we also show the relation between η_{turb} and η_{λ}, i.e.,

η_{turb} ≈ 1.6 η_{λ},  (37)

which shows the validity of the estimate of u_{rms}^{max} in terms of ṽ_{λ}. Equation (36) is useful for estimating the properties of magnetic field strength and coherence length at later times. Therefore, we conclude that the numerical results support, at least for a moderate scale separation, 1 < |μ_{50}|/k_{0} ≲ O(10), the theoretical estimate for the evolution of the characteristic scales given in Sec. II E

FIG. 8. Comparison of ⟨μ_{M}^{±}⟩ for Run II (red lines), Run J (blue lines), and Run G (orange lines). The times η_{ACC} are marked by the correspondingly colored filled symbol at the crossing points of the extrapolated η^{−2/3} decay law with the initially constant values, indicated by dotted line. The η^{−4/3} decay law ⟨μ_{M}^{−}⟩ is shown as the dashed-dotted line.

with a more accurate determination of the time of the onset of the scaling evolution, Eq. (36).

D. Evolution of ⟨μ_{5}⟩ and ⟨μ_{M}^{±}⟩

We now discuss how the chirality of the system evolves. Using Eqs. (18) and (19), we divide the magnetic helicity into ⟨μ_{M}^{+}⟩ and ⟨μ_{M}^{−}⟩. The typical evolution of ⟨μ_{5}⟩ and ⟨μ_{M}^{±}⟩ is as follows. (i) ⟨μ_{5}⟩ and ⟨μ_{M}^{+}⟩ stay constant until the time η = η_{μ_{M}^{+}}, when the ACC commences exhibiting a power law decay. (ii) ⟨μ_{M}^{−}⟩ grows until the time η = η_{μ_{M}^{−}} and then decays.

As discussed in Sec. II E, the decay of ⟨μ_{5}⟩ and ⟨μ_{M}^{+}⟩ due to the ACC is expected to be like η^{−2/3}. In Fig. 8, we have overplotted the asymptotic η^{−2/3} decay laws of magnetic helicity with results of some of the representative numerical runs (Runs II, J, and G), which clearly shows that the numerical results support the theoretical prediction.

The decay of ⟨μ_{M}^{−}⟩ is faster than that of ⟨μ_{5}⟩ and ⟨μ_{M}^{+}⟩ and follows an approximate η^{−4/3} law, resulting in a decay of the ratio ⟨μ_{M}^{−}⟩/⟨μ_{M}^{+}⟩ ∝ η^{−2/3}. Therefore, unless ⟨μ_{M}^{−}⟩ becomes comparable to ⟨μ_{M}^{+}⟩ when the grow stops, a complete cancellation between ⟨μ_{M}^{−}⟩ and ⟨μ_{M}^{+}⟩ never occurs.

The production of ⟨μ_{M}^{−}⟩ is expected to be a result of the CPI. We now address the question of how much ⟨μ_{M}^{−}⟩ is being produced and what its maximum value depends on. Figure 8 shows that ⟨μ_{M}^{−}⟩ is generally rather small, and at least for μ_{50}/k_{0} ≲ 20 there is always a strong imbalance between |⟨μ_{M}^{+}⟩| and |⟨μ_{M}^{−}⟩|, which never enters a phase with a near-complete cancellation.

To see whether this is related to the value of the conductivity, we compare simulations with different values of σ. It turns out that runs with smaller magnetic diffusivity (σ^{−1} = 10^{−4}) result in an even larger imbalance, while those with a larger diffusivity (σ^{−1} = 5×10^{−4}) have a smaller imbalance; compare Runs II+, II, and II− in Table II.

Before closing this section, let us comment on another trend in the numerical runs we conducted regarding the absence of a near-complete cancellation between ⟨μ_{M}^{−}⟩ and ⟨μ_{M}^{+}⟩. For Runs III–VI, the ratio |⟨μ_{M}^{−}⟩|^{max}/|μ_{50}| becomes rather large. This could be due to the very large scale separation of k_{0} and |μ_{50}|. This suggests a possibility that the CPI completes the cancellation between the magnetic helicity and chirality immediately. However, the positive and negative helicity modes are distributed at separate length scales with the negative ones sitting at higher length scales and the latter receives a stronger magnetic diffusion. Therefore we expect the cancellation not to be complete and that the two helicity modes decay with a power-law decay, not an exponential one, though the scaling index can be different from −2/3. In order to investigate the evolution of the system in such extreme cases, |μ_{50}|/k_{0} ≫ O(10), we need to have a sufficiently large box size to realize the corresponding scale separation. The detailed study is left for future study.

E. Onset of ACC

In Fig. 9, we show the dependence of η_{μ_{M}^{+}} and η_{μ_{M}^{−}} on η_{ACC} = σ/|μ_{50}k_{0}| (for the case k_{0} ≪ |μ_{50}|; see Eq. (26)). It turns out that η_{μ_{M}^{+}} increases with η_{ACC} = σ/|μ_{50}k_{0}| such that

η_{μ_{M}^{+}} ≈ 0.2 η_{ACC} = 0.2 σ/|μ_{50}k_{0}|  (38)

provides a good description to the data, which supports the discussion in Sec. II E, at least for a mild hierarchy |μ_{50}| ≳ k_{0}. Furthermore, η_{μ_{M}^{−}} shows an approximately linear dependence on η_{ACC}. This is reasonable because the CPI becomes ineffective when the ACC onsets such that ⟨μ_{M}^{−}⟩ is no longer amplified by the CPI after that.

F. The scale ratios k_{±}/k_{I} and k_{II}/k_{I}

We also mention another observation in the case with k_{0} ≪ |μ_{50}|. At late times, the scale ratios k_{±}/k_{I} and k_{II}/k_{I} reach values that are approximately constant in time. It is about 10 in the case of Run O, i.e., equal to the initial scale separation, |μ_{50}|/k_{0} = 10. One might have expected the scale ratios to increase with |μ_{50}|/k_{0}. However, in all other cases, this ratio is smaller. Some of this might also be caused by one of the two scale separation constraints not being well enough obeyed, although the counter-intuitive trend remains surprising.

FIG. 9. Dependence of η_{μ_{M}^{+}} and η_{μ_{M}^{−}} on (a) σ/|μ_{50}k_{0}| as well as the geometric means of η_{CPI} and (b) either η_{turb} or (c) η_{λ}.

In Fig. 10, we show the ratios k_{±}/k_{I} and k_{II}/k_{I} versus |μ_{50}|/k_{0} = (η_{diff}/η_{CPI})^{1/2}. The two insets give separately the dependencies on 1/η_{CPI}, showing an η_{CPI}^{−0.4} behavior, and on η_{diff}, with a ∝ η_{diff}^{1/2} behavior. We see that k_{±}/k_{I} and k_{II}/k_{I} decrease both with 1/η_{CPI} and with η_{diff}, giving a combined dependence on just the ratio |μ_{50}|/k_{0}. Thus, we see that, somewhat unexpectedly, large |μ_{50}| and small k_{0} tend to be detrimental to producing large scale ratios.

G. Effect of chirality-flipping

The simulations discussed so far had Γ = 0 and they resulted in a final state where ⟨μ_{5}⟩ and ⟨μ_{M}⟩ vanish at late times. As discussed in the introduction, spin flipping could prematurely lead to a vanishing ⟨μ_{5}⟩, which would imply that the decay of ⟨μ_{M}⟩ would slow down and level off at a value away from zero. To study this quantitatively, we show in Fig. 11 the evolution of ⟨μ_{5}⟩, ⟨μ_{M}⟩, and ⟨μ_{5}⟩+⟨μ_{M}⟩ for Run F with η_{flip} = 100 and Γ = 10^{−2} either for the rest of the run or only until η_{off} = 10^{3} (Run F in Table II).

First, we study a case where spin flipping acts permanently (after η = η_{flip}), which is shown in Fig. 11 as

FIG. 10. Dependence of k_{±}/k_{I} and k_{II}/k_{I} on η_{CPI}, showing an η_{CPI}^{−0.4} behavior. on η_{diff}, showing a ∝ η_{diff}^{1/2} behavior.

solid lines. We see that |⟨μ_{5}⟩| begins to decrease rapidly to zero after η_{flip} = 100. This slows down the decay of ⟨μ_{M}⟩, which then declines at a much smaller rate. Qualitatively similar behaviors are also seen for smaller values of Γ. In all cases, we see that ⟨μ_{5}⟩ + ⟨μ_{M}⟩ evolves away from zero. This is because the total chirality is then no longer conserved. The decay of ⟨μ_{M}⟩ is understood by magnetic diffusion. Thus we expect the decrease to slow down for a larger scale separation between the magnetic diffusion scale and k_{I}.

Next, it is also of interest to study a case where spin flipping acts only for a certain time interval and is then turned off again at η = η_{off}. This case is shown in Fig. 11 as dashed lines. We see that, when Γ = 0 after η_{off} = 10^{3}, the sum ⟨μ_{5}⟩ + ⟨μ_{M}⟩ is strictly constant and away from zero. This is in contrast to the case with permanently nonvanishing Γ, where the sum continues to decrease slowly. The constancy of the total chirality leads to the behavior that ⟨μ_{M}⟩ stops to decline rapidly at a larger value. Furthermore, during that time, some of the magnetic helicity decays due to the magnetic diffusion and is temporarily converted back into fermion chirality through the total chirality conservation; see the small increase of ⟨μ_{5}⟩ with a positive maximum at η ≈ 4000 in Fig. 11. Later, however, this excess fermion chirality gets converted back into magnetic fields, which explains the slight uprise of ⟨μ_{M}⟩ near the end of the simulation. Indeed, this process is similar to the one seen in Refs. [14, 15]. This is natural because after the decay of ⟨μ_{5}⟩ the setup becomes very similar to the ones in these studies.

In Fig. 12, we show η_{±}^{(i)} and η_{I} in the presence of spin flipping. The results suggests that the η^{−4/9} decay changes into the faster η^{−2/3} decay. Spin flipping brings ⟨μ_{5}⟩ close to zero. This process stops or slows down the decline of magnetic helicity, which therefore remains positive. At late times, ⟨μ_{5}⟩, which was originally negative, now becomes positive and settles at a value of around

FIG. 11. Evolution of ⟨μ_{M}⟩ (red), ⟨μ_{5}⟩ (blue), and their sum (black) for Run F with η_{flip} = 100 and Γ = 10^{−2} either for the rest of the run or only until η_{off} = 10^{3}.

FIG. 12. k_{±}^{(i)} and k_{I} for Run F with spin flipping, η_{flip} = 100 and Γ = 10^{−2} for the rest of the run. As in Fig. 4, the green dashed line shows ⟨μ_{5}⟩ and the green dotted line shows μ_{5}^{rms}.

⟨μ_{5}⟩ ≈ k_{1}. This is because at later times the positive chirality induced due to the helicity decay by the magnetic diffusion through the chiral anomaly is balanced by the erasure of the chirality through the CME [14, 15], similar to the baryon asymmetry through the magnetic helicity decay much before the electroweak phase transition [19, 20, 24].

The sign of the final value of ⟨μ_{5}⟩ is determined by the magnetic helicity after the decay of ⟨μ_{5}⟩ due to the onset of spin flipping. In the cases presented above, the sign of the magnetic helicity at the time of the onset of spin flipping was positive and thus the chiral chemical potential at later times was also positive. If the initial magnetic field is weaker and the total chirality being negative (see App. A), the sign of the final value of ⟨μ_{5}⟩ can stay negative.

Our runs show that spin flipping can lead to a signifi-

TABLE III. Empirical values of β for cases with |μ_{50}| < k_{0}. For a given value of λ, the values of v_{A0} followed from the requirement that the total chirality vanishes. The resulting maximum rms velocity u_{rms}^{max} is listed for completeness.

Run −μ_{50} λ  v_{μ}  ṽ_{λ} v_{A0} u_{rms}^{max} SLD β η_{I}
P 0.1 500 2 × 10^{−5} 0.014 0.026 0.008 no 0.33 160
Q 0.1 50 2 × 10^{−5} 0.045 0.081 0.028 no 0.15 50
R 0.1 5 2 × 10^{−5} 0.141 0.257 0.076 yes 0.05 27
S 0.5 500 10^{−4} 0.032 0.057 0.019 no 0.33 70

cant increase of the fraction of the magnetic helicity that can be preserved in spite of the fact that the system has vanishing total chirality. This also reduces the total energy density dissipation of the system. In the absence of spin flipping, both magnetic helicity and chiral chemical potential would approach zero, so there would be no magnetic helicity available for successful baryogenesis. In the real Universe, however, spin flipping due to the electron Yukawa interaction, which really violates the (total) chirality conservation, inevitably acts at T ≲ 10^{2} TeV [45, 46], and hence magnetic helicity survives more or less at the electroweak phase transition.

In Fig. 12, we see an interval between the onset of spin flipping, η = η_{flip} = 10^{2}, and the onset of the η^{−2/3} scaling evolution of k_{I}, η ∼ 6 × 10^{3}, which marks the real onset of the evolution with (pure) magnetic helicity conservation. For a rough estimate of the magnetic field evolution, however, we shall practically use η_{flip} as the switching time between the adapted Hosking integral conservation and the (pure) magnetic helicity conservation.

H. Cases with initially small |μ_{50}|/k_{0}

In all the cases considered so far, we assumed |μ_{50}|/k_{0} > 1. We now consider the opposite case and discuss runs with μ_{50} = −0.1, keeping still k_{0} = 1, so |μ_{50}|/k_{0} = 0.1 (Runs P, Q and R), and also a run with μ_{50} = −0.5 and k_{0} = 1 (Run S). To prevent the magnetic field from being too weak, while still preserving vanishing total chirality, we now decrease the value of λ and choose λ = 500, 50, and 5 for Runs P (and S), Q, and R, respectively. All the runs end at η ∼ 10^{4}. The parameters of these runs are summarized in Table III. The electric conductivity is taken as σ^{−1} = 2 × 10^{−4}. In all our cases, the system does not even reach η_{CPI}, because of the small value of |μ_{50}|, which enters the CPI time inverse quadratically; see Eq. (10).

Smaller values of λ correspond to to larger magnetic fields. We see that this also leads to a gradual decrease of the scaling index of the envelope of the magnetic energy spectrum, β, toward zero. For a given value of β, we expect that the scaling indices q, p, and r, which are those

TABLE IV. Possible combinations of q = 2/(β + 3), p = 2(1 − q), and r = p − q in the range 0 ≤ β ≤ 1/3.

\begin{matrix}β & q & p & r \\ 3/2 = 1.50 & 4/9 ≈ 0.44 & 10/9 ≈ 1.11 & 2/3 ≈ 0.67 \\ 1/3 ≈ 0.33 & 3/5 = 0.60 & 4/5 = 0.80 & 1/5 = 0.20 \\ 0.15 & 0.63 & 0.73 & 0.10 \\ 0.05 & 0.66 & 0.69 & 0.03 \\ 0 & 2/3 = 0.67 & 2/3 = 0.67 & 0 \end{matrix}

FIG. 13. Time dependence of E_{M} (black), ξ_{M}^{−1} (orange), H_{M} (red), and −2⟨μ_{5}⟩/λ (blue), for Runs P (a) and R (b).

of the evolution of the magnetic coherence length, energy density, and helicity, respectively (ξ_{M} ∝ η^{q}, E_{M} ∝ η^{−p}, and H_{M} ∝ η^{−r}), are given as q = 2/(β + 3), p = 2(1 − q), and r = p − q. In Fig. 13(a), we see that for Run P the exponents in agree reasonably well with those expected for β = 1/3. In Fig. 13(b), we also show the results for Run R, where λ is a hundred times smaller and the magnetic field ten times stronger. Now the value of β is very small (about 0.05), corresponding to q = 0.66, p = 0.69, and r = 0.03.

In Table IV, we list several combinations of the expected scaling indices q, p, and r for 0 ≤ β ≤ 3/2. Interestingly, in the range 0 ≤ β ≤ 1/3, the values of q and p do not vary much in this range, especially compared to the case for the evolution with the (adapted) Hosking integral conservation, β = 3/2, so if they do not agree with

FIG. 14. pq diagram for Runs P (red symbols), Q (orange symbols), and R (blue symbols) at times t = 700, 1000, 1500, 2200, 3200, 4600, 6800, 10^{4}, 1.5 × 10^{4}, 1.5 × 10^{4}, 2.2 × 10^{4}, and 3.2 × 10^{4}, corresponding to symbols of increasing size. The solid line denotes the scale-invariance line p = 2(1 − q), the dashed line the empirical β = 1/3 line, and the dasheddotted line is the resulting r = 1/5 line for the magnetic helicity decay. We also show cases with stronger magnetic field strength for Runs Q and R, where the solutions evolve along β ≈ 0.15 and β ≈ 0.05, respectively. Toward the end of the runs, the finite size effects of the domain begin to affect the solution. The dotted line denotes the β = 0 line for magnetic helicity conservation and is shown for comparison.

FIG. 15. Comparison of k_{I} (red), k_{±} (orange), and k_{II} (blue) for Run S with μ_{50} = −0.5. The dashed-dotted line indicates η^{−4/9}. The green dashed line shows ⟨μ_{5}⟩ and the green dotted line shows the rms value μ_{5}^{rms}.

those from the simulations, the discrepancy cannot easily be resolved by changing the value of β within reasonable limits. Note that β = 0 is expected if the evolution is governed by (pure) magnetic helicity conservation.

In the corresponding pq diagram Fig. 14, we see that Run P approaches the scale-invariance line p = 2(1 − q) along the line r = 1/5. At the intersection, we have q = 3/5 and p = 4/5. However, for Runs Q and R with stronger magnetic field strengths, here measured in terms of the initial Alfvén speed v_{A0} = ṽ_{λ} = B_{0}/\sqrt{ρ_{0}}, the solution approaches the β = 0 line, which suggests better conservation of magnetic helicity. Note that near the end of those runs, the data points may not be reliable because of the finite size of the domain. In addition, because of the stronger magnetic field, the Alfvén time is shorter and therefore k_{I} reaches k_{1} more quickly. In any case, it is likely that for small |μ_{50}|/k_{0} we see an intermediate stage of the evolution of the system where the magnetic helicity and chirality are temporarily conserved individually as discussed in Sec. II E. This is supported by the fact that the theoretically predicted time of the onset of ACC (or even CPI) comes much later than near the end of the run; see Eqs. (10) and (31). In other words, from the present simulation results, we cannot distinguish the two possibilities of the theoretical prediction for the later evolution in the case k_{0} > |μ_{50}|; see steps 3 or 3′ given in Sec. II E.

For all the simulations where initially |μ_{50}| < k_{0}, we find that μ_{5} decays more slowly than k_{I}; see Fig. 15 for Run S, as an example, where we see the crossing of μ_{5} and k_{I}. The same is also seen for |μ_{50}| = 0.1, but then the crossing of |μ_{5}| and k_{I} is less prominent. Again, these observations suggest that the magnetic helicity-conserving phase is an intermediate one before the solution resumes the decay governed by the adapted Hosking integral, as discussed in Sec. II E.

The time evolution of the magnetic energy and helicity spectra for Run P are given in Fig. 16(a). We can see that a negative magnetic helicity part of the spectrum still emerges, again only at large wave numbers, although now much later. This means that |⟨μ_{M}^{−}⟩| is induced by the CPI, but it stays extremely small. However, the time is still before η_{CPI} and hence it is likely that we see the early stages of the CPI before its much stronger amplification. Furthermore, |⟨μ_{M}^{+}⟩| does not decay much during the time of the run; see also Fig. 13(b). This can easily be understood by the fact that η_{CPI} becomes very large in this run. Many other features of the magnetic field evolution remain superficially similar to the limit of large |μ_{50}|/k_{0}. One still sees inverse cascading of positive magnetic helicity.

In Fig. 16(b), we show the total helicity variance spectrum Sp(h_{tot}). We clearly see that the spectral slope changes from k^{2} to k^{4} as time goes on. Again, this suggests that the adapted Hosking integral, as defined in Ref. [55], is not conserved.

Looking at the scaling of the magnetic helicity correlation function in Fig. 17, we see that there is a plateau only for very small values of R, which are generally smaller than ξ_{M} (marked by colored symbols). This shows that the conservation of the Hosking integral is now superseded by the conservation of the magnetic helicity. This can be understood by the fact that the CME is inefficient at k_{I} during the time of the run [see Eq. (31)], and hence the magnetic helicity and fermion chirality would be individually relatively well conserved quantities.

Regarding the conserved quantity for runs in the limit

FIG. 16. (a) Magnetic energy and (b) total helicity variance spectra at t = 31 (dashed), 100 (solid), 316 (dotted), 10^{3} (blue), 3.16 × 10^{3} (green), 10^{4} (orange), and 3.16 × 10^{4} (red) for Run P. In (a), note that the E_{M}(k, t) evolve underneath the envelope k^{1/3}, and the upward arrow indicates the sense of time. In (b), the slopes k^{4} and k^{−4} have been indicated and the inset compares Sp(2μ_{5}/λ) (solid) with Sp(h_{tot}) (dotted) at the last time.

of small |μ_{50}|/k_{0}, we can say that, in spite of vanishing total chirality, the Hosking integral is here not conserved, because the magnetic energy now peaks at scales where the CME is not effective during the time of the run and the magnetic helicity is conserved. As a result, the net chirality is no longer random, but systematically of positive sign. The subinertial range of the magnetic helicity variance begins to be dominated by a k^{4} spectrum, which suggests that the Hosking integral in the expansion Sp(h_{tot}) = I_{H}k^{2}/2π^{2} + O(k^{4}) is now subdominant.

To summarize, these runs are consistent with the theoretical prediction in Sec. II E, up to the intermediate stage (step 2), although a moderate violation of helicity conservation has been seen for Run P. Note that Run P has a larger value of λ, which makes the theoretically predicted η_{ACC} smaller [see Eq. (31)], so that an earlier transition to the evolution with adapted Hosking integral conservation is expected. For an analytic estimate of the evolution of the system in the next section, we shall use the theoretical prediction discussed in Sec. II E. Namely, the system is frozen until η = η_{λ} and then evolves with the usual inverse cascade for the helical magnetic field as an intermediate stage. We choose step 3′ for the onset of the ACC, which would be more realistic, though at present we do not have any support from numerical results. That is, at η = η_{5dec} the system starts to evolve

FIG. 17. (a) I_{H}(R, t) versus R for different times t_{∗} (indicated by the same colors/line styles, and (b) I_{H}(R, t) versus t (normalized) for R = ξ_{M}(t_{∗}) marked by the four colors for Run P. In (a), the four colored symbols indicate the positions of k_{0}ξ_{M}(t_{∗}), and in (b), the time dependencies are plotted for those R = ξ_{M}(t_{∗}).

with a decay law determined by the conservation of the adapted Hosking integral.

IV. APPLICATION TO THE EARLY UNIVERSE

A. From QED to the Standard Model

Now we investigate the impact of our findings in the previous sections on the cosmology of the early Universe, especially, baryogenesis. Up to here, we focused on a QED-like theory. Thus, we first would like to clarify its relation to the dynamics in the early Universe. The SM involves the right-handed leptons e_{Rf}, the left-handed lepton doublets ℓ_{Lf}, the right-handed up- and down-type quarks, u_{Rf} and d_{Rf}, and the left-handed quark doublets q_{Lf} with the flavor index running through f = 1, 2, 3, alongside the scalar Higgs doublet Φ, which are in total 16 species. On top of this, we have gauge interactions of U(1)_{Y} × SU(2)_{L} × SU(3)_{C}. It is not obvious why we can reduce this complicated system to chiral MHD based on a QED-like theory like the one introduced in Sec. II.

What we are interested in here is the slow dynamics at long wave lengths compared to interactions among particles. The key idea for the reduction is to assume the equilibration of fast interactions and to keep only the slow variables. The hypermagnetic field of U(1)_{Y} with a correlation length much larger than the mean free path of the particles stands out as a slow variable because the magnetic flux cannot be cut thanks to the absence of monopoles. This feature does not hold for non-Abelian gauge fields because they are charged under their own gauge group. We also need the chiral chemical potential, since it is related to the magnetic field via the anomaly equation. Apart from these two fundamental building blocks, we can coarse-grain the microscopic properties of all particle in the form of transport coefficients such as the diffusion constant and the electric conductivity, besides macroscopic quantities such as the pressure, energy density, and velocity field. In this way, one may see that the system can be reduced to chiral MHD as far as the slow and long-wave dynamics is concerned.

Still, one might wonder why we can just focus on one particular chiral chemical potential, as in Eq. (3), since we have 15 chiral fermion species in the SM. To illustrate this, let us focus on the temperature right above 10^{5} GeV, where the electron Yukawa interaction is not efficient compared to the cosmic expansion, but other interactions are fast enough. In this case, the chiral chemical potential for the right-handed electron, μ̃_{e}, should be counted as a slow variable, as it is directly related to the hypermagnetic field via the anomaly equation. On the other hand, other chiral chemical potentials are subject to fast SM interactions, which provides 11 nontrivial constraints among them. Recalling that the SM has four conserved charges, hypercharge Y and the flavored baryon-minus-lepton numbers B/3 − L_{f} with f = 1, 2, 3, one may immediately see that the remaining 15 chemical potentials can be expressed as a function of μ̃_{e} by solving 11 + 4 constraints. The chiral chemical potential μ̃_{5} originates from the generalized Ohm’s law,

J_{Y} = σ_{Y} E_{Y} + \frac{2α_{Y}}{π} μ̃_{5} B_{Y} ,  (39)

where α_{Y} is now the U(1)_{Y} fine-structure constant and σ_{Y} the hyperelectric conductivity of the plasma. In the following, we will work with the α_{Y} value around the electroweak scale, α_{Y} ≃ 0.01, and neglect its renormalization group running when considering the dynamics of the hypermagnetic field at high energies. Also, note that, in this section, we set ℏ = c = k_{B} = 1, and all quantities are physical rather than comoving, unless explicitly stated otherwise. For the SM U(1)_{Y}, at T ∼ 10^{5···6}GeV, one may express this μ̃_{5} as a summation of chiral chemical potentials for the SM fermions as [12, 49]

μ̃_{5} = ∑_{i=1}^{15} ϵ_{i} g_{i} y_{i}^{2} \frac{μ̃_{i}}{2} = \frac{711}{481} \frac{μ̃_{e}}{2} ,  (40)

where i runs over all SM fermions, ϵ_{i} = ± for right- and left-handed fermions, respectively, g_{i} counts internal degrees of freedom, and y_{i} is the hypercharge of fermion species i. In the second equality, we inserted the solution of the 15 constraint equations mentioned above. We now see that, up to the O(1) coefficient of 711/481, one chiral chemical potential suffices to describe the system. In higher temperature regimes, we will have additionally more slow variables that enter the expression of μ̃_{5}, but it is still written as a linear combination of their chemical potentials with O(1) coefficients. It still holds that the evolution of the system is described by chiral MHD as discussed in Sec. II with μ̃_{5} being evaluated accordingly.

B. Baryogenesis

After these general remarks, let us now turn to the implications of our analysis for the generation of the baryon asymmetry of the Universe. We are primarily interested in the scenario of baryogenesis from decaying hypermagnetic helicity [19, 20, 24–26], which assumes the presence of a strongly helical hypermagnetic field during the radiation-dominated era in the early Universe. This scenario is based on the observation that the helicity stored in the hypermagnetic field decays at the time of the electroweak phase transition, not because of some exotic helicity-violating interactions, but simply because hypermagnetic helicity is converted to magnetic helicity. This decay of hypermagnetic helicity then sources a baryon asymmetry via the chiral anomaly of the baryonnumber current.

One possibility to generate the helical hypermagnetic field required for baryogenesis consists of axion inflation featuring a Chern–Simons coupling to U(1)_{Y}. Such a model leads to the nonperturbative production of hypermagnetic gauge fields in combination with charge asymmetries for the 15 chiral SM fermion species [48, 49],

n_{i} − n̄_{i} = \frac{1}{6} g_{i} μ̃_{i} T^{2} = −ϵ_{i} g_{i} y_{i}^{2} \frac{α_{Y}}{2π} h_{Y} + ⋯, (41)

where the ellipsis represents all other SM contributions, which, however, can safely be neglected during inflation.^{4} Furthermore, h_{Y} in Eq. (41) is the physical helicity density, which we define in terms of the comoving vector potential A_{Y,com}, comoving hypermagnetic field B_{Y,com}, and scale factor a,

h_{Y} = \frac{1}{a^{3}} ⟨A_{Y,com} · B_{Y,com}⟩ ,  (42)

where the angle brackets now stand for a double average including the spatial average and the quantum mechanical expectation value during inflation. From Eq. (41), we can read off the fermion chemical potentials at the end of inflation in terms of the helicity density at the end of inflation. Inserting this result into Eq. (40), we obtain the chiral chemical potential at the end of inflation,

\frac{μ̃_{5}}{T} = −\frac{c_{5}}{2} 6 χ ,  c_{5} = \frac{95}{18} ,  (43)

4 The top-quark Yukawa interaction would be a possible exception; see the discussion in footnote 5 of Ref. [49] for more details.

where the dimensionless yield parameter χ quantifies the amount of CP violation during axion inflation [49],

χ = \frac{α_{Y}}{2π} \frac{h_{Y}}{T^{3}} .  (44)

Here, we assume instantaneous reheating. The same coefficient c_{5} was found in Ref. [49]; in total, the expression for μ̃_{5} in Eq. (43) is, however, smaller than the one in Ref. [49] by a factor of ^{1}/2 because, in the present paper, we include a factor of ^{1}/2 in Eq. (40).

The fermion asymmetries generated during axion inflation are consistent with the chiral anomalies of the respective fermion currents. In fact, it is straightforward to generalize the conversion law in Eq. (8) to the early Universe. To see this, let us rewrite Eq. (43) as follows,

μ_{5} + \frac{3 c_{5}}{2} (\frac{2α_{Y}}{π} )^{2} \frac{1}{2a^{3}T^{2}} ⟨A_{Y,com} · B_{Y,com}⟩ = 0 , (45)

where we used μ_{5} = (2α_{Y} /π) μ̃_{5}. Then, introducing

λ_{Y} = 3 (\frac{2α_{Y}}{πaT} )^{2} ,  (46)

we obtain the relation

μ_{5} + \frac{c_{5}}{2} μ_{M}^{Y} = 0 ,  (47)

where

μ_{M}^{Y} = \frac{1}{2a}λ_{Y} ⟨A_{Y,com} · B_{Y,com}⟩ .  (48)

As the temperature in the early Universe decreases, more and more SM interactions reach chemical equilibrium. This includes the SM Yukawa interactions, which violate parity and hence render the coefficient c_{5} in Eq. (43) a time-dependent quantity [49]. During axion inflation, c_{5} assumes its maximal value, c_{5} = 95/18 ≃ 5.3, before it then decreases down to c_{5} = 711/481 ≃ 1.5 at temperatures of a few 100 TeV [see Eq. (40)]. This change in c_{5} is reflected in a changing value of the chiral chemical potential μ_{5}, which is always given by μ_{5} = −c_{5}/2 μ_{M}^{Y} according to Eq. (47), with μ_{M}^{Y} remaining constant until the onset of CPI, ACC or electroweak phase transition. At T ≲ 10^{5} GeV, c_{5} and hence μ_{5} vanish because all SM interactions have reached chemical equilibrium.

The CP asymmetry parameter χ in Eq. (44) controls the outcome of baryogenesis from helicity decay. That is, if no CPI or ACC takes place before the onset of spin flipping, the decay of hypermagnetic helicity around the electroweak phase transition results in a present-day baryon asymmetry (quantified in terms of the baryon-to-photon ratio) that is fully controlled by χ [49],

η_{B}^{0} ≡ \frac{n_{B}^{0}}{n_{γ}^{0}} ≃ 0.15 c_{B}^{dec}χ ,  (49)

where n_{γ} = 2 ζ(3)T^{3}/π^{2} and the superscript 0 indicates that a quantity is evaluated at the present time. Here, the coefficient c_{B}^{dec} has a theoretical uncertainty of possibly two orders of magnitude [26]. In the following, we will work with the representative value c_{B}^{dec} = 0.05 [49, 50], which implies that χ values of the order of χ ∼ 10^{−7} are necessary to reproduce the observed baryon asymmetry, η_{B}^{obs} ≃ 6.1 × 10^{−10} [74, 75]. Meanwhile, the parameter χ also allows us to evaluate the ratio of k_{0} and μ_{5} at the end of axion inflation. Specifically, if we estimate the comoving peak wave number k_{0} in terms of the comoving wave number that enters the Hubble horizon at the end of reheating, k_{rh} = a_{rh}H_{rh} [49], we find

\frac{|μ_{5}|}{k_{rh}/a_{rh}} = \frac{6 α_{Y} c_{5} χ}{π} \frac{T_{rh}}{H_{rh}} = \frac{6 α_{Y} c_{5} χ}{π} \frac{M_{∗}}{T_{rh}}
∼ 10^{−4} (\frac{χ}{10^{−7}} ) (\frac{10^{14}GeV}{T_{rh}} ) ,  (50)

where M_{∗} = (90/π^{2}/g_{∗})^{1/2} M_{Pl} ≃ 7.1 × 10^{17} GeV is the reduced Planck mass, M_{Pl} ≃ 2.4 × 10^{18} GeV, rescaled by the effective number of relativistic degrees of freedom in the Standard Model plasma, g_{∗} = 427/4. Axion inflation typically results in small values of the χ parameter (e.g., χ ∼ 10^{−7}; see above) and large values of the reheating temperature (e.g., T_{rh} ∼ 10^{14} GeV; see Ref. [49]), which puts us in the parametric regime where |μ_{5}| ≪ k_{rh}/a_{rh} at the end of axion inflation. Moreover, smaller values of T_{rh} typically result in smaller values of χ, following the scaling relation χ ∝ (T_{rh}/M_{∗})^{3} [49], which means that the opposite hierarchy, |μ_{5}| ≫ k_{rh}/a_{rh}, cannot simply be obtained by considering a smaller reheating temperature.

For |μ_{5}| ≪ k_{rh}/a_{rh}, the chiral chemical potential eventually becomes larger than the peak wave number of the hypermagnetic energy spectrum. This is because the peak momentum decreases via the inverse cascade, while the chiral chemical potential is approximately conserved. As mentioned above Eq. (34), we expect a large hierarchy between μ_{5} and k_{I} at η_{CPI} for the SM plasma,

\frac{|μ_{5}|}{k_{I}/a}|_{CPI} ≃ (\frac{2σ_{Y}^{2}}{ρ̄λ_{Y}} )^{1/3}  (51)

∼ 2 × 10^{2} (\frac{106.75}{g_{∗}})^{1/3} (\frac{0.01}{α_{Y}})^{2/3} , (52)

where we evaluate the hyperelectric conductivity σ_{Y} as σ_{Y} = c_{σ_{Y}} T, with c_{σ_{Y}} ∼ 10^{2} [71, 72], the average radiation energy density ρ̄ as ρ̄ = cρ̄ T^{4}, with cρ̄ = π^{2}g_{∗}/30, and the parameter λ_{Y} as λ_{Y} = c_{λ_{Y}} /(aT)^{2} with c_{λ_{Y}} = 12 α_{Y}^{2} /π^{2}. The net chiral chemical potential may start to decay after some duration of CPI as given in Eq. (34)

η_{5dec} = c_{A} η_{CPI} = c_{A} [\frac{1}{a} \frac{σ_{Y}}{μ_{5}^{2}} ]_{rh} ,  (53)

where the factor of a^{−1} follows from the mass dimension of the factor σ_{Y} /μ_{5}^{2}.

The time η_{5dec} marks the onset of the net chirality decay via the CPI and needs to be compared to the time η_{sf} = 1/k_{sf} when spin flipping for left- and right-handed electrons becomes efficient, where k_{sf} = a_{sf}H_{sf} is the comoving horizon scale at η = η_{sf}. Using Eq. (53), together with T_{sf} = 10^{5} GeV for the electron Yukawa interaction in the SM [45, 46], we obtain

\frac{η_{5dec}}{η_{sf}} ∼ (\frac{g_{∗}}{106.75})^{1/2} (\frac{0.01}{α})^{2} (\frac{c_{A}}{10} ) (\frac{10^{−4}}{χ} )^{2} . (54)

If chirality decay occurs before the onset of spin flipping, i.e., χ ≳ 10^{−4}, we expect the hypermagnetic helicity to decay simultaneously. Therefore, the correct amount of baryon asymmetry would be obtained for such a large χ ≫ 10^{−7} that would overproduce baryon (49) at a first glance, since the reduction of the hypermagnetic helicity because of the CPI might counteract this overproduction of baryon number. To see this explicitly, one may use the ratio in Eq. (54) to introduce a dilution factor, Δ = (η_{5dec}/η_{sf})^{q_{5}} with q_{5} being the scaling index introduced in Eq. (35), that multiplies the naive baryon asymmetry in Eq. (49) whenever the CPI should occur before the onset of spin flipping,

η_{B}^{0} ≃ min {1, Δ} × 0.15 c_{B}^{dec}χ .  (55)

By requiring η_{B}^{0} ≃ 6.1 × 10^{−10}, we can estimate the size of the CP violation χ required to obtain the observed baryon-to-photon ratio. For q_{5} = ^{2}/3, we need χ ∼ 10^{6}. Such large χ values are extremely difficult, if not impossible, to realize in realistic models of axion inflation. For q_{5} = ^{4}/3, we need instead χ ∼ 10^{−2}, which is still large but not an unrealistic value for axion inflation. In this way, the result for the final baryon asymmetry strongly depends on the scaling index q_{5}, whose precise value is, however, beyond the scope of our current simulations.

To sum up, even if we initially start from |μ_{5}| ≪ k_{rh}/a_{rh}, the system eventually ends up with the large hierarchy of |μ_{5}| ≫ k_{I}/a. Moreover, if χ is extremely large, we can have initially |μ_{5}| ≫ k_{I}/a. These cases might display similar dynamics as our Runs III–VI, which we commented on at the end of Sec. III D. The decay law for the magnetic helicity may then be different from the η^{−2/3} behavior that we typically find for ACC, which strongly affects the outcome of baryogenesis. As already stated in Section III D, we leave a more detailed study of this more exceptional case for future work.

V. CONCLUSIONS

We have performed numerical simulations of chiral MHD with zero initial total chirality for a range of parameters to determine the dependence of characteristic time and scale ratios, which are well explained by the analytical estimate in Sec. II E. Namely, they are consistent with the scaling evolution, ξ_{M} ∝ η^{4/9}, E_{M} ∝ η^{−10/9}, and ⟨μ_{5}⟩ ∝ η^{−2/3}, derived from the conservation of the adapted Hosking integral [55], and also the time scale of the onset of this scaling evolution, η_{ACC}; see Eqs. (26) and (31). Our numerical simulations also assess the possibility of artifacts resulting from insufficient scale separation. A particularly important constraint is a sufficiently large size of the computational domain (small k_{1}), which is needed to obtain the expected η^{4/9} scaling of the correlation length. When this constraint is not obeyed, the scaling is closer to η^{1/3}. The second constraint of a sufficiently large Nyquist wave number is important to obtain the correct values of the scale ratio of the positive and negative magnetic helicity peaks, i.e., k_{II}/k_{I}. Somewhat surprisingly, this ratio scales inversely with the initial scale separation between the scale of the magnetic field and the CPI scale. Increased values of ⟨μ_{M}^{−}⟩, which characterize the strength of the CPI, are obtained when σ is small or |μ_{50}| is large and therefore the coupling to the CPI is more efficient.

In the absence of spin flipping, even the slightest initial imbalance will amplify as the magnetic energy decays; see Appendix A. On long time scales, this eventually leads to a fully helical state, although simulations of this are at present unable to demonstrate this conclusively owing to the finite size of the computational domain. Spin flipping is another mechanism that can produce an imbalance between magnetic helicity and fermion chirality. In any case, however, the finally available magnetic energy and helicity densities are always limited by the finiteness of the initial total chirality imbalance. For η < η_{ACC}, when the chiral magnetic effect is not effective at the peak scale, magnetic helicity conservation governs the decay of magnetic energy and the Hosking integral does not play a role.

We also discussed the implications of our findings for the generation of the baryon asymmetry of the Universe, in particular, the scenario of baryogenesis from helicity decay. The final baryon asymmetry in this scenario is controlled by a dimensionless yield parameter χ that quantifies the helicity density produced in the very early Universe, for instance, during a stage of axion inflation. In previous work, it was shown how the observed baryon asymmetry can be generated from helicity decay at the time of the electroweak phase transition for a specific χ value, χ_{0} ∼ 10^{−7}; see Eq. (49) and Ref. [49]. The situation at larger χ values, however, remained unclear. At χ ≫ 1, one may have anticipated either (A) the overproduction of baryon number or (B) catastrophic helicity erasure by the chiral plasma instability and consequently no baryon asymmetry at all.

Thanks to the analysis in this paper, we now extend the previous analysis to the case where the decay of the hypermagnetic helicity occurs before the spin flipping for electrons becomes efficient, i.e., χ ≳ 10^{−4}. For typical parameters for axion inflation, although we initially have |μ_{5}| ≪ k_{rh}/a_{rh}, the opposite large hierarchy is eventually realized |μ_{5}| ≫ k_{I}/a_{I}. Therefore the required value of χ to reproduce the observed baryon asymmetry depends on the value of the scaling index q_{5} in the large hierarchy regime. Further studies for the large hierarchy are indispensable for understanding the outcome of baryogenesis.

Data availability—The source code used for the simulations of this study, the Pencil Code, is freely available from Ref. [66]. The simulation setups and the corresponding data are freely available from Ref. [76].

ACKNOWLEDGMENTS

We thank V. Domcke for fruitful discussions and J. Warnecke for his work on the implementation of SLD, which is used in some of the simulations. Support through the grant 2019-04234 from the Swedish Research Council (Vetenskapsrådet) (AB), Grant-in-Aid for Scientific Research No. (C) JP19K03842 from the JSPS KAKENHI (KK), MEXT Leading Initiative for Excellent Young Researchers No. JPMXS0320200430 (KM), Grant-in-Aid for Young Scientists No. JP22K14044 from the JSPS KAKENHI (KM), and the grant 185863 from the Swiss National Science Foundation (JS) are gratefully acknowledged. We acknowledge the allocation of computing resources provided by the Swedish National Infrastructure for Computing (SNIC) at the PDC Center for High Performance Computing Stockholm and Linköping.

Appendix A: Behavior in imbalanced chirality decay

In Sect. V, we emphasized that even the slightest initial imbalance between magnetic helicity and fermion chirality will amplify as the magnetic energy decays. It is therefore important to remember that the dynamics discussed in this paper is specific to the case of balanced chirality, which is arguably also the most generic case. We know that the decay of magnetic energy and the increase of the correlation length follow a different behavior in the completely imbalanced case compared to the unbalanced one. We now discuss the behavior for the mildly imbalanced case. Here, we show that there is a tendency for the system to approach the behavior of a completely imbalanced one.

We discuss two runs, Run A where the initial ⟨μ_{M}⟩ is enhanced by 20% compared with |⟨μ_{5}⟩|, and Run B where it is decreased by 20%. Apart from that, the runs are the same as Run O, i.e., the run discussed in Ref. [55]. In Run A, where the magnetic helicity is weaker than in Run O, the CPI becomes dominant and overcompensates the magnetic helicity. The net chirality is then negative. Eventually, the sign of the magnetic helicity changes and all the remaining fermion chirality is converted to magnetic fields with negative helicity; see Fig. 18, where we show the magnetic energy E_{M}(k, t) and the normalized magnetic helicity spectra kH_{M}(k)/2 for Run A at times η = 32, 320, 1000, 3200, 10,000, and 32,000. We see that k|H_{M}(k)|/2 approaches E_{M}(k) near the maximum. In view of the spectral realizability condition, Eq. (16), this means that the magnetic field is fully helical. Away

FIG. 18. Magnetic energy (solid lines) and normalized helicity spectra kH_{M}(k)/2 (dotted lines with red and blue symbols for positive and negative helicity spectra, respectively) for Run A at times η = 32, 100, 320, 1000, 3200, 15,000, 22,000, and 32,000.

FIG. 19. Same as Fig. 3, but for Run B at times η = 32, 100, 320, 1000, 3200, 10,000, and 32,000.

from the maxima, the inequality is no longer saturated, but this is a typical effect in all turbulent flows where the current helicity spectrum shows a Kolmogorov-type spectrum, making the magnetic helicity spectrum therefore steeper than what could still be allowed by the spectral realizability condition [77].

On the other hand, when the fermion chirality is weak (Run B), the usual inverse magnetic cascade quickly gets established; see Fig. 19. In either case, the fermion chirality gets ultimately converted into magnetic helicity. It just takes a little longer than when the magnetic helicity is initially weak. At the end, however, the usual inverse cascade for a fully helical magnetic field commences. The sign of magnetic helicity can be positive or negative, depending on the sign of the initial total chirality.

To illustrate how the decay laws change when magnetic helicity and fermion chirality no longer balance, we plot in Fig. 20 the time dependencies of E_{M}, ξ_{M}, H_{M}, and −2⟨μ_{5}⟩/λ, for (a) Run A with 20% smaller and (b)

FIG. 20. Time dependence of E_{M} (black), ξ_{M}^{−1} (orange), H_{M} (red), and −2⟨μ_{5}⟩/λ (blue), for (a) Run A with smaller and (b) Run B with larger magnetic helicity than in the balanced case. Dashed lines indicate negative values; at late times −2⟨μ_{5}⟩/λ changes sign in (a), and H_{M} changes sign in (b). In (a) the dotted line denotes the η^{−1/2} log(η/η_{log}) scaling of Ref. [15] with η_{log} = 3.

Run B with 20% larger magnetic helicity than in the balanced case. In both cases, we see a tendency of the decays of E_{M} and ξ_{M}^{−1} to slow down while those of H_{M} and −2⟨μ_{5}⟩/λ follow separate evolutions. Especially in the case of Run B, where the magnetic helicity dominates of the fermion chirality, we see a tendency toward a E_{M} ∝ ξ^{−1} ∝ t^{−2/3} as well as H_{M} = const evolution, as expected from magnetic helicity conservation.

(2009), arXiv:0906.5044 [hep-th]. (2011), arXiv:1011.5107 [hep-th].

[1] D. T. Son and P. Surowka, Phys. Rev. Lett. 103, 191601 [2] Y. Neiman and Y. Oz, J. High. Energy Phys. 03, 023 [3] A. Boyarsky, J. Fröhlich, and O. Ruchayskiy, Phys. Rev.

Lett. 108, 031301 (2012), arXiv:1109.3350 [astro-ph.CO].

[4] A. Boyarsky, J. Fröhlich, and O. Ruchayskiy, Phys. Rev.

D 92, 043004 (2015), arXiv:1504.04854 [hep-ph].

[5] P. Pavlović, N. Leite, and G. Sigl, Phys. Rev. D 96,

023504 (2017), arXiv:1612.07382 [astro-ph.CO].

[6] I. Rogachevskii, O. Ruchayskiy, A. Boyarsky, J. Fröhlich, N. Kleeorin, A. Brandenburg, and J. Schober, Astrophys. J. 846, 153 (2017), arXiv:1705.00378 [physics.plasm-ph].

[7] A. Brandenburg, J. Schober, I. Rogachevskii, T. Kahniashvili, A. Boyarsky, J. Fröhlich, O. Ruchayskiy, and N. Kleeorin, Astrophys. J. 845, L21 (2017).

[8] J. Schober, I. Rogachevskii, A. Brandenburg, A. Boyarsky, J. Fröhlich, O. Ruchayskiy, and N. Kleeorin, Astrophys. J. 858, 124 (2018).

[9] L. Del Zanna and N. Bucciantini, Mon. Not. R. Astron. Soc. 479, 657 (2018), arXiv:1806.07114 [astro-ph.HE].

[10] S. L. Adler, Phys. Rev. 177, 2426 (1969).

[11] J. S. Bell and R. Jackiw, Nuovo Cim. A 60, 47 (1969).

[12] M. Joyce and M. Shaposhnikov, Phys. Rev. Lett. 79, 1193 (1997).

[13] Y. Akamatsu and N. Yamamoto, Phys. Rev. Lett. 111, 052002 (2013), arXiv:1302.2125 [nucl-th].

[14] Y. Hirono, D. E. Kharzeev, and Y. Yin, Phys. Rev. D 92, 125031 (2015), arXiv:1509.07790 [hep-th].

[15] J. Schober, T. Fujita, and R. Durrer, Phys. Rev. D 101, 103028 (2020), arXiv:2002.09501 [physics.plasm-ph].

[16] C. Manuel and J. M. Torres-Rincon, Phys. Rev. D 92, 074018 (2015), arXiv:1501.07608 [hep-ph].

[17] V. Domcke, Y. Ema, K. Mukaida, and R. Sato, J. High Energy Phys. 2019, 111 (2019), arXiv:1812.08021 [hep- -ph].

[18] V. Domcke, Y. Ema, and K. Mukaida, J. High. Energy Phys. 2020, 55 (2020), arXiv:1910.01205 [hep-ph].

[19] M. Giovannini and M. E. Shaposhnikov, Phys. Rev. D 57, 2186 (1998), arXiv:hep-ph/9710234.

[20] M. Giovannini and M. E. Shaposhnikov, Phys. Rev. Lett. 80, 22 (1998), arXiv:hep-ph/9708303.

[21] M. Dvornikov and V. B. Semikoz, J. Cosmol. Astropart. Phys. 02, 040 (2012), [Erratum: J. Cosmol. Astropart. Phys.,08, E01 (2012)], arXiv:1111.6876 [hep-ph].

[22] H. Tashiro, T. Vachaspati, and A. Vilenkin, Phys. Rev. D 86, 105033 (2012), arXiv:1206.5549 [astro-ph.CO].

[23] A. J. Long, E. Sabancilar, and T. Vachaspati, J. Cosmol. Astropart. Phys. 02, 036 (2014), arXiv:1309.2315 [astro-ph.CO].

[24] T. Fujita and K. Kamada, Phys. Rev. D 93, 083520 (2016), arXiv:1602.02109 [hep-ph].

[25] K. Kamada and A. J. Long, Phys. Rev. D 94, 063501 (2016), arXiv:1606.08891 [astro-ph.CO].

[26] K. Kamada and A. J. Long, Phys. Rev. D 94, 123509 (2016), arXiv:1610.03074 [hep-ph].

[27] P. Pavlović, N. Leite, and G. Sigl, Phys. Rev. D 96, 023504 (2017), arXiv:1612.07382 [astro-ph.CO].

[28] J. Schober, I. Rogachevskii, and A. Brandenburg, Phys. Rev. Lett. 128, 065002 (2022), arXiv:2107.12945 [physics.plasm-ph].

[29] J. Schober, I. Rogachevskii, and A. Brandenburg, Phys. Rev. D 105, 043507 (2022), arXiv:2107.13028 [physics.plasm-ph].

[30] J. García-Bellido, D. Grigoriev, A. Kusenko, and M. Shaposhnikov, Phys. Rev. D 60, 123504 (1999), arXiv:hep-ph/9902449 [hep-ph].

[31] K. Kamada, Phys. Rev. D 97, 103506 (2018), arXiv:1802.03055 [hep-ph].

[32] V. Domcke, K. Kamada, K. Mukaida, K. Schmitz, and M. Yamada, Phys. Rev. Lett. 126, 201802 (2021), arXiv:2011.09347 [hep-ph].

[33] M. Yoshimura, Phys. Rev. Lett. 41, 281 (1978), [Erra-

tum: ibid. 42, 746 (1979)].

[34] S. Dimopoulos and L. Susskind, Phys. Rev. D 18, 4500 (1978).

[35] D. Toussaint, S. Treiman, F. Wilczek, and A. Zee, Phys. Rev. D 19, 1036 (1979).

[36] S. Weinberg, Phys. Rev. Lett. 42, 850 (1979).

[37] S. M. Barr, G. Segre, and H. Weldon, Phys. Rev. D 20, 2494 (1979).

[38] A. Ohnishi and N. Yamamoto, (2014), arXiv:1402.4760 [astro-ph.HE].

[39] D. Grabowska, D. B. Kaplan, and S. Reddy, Phys. Rev. D 91, 085035 (2015), arXiv:1409.3602 [hep-ph].

[40] Y. Masada, K. Kotake, T. Takiwaki, and N. Yamamoto, Phys. Rev. D 98, 083018 (2018), arXiv:1805.10419 [astro-ph.HE].

[41] M. Dvornikov, V. B. Semikoz, and D. D. Sokoloff, Phys. Rev. D 101, 083009 (2020), arXiv:2001.08139 [astro-ph.HE].

[42] J. Matsumoto, N. Yamamoto, and D.-L. Yang, Phys. Rev. D 105, 123029 (2022), arXiv:2202.09205 [astro-ph.HE].

[43] K. Kamada, N. Yamamoto, and D.-L. Yang, Prog. Part. Nucl. Phys. 129, 104016 (2023), arXiv:2207.09184 [astro-ph.CO].

[44] V. Domcke and K. Mukaida, J. Cosmol. Astropart. Phys. 2018, 020 (2018).

[45] B. A. Campbell, S. Davidson, J. R. Ellis, and K. A. Olive, Phys. Lett. B 297, 118 (1992), arXiv:hep-ph/9302221.

[46] D. Bödeker and D. Schröder, J. Cosmol. Astropart. Phys. 05, 010 (2019), arXiv:1902.07220 [hep-ph].

[47] D. Jiménez, K. Kamada, K. Schmitz, and X.- J. Xu, J. Cosmol. Astropart. Phys. 12, 011 (2017), arXiv:1707.07943 [hep-ph].

[48] V. Domcke, B. von Harling, E. Morgante, and K. Mukaida, J. Cosmol. Astropart. Phys. 2019, 032 (2019), arXiv:1905.13318 [hep-ph].

[49] V. Domcke, K. Kamada, K. Mukaida, K. Schmitz, and M. Yamada, J. High. Energy Phys. 01, 053 (2023), arXiv:2210.06412 [hep-ph].

[50] K. Kamada, F. Uchida, and J. Yokoyama, J. Cosmol. Astropart. Phys. 04, 034 (2021), arXiv:2012.14435 [astro-ph.CO].

[51] T. Hatori, J. Phys. Soc. Jpn 53, 2539 (1984).

[52] D. Biskamp and W.-C. Müller, Phys. Rev. Lett. 83, 2195 (1999), arXiv:physics/9903028 [physics.flu-dyn].

[53] T. Kahniashvili, A. G. Tevzadze, A. Brandenburg, and A. Neronov, Phys. Rev. D 87, 083007 (2013), arXiv:1212.0596 [astro-ph.CO].

[54] A. Brandenburg and T. Kahniashvili, Phys. Rev. Lett. 118, 055102 (2017), arXiv:1607.01360 [physics.flu-dyn].

[55] A. Brandenburg, K. Kamada, and J. Schober, Phys. Rev. Research 5, L022028 (2023), arXiv:2302.00512 [physics.plasm-ph].

[56] D. N. Hosking and A. A. Schekochihin, Phys. Rev. X 11, 041005 (2021), arXiv:2012.01393 [physics.flu-dyn].

[57] D. N. Hosking and A. A. Schekochihin, (2022), arXiv:2203.03573 [astro-ph.CO].

[58] H. Zhou, R. Sharma, and A. Brandenburg, J. Plasma Phys. 88, 905880602 (2022), arXiv:2206.07513 [physics.plasm-ph].

[59] A. Brandenburg, J. Plasma Phys. 89, 175890101 (2023), arXiv:2211.14197 [physics.plasm-ph].

[60] P. Goldreich and A. Reisenegger, Astrophys. J. 395, 250

(1992).

[61] A. Vilenkin, Phys. Rev. D 22, 3080 (1980).

[62] A. Yu. Alekseev, V. V. Cheianov, and J. Fröhlich, Phys. Rev. Lett. 81, 3503 (1998), arXiv:cond-mat/9803346 [cond-mat].

[63] K. Fukushima, D. E. Kharzeev, and H. J. Warringa, Phys. Rev. D78, 074033 (2008), arXiv:0808.3382 [hep- -ph].

[64] A. Brandenburg, K. Enqvist, and P. Olesen, Phys. Rev. D 54, 1291 (1996), arXiv:astro-ph/9602031 [astro-ph].

[65] A. Brandenburg, Y. He, T. Kahniashvili, M. Rheinhardt, and J. Schober, Astrophys. J. 911, 110 (2021), arXiv:2101.08178 [astro-ph.CO].

[66] Pencil Code Collaboration, A. Brandenburg, A. Johansen, P. Bourdin, W. Dobler, W. Lyra, M. Rheinhardt, S. Bingert, N. Haugen, A. Mee, F. Gent, N. Babkovskaia, C.-C. Yang, T. Heinemann, B. Dintrans, D. Mitra, S. Candelaresi, J. Warnecke, P. Käpylä, A. Schreiber, P. Chatterjee, M. Käpylä, X.-Y. Li, J. Krüger, J. Aarnes, G. Sarson, J. Oishi, J. Schober, R. Plasson, C. Sandin, E. Karchniwy, L. Rodrigues, A. Hubbard, G. Guerrero, A. Snodin, I. Losada, J. Pekkilä, and C. Qian, J. Open Source Softw. 6, 2807 (2021).

[67] M. Rempel, Astrophys. J. 789, 132 (2014),

arXiv:1405.6814 [astro-ph.SR].

[68] See Fig. 7 of the Pencil Code Newsletter 2022/2;

                                                        https://github.com/pencil-code/website.

[69] H. K. Moffatt, Magnetic field generation in electrically conducting fluids (1978).

[70] R. Durrer and C. Caprini, J. Cosmol. Astropart. Phys. 2003, 010 (2003), arXiv:astro-ph/0305059 [astro-ph].

[71] G. Baym and H. Heiselberg, Phys. Rev. D 56, 5254 (1997), arXiv:astro-ph/9704214.

[72] P. B. Arnold, G. D. Moore, and L. G. Yaffe, J. High. Energy Phys. 11, 001 (2000), arXiv:hep-ph/0010177.

[73] A. Brandenburg, T. Kahniashvili, and A. G. Tevzadze, PhRvL 114, 075001 (2015), arXiv:1404.2238.

[74] N. Aghanim et al. (Planck), Astron. Astrophys. 641, A6 (2020), [Erratum: Astron.Astrophys. 652, C4 (2021)], arXiv:1807.06209 [astro-ph.CO].

[75] R. L. Workman et al. (Particle Data Group), Prog. Theor. Exp. Phys. 2022, 083C01 (2022).

[76] See  the  Supplemental  Material  on

                                                        http://norlx65.nordita.org/~brandenb/projects/0Chir/

for data sets and input files for runs with the Pencil Code.

[77] A. Brandenburg and K. Subramanian, Astron. Astrophys. 439, 835 (2005), astro-ph/0504222.