NORDITA-2023-0XX

Relic Gravitational Waves from the Chiral Plasma Instability in the Standard Cosmological Model

Axel Brandenburg,^{1, 2, 3, 4}^{,} ^{∗} Emma Clarke,^{4}^{,} ^{†} Tina Kahniashvili,^{4, 3, 5}^{,} ^{‡} Andrew J. Long,^{6}^{,} ^{§} and Guotong Sun^{4}^{,} ^{¶}

1Nordita, KTH Royal Institute of Technology and Stockholm University, 10691 Stockholm, Sweden 2Department of Astronomy, AlbaNova University Center, Stockholm University, 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, PA 15213, USA 5Abastumani Astrophysical Observatory, Tbilisi, GE-0179, Georgia

6Department of Physics and Astronomy, Rice University, 6100 Main St., Houston, TX 77005, USA (Dated: July 4, 2023)

In the primordial plasma, at temperatures above the scale of electroweak symmetry breaking, the presence of chiral asymmetries is expected to induce the development of helical hypermagnetic fields through the phenomenon of chiral plasma instability. It results in magnetohydrodynamic turbulence due to the high conductivity and low viscosity and sources gravitational waves that survive in the universe today as a stochastic polarized gravitational wave background. In this article, we show that this scenario only relies on Standard Model physics, and therefore the observable signatures, namely the relic magnetic field and gravitational background, are linked to a single parameter controlling the initial chiral asymmetry. We estimate the magnetic field and gravitational wave spectra, and validate these estimates with 3D numerical simulations.

I. INTRODUCTION

The excess of matter over antimatter on cosmological scales in the universe today is well measured but its origin is not yet established. In studies of early universe cosmology, it is typically assumed that the matter-antimatter asymmetry arose dynamically in the first fractions of a second after the Big Bang through a process called baryogenesis [1]. In addition to creating the baryon asymmetry, e.g., the excess of nuclei over antinuclei, baryogenesis may have created other (possibly unstable) particle asymmetries as well; a few examples include lepton asymmetry [2], Higgs asymmetry [3], neutrino asymmetry [4, 5], and right-chiral electron asymmetry [6]. Some of these are examples of chiral asymmetries, n_{5} = n_{R} − n_{L}, namely an excess (or deficit) of right-chiral particles and antiparticles over their left-chiral partners. A particular linear combination of various particle asymmetries, which we call the hypercharge-weighted chiral asymmetry, has attracted interest because of its connections with primordial magnetogenesis [7] through a phenomenon known as the chiral plasma instability [8].

The primordial magnetic field may survive in the universe today as an intergalactic magnetic field, thereby opening a pathway to test this scenario [9, 10]. In addition, the primordial magnetic field and its interaction with the turbulent plasma are expected to source gravitational radiation; see Ref. [11] for pioneering work and

∗Electronic address: brandenb@nordita.org †Electronic address: emmaclar@andrew.cmu.edu ‡Electronic address: tinatin@andrew.cmu.edu §Electronic address: andrewjlong@rice.edu ¶Electronic address: guotongs@andrew.cmu.edu

Ref. [12] for numerical simulations of the gravitational waves induced by the primordial magnetic field originating from the chiral plasma instability. In this work, we investigate the gravitational wave signatures of a primordial hypercharge-weighted chiral asymmetry via the chiral plasma instability.

Contrary to earlier numerical simulations, we study here a parameter regime that is more realistic in various respects. The resulting gravitational wave energy from our simulations confirm the scaling with the sixth power of the chiral chemical potential and the fifth power of the inverse square root of the chiral dilution parameter, which can be combined into a single parameter, as already found previously [12].

II. DESCRIPTION OF THE MODEL

We consider the primordial Standard Model plasma at temperatures T > ∼ 100 TeV in the phase of unbroken electroweak symmetry. We remain agnostic as to the physics of baryogenesis, but assume that a nonzero hyperchargeweighted chiral asymmetry is present in the plasma initially. We study the growth of an initially vanishingly small hypermagnetic field via the chiral plasma instability and calculate the resulting gravitational wave radiation. The hypermagnetic field generated by the chiral plasma instability is always maximally helical and therefore also leads to the production of maximally circularly polarized gravitational waves. The present work is conceptually different from that of Refs. [13, 14], where a helical magnetic field was present initially such that the net chirality of the system was balanced to zero by a fermion chirality of opposite sign.

One appealing aspect of our approach is its minimalism: we only assume Standard Model particle physics and the standard cosmological model after reheating. Our only free parameter is the initial hypercharge-weighted chiral asymmetry, which presumably arises from physics beyond the Standard Model. We work in the Lorentz-Heaviside unit system with h̄ = c = k_{B} = 1. We account for the cosmological expansion using an Friedmann-Lemaître-Robertson-Walker metric with dimensionless scale factor a(t) and set a_{0} = 1 today. Unless otherwise specified, all dimensionful variables are comoving; this includes conformal time dt = dt_{phys}/a, comoving magnetic field B = a^{2}B_{phys}, comoving magnetic correlation length ξ_{M} = ξ_{M,phys}/a, comoving temperature T = aT_{phys}, comoving wavenumber k = ak_{phys}, comoving Hubble parameter H = aH_{phys} (with H ≡ (da/dt)/a), comoving energy density of any relativistic component (including frozen-in magnetic fields, gravitational waves, etc) E = E_{phys}a^{4} for the relativistic plasma, and so on. We denote Newton’s gravitational constant by G, the Planck mass by M_{Pl} = 1/\sqrt{G} = 1.2 × 10^{19} GeV, the physical Hubble constant by H_{phys,0} = 100h_{0} km/sec/Mpc, and the critical energy density today by E_{cr} = 3H_{0}^{2}/(8πG). We use the subscript “CPI” to denote the time when the chiral plasma instability (CPI) develops. Assuming that the plasma’s entropy density is conserved between the CPI epoch and today leads to the relation g_{∗S,cpi}a_{cpi}^{3}T_{phys,cpi}^{3} = g_{∗S,0}a_{0}^{3}T_{phys,0}^{3}. Taking g_{∗S,0} = 3.91 and T_{phys,0} = 0.234 meV gives

\frac{a_{cpi}}{a_{0}} = (8 × 10^{−19})(\frac{g_{∗S,cpi}}{106.75})^{−1/3}(\frac{T_{phys,cpi}}{100 TeV} )^{−1} . (1)

We fiducialize the effective number of relativistic degrees of freedom during the CPI epoch to g_{∗S,cpi} = 106.75, which is the expected value for Standard Model cosmology at temperatures above 100 GeV. We fiducialize the physical plasma temperature at the CPI epoch to T_{phys,cpi} = 100 TeV, and the expected value varies according to Eq. (2).

a. Chiral magnetic effect. The chiral plasma instability and chiral magnetic effect (CME) were first studied in the context of a relativistic electron-positron plasma described by quantum electrodynamics (QED). Although chirality is conserved at the classical level for massless electrons, chirality is broken in the quantum theory and this is expressed by the Adler-Bell-Jackiw axial anomaly [15, 16]. A manifestation of the anomalous chiral symmetry is the CME [17]: in a QED plasma that possesses a chiral asymmetry, a magnetic field induces a proportional current. The CME corresponds to an anomalous contribution to the electric current density J(x, t) = μ_{5}(t)B(x, t) where μ_{5} = 2αμ̃_{5}/π is proportional to the chiral chemical potential μ̃_{5}, α = e^{2}/4π ≈ 1/137 is the electromagnetic fine structure constant, and B is the magnetic field. Implications of the CME for a turbulent QED plasma have been studied extensively with a combination of analytical techniques and numerical simulations [12, 14, 18–21]; see also Ref. [22] for a recent review article.

b. Adaptation to hypercharge. The formalism used to study the CME in QED is easily adapted to the hypercharge sector of the Standard Model for a plasma in the phase of unbroken electroweak symmetry at temperatures T_{phys} > ∼ 100 GeV. The quantity of interest is the hypercharge-weighted chiral chemical potential μ̃_{Y,5}(t), which is given by μ̃_{Y,5}(t) = ∑_{i} ε_{i}g_{i}Y_{i}^{2}μ̃_{i}(t), where the sum runs over all Standard Model particle species (indexed by i), ε_{i} = ±1 for right/left-chiral particles (and 0 otherwise), g_{i} is a multiplicity factor (counting color, spin, etc), Y_{i} is the hypercharge of species i, and μ̃_{i} is the chemical potential that parameterizes the asymmetry (excess of particles over antiparticle partners) in species i via n_{i} ∝ μ̃_{i}T^{2}; see Refs. [23, 24] for additional details. c. Chiral plasma instability. In the presence of a chiral asymmetry, the equations of magnetohydrodynamics (MHD) are modified due to the CME, and the new equations exhibit a tachyonic instability toward the growth of long-wavelength modes of the magnetic field, which is known as the chiral plasma instability [8]. To illustrate the instability in the hypercharge sector of the primordial plasma, we present the evolution equation for the hyper-magnetic field assuming negligible plasma velocity: Ḃ_{Y} = η_{Y} ∇^{2}B_{Y} + (2α_{Y} μ̃_{Y,5}/π)η_{Y} ∇ × B_{Y}. Here and below, dots represent partial derivatives with respect to conformal time, B_{Y} (x, t) is the hypermagnetic field, η_{Y} = 1/σ_{Y} is the hypermagnetic diffusivity, σ_{Y} is the hypercharge conductivity, and α_{Y} = g^{′2}/4π ≈ 0.01 is the hypercharge fine structure constant. Long-wavelength modes of the hypermagnetic field with wave number k < k_{cpi} = 2α_{Y} |μ̃_{Y,5}(t)|/π experience a tachyonic instability in one of the two circular polarization modes, and their amplitude increases exponentially ∝ exp(t/t_{cpi}). The fastest growing modes have k = k_{cpi}/2, and for these modes t_{cpi} = 4/η_{Y} k_{cpi}^{2} = π^{2}/η_{Y} α_{Y}^{2} |μ̃_{Y,5}|^{2}. Assuming a radiation-dominated cosmology with g_{∗} = 106.75, the physical plasma temperature at this time is

T_{phys,cpi} = (70 TeV)(\frac{η_{Y}}{0.01T^{−1}} )(\frac{|μ̃_{Y,5}|/T}{10^{−3}} )^{2} . (2)

In other words, although the chiral asymmetry may be present in the plasma from a very early time, its effect on the hypermagnetic field does not develop until (possibly much) later when the age of the universe is comparable to conformal time t_{cpi} and the plasma has cooled to temperature T_{cpi}. Reducing the magnitude of the chiral asymmetry, i.e., assuming a smaller |μ̃_{Y,5}|/T initially, delays the onset of the chiral plasma instability.

d. Chiral asymmetry erasure. In a relativistic electron-positron plasma described by the theory of QED, the electromagnetic charge is exactly conserved and the chiral charge is approximately conserved. The chiral charge changes in a scattering that converts rightchiral particles into left-chiral particles, or vice versa, and the rate for such ‘spin-flip’ scatterings is proportional the squared electron mass (m/T)^{2}. Although the chiral charge is not exactly conserved, it is important to recognize that it is approximately conserved on time scales that are small compared to the inverse slip-flip rate. Similarly, the hypercharge-weighted chiral asymmetry is eventually driven to zero by scatterings involving the Yukawa couplings; the most relevant processes are Higgs decays and inverse decays with right-chiral electrons. The rate for these chirality-changing reactions is Γ_{f} ≈ 10^{−2}y_{e}^{2}T, with y_{e} electron Yukawa coupling y_{e} = \sqrt{2}m_{e}/v ≃ 3 × 10^{−6} and assuming a standard radiationdominated cosmology, these reactions come into equilibrium when the plasma cools to a physical temperature of (T/a)_{f} ≃ 80 TeV [25]. To ensure that the chiral plasma instability develops before the hypercharge-weighted chiral asymmetry is erased by Higgs decays and inverse decays, it is necessary to have |μ̃_{Y,5}|/T > 10^{−3}. For reference, the observed baryon asymmetry of the universe today corresponds to a much smaller chemical potential of μ̃_{B}/T ≈ 10^{−8}, but it is not unusual for large chemical potentials to be generated during the course of baryogenesis. New physics such as a matter-dominated phase or an injection of e_{R} asymmetry can change the temperature of chiral asymmetry erasure; for an example, see Ref. [26].

e. Magnetogenesis. As the chiral plasma instability develops, the growing helical hypermagnetic field is accompanied by a depletion of the hypercharge-weighted chiral asymmetry. This is because the hyperchargeweighted chiral number density n_{Y,5} = μ̃_{Y,5}T^{2}/6 and the hypermagnetic helicity H_{B} are linked by the chiral anomaly, which imposes ṅ_{Y,5} ∝ −α_{Y} Ḣ_{M}/π [7]. If the chiral plasma instability shuts off after the hypercharge-weighted chiral asymmetry depletes by an order one factor, the hypermagnetic helicity can be estimated as H_{M,cpi} ∼ π|μ̃_{Y,5}|T_{cpi}^{2}/6α_{Y}. The coherence length and field strength are estimated as ξ_{M,cpi} ≈ 2π/(k_{cpi}/2) and B_{cpi} ≈ \sqrt{H_{M}/ξ_{M,cpi}}, which gives ξ_{M,cpi} ≈ (5 × 10^{5} cm)(|μ̃_{Y,5}|/10^{−3}T)^{−1} and B_{cpi} ≈ (5 × 10^{−11} G) (|μ̃_{Y,5}|/10^{−3}T). If the magnetic field evolves according to the inverse cascade scaling (in the fully helical case), ξ_{M} ∝ t^{2/3} and B ∝ t^{−1/3} [27], until recombination, then the physical coherence length and field strength today (assuming a frozen-in magnetic field and neglecting MHD dynamics at late epochs, after re-ionization) are expected to be on the order of

ξ_{M,phys,0} = (1 × 10^{−3} pc)(\frac{η_{Y}}{0.01T^{−1}} )^{2/3}(\frac{|μ̃_{Y,5}|/T}{10^{−3}} )^{1/3} ,

B_{phys,0} = (7 × 10^{−16} G)(\frac{η_{Y}}{0.01T^{−1}} )^{−1/3}(\frac{|μ̃_{Y,5}|/T}{10^{−3}} )^{1/3} .

(3)

A larger chiral asymmetry leads to a stronger magnetic field on larger length scales today.

f. Gravitational wave generation. The anisotropic stress of the growing hypermagnetic field provides a source of gravitational wave radiation [11]. As the chiral plasma instability develops, most of the magnetic energy is carried by the modes with coherence length ξ_{M,cpi} (i.e., the magnetic energy is characterized by the spectrum at the peak at wave number k_{I} ≃ 1/ξ_{M,cpi}). As long as the magnetic field is still growing, however, the induced gravitational wave spectrum peaks at the characteristic wave number k = 2/t_{cpi} = k_{cpi}^{2}/(2σ_{Y}) [12]. For k_{cpi}/(2σ_{Y}) < 1, this wave number is below the cutoff wave number for gravitational waves, k_{cpi}. Above this wave number, very little gravitational wave energy is produced by the chiral plasma instability [12]. The gravitational wave cutoff frequency is f_{GW} ≃ 2k_{I}/(2π) ≃ 2/ξ_{M,cpi} (the factor “2” is due to the quadratic nature of the source). Since the gravitational waves’ comoving frequency remains constant, the physical frequency today corresponds to f_{GW,0} = 2/ξ_{M,cpi}. Once the CPI stops and μ̃_{Y,5} becomes depleted, the low wave number part of the gravitational wave spectrum becomes shallower and the peak moves toward smaller wave numbers.

The energy density carried by the gravitational waves is estimated as E_{GW} ∼ (G/2π)a_{cpi}^{−2}ξ_{M,cpi}^{2}B_{cpi}^{4}. This estimate follows from writing the energy density as E_{GW}(x, t) = [ḣ_{ij}(x, t)]^{2}/(32πG) where h_{ij}(x, t) is the transverse and traceless tensor mode of the metric perturbations, and using the gravitational wave equation ∂_{t}^{2}h_{ij} − ∇^{2}h_{ij} = 16πGT_{ij}/a where T_{ij} ∼ B^{2}/2 is the transverse and traceless part of the anisotropic part of the magnetic field stress-energy tensor, to estimate the field amplitude [28, 29]. After volume averaging, we define Ω_{GW} = E_{GW}/E_{cr} to be the gravitational wave energy fraction today.

Numerical estimates give

\begin{aligned}f_{GW,0} & = (1 × 10^{5} Hz)(\frac{|μ̃_{Y,5}|/T}{10^{−3}} ), \\ Ω_{GW}h_{0}^{2} & = (7 × 10^{−39})(\frac{η_{Y}}{0.01T^{−1}} )^{2}(\frac{|μ̃_{Y,5}|/T}{10^{−3}} )^{6} . \end{aligned} (4)

A larger hypercharge-weighted chiral asymmetry moves the peak of the gravitational wave spectrum to higher frequencies (since the chiral plasma instability develops earlier) and increases the gravitational wave strength. For reference, the LIGO-Virgo-KAGRA gravitational wave interferometer array is sensitive to a stochastic gravitational wave background at the level of ∼ 10^{−7} for frequencies of ∼ 10–100 Hz [30]. The future space-based interferometer LISA will push this sensitivity down to ∼ 10^{−12} at frequencies of ∼ 1 − 10 mHz [31–33]. Various strategies for probing higher-frequency gravitational waves, even up to the GHz band, have been explored in recent years; see Ref. [34] for a review of these activities. Nevertheless, a detection of gravitational wave radiation at the level expected here, even for |μ̃_{Y,5}|/T ≈ 1, seems far out of reach.

g. Baryon number overproduction. The presence of a helical hypermagnetic field in the early universe is expected to give rise to a baryon asymmetry [23, 24, 35]. This is because time-varying hypermagnetic helicity sources baryon and lepton number through the electroweak anomaly [36]. Specifically, the conversion of a hypermagnetic field into an electromagnetic field at the electroweak epoch at T_{phys} ≈ 100 GeV sources baryon number after the electroweak sphaleron has gone out of equilibrium, leading to a boost in the baryon asymmetry [24].

The baryon number can easily be over-produced if the magnetic field strength is too large. Avoidance of this baryon-number overproduction imposes an upper bound of |μ̃_{Y,5}|/T < ∼ 10^{−2} [37]. This bound is somewhat uncertain as the baryon production calculation depends on a detailed modeling of magnetic field evolution at the Standard Model electroweak crossover [24], which is not well understood.

III. NUMERICAL SIMULATIONS

In order to validate the preceding estimates, we have performed three-dimensional numerical simulations using the Pencil Code [38]. These simulations allow us to study the growth and evolution of the magnetic field during the chiral plasma instability and to evaluate the spectrum of the resulting gravitational wave radiation.

We model the Standard Model matter and radiation as a single component plasma of charged particles interacting with the hypermagnetic field. Several properties of the plasma are relevant to the evolution: the magnetic diffusivity (for simplicity here and below we suppress the subscript “Y ”) η(t) = 1/σ(t), the kinematic viscosity ν(t), the chiral diffusion coefficient D_{5}(t), the chiral depletion parameter λ(t), and the chiral chemical potential μ_{50} ≡ μ_{5}(x, 0) = 2αμ̃_{5}/π that enters as an initial condition. One can calculate σ, η, ν, and D_{5} from first principles using Standard Model particle physics. The hypercharge conductivity is predicted to be σ ∼ T/α ≈ 100 T [39] implying η ≈ 0.01T^{−1}, and we assume for simplicity η = ν = D_{5}. The chiral depletion parameter λ arises from the Standard Model chiral anomalies, and past studies have obtained the prediction λ = 192 α^{2}/T^{2} ≃ 0.02 T^{−2} [21, 40]. The initial chiral chemical potential can be written as μ_{50} ≈ (6 × 10^{−6} T)(μ̃_{5}/T/10^{−3}) by fiducializing to μ̃_{5}/T = 10^{−3}.

Given the limited dynamic range of numerical simulations, it is not possible to set the parameters, η, ν, D_{5}, λ, and μ_{50}, equal to the Standard Model predictions. Instead we consider sets of simulations with different parameters. They can be distinguished by the relative ordering of the characteristic quantities v_{λ} = μ_{50}/(\overline{ρ}λ)^{1/2} and v_{μ} = μ_{50}η, where \overline{ρ} is the initial volume-averaged energy density. We consider runs in regimes I (where v_{λ} > v_{μ}) and II (where v_{λ} < v_{μ}).

The simulations solve a coupled system of partial differential equations that account for MHD and the CME [21] to determine the evolution of the magnetic field B(x, t), the energy density of the plasma ρ(x, t), the plasma velocity u(x, t), and the chiral chemical potential μ_{5}(x, t). In the following, we solve the following set of equations [40]

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

\frac{Du}{Dt} = \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,  (7)
\frac{∂ ln ρ}{∂t} = −\frac{4}{3} (∇ · u + u · ∇ ln ρ)
+\frac{1}{ρ} [u · (J × B) + ηJ^{2}] ,  (8)

where S_{ij} = (∂_{j}u_{i} + ∂_{i}u_{j})/2 − δ_{ij}∇ · u/3 are the components of the rate-of-strain tensor. We solve Eqs (5)– (7) using the Pencil Code, which is a massively parallel MHD code using sixth-order finite differences and a third-order time stepping scheme [41].

The linearized gravitational wave equations are solved in wave number space [29],

\frac{∂^{2}}{∂t^{2}} h̃_{+/×}(k, t) + k^{2}h̃_{+/×}(k, t) = \frac{6}{t}T̃_{+/×}(k, t), (9)

where h̃_{+/×} = e_{ij}^{+/×}(P_{il}P_{jm} − \frac{1}{2}P_{ij}P_{lm}) h̃_{lm}(k, t) are the Fourier-transformed + and × modes of h, with e_{ij}^{+}(k) = e_{i}^{1}e_{j}^{1} − e_{i}^{2}e_{j}^{2} and e_{ij}^{×}(k) = e_{i}^{1}e_{j}^{2} + e_{i}^{2}e_{j}^{1} being the linear polarization basis, e^{1} and e^{2} are unit vectors perpendicular to k and perpendicular to each other, and P_{ij}(k) = δ_{ij} − k_{i}k_{j} is the projection operator. T̃_{+/×} are defined analogously and normalized by the critical density. We solve Eq. (9) accurate to second order in the time step and use 1024^{3} mesh points in all of our calculations. Our initial conditions have a weak seed magnetic field and vanishing plasma velocities, and the chiral chemical potential is homogeneous and equal to the value given above. At each time step, we calculate the spectrum of gravitational wave radiation by solving the gravitational wave equation sourced by the stress-energy of the plasma and magnetic field; see Ref. [29] for details regarding our computational approach.

In Table I, we summarize our results, where η, λ, μ_{50}, and the smallest wave number, k_{1}. We consider two series of runs that we refer to as X and Y. They are subdivided further into Runs X1–X4 and Runs Y1–Y3. Our runs of series X have increasing values of v_{μ} and cross from regime Y (for Run I1) into regime I (for Run I4). For the runs of series X, we take η = ν = D_{5} = 5 × 10^{−11}/H_{∗}, \overline{ρ}λ = 10^{20} H_{∗}^{2}, and μ_{50} = 10^{6} H_{∗}. We also give the wave number corresponding to the outer scale of the k^{−2} inverse cascade range, k_{λ} ≈ (E_{∗}^{1/2}t_{∗}/l_{∗}^{5/2})4μ_{50}ηλ^{1/2}, as well as the efficiency of gravitational wave production,

q = (k_{peak}/H_{∗}) \sqrt{E_{GW}^{sat}E_{cr}}/ E_{M}^{max},  (10)

TABLE I: Summary of Runs discussed in this paper. Runs B1, B10, A1, and A12 of Ref. [12] are included for comparison. In the last row, theoretically expected values are listed where η_{2} = η_{Y} /(0.01T^{−1}) and μ_{3} = μ_{Y,5}/10^{−3}T.

Run ηH_{∗} (\overline{ρ}λ)^{1/2}/H_{∗} μ_{50}/H_{∗}  v_{μ}  v_{λ}  ημ_{50}^{2}/H_{∗} k_{1}/H_{∗} E_{M}^{max}/E_{cr}  E_{GW}^{sat}/E_{cr}  q
B1 1 × 10^{−6} 2 × 10^{4}  10^{4}  1 × 10^{−2}  5 × 10^{−1} 1 × 10^{2} 1 × 10^{2} 1.6 × 10^{−2} 4.7 × 10^{−12} 0.027
B10 1 × 10^{−3} 2 × 10^{4}  10^{4}  1 × 10^{1}  5 × 10^{−1} 1 × 10^{5} 1 × 10^{2} 6.0 × 10^{−2} 6.0 × 10^{−9}  12
A1 1 × 10^{−6} 5 × 10^{4}  10^{4}  1 × 10^{−2}  2 × 10^{−1} 1 × 10^{2} 1 × 10^{2} 4.6 × 10^{−3} 8.9 × 10^{−14} 0.032
A12 5 × 10^{−3} 5 × 10^{4}  10^{4}  5 × 10^{1}  2 × 10^{−1} 5 × 10^{5} 5 × 10^{1} 9.2 × 10^{−3} 3.0 × 10^{−10} 18
X1 5 × 10^{−8}  10^{10}  10^{6}  5 × 10^{−2}  1 × 10^{−4} 5 × 10^{4} 5 × 10^{3} 2.4 × 10^{−9} 8.8 × 10^{−31} 0.39
X2 5 × 10^{−9}  10^{10}  10^{6}  5 × 10^{−3}  1 × 10^{−4} 5 × 10^{3} 5 × 10^{3} 2.4 × 10^{−9} 1.6 × 10^{−30} 0.53
X3 5 × 10^{−10}  10^{10}  10^{6}  5 × 10^{−4}  1 × 10^{−4} 5 × 10^{2} 5 × 10^{3} 2.4 × 10^{−9} 1.1 × 10^{−30} 0.44
X4 5 × 10^{−11}  10^{10}  10^{6}  5 × 10^{−5}  1 × 10^{−4} 5 × 10^{1} 5 × 10^{3} 2.3 × 10^{−9} 3.1 × 10^{−31} 0.12
Y1 5 × 10^{−8} 7 × 10^{11}  10^{6}  5 × 10^{−2}  1 × 10^{−6} 5 × 10^{4} 5 × 10^{3} 4.9 × 10^{−13} 3.6 × 10^{−38} 0.39
Y2 5 × 10^{−8} 7 × 10^{11}  10^{6}  5 × 10^{−2}  1 × 10^{−6} 5 × 10^{4} 2 × 10^{3} 4.4 × 10^{−13} 3.2 × 10^{−37} 1.3
Y3 5 × 10^{−8} 7 × 10^{11}  10^{6}  5 × 10^{−2}  1 × 10^{−6} 5 × 10^{4} 1 × 10^{3} 3.3 × 10^{−13} 6.9 × 10^{−37} 2.5
exp 10^{−15}η_{2} 6 × 10^{12} 5 × 10^{7}μ_{3} 6 × 10^{−8}η_{2}μ_{3} 8 × 10^{−6}μ_{3} 3η_{2}μ_{3}^{2}  — 6 × 10^{−15}μ_{3}^{2} 7 × 10^{−39}η_{2}^{2}μ_{3}^{6} —

where we estimate k_{peak} = k_{μ} min(1, v_{μ}/v_{λ}). This means that k_{peak} = k_{μ} when v_{μ} > v_{λ} (regime II) and k_{peak} = k_{λ}/4 when v_{μ} < v_{λ} (regime I); see Ref. [21].

The evolution of the magnetic and gravitational wave energy spectra for Run I4 is shown in Fig. 1. The energy densities may be written as E = ∫_{0}^{∞}dk E(k) where k is the wavenumber and E(k) is the energy spectrum. For the Run I4 parameters, the instability length scale corresponds to a wavenumber of k_{cpi} = 10^{6}H_{∗}, which agrees with the wave number above which E_{GW}(k) drops sharply. The instability time scale normalized to the Hubble time is t_{cpi}/t_{H} = 0.08 (with t_{H} = H_{∗}^{−1} being the Hubble time) which is about 100 times longer than the time step. The magnetic energy spectrum grows initially for modes with k ≈ k_{cpi}/2 = 5 × 10^{5}H_{∗} (see the upper set of dotted lines in Fig. 1). Later, the peak evolves to smaller k with an inverse cascade scaling, which is consistent with earlier simulations [21]. The generated magnetic field is then maximally helical; see Fig. 8(b) of Ref. [42].

The gravitational wave energy spectra grow in time as long as the magnetic energy has not yet reached its maximum. In this phase, as discussed above, the gravitational wave spectrum is expected to peak at the characteristic wave number k = 2/t_{cpi} = ηk_{cpi}^{2}/2 = 25 H_{∗}, which is here much smaller than k_{cpi} = 10^{6}H_{∗}, but larger than the horizon wave number, k = H_{∗}. When the magnetic energy density has reached its maximum value, the gravitational wave spectrum has nearly saturated and is approximately independent of k for k < k_{cpi}/2. In principle, it is possible to have a declining k^{−2} spectrum in the range ηk_{cpi}^{2}/2 ≤ k ≤ k_{cpi}, but this is only seen in our models with larger diffusivity. The absence of a k^{−2} subrange in the gravitational wave spectrum could also be an artifact of insufficient numerical resolution. In any case, once the gravitational wave spectrum saturates, we would expect the development of a flat (E_{GW} ∝ k^{0}) spec-

FIG. 1: Spectra (per linear wave number interval) of magnetic energy E_{M}(k, t) (upper curves) and gravitational wave energy E_{GW}(k, t) (lower curves) from the chiral plasma instability and turbulent MHD evolution for Run I, where μ_{5} = 10^{6} H_{∗}, \overline{ρ}λ = 10^{20}/H_{∗}, η = 5×10^{−11} H_{∗} which implies v_{λ} = 10^{−4} and v_{μ} = 5×10^{−5} (corresponding to regime I). The solid curves are for tH_{∗} = 2.98, when E_{M} is maximum. The dotted curves are for tH_{∗} = 2.41 (black), 2.56 (red), and 2.71 (orange), before E_{M} is maximum, while the dashed curves are for tH_{∗} = 3.66 (blue) and tH_{∗} = 5.37 (green), when E_{M} is decaying.

trum. Such a flat spectrum is expected to extend all the way to the horizon wave number k = H_{∗} [43–45]. Since the smallest wave number that is resolved in our simulations is already larger than the horizon wave number, we expect that the gravitational wave energy is underestimated by our simulations by an amount that depends on the value of the smallest resolved wave number, k_{1}. In Fig. 2, we show gravitational wave spectra for a few runs with smaller values of the minimum wave number in the simulations. We see that the spectra remain nearly flat, but the spectra are also becoming more irregular at

FIG. 2: Comparison of h_{0}^{2}Ω_{GW}(f) versus f for runs with k_{1}/H_{∗} = 10^{3} (blue), 2 × 10^{3} (orange), and 5 × 10^{3} (red), with λ = 49 × 10^{22}H_{∗}^{2}/E_{cr,∗} and η = 5 × 10^{−8}H_{∗}^{−1} with H_{0} = 100 h_{0} km s^{−1} Mpc^{−1}. A run with with λ = 10^{20}H_{∗}^{2}/E_{cr,∗} and again η = 5×10^{−8}H_{∗}^{−1} is shown as the black line for reference.

large wave numbers. This is likely an artifact of insufficient numerical resolution. We also see that most of the gravitational wave energy is at frequencies below about 1 kHz, but this value would increase with increasing values of μ_{50}, beyond the value of 10^{6}H_{∗} adopted here.

Earlier work showed that E_{GW}^{sat} grows approximately linearly with η and was proportional to (\overline{ρ}λ)^{−5/2}, which leads to the combined dependence [12]

E_{GW}^{sat}/E_{cr} ≈ 6 × 10^{−8} v_{λ}^{5}v_{μ},  (11)

which implies that E_{GW}^{sat} ∝ μ_{50}^{6}. In Fig. 3, we plot E_{GW}^{sat} versus v_{λ}^{5}v_{μ} for Runs X1–X4 and Y1–Y3. We see that Eq. (11) agrees reasonably well with our numerical data. Compared with the runs of Ref. [12], the new one in Fig. 1 has much smaller values of v_{λ} (here v_{λ} = 10^{−4} instead of 0.5 for the old runs of Series B) and v_{μ} (here v_{μ} = 5 × 10^{−5} instead of 10^{−2}, which was their smallest value). This has been achieved by having k_{cpi} much larger (here 10^{6} instead of 10^{4}, for example). This also means that we have to choose a correspondingly larger value of the minimum wave number, k_{1}. This means that we consistently underestimate the value of E_{GW}, because the contribution from low values of k are not resolved.

IV. CONCLUSIONS

Our estimates of the key variables are summarized in Fig. 4. Since we assume Standard Model particles and interactions, as well as a standard cosmology with radiation domination at temperatures T > 100 TeV, the observables depend only on the single dimensionless parameter |μ_{Y,5}|/T, which controls the size of the initial hypercharge-weighted chiral asymmetry. To ensure that the instability develops before the chiral asymmetry is

FIG. 3: Dependence of E_{GW}^{sat} on v_{λ}^{5}v_{μ} for our runs of Series X (red) and Y (blue), as well as Series A (black) and B (orange) of Ref. [12].

washed out by reactions such as Higgs decays and inverse decays, we need |μ_{Y,5}|/T > ∼ 10^{−3}. On the other hand, to avoid over-producing the baryon asymmetry we need |μ_{Y,5}|/T < ∼ 10^{−2}. This leaves an approximately onedecade wide window of viable parameter space. The predicted magnetic field strength today, assuming inverse cascade scaling from production until recombination, is at the level of 10^{−15} Gauss. An intergalactic magnetic field at this level is barely strong enough to explain observations of distant TeV blazars, which provide evidence for a nonzero intergalactic magnetic field at the level > ∼ 10^{−16} Gauss [9]. The same magnetic field may help to explain the origin of galactic magnetic fields by providing a seed for the galactic dynamo. The strength of the gravitational wave signal is expected to depend strongly on the value of |μ_{Y,5}|/T, going as its sixth power. The typical frequency of these this signal is expected to fall near ∼ GHz, putting it a frequency band that is being targeted by several recently-proposed probes of highfrequency gravitational wave radiation. However, within the viable window, the gravitational wave signal is likely far too weak for detection.

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

Acknowledgements—We thank Kohei Kamada, Sayan Mandal, and Jonathan Stepp for useful discussions and comments. A.J.L. is grateful to the Mainz Institute for Theoretical Physics (MITP) of the DFG Cluster of Excellence PRISMA+ (Project ID 39083149), for its hospitality during the completion of this work. Support through the NASA ATP award 80NSSC22K0825, the Swedish Research Council, grant 2019-04234, and Shota Rustaveli GNSF (grant FR/18-1462) are gratefully acknowledged. We acknowledge the allocation of computing resources provided by the Swedish National Alloca-

FIG. 4: Summary figure. AL: I have to change μ_{Y,5} to μ̃_{Y,5}

tions Committee at the Center for Parallel Computers at the Royal Institute of Technology in Stockholm. E.C. acknowledges the Pake Fellowship, and G.S. acknowledges support from the Undergraduate Research Office in the form of a Summer Undergraduate Research Fellowships at Carnegie Mellon University.

[1] E. W. Kolb and S. Wolfram, Nucl. Phys. B 172, 224 (1980), [Erratum: Nucl.Phys.B 195, 542 (1982)].

[2] M. Fukugita and T. Yanagida, Phys. Lett. B174, 45 (1986).

[3] G. Servant and S. Tulin, (2013), arXiv:1304.3464.

[4] K. Dick, M. Lindner, M. Ratz, and D. Wright, Phys.Rev.Lett. 84, 4039 (2000), arXiv:hep-ph/9907562.

[5] H. Murayama and A. Pierce, Phys. Rev. Lett. 89, 271601 (2002), arXiv:hep-ph/0206177.

[6] B. A. Campbell, S. Davidson, J. R. Ellis, and K. A. Olive, Phys.Lett. B256, 484 (1991).

[7] M. Joyce and M. E. Shaposhnikov, Phys.Rev.Lett. 79, 1193 (1997), arXiv:astro-ph/9703005.

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

[9] A. Neronov and I. Vovk, Science 328, 73 (2010), arXiv:1006.3504.

[10] T. Vachaspati, Rept. Prog. Phys. 84, 074901 (2021), arXiv:2010.10525.

[11] D. V. Deryagin, D. Y. Grigoriev, V. A. Rubakov, and M. V. Sazhin, Mod. Phys. Lett. A 1, 593 (1986).

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

[13] A. Brandenburg, K. Kamada, and J. Schober, Phys. Rev. Res. 5, L022028 (2023), arXiv:2302.00512.

[14] A. Brandenburg, K. Kamada, K. Mukaida, K. Schmitz, and J. Schober, (2023), arXiv:2304.06612.

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

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

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

[18] A. Boyarsky, J. Fröhlich, and O. Ruchayskiy, Phys. Rev. Lett. 108, 031301 (2012), arXiv:1109.3350.

[19] A. Boyarsky, O. Ruchayskiy, and M. Shaposhnikov, Phys. Rev. Lett. 109, 111602 (2012), arXiv:1204.3604.

[20] A. Boyarsky, J. Fröhlich, and O. Ruchayskiy, Phys. Rev. D 92, 043004 (2015), arXiv:1504.04854.

[21] A. Brandenburg et al., Astrophys. J. Lett. 845, L21 (2017), arXiv:1707.03385.

[22] K. Kamada, N. Yamamoto, and D.-L. Yang, Prog. Part.

Nucl. Phys. 129, 104016 (2023), arXiv:2207.09184.

[23] K. Kamada and A. J. Long, Phys. Rev. D 94, 063501 (2016), arXiv:1606.08891.

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

[25] D. Bödeker and D. Schröder, JCAP 05, 010 (2019), arXiv:1902.07220.

[26] M.-C. Chen, S. Ipek, and M. Ratz, Phys. Rev. D 100, 035011 (2019), arXiv:1903.06211.

[27] T. Hatori, JPSJ 53, 2539 (1984).

[28] G. Gogoberidze, T. Kahniashvili, and A. Kosowsky, Phys. Rev. D 76, 083002 (2007), arXiv:0705.1733.

[29] A. Roper Pol, A. Brandenburg, T. Kahniashvili, A. Kosowsky, and S. Mandal, Geophys. Astrophys. Fluid Dynamics 114, 130 (2020), arXiv:1807.05479.

[30] LIGO Scientific, Virgo, B. P. Abbott et al., Phys. Rev. Lett. 118, 121101 (2017), arXiv:1612.02029, [Erratum: Phys.Rev.Lett. 119, 029901 (2017)].

[31] N. Bartolo et al., JCAP 11, 034 (2018), arXiv:1806.02819.

[32] J. Baker et al., (2019), arXiv:1907.06482.

[33] C. Caprini et al., JCAP 11, 017 (2019), arXiv:1906.09244.

[34] N. Aggarwal et al., Living Rev. Rel. 24, 4 (2021), arXiv:2011.12414.

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

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

[37] V. Domcke, K. Kamada, K. Mukaida, K. Schmitz, and M. Yamada, (2022), arXiv:2208.03237.

[38] The pencil code. doi:10.5281/zenodo.2315093. https://github.com/pencil-code.

[39] P. B. Arnold, G. D. Moore, and L. G. Yaffe, JHEP 11, 001 (2000), arXiv:hep-ph/0010177.

[40] I. Rogachevskii et al., Astrophys. J. 846, 153 (2017), arXiv:1705.00378.

[41] Pencil Code, A. Brandenburg et al., J. Open Source Softw. 6, 2807 (2021), arXiv:2009.08231.

[42] A. Brandenburg, Y. He, and R. Sharma, Astrophys. J.

922, 192 (2021), arXiv:2107.12333.

[43] A. Roper Pol, S. Mandal, A. Brandenburg, T. Kahniashvili, and A. Kosowsky, Phys. Rev. D 102, 083512 (2020), arXiv:1903.08585.

[44] A. Roper Pol, C. Caprini, A. Neronov, and D. Semikoz, Phys. Rev. D 105, 123502 (2022), arXiv:2201.05630.

[45] R. Sharma and A. Brandenburg, Phys. Rev. D 106, 103536 (2022), arXiv:2206.00055.

                                                    [46] http://norlx65.nordita.org/~brandenb/projects/GWfromSM/

contains data sets and run directories for the runs of Table I.