# Leading-order nonlinear gravitational waves from reheating magnetogeneses

Yutong He<EMAIL_ADDRESS>Nordita, KTH Royal Institute of Technology and
Stockholm University,Hannes Alfvéns väg 12, 10691 Stockholm, Sweden
Department of Astronomy, AlbaNova University Center, Stockholm University,
10691 Stockholm, Sweden  Alberto Roper Pol<EMAIL_ADDRESS>Université
de Paris, CNRS, Astrophysique et Cosmologie, Paris, F-75013, France  School of
Natural Sciences and Medicine, Ilia State University, 3-5 Cholokashvili Ave,
Tbilisi, GE-0194, Georgia  Axel Brandenburg<EMAIL_ADDRESS>Nordita, KTH
Royal Institute of Technology and Stockholm University,Hannes Alfvéns väg 12,
10691 Stockholm, Sweden  Department of Astronomy, AlbaNova University Center,
Stockholm University, 10691 Stockholm, Sweden  School of Natural Sciences and
Medicine, Ilia State University, 3-5 Cholokashvili Ave, Tbilisi, GE-0194,
Georgia  McWilliams Center for Cosmology and Department of Physics, Carnegie
Mellon University, 5000 Forbes Ave, Pittsburgh, PA 15213, USA

###### Abstract

We study the leading-order nonlinear gravitational waves (GWs) produced by an
electromagnetic (EM) stress in reheating magnetogenesis scenarios. Both
nonhelical and helical magnetic fields are considered. By numerically solving
the linear and leading-order nonlinear GW equations, we find that the GW
energy from the latter is usually larger. We compare their differences in
terms of the GW spectrum and parameterize the GW energy difference due to the
nonlinear term, $\Delta\mathcal{E}_{\rm GW}$, in terms of EM energy
$\mathcal{E}_{\rm EM}$ as $\Delta\mathcal{E}_{\rm
GW}=(\tilde{p}\mathcal{E}_{\rm EM}/k_{*})^{3}$, where $k_{*}$ is the
characteristic wave number, $\tilde{p}=0.84$ and $0.88$ are found in the
nonhelical and helical cases, respectively, with reheating around the QCD
energy scale, while $\tilde{p}=0.45$ is found at the electroweak energy scale.
We also compare the polarization spectrum of the linear and nonlinear cases
and find that adding the nonlinear term usually yields a decrease in the
polarization that is proportional to the EM energy density. We parameterize
the fractional polarization suppression as $|\Delta\mathcal{P}_{\rm
GW}/\mathcal{P}_{\rm GW}|=\tilde{r}\mathcal{E}_{\rm EM}/k_{*}$ and find
$\tilde{r}=1.2\times 10^{-1}$, $7.2\times 10^{-4}$, and $3.2\times 10^{-2}$
for the helical cases with reheating temperatures $T_{\rm r}=300{\rm\,TeV}$,
$8{\rm\,GeV}$, and $120{\rm\,MeV}$, respectively. Prospects of observation by
pulsar timing arrays, space-based interferometers, and other novel detection
proposals are also discussed.

gravitational waves, nonlinear memory, reheating, magnetogenesis

††preprint: NORDITA-2021-080

##  I Introduction

The recent detections of gravitational wave (GW) events produced by the
collision of compact binary objects by the LIGO-Virgo collaboration mark the
beginning of the rich field of GW astronomy Abbott  _et al._ (2016a, b, 2017).
Over 50 events have been detected so far by Advanced LIGO and Advanced Virgo,
which are collected in the GW transient catalogs Abbott  _et al._ (2019,
2021). These events yield propagating oscillatory strains of the metric
tensor, which were predicted by linearized general relativity (GR). On the
other hand, a different type of perturbations in the spacetime metric yielding
non-oscillatory strains have been predicted by nonlinear GR, although they
have not yet been detected. This effect, known as GW memory, is manifested by
a permanent displacement of freely-falling test masses after the passage of
GWs Zel’dovich and Polnarev (1974); Turner (1978); Epstein (1978); Braginskii
and Thorne (1987); Christodoulou (1991); Blanchet and Damour (1992); Thorne
(1992). It is a consequence of the facts that GR admits highly-degenerate
Bondi-Metzner-Sachs (BMS) vacua Bondi  _et al._ (1962); Sachs (1962) and that
the initial and final metrics, although both flat, differ by a BMS
supertranslation induced by the passing GWs. Theoretically, it has been shown
that the GW memory effect is equivalent to the temporal Fourier transform of
the Weinberg soft graviton theorem Strominger and Zhiboedov (2016); He  _et
al._ (2015), which is a universal formula relating to each other two
scattering matrices that differ by a zero-energy (soft) graviton Weinberg
(1965).

Recently, the memory effect has seen a pick-up in interest in the GW community
extensively due to the expected improvement in the coming years in sensitivity
by upgraded GW detectors. In particular, the advanced LIGO-Virgo-KAGRA network
will be able to detect an ensemble of memory signals from binary mergers Lasky
_et al._ (2016); McNeill  _et al._ (2017); Hübner  _et al._ (2020); Boersma
_et al._ (2020). On the other hand, space-based interferometers like the Laser
Interferometer Space Antenna (LISA), planned to be launched in 2034 Amaro-
Seoane  _et al._ (2017), will be capable of storing the permanent
displacements since the test masses are actually free-falling, unlike in
ground-based detectors Favata (2009a). It is expected that LISA will detect
compact binaries memory signals individually Favata (2009a); Islo  _et al._
(2019); Yang and Martynov (2018). In the lower-frequency end, pulsar timing
arrays (PTAs) are also considered as potential memory detectors and certain
constraints have been put on the detectable signals Wang  _et al._ (2015);
Arzoumanian  _et al._ (2015); Madison  _et al._ (2017); Aggarwal  _et al._
(2020). Astrophysical sources of the memory effect have been studied in binary
black holes Favata (2009a); Pollney and Reisswig (2011); Khera  _et al._
(2021); Ebersold and Tiwari (2020); Islo  _et al._ (2019), supernova neutrinos
Turner (1978); Epstein (1978); Mukhopadhyay  _et al._ (2021), and gamma-ray
burst jets Burrows and Hayes (1996); Segalis and Ori (2001); Sago  _et al._
(2004); Hiramatsu  _et al._ (2005). As for cosmological sources, cosmic
strings have been considered Aurrekoetxea  _et al._ (2020); Jenkins and
Sakellariadou (2021). In particular, the memory strain from compact binary
objects has been expected to be only one order of magnitude weaker than the
standard GW strain Favata (2009a); Khera  _et al._ (2021). As a test of
gravity theories, the memory effect has been studied as a result of massive
gravity Kilicarslan and Tekin (2019), scalar-tensor theories Hou (2021); Hou
and Zhu (2021); Koyama (2020), as well as in the strong-field regime Shore
(2018). In addition, novel phenomena, different from the conventional
displacement memory effect, have been proposed. These are the velocity
Grishchuk and Polnarev (1989); Zhang  _et al._ (2017, 2018); Divakarla and
Whiting (2021) and spin memory effect Pasterski  _et al._ (2016); Nichols
(2017); Mitman  _et al._ (2020); Tahura  _et al._ (2021), where the freely-
falling particles would acquire permanent relative velocities and spins,
respectively, after the passage of GWs.

However, despite the plethora of literature on GW memory effect to date, to
the best of our knowledge, none has explored the effect with continuous
sources on cosmological scales, which are ubiquitous in the early universe.
With the recent advance of numerical simulations, there has been significant
progress in calculating the GW backgrounds (GWBs) from distributed sources
such as homogeneous hydrodynamic and hydromagnetic turbulence in the early
universe Roper Pol  _et al._ (2020a, b). We are interested in studying the
leading-order nonlinear term producing gravitational radiation, for which we
require strong sources. For this reason, we consider in this work the
production of GWs by reheating magnetogeneses, which yield electromagnetic
(EM) strengths that depend only on the initial field and can grow several
orders of magnitude Brandenburg and Sharma (2021); Brandenburg  _et al._
(2021a). For simplicity, we defer the study of the nonlinear effect produced
by hydromagnetic turbulence from cosmological phase transitions to future
work, since it requires a fully relativistic framework to reach plasma and/or
Alfvén velocities near speed of light, as expected for very strong sources.
There have been recent analytical works on nonhelical and helical
magnetogeneses during the reheating era Sharma  _et al._ (2017, 2018), which
circumvent known difficulties such as the backreaction and strong coupling
problems Demozzi  _et al._ (2009). Numerical simulations of GWs from these
magnetogeneses have also been performed Brandenburg and Sharma (2021);
Brandenburg  _et al._ (2021a) using the Pencil Code Collaboration (2021).

This paper is organized as follows. In Sec. II, we introduce the
magnetogenesis models that will be used to generate the EM sources of GWs. In
Sec. III, we discuss the relevant GW equation at the leading order beyond the
standard linear case. In Sec. IV, we present the simulation parameters,
discuss the results, and study the potential detectability of the sources
considered. Finally, we conclude in Sec. V.

We use the $(-\\!\\!+\\!\\!++)$ metric signature, set $c=1$, and normalize the
critical energy density by the end of reheating $\eta_{*}$ to be unity, i.e.,
$\rho_{\rm crit}(\eta_{*})=1$, with $\eta$ denoting conformal time. We also
adopt the convention that, for a function $F$, the Fourier transformation and
its inverse read $\tilde{F}(\eta,\bm{k})=\int
F(\eta,\bm{x})e^{-i\bm{k}\cdot\bm{x}}{\rm\,d}^{3}\bm{x}$ and
$F(\eta,\bm{x})=\int\tilde{F}(\eta,\bm{k})e^{i\bm{k}\cdot\bm{x}}{\rm\,d}^{3}\bm{x}/(2\pi)^{3}$,
respectively. For long expressions in Fourier space, we also adopt the
notation ${\cal F}(F(\eta,\bm{x}))=\tilde{F}(\eta,\bm{k})$. In addition, Greek
letters are employed for space-time indices running from 0 to 3 and Latin
letters are used for space indices 1 to 3. Finally, variables by default refer
to linear solutions but the ones with ‘nlin’ superscript denote nonlinear
solutions.

##  II Reheating magnetogenesis

A promising mechanism that can produce helical and nonhelical relic magnetic
fields today is reheating magnetogenesis. It is made possible by the coupling
of the inflaton field and the Maxwell Lagrangian $\mathcal{L}^{\rm
EM}=F_{\mu\nu}F^{\mu\nu}$ Ratra (1992), where
$F_{\mu\nu}\equiv\partial_{\mu}A_{\nu}-\partial_{\nu}A_{\mu}$ is the Faraday
tensor in terms of the EM four-potential $A_{\mu}$. However, due to the
conformal invariance of $\mathcal{L}^{\rm EM}$, the EM amplitudes decay
rapidly as the square of the scale factor during inflation. Therefore,
violation of conformal invariance is required to obtain a growing EM field
Turner and Widrow (1988); Dolgov (1993). This can be done, in the nonhelical
case, via a coupling term of the form $\mathcal{L}^{\rm
nhel}=f^{2}F_{\mu\nu}F^{\mu\nu}$, with $f$ being a time-dependent function
that scales as

$f(a)\begin{cases}\propto a^{\alpha}\;\;&(\text{during inflation}),\\\ \propto
a^{-\beta}\;\;&(\text{during reheating}),\\\ =1\;\;&(\text{during
RD}),\end{cases}$ (1)

where RD refers to the radiation-dominated era, which starts at the end of
reheating, $a$ is the scale factor, $\alpha$ is chosen to be 2, to avoid the
backreaction problem Sharma  _et al._ (2017), or 1, which enables reheating
temperatures above the electroweak (EW) scale Brandenburg  _et al._ (2021a)
(see Sec. IV), and $\beta>0$ parameterizes the reheating temperature $T_{\rm
r}$. The detailed relations between $T_{\rm r}$ and $\beta$ can be found in
Ref. Sharma  _et al._ (2017, 2018). Note that, in order to realize the
standard electromagnetism, $f(a)=1$ is needed by the onset of radiation era,
which we normalize to $\eta_{*}=1$ Roper Pol  _et al._ (2020a).

On the other hand, helical magnetic fields can also be generated by a modified
coupling term of the form $\mathcal{L}^{\rm hel}=\gamma
f^{2}F_{\mu\nu}\tilde{F}^{\mu\nu}$, where $\gamma\neq 0$ is a constant, $f$
obeys the same scaling relations as in Eq. (1), and
$\tilde{F}^{\mu\nu}=\epsilon^{\mu\nu\alpha\beta}F_{\alpha\beta}/2$ is the dual
of the Faraday tensor, with $\epsilon^{\mu\nu\alpha\beta}$ being the fully
antisymmetric tensor.

In addition, recall that the conformal time $\eta$ relates to the scale factor
$a$ as $a\propto\eta^{2/(1+3w)}$, where $w\equiv p/\rho$ is the equation of
state of the universe. During the radiation era with $w=1/3$, we have
$a\propto\eta$. During the matter-dominated reheating era, $w=0$ gives
$a\propto\eta^{2}$. Since we adopt the normalization such that $a_{*}=1$ marks
the onset of the radiation era at $\eta_{*}=1$ Roper Pol  _et al._ (2020a), we
have the following relations Sharma  _et al._ (2017)

$a=\begin{cases}(\eta+1)^{2}/4\;\;&(\text{during reheating}),\\\
\eta\;\;&(\text{during RD}).\end{cases}$ (2)

###  II.1 Equations of motion for the four-potential

The EM vector potential $\tilde{\bm{A}}$ can split into $\tilde{A}_{+}$ and
$\tilde{A}_{-}$ polarization modes in Fourier space as Varshalovich  _et al._
(1988)

$\tilde{\bm{A}}(\eta,\bm{k})=\tilde{A}_{+}(\eta,\bm{k})\tilde{\bm{e}}_{+}(\bm{k})+\tilde{A}_{-}(\eta,\bm{k})\tilde{\bm{e}}_{-}(\bm{k}),$
(3)

where

$\tilde{\bm{e}}_{\pm}(\bm{k})=-\frac{i}{\sqrt{2}}[\tilde{\bm{e}}_{1}(\bm{k})\pm\tilde{\bm{e}}_{2}(\bm{k})]$
(4)

is the polarization basis satisfying $i\bm{k}\times\tilde{\bm{e}}_{\pm}=\pm
k\tilde{\bm{e}}_{\pm}$. Here $k=|\bm{k}|$ is the wave number, and
$\tilde{\bm{e}}_{1}(\bm{k})$ and $\tilde{\bm{e}}_{2}(\bm{k})$ are unit vectors
orthogonal to $\bm{k}$ and to each other.

In Fourier space, the relevant equations for the four-potential modes
$\tilde{A}_{\pm}$ during reheating take the form Sharma  _et al._ (2018);
Brandenburg  _et al._ (2021a)

$\tilde{\mathcal{A}}_{\pm}^{\prime\prime}+\left(\bm{k}^{2}\pm 2\gamma
k\frac{f^{\prime}}{f}-\frac{f^{\prime\prime}}{f}\right)\tilde{\mathcal{A}}_{\pm}=0,$
(5)

where $\tilde{\mathcal{A}}_{\pm}\equiv f\tilde{A}_{\pm}$ is the scaled four-
potential and
$\tilde{\mathcal{A}}_{\pm}^{\prime\prime}\equiv\partial^{2}\tilde{\mathcal{A}}_{\pm}/\partial\eta^{2}$
is the second derivative with respect to conformal time $\eta$. Here $\gamma$
is a constant, and the terms involving the coupling function $f$, scale factor
$a$, and their derivatives can be expressed using Eqs. (1) and (2), yielding

$\frac{a^{\prime\prime}}{a}=\frac{2}{(\eta+1)^{2}}\;,\;\frac{f^{\prime}}{f}=-\frac{2\beta}{\eta+1}\;,\;\frac{f^{\prime\prime}}{f}=\frac{2\beta(2\beta+1)}{(\eta+1)^{2}}.$
(6)

In the helical basis, helicity is proportional to the difference between the
$+$ and $-$ modes. This means that
$\tilde{\mathcal{A}}_{+}=\tilde{\mathcal{A}}_{-}$ in the nonhelical case,
where $\gamma=0$ and the modes are governed by $f^{\prime\prime}/f$ only. In
the helical case, $\gamma=1$ and both $f^{\prime}/f$ and $f^{\prime\prime}/f$
play a role in the solution. The helical term in Eq. (5), proportional to
$\gamma$, leads to a difference in the growth rates between
$\tilde{\mathcal{A}}_{+}$ and $\tilde{\mathcal{A}}_{-}$. In our simulations,
since there is enough time for the field to grow over many orders of
magnitude, we consider only the dominant polarization mode in practice.

##  III GW Equations

###  III.1 Leading-order nonlinear equation

Given the metric $g_{\mu\nu}=\bar{g}_{\mu\nu}+h_{\mu\nu}$, where
$\bar{g}_{\mu\nu}$ is the background metric tensor and $|h_{\mu\nu}|\ll 1$
correspond to high-frequency (compared to the slowly varying background) small
perturbations, then the Einstein field equations (EFEs) can be expanded to
arbitrary orders in $h_{\mu\nu}$, following Refs. Isaacson (1968a, b); see Ch.
20.3 of Ref. Misner  _et al._ (1973) and Ch. 1 of Ref. Maggiore (2007) for
textbook references.

In general, EFEs read $G_{\mu\nu}=\kappa T_{\mu\nu}$. Here $T_{\mu\nu}$ is the
stress-energy tensor of the source, $G_{\mu\nu}\equiv
R_{\mu\nu}-\frac{1}{2}g_{\mu\nu}R$ is the Einstein tensor in terms of the
Ricci tensor $R_{\mu\nu}$ and Ricci scalar $R$, and $\kappa\equiv 8\pi G_{\rm
N}$, where $G_{\rm N}$ is Newton’s gravitational constant. As we are
interested in nonlinear GWs, we expand $G_{\mu\nu}$ up to, at least, quadratic
order in $h_{\mu\nu}$, giving

$G_{\mu\nu}=\bar{G}_{\mu\nu}+G_{\mu\nu}^{(1)}+G_{\mu\nu}^{(2)}+{\cal
O}(h^{3})=\kappa T_{\mu\nu},$ (7)

where the superscripts (1) and (2) denote $\mathcal{O}(h)$ and
$\mathcal{O}(h^{2})$, respectively, and $\bar{G}_{\mu\nu}$ depends on
$\bar{g}_{\mu\nu}$. The equation for the $\mathcal{O}(h)$ terms can be recast
as

$\displaystyle G_{\mu\nu}^{(1)}=\kappa(T_{\mu\nu}+t_{\mu\nu}),$ (8)

where we have introduced the effective stress of GWs, generally referred to as
the Landau-Lifshitz (LL) pseudotensor Landau and Lifshitz (1962),
$t_{\mu\nu}=-G_{\mu\nu}^{(2)}/\kappa$, with $G_{\mu\nu}^{(2)}$ satisfying

$G_{\mu\nu}^{(2)}=R_{\mu\nu}^{(2)}-\frac{1}{2}\bar{g}_{\mu\nu}R^{(2)}.$ (9)

The RHS of Eq. (8) shows the two terms sourcing GWs. Rewriting with the trace-
reversed perturbations $\bar{h}_{\mu\nu}\equiv
h_{\mu\nu}-\frac{1}{2}\bar{g}_{\mu\nu}h$, with $h=\bar{g}^{\mu\nu}h_{\mu\nu}$,
and applying the harmonic gauge $\partial^{\mu}\bar{h}_{\mu\nu}=0$, the
linearized Einstein tensor becomes
$G_{\mu\nu}^{(1)}=-\frac{1}{2}\Box\bar{h}_{\mu\nu}$, where
$\Box\equiv\partial_{\rho}\partial^{\rho}=-\partial_{t}^{2}+\bm{\nabla}^{2}$
is the d’Alembert operator. Therefore, the linearized EFEs, i.e., omitting
terms of order ${\cal O}(h^{2})$ and higher, correspond to the usual GW
equation:

$\Box\bar{h}_{\mu\nu}=-2\kappa T_{\mu\nu},$ (10)

while the expansion up to second order includes the additional leading-order
nonlinear term $t_{\mu\nu}$. Since we are interested in the resulting
propagation of GWs at the location of the observer, we can assume a flat
background space-time away from the source.

The second-order Ricci tensor is Maggiore (2007)

$\displaystyle R_{\mu\nu}^{(2)}$ $\displaystyle=$
$\displaystyle{\textstyle{1\over 2}}\big{(}{\textstyle{1\over
2}}\partial_{\mu}h_{\alpha\beta}\partial_{\nu}h^{\alpha\beta}+h^{\alpha\beta}\partial_{\mu}\partial_{\nu}h_{\alpha\beta}-h^{\alpha\beta}\partial_{\nu}\partial_{\beta}h_{\alpha\mu}$
(11) $\displaystyle-$ $\displaystyle
h^{\alpha\beta}\partial_{\mu}\partial_{\beta}h_{\alpha\nu}+h^{\alpha\beta}\partial_{\alpha}\partial_{\beta}h_{\mu\nu}+\partial^{\beta}h^{\alpha}_{\
\nu}\partial_{\beta}h_{\alpha\mu}$ $\displaystyle-$
$\displaystyle\partial^{\beta}h^{\alpha}_{\
\nu}\partial_{\alpha}h_{\beta\mu}-\partial_{\beta}h^{\alpha\beta}\partial_{\nu}h_{\alpha\mu}+\partial_{\beta}h^{\alpha\beta}\partial_{\alpha}h_{\mu\nu}$
$\displaystyle-$
$\displaystyle\partial_{\beta}h^{\alpha\beta}\partial_{\mu}h_{\alpha\nu}-{\textstyle{1\over
2}}\partial^{\alpha}h\,\partial_{\alpha}h_{\mu\nu}+{\textstyle{1\over
2}}\partial^{\alpha}h\,\partial_{\nu}h_{\alpha\mu}$ $\displaystyle+$
$\displaystyle{\textstyle{1\over
2}}\partial^{\alpha}h\partial_{\mu}h_{\alpha\nu}\bigr{)}.$

We choose the traceless and harmonic gauge, such that $h=0$ (and hence,
$\bar{h}_{\mu\nu}=h_{\mu\nu}$, so we can drop the overbar from now on) and
$\partial^{\mu}h_{\mu\nu}=0$.

The energy density carried by GWs corresponds to the curvature they induce on
the background so we want to identify the terms of the LL pseudotensor that
contribute to the background metric by averaging over the high frequencies
that characterize the strain perturbations $h_{\mu\nu}$ compared to the
background $\bar{g}_{\mu\nu}$. Therefore, Eq. (11) gets simplified to (see
Ref. Maggiore (2007) for details)

$\langle t_{\mu\nu}\rangle=-\frac{1}{\kappa}\langle
R_{\mu\nu}^{(2)}\rangle=\frac{1}{4\kappa}\langle\partial_{\mu}h_{\alpha\beta}\partial_{\nu}h^{\alpha\beta}\rangle,$
(12)

where the angle brackets refer to the aforementioned average. If we require
$k^{i}\tilde{h}_{i}=0$, then the strains are expressed in the traceless-
transverse (TT) gauge, such that the gauge-invariant 00-component corresponds
to the GW energy density Arnowitt  _et al._ (1961),

$\langle t_{00}\rangle=\frac{1}{4\kappa}\langle\partial_{t}h_{ij}^{\rm
TT}\partial_{t}h_{ij}^{\rm TT}\rangle,$ (13)

and the spatial components are

$\langle t_{ij}\rangle=\frac{1}{4\kappa}\langle\partial_{i}h_{lm}^{\rm
TT}\partial_{j}h_{lm}^{\rm TT}\rangle.$ (14)

In the present context of continuous sources during cosmological epochs, we
consider the leading-order nonlinear GW equation, given in Eq. (8), and choose
the term that survives the average over high frequencies from all the
nonlinear terms in Eq. (11), which corresponds to the memory effect Will and
Wiseman (1996); Favata (2009a, 2010),

$\Box h_{ij}^{\rm nlin}=-2\kappa(T_{ij}+t_{ij}),$ (15)

where the “nlin” superscript refers to the GW strains sourced by the nonlinear
terms up to leading-order and

$t_{ij}=\frac{1}{4\kappa}\partial_{i}h_{lm}^{\rm TT}\partial_{j}h_{lm}^{\rm
TT}.$ (16)

This is the equation, adapted to an expanding universe, that we use in our
study of the early universe in Sec. III.3. We continue using the term “memory
effect,” even though there may not be any obvious permanent displacement of
freely-falling test masses due to the GWB.

###  III.2 Polarization modes of nonlinear solution

We express the linear polarization modes of the GW strain in Fourier space as
$\tilde{h}_{+,\times}(t,\bm{k})=\tilde{e}_{+,\times}^{\,ij}(\bm{k})\tilde{h}_{ij}(t,\bm{k})$
Varshalovich  _et al._ (1988). The $+$ and $\times$ modes are obtained as a
result of the TT projection, and $\tilde{e}_{+,\times}^{\,ij}$ are the basis
tensors,

$\tilde{e}_{+}^{\,ij}=\tilde{e}_{i}^{1}\tilde{e}_{j}^{1}-\tilde{e}_{i}^{2}\tilde{e}_{j}^{2},\quad\tilde{e}_{\times}^{\,ij}=\tilde{e}_{i}^{1}\tilde{e}_{j}^{2}+\tilde{e}_{i}^{2}\tilde{e}_{j}^{1},$
(17)

related to the helical polarization basis, given in Eq. 4 for a vector field,
as $\tilde{e}^{\,ij}_{\pm}=(\tilde{e}_{+}^{\,ij}\pm
i\tilde{e}^{\,ij}_{\times})/\sqrt{2}$ Caprini  _et al._ (2004). In the same
way, we define $\tilde{T}_{+,\times}$ and $\tilde{t}_{+,\times}$ that source
the corresponding GW mode.

To study the polarization modes of $t_{ij}$, we take the WKB or eikonal
approximation, i.e., we approximate the solution $\tilde{h}_{\mu\nu}={\cal
C}_{\mu\nu}e^{ik_{\alpha}x^{\alpha}}\approx{\cal C}_{\mu\nu}e^{i\phi}$, being
${\cal C}_{\mu\nu}$ the amplitudes of the perturbations, and $\phi$ the
constant GW phase. The effective stress tensor of GWs, given in Eq. (16), can
be expressed as [see Eq. (4.5) of Ref. Isaacson (1968b)]

$t_{\mu\nu}^{\rm\,WKB}=\frac{{\cal
C}^{2}}{4\kappa}k_{\mu}k_{\nu}\sin^{2}\phi,$ (18)

where ${\cal C}^{2}={\cal C}^{\mu\nu}{\cal C}_{\mu\nu}$ and
$k^{\mu}=(\omega=k,k^{i})$, with $k^{i}$ being the characteristic wave number
and $\omega$ the angular frequency, or, equivalently, as

$t_{ij}^{\rm\,WKB}=t_{00}\hat{k}_{i}\hat{k}_{j}=\frac{1}{4\kappa}\bigl{(}\partial_{t}h_{ab}^{\rm
TT}\partial_{t}h_{ab}^{\rm TT}\bigr{)}\hat{k}_{i}\hat{k}_{j}.$ (19)

This expression is usually used in analytical studies of the memory effect;
see, e.g., Refs. Favata (2009b); Jenkins and Sakellariadou (2021). Note that
in our present work, we source the GWs directly from $t_{\mu\nu}$, such that
the WKB approximation is not necessary. However, Eq. (19) is useful since it
allows us to show that only one polarization mode is sourced by the effective
stress tensor. To see this, we Fourier transform Eq. (19) to get

$\tilde{t}_{ij}^{\rm\,WKB}=\frac{1}{4\kappa}{\cal
F}\bigl{(}\partial_{t}h_{ab}^{\rm TT}\partial_{t}h_{ab}^{\rm TT}\bigr{)}*{\cal
F}(\hat{k}_{i}\hat{k}_{j}),$ (20)

where $\cal F$ refers to the Fourier transform and is equivalent to previous
notation, which uses a tilde. One can show that ${\cal
F}(\hat{k}_{i}\hat{k}_{j})\propto\delta_{ij}$ and, hence,
$\tilde{t}_{ij}^{\rm\,WKB}$ only has one non-zero mode, since $\delta^{+}=1$
and $\delta^{\times}=0$. Hence, the GW Eq. (15) can be expressed as

$\Box\tilde{h}^{\rm
nlin}_{+}=-2\kappa(\tilde{T}_{+}+\tilde{t}_{+}),\quad\Box\tilde{h}^{\rm
nlin}_{\times}=-2\kappa\tilde{T}_{\times},$ (21)

and the polarization modes of the nonlinear strains are $\tilde{h}^{\rm
nlin}_{+}=\tilde{h}_{+}+\Delta\tilde{h}_{+}$ and $\tilde{h}^{\rm
nlin}_{\times}=\tilde{h}_{\times}$, where $\tilde{h}_{+,\times}$ are the
linear solutions sourced by $\tilde{T}_{+,\times}$ and $\Delta\tilde{h}_{+}$
is the additional strain sourced by $\tilde{t}_{+}$.

###  III.3 EM-sourced GWs in the early universe

For the early universe, we adopt conformal time $\eta$ such that
$a{\rm\,d}\eta={\rm\,d}t$ and set $\eta_{*}=H_{*}^{-1}$ as the end of
reheating time, where $H_{*}=\sqrt{\kappa\rho_{\rm crit}(\eta_{*})/3}$ is the
Hubble parameter at time $\eta_{*}$ with critical energy density $\rho_{\rm
crit}(\eta_{*})$.

We normalize the conformal time $\bar{\eta}=\eta/\eta_{*}$, such that
$-1<\bar{\eta}<1$ corresponds to the reheating era, and we use comoving and
normalized wave vector $\bar{\bm{k}}=\bm{k}/(aH_{*})$. Since we study the GW
relics at the end of reheating, i.e., at the onset of radiation-dominated era,
we approximate $\rho_{\rm crit}$ by the radiation energy density as $\rho_{\rm
crit}(\eta_{*})\simeq\mathcal{E}_{\rm rad}(\eta_{*})=\pi^{2}g_{\rm
r}(T_{*})k_{\rm B}^{4}T_{*}^{4}/(30\hbar^{3})$, where $g_{\rm r}(T_{*})$ is
the number of relativistic degrees of freedom at temperature $T_{*}$, $k_{\rm
B}$ is the Boltzmann constant, and $\hbar$ is the reduced Planck constant. We
also scale the strain $\bar{h}_{ij}=ah_{ij}$ and normalize the comoving stress
of the source $\bar{T}_{ij}=a^{4}T_{ij}/\rho_{\rm crit}$. Then dropping all
the overbars, the resulting GW equation Roper Pol  _et al._ (2020a), with the
leading-order nonlinear term added, reads111 In Eq. (21), we set
$\tilde{t}_{\times}=0$ based on the eikonal approximation. However, we compute
$\tilde{t}_{\times}$ in the numerical simulations, so we keep this term in Eq.
(22).

$\Bigl{(}\tilde{h}^{\rm
nlin}_{+,\times}\Bigr{)}^{\prime\prime}+\biggl{(}k^{2}-\frac{a^{\prime\prime}}{a}\biggr{)}\tilde{h}_{+,\times}^{\rm
nlin}=\frac{6}{a}\bigl{(}\tilde{T}_{+,\times}+\tilde{t}_{+,\times}\bigr{)},$
(22)

where $a^{\prime\prime}/a$ is given in Eq. (6). See Refs. Maggiore (2007);
Caprini and Figueroa (2018) for details and Ref. Roper Pol  _et al._ (2020a)
for the implementation in the Pencil Code Collaboration (2021). The EM stress
components in physical space are assembled as

$T_{ij}=f^{2}(B_{i}B_{j}+E_{i}E_{j}),$ (23)

where $B_{i}$ and $E_{i}$ are components of the magnetic field
$\bm{B}=\bm{\nabla}\times\bm{A}$ and electric field
$\bm{E}=-\partial\bm{A}/\partial\eta$.

In connection with potential observations, we obtain the present day values of
the GW energy and helicity spectrum in the forms $h_{0}^{2}\Omega_{\rm GW}$
and $h_{0}^{2}\Xi_{\rm GW}$, which are independent of the uncertainty in the
present day Hubble parameter
$H_{0}=100\,h_{0}{\rm\,km}{\rm\,s}^{-1}{\rm\,Mpc}^{-1}$, with $h_{0}\approx
0.7$ Aghanim  _et al._ (2020),

$\displaystyle\Omega_{\rm GW}(f)$ $\displaystyle=$
$\displaystyle\biggl{(}\frac{H_{*}}{H_{0}}\biggr{)}^{2}\biggl{(}\frac{a_{*}}{a_{0}}\biggr{)}^{4}kE_{\rm
GW}(k),$ (24) $\displaystyle\Xi_{\rm GW}(f)$ $\displaystyle=$
$\displaystyle\biggl{(}\frac{H_{*}}{H_{0}}\biggr{)}^{2}\biggl{(}\frac{a_{*}}{a_{0}}\biggr{)}^{4}kH_{\rm
GW}(k),$ (25)

where $H_{*}$ and $a_{*}$ are the Hubble parameter and scale factor at the end
of reheating, and $H_{0}$ and $a_{0}$ are their present day values. Here
$f=a_{*}H_{*}k/(2\pi a_{0})$ is the physical frequency today. The GW spectra
$E_{\rm GW}(k)$ and $H_{\rm GW}(k)$ are computed from the time derivatives of
the strains as

$\displaystyle E_{\rm GW}(k)$
$\displaystyle\,=\int\Bigl{(}\bigl{|}\dot{\tilde{h}}_{+}\bigr{|}^{2}+\bigl{|}\dot{\tilde{h}}_{\times}\bigr{|}^{2}\Bigr{)}k^{2}{\rm\,d}\Omega_{k},$
(26) $\displaystyle H_{\rm GW}(k)$ $\displaystyle\,=\int 2\,{\rm
Im}\,\Bigl{(}\dot{\tilde{h}}_{+}\dot{\tilde{h}}_{\times}^{*}\Bigr{)}k^{2}{\rm\,d}\Omega_{k},$
(27)

where $\Omega_{k}$ is the solid angle of the shell with size $k$.

###  III.4 Order-of-magnitude estimates

####  III.4.1 Effects on the Maxwell equation

Nonlinear effects at the order $\mathcal{O}(h^{2})$ can potentially enter GW
production at the order $\mathcal{O}(h)$ in two ways. The first way, as has
been shown in Sec. III.1, is to have a term $t_{ij}=\mathcal{O}(h^{2})$, given
in Eq. (16), directly on the RHS of the standard GW equation, i.e., Eq. (15),
thereby acting as an additional source. The second way is to see how
$\mathcal{O}(h^{2})$ changes the regular EM source and whether in turn the
modified EM source could also affect GW productions at the linear order
$\mathcal{O}(h)$.

For the latter, we use the first of Maxwell’s equations in curved spacetime

$\nabla_{\mu}F^{\mu\nu}=\frac{1}{\sqrt{-g}}\partial_{\mu}\bigl{(}\sqrt{-g}F^{\mu\nu}\bigr{)}=\mu_{0}j^{\nu}=0.$
(28)

We have set the current to vanish during the reheating era, where the relevant
Maxwell’s equations are in vacuum. This can be rewritten as

$\displaystyle\partial_{\mu}F^{\mu\nu}$ $\displaystyle=$
$\displaystyle-\dfrac{\partial_{\mu}\sqrt{-g}}{\sqrt{-g}}F^{\mu\nu}$ (29)
$\displaystyle=$
$\displaystyle-\biggl{(}\dfrac{1}{2}\partial_{\mu}(\bar{g}^{\rho\sigma}h_{\rho\sigma})+\mathcal{O}(h^{2})\biggr{)}F^{\mu\nu},$
$\displaystyle\Box A^{\nu}$ $\displaystyle=$
$\displaystyle\Big{(}\mathcal{O}(h)+\mathcal{O}(h^{2})\Big{)}\mathcal{O}(A),$
(30)

where we have used

$\frac{\partial_{\mu}\sqrt{-g}}{\sqrt{-g}}=\frac{1}{2}\bar{g}^{\rho\sigma}\partial_{\mu}g_{\rho\sigma}=\frac{1}{2}\partial_{\mu}(\bar{g}^{\rho\sigma}h_{\rho\sigma})+\mathcal{O}(h^{2}),$
(31)

and applied the Lorenz gauge condition $\partial_{\mu}A^{\mu}=0$, such that
$\partial_{\mu}F^{\mu\nu}=\Box A^{\nu}$.

Therefore, Eq. (30) shows that the leading-order nonlinearity in the GWs, of
the order $\mathcal{O}(h^{2})$, does modify the four-potential, which
registers a difference

$\Delta A=\mathcal{O}(h^{2})\mathcal{O}(A).$ (32)

However, the result of such additional EM amplitude $\mathcal{E}_{\rm EM}$,
which is proportional to $A^{2}$; see Eq. (23), introduces a modification to
the production of GWs at $h={\cal O}(\mathcal{E}_{\rm
EM})=\mathcal{O}(A^{2})$; see Eq. (22), of larger order than quadratic in $h$.
For this reason, we only consider the additional GWs produced by their self-
backreaction at $t_{ij}=\mathcal{O}(h^{2})$, rather than the backreaction from
the vector potential at $\mathcal{O}((\Delta
A)^{2})=\mathcal{O}(h^{4})\mathcal{O}(A^{2})=\mathcal{O}(h^{5})$.

####  III.4.2 GW energy

The standard GW equation is linear in the strain and its sourcing energy,
i.e., $h=\mathcal{O}(\mathcal{E}_{\rm EM})$. Since GW energy is quadratic in
strain, i.e., $\mathcal{E}_{\rm GW}=\mathcal{O}(h^{2})$, we have a quadratic
relation

$\mathcal{E}_{\rm GW}=\mathcal{O}(\mathcal{E}_{\rm EM}^{2}),$ (33)

which has been verified by a number of earlier simulations Roper Pol  _et al._
(2020b); Brandenburg  _et al._ (2021b, c); Roper Pol  _et al._ (2021). Adding
nonlinearities of the order $\mathcal{O}(h^{2})$ to the source means that
$h^{\rm nlin}=h+\mathcal{O}(h^{2})$. Therefore, the nonlinear GW energy
relates to the source energies as

$\displaystyle\mathcal{E}_{\rm GW}^{\rm nlin}$
$\displaystyle=\bigl{[}(\mathcal{E}_{\rm
GW})^{1/2}+\mathcal{O}(h^{2})\bigr{]}^{2}$ $\displaystyle=\mathcal{E}_{\rm
GW}+\mathcal{E}_{\rm GW}^{1/2}\mathcal{O}(h^{2})+{\cal O}(h^{4})$
$\displaystyle=\mathcal{E}_{\rm GW}+\mathcal{O}(\mathcal{E}_{\rm
EM}^{3})+{\cal O}(\mathcal{E}_{\rm EM}^{4}),$ (34)

which implies that the leading-order energy difference $\Delta\mathcal{E}_{\rm
GW}\equiv\mathcal{E}_{\rm GW}^{\rm nlin}-\mathcal{E}_{\rm GW}$ should be
proportional to the cube of the sourcing energy, i.e.,

$\Delta\mathcal{E}_{\rm GW}=\mathcal{O}(\mathcal{E}_{\rm EM}^{3}).$ (35)

This relation will be verified and quantified numerically in Sec. IV.3.

####  III.4.3 GW polarization

The GW polarization is defined to be Kahniashvili  _et al._ (2005)

$\mathcal{P}_{\rm GW}(k)=\frac{H_{\rm GW}(k)}{E_{\rm GW}(k)}.$ (36)

The differences in GW spectra due to the nonlinear solution; see Eq. (21), are

$\displaystyle E_{\rm GW}^{\rm nlin}(k)-E_{\rm GW}(k)\propto$
$\displaystyle\,\Delta\tilde{h}_{+}\tilde{h}_{+}^{*}+\tilde{h}_{+}\Delta\tilde{h}_{+}^{*}$
$\displaystyle=$ $\displaystyle\,{\cal O}(h^{3})={\cal O}(\mathcal{E}_{\rm
EM}^{3}),$ (37) $\displaystyle H_{\rm GW}^{\rm nlin}(k)-H_{\rm GW}(k)\propto$
$\displaystyle\,\Delta\tilde{h}_{+}\tilde{h}_{\times}^{*}-\Delta\tilde{h}_{+}^{*}\tilde{h}_{\times}$
$\displaystyle=$ $\displaystyle\,{\cal O}(h^{3})={\cal O}(\mathcal{E}_{\rm
EM}^{3}).$ (38)

Note that Eq. (37) is consistent with Eq. (34). We can now express the
relative differences in the GW polarization due to the addition of the
nonlinear term (omitting the $k$ dependence to simplify the notation) as

$\displaystyle\frac{\mathcal{P}_{\rm GW}-\mathcal{P}_{\rm GW}^{\rm
nlin}}{\mathcal{P}_{\rm GW}}=$ $\displaystyle\,1-\frac{H_{\rm GW}^{\rm
nlin}E_{\rm GW}}{H_{\rm GW}E_{\rm GW}^{\rm nlin}}$ $\displaystyle=$
$\displaystyle\,1-\frac{1+H_{\rm GW}^{-1}\,{\cal O}(h^{3})}{1+E_{\rm
GW}^{-1}\,{\cal O}(h^{3})}$ $\displaystyle=$ $\displaystyle\,1-(1+H_{\rm
GW}^{-1}\,{\cal O}(h^{3}))(1-E_{\rm GW}^{-1}\,{\cal O}(h^{3}))$
$\displaystyle=$ $\displaystyle\,E_{\rm GW}^{-1}\,{\cal O}(h^{3})-H_{\rm
GW}^{-1}\,{\cal O}(h^{3})$ $\displaystyle =$ $\displaystyle\,{\cal
O}(h)-\mathcal{P}_{\rm GW}^{-1}{\cal O}(h)$ $\displaystyle=$
$\displaystyle\,{\cal O}(h)={\cal O}(\mathcal{E}_{\rm EM}).$ (39)

We then expect a decrease in polarization due to the nonlinear contributions
proportional to the EM energy density and to the polarization obtained in the
linearized approach, which will be verified in Sec. IV.2. We find that the
sign of the difference and its exact value depend on the balance between the
two terms of order ${\cal O}(h)$, which correspond to those given in Eqs. (37)
and (38) divided by $H_{\rm GW}(k)$ and $E_{\rm GW}(k)$, respectively.

##  IV Numerical simulations

Tables 1 and 2 summarize the relevant simulation parameters and their final
output. In this work, all simulations are performed on $n^{3}=512^{3}$ mesh
points with the smallest wave number $k_{1}=1$, except in series E, where
$k_{1}=0.2$. The EM energy spectra peak at wave numbers $k_{*}$ within the
simulation domain, i.e., $k_{1}<k_{*}<k_{\rm Ny}$, where $k_{\rm Ny}=nk_{1}/2$
is the Nyquist wave number. We conduct five series of simulations: two
nonhelical series (A and B) and three helical series (C, D, and E).

Table 1: Choice of model parameters, initial magnetic spectrum for $k\leq k_{*}$, and energy dilution factors. Series | $T_{\rm r}[{\rm\,GeV}]$ | $\gamma$ | $\beta$ | $g_{\rm r}$ | $E_{\rm M}(k)$ | $(H_{*}/H_{0})^{2}(a_{*}/a_{0})^{4}$
---|---|---|---|---|---|---
A | 100 | 0 | 7.3 | 106 | $\propto k^{3}$ | $1.6\times 10^{-5}$
B | 0.15 | 0 | 2.7 | 15 | $\propto k^{3}$ | $3.1\times 10^{-5}$
C | 8 | 1 | 7.3 | 86 | $\propto k^{3}$ | $1.7\times 10^{-5}$
D | 0.12 | 1 | 2.7 | 20 | $\propto k^{3}$ | $2.8\times 10^{-5}$
E | $3\times 10^{5}$ | 1 | 1.7 | 106 | $\propto k^{5}$ | $1.6\times 10^{-5}$
Table 2: Summary of simulations and relevant quantities. A’, B’, C’, D’, and E’ runs adopt the same parameters as A, B, C, D, and E, respectively, but are run up to $\eta=10$ and we show their saturated values. For the rest of the runs, the values correspond to $\eta=1$. In all cases, $k_{*}$ is shown at $\eta=1$. Runs |  $B_{0}$ | $\mathcal{E}_{\rm EM}$ | $k_{*}(1)$ |  $\mathcal{E}_{\rm GW}$ | $\hskip 5.69054pt\Delta\mathcal{E}_{\rm GW}$ | $\Delta\mathcal{E}_{\rm GW}/\mathcal{E}_{\rm GW}^{\rm nlin}$ |  $\mathcal{P}_{\rm GW}^{\rm nlin}$ | $\Delta\mathcal{P}_{\rm GW}$ | $\Delta\mathcal{P}_{\rm GW}/\mathcal{P}_{\rm GW}$
---|---|---|---|---|---|---|---|---|---
A1 | $3.3\times 10^{-19}$ | 0.02 | $7.5$ | $1.2\times 10^{-5}$ | $2.7\times 10^{-10}$ | $2.3\times 10^{-5}$ | $-0.108$ | $1.8\times 10^{-7}$ | $-1.6\times 10^{-6}$
A1’ | $3.3\times 10^{-19}$ | 0.02 | $7.5$ | $7.7\times 10^{-6}$ | $1.6\times 10^{-9}$ | $2.1\times 10^{-4}$ | $-0.108$ | $2.5\times 10^{-4}$ | $-2.3\times 10^{-3}$
A2 | $7.5\times 10^{-19}$ | 0.1 | $7.5$ | $3.2\times 10^{-4}$ | $3.8\times 10^{-8}$ | $1.2\times 10^{-4}$ | $-0.108$ | $9.2\times 10^{-7}$ | $-8.6\times 10^{-6}$
A2’ | $7.5\times 10^{-19}$ | 0.1 | $7.5$ | $2.1\times 10^{-4}$ | $2.7\times 10^{-7}$ | $1.3\times 10^{-3}$ | $-0.109$ | $1.3\times 10^{-3}$ | $-1.2\times 10^{-2}$
A3 | $2.4\times 10^{-18}$ | 1 | $7.5$ | $3.4\times 10^{-2}$ | $4.0\times 10^{-5}$ | $1.2\times 10^{-3}$ | $-0.108$ | $9.6\times 10^{-6}$ | $-8.9\times 10^{-5}$
A4 | $7.5\times 10^{-18}$ | 10 | $7.5$ | $3.2\times 10^{0}$ | $3.8\times 10^{-2}$ | $1.2\times 10^{-2}$ | $-0.108$ | $1.0\times 10^{-4}$ | $-9.3\times 10^{-4}$
B1 | $2.7\times 10^{-7}$ | 0.02 | $2.9$ | $1.1\times 10^{-4}$ | $9.5\times 10^{-8}$ | $8.7\times 10^{-4}$ | $-0.405$ | $6.2\times 10^{-5}$ | $-1.5\times 10^{-4}$
B1’ | $2.7\times 10^{-7}$ | 0.02 | $2.9$ | $7.9\times 10^{-5}$ | $1.8\times 10^{-7}$ | $2.2\times 10^{-3}$ | $-0.380$ | $1.8\times 10^{-3}$ | $-4.8\times 10^{-3}$
B2 | $6.0\times 10^{-7}$ | 0.1 | $2.9$ | $2.6\times 10^{-3}$ | $1.1\times 10^{-5}$ | $4.3\times 10^{-3}$ | $-0.406$ | $3.0\times 10^{-4}$ | $-7.5\times 10^{-4}$
B2’ | $6.0\times 10^{-7}$ | 0.1 | $2.9$ | $1.9\times 10^{-3}$ | $2.9\times 10^{-5}$ | $1.5\times 10^{-2}$ | $-0.386$ | $8.2\times 10^{-3}$ | $-2.2\times 10^{-2}$
B3 | $1.9\times 10^{-6}$ | 1 | $2.9$ | $2.7\times 10^{-1}$ | $1.2\times 10^{-2}$ | $4.4\times 10^{-2}$ | $-0.408$ | $2.8\times 10^{-3}$ | $-6.9\times 10^{-3}$
B4 | $6.0\times 10^{-6}$ | 10 | $2.9$ | $2.6\times 10^{1}$ | $1.9\times 10^{1}$ | $4.2\times 10^{-1}$ | $-0.411$ | $5.5\times 10^{-3}$ | $-1.4\times 10^{-2}$
C1 | $1.7\times 10^{-24}$ | 0.02 | $17.8$ | $2.4\times 10^{-6}$ | $7.2\times 10^{-10}$ | $3.0\times 10^{-4}$ | $0.942$ | $3.9\times 10^{-5}$ | $4.1\times 10^{-5}$
C1’ | $1.7\times 10^{-24}$ | 0.02 | $17.8$ | $3.7\times 10^{-6}$ | $1.1\times 10^{-10}$ | $2.9\times 10^{-5}$ | $0.971$ | $-1.0\times 10^{-6}$ | $-1.0\times 10^{-6}$
C2 | $3.9\times 10^{-24}$ | 0.1 | $17.8$ | $6.7\times 10^{-5}$ | $1.0\times 10^{-7}$ | $1.6\times 10^{-3}$ | $0.942$ | $2.0\times 10^{-4}$ | $2.2\times 10^{-4}$
C2’ | $3.9\times 10^{-24}$ | 0.1 | $17.8$ | $1.0\times 10^{-4}$ | $1.7\times 10^{-8}$ | $1.7\times 10^{-4}$ | $0.971$ | $3.0\times 10^{-6}$ | $3.1\times 10^{-6}$
C3 | $1.2\times 10^{-23}$ | 1 | $17.8$ | $6.0\times 10^{-3}$ | $9.0\times 10^{-5}$ | $1.5\times 10^{-2}$ | $0.940$ | $2.0\times 10^{-3}$ | $2.1\times 10^{-3}$
C4 | $3.9\times 10^{-23}$ | 10 | $17.8$ | $6.7\times 10^{-1}$ | $1.2\times 10^{-1}$ | $1.5\times 10^{-1}$ | $0.918$ | $2.4\times 10^{-2}$ | $2.6\times 10^{-2}$
D1 | $2.4\times 10^{-9}$ | 0.02 | $6.7$ | $1.6\times 10^{-5}$ | $4.1\times 10^{-8}$ | $2.5\times 10^{-3}$ | $0.944$ | $2.8\times 10^{-4}$ | $3.0\times 10^{-4}$
D1’ | $2.4\times 10^{-9}$ | 0.02 | $6.7$ | $2.7\times 10^{-5}$ | $1.6\times 10^{-8}$ | $6.0\times 10^{-4}$ | $0.972$ | $8.0\times 10^{-5}$ | $8.2\times 10^{-5}$
D2 | $5.4\times 10^{-9}$ | 0.1 | $6.7$ | $4.2\times 10^{-4}$ | $5.4\times 10^{-6}$ | $1.3\times 10^{-2}$ | $0.943$ | $1.5\times 10^{-3}$ | $1.5\times 10^{-3}$
D2’ | $5.4\times 10^{-9}$ | 0.1 | $6.7$ | $6.9\times 10^{-4}$ | $2.3\times 10^{-6}$ | $3.4\times 10^{-3}$ | $0.972$ | $6.2\times 10^{-4}$ | $6.4\times 10^{-4}$
D3 | $1.7\times 10^{-8}$ | 1 | $6.7$ | $4.1\times 10^{-2}$ | $5.8\times 10^{-3}$ | $1.2\times 10^{-1}$ | $0.927$ | $1.7\times 10^{-2}$ | $1.8\times 10^{-2}$
D4 | $5.4\times 10^{-8}$ | 10 | $6.7$ | $4.2\times 10^{0}$ | $2.4\times 10^{1}$ | $8.5\times 10^{-1}$ | $0.672$ | $2.7\times 10^{-1}$ | $2.9\times 10^{-1}$
E1 | $4.5\times 10^{-6}$ | 0.02 | $4.3$ | $5.3\times 10^{-5}$ | $3.4\times 10^{-7}$ | $6.3\times 10^{-3}$ | $0.909$ | $8.1\times 10^{-4}$ | $8.9\times 10^{-4}$
E1’ | $4.5\times 10^{-6}$ | 0.02 | $4.3$ | $9.4\times 10^{-5}$ | $3.4\times 10^{-7}$ | $3.6\times 10^{-3}$ | $0.941$ | $4.6\times 10^{-4}$ | $4.8\times 10^{-4}$
E2 | $1.0\times 10^{-5}$ | 0.1 | $4.3$ | $1.3\times 10^{-3}$ | $4.2\times 10^{-5}$ | $3.2\times 10^{-2}$ | $0.905$ | $4.4\times 10^{-3}$ | $4.9\times 10^{-3}$
E2’ | $1.0\times 10^{-5}$ | 0.1 | $4.3$ | $2.3\times 10^{-3}$ | $4.5\times 10^{-5}$ | $1.9\times 10^{-2}$ | $0.938$ | $3.3\times 10^{-3}$ | $3.5\times 10^{-3}$
E3 | $3.2\times 10^{-5}$ | 1 | $4.3$ | $1.4\times 10^{-1}$ | $7.3\times 10^{-2}$ | $3.5\times 10^{-1}$ | $0.805$ | $1.0\times 10^{-1}$ | $1.1\times 10^{-1}$

In Table 1, we list the relevant model parameters for all series of the
nonhelical ($\gamma=0$) and helical ($\gamma=1$) magnetogenesis scenarios
considered in the present work. Note that, for the same values of $\beta$,
despite the different corresponding temperatures $T_{\rm r}$, due to helicity,
the energy dilution factors $(H_{*}/H_{0})^{2}(a_{*}/a_{0})^{4}$ remain
similar. This is because the dilution factor is proportional to $g_{\rm
r}^{-1/3}$.

There are five variants for series A, B, C, and D each, corresponding to final
EM energy $\mathcal{E}_{\rm EM}\in\\{0.02,0.1,1,10\\}$, and four variants for
series E, where $\mathcal{E}_{\rm EM}=10$ is absent since the GW solution
diverged in this case. In reality, we expect the EM energy density to be
$\mathcal{E}_{\rm EM}\lesssim 0.1$, which is generally believed to be the
upper limit on the additional relativistic components based on the abundance
of light elements imposed by the constraints of Big Bang Nucleosynthesis (BBN)
Shvartsman (1969); Grasso and Rubinstein (1996). However, we also take
unrealistically large values in the present work to verify the scaling of the
memory effect.

The current study mainly focuses on the nonlinear features of GWs during the
reheating era, when the sourcing magnetogenesis takes place concurrently.
However, we are also interested in the resulting nonlinear solutions up to our
present time to study their potential detectability. For this reason, we
choose additional runs per series with $\mathcal{E}_{\rm EM}=0.02$ and $0.1$,
and evolve them beyond the end of reheating at $\eta=1$. The simulation is
then continued up to $\eta=10$, but with the EM source being turned off at
$\eta=1$. In these longer runs, we observe that the resulting GW energy
density, polarization, and spectra, become oscillatory around stationary
solutions, just as in previous numerical simulations Roper Pol  _et al._
(2020b, 2021). Hence, we average the results over the oscillations in time
between $\eta=2$ and $10$ and assume the GW energy density to evolve
proportional to $a^{-4}$ up to our present time. The new runs are denoted as
A’, B’, C’, D’, and E’. To avoid any additional discontinuities at $\eta=1$,
we continue using $a=(\eta+1)^{2}/4$, i.e., no attempt is made to model the
subsequent radiation-dominated era. Note that $a^{\prime\prime}/a$ is 0.5 at
$\eta=1$ and decreases to zero as $\eta^{-2}$. All these runs produce a
certain change of the GW energy, but it is only weakly related to the effect
of $t_{ij}$.

For each variant, with fixed $\beta$ and $\mathcal{E}_{\rm EM}$, we then
obtain the total GW energy density $\mathcal{E}_{\rm GW}$ and polarization222
The total GW energy density and polarization are computed from the GW spectra
as $\mathcal{E}_{\rm GW}=\int E_{\rm GW}(k){\rm\,d}k$ and $\mathcal{P}_{\rm
GW}=\int H_{\rm GW}(k){\rm\,d}k/\int E_{\rm GW}(k){\rm\,d}k$.
$\mathcal{P}_{\rm GW}$ for both the linear and nonlinear cases. Table 2
summarizes the relevant simulation parameters and their output at $\eta=1$,
where the nonlinear effects are represented by $\Delta\mathcal{E}_{\rm
GW}\equiv\mathcal{E}_{\rm GW}^{\rm nlin}-\mathcal{E}_{\rm GW}$ and
$\Delta\mathcal{P}_{\rm GW}\equiv\mathcal{P}_{\rm GW}-\mathcal{P}_{\rm
GW}^{\rm nlin}$. For the continued runs, we show their values averaged between
$\eta=2$ and 10.

In the following, we introduce the relevant model parameters and initial
conditions in Sec. IV.1 and we present the nonlinear effects to the GW and
polarization spectra for different initial fields $B_{0}$ and reheating
parameters $\beta$ in Sec. IV.2. Their relation with the sourcing energy
$\mathcal{E}_{\rm EM}$ is parameterized in Sec. IV.3, their long-term
evolution is shown in Sec. IV.4, and the possibility for detection is
considered in Sec. IV.5. Only the results of Sec. IV.5 correspond to
observables at the present time, while previous results are shown in the
scaled and normalized units presented in Secs. II and III.

###  IV.1 Choice of parameters

Here we briefly explain the choice of parameters used for this work. The
values of $\alpha$ and $\beta$, governing the evolution of the coupling
function $f$; see Eq. (1), and the magnetic field initial conditions are all
chosen in line with Refs. Brandenburg and Sharma (2021); Brandenburg  _et al._
(2021a). For series A and B, we choose $\alpha=2$ to avoid the backreaction
problem Sharma  _et al._ (2017), and $\beta\in\\{7.3,2.7\\}$, corresponding to
reheating temperatures of $T_{\rm r}\in\\{100,0.15\\}{\rm\,GeV}$, i.e., the
energy levels of EW and QCD phase transitions. For helical scenarios, series C
and D carry the same value of $\alpha$ and $\beta$ values are maintained for
comparison with series A and B, although now they correspond to $T_{\rm
r}\in\\{8,0.12\\}{\rm\,GeV}$. The helical series E takes $\alpha=1$ to enable
a higher-than-EW reheating temperature and $\beta=1.7$, which corresponds to
$T_{\rm r}=3\times 10^{5}{\rm\,GeV}$; see Ref. Brandenburg  _et al._ (2021a)
for details. For series A, B, C, and D, the initial conditions at
$\eta=\eta_{\rm ini}=-0.99$ for Eq. (5) are random, Gaussian-distributed,
magnetic fields with strengths $B_{0}\ll 1$ and energy spectra

$E_{\rm M}(k)\propto\begin{cases}k^{3}\;\;&(k\leq k_{*}(\eta_{\rm ini})),\\\
k^{1-4\beta}\;\;&(k>k_{*}(\eta_{\rm ini})),\end{cases}$ (40)

and

$k_{*}(\eta)=\frac{2\beta}{\eta+1}\Bigg{(}\gamma+\sqrt{1+\gamma^{2}+\frac{1}{2\beta}}\Bigg{)}$
(41)

indicates the peak of the EM energy spectrum. For series E, we initialize with
$E_{\rm M}(k)\propto k^{5}$ for all wave numbers. However, Eq. (41) is still
valid for runs in series E at the end of reheating as they develop
$k_{*}(\eta)$ during their evolution.

###  IV.2 Dependence on $B_{0},\beta$, and helicity

To observe the roles of the initial field $B_{0}$, the reheating parameter
$\beta$, and the helicity $\gamma$, we show the GW energy density spectrum
$E_{\rm GW}(k)$ in Fig. 1, where each row of panels corresponds to one of the
series A to E. In the left column, the faintest curves correspond to
$\mathcal{E}_{\rm EM}=0.02$ whereas the darkest shades correspond to the
unrealistically extreme case of $\mathcal{E}_{\rm EM}=10$. Linear and
nonlinear solutions are shown as solid and dashed curves, respectively. Figure
1 shows the GW spectra obtained at $\eta=1$.

Figure 1: Left panels: energy spectra $E_{\rm GW}(k)$ for $\mathcal{E}_{\rm
EM}=0.02$ (faintest), 0.1, 1, and 10 (darkest) at $\eta=1$. Linear and
nonlinear solutions are in solid and dashed curves, respectively. Right
panels: differences in the energy spectra $|\Delta E_{\rm GW}(k)|$ with
$y$-axis to the left and ratio of the two spectra $|\Delta E_{\rm
GW}(k)|/E_{\rm GW}^{\rm nlin}(k)$ in faint dotted curves with $y$-axis to the
right.

From the left column of Fig. 1, we see that all linear solutions exhibit a
shallow spectrum at small $k$, with empirical slopes $\sim\\!k$ and
$\sim\\!k^{3/2}$ for nonhelical and helical cases, respectively. Around the GW
spectral peak, located at $\sim\\!2k_{*}$, the spectra drop by many orders of
magnitude, although the rate of such sharp drops varies – it is quicker for
larger values of $\beta$ and for helical cases, which is in agreement with
Ref. Brandenburg  _et al._ (2021a).

We then see that all the nonlinear solutions, shown as dashed curves in the
left column of Fig. 1, closely follow the linear ones, but depart at the steep
slopes when the linear spectra start to drop significantly. These departures
are in the form of an upward “bulge”, indicating that stronger GWs are
produced at the nonlinear level. The differences diminish towards larger wave
numbers. Helicity plays a visible role of extending the range of the bulges.
This can be seen by comparing Series A and C in Fig. 1, where the bulges span
over only a fraction of the intermediate $k$ range around $k\sim 10$ without
helicity but extend down to the largest wave numbers with helicity. This is
more pronounced in Series C, D, and E, where the top bulge now hovers above
its linear counterpart for almost the entire $k$ range.

In all cases, larger initial $B_{0}$ produce larger nonlinear effects. Since
the nonlinear features are not clearly visible in $E_{\rm GW}(k)$ for all
values of $B_{0}$, we show the difference in GW energy density between the
linear and nonlinear results, $\Delta E_{\rm GW}\equiv E_{\rm GW}^{\rm
nlin}-E_{\rm GW}$, and the relative difference, $\Delta E_{\rm GW}/E_{\rm
GW}^{\rm nlin}$, in the right panels of Fig. 1. We see that the bulges are
indeed the features that stand out, as they give the largest ratios. However,
whether they contribute significantly to the overall energy difference varies
among the series. In series D and E, as the bulges occur when the linear
spectra are shallow, the nonlinear effects contribute significantly to the
overall GW energy density. On the other hand, spectra for series A, B, and C
at small $k$ display no visible bulges, yet register relatively large values
of the difference $\Delta E_{\rm GW}$. In these runs, the bulge only appears
where the spectra are already dropping significantly, such that their effect
to the total GW energy density is negligible. At most wave numbers, we observe
the nonlinear spectrum to be larger than the linear one, as expected. However,
at large wave numbers, we find some values of $k$ at which the opposite
occurs, in particular in the nonhelical series A and B. We find that the time
evolution of the GW modes shows a growing phase followed by an oscillatory
one, as was previously found in other numerical simulations Roper Pol  _et
al._ (2020b, 2021). It is during the oscillatory phase that the linear
solution can become larger than the nonlinear one due to the fact that, even
though the latter has a larger oscillatory amplitude, there is a slight shift
in the phase of the oscillations, allowing the former to acquire larger values
between the minimum and the maximum of the oscillations.

Figure 2: GW polarization spectra for $\beta=7.3$ (top panel), $\beta=2.7$
(middle panel), and $\beta=1.7$ (bottom panel) at $\eta=1$. Solid and dashed
curves are linear and nonlinear solutions. Fainter curves indicate smaller
values of $\mathcal{E}_{\rm EM}$. In series D and E, the dashed curves for
$\mathcal{E}_{\rm EM}=0.02$, $0.1$, and 1 overlap almost entirely with the
solid piece and with each other.

In Fig. 2, we show the degree of circular polarization $\mathcal{P}_{\rm
GW}(k)$, defined in Eq. (36), for series of helical runs C, D, and E, with
$\beta=7.3$, $2.7$, and $1.7$, respectively. We observe that the nonlinearity
produces weaker polarization where the nonlinear strains are larger, which
occurs at large wave numbers, $k\sim 100$, i.e., where the bulges occur; see
series C and D in Fig. 1. Note that the differences in polarization seem to
affect all variants of series C with $\beta=7.3$, but are almost negligible
for $\beta=2.7$ apart from D4 with $\mathcal{E}_{\rm EM}=10$. This is also
comparable to the fractional effects $\Delta E_{\rm GW}(k)/E_{\rm GW}^{\rm
nlin}(k)$ shown in series C and D of Fig. 1, where all variants of series C
register significant nonlinear effects, but among series D only D4 does.
Similarly to series D, the runs of series E do not show a large decrease in
the polarization, only visible for E3, which corresponds to the largest EM
energy density considered in these runs, i.e., $\mathcal{E}_{\rm EM}=1$. This
is expected since the decrease in polarization is induced by the increased
ratio in the GW spectrum as shown in Eq. (39).

###  IV.3 Scaling with $\mathcal{E}_{\rm EM}$

In Eq. (35), we have alluded to the expected order-of-magnitude relation
between $\Delta\mathcal{E}_{\rm GW}$ and $\mathcal{E}_{\rm EM}$. Here we
present the empirical findings of such a relation.

We show in Fig. 3 the GW production efficiency by comparing the GW energy
reached at the end of reheating $\eta_{*}$, $\mathcal{E}_{\rm GW}$ (top
panel), and the energy difference, $\Delta\mathcal{E}_{\rm GW}$ (middle
panel), against the maximum sourcing energy, also reached at $\eta_{*}$,
$\mathcal{E}_{\rm EM}$. The gray shaded regions in both panels correspond to
values of EM energy densities larger than the upper bound imposed by BBN
Shvartsman (1969); Grasso and Rubinstein (1996); Kahniashvili  _et al._
(2011). We find that the quadratic relation $\mathcal{E}_{\rm
GW}\propto\mathcal{E}_{\rm EM}^{2}$ and the cubic relation
$\Delta\mathcal{E}_{\rm GW}\propto\mathcal{E}_{\rm EM}^{3}$ hold at all values
of $\mathcal{E}_{\rm EM}$ considered.

Figure 3: GW energy density $\mathcal{E}_{\rm GW}$ (top panel), energy
difference $\Delta\mathcal{E}_{\rm GW}$ (middle panel), and relative
polarization suppression $|\Delta\mathcal{P}_{\rm GW}/\mathcal{P}_{\rm GW}|$
(bottom panel) against sourcing energy density $\mathcal{E}_{\rm EM}$ at the
end of reheating. We show series A (blue hollow circles), series B (green
hollow circles), series C (yellow dots), series D (red dots), and series E
(purple dots); see Table 2. Dashed lines and crosses are the fittings
corresponding to the saturated values of the extended runs A’, B’, C’, D’, and
E’. Shaded areas indicate $\mathcal{E}_{\rm EM}>0.1$.

Then, we denote the efficiency coefficients as $q$ and $p$ and rewrite the
order-of-magnitude expressions in Eqs. (33) and (35) as

$\mathcal{E}_{\rm GW}=(q\mathcal{E}_{\rm EM})^{2},\;\Delta\mathcal{E}_{\rm
GW}=(p\mathcal{E}_{\rm EM})^{3}.$ (42)

We obtain the pairs of values $(q,p)$ that fit the runs of Table 2 for each of
the series A to E. These are shown in Table 3 and in the legend of Fig. 3. We
see that runs with $\beta=2.7$ overall produce GWs more efficiently than those
with $\beta=7.3$, with or without helicity. The helical runs with $\beta=1.7$
have even larger efficiency than the other helical runs. Hence, smaller values
of $\beta$ seem to induce GWs with larger efficiencies.

Table 3: Empirical findings of the coefficients $(q,p)$, $(\tilde{q},\tilde{p})$, $r$, and $\tilde{r}$ for the runs of Table 2. A, B, C, D, and E refer to the values at the end of reheating while A’, B’, C’, D’, and E’ refer to the saturated values.  Runs | $(q,p)$ | $(\tilde{q},\tilde{p})$ | $r$ | $\tilde{r}$
---|---|---|---|---
A | $(0.18,0.03)$ | $(1.36,0.25)$ | $8.5\times 10^{-5}$ | $6.4\times 10^{-4}$
A’ | $(0.14,0.06)$ | $(1.08,0.45)$ | $1.2\times 10^{-1}$ | $9.0\times 10^{-1}$
B | $(0.52,0.23)$ | $(1.53,0.68)$ | $7.0\times 10^{-3}$ | $2.1\times 10^{-2}$
B’ | $(0.44,0.29)$ | $(1.29,0.84)$ | $2.3\times 10^{-1}$ | $6.7\times 10^{-1}$
C | $(0.08,0.05)$ | $(1.42,0.84)$ | $2.0\times 10^{-3}$ | $3.6\times 10^{-2}$
C’ | $(0.10,0.03)$ | $(1.79,0.46)$ | $4.0\times 10^{-5}$ | $7.2\times 10^{-4}$
D | $(0.21,0.18)$ | $(1.41,1.20)$ | $1.5\times 10^{-2}$ | $1.0\times 10^{-1}$
D’ | $(0.27,0.13)$ | $(1.78,0.88)$ | $4.8\times 10^{-3}$ | $3.2\times 10^{-2}$
E | $(0.36,0.37)$ | $(1.54,1.58)$ | $4.5\times 10^{-2}$ | $1.9\times 10^{-1}$
E’ | $(0.48,0.36)$ | $(2.06,1.53)$ | $2.8\times 10^{-2}$ | $1.2\times 10^{-1}$

In some earlier parameterizations of GW efficiencies Roper Pol  _et al._
(2020b); Brandenburg  _et al._ (2021); Roper Pol  _et al._ (2021), it was
found advantageous to take the dependence on the typical wave number $k_{*}$
of Eq. (41) into account via $\tilde{q}\equiv qk_{*}$. Analogously, we might
write $\tilde{p}=k_{*}p$. Then Eq. (42) is rewritten as

$\mathcal{E}_{\rm GW}=(\tilde{q}\mathcal{E}_{\rm
EM}/k_{*})^{2},\;\Delta\mathcal{E}_{\rm GW}=(\tilde{p}\mathcal{E}_{\rm
EM}/k_{*})^{3}.$ (43)

The corresponding pairs $(\tilde{q},\tilde{p})$ are also listed in Table 3. We
see that, for both nonhelical and helical cases, $\tilde{q}$ and $\tilde{p}$
are indeed justified parameterizations, as they vary less across different
values of $\beta$ than $q$ and $p$.

In an analogous manner to the dependences of $\mathcal{E}_{\rm GW}$ and
$\Delta\mathcal{E}_{\rm GW}$ on $\mathcal{E}_{\rm EM}$, we also found in Eq.
(39) an expected linear dependence of the relative ratio of the decrease in
polarization due to the nonlinear contributions, i.e., $\Delta\mathcal{P}_{\rm
GW}/\mathcal{P}_{\rm GW}\propto\mathcal{E}_{\rm EM}$. This relation is also
found in the numerical simulations; see bottom panel of Fig. 3, especially for
lower values of $\mathcal{E}_{\rm EM}$. In this case, we define the
polarization suppression coefficients $r$ and $\tilde{r}\equiv rk_{*}$ such
that

$\Big{|}\frac{\Delta\mathcal{P}_{\rm GW}}{\mathcal{P}_{\rm
GW}}\Big{|}=r\,\mathcal{E}_{\rm EM},\quad\Big{|}\frac{\Delta\mathcal{P}_{\rm
GW}}{\mathcal{P}_{\rm GW}}\Big{|}=\tilde{r}\,\mathcal{E}_{\rm EM}/k_{*},$ (44)

and compute their values by fitting the results from the simulations, listed
in Table 2. The resulting values are shown in Table 3. Here $\tilde{r}$ is
also a somewhat better choice of parameterization than $r$.

###  IV.4 GW evolution after the end of reheating

Figure 4: Time evolution of the nonlinear GW energy density (upper panel) and
its difference with the linear one (bottom panel). We show runs A2’ (blue),
B2’ (green), C2’ (orange), D2’ (red), and E2’ (purple), all corresponding to
$\mathcal{E}_{\rm EM}=0.1$. The horizontal lines are the averaged values over
times larger than $\eta=2$ and the values at $\eta=1$. The vertical lines
correspond to $\eta=1$, i.e., the end of reheating, and $\eta=2$. The zoomed-
in plot (upper panel) shows the difference between nonlinear (dashed) and
linear solution for one of the runs. Positive and negative values in the
difference (lower panel) are shown using solid and dotted lines, respectively,
in the nonhelical runs, while helical runs only present positive values.

The long-term time evolution of GWs, in the absence of the EM source are
studied with runs A’, B’, C’, D’, and E’. The results for the runs with
$\mathcal{E}_{\rm EM}=0.1$ can be seen in Fig. 4, where both $\mathcal{E}_{\rm
GW}^{\rm nlin}$ (top) and $\Delta\mathcal{E}_{\rm GW}$ (bottom) become
oscillatory around a stationary value333We could analogously refer to the
linear energy density $\mathcal{E}_{\rm GW}$ since the difference between the
nonlinear and linear is at least 1.5 orders of magnitude smaller. after
$\eta\sim 2$. The resulting values of averaging over oscillations (i.e.,
between $\eta=2$ and 10) are recorded as the saturated values in Table 2.
Additionally, the parameterization coefficients are obtained analogously and
shown in Table 3. We observe that the saturated values of the nonlinear GW
energy density become slightly larger than those at the end of reheating in
the helical cases, while the opposite is found in the nonhelical runs. On the
other hand, the difference between nonlinear and linear solutions becomes
larger for nonhelical runs and smaller for helical ones. This is probably
related to the fact that the nonhelical runs present larger oscillations of
the solutions over time already in the stationary phase. As discussed in Sec.
IV.2, the linear solution might become larger than the nonlinear one at some
exceptional times due to a different phase in the time oscillations of the GW
modes, even though the latter presents larger oscillatory amplitudes. This
behavior can similarly occur in the total GW energy density, integrated over
wave numbers, as it is seen in Fig. 4 (see zoomed-in plot of the upper panel).
We observe that between the minimum and maximum of the oscillations, due to
the shift in phase, the nonlinear solution is retarded in time, allowing
smaller values with respect to the linear one. The negative values of the
differences (shown in the lower panel) always occur at this moment of the
oscillations, but the opposite is not necessarily true.

Besides the total GW energy density, the individual GW modes also enter an
oscillatory phase around stationary values. In Fig. 5, we show the linear and
nonlinear spectra, comparing those obtained at $\eta=1$ with the saturated
values. In all cases, the GW spectra decrease after $\eta=1$ at low wave
numbers, while they increase around the spectral peak. The latter increase is
less pronounced in the nonhelical runs, which then explains why in these runs,
the total saturated GW energy density becomes smaller than its value at
$\eta=1$. In the helical runs, the increase around the spectral peak
compensates for the decrease at low wave numbers, hence leading to a total
increase of the saturated values.

The differences between the linear and nonlinear GW spectra are not noticeable
by eye at $\eta=1$; see runs with $\mathcal{E}_{\rm EM}=0.1$ in Fig. 1.
However, the saturated values of the spectra show a pronounced bulge at large
wave numbers for the nonlinear solutions of runs A2’ and C2’, i.e., the runs
with $\beta=7.3$. Since this feature was observed in Fig. 1 for larger values
of $\mathcal{E}_{\rm EM}$, we can imply that it is enhanced by letting the GW
equation evolve for longer times after the end of reheating.

Figure 5: Nonlinear GW energy spectra. We show runs A2’ (blue), B2’ (green),
C2’ (orange), D2’ (red), and E2’ (purple), all corresponding to
$\mathcal{E}_{\rm EM}=0.1$. The solid thin lines correspond to the spectra at
the end of reheating, while the dashed lines to the saturated values of the
spectra. The inset shows the linear GW spectra near the spectral peak, where
both linear and nonlinear solutions are almost identical.

The polarization spectra of the helical runs is shown in Fig. 6. As in Fig. 2,
series D and E do not show a significant difference between linear and
nonlinear polarizations for $\mathcal{E}_{\rm EM}=0.1$. A noticeable
difference appears in series C, which already showed a strong decrease in
polarization at the end of reheating, from 1 down to values below 0.4 at wave
numbers between 50 and 150, which corresponds to the position of the GW energy
density bulge. In this case, the saturated polarization shows an even stronger
decrease, going all the way down to zero, i.e., becoming unpolarized, in the
same range of wave numbers. This again coincides with the appearance of a
stronger bulge in the GW energy density; see Fig. 5.

Figure 6: GW polarization spectra. We show helical runs C2’ (orange), D2’
(red), and E2’ (purple), all corresponding to $\mathcal{E}_{\rm EM}=0.1$. We
compare the linear (dotted) and nonlinear (dashed) saturated polarization
spectra with those at the end of reheating (thin solid curves).

###  IV.5 Comparison with observational limits

Table 2 and Fig. 1 have demonstrated that the nonlinear effects are
subdominant, and their spectra at the end of reheating are virtually
indistinguishable from the linear ones when $\mathcal{E}_{\rm EM}$ is smaller
than unity. This is still true when we consider the spectra at later times,
with the exception of the appearance of a bulge below the spectral peaks in
the cases with $\beta=7.3$. However, these features appear around values much
smaller than those near the spectral peak, making them difficult to be
observed in the future. We explore whether it is sensible to consider the
detection prospects of such nonlinear features, as well as the detectability
of the signals considered in the present work. To do this, we convert the
direct simulation results into present day observable GW energy density
$\Omega_{\rm GW}$ and helicity $\Xi_{\rm GW}$ spectra in terms of physical
frequency $f$, discussed in Sec. III.3. Figure 7 presents such estimated GW
energy densities

Figure 7:  GW energy $h_{0}^{2}\Omega_{\rm GW}$ (left) and helicity
$h_{0}^{2}\,\Xi_{\rm GW}(f)$ (right) spectra scaled to the present day, in
units of physical frequency $f$: run A2’ is in blue, B2’ in green, C2’ in
orange, D2’ in red, and E1’ and E2’ in purple. Solid curves are linear
solutions and dashed-dotted curves are the differences between nonlinear and
linear solutions. Dotted curves correspond to negative values of the
differences. We show the $2\sigma$ confidence region for the 30-frequency
power law fit for a common-process spectrum reported in the NANOGrav 12.5-year
dataset (blue shaded wedge) Arzoumanian  _et al._ (2020), as well as the
expected sensitivity curve for the SKA (black line) Moore  _et al._ (2015).
LISA instrument (dash-dotted line) Robson  _et al._ (2019) and power law
sensitivity (PLS) (solid line) Caprini  _et al._ (2019); Roper Pol  _et al._
(2021) are shown in green. The PLS of DECIGO Kawamura  _et al._ (2020) and of
the $\mu$Ares concept, proposed within the ESA Voyage 2050 programme Sesana
_et al._ (2021), are shown in dark red and cyan, respectively, assuming 4
years of observations. The PLSs estimated for astrometry methods such as Gaia
and Theia Garcia-Bellido  _et al._ (2021); for atomic interferometer projects
AION ($2{\rm\,km}$ design) Badurina  _et al._ (2020) and AEDGE El-Neaj  _et
al._ (2020), and using binary resonance forecasts by 2038 from binary
millisecond pulsars (MSPs) and lunar laser-ranging (LRR) measurements are
shown in dark blue, brown, and violet, respectively. Regarding the
detectability of GW helicity (right panel), we show the PLS of LISA to a
polarized GWB using the dipole response function induced by our proper motion
Domcke  _et al._ (2020); Ellis  _et al._ (2020); Roper Pol  _et al._ (2021)
and the one obtained using the LISA–Taiji network Orlando  _et al._ (2021);
Seto (2020); Roper Pol  _et al._ (2021). All PLS curves assume an SNR of 10.

$h_{0}^{2}\Omega_{\rm GW}(f)$ and $h_{0}^{2}\Xi_{\rm GW}(f)$ from runs A2’,
B2’, C2’, D2’, and E2’, listed in Table 2, which correspond to the upper bound
$\mathcal{E}_{\rm EM}=0.1$ imposed by the BBN limit Shvartsman (1969); Grasso
and Rubinstein (1996); Kahniashvili  _et al._ (2011), and E1’, with smaller
$\mathcal{E}_{\rm EM}=0.02$. We compare our numerical results with the
sensitivities of some current and planned GW detectors: space-based GW
detectors, like the Laser Interferometer Space Antenna (LISA) Amaro-Seoane
_et al._ (2017), the DECi-hertz Interferometer Gravitational wave Observatory
(DECIGO) Seto  _et al._ (2001), and the $\mu$Ares concept proposed within the
ESA Voyage 2050 programme Sesana  _et al._ (2021); PTAs, like the North
American Nanohertz Observatory for Gravitational Waves (NANOGrav) Arzoumanian
_et al._ (2020) and the Square Kilometer Array (SKA) sensitivity Moore  _et
al._ (2015); proposed astrometry methods using Gaia Moore  _et al._ (2017) and
Theia Boehm  _et al._ (2017); atomic interferometer designs, like the
terrestrial Atom Interferometer Observatory and Network (AION) Badurina  _et
al._ (2020) and the space-based Atomic Experiment for Dark Matter and Gravity
Exploration (AEDGE) El-Neaj  _et al._ (2020), and the proposed detections
based on binary resonance Blas and Jenkins (2021a, b).

For both helical and nonhelical magnetogeneses, $\beta=2.7$ (with $T_{\rm
r}\approx 100{\rm\,MeV}$) can produce GWs within the detection range of PTAs,
or astrometry methods using Gaia or Theia. However, since a primordial GWB is
not a transient event and will not allow detectors to compare linear and
nonlinear signals, very detailed spectral templates might be needed in order
to extract specific features regarding nonlinear GWs in future PTA data. The
helical run D2’ produces a linear GW solution that would be compatible with
the observations reported by NANOGrav, while the helical run B2’ would lead to
a larger GW signal, which has not been observed by current PTA measurements.
Therefore, $\mathcal{E}_{\rm EM}<0.1$ is required for nonhelical signals with
an end-of-reheating temperature around QCD. In both cases, with the
improvement of PTA measurements in the coming years, as well as astrometry
methods, we might be able to detect the corresponding nonlinear GW signals.

On the other hand, the runs with $\beta=7.3$, corresponding to $T_{\rm r}=100$
and $8{\rm\,GeV}$ for the nonhelical and helical cases, respectively, peak
around $\sim\\!10^{-5}{\rm\,Hz}$, which is located in between the frequencies
reached by PTA experiments and space-based GW detectors. There are some
proposals to detect GWs in this range of frequencies, like methods based on
binary resonance, which have already been used to put some upper constraints
on GWBs Blas and Jenkins (2021a, b), and next-generation of space-borne GW
detectors, like the concept proposed within the ESA Voyage 2050 programme,
$\mu$-Ares Sesana  _et al._ (2021). The former method, using forecasts of
lunar-laser ranging data available by 2038 (see Fig. 7) would allow to probe
GWs at $\lesssim$ 1 $\mu$Hz. However, the resulting signals from A2’ and C2’
runs are peaking at 10 and 5 $\mu$Hz, respectively. The expected sensitivity
of $\mu$Ares would allow to detect this type of signals with very large
signal-to-noise ratios (SNRs), as well as their nonlinear contribution
(especially for the nonhelical run A2’).

At larger energy scales, with $T_{\rm r}=3\times 10^{5}{\rm\,GeV}$, the
linearized contribution to the GW signal of the run E2’, with
$\mathcal{E}_{\rm EM}=0.1$, can be detectable by LISA with a large SNR $\sim
9500$, peaking at $6\times 10^{-2}{\rm\,Hz}$. In addition, the memory effect
can also be detectable with a SNR $\sim 580$. For lower EM energy density
$\mathcal{E}_{\rm EM}=0.02$ (E1’) we still find a detectable signal, with SNRs
of 380 and 4 for the linear and nonlinear contributions, respectively. Next-
generation planned space-based GW detectors like DECIGO and Big Bang Observer
(BBO) Crowder and Cornish (2005) (the latter is not shown in Fig. 7 but it is
expected to have a slightly lower sensitivity curve than DECIGO in a similar
range of frequencies) will improve even more the detectability of this type of
signals and their memory effect contributions. In addition, atomic
interferometry projects, like AION and AEDGE, will be able to detect the GW
signal of E2’. While the former will cover only the high-frequency range of
the GW signal and its nonlinear contribution, the sensitivity of the latter
will allow to probe the full signal.

Alternatively, the polarization features shown in Fig. 2 might offer a
different approach for observation. While PTAs are not capable of detecting
circular polarization of GWBs Belgacem and Kamionkowski (2020), different
methods to detect polarization have been proposed using space-based GW
detectors like LISA and Taiji Ruan  _et al._ (2020). For example, the
anisotropies induced in a polarized isotropic GWB by our proper motion produce
a dipolar response to the GW signal, which allows to detect their polarization
Seto (2006, 2007), as studied in Ref. Domcke  _et al._ (2020) for LISA and in
Refs. Seto (2020); Orlando  _et al._ (2021) for the LISA–Taiji network. The
detectability of polarized GWBs produced by primordial helical magnetic fields
originated or present at the EW scale, following this approach, has been
recently studied in Ref. Roper Pol  _et al._ (2021). We show in Fig. 7 the
helical GW spectra produced by our runs and compare them to the sensitivities
of LISA (via the dipole response) and the LISA–Taiji network. Only the runs of
series E’ peak within the range of frequencies where LISA can detect GW
signals. According to Ref. Roper Pol  _et al._ (2021), the SNR of 9500, which
has been obtained for run E2’, corresponds to an SNR of a fully polarized GWB
of $16$ assuming a flat spectrum. In our case, this run has an SNR of 8, which
is slightly below the commonly used value of 10 to claim the signal to be
detectable Caprini  _et al._ (2019). The difference in helicity due to the
memory effect corresponds to a signal (see Fig. 7) with an SNR of
$\sim\\!0.25$. By considering the LISA–Taiji network, the polarized SNR of a
flat spectrum would be 3000 following Ref. Roper Pol  _et al._ (2021). In our
specific case, the spectral shape of run E2’ yields a smaller value of 360.
The difference in helicity corresponds to an SNR of 12. For the run E1’, with
$\mathcal{E}_{\rm EM}=0.02$, we find SNRs of 18 and 0.18 for the linear and
nonlinear contributions, respectively. Hence, the polarization spectra
corresponding to the series E’, with EM energies between 0.02 and 0.1, could
be detectable, and the nonlinear contribution could be detectable around the
$\mathcal{E}_{\rm EM}=0.1$ case.

The amplitude of the polarization spectrum decreases when we include the
leading-order nonlinear term in the GW production, as shown in Fig. 2 and
quantified, using the total polarization, by the suppression coefficient $r$
(see Table 3 and bottom panel of Fig. 3). However, this suppression is only
significant for large values of the EM energy density $\mathcal{E}_{\rm
EM}>1$. Hence, the memory effect does not affect the prospects of polarization
detectability for $\mathcal{E}_{\rm EM}\leq 0.1$, which is the upper bound
imposed by the BBN limit.

##  V Conclusions

In this work, we take the GW equation beyond the standard linear order and
study the consequent spectral features due to such nonlinearities in the
context of cosmological GWBs driven by inflationary-generated primordial
magnetic fields during reheating era. Instead of providing a stationary piece
to the final strain for localized transient events such as binary black hole
mergers, the displacement memory effect in the context of a primordial
continuous source manifests in an overall boost in the GW energy and helicity.
In terms of energy spectra $E_{\rm GW}(k)$, the nonlinear boosts closely
follow the spectral shapes of the linear solutions, but become more evident by
bulging on top of the linear solutions, when the latter start decreasing
sharply. The bulges are also extended over a longer range of $k$ in the
presence of helicity. We find that the boost of GW energy density due to the
inclusion of the leading-order nonlinear term occurs at all wave numbers of
the spectrum and we predict analytically that it is due to terms of order
${\cal O}(h^{3})$; see Eqs. (34) and (37). However, we find a small shift in
the phase of the GW mode oscillations due to the nonlinear term. Hence, at
some specific times, this shift can compensate the nonlinear boost and smaller
values of the nonlinear solution can be found.

Another nonlinear aspect related to helicity can be seen in the polarization
spectra $\mathcal{P}_{\rm GW}(k)$ where, for $\beta=7.3$, the polarization is
weakened approximately at the $k$ range where the bulges occur. However, this
weakening is not visible for $\beta=2.7$, until the unrealistically large
$\mathcal{E}_{\rm EM}=10$ is reached. Similarly, in the case with $\beta=1.7$,
we find a slight decrease in the spectral polarization, which is only
significant in the case with $\mathcal{E}_{\rm EM}=1$. In general, we find
that the addition of the leading-order nonlinear term suppresses polarization
at all wave numbers and we show analytically that this decrease depends on the
polarization of the linearized solution and that it is proportional to terms
of the order ${\cal O}(h)$. This occurs because the nonlinear boost is, in
general, stronger in the GW energy density than in its helicity.

In terms of the total energy, the relation $\Delta\mathcal{E}_{\rm
GW}=(p\mathcal{E}_{\rm EM})^{3}$ is found. Together with the quadratic
relation $\mathcal{E}_{\rm GW}=(q\mathcal{E}_{\rm EM})^{2}$, found in earlier
similar numerical computations, the coefficients $q$ and $p$ are determined
numerically and given in Table 3. Furthermore, the scaled coefficients
$\tilde{p}\equiv k_{*}p$ and $\tilde{q}\equiv k_{*}q$ are found to be a better
parameterization, presenting less variation over different scenarios, which is
consistent with the scaling $\mathcal{E}_{\rm GW}=(q\mathcal{E}_{\rm
EM}/k_{*})^{2}$, also obtained in previous numerical computations. Similarly,
regarding the total GW polarization, we find the linear relation
$\Delta\mathcal{P}_{\rm GW}=r\mathcal{P}_{\rm GW}\mathcal{E}_{\rm EM}$, where
we define the polarization suppression coefficients $r$, given in Table 3.

Finally, the detection prospects of such nonlinear GWs are discussed. We find
that PTAs and astrometry methods are relatively promising in detecting certain
low end-of-reheating temperature magnetogenesis-induced GWs. However, it would
be rather challenging to detect nonlinear features directly, since they
correspond to a second-order addition to the linear GWBs. Hence, very accurate
models for the linear contribution are needed in order to remove it and
recover the nonlinear part. The signals produced considering end-of-reheating
temperatures between 8 and $100{\rm\,GeV}$ are the most challenging to detect,
due to the gap of proposed GW detectors around the $\mu$Hz band. The $\mu$Ares
concept proposed within the ESA Voyage 2050 programme Sesana  _et al._ (2021)
could allow to detect this type of signals. Space-borne GW detectors like
LISA, Taiji, and DECIGO, and atomic interferometers, also seem promising
candidates to detect higher temperature inflationary scenarios with expected
large SNR for the signals produced by helical magnetic fields. In addition,
the combination of a network of space detectors, e.g., LISA and Taiji, can
provide measurements of the polarization of such signals. Again, the specific
features of the nonlinear contribution are of second-order, so very accurate
models and predictions of the linear counterparts are required. Overall,
further studies are required for a more detailed understanding of nonlinear
GWs themselves, as well as the parameter space of various GW generation
processes.

## Data availability

The source code used for the simulations of this study, the Pencil Code, is
freely available Collaboration (2021). The simulation data are also available
at Ref. He  _et al._ . The calculations, the simulation data, and the routines
generating the plots are also available at
https://github.com/AlbertoRoper/GW_turbulence/tree/master/memory_effect.

###### Acknowledgements.

Support through the grant 2019-04234 from the Swedish Research Council
(Vetenskapsrådet) and the Shota Rustaveli National Science Foundation (SRNSF)
of Georgia (grant FR/18-1462) are gratefully acknowledged. ARP is supported by
the French National Research Agency (ANR) project MMUniverse
(ANR-19-CE31-0020). 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 the A9 allocation at
the Occigen supercomputer provided by the Grand Équipement National de Calcul
Intensif (GENCI) to the project “Opening new windows on Early Universe with
multi-messenger astronomy”.

## References

  * Abbott  _et al._ (2016a) B. P. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. Lett. 116, 061102 (2016a), arXiv:1602.03837 [gr-qc] .
  * Abbott  _et al._ (2016b) B. P. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. Lett. 116, 241103 (2016b), arXiv:1606.04855 [gr-qc] .
  * Abbott  _et al._ (2017) B. P. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. Lett. 119, 161101 (2017), arXiv:1710.05832 [gr-qc] .
  * Abbott  _et al._ (2019) B. P. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. X 9, 031040 (2019), arXiv:1811.12907 [astro-ph.HE] .
  * Abbott  _et al._ (2021) R. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. X 11, 021053 (2021), arXiv:2010.14527 [gr-qc] .
  * Zel’dovich and Polnarev (1974) Y. B. Zel’dovich and A. G. Polnarev, Sov. Astron. 18, 17 (1974).
  * Turner (1978) M. S. Turner, Nature 274, 565 (1978).
  * Epstein (1978) R. Epstein, Astrophys. J. 223, 1037 (1978).
  * Braginskii and Thorne (1987) V. B. Braginskii and K. S. Thorne, Nature 327, 123 (1987).
  * Christodoulou (1991) D. Christodoulou, Phys. Rev. Lett. 67, 1486 (1991).
  * Blanchet and Damour (1992) L. Blanchet and T. Damour, Phys. Rev. D 46, 4304 (1992).
  * Thorne (1992) K. S. Thorne, Phys. Rev. D 45, 520 (1992).
  * Bondi  _et al._ (1962) H. Bondi, M. G. J. van der Burg, and A. W. K. Metzner, Proc. Roy. Soc. Lond. A 269, 21 (1962).
  * Sachs (1962) R. K. Sachs, Proc. Roy. Soc. Lond. A 270, 103 (1962).
  * Strominger and Zhiboedov (2016) A. Strominger and A. Zhiboedov, JHEP 01, 086 (2016), arXiv:1411.5745 [hep-th] .
  * He  _et al._ (2015) T. He, V. Lysov, P. Mitra, and A. Strominger, JHEP 05, 151 (2015), arXiv:1401.7026 [hep-th] .
  * Weinberg (1965) S. Weinberg, Phys. Rev. 140, B516 (1965).
  * Lasky  _et al._ (2016) P. D. Lasky, E. Thrane, Y. Levin, J. Blackman, and Y. Chen, Phys. Rev. Lett. 117, 061102 (2016), arXiv:1605.01415 [astro-ph.HE] .
  * McNeill  _et al._ (2017) L. O. McNeill, E. Thrane, and P. D. Lasky, Phys. Rev. Lett. 118, 181103 (2017), arXiv:1702.01759 [astro-ph.IM] .
  * Hübner  _et al._ (2020) M. Hübner, C. Talbot, P. D. Lasky, and E. Thrane, Phys. Rev. D 101, 023011 (2020), arXiv:1911.12496 [astro-ph.HE] .
  * Boersma  _et al._ (2020) O. M. Boersma, D. A. Nichols, and P. Schmidt, Phys. Rev. D 101, 083026 (2020), arXiv:2002.01821 [astro-ph.HE] .
  * Amaro-Seoane  _et al._ (2017) P. Amaro-Seoane _et al._ (LISA), (2017), arXiv:1702.00786 [astro-ph.IM] .
  * Favata (2009a) M. Favata, Astrophys. J. Lett. 696, L159 (2009a), arXiv:0902.3660 [astro-ph.SR] .
  * Islo  _et al._ (2019) K. Islo, J. Simon, S. Burke-Spolaor, and X. Siemens, (2019), arXiv:1906.11936 [astro-ph.HE] .
  * Yang and Martynov (2018) H. Yang and D. Martynov, Phys. Rev. Lett. 121, 071102 (2018), arXiv:1803.02429 [gr-qc] .
  * Wang  _et al._ (2015) J. B. Wang _et al._ , Mon. Not. Roy. Astron. Soc. 446, 1657 (2015), arXiv:1410.3323 [astro-ph.GA] .
  * Arzoumanian  _et al._ (2015) Z. Arzoumanian _et al._ (NANOGrav), Astrophys. J. 810, 150 (2015), arXiv:1501.05343 [astro-ph.GA] .
  * Madison  _et al._ (2017) D. R. Madison, D. F. Chernoff, and J. M. Cordes, Phys. Rev. D 96, 123016 (2017), arXiv:1710.04974 [astro-ph.GA] .
  * Aggarwal  _et al._ (2020) K. Aggarwal _et al._ (NANOGrav), Astrophys. J. 889, 38 (2020), arXiv:1911.08488 [astro-ph.HE] .
  * Pollney and Reisswig (2011) D. Pollney and C. Reisswig, Astrophys. J. Lett. 732, L13 (2011), arXiv:1004.4209 [gr-qc] .
  * Khera  _et al._ (2021) N. Khera, B. Krishnan, A. Ashtekar, and T. De Lorenzo, Phys. Rev. D 103, 044012 (2021), arXiv:2009.06351 [gr-qc] .
  * Ebersold and Tiwari (2020) M. Ebersold and S. Tiwari, Phys. Rev. D 101, 104041 (2020), arXiv:2005.03306 [gr-qc] .
  * Mukhopadhyay  _et al._ (2021) M. Mukhopadhyay, C. Cardona, and C. Lunardini, (2021), arXiv:2105.05862 [astro-ph.HE] .
  * Burrows and Hayes (1996) A. Burrows and J. Hayes, Phys. Rev. Lett. 76, 352 (1996), arXiv:astro-ph/9511106 .
  * Segalis and Ori (2001) E. B. Segalis and A. Ori, Phys. Rev. D 64, 064018 (2001), arXiv:gr-qc/0101117 .
  * Sago  _et al._ (2004) N. Sago, K. Ioka, T. Nakamura, and R. Yamazaki, Phys. Rev. D 70, 104012 (2004), arXiv:gr-qc/0405067 .
  * Hiramatsu  _et al._ (2005) T. Hiramatsu, K. Kotake, H. Kudoh, and A. Taruya, Mon. Not. Roy. Astron. Soc. 364, 1063 (2005), arXiv:astro-ph/0509787 .
  * Aurrekoetxea  _et al._ (2020) J. C. Aurrekoetxea, T. Helfer, and E. A. Lim, Class. Quant. Grav. 37, 204001 (2020), arXiv:2002.05177 [gr-qc] .
  * Jenkins and Sakellariadou (2021) A. C. Jenkins and M. Sakellariadou, (2021), arXiv:2102.12487 [gr-qc] .
  * Kilicarslan and Tekin (2019) E. Kilicarslan and B. Tekin, Eur. Phys. J. C 79, 114 (2019), arXiv:1805.02240 [gr-qc] .
  * Hou (2021) S. Hou, Astron. Nachr. 342, 96 (2021), arXiv:2011.02087 [gr-qc] .
  * Hou and Zhu (2021) S. Hou and Z.-H. Zhu, JHEP 01, 083 (2021), arXiv:2005.01310 [gr-qc] .
  * Koyama (2020) K. Koyama, Phys. Rev. D 102, 021502 (2020), arXiv:2006.15914 [gr-qc] .
  * Shore (2018) G. M. Shore, JHEP 12, 133 (2018), arXiv:1811.08827 [gr-qc] .
  * Grishchuk and Polnarev (1989) L. P. Grishchuk and A. G. Polnarev, Sov. Phys. JETP 69, 653 (1989).
  * Zhang  _et al._ (2017) P. M. Zhang, C. Duval, G. W. Gibbons, and P. A. Horvathy, Phys. Lett. B 772, 743 (2017), arXiv:1704.05997 [gr-qc] .
  * Zhang  _et al._ (2018) P. M. Zhang, C. Duval, G. W. Gibbons, and P. A. Horvathy, JCAP 05, 030 (2018), arXiv:1802.09061 [gr-qc] .
  * Divakarla and Whiting (2021) A. K. Divakarla and B. F. Whiting, (2021), arXiv:2106.05163 [gr-qc] .
  * Pasterski  _et al._ (2016) S. Pasterski, A. Strominger, and A. Zhiboedov, JHEP 12, 053 (2016), arXiv:1502.06120 [hep-th] .
  * Nichols (2017) D. A. Nichols, Phys. Rev. D 95, 084048 (2017), arXiv:1702.03300 [gr-qc] .
  * Mitman  _et al._ (2020) K. Mitman, J. Moxon, M. A. Scheel, S. A. Teukolsky, M. Boyle, N. Deppe, L. E. Kidder, and W. Throwe, Phys. Rev. D 102, 104007 (2020), arXiv:2007.11562 [gr-qc] .
  * Tahura  _et al._ (2021) S. Tahura, D. A. Nichols, and K. Yagi, (2021), arXiv:2107.02208 [gr-qc] .
  * Roper Pol  _et al._ (2020a) A. Roper Pol, A. Brandenburg, T. Kahniashvili, A. Kosowsky, and S. Mandal, Geophys. Astrophys. Fluid Dynamics 114, 130 (2020a), arXiv:1807.05479 [physics.flu-dyn] .
  * Roper Pol  _et al._ (2020b) A. Roper Pol, S. Mandal, A. Brandenburg, T. Kahniashvili, and A. Kosowsky, Phys. Rev. D 102, 083512 (2020b), arXiv:1903.08585 [astro-ph.CO] .
  * Brandenburg and Sharma (2021) A. Brandenburg and R. Sharma, Astrophys. J. 920, 26 (2021), arXiv:2106.03857 [astro-ph.CO] .
  * Brandenburg  _et al._ (2021a) A. Brandenburg, Y. He, and R. Sharma, (2021a), arXiv:2107.12333 [astro-ph.CO] .
  * Sharma  _et al._ (2017) R. Sharma, S. Jagannathan, T. R. Seshadri, and K. Subramanian, Phys. Rev. D 96, 083511 (2017), arXiv:1708.08119 [astro-ph.CO] .
  * Sharma  _et al._ (2018) R. Sharma, K. Subramanian, and T. R. Seshadri, Phys. Rev. D 97, 083503 (2018), arXiv:1802.04847 [astro-ph.CO] .
  * Demozzi  _et al._ (2009) V. Demozzi, V. Mukhanov, and H. Rubinstein, JCAP 08, 025 (2009), arXiv:0907.1030 [astro-ph.CO] .
  * Collaboration (2021) P. C. Collaboration, J. Open Source Softw. 6, 2807 (2021), arXiv:2009.08231 [astro-ph.IM] .
  * Ratra (1992) B. Ratra, Astrophys. J. Lett. 391, L1 (1992).
  * Turner and Widrow (1988) M. S. Turner and L. M. Widrow, Phys. Rev. D 37, 2743 (1988).
  * Dolgov (1993) A. Dolgov, Phys. Rev. D 48, 2499 (1993), arXiv:hep-ph/9301280 .
  * Varshalovich  _et al._ (1988) D. A. Varshalovich, A. N. Moskalev, and V. K. Khersonsky,  _Quantum Theory of Angular Momentum: Irreducible Tensors, Spherical Harmonics, Vector Coupling Coefficients, 3nj Symbols_ (World Scientific, Singapore, 1988).
  * Isaacson (1968a) R. A. Isaacson, Phys. Rev. 166, 1263 (1968a).
  * Isaacson (1968b) R. A. Isaacson, Phys. Rev. 166, 1272 (1968b).
  * Misner  _et al._ (1973) C. W. Misner, K. S. Thorne, and J. A. Wheeler,  _Gravitation_ (W. H. Freeman, San Francisco, 1973).
  * Maggiore (2007) M. Maggiore,  _Gravitational Waves. Vol. 1: Theory and Experiments_ , Oxford Master Series in Physics (Oxford University Press, 2007).
  * Landau and Lifshitz (1962) L. D. Landau and E. M. Lifshitz,  _The classical theory of fields_ (Pergamon Press, Oxford, 1962).
  * Arnowitt  _et al._ (1961) R. L. Arnowitt, S. Deser, and C. W. Misner, Phys. Rev. 121, 1556 (1961).
  * Will and Wiseman (1996) C. M. Will and A. G. Wiseman, Phys. Rev. D 54, 4813 (1996), arXiv:gr-qc/9608012 .
  * Favata (2010) M. Favata, Class. Quant. Grav. 27, 084036 (2010), arXiv:1003.3486 [gr-qc] .
  * Caprini  _et al._ (2004) C. Caprini, R. Durrer, and T. Kahniashvili, Phys. Rev. D 69, 063006 (2004), arXiv:astro-ph/0304556 .
  * Favata (2009b) M. Favata, Phys. Rev. D 80, 024002 (2009b), arXiv:0812.0069 [gr-qc] .
  * Caprini and Figueroa (2018) C. Caprini and D. G. Figueroa, Class. Quant. Grav. 35, 163001 (2018), arXiv:1801.04268 [astro-ph.CO] .
  * Aghanim  _et al._ (2020) N. Aghanim _et al._ (Planck), Astron. Astrophys. 641, A6 (2020), [Erratum: Astron.Astrophys. 652, C4 (2021)], arXiv:1807.06209 [astro-ph.CO] .
  * Brandenburg  _et al._ (2021b) A. Brandenburg, Y. He, T. Kahniashvili, M. Rheinhardt, and J. Schober, Astrophys. J. 911, 110 (2021b), arXiv:2101.08178 [astro-ph.CO] .
  * Brandenburg  _et al._ (2021c) A. Brandenburg, E. Clarke, Y. He, and T. Kahniashvili, Phys. Rev. D 104, 043513 (2021c), arXiv:2102.12428 [astro-ph.CO] .
  * Roper Pol  _et al._ (2021) A. Roper Pol, S. Mandal, A. Brandenburg, and T. Kahniashvili, (2021), arXiv:2107.05356 [gr-qc] .
  * Kahniashvili  _et al._ (2005) T. Kahniashvili, G. Gogoberidze, and B. Ratra, Phys. Rev. Lett. 95, 151301 (2005), arXiv:astro-ph/0505628 .
  * Shvartsman (1969) V. F. Shvartsman, Pisma Zh. Eksp. Teor. Fiz. 9, 315 (1969).
  * Grasso and Rubinstein (1996) D. Grasso and H. R. Rubinstein, Phys. Lett. B 379, 73 (1996), arXiv:astro-ph/9602055 .
  * Kahniashvili  _et al._ (2011) T. Kahniashvili, A. G. Tevzadze, and B. Ratra, Astrophys. J. 726, 78 (2011), arXiv:0907.0197 [astro-ph.CO] .
  * Brandenburg  _et al._ (2021) A. Brandenburg, G. Gogoberidze, T. Kahniashvili, S. Mandal, A. R. Pol, and N. Shenoy, Classical and Quantum Gravity 38, 145002 (2021), arXiv:2103.01140 [gr-qc] .
  * Arzoumanian  _et al._ (2020) Z. Arzoumanian _et al._ (NANOGrav), Astrophys. J. Lett. 905, L34 (2020), arXiv:2009.04496 [astro-ph.HE] .
  * Moore  _et al._ (2015) C. J. Moore, R. H. Cole, and C. P. L. Berry, Class. Quant. Grav. 32, 015014 (2015), arXiv:1408.0740 [gr-qc] .
  * Robson  _et al._ (2019) T. Robson, N. J. Cornish, and C. Liu, Class. Quant. Grav. 36, 105011 (2019), arXiv:1803.01944 [astro-ph.HE] .
  * Caprini  _et al._ (2019) C. Caprini, D. G. Figueroa, R. Flauger, G. Nardini, M. Peloso, M. Pieroni, A. Ricciardone, and G. Tasinato, JCAP 11, 017 (2019), arXiv:1906.09244 [astro-ph.CO] .
  * Kawamura  _et al._ (2020) S. Kawamura _et al._ , (2020), arXiv:2006.13545 [gr-qc] .
  * Sesana  _et al._ (2021) A. Sesana _et al._ , Exper. Astron. 51, 1333 (2021), arXiv:1908.11391 [astro-ph.IM] .
  * Garcia-Bellido  _et al._ (2021) J. Garcia-Bellido, H. Murayama, and G. White, (2021), arXiv:2104.04778 [hep-ph] .
  * Badurina  _et al._ (2020) L. Badurina _et al._ , JCAP 05, 011 (2020), arXiv:1911.11755 [astro-ph.CO] .
  * El-Neaj  _et al._ (2020) Y. A. El-Neaj _et al._ (AEDGE), EPJ Quant. Technol. 7, 6 (2020), arXiv:1908.00802 [gr-qc] .
  * Domcke  _et al._ (2020) V. Domcke, J. Garcia-Bellido, M. Peloso, M. Pieroni, A. Ricciardone, L. Sorbo, and G. Tasinato, JCAP 05, 028 (2020), arXiv:1910.08052 [astro-ph.CO] .
  * Ellis  _et al._ (2020) J. Ellis, M. Fairbairn, M. Lewicki, V. Vaskonen, and A. Wickens, JCAP 10, 032 (2020), arXiv:2005.05278 [astro-ph.CO] .
  * Orlando  _et al._ (2021) G. Orlando, M. Pieroni, and A. Ricciardone, JCAP 03, 069 (2021), arXiv:2011.07059 [astro-ph.CO] .
  * Seto (2020) N. Seto, Phys. Rev. Lett. 125, 251101 (2020), arXiv:2009.02928 [gr-qc] .
  * Seto  _et al._ (2001) N. Seto, S. Kawamura, and T. Nakamura, Phys. Rev. Lett. 87, 221103 (2001), arXiv:astro-ph/0108011 .
  * Moore  _et al._ (2017) C. J. Moore, D. P. Mihaylov, A. Lasenby, and G. Gilmore, Phys. Rev. Lett. 119, 261102 (2017), arXiv:1707.06239 [astro-ph.IM] .
  * Boehm  _et al._ (2017) C. Boehm _et al._ (Theia), (2017), arXiv:1707.01348 [astro-ph.IM] .
  * Blas and Jenkins (2021a) D. Blas and A. C. Jenkins, (2021a), arXiv:2107.04063 [gr-qc] .
  * Blas and Jenkins (2021b) D. Blas and A. C. Jenkins, (2021b), arXiv:2107.04601 [astro-ph.CO] .
  * Crowder and Cornish (2005) J. Crowder and N. J. Cornish, Phys. Rev. D 72, 083005 (2005), arXiv:gr-qc/0506015 .
  * Belgacem and Kamionkowski (2020) E. Belgacem and M. Kamionkowski, Phys. Rev. D 102, 023004 (2020), arXiv:2004.05480 [astro-ph.CO] .
  * Ruan  _et al._ (2020) W.-H. Ruan, Z.-K. Guo, R.-G. Cai, and Y.-Z. Zhang, Int. J. Mod. Phys. A 35, 2050075 (2020), arXiv:1807.09495 [gr-qc] .
  * Seto (2006) N. Seto, Phys. Rev. Lett. 97, 151101 (2006), arXiv:astro-ph/0609504 .
  * Seto (2007) N. Seto, Phys. Rev. D 75, 061302 (2007), arXiv:astro-ph/0609633 .
  * (108) Y. He, A. Brandenburg, and A. Roper Pol,  10.5281/zenodo.5524454.