# Relic Cosmological Vector Fields and Inflationary Gravitational Waves

Avery J. Tishue<EMAIL_ADDRESS>Robert R. Caldwell  Department
of Physics and Astronomy, Dartmouth College, 6127 Wilder Laboratory, Hanover,
NH 03755

###### Abstract

We show that relic vector fields can significantly impact a spectrum of
primordial gravitational waves in the post-inflationary era. We consider a
triplet of U(1) fields in a homogeneous, isotropic configuration. The
interaction between the gravitational waves and the vector fields, from the
end of reheating to the present day, yields novel spectral features. The
amplitude, tilt, shape, and net chirality of the gravitational wave spectrum
are shown to depend on the abundance of the electric- and magnetic-like vector
fields. Our results show that even a modest abundance can have strong
implications for efforts to detect the imprint of gravitational waves on the
cosmic microwave background polarization. We find that a vector field
comprising less than $2\%$ of the energy density during the radiation
dominated era can have a greater than order unity effect on the predicted
inflationary gravitational wave spectrum.

##  I Introduction

The inflationary paradigm is widely regarded as the leading framework to
describe the early stages of the Hot Big Bang cosmology Guth (1981); Linde
(1982); Albrecht and Steinhardt (1982). Indeed, the inflationary prediction of
a nearly scale-invariant spectrum of primordial adiabatic density
perturbations Mukhanov and Chibisov (1981); Hawking (1982); Guth and Pi
(1982); Starobinsky (1982); Bardeen et al. (1983) underpins the successes of
the Standard Cosmological Model. Inflation also predicts a nearly scale-
invariant spectrum of primordial gravitational waves (GWs) Grishchuk (1975);
Starobinsky (1979); Rubakov et al. (1982); Fabbri and Pollock (1983); Abbott
and Wise (1984). Detection of the unique signature of these long wavelength
GWs in the polarization pattern of the cosmic microwave background (CMB) would
provide strong confirmation of inflation Kamionkowski et al. (1997); Seljak
and Zaldarriaga (1997). Consequently, there is an enormous experimental effort
to detect the polarization signal Ade et al. (2017, 2018, 2019); Akrami et al.
(2020); Tristram et al. (2020); Abazajian et al. (2016); Hazumi et al. (2019).
Similarly, there is a concerted theoretical effort to connect such a signal to
the physics of inflation Spergel and Zaldarriaga (1997); Kamionkowski and
Kosowsky (1998); Knox and Song (2002); Baumann et al. (2015); Kamionkowski and
Kovetz (2016). Yet, despite these successes and prospects, there is no leading
theoretical model of inflation, and investigations into the fundamental
physics origin of the inflationary era is an ongoing concern. Most efforts
focus on building a sufficiently flat potential for slow-roll scalar field
evolution. However, Planck-scale quantum corrections make it a challenge to
obtain viable inflation with an appreciable level of GWs Lyth (1997).

A new approach, consisting of a method to achieve slow roll inflation in a
steep potential, has recently garnered attention. In this scenario, a coupling
between an axion-like inflaton and the vacuum expectation value (vev) of a
vector gauge field, through the Chern-Simons term, acts like a brake on the
inflaton evolution Anber and Sorbo (2010); Adshead and Wyman (2012). These
scenarios can produce sufficient inflation without resorting to fine tuning or
Planck-scale physics. Variations on the basic model have been explored
Maleknejad and Sheikh-Jabbari (2013); Maleknejad et al. (2013, 2018); Adshead
et al. (2013a, 2016a); Adshead and Sfakianakis (2017); Dimastrogiovanni et al.
(2018); Domcke et al. (2019a, b), including the impact on reheating Adshead et
al. (2015, 2017, 2018); Adshead et al. (2020a, b). Models such as chromo-
natural inflation Adshead and Wyman (2012) organize the gauge field under
SU(2) to maintain isotropy and spatial homogeneity, though other group
structures are feasible. But all such models share a common feature: the
tensor shear due to the gauge field vev strongly affects GW amplification and
evolution during inflation Anber and Sorbo (2012); Adshead et al. (2013b);
Namba et al. (2013); Maleknejad (2016a, b); Dimastrogiovanni et al. (2017);
Fujita et al. (2018); Caldwell and Devulder (2018); Agrawal et al. (2018a, b);
Fujita et al. (2019); Watanabe and Komatsu (2020). The upshot is that these
models overwrite the standard inflationary expressions that connect the
primordial tensor amplitude to the scale of inflation or the tensor tilt to
the slow roll parameters. A significant piece of this story is missing,
however; the gauge field vev may be expected to survive inflation and
reheating Adshead et al. (2015, 2017); Adshead et al. (2020b), and further
transform the GW spectrum.

There are various scenarios in which relic vector fields emerge from inflation
Davis et al. (2001); Dimopoulos et al. (2002); Demozzi et al. (2009); Adshead
et al. (2016b); Kandus et al. (2011). In the simplest cases, a coupling
between the inflaton or spectator scalar and a U(1) field pumps energy into
the vector field, producing a spectrum of long-wavelength vector
perturbations. The cumulative effect of superhorizon modes contributes an
effective homogeneous background, consisting of both electric and magnetic
components Anber and Sorbo (2010); Domcke et al. (2019a).

In this paper, we investigate the effect of the relic vector fields on the
primordial GW spectrum during the post-inflationary era. For simplicity, we
consider a toy model consisting of a triplet of U(1) fields with vevs that are
arranged in a spatially homogeneous and isotropic configuration. Such a
configuration may occur as a central feature of inflation with Abelian fields
Anber and Sorbo (2010) but may also arise from non-Abelian gauge field
inflation, in the weakly-coupled regime Domcke et al. (2019a), or as a
consequence of the post-inflationary thermalization process Davis et al.
(2001); Dimopoulos et al. (2002); Adshead et al. (2015). The GW - gauge field
system is parameterized by the amplitude of the electric- and magnetic-like
field strengths, with constant fractional energy density during the radiation
era. We show that the spectrum amplitude, tilt, shape, and net chirality are
transformed. This is a new twist on cosmic archaeology Cui et al. (2018): the
imprinting of a non-standard expansion history Allahverdi et al. (2021) or
particle physics history Seto and Yokoyama (2003); Boyle and Steinhardt
(2008); Watanabe and Komatsu (2006); Nakayama et al. (2008); Kuroyanagi et al.
(2011); Saikawa and Shirai (2018); Caldwell et al. (2019); Figueroa and Tanin
(2019) on a GW spectrum. Through examples, we illustrate the consequences for
CMB probes of primordial B-modes and efforts to directly detect a stochastic
GW background. In the context of these scenarios, the strong amplification or
suppression of a primordial GW spectrum has significant observational
consequences.

The paper is organized as follows. In Sec. II we introduce the model and
conventions. We present our results in Sec. III and interpretations in Sec.
IV. Calculation details are presented in the appendices.

##  II Model and Conventions

The vector field sector of this model is composed of a triplet of classical
U(1) fields, $A_{\mu}^{(i)}$, with each copy indexed by $i=1-3$. Hence, after
the relaxation and decay of the inflaton through the epoch of reheating, all
remaining particles and fields are given by the Lagrangian density

$\mathcal{L}=\frac{1}{2}M_{P}^{2}R-\frac{1}{4}F_{\mu\nu}^{(i)}F_{(i)}^{\mu\nu}+\mathcal{L}_{rm}$
(1)

where $\mathcal{L}_{rm}$ is the Lagrangian for matter and radiation. Here and
throughout, the reduced Planck mass is $M_{P}=1/\sqrt{8\pi G}$. We assume the
fermions associated with the vector fields are not present, having been
diluted by inflation or simply too heavy to reach equilibrium with the thermal
radiation. As usual, the field strength tensor for the $i^{\rm{th}}$ field is
$F_{\mu\nu}^{(i)}=\partial_{\mu}A_{\nu}^{(i)}-\partial_{\nu}A_{\mu}^{(i)}$ and
the associated stress-energy tensor is
$T_{\mu\nu}^{(i)}=F^{(i)\alpha}_{\mu}F_{\nu\alpha}^{(i)}-g_{\mu\nu}F^{(i)2}/4$.
The total stress energy of the system is the sum over the three,
$T_{\mu\nu}=\sum_{i}T_{\mu\nu}^{(i)}$. To proceed, both the metric and the
gauge fields are split into a background piece and a linear perturbation:
$g_{\mu\nu}=\eta_{\mu\nu}+h_{\mu\nu}$ and
$A_{\mu}^{(i)}=\bar{A}_{\mu}^{(i)}+\delta A_{\mu}^{(i)}$.

###  II.1 Background

We consider a homogeneous, isotropic spacetime with line element
$ds^{2}=a^{2}(\tau)(-d\tau^{2}+d\vec{x}^{2})$. The conformal time is $\tau$
and the expansion scale factor $a$ is normalized so that $a_{0}=1$ today. The
background metric evolves according to the Friedmann equation,

$H^{2}\equiv\left(\frac{a^{\prime}}{a^{2}}\right)^{2}=H_{0}^{2}\biggl{[}\Omega_{M}\left(\frac{a_{0}}{a}\right)^{3}+\Omega_{R}\left(\frac{a_{0}}{a}\right)^{4}+\Omega_{\Lambda}\biggr{]}$
(2)

where $\Omega_{M}+\Omega_{R}+\Omega_{\Lambda}=1$, and $\prime=d/d\tau$
indicates differentiation with respect to conformal time.

In accordance with the symmetries of the RW spacetime, the vector fields must
preserve homogeneous and isotropic stress-energy at the background level; the
background fields are set up in a “flavor-space” locked configuration,

$F^{(i)}_{j0}=E_{0}\delta^{i}_{j},\qquad F^{(i)}_{kj}=B_{0}\epsilon^{ijk}$ (3)

with constants $E_{0},B_{0}$ representing the vevs of the fields. We stress
that these are not the electromagnetic fields of the Standard Model. For each
flavor of the U(1) there is an “electric” field and co-parallel “magnetic”
field; each flavor is aligned with one of the principal spatial directions. We
note that the two triads must be so aligned, or anti-aligned, in order for the
stress-energy to be isotropic. Put another way, if the electric and magnetic
fields of a given flavor were not aligned, there would be a non-zero Poynting
vector that would break the rotational symmetry. Later, we also consider the
case that the electric- and magnetic-like fields are associated with different
flavors, requiring two triplets of U(1) fields.

With this ansatz the energy density in the background fields is computed via
$\rho_{\bar{A}}=u^{\mu}u^{\nu}T_{\mu\nu}=(3/2)(B_{0}^{2}+E_{0}^{2})/a^{4}$
where $u^{\mu}$ is the four velocity of an observer at rest with respect to
the cosmic frame. The pressure in the $i=x,\,y,\,$ or $z$ direction is
$p_{\bar{A},(i)}=e^{\mu}_{(i)}e^{\nu}_{(i)}T_{\mu\nu}=(1/2)(B_{0}^{2}+E_{0}^{2})/a^{4}$,
the same in all directions, where $e^{\nu}_{(i)}$ are a set of mutually
orthogonal basis vectors. The equation of state is $w=1/3$, but these
stationary fields are not thermal radiation. To parameterize the ratio of the
U(1) energy density to Standard Model radiation (SMR) energy density, we
introduce $R_{B},\,R_{E}$,

$R_{B}\equiv\frac{1}{2}\frac{B_{0}^{2}}{\Omega_{SMR}M_{P}^{2}H_{0}^{2}},\quad
R_{E}\equiv\frac{1}{2}\frac{E_{0}^{2}}{\Omega_{SMR}M_{P}^{2}H_{0}^{2}},$ (4)

where $\Omega_{R}=\Omega_{SMR}(1+R_{B}+R_{E})$.

###  II.2 Perturbations

Our focus in this scenario is the interaction of GWs with the vector fields,
so we begin by considering transverse, traceless, synchronous metric
perturbations, $g_{\mu\nu}=\eta_{\mu\nu}+h_{\mu\nu}$. It will be sufficient to
consider a GW propagating in the $+\hat{z}$ direction,

$\displaystyle h_{\mu\nu}=a^{2}(\tau)\begin{pmatrix}0&0&0&0\\\
0&h_{+}(\tau,z)&h_{\times}(\tau,z)&0\\\
0&h_{\times}(\tau,z)&-h_{+}(\tau,z)&0\\\ 0&0&0&0\end{pmatrix},$ (5)

where we allow both $+$ and $\times$ polarizations.

Similarly, we consider linear U(1) perturbations, $\delta A_{\mu}^{(i)}$,
propagating in the $+\hat{z}$ direction. These wavelike perturbations can be
expressed in terms of two independent polarizations, whereby

$\displaystyle\delta A_{\mu}^{(1)}$
$\displaystyle=\frac{1}{2}M_{P}\left(0,w_{+}(\tau,z),w_{\times}(\tau,z),0\right)$
(6) $\displaystyle\delta A_{\mu}^{(2)}$
$\displaystyle=\frac{1}{2}M_{P}\left(0,w_{\times}(\tau,z),-w_{+}(\tau,z),0\right)$
(7) $\displaystyle\delta A_{\mu}^{(3)}$ $\displaystyle=\left(0,0,0,0\right).$
(8)

The factor of $M_{P}$ makes $w_{+,\times}$ dimensionless, and the factor of
$1/2$ is introduced to make the equations of motion more symmetrical.

The equations of motion are obtained from the perturbed field equations,
$G_{\mu\nu}=M_{P}^{-2}T_{\mu\nu}$ for GWs, and
$\nabla_{\mu}F_{(i)}^{\mu\nu}=0$ for the vector field excitations. We find
that the $+$ and $\times$ polarizations mix, but by switching to a circularly
polarized basis,

$\displaystyle h_{R,L}$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{2}}(h_{+}\pm ih_{\times})$ (9) $\displaystyle
w_{R,L}$ $\displaystyle=$ $\displaystyle\frac{1}{\sqrt{2}}(w_{+}\pm
iw_{\times}),$ (10)

the wave equations for the two polarizations separate. For further
convenience, we define $u\equiv ah$ and move to Fourier space. The coupled
equations of motion are

$\displaystyle
u_{\sigma}^{\prime\prime}+\left(k^{2}-\frac{a^{\prime\prime}}{a}+2\frac{B_{0}^{2}-E_{0}^{2}}{a^{2}M_{P}^{2}}\right)u_{\sigma}$
(11)
$\displaystyle\qquad=\frac{2}{aM_{P}}\left[E_{0}w_{\sigma}^{\prime}-\sigma
kB_{0}w_{\sigma}\right],$ (12) $\displaystyle
w_{\sigma}^{\prime\prime}+k^{2}w_{\sigma}$ (13)
$\displaystyle\qquad=\frac{2}{aM_{P}}\left[E_{0}\left(\frac{a^{\prime}}{a}u_{\sigma}-u_{\sigma}^{\prime}\right)-\sigma
kB_{0}u_{\sigma}\right],$ (14)

where $\sigma=\pm 1$ denotes the right (left) circular polarization. We ignore
the effects of photon and neutrino free-streaming and details of the thermal
history of the radiation fluid on the GW background Watanabe and Komatsu
(2006).

The above system of equations displays significant new features affecting the
propagation of GWs in the cosmological spacetime. To begin, there now appears
an effective mass term in the $u$ equation. This term plays an important role
during the radiation-dominated era, when $a^{\prime\prime}/a$ is negligible,
leading to suppression (enhancement) of superhorizon modes when $B_{0}^{2}$ is
greater (less) than $E_{0}^{2}$. Hence, a blue (red) tilt will be induced in a
scale-invariant spectrum. Next, the GWs couple to the wave excitations of the
U(1) fields. Interconversion between the two, $u$ and $w$, will deplete a
spectrum of GWs across a range of frequencies. This will imprint periodic
features onto an otherwise power-law spectrum. Last, the equations of motion
differ in chirality. In the presence of both $E_{0}$ and $B_{0}$ fields,
right-circular polarized waves propagate differently than left-circular.
Consequently, the processed spectrum will exhibit an excess handedness as a
function of frequency. Together, these three features will transform a
primordial spectrum and leave interesting targets for observation and
experiment.

Finally, we note that scalar perturbations of the flavor-space locked
configuration may be equivalently described as perturbations of an $E$\- and
$B$-fluid, both with $w=\delta p/\delta\rho=1/3$ and a scalar shear $\sigma$.
Details are given in Appendix A. The equations of motion are identical to the
color electrodynamics case previously studied Bielefeld and Caldwell (2015).
There, it was shown that the scalar perturbations are stable and leave
negligible imprint on the CMB, given the bounds on the relic energy density.
Hence, we will not consider the scalar perturbations further in our analysis.

###  II.3 Stochastic Gravitational Wave Background

The observable of interest that displays these new features is the GW spectrum
$\Omega_{GW}$. This quantity is defined as the energy density in GWs,
$\rho_{GW}$, per logarithmic frequency interval, $\ln k$, normalized by the
critical density, $\rho_{c}$ Watanabe and Komatsu (2006),

$\displaystyle\Omega_{GW}(k,\tau)$
$\displaystyle\equiv\frac{1}{\rho_{c}}\frac{d\rho_{GW}}{d\ln k}$ (15)
$\displaystyle=\frac{1}{12(Ha)^{2}}\sum_{\sigma}\Delta^{2}_{h_{\sigma},\rm{prim}}\left|\mathcal{T}_{h_{\sigma},k}^{\prime}(\tau)\right|^{2}.$
(16)

In the second equality we introduce the GW transfer function $\mathcal{T}$,

$\mathcal{T}_{h_{\sigma},k}(\tau)\equiv
h_{k,\sigma}(\tau)/h_{k,\sigma}^{\rm{prim}},$ (17)

which connects the initial GW amplitude to the later, processed amplitude. The
primordial GW power spectrum is

$\Delta_{h_{\sigma},\mathrm{prim}}^{2}(k)\equiv
2\frac{k^{3}}{2\pi^{2}}\left|h_{\sigma,k,\mathrm{prim}}\right|^{2}$ (18)

for each polarization. In the standard case of slow-roll inflation,
$|h_{\sigma,k,{\rm prim}}|^{2}=H_{I}^{2}/k^{3}M_{P}^{2}$ where $H_{I}$ is the
Hubble expansion rate during inflation at horizon exit of a mode with comoving
wavenumber $k$.

The spectrum of U(1) excitations can be described similarly,

$\displaystyle\Omega_{\delta A}(k,\tau)$
$\displaystyle\equiv\frac{1}{\rho_{c}}\frac{d\rho_{\delta A}}{d\ln k}$ (19)
$\displaystyle=\frac{k^{3}}{12\pi^{2}(Ha^{2})^{2}}\sum_{\sigma}\left|w_{\sigma,k}^{\prime}(\tau)\right|^{2}$
(20)
$\displaystyle=\frac{1}{12(Ha^{2})^{2}}\sum_{\sigma}\Delta^{2}_{h_{\sigma},\rm{prim}}\left|\mathcal{T}_{w_{\sigma},k}^{\prime}(\tau)\right|^{2}.$
(21)

Here, we cast the U(1) spectrum in terms of the primordial GW power spectrum,
and define the $w$ transfer function $\mathcal{T}_{w_{\sigma},k}(\tau)\equiv
w_{k,\sigma}(\tau)/h_{k,\sigma}^{\rm{prim}}$, also normalized by the
primordial GW amplitude.

We note that the calculation of the energy densities invokes spatial averaging
over lengths much larger than the wavelength. Although this procedure makes
sense for waves that are well within the horizon, we nevertheless apply the
above expressions to wavelengths approaching the Hubble radius. Furthermore,
at high frequencies we apply a time-averaging filter to obtain the envelope of
the spectrum, $\Omega_{GW,\mathrm{env}}\equiv
2\langle\Omega_{GW}\rangle_{\tau}$, where the time average
$\langle\dots\rangle_{\tau}$ is to be taken over a period of oscillation,
$\tau=2\pi/k$. In subsequent figures, we plot the full oscillatory spectrum
$\Omega_{GW}$ up to $k\simeq k_{eq}$. At higher wavenumbers we plot only the
envelope.

###  II.4 Initial Conditions and Parameter Choices

We prepare the background spacetime in a radiation-dominated universe
following reheating. We use cosmological parameters
$(H_{0},\,\Omega_{M},\,\Omega_{SMR})=\left(67.66\,\mathrm{km}\,\mathrm{s}^{-1}\,\mathrm{Mpc}^{-1},\,0.3111,\,9.138\times
10^{-5}\right)$ obtained from the Planck 2018 TT, TE, EE, lowE, lensing, BAO
best fit values Aghanim et al. (2020). The energy density of the U(1) fields
is set by choosing values of $R_{B},\,R_{E}$. This additional, relativistic
energy density is constrained by Big Bang Nucleosynthesis (BBN) to add no more
than $\Delta N_{\rm eff}<0.43$ (2$\sigma$) extra massless neutrinos Fields et
al. (2020). CMB constraints on extra species, which require some assumptions
about the clustering properties, yield $\Delta N_{\rm eff}<0.3$ (95% CL)
Aghanim et al. (2020). A joint BBN-CMB analysis gives $\Delta N_{\rm
eff}<0.168$ (95% CL) Fields et al. (2020). We will use this tighter bound to
set a limit on the energy density. Hence, expressing $R_{B},\,R_{E}$ in terms
of the effective number of neutrino species,

$R_{B}+R_{E}=\frac{\frac{7}{8}\left(\frac{4}{11}\right)^{4/3}\Delta N_{\rm
eff}}{1+\frac{7}{8}\left(\frac{4}{11}\right)^{4/3}N_{\rm eff}},$ (22)

where $N_{\rm eff}=3.046$, we obtain the upper limit
$R_{B}+R_{E}<R_{max}=0.0225$ (95% CL).

We assume the existence of a spectrum of primordial GWs generated by slow roll
inflation. We set the spectrum amplitude at the CMB pivot scale,
$k_{*}=0.05\,\mathrm{Mpc}^{-1}$, to be consistent with bounds on the tensor-
to-scalar ratio $r_{*}$ and the normalization of the scalar spectrum,
$A_{S}=2.101^{+0.031}_{-0.034}\times 10^{-9}$ Aghanim et al. (2020). A Planck
2018 analysis yields an upper bound $r_{*}<0.06$ (95% CL) Akrami et al.
(2020), whereas a subsequent study places the upper limit at $r_{*}<0.044$
(95% CL) Tristram et al. (2020). To be concrete, we will use $r_{*}=0.04$ to
set the spectrum amplitude, $\Delta^{2}_{h_{\sigma},{\rm
prim}}(k_{*})=H_{I}^{2}(k_{*})/M_{P}^{2}\pi^{2}$, whereby summing over
polarizations we obtain $H_{I}(k_{*})\approx 2.0\times 10^{-5}M_{P}$. The
amplitude at other wavelengths, ranging from the present-day Hubble radius
down to the Hubble radius at the end of reheating, depends on the details of
both the inflationary and reheating scenarios. For simplicity, in this work,
we make the gross assumptions that the primordial spectral tilt is vanishingly
small, $n_{T}\ll 1$, and that reheating is instantaneous. Hence, we initialize
the system at the end of reheating, at conformal time
$\tau_{RH}=1/(a_{RH}H_{RH})$, where
$a_{RH}=(\Omega_{R}H_{0}^{2}/H_{RH}^{2})^{1/4}$ and $H_{RH}=H_{I}$. For
simplicity, we ignore the change in number of degrees of freedom of the
radiation fluid. For the perturbations $u_{\sigma,k}$ and $w_{\sigma,k}$, we
focus our attention on modes that are outside of the horizon at the end of
reheating. This sets the maximum frequency that we consider
$k_{RH}=a_{RH}H_{RH}$. We initialize the perturbations so that the GW is
frozen outside the horizon, as it would be in the absence of the relic vector
fields, while the vector field perturbation is frozen at zero. The transfer
functions cast in terms of the variables $h$, or $u$, and $w$ have initial
conditions

$\displaystyle(\mathcal{T}_{h_{\sigma},k},\mathcal{T}_{h_{\sigma},k}^{\prime},\mathcal{T}_{w_{\sigma},k},\mathcal{T}_{w_{\sigma},k}^{\prime})_{RH}=$
$\displaystyle\,(1,0,0,0)_{RH},$ (23)
$\displaystyle(\mathcal{T}_{u_{\sigma},k},\mathcal{T}_{u_{\sigma},k}^{\prime},\mathcal{T}_{w_{\sigma},k},\mathcal{T}_{w_{\sigma},k}^{\prime})_{RH}=$
$\displaystyle\,(1,aH,0,0)_{RH}.$ (24)

We are now prepared to evolve the system of equations.

###  II.5 Numerical Evolution in the High Frequency Regime

To compute the GW spectral energy density across a broad range of frequencies
we must evolve Eqs. (2), (12) and (14) for the transfer functions of $u,\,w$.
For each Fourier mode $k$, we evolve the system of equations from the end of
reheating, $\tau_{RH}$, to today, $\tau_{0}$. We sample approximately 27
decades in wavenumber, $k\in[10^{-5},\,10^{22}]~{}\mathrm{Mpc}^{-1}$, sampling
linearly for $k\lesssim k_{eq}$ and logarithmically for greater values. Each
point in the $(k,\Omega_{GW}(k,\tau_{0}))$ plane thus corresponds to a single
mode evolution of the equations of motion. Accurately tracking the oscillatory
behavior of every mode once it enters the horizon presents a computational
challenge, so some comments on the procedure are in order.

Solving the system of equations exactly must be done numerically. However, the
numerical integration becomes computationally unfeasible once the modes are
deep within the horizon, $k\gg\mathcal{H}$, due to the small step size
required to resolve the high frequency oscillations of the GW and U(1) modes.
To circumvent this issue, we evolve each mode numerically from $\tau_{RH}$
until it is deep within the horizon at some time $\tau_{m}\sim 10^{4}/k$. At
this point we match the numerical solution to an approximate analytic
solution. In this high frequency limit, $k\tau\gg 1$, a WKB approximation is
made, setting $u(\tau)=\mathrm{Re}\left[U(\tau)e^{ik(\tau-\tau_{m})}\right]$
and assuming several conditions: (1) the mode is deep within the horizon,
$k\gg\mathcal{H}$; (2) the envelope function $U$ varies slowly compared to
$k$, $\left|U^{\prime}\right|\ll k\left|U\right|$; (3) the
$a^{\prime\prime}/a$ term is negligible (e.g. in a radiation-dominated
background). This analytic solution captures the remaining evolution from
$\tau_{m}$ to $\tau_{0}$, whereby the envelope function $U(\tau)$ is

$U(\tau)=U_{+}e^{i\omega\int_{\tau_{m}}^{\tau}d\tau^{\prime}/a(\tau^{\prime})}+U_{-}e^{-i\omega\int_{\tau_{m}}^{\tau}d\tau^{\prime}/a(\tau^{\prime})}$
(25)

and $\omega=\sqrt{({B}_{0}^{2}+{E}_{0}^{2})/M_{P}^{2}}$ (with an analogous
solution for $w(\tau)$). The full details of this high-frequency analytic
solution and matching technique are presented in Appendix B.

Here we briefly consider the standard case, without the U(1) fields, to
illustrate this procedure. The equation of motion for the GW is

$\displaystyle
u^{\prime\prime}_{\mathrm{std}}+\left(k^{2}-\frac{a^{\prime\prime}}{a}\right)u_{\mathrm{std}}=0.$
(26)

The numerical solution requires prohibitively small step sizes once the mode
is deep within the horizon, but this can be easily circumvented with the
analytic approximations described above. Once a mode is deep inside the
horizon, we use the high frequency ansatz for $u(\tau)$ to find a matching
solution to the exact evolution, valid for $\tau>\tau_{m}$. In this standard
case, the envelope function $U$ is a constant. The resulting spectral energy
density, shown in Fig. 1, is easily characterized. For superhorizon modes
$k<a_{0}H_{0}$, $\Omega_{GW}\propto k^{4}$; for modes that entered the horizon
since radiation-matter equality, $a_{0}H_{0}<k<k_{eq}$, $\Omega_{GW,{\rm
env}}\propto k^{-2}$; and for modes that entered the horizon prior,
$k>k_{eq}$, the envelope is a constant $\Omega_{GW,{\rm
env}}\simeq\Omega_{R}H_{I}^{2}/12\pi^{2}M_{P}^{2}$.

In more general cases, the envelope functions for $u,\,w$ satisfy coupled
first-order differential equations; the values of $u,\,w$ at the matching
point are used to provide boundary conditions for the envelope solution. The
analytic solution then captures the evolution from $\tau_{m}$ to $\tau_{0}$.

Figure 1: A comparison of $\Omega_{GW}$ between the standard case (dark blue,
short dashed) and this model in a) a B-dominant case, upper panel,
$(R_{B},R_{E})=(10^{-2},10^{-3})$; and b) an E-dominant case, lower panel,
$(R_{B},R_{E})=(10^{-3},10^{-2})$. Both polarizations in this model are shown:
$\sigma=1$ (light red, long dashed) and $\sigma=-1$ (light blue, solid).

##  III Results

We compute the present-day GW spectrum $\Omega_{GW}$ in a variety of cases, to
demonstrate the impact of the U(1) gauge fields. We begin with a B-dominant
case, $(R_{B},R_{E})=(10^{-2},10^{-3})$, and an E-dominant case,
$(R_{B},R_{E})=(10^{-3},10^{-2})$, as illustrated in the upper and lower
panels respectively of Fig. 1. Both circular polarizations are shown. All
curves are generated assuming a standard inflationary scenario with the same
inflationary scale $H_{I}\approx 2.0\times 10^{-5}M_{P}$ and a flat primordial
GW spectrum, $n_{T}=0$. Upon inspection, it is seen that the GW spectrum
acquires three new features compared to the standard case: a tilt, a net
circular polarization, and oscillations with logarithmic frequency. Let us now
consider each feature in turn.

Figure 2: The behavior of the tilt functions $P_{B}$ (dark red, dot-dashed)
and $P_{E}$ (green, solid) at $k=k_{eq}=k_{piv}$ across a range of background
values $R_{X}$. In the B-dominant case the suppression factor reaches a
minimal value of $P_{B}\approx 10^{-5}$ for $R_{B}=R_{max}$. The E-dominant
case reaches a maximal value of $P_{E}\approx 4$ for $R_{E}=R_{max}$.

###  III.1 Tilt

The tilt of the spectrum is a consequence of the effective mass term
contributed by the U(1) background fields, $E_{0}$ and $B_{0}$. The mass term,
seen in Eq. (12), modifies the dispersion relation during the radiation era,
so that low frequency, long wavelength modes are no longer frozen outside the
horizon. Rather, superhorizon modes are amplified for $E_{0}>B_{0}$, and
suppressed for $B_{0}>E_{0}$. Because modes with lower frequencies enter the
horizon later, this modified superhorizon evolution has a larger cumulative
effect for longer wavelengths. Hence, $E_{0}>B_{0}$ will impart a red tilt to
modes $k\gtrsim k_{eq}$, and a blue tilt for $B_{0}>E_{0}$. The magnitude of
the enhancement or suppression $P_{X}$ for $X=B,\,E$ at wavenumber $k_{eq}$ is
illustrated in Fig. 2. The maximum enhancement is a factor of $4$ in the
E-dominant case, whereas the B-dominant suppression can push the spectrum down
by orders of magnitude. Because this phenomenon produces a constant offset for
$k<k_{piv}\approx k_{eq}$, the enhancement or suppression affects wavelengths
relevant for the CMB.

Figure 3: Demonstration of the approximation of the modified GW spectrum using
the tilt function, Eqs. (27-28), in a scenario with
$(R_{B},R_{E})=(R_{max},0)$. The approximation (orange, long dashed) of the
spectrum as the product of the standard spectrum, $\Omega_{GW,\rm{std}}$ (dark
blue, short dashed), and a tilt function, $P_{B}$, reliably predicts the tilt
of the numerical result (light blue, solid).

The degree of tilt can be predicted by solving the equations of motion in the
limit $k\rightarrow 0$ and taking the ratio of the transfer function amplitude
at horizon entry $k\tau_{entry}=1$ in the modified and standard cases,
$P=\left|\mathcal{T}_{u}(\tau_{entry})/\mathcal{T}_{u,\rm{std}}(\tau_{entry})\right|^{2}$.
This ratio $P$ describes the frequency-dependent modification of the standard
spectrum produced by the gauge fields. Here we quote the result for the simple
case where one of the background fields is zero, since the tilt is only
significant when one background field dominates over the other. In these
cases, momentarily setting aside other effects, the newly tilted spectrum can
be well-approximated as a product of the standard spectrum and a tilt function

$\Omega_{GW,X}\simeq\Omega_{GW,\mathrm{std}}P_{X}(k).$ (27)

The tilt functions in this case are given by piecewise continuous functions

$\displaystyle P_{B}(k)$ $\displaystyle={\cal A}_{R}^{2}\left(\frac{{\rm
max}(k,k_{piv})}{k_{RH}}\right)^{n_{B}}$ (28) $\displaystyle P_{E}(k)$
$\displaystyle=\left[2+{\cal A}_{R}\left(\frac{{\rm
max}(k,k_{piv})}{k_{RH}}\right)^{\frac{1}{2}n_{E}}\right]^{2},$ (29)

where $k_{piv}$ is a pivot scale, $n_{B}=1-\sqrt{1-16R_{B}}$,
$n_{E}=1-\sqrt{1-16R_{E}}$, $k_{RH}=a_{RH}H_{RH}$, and

${\cal
A}_{R}=\frac{R_{B}-R_{E}}{R_{B}+R_{E}}\frac{1+\sqrt{1-16(R_{B}+R_{E})}}{2\sqrt{1-16(R_{B}+R_{E})}}.$
(30)

See Appendix C for details of this calculation. The modification to the
standard case is a power law for $k_{piv}<k<k_{RH}$, while for $k<k_{piv}$ it
becomes a constant offset. Intuitively $k_{piv}$ should be close to $k_{eq}$,
because after radiation-matter equality, the U(1) fields dilute rapidly
compared to the background and no longer affect the superhorizon evolution of
the GWs. Hence all modes that enter the horizon after equality have the same
modification to their amplitude. A comparison of the above model to the exact
spectrum is shown in Fig. 3. We find good agreement between the approximate
tilt functions and the numerical results when using $k_{piv}=k_{eq}$.

###  III.2 Chirality

Starting from a GW spectrum with equal amplitude left- and right-circular
polarizations, the presence of the gauge field vevs in this scenario will
cause a net chirality to develop over time. This effect requires the presence
of both background gauge fields $E_{0}$ and $B_{0}$, and is relevant as the
wave enters the horizon. The necessity of $B_{0}$ is an obvious consequence of
the axial nature of the magnetic field. However, $E_{0}$ is also necessary, if
only to provide a reference against which left and right can be defined. This
can be seen by recognizing that when $E_{0}=0$, the equations of motion for
the two polarizations differ only by a minus sign, $u_{R}=-u_{L}$, which has
no consequence. Furthermore, in both the $k\rightarrow 0$ and $k\gg{\cal H}$
limits, for wavenumbers that are far from the effective Hubble scale
introduced by the gauge fields, the polarization $\sigma$ drops from the
dynamical equations.

Figure 4: The squared GW transfer function $\mathcal{T}_{h}^{2}$ at
recombination with $(R_{B},R_{E})/0.02=(1-10^{-3},10^{-3})$ for both the
$\sigma=-1$ (light blue, solid) and $\sigma=1$ (light red, long dashed)
polarizations, compared to the standard case (dark blue, short dashed).

Interestingly, considering the cases in Fig. 1, the maximum chirality in the
B-dominant case can be quite significant, with the amplitude of the two
polarizations differing by over an order of magnitude at some frequencies,
whereas the chirality is weaker when the dominance is reversed; the chirality
depends on the comparison of $B_{0}$ relative to $E_{0}$. To further
characterize the chirality in this model, in Fig. 4 we compare the GW transfer
function at CMB formation in a B-dominant example to the standard case. The
suppressed GW amplitude and the net circular polarization are both
consequences of the interactions between the U(1) gauge fields and the
primordial GWs. The $\sigma=\pm 1$ polarization enters the horizon later
(earlier) and is relatively suppressed (amplified). In this scenario the
resultant B mode spectrum of the CMB would have a net circular polarization.

###  III.3 Oscillations

The oscillation of the spectral energy density over broad frequency ranges is
best understood by considering the spectra of both GWs and the U(1)
perturbations $\delta A$ together. Fig. 5 makes clear that the oscillations in
$\Omega_{GW}$ are complementary to those in $\Omega_{\delta A}$: as one rises,
the other falls, and the sum yields a flat spectrum. This feature is a
consequence of the interconversion between the GWs and the U(1) perturbations
in the presence of the background $E_{0}$ and $B_{0}$ fields. The phenomenon
is referred to alternately as the Gertsenshteyn effect Gertsenshteyn (1962),
photon-graviton conversion Poznanin (1969); Boccaletti et al. (1970);
Zeldovich (1974), and more generally as GW - gauge field oscillations Caldwell
et al. (2016). We now discuss the main idea underpinning this effect.

The interconversion is seen in the amplitude of the GW and U(1) perturbations
during the radiation era. Subhorizon, the modes oscillate rapidly with
frequency $k$, modulated by the more slowly varying envelope

$e^{i\omega\int^{\tau}_{\tau_{m}}d\tau^{\prime}/a(\tau^{\prime})}.$ (31)

The lower bound of integration is set by horizon entry $\tau\sim k^{-1}$. The
upper bound is set by radiation-matter equality, when the interconversion
process effectively turns off. Hence, modes that enter the horizon earlier
(later) will accumulate more (less) phase modulation. Because the scale factor
evolves as $a\propto\tau$ in the radiation era, the phase acquires a
logarithmic dependence on wavenumber. This slow modulation survives the time-
averaging in the evaluation of the spectral energy density, and imparts the
log-scale modulation seen in Fig. 5. A detailed derivation is given in
Appendix B.

Figure 5: The spectral energy densities of the GW (light red, long dashed) and
U(1) perturbations (dark yellow, dot-dashed) for the right-handed ($\sigma=1$)
polarizations in an $R_{B}=R_{E}=R_{max}/2$ case. For comparison we also show
$\Omega_{GW,\rm{std}}$ (dark blue, short dashed) and the sum of the GW and
U(1) spectral energy densities (light blue, solid); note these two curves are
nearly on top of one another. Figure 6: Demonstration of the approximation of
the oscillatory GW spectrum using the fitting function, Eq. (32), in a
scenario with $R_{B}=R_{E}=R_{max}/2$. Using a pivot scale $k_{piv}=k_{eq}/6$,
the approximation (orange, long dashed) of the spectrum as the product
$\cos^{2}[\Phi]\times\Omega_{GW,\mathrm{std}}$ reliably predicts the shape of
the numerical result (light blue, solid). The standard case (dark blue, short
dashed) is included for comparison. Figure 7: The behavior of the $\sigma=-1$
(light blue, solid) and $\sigma=1$ (light red, long dashed) polarizations of
the GW spectra in increasingly B-dominant cases for $R_{tot}=0.02$.

This analysis allows one to predict the oscillatory shape of the spectra by
relating the background U(1) energy density in the model to the frequency
dependence. The new oscillatory contribution to $\Omega_{GW}$ can be
approximated simplistically by multiplying the $k\gtrsim k_{eq}$ part of the
otherwise power-law spectrum by $\cos^{2}\Phi$ where

$\Phi=\sqrt{2(R_{E}+R_{B})}\ln\left({k}/{k_{piv}}\right),$ (32)

and $k_{piv}$ is a new pivot scale that should also be close to $k_{eq}$ for
the same reasons as in Section III.1. Because of the simplistic nature of this
approximation, the pivot scale $k_{piv}$ varies slightly with $R_{E},R_{B}$.
When $R_{E}\approx R_{B}$, using $k_{piv}/k_{eq}=1/6$ in Eq. (32) provides an
excellent description of the oscillations. In the limiting case where
$R_{B}\rightarrow 0$, we find good results using the empirical relation
$k_{piv}/k_{eq}=a(R_{E}/b)^{c}$ for $(a,b,c)=(1/5,1/50,-9/50)$. When
$R_{E}\rightarrow 0$, $k_{piv}/k_{eq}=a(R_{B}/b)^{c}$ for
$(a,b,c)=(1/18,1/50,1/20)$ gives good results. Setting the phase $\Phi$ to an
even (odd) multiple of $\pi/2$ gives the frequency for a peak (dip). Adjacent
extrema at $k_{1}$ and $k_{2}$ are thus related by
$\ln(k_{1}/k_{2})=(\pi/2)/\sqrt{2(R_{E}+R_{B})}$. Fig. 6 shows a comparison of
this simple prediction with the detailed numerical calculation.

The simple cosine model of Eq. (32) is effective at capturing the locations of
the oscillatory features, but the full structure of the GW spectrum can be far
richer in the presence of chirality. We provide an example in the four panels
of Fig. 7, in which we fix $R_{tot}=R_{B}+R_{E}=0.02$ and progressively lower
the fraction of $R_{E}/R_{tot}$. We observe that the frequency dependence of
the GW spectra is complicated, and the spectra are substantially different for
the two polarizations. The $\sigma=1$ polarization spectrum is generally lower
in amplitude than for $\sigma=-1$; the underlying reason is that for a given
k-mode, the $\sigma=1$ transfer function begins to oscillate sooner as it
approaches horizon entry. Furthermore, the $\sigma=1$ spectrum appears to have
one deep minimum, like a beat frequency, that progressively shifts to lower
frequency minima as $R_{E}/R_{tot}$ is lowered. The behaviour displayed in the
figure clearly shows the limitations of the $\cos^{2}(\Phi)$ model when
significant chirality is present. Finally, it is only once $R_{E}/R_{tot}$ is
lowered below $10^{-4}$ that the two polarization spectra begin to converge.
Roughly speaking, for larger (smaller) $R_{tot}$, the trend towards
convergence occurs at smaller (larger) values of $R_{E}/R_{tot}$.

###  III.4 Sextet Model

We extend our analysis by considering a scenario in which the electric- and
magnetic-like fields originate from distinct gauge groups. Specifically we
consider a model variation in which the electric and magnetic backgrounds are
distributed among two separate U(1) triplets. This model, which we refer to as
a sextet, is very similar to the triplet model, so we keep our discussion
brief. Calculation details are presented in Appendix D.

The sextet and triplet models have the same equations of motion in the
$k\rightarrow 0$ limit, so the superhorizon behavior, and therefore the tilt
in the GW spectrum, is unchanged. The models are also the same in the limit
that either background, $E_{0}$ or $B_{0}$, goes to zero. However, the sextet
model exhibits no handedness, so the excess circular polarization in the
transfer function in Fig. 4, and the rich features seen in the upper panels of
Fig. 7, are absent. Rather, the unpolarized GW spectral energy density in the
sextet model generally behaves like an intermediate of the enhanced and
suppressed polarizations seen in the triplet model. We demonstrate this in the
upper panel of Fig. 8.

Figure 8: Upper panel: The GW spectrum for a single polarization of the sextet
model (green, dot-dashed) compared to that of the $\sigma=-1$ (light blue,
solid) and $\sigma=1$ (light red, long dashed) polarizations in the triplet
model, for $R_{E}=R_{B}=R_{max}/2$. The standard case (dark blue, short
dashed) is included for comparison. Lower panel: Comparison of spectra in the
sextet model for a single polarization of GW (green, dot-dashed), $w$
excitation (orange, long dashed), and $y$ excitation (light blue, solid), for
$R_{E}=R_{B}=R_{max}/2$.

Furthermore, the sextet model includes two independent wavelike gauge field
excitations that mix with the GWs, so the structure of the gauge field spectra
can be more complicated. Importantly, the oscillations in the GW spectrum
envelope in the two models are described equally well by the approximation Eq.
(32), excepting strongly chiral cases in the triplet model. As seen in the
lower panel of Fig. 8, the oscillations in the two gauge field excitations are
together complementary to the oscillations in the GW spectrum.

##  IV Discussion

In the present work we have investigated a toy model in which a homogeneous,
isotropic configuration of relic U(1) vector fields interact with primordial
GWs throughout the radiation-dominated portion of the post-inflationary
Universe. For simplicity, and to isolate the effects of the relics after
inflation, we have assumed instantaneous reheating and that the primordial GW
spectrum is generated from a standard inflationary scenario, with $n_{T}=0$.
However, it is straightforward to superimpose the effects of the relics on any
other primordial spectrum.

We have focused our study on the imprint left on the GW spectral energy
density, $\Omega_{GW}$, finding that the relic fields impart to the spectrum
three distinct features: a tilt, a net circular polarization, and oscillations
across the high frequency ($k\gtrsim k_{eq}$) spectrum.

In particular, the tilt, induced by the U(1) background and demonstrated in
Fig. 9, is our primary result. Compared to previous work, e.g. Adshead and
Wyman (2012); Maleknejad and Sheikh-Jabbari (2013); Bielefeld and Caldwell
(2015), which explored a gauge vev that produces a background “electric”
field, the setup presented here contains a new ingredient at the background
level in the form of an independent axial “magnetic” field. This new
ingredient allows the GW effective mass term to be positive, altering the
dispersion relation and producing the tilt described in Section III.1. In the
maximal scenario, the $R_{B}\gg R_{E}$ case permits a much stronger tilt than
in the reverse case, and so the suppression when $R_{B}$ dominates can be
quite significant. The consequences of such a scenario are noteworthy.

Broadly speaking this scenario is one of several in which the additional
ingredients complicate the relationship between the present day tensor-to-
scalar ratio $r$ and the inflationary scale $H_{I}$. The novel element here is
that this effect is a consequence of dynamics after inflation, rather than
during it. We have demonstrated that this tilt can reduce the amplitude of
$\Omega_{GW}$ by up to $\sim 5$ orders of magnitude at CMB frequencies
$k\lesssim k_{eq}$, all while keeping $H_{I}$ fixed. This is important for
inflation model building because it implies that models that predict $r$
values previously considered too large may instead be within observational
bounds due to this post-inflationary effect.

Figure 9: GW spectral density versus wavenumber and frequency are shown,
relative to the power-law sensitivity curves for several proposed or
futuristic GW observatories. The GW spectra correspond to a B-dominant case,
$(R_{B},R_{E})=(R_{max},0)$, with $H_{I}=2\times 10^{-5}M_{P}$ (light red,
long dashed) and $H_{I}=2M_{P}$ (light blue, solid). This scenario represents
the maximal modification to the standard case (dark blue, short dashed) in
which $(R_{E},R_{B})=(0,0)$, given the energy density bounds on the relic
cosmic vector fields. Any power-law spectral density that crosses a
sensitivity curve is detectable at $\mathrm{SNR}=5$ in a
$t_{obs}=4\,\mathrm{yr}$ experiment.

One illustrative example is single field slow-roll inflation with a quadratic
potential. Another example is chromonatural inflation Adshead and Wyman
(2012). Both of these models overproduce GWs for acceptable values of the
scalar spectral index, $n_{s}$. For the free massive field, $r\sim 0.133$.
However, $R_{B}\gtrsim 10^{-1}R_{max}$ would be sufficient to drop the
amplitude at CMB scales down below $r\lesssim 0.04$. For the case of
chromonatural inflation, the tensor-to-scalar ratio is estimated to be
$r\sim{\cal O}(10^{2})$. In this case, a larger suppression factor with
$R_{B}\sim 3R_{max}/4$ could bring the model within the $1\sigma$ contour in
the $n_{s}-r$ plane Ade et al. (2018); Akrami et al. (2020). Fig. 2 and Eq.
(28) show the range of suppression possible in this model. In numbers, the
range $R_{B}\in\left[5\times 10^{-3},10^{-2}\right]$ roughly corresponds to a
suppression factor of $P_{B}\in\left[10^{-1},10^{-2}\right]$, which may be
enough to revitalize a variety of models with otherwise too strong GW
backgrounds. These arguments also apply to models of inflation with axion
spectator fields.

The post-inflationary effect of the vector fields on $\Omega_{GW}$ has
implications for CMB probes of B-modes as well as direct detection by future
GW observatories. Future CMB experiments such as the Simons Observatory Ade et
al. (2019), LiteBIRD Hazumi et al. (2019), and CMB-S4 Abazajian et al. (2016),
with sensitivities approaching $\delta r\sim 10^{-3}$ can constrain the
$R_{E}-R_{B}$ parameter space in multiple ways. A bound on, or detection of,
primordial B-modes may be interpreted as a joint bound on a model of inflation
and the presence of post-inflationary vector fields. Forthcoming CMB
experiments are also expected to improve measurement of $\Delta
N_{\mathrm{eff}}$. The GW chirality catalyzed by the vector fields, however,
may prove too weak to be detectable with the CMB Gluscevic and Kamionkowski
(2010); Gerbino et al. (2016); Thorne et al. (2018).

Future GW detectors may find a new target within these scenarios. For
illustrative purposes only, in Fig. 9 we show the power-law sensitivity curves
Thrane and Romano (2013); Smith and Caldwell (2019); Schmitz (2021) for the
Laser Interferometer Space Antenna (LISA) Amaro-Seoane et al. (2017), Cosmic
Explorer (CE) Abbott et al. (2017), Einstein Telescope (ET) Hild et al.
(2011); Punturo  (assuming a single Michelson interferometer), and the
futuristic Big Bang Observer (BBO) Crowder and Cornish (2005); Harry et al.
(2006). Under optimal conditions, one of these detectors could be sensitive to
the tilt, chirality, and oscillatory features imprinted in the stochastic GW
background.

###### Acknowledgements.

This work is supported in part by U.S. Department of Energy Award No. DE-
SC0010386.

##  Appendix A Scalar Perturbations

In this appendix we compute the scalar perturbations in the U(1) vector field
model in the background flavor-space locked configuration, Eq. (3). For this
analysis both the metric and U(1) fields are perturbed with only scalar
degrees of freedom, as linear perturbations in a scalar-vector-tensor
decomposition ensure the scalar, vector, and tensor modes do not mix. The U(1)
freedom in this model allows us to perform a coordinate rotation such that the
Fourier vector points along the $+z$ axis, so the perturbations will be
functions of $(\tau,z)$ only. The metric in a flat, RW Universe is then split
into a background metric and linear perturbation,
$g_{\mu\nu}=\eta_{\mu\nu}+\delta g_{\mu\nu}$. In a gauge with smooth spatial
sections, the flat-slicing gauge, the background metric and perturbations,

$\displaystyle\eta_{00}$ $\displaystyle=-a^{2}(\tau)$ (33)
$\displaystyle\eta_{ij}$ $\displaystyle=a^{2}(\tau)\delta_{ij}$ (34)
$\displaystyle\delta g_{00}$ $\displaystyle=2a^{2}(\tau)\Phi_{S}(\tau,z)$ (35)
$\displaystyle\delta g_{0i}$ $\displaystyle=a^{2}(\tau)\partial_{z}b(\tau,z),$
(36)

give the corresponding line element $ds^{2}=g_{\mu\nu}dx^{\mu}dx^{\nu}$ as

$ds^{2}=a^{2}(\tau)\left[-(1-2\Phi_{S})d\tau^{2}-2\partial_{z}b\,dzd\tau+\delta_{ij}\,dx^{i}dx^{j}\right].$
(37)

Similarly the U(1) fields are linearly perturbed,
$A^{(i)}_{\mu}=\bar{A}^{(i)}_{\mu}+\delta A^{(i)}_{\mu}$. In analogy with the
SU(2) case Maleknejad and Sheikh-Jabbari (2011) the U(1) triplet perturbations
can be written

$\displaystyle\delta A^{(i)}_{0}$
$\displaystyle=\delta^{(i)k}\partial_{k}\dot{Y}+\delta^{(i)j}u_{j}$ (38)
$\displaystyle\delta A^{(i)}_{j}$
$\displaystyle=\delta^{(i)}_{\phantom{j}j}Q+\delta^{(i)k}\partial_{jk}M+\varepsilon^{(i)k}_{\phantom{j}j}\partial_{k}P+\delta^{(i)k}\partial_{j}v_{k}+\varepsilon^{(i)k}_{\phantom{j}j}w_{k}+\delta^{(i)k}t_{jk}.$
(39)

We consider only scalar perturbations,

$\displaystyle\delta A^{(i)}_{0}$
$\displaystyle=\delta^{(i)k}\partial_{k}\dot{Y}$ (40) $\displaystyle\delta
A^{(i)}_{j}$
$\displaystyle=\delta^{(i)}_{\phantom{j}j}Q+\delta^{(i)k}\partial_{jk}M+\varepsilon^{(i)k}_{\phantom{j}j}\partial_{k}P,$
(41)

where the overdot denotes differentiation with respect to cosmic time $t$.
Imposing the $(\tau,z)$-only dependence of the perturbations and defining
$(W,\delta G,\delta m,\delta C)\equiv(\dot{Y},Q,-\partial_{z}^{2}M,P)$ gives

$\displaystyle\delta A^{(1)}_{\mu}$ $\displaystyle=\left(0,\delta
G(\tau,z),\partial_{z}\delta C(\tau,z),0\right)$ (42) $\displaystyle\delta
A^{(2)}_{\mu}$ $\displaystyle=\left(0,-\partial_{z}\delta C(\tau,z),\delta
G(\tau,z),0\right)$ (43) $\displaystyle\delta A^{(3)}_{\mu}$
$\displaystyle=\left(\partial_{z}W(\tau,z),0,0,\delta G(\tau,z)-\delta
m(\tau,z)\right)$ (44)

which includes all possible independent scalar perturbations to the U(1)
fields. Stress-energy conservation, $\nabla_{\mu}T^{\mu}_{\phantom{a}\nu}=0$,
and the free U(1) equations of motion, $\nabla_{\mu}F^{\mu\nu,(i)}=0$, give
three equations of motion (second order) and one constraint (first order),

$\displaystyle\delta G^{\prime\prime}+k^{2}\delta
G-k^{2}E_{0}b-E_{0}\Phi_{S}^{\prime}$ $\displaystyle=0$ (45)
$\displaystyle\delta C^{\prime\prime}+k^{2}\delta
C+B_{0}\Phi_{S}-B_{0}b^{\prime}$ $\displaystyle=0$ (46) $\displaystyle\delta
m^{\prime\prime}+k^{2}\delta G-k^{2}E_{0}b-k^{2}W^{\prime}$ $\displaystyle=0$
(47) $\displaystyle\delta m^{\prime}-\delta G^{\prime}+E_{0}\Phi_{S}-k^{2}W$
$\displaystyle=0.$ (48)

With these equations we can obtain the fluid variables through the perturbed
stress energy tensor Ma and Bertschinger (1995), giving

$\displaystyle\bar{\rho}\sigma$
$\displaystyle=\frac{k^{2}}{a^{4}}(-E_{0})W+\frac{E_{0}}{a^{4}}\delta
m^{\prime}+\frac{k^{2}}{a^{4}}B_{0}\delta C$ (49)
$\displaystyle(\bar{\rho}+\bar{P})\theta$
$\displaystyle=2\frac{k^{2}}{a^{4}}(-E_{0})\delta
G^{\prime}+2\frac{k^{2}}{a^{4}}B_{0}(\delta C^{\prime}-B_{0}b)$ (50)
$\displaystyle\delta\rho$ $\displaystyle=\frac{-E_{0}}{a^{4}}(3\delta
G^{\prime}-\delta
m^{\prime})-\frac{k^{2}}{a^{4}}E_{0}W+3\frac{E_{0}^{2}}{a^{4}}\Phi_{S}-2\frac{k^{2}}{a^{4}}B_{0}\delta
C$ (51)

where the energy density and pressure at background and linear order satisfy
the equation of state of radiation, $\bar{P}/\bar{\rho}=1/3=\delta
P/\delta\rho$. The evolution of the fluid variables is given by

$\displaystyle\delta^{\prime}$
$\displaystyle=-\frac{4}{3}\left(\theta+k^{2}b\right)$ (52)
$\displaystyle\theta^{\prime}$
$\displaystyle=k^{2}\left(\frac{1}{4}\delta-\Phi_{S}-\sigma\right)$ (53)

where $\delta\equiv\delta\rho/\bar{\rho}$.

Now we transform to the conformal Newtonian (cN) gauge and express the fluid
equations in that gauge. The metric in the cN gauge is

$\displaystyle g_{00}$ $\displaystyle=-a^{2}(1+2\psi_{cN})$ (54)
$\displaystyle g_{0i}$ $\displaystyle=0$ (55) $\displaystyle g_{ij}$
$\displaystyle=a^{2}(1-2\phi_{cN})\delta_{ij}.$ (56)

The cN gauge (with coordinates $\hat{x}^{\mu}$) and our gauge (with
coordinates $x^{\mu}$) are related by a simple coordinate change
$x^{0}=\hat{x}^{0}+b$, with no transformation of the spatial coordinates. This
yields the relationships Ma and Bertschinger (1995) between the metric
perturbations in the two gauges,

$\displaystyle\mathcal{H}b$ $\displaystyle=-\phi_{cN}$ (57)
$\displaystyle-\Phi_{S}+\mathcal{H}b+b^{\prime}$ $\displaystyle=\psi_{cN}.$
(58)

Under coordinate transformation, the stress energy tensor and thus the fluid
variables also transform. It is simplest to proceed by calculating the gauge-
invariant (gi) perturbations, which are equivalent to those in the cN gauge
Mukhanov et al. (1992). With this we can compute the gauge invariant
(equivalently, cN) fluid variables,

${}^{(gi)}\delta$
$\displaystyle=\delta_{cN}=\frac{3E_{0}^{2}\Phi_{S}-k^{2}E_{0}W-E_{0}(3\delta
G^{\prime}-\delta m^{\prime})-2B_{0}k^{2}\delta
C}{3(E_{0}^{2}+B_{0}^{2})/2}+4\phi_{cN}$ (59) ${}^{(gi)}\theta$
$\displaystyle=\theta_{cN}=\frac{-k^{2}E_{0}\delta
G+k^{2}E_{0}b+k^{2}B_{0}\delta C^{\prime}}{E_{0}^{2}+B_{0}^{2}}$ (60)
${}^{(gi)}\sigma$ $\displaystyle=\sigma_{cN}=\frac{-k^{2}E_{0}W+E_{0}\delta
m^{\prime}+k^{2}B_{0}\delta C}{3(E_{0}^{2}+B_{0}^{2})/2},$ (61)

where one can replace and $b$ and $\Phi_{S}$ with cN variables using Eqs. (57)
and (58). These cN fluid variables obey

${}^{(gi)}\delta^{\prime}$
$\displaystyle=\delta_{cN}^{\prime}=-\frac{4}{3}\theta_{cN}+4\phi_{cN}^{\prime}$
(62) ${}^{(gi)}\theta^{\prime}$
$\displaystyle=\theta_{cN}^{\prime}=k^{2}\frac{1}{4}\delta_{cN}-k^{2}\sigma_{cN}+k^{2}\psi_{cN}$
(63) ${}^{(gi)}\sigma^{\prime}$
$\displaystyle=\sigma_{cN}^{\prime}=\frac{2}{3}\theta_{cN}.$ (64)

These are analogous to the first two equations in (63) of Ref. Ma and
Bertschinger (1995) with a modified shear component, in accordance with what
was found in the color electrodynamics case Bielefeld and Caldwell (2015).

##  Appendix B High Frequency GW-U(1) Solutions and Computation of
$\Omega_{GW}$

In this section we develop an approximate analytic solution to the equations
of motion (12-14) in a high frequency limit. This solution is used to compute
$\Omega_{GW}$ given the computational challenges described in Section II.5.
This solution also demonstrates the origins of the GW envelope, Eq. (25), and
the consequent oscillations in $\Omega_{GW}$ discussed in Section III.3. The
equations of motion are

$\displaystyle
u_{\sigma}^{\prime\prime}+\left(k^{2}-\frac{a^{\prime\prime}}{a}+2\frac{(B_{0}^{2}-E_{0}^{2})}{a^{2}M_{P}^{2}}\right)u_{\sigma}=\frac{2}{aM_{P}}\left[E_{0}w_{\sigma}^{\prime}-\sigma
kB_{0}w_{\sigma}\right]$ (65) $\displaystyle
w_{\sigma}^{\prime\prime}+k^{2}w_{\sigma}=\frac{2}{aM_{P}}\left[E_{0}\left(\frac{a^{\prime}}{a}u_{\sigma}-u_{\sigma}^{\prime}\right)-\sigma
kB_{0}u_{\sigma}\right].$ (66)

Let $u_{k}(\tau)=\mathrm{Re}\left[U_{k}(\tau)e^{ik(\tau-\tau_{m})}\right]$ and
likewise for $w$. Here $\tau_{m}\sim 10^{4}/k$ is a constant matching time at
which we match the high frequency solution to the numerical solution. Dropping
subscripts for now, the equations of motion become

$\displaystyle
U^{\prime\prime}+2ikU^{\prime}-\frac{a^{\prime\prime}}{a}U+2\frac{B_{0}^{2}-E_{0}^{2}}{M_{P}^{2}a^{2}}U=\frac{2}{aM_{P}}\left[E_{0}W^{\prime}+ikE_{0}W-\sigma
kB_{0}W\right]$ (67) $\displaystyle
W^{\prime\prime}+2ikW^{\prime}=\frac{2}{aM_{P}}\left[-E_{0}U^{\prime}-ikE_{0}U+E_{0}\mathcal{H}U-\sigma
kB_{0}U\right].$ (68)

We are searching for a high-frequency solution, so assume the following four
conditions: $k\left|U\right|\gg\left|U^{\prime}\right|$, $k\gg\mathcal{H}$,
$k\left|U^{\prime}\right|\gg\left|(a^{\prime\prime}/a)U\right|$, and
$k\left|U^{\prime}\right|\gg\left|U(B_{0}^{2}-E_{0}^{2})/(a^{2}M_{P}^{2})\right|$.
The first condition assumes the envelope functions $U,W$ vary appreciably only
on time scales much larger than $1/k$. The second condition assumes the mode
is deep within the horizon. The third condition is easily satisfied during
radiation domination, during which the gauge radiation terms are most
significant. Using the definition in Eq. (4), the last condition can be
expressed

$\frac{a^{2}}{2\left|R_{B}-R_{E}\right|\Omega_{SMR}H_{0}^{2}}\frac{k\left|U^{\prime}\right|}{\left|U\right|}\gg
1$ (69)

which we will show shortly is valid in the regime of interest. When the above
four conditions hold the equations of motion simplify to

$\displaystyle U^{\prime}$ $\displaystyle=\frac{1}{aM_{P}}\left[E_{0}W+i\sigma
B_{0}W\right]$ (70) $\displaystyle W^{\prime}$
$\displaystyle=-\frac{1}{aM_{P}}\left[E_{0}U-i\sigma B_{0}U\right].$ (71)

Defining a new time variable $dx\equiv d\tau/a$ then yields equations of
motion for simple harmonic oscillators in $x$, giving

$U=U_{+}e^{i\omega x}+U_{-}e^{-i\omega x}$ (72)

where $\omega=\sqrt{(E_{0}^{2}+B_{0}^{2})/M_{P}^{2}}$. We can now return to
the condition in Eq. (69) which becomes

$\frac{ak}{\Omega_{SMR}^{1/2}H_{0}}\frac{\sqrt{2(R_{B}+R_{E})}}{2\left|R_{B}-R_{E}\right|}\frac{|U_{+}e^{i\omega
x}-U_{-}e^{-i\omega x}|}{|U_{+}e^{i\omega x}+U_{-}e^{-i\omega x}|}\gg 1$ (73)

which is satisfied for all relevant modes. This is easily seen by recognizing
that the fraction containing $R_{E},R_{B}$ is always larger than unity, and
the fraction containing $k/H_{0}$ is roughly $k\tau$ during the radiation era,
which is large by assumption $k\gg\mathcal{H}$.

Next we detail how this analytic solution is used to compute the spectral
energy densities. The high frequency solution is matched at time $\tau_{m}$ to
the numerical solution, so it is convenient to redefine the coefficients
$U_{+},U_{-}$ so that

$U=(U_{+c}+iU_{+s})e^{-ik\tau_{m}-i\omega x_{m}}e^{i\omega
x}+(U_{-c}+iU_{-s})e^{-ik\tau_{m}+i\omega x_{m}}e^{-i\omega x},$ (74)

where $(U_{+c},U_{+s},U_{-c},U_{-s})\in\mathbb{R}$. The above choice yields
the full high frequency solution

$\displaystyle
u_{k}(\tau)=U_{+c}\cos(\delta_{+})-U_{+s}\sin(\delta_{+})+U_{-c}\cos(\delta_{-})-U_{-s}\sin(\delta_{-})$
(75)

where we have defined the argument $\delta_{\pm}\equiv
k(\tau-\tau_{m})\pm\omega(x-x_{m})$, so

$\frac{d\delta_{\pm}}{d\tau}=k\pm\frac{\omega}{a(\tau)}\equiv k\pm f(\tau).$
(76)

Note $k\gg\mathcal{H}$ implies $k\gg f$. The solution for $W$ is obtained with
Eq. (70) and inserting the solution for $U$, giving

$\displaystyle
w=W_{c+}\cos(\delta_{+})-W_{s+}\sin(\delta_{+})+W_{c-}\cos(\delta_{-})-W_{s-}\sin(\delta_{-})$
(77)

with the $W_{s,c,\pm}$ coefficients

$\displaystyle W_{c+}$ $\displaystyle\equiv\frac{\sigma
B_{0}U_{+c}-E_{0}U_{+s}}{M_{P}\omega}$ (78) $\displaystyle W_{s+}$
$\displaystyle\equiv\frac{\sigma B_{0}U_{+s}+E_{0}U_{+c}}{M_{P}\omega}$ (79)
$\displaystyle W_{c-}$ $\displaystyle\equiv-\frac{(\sigma
B_{0}U_{-c}-E_{0}U_{-s})}{M_{P}\omega}$ (80) $\displaystyle W_{s-}$
$\displaystyle\equiv-\frac{(\sigma B_{0}U_{-s}+E_{0}U_{-c})}{M_{P}\omega}.$
(81)

The four independent coefficients then are found by solving a system of four
linear equations

$\displaystyle u_{k}(\tau=\tau_{m})$ $\displaystyle\equiv u_{m}=U_{+c}+U_{-c}$
(82) $\displaystyle u_{k}^{\prime}(\tau=\tau_{m})$ $\displaystyle\equiv
du_{m}=-(k+f_{m})U_{+s}-(k-f_{m})U_{-s}$ (83) $\displaystyle
w_{k}(\tau=\tau_{m})$ $\displaystyle\equiv w_{m}=\frac{(\sigma
B_{0}U_{+c}-E_{0}U_{+s})}{M_{P}\omega}-\frac{(\sigma
B_{0}U_{-c}-E_{0}U_{-s})}{M_{P}\omega}$ (84) $\displaystyle
w_{k}^{\prime}(\tau=\tau_{m})$ $\displaystyle\equiv
dw_{m}=-(k+f_{m})\frac{(\sigma
B_{0}U_{+s}+E_{0}U_{+c})}{M_{P}\omega}+(k-f_{m})\frac{(\sigma
B_{0}U_{-s}+E_{0}U_{-c})}{M_{P}\omega}$ (85)

for $f_{m}\equiv f(\tau_{m})$, with solutions

$\displaystyle U_{+c}$ $\displaystyle=\frac{-\bar{E}_{0}k\omega
dw_{m}+\bar{E}_{0}^{2}k(k-f_{m})u_{m}+\sigma\bar{B}_{0}\bar{E}_{0}f_{m}du_{m}+(k^{2}-f_{m}^{2})\sigma\bar{B}_{0}(\omega
w_{m}+\sigma\bar{B}_{0}u_{m})}{2k^{2}\omega^{2}-2f_{m}^{2}\bar{B}_{0}^{2}}$
(86) $\displaystyle U_{+s}$
$\displaystyle=\frac{(f_{m}-k)\left[\bar{E}_{0}k\omega
w_{m}+\sigma\bar{B}_{0}\omega
dw_{m}+\sigma\bar{B}_{0}\bar{E}_{0}f_{m}u_{m}\right]-du_{m}\left[(k-f_{m})\bar{B}_{0}^{2}+k\bar{E}_{0}^{2}\right]}{2k^{2}\omega^{2}-2f_{m}^{2}\bar{B}_{0}^{2}}$
(87) $\displaystyle U_{-c}$ $\displaystyle=\frac{\bar{E}_{0}k\omega
dw_{m}+\bar{E}_{0}^{2}k(k+f_{m})u_{m}-\sigma\bar{B}_{0}\bar{E}_{0}f_{m}du_{m}+(k^{2}-f_{m}^{2})\sigma\bar{B}_{0}\left[\sigma\bar{B}_{0}u_{m}-\omega
w_{m}\right]}{2k^{2}\omega^{2}-2f_{m}^{2}\bar{B}_{0}^{2}}$ (88) $\displaystyle
U_{-s}$
$\displaystyle=\frac{-\left[\bar{E}_{0}^{2}k+(k+f_{m})\bar{B}_{0}^{2}\right]du_{m}+(k+f_{m})\left[\bar{E}_{0}k\omega
w_{m}+\sigma\bar{B}_{0}\omega
dw_{m}+\sigma\bar{B}_{0}\bar{E}_{0}f_{m}u_{m}\right]}{2k^{2}\omega^{2}-2f_{m}^{2}\bar{B}_{0}^{2}},$
(89)

where $(\bar{E}_{0},\bar{B}_{0})\equiv\tfrac{1}{M_{P}}(E_{0},B_{0})$ have been
defined for notational compactness.

For the purposes of computing the spectral energy densities, we need the time
average of $u^{\prime}{}^{2}$ (evaluated at $\tau_{0}$); we have

$u^{\prime}_{k}=-(k+f)\left[U_{+c}\sin(\delta_{+})+U_{+s}\cos(\delta_{+})\right]-(k-f)\left[U_{-c}\sin(\delta_{-})+U_{-s}\cos(\delta_{-})\right].$
(90)

The trigonometric functions in $\delta_{\pm}$ can be expanded into products of
trigonometric functions in $k\tau$ and $\omega x$. For all modes of interest
we have $k\gg f$ and $k\gg\mathcal{H}$, so $\omega x$ and $f(\tau)$ are
effectively constant over the averaging period $2\pi/k$. The time average can
then be safely taken only over the functions in $k\tau$, giving

$\displaystyle\langle u^{\prime}_{k}{}^{2}\rangle_{\tau}=\frac{1}{2}$
$\displaystyle\biggl{\\{}\bigl{[}k^{2}((U_{+c}+U_{-c})^{2}+(U_{+s}+U_{-s})^{2})+2fk(U_{+c}^{2}+U_{+s}^{2}-U_{-c}^{2}-U_{-s}^{2})$
$\displaystyle+f^{2}((U_{+c}-U_{-c})^{2}+(U_{+s}-U_{-s})^{2})\bigr{]}\cos^{2}\left[\omega(x-x_{m})\right]$
$\displaystyle+\bigl{[}k^{2}((U_{+c}-U_{-c})^{2}+(U_{+s}-U_{-s})^{2})+2fk(U_{+c}^{2}+U_{+s}^{2}-U_{-c}^{2}-U_{-s}^{2})$
$\displaystyle+f^{2}((U_{+c}+U_{-c})^{2}+(U_{+s}+U_{-s})^{2}))\bigr{]}\sin^{2}\left[\omega(x-x_{m})\right]$
$\displaystyle+4(k^{2}-f^{2})(U_{+c}U_{-s}-U_{+s}U_{-c})\sin\left[\omega(x-x_{m})\right]\cos\left[\omega(x-x_{m})\right]\biggr{\\}}$
(91)

and an identical expression for $\langle w^{\prime}_{k}{}^{2}\rangle$ with the
replacement
$\left\\{U_{+c},U_{+s},U_{-c},U_{-s}\right\\}\rightarrow\left\\{W_{c+},W_{s+},W_{c-},W_{s-}\right\\}$.
Eq. (91) and the $w$ analogue are what we use in practice to compute the
spectral energy densities $\Omega_{GW}$ and $\Omega_{\delta A}$ in Eqs. (16)
and (21) respectively. The trigonometric functions in $\omega x$ that appear
in Eq. (91) are the cause of the oscillations described in Sec. III.3, and for
a mode that enters in the radiation dominated era, $x-x_{m}$ has a $\ln(k)$
contribution. These provide the justification for the approximate fitting
formula in Eq. (32), although it is clear that the frequency dependence of Eq.
(91) is complicated.

##  Appendix C Modified Superhorizon GW Evolution

In this appendix we solve the equations of motion (12-14) in the low-frequency
limit to examine the effect of the U(1) fields on superhorizon GW evolution.
In the $k=0$ limit the equations of motion become

$\displaystyle
u^{\prime\prime}-\frac{a^{\prime\prime}}{a}u+2\frac{B_{0}^{2}-E_{0}^{2}}{M_{pl}^{2}a^{2}}u=\frac{2E_{0}}{aM_{P}}\left[w^{\prime}\right]$
(92) $\displaystyle
w^{\prime\prime}=\frac{2E_{0}}{aM_{P}}\left[-u^{\prime}+\mathcal{H}u\right].$
(93)

Dividing both equations by $u_{k,\mathrm{prim}}$ recasts the equations of
motion in term of the transfer functions $\mathcal{T}_{u},\mathcal{T}_{w}$. We
again focus on the radiation era during which the background U(1) terms
present here are most relevant, so
$a(\tau)=a_{RH}(\tau/\tau_{RH})=\Omega_{R}^{1/2}H_{0}\tau$ (see Section II.4).
Then the equations of motion can be compactly written,

$\displaystyle\mathcal{T}_{u}^{\prime\prime}+2\frac{\bar{B}_{0}^{2}-\bar{E}_{0}^{2}}{\Omega_{R}H_{0}^{2}\tau^{2}}\mathcal{T}_{u}=\frac{2\bar{E}_{0}}{\Omega_{R}^{1/2}H_{0}\tau}\left[\mathcal{T}_{w}^{\prime}\right]$
(94)
$\displaystyle\mathcal{T}_{w}^{\prime\prime}=-\frac{2\bar{E}_{0}}{\Omega_{R}^{1/2}H_{0}}\frac{d}{d\tau}\left[\mathcal{T}_{u}\tau^{-1}\right],$
(95)

which are easily decoupled and solved by integrating the second equation of
motion and inserting it into the first, keeping in mind the initial conditions
in Eq. (24). This gives the $\mathcal{T}_{u}$ equation of motion

$\tau_{RH}^{2}\mathcal{T}_{u}^{\prime\prime}+2\frac{(\bar{B}_{0}^{2}+\bar{E}_{0}^{2})}{\Omega_{R}H_{0}^{2}(\tau/\tau_{RH})^{2}}\mathcal{T}_{u}-\frac{4\bar{E}_{0}^{2}}{\Omega_{R}H_{0}^{2}(\tau/\tau_{RH})}=0.$
(96)

Dropping the negligible $\mathcal{O}(R_{X}^{2})$ difference between
$(E_{0}^{2},B_{0}^{2})/\Omega_{R}$ and $(E_{0}^{2},B_{0}^{2})/\Omega_{SMR}$,
this yields simple solutions,

$\displaystyle\mathcal{T}_{u}$
$\displaystyle=\mathcal{B}_{R}\left(\frac{\tau}{\tau_{RH}}\right)^{n_{-}}+\mathcal{A}_{R}\left(\frac{\tau}{\tau_{RH}}\right)^{n_{+}}+\frac{2R_{E}}{R_{B}+R_{E}}\frac{\tau}{\tau_{RH}}$
(97) $\displaystyle\mathcal{T}_{w}$
$\displaystyle=-2\sqrt{2R_{E}}\left[\frac{R_{B}-R_{E}}{R_{B}+R_{E}}\left(\frac{\tau}{\tau_{RH}}-1\right)+\frac{\mathcal{B}_{R}}{n_{-}}\left(\frac{\tau}{\tau_{RH}}\right)^{n_{-}}+\frac{\mathcal{A}_{R}}{n_{+}}\left(\frac{\tau}{\tau_{RH}}\right)^{n_{+}}\right]$
(98)

for $n_{\pm}=\frac{1}{2}\Bigl{(}1\pm\sqrt{1-16(R_{B}+R_{E})}\Bigr{)}$ and

$\displaystyle\mathcal{B}_{R}=\frac{R_{E}-R_{B}}{R_{B}+R_{E}}\frac{1-\sqrt{1-16(R_{B}+R_{E})}}{2\sqrt{1-16(R_{B}+R_{E})}}$
$\displaystyle=\frac{R_{E}-R_{B}}{R_{B}+R_{E}}\frac{n_{-}}{\sqrt{1-16(R_{B}+R_{E})}}$
(99)
$\displaystyle\mathcal{A}_{R}=\frac{R_{B}-R_{E}}{R_{B}+R_{E}}\frac{1+\sqrt{1-16(R_{B}+R_{E})}}{2\sqrt{1-16(R_{B}+R_{E})}}$
$\displaystyle=\frac{R_{B}-R_{E}}{R_{B}+R_{E}}\frac{n_{+}}{\sqrt{1-16(R_{B}+R_{E})}},$
(100)

where $|\mathcal{A}_{R}|>|\mathcal{B}_{R}|$. In the standard case without the
gauge fields, the GW transfer function equation of motion (also for
$k\rightarrow 0$ and in the radiation era) is
$\mathcal{T}_{u,\rm{std}}^{\prime\prime}=0$ with solution
$\mathcal{T}_{u,\rm{std}}=\tau/\tau_{RH}$. Therefore in the presence of the
gauge fields the GW amplitude is altered while outside the horizon by a factor

$\frac{\mathcal{T}_{u}(\tau)}{\mathcal{T}_{u,\rm{std}}(\tau)}=\mathcal{B}_{R}\left(\frac{\tau}{\tau_{RH}}\right)^{\frac{1}{2}\left(-1-\sqrt{1-16(R_{B}+R_{E})}\right)}+\mathcal{A}_{R}\left(\frac{\tau}{\tau_{RH}}\right)^{\frac{1}{2}\left(-1+\sqrt{1-16(R_{B}+R_{E})}\right)}+\frac{2R_{E}}{R_{B}+R_{E}}$
(101)

before the mode enters the horizon at $\tau_{entry}=1/k$. Lower frequency
modes enter later and accumulate a larger modification, so this effect is
frequency dependent. Taking a piecewise approximation in which this
$k\rightarrow 0$ description is valid until $\tau=\tau_{entry}$, after which
point the amplitude modulation shuts off and subhorizon oscillation begins,
this modulated amplitude means that $\Omega_{GW}$ is different by a factor

$\displaystyle\left|\mathcal{T}_{u_{k}}(\tau_{entry})/\mathcal{T}_{u_{k},\rm{std}}(\tau_{entry})\right|^{2}=\biggl{[}$
$\displaystyle\mathcal{B}_{R}\left(\frac{k}{k_{RH}}\right)^{\frac{1}{2}\left(1+\sqrt{1-16(R_{B}+R_{E})}\right)}+\mathcal{A}_{R}\left(\frac{k}{k_{RH}}\right)^{\frac{1}{2}\left(1-\sqrt{1-16(R_{B}+R_{E})}\right)}+\frac{2R_{E}}{R_{B}+R_{E}}\biggr{]}^{2}$
(102)

where we have restored the $k$-subscript and used $\tau_{RH}=1/k_{RH}$. The
power laws in $k$ give the GW spectrum the new tilt seen in this model. When
the background fields are absent, the standard case is recovered, and there is
no additional tilt. The tilt also disappears for $R_{E}\approx R_{B}$ so it is
useful to consider the limiting cases when one background field is zero,

$\displaystyle P_{B}$
$\displaystyle=\left|\mathcal{T}_{u_{k}}(\tau_{entry})/\mathcal{T}_{u_{k},\rm{std}}(\tau_{entry})\right|^{2}_{R_{E}=0}=\biggl{[}\mathcal{B}_{R}\left(\frac{k}{k_{RH}}\right)^{\frac{1}{2}\left(1+\sqrt{1-16R_{B}}\right)}+\mathcal{A}_{R}\left(\frac{k}{k_{RH}}\right)^{\frac{1}{2}\left(1-\sqrt{1-16R_{B}}\right)}\biggr{]}^{2}$
(103) $\displaystyle P_{E}$
$\displaystyle=\left|\mathcal{T}_{u_{k}}(\tau_{entry})/\mathcal{T}_{u_{k},\rm{std}}(\tau_{entry})\right|^{2}_{R_{B}=0}=\biggl{[}2+\mathcal{B}_{R}\left(\frac{k}{k_{RH}}\right)^{\frac{1}{2}\left(1+\sqrt{1-16R_{E}}\right)}+\mathcal{A}_{R}\left(\frac{k}{k_{RH}}\right)^{\frac{1}{2}\left(1-\sqrt{1-16R_{E}}\right)}\biggr{]}^{2}.$
(104)

We always have $k<k_{RH}$, so the $\mathcal{A}_{R}$ term is more significant;
dropping the subdominant $\mathcal{B}_{R}$ term then gives

$\displaystyle P_{B}$
$\displaystyle\simeq\biggl{[}\mathcal{A}_{R}\left(\frac{k}{k_{RH}}\right)^{\frac{1}{2}\left(1-\sqrt{1-16R_{B}}\right)}\biggr{]}^{2}$
(105) $\displaystyle P_{E}$
$\displaystyle\simeq\biggl{[}2+\mathcal{A}_{R}\left(\frac{k}{k_{RH}}\right)^{\frac{1}{2}\left(1-\sqrt{1-16R_{E}}\right)}\biggr{]}^{2}.$
(106)

Hence at low frequencies the E-enhancement is limited to roughly a factor of 4
while the B-suppression can be far more significant. Once radiation domination
gives way to matter domination at $\tau_{eq}$, the $E_{0},B_{0}$ background
rapidly dilutes away and the modification to the superhorizon evolution ends.
Thus all modes $k\lesssim k_{eq}$ accumulate the same maximum modification.
This is why the tilt for $k\lesssim k_{eq}$ becomes a constant offset. Hence
we replace $k\rightarrow\mathrm{min}(k,k_{eq})$ in Eqs. (105-106).

##  Appendix D Sextet Model

In this section we extend the model to incorporate six U(1) fields, in which
each field obtains only a single electric or magnetic background. The
corresponding ansatz is

$F^{(i)}_{j0}=E_{0}\delta^{i}_{j},\qquad
F^{(i^{\prime})}_{kj}=B_{0}\epsilon^{(i^{\prime}-3)jk}$ (107)

where charge flavor $i=1,2,3$ is mapped to Cartesian directions $x,y,z$,
whereas charge flavor $i^{\prime}=4,5,6$ is separately mapped to $x,y,z$. That
is, the electric- and magnetic-like fields are associated with different
triplets of U(1) fields. This field configuration gives identical background
behavior to the triplet case, but the linear wavelike perturbations

$\displaystyle\delta A_{\mu}^{(1)}$
$\displaystyle=\frac{M_{P}}{2}\left(0,w_{+}(\tau,z),w_{\times}(\tau,z),0\right)$
(108) $\displaystyle\delta A_{\mu}^{(2)}$
$\displaystyle=\frac{M_{P}}{2}\left(0,w_{\times}(\tau,z),-w_{+}(\tau,z),0\right)$
(109) $\displaystyle\delta A_{\mu}^{(4)}$
$\displaystyle=\frac{M_{P}}{2}\left(0,y_{+}(\tau,z),y_{\times}(\tau,z),0\right)$
(110) $\displaystyle\delta A_{\mu}^{(5)}$
$\displaystyle=\frac{M_{P}}{2}\left(0,y_{\times}(\tau,z),-y_{+}(\tau,z),0\right)$
(111)

change the gravitational wave dynamics. The equations of motion

$\displaystyle
u_{\sigma}^{\prime\prime}+\left(k^{2}-\frac{a^{\prime\prime}}{a}+2\frac{(B_{0}^{2}-E_{0}^{2})}{a^{2}M_{P}^{2}}\right)u_{\sigma}=\frac{2}{aM_{P}}\left[E_{0}w_{\sigma}^{\prime}-\sigma
kB_{0}y_{\sigma}\right]$ (112) $\displaystyle
w_{\sigma}^{\prime\prime}+k^{2}w_{\sigma}=\frac{2}{aM_{P}}\left[E_{0}\left(\frac{a^{\prime}}{a}u_{\sigma}-u_{\sigma}^{\prime}\right)\right]$
(113) $\displaystyle
y_{\sigma}^{\prime\prime}+k^{2}y_{\sigma}=\frac{2}{aM_{P}}\left[-\sigma
kB_{0}u_{\sigma}\right]$ (114)

demonstrate different features than the triplet case. The equations of motion
are invariant under $(\sigma,y)\rightarrow-(\sigma,y)$ or
$(\sigma,u,w)\rightarrow-(\sigma,u,w)$, neither of which modify the spectral
energy densities. Hence the spectral energy densities do not exhibit an excess
handedness. Furthermore, the sextet and triplet models have the same equations
of motion the limit that either $E_{0}$ or $B_{0}$ vanish, and in the
superhorizon $k\rightarrow 0$ limit.

An analogous high frequency solution exists for this system, giving the
following system,

$\displaystyle\frac{dU}{dx}$ $\displaystyle=\bar{E}_{0}W+i\sigma\bar{B}_{0}Y$
(115) $\displaystyle\frac{dW}{dx}$ $\displaystyle=-\bar{E}_{0}U$ (116)
$\displaystyle\frac{dY}{dx}$ $\displaystyle=i\sigma\bar{B}_{0}U.$ (117)

The GW equation of motion can again be cast in terms of a harmonic oscillator
in $x$ with frequency $\omega$, but the gauge field system is changed,

$\displaystyle U$ $\displaystyle=U_{+}e^{i\omega x}+U_{-}e^{-i\omega x}$ (118)
$\displaystyle W$
$\displaystyle=W_{0}e^{-ik\tau_{m}}+i\frac{\bar{E}_{0}}{\omega}U_{+}e^{i\omega
x}-i\frac{\bar{E}_{0}}{\omega}U_{-}e^{-i\omega x}$ (119) $\displaystyle Y$
$\displaystyle=i\frac{\bar{E}_{0}}{\sigma\bar{B}_{0}}W_{0}e^{-ik\tau_{m}}+\sigma\frac{\bar{B}_{0}}{\omega}U_{+}e^{i\omega
x}-\sigma\frac{\bar{B}_{0}}{\omega}U_{-}e^{-i\omega x}$ (120)

where we have three complex coefficients for our six initial conditions. The
corresponding solution to the coupled system is

$\displaystyle u$
$\displaystyle=U_{+c}\cos(\delta_{+})-U_{+s}\sin(\delta_{+})+U_{-c}\cos(\delta_{-})-U_{-s}\sin(\delta_{-})$
(121) $\displaystyle w$
$\displaystyle=W_{1}\cos(k(\tau-\tau_{m}))-W_{2}\sin(k(\tau-\tau_{m}))+\frac{\bar{E}_{0}}{\omega}\left[-U_{+c}\sin(\delta_{+})-U_{+s}\cos(\delta_{+})+U_{-c}\sin(\delta_{-})+U_{-s}\cos(\delta_{-})\right]$
(122) $\displaystyle y$
$\displaystyle=\frac{\bar{E}_{0}}{\sigma\bar{B}_{0}}\left[-W_{2}\cos(k(\tau-\tau_{m}))-W_{1}\sin(k(\tau-\tau_{m}))\right]+\frac{\sigma\bar{B}_{0}}{\omega}\left[U_{+c}\cos(\delta_{+})-U_{+s}\sin(\delta_{+})-U_{-c}\cos(\delta_{-})+U_{-s}\sin(\delta_{-})\right]$
(123)

which yields the following four matching conditions

$\displaystyle u_{m}$ $\displaystyle=U_{+c}+U_{-c}$ (124) $\displaystyle
du_{m}$ $\displaystyle=-U_{+s}(k+f_{m})-U_{-s}(k-f_{m})$ (125) $\displaystyle
w_{m}$
$\displaystyle=W_{1}+\frac{\bar{E}_{0}}{\omega}\left[U_{-s}-U_{+s}\right]$
(126) $\displaystyle dw_{m}$
$\displaystyle=-kW_{2}+\frac{\bar{E}_{0}}{\omega}\left[-U_{+c}(k+f_{m})+U_{-c}(k-f_{m})\right]$
(127) $\displaystyle y_{m}$
$\displaystyle=-\frac{\bar{E}_{0}}{\sigma\bar{B}_{0}}W_{2}+\frac{\sigma\bar{B}_{0}}{\omega}\left[U_{+c}-U_{-c}\right]$
(128) $\displaystyle dy_{m}$
$\displaystyle=-\frac{\bar{E}_{0}}{\sigma\bar{B}_{0}}kW_{1}+\frac{\sigma\bar{B}_{0}}{\omega}\left[-U_{+s}(k+f_{m})+U_{-s}(k-f_{m})\right]$
(129)

whose solutions give the coefficients

$\displaystyle U_{+c}$
$\displaystyle=\frac{\bar{E}_{0}^{2}(k-f_{m})u_{m}-\omega\bar{E}_{0}dw_{m}+k\bar{B}_{0}\left(\bar{B}_{0}u_{m}+\sigma\omega
y_{m}\right)}{2k\omega^{2}}$ (130) $\displaystyle U_{-c}$
$\displaystyle=\frac{\bar{E}_{0}^{2}(k+f_{m})u_{m}+\omega\bar{E_{0}}dw_{m}+k\bar{B}_{0}\left(\bar{B}_{0}u_{m}-\sigma\omega
y_{m}\right)}{2k\omega^{2}}$ (131) $\displaystyle U_{+s}$
$\displaystyle=\frac{-\bar{B}_{0}^{2}(k-f_{m})du_{m}-k\bar{E}_{0}^{2}du_{m}-\omega(k-f_{m})(k\bar{E}_{0}w_{m}+\sigma\bar{B}_{0}dy_{m})}{2k^{2}\bar{E}_{0}^{2}+2\bar{B}_{0}^{2}(k^{2}-f_{m}^{2})}$
(132) $\displaystyle U_{-s}$
$\displaystyle=\frac{-\bar{B}_{0}^{2}(k+f_{m})du_{m}-k\bar{E}_{0}^{2}du_{m}+\omega(k+f_{m})(k\bar{E}_{0}w_{m}+\sigma\bar{B}_{0}dy_{m})}{2k^{2}\bar{E}_{0}^{2}+2\bar{B}_{0}^{2}(k^{2}-f_{m}^{2})}$
(133) $\displaystyle W_{1}$ $\displaystyle=\frac{\bar{B}_{0}\left[-\sigma
k\bar{E}_{0}\omega
dy_{m}+f_{m}\bar{B}_{0}\bar{E}_{0}du_{m}+(k^{2}-f_{m}^{2})\omega\bar{B}_{0}w_{m}\right]}{\omega\left(k^{2}\bar{E}_{0}^{2}+(k^{2}-f_{m}^{2})\bar{B}_{0}^{2}\right)}$
(134) $\displaystyle W_{2}$ $\displaystyle=\frac{-\bar{B}_{0}\left[\sigma
k\bar{E}_{0}\omega
y_{m}+f_{m}\bar{B}_{0}\bar{E}_{0}u_{m}+\omega\bar{B}_{0}dw_{m}\right]}{k\omega^{3}}.$
(135)

For computing the spectral energy densities, we still have Eq. (91) for the
GW, although the coefficients $U_{+c}$, etc. have changed. The spectral energy
densities for the gauge field system, on the other hand, are slightly
modified. It is convenient to define a new set of coefficients,

$\displaystyle W_{+c}$ $\displaystyle=-\frac{\bar{E}_{0}}{\omega}U_{+s}\qquad
W_{+s}=\frac{\bar{E}_{0}}{\omega}U_{+c}$ (136) $\displaystyle W_{-c}$
$\displaystyle=\frac{\bar{E}_{0}}{\omega}U_{-s}\qquad
W_{-s}=-\frac{\bar{E}_{0}}{\omega}U_{-c}$ $\displaystyle Y_{+c}$
$\displaystyle=\frac{\sigma\bar{B}_{0}}{\omega}U_{+c}\qquad
Y_{+s}=\frac{\sigma\bar{B}_{0}}{\omega}U_{+s}$ $\displaystyle Y_{-c}$
$\displaystyle=-\frac{\sigma\bar{B}_{0}}{\omega}U_{-c}\qquad
Y_{-s}=-\frac{\sigma\bar{B}_{0}}{\omega}U_{-s}$ $\displaystyle Y_{1}$
$\displaystyle=-\frac{\bar{E}_{0}}{\sigma\bar{B}_{0}}W_{2}\qquad
Y_{2}=\frac{\bar{E}_{0}}{\sigma\bar{B}_{0}}W_{1}$

which casts the solutions into similar forms,

$\displaystyle u$
$\displaystyle=U_{+c}\cos(\delta_{+})-U_{+s}\sin(\delta_{+})+U_{-c}\cos(\delta_{-})-U_{-s}\sin(\delta_{-})$
(137) $\displaystyle w$
$\displaystyle=W_{1}\cos(k(\tau-\tau_{m}))-W_{2}\sin(k(\tau-\tau_{m}))+W_{+c}\cos(\delta_{+})-W_{+s}\sin(\delta_{+})+W_{-c}\cos(\delta_{-})-W_{-s}\sin(\delta_{-})$
(138) $\displaystyle y$
$\displaystyle=Y_{1}\cos(k(\tau-\tau_{m}))-Y_{2}\sin(k(\tau-\tau_{m}))+Y_{+c}\cos(\delta_{+})-Y_{+s}\sin(\delta_{+})+Y_{-c}\cos(\delta_{-})-Y_{-s}\sin(\delta_{-}).$
(139)

These resemble the solutions from the triplet scenario, except the gauge field
excitations have acquired homogeneous solutions. These homogeneous solutions
contribute additional terms to the time average of the squared gauge field
mode functions,

$\displaystyle\langle{w^{\prime}}_{k}^{2}\rangle=$
$\displaystyle\langle{w^{\prime}}_{u}^{2}\rangle+\frac{k^{2}(W_{1}^{2}+W_{2}^{2})}{2}+2\frac{kW_{1}}{2}\bigl{[}(k+f)W_{+c}\cos(\omega(x-x_{m}))-(k+f)W_{+s}\sin(\omega(x-x_{m}))$
$\displaystyle+(k-f)W_{-c}\cos(\omega(x-x_{m}))+(k-f)W_{-s}\sin(\omega(x-x_{m}))\bigr{]}+2\frac{kW_{2}}{2}\bigl{[}(k+f)W_{+c}\sin(\omega(x-x_{m}))$
$\displaystyle+(k+f)W_{+s}\cos(\omega(x-x_{m}))-(k-f)W_{-c}\sin(\omega(x-x_{m}))+(k-f)W_{-s}\cos(\omega(x-x_{m}))\bigr{]}$
(140)

where $\langle w_{u}^{\prime 2}\rangle$ schematically indicates usage of the
expression for $u$, Eqn. (91), with the coefficient replacement
$\left\\{U_{+c},U_{+s},U_{-c},U_{-s}\right\\}\rightarrow\left\\{W_{c+},W_{s+},W_{c-},W_{s-}\right\\}$.

## References

  * Guth (1981) A. H. Guth, Phys. Rev. D 23, 347 (1981).
  * Linde (1982) A. D. Linde, Phys. Lett. B 108, 389 (1982).
  * Albrecht and Steinhardt (1982) A. Albrecht and P. J. Steinhardt, Phys. Rev. Lett. 48, 1220 (1982).
  * Mukhanov and Chibisov (1981) V. F. Mukhanov and G. V. Chibisov, JETP Lett. 33, 532 (1981).
  * Hawking (1982) S. W. Hawking, Phys. Lett. B 115, 295 (1982).
  * Guth and Pi (1982) A. H. Guth and S. Y. Pi, Phys. Rev. Lett. 49, 1110 (1982).
  * Starobinsky (1982) A. A. Starobinsky, Phys. Lett. B 117, 175 (1982).
  * Bardeen et al. (1983) J. M. Bardeen, P. J. Steinhardt, and M. S. Turner, Phys. Rev. D 28, 679 (1983).
  * Grishchuk (1975) L. P. Grishchuk, Sov. Phys. JETP 40, 409 (1975).
  * Starobinsky (1979) A. A. Starobinsky, JETP Lett. 30, 682 (1979).
  * Rubakov et al. (1982) V. A. Rubakov, M. V. Sazhin, and A. V. Veryaskin, Phys. Lett. B 115, 189 (1982).
  * Fabbri and Pollock (1983) R. Fabbri and M. d. Pollock, Phys. Lett. B 125, 445 (1983).
  * Abbott and Wise (1984) L. F. Abbott and M. B. Wise, Nucl. Phys. B 244, 541 (1984).
  * Kamionkowski et al. (1997) M. Kamionkowski, A. Kosowsky, and A. Stebbins, Phys. Rev. Lett. 78, 2058 (1997), eprint astro-ph/9609132.
  * Seljak and Zaldarriaga (1997) U. Seljak and M. Zaldarriaga, Phys. Rev. Lett. 78, 2054 (1997), eprint astro-ph/9609169.
  * Ade et al. (2017) P. A. R. Ade et al. (POLARBEAR), Astrophys. J. 848, 121 (2017), eprint 1705.02907.
  * Ade et al. (2018) P. A. R. Ade et al. (BICEP2, Keck Array), Phys. Rev. Lett. 121, 221301 (2018), eprint 1810.05216.
  * Ade et al. (2019) P. Ade et al. (Simons Observatory), JCAP 02, 056 (2019), eprint 1808.07445.
  * Akrami et al. (2020) Y. Akrami et al. (Planck), Astron. Astrophys. 641, A10 (2020), eprint 1807.06211.
  * Tristram et al. (2020) M. Tristram et al. (2020), eprint 2010.01139.
  * Abazajian et al. (2016) K. N. Abazajian et al. (CMB-S4) (2016), eprint 1610.02743.
  * Hazumi et al. (2019) M. Hazumi et al., J. Low Temp. Phys. 194, 443 (2019).
  * Spergel and Zaldarriaga (1997) D. N. Spergel and M. Zaldarriaga, Phys. Rev. Lett. 79, 2180 (1997), eprint astro-ph/9705182.
  * Kamionkowski and Kosowsky (1998) M. Kamionkowski and A. Kosowsky, Phys. Rev. D 57, 685 (1998), eprint astro-ph/9705219.
  * Knox and Song (2002) L. Knox and Y.-S. Song, Phys. Rev. Lett. 89, 011303 (2002), eprint astro-ph/0202286.
  * Baumann et al. (2015) D. Baumann, D. Green, and R. A. Porto, JCAP 01, 016 (2015), eprint 1407.2621.
  * Kamionkowski and Kovetz (2016) M. Kamionkowski and E. D. Kovetz, Ann. Rev. Astron. Astrophys. 54, 227 (2016), eprint 1510.06042.
  * Lyth (1997) D. H. Lyth, Phys. Rev. Lett. 78, 1861 (1997), eprint hep-ph/9606387.
  * Anber and Sorbo (2010) M. M. Anber and L. Sorbo, Phys. Rev. D 81, 043534 (2010), eprint 0908.4089.
  * Adshead and Wyman (2012) P. Adshead and M. Wyman, Phys. Rev. Lett. 108, 261302 (2012), eprint 1202.2366.
  * Maleknejad and Sheikh-Jabbari (2013) A. Maleknejad and M. M. Sheikh-Jabbari, Phys. Lett. B 723, 224 (2013), eprint 1102.1513.
  * Maleknejad et al. (2013) A. Maleknejad, M. M. Sheikh-Jabbari, and J. Soda, Phys. Rept. 528, 161 (2013), eprint 1212.2921.
  * Maleknejad et al. (2018) A. Maleknejad, M. Noorbala, and M. M. Sheikh-Jabbari, Gen. Rel. Grav. 50, 110 (2018), eprint 1208.2807.
  * Adshead et al. (2013a) P. Adshead, E. Martinec, and M. Wyman, JHEP 09, 087 (2013a), eprint 1305.2930.
  * Adshead et al. (2016a) P. Adshead, E. Martinec, E. I. Sfakianakis, and M. Wyman, JHEP 12, 137 (2016a), eprint 1609.04025.
  * Adshead and Sfakianakis (2017) P. Adshead and E. I. Sfakianakis, JHEP 08, 130 (2017), eprint 1705.03024.
  * Dimastrogiovanni et al. (2018) E. Dimastrogiovanni, M. Fasiello, R. J. Hardwick, H. Assadullahi, K. Koyama, and D. Wands, JCAP 11, 029 (2018), eprint 1806.05474.
  * Domcke et al. (2019a) V. Domcke, B. Mares, F. Muia, and M. Pieroni, JCAP 04, 034 (2019a), eprint 1807.03358.
  * Domcke et al. (2019b) V. Domcke, B. von Harling, E. Morgante, and K. Mukaida, JCAP 10, 032 (2019b), eprint 1905.13318.
  * Adshead et al. (2015) P. Adshead, J. T. Giblin, T. R. Scully, and E. I. Sfakianakis, JCAP 12, 034 (2015), eprint 1502.06506.
  * Adshead et al. (2017) P. Adshead, J. T. Giblin, and Z. J. Weiner, Phys. Rev. D 96, 123512 (2017), eprint 1708.02944.
  * Adshead et al. (2018) P. Adshead, J. T. Giblin, and Z. J. Weiner, Phys. Rev. D 98, 043525 (2018), eprint 1805.04550.
  * Adshead et al. (2020a) P. Adshead, J. T. Giblin, M. Pieroni, and Z. J. Weiner, Phys. Rev. Lett. 124, 171301 (2020a), eprint 1909.12843.
  * Adshead et al. (2020b) P. Adshead, J. T. Giblin, M. Pieroni, and Z. J. Weiner, Phys. Rev. D 101, 083534 (2020b), eprint 1909.12842.
  * Anber and Sorbo (2012) M. M. Anber and L. Sorbo, Phys. Rev. D 85, 123537 (2012), eprint 1203.5849.
  * Adshead et al. (2013b) P. Adshead, E. Martinec, and M. Wyman, Phys. Rev. D 88, 021302 (2013b), eprint 1301.2598.
  * Namba et al. (2013) R. Namba, E. Dimastrogiovanni, and M. Peloso, JCAP 11, 045 (2013), eprint 1308.1366.
  * Maleknejad (2016a) A. Maleknejad, JHEP 07, 104 (2016a), eprint 1604.03327.
  * Maleknejad (2016b) A. Maleknejad (2016b), eprint 1612.05701.
  * Dimastrogiovanni et al. (2017) E. Dimastrogiovanni, M. Fasiello, and T. Fujita, JCAP 01, 019 (2017), eprint 1608.04216.
  * Fujita et al. (2018) T. Fujita, R. Namba, and Y. Tada, Phys. Lett. B 778, 17 (2018), eprint 1705.01533.
  * Caldwell and Devulder (2018) R. R. Caldwell and C. Devulder, Phys. Rev. D 97, 023532 (2018), eprint 1706.03765.
  * Agrawal et al. (2018a) A. Agrawal, T. Fujita, and E. Komatsu, Phys. Rev. D 97, 103526 (2018a), eprint 1707.03023.
  * Agrawal et al. (2018b) A. Agrawal, T. Fujita, and E. Komatsu, JCAP 06, 027 (2018b), eprint 1802.09284.
  * Fujita et al. (2019) T. Fujita, E. I. Sfakianakis, and M. Shiraishi, JCAP 05, 057 (2019), eprint 1812.03667.
  * Watanabe and Komatsu (2020) Y. Watanabe and E. Komatsu (2020), eprint 2004.04350.
  * Davis et al. (2001) A.-C. Davis, K. Dimopoulos, T. Prokopec, and O. Tornkvist, Phys. Lett. B 501, 165 (2001), eprint astro-ph/0007214.
  * Dimopoulos et al. (2002) K. Dimopoulos, T. Prokopec, O. Tornkvist, and A. C. Davis, Phys. Rev. D 65, 063505 (2002), eprint astro-ph/0108093.
  * Demozzi et al. (2009) V. Demozzi, V. Mukhanov, and H. Rubinstein, JCAP 08, 025 (2009), eprint 0907.1030.
  * Adshead et al. (2016b) P. Adshead, J. T. Giblin, T. R. Scully, and E. I. Sfakianakis, JCAP 10, 039 (2016b), eprint 1606.08474.
  * Kandus et al. (2011) A. Kandus, K. E. Kunze, and C. G. Tsagas, Phys. Rept. 505, 1 (2011), eprint 1007.3891.
  * Cui et al. (2018) Y. Cui, M. Lewicki, D. E. Morrissey, and J. D. Wells, Phys. Rev. D 97, 123505 (2018), eprint 1711.03104.
  * Allahverdi et al. (2021) R. Allahverdi et al., Open J. Astrophys. 4 (2021), eprint 2006.16182.
  * Seto and Yokoyama (2003) N. Seto and J. Yokoyama, J. Phys. Soc. Jap. 72, 3082 (2003), eprint gr-qc/0305096.
  * Boyle and Steinhardt (2008) L. A. Boyle and P. J. Steinhardt, Phys. Rev. D 77, 063504 (2008), eprint astro-ph/0512014.
  * Watanabe and Komatsu (2006) Y. Watanabe and E. Komatsu, Phys. Rev. D 73, 123515 (2006), eprint astro-ph/0604176.
  * Nakayama et al. (2008) K. Nakayama, S. Saito, Y. Suwa, and J. Yokoyama, Phys. Rev. D 77, 124001 (2008), eprint 0802.2452.
  * Kuroyanagi et al. (2011) S. Kuroyanagi, K. Nakayama, and S. Saito, Phys. Rev. D 84, 123513 (2011), eprint 1110.4169.
  * Saikawa and Shirai (2018) K. Saikawa and S. Shirai, JCAP 05, 035 (2018), eprint 1803.01038.
  * Caldwell et al. (2019) R. R. Caldwell, T. L. Smith, and D. G. E. Walker, Phys. Rev. D 100, 043513 (2019), eprint 1812.07577.
  * Figueroa and Tanin (2019) D. G. Figueroa and E. H. Tanin, JCAP 08, 011 (2019), eprint 1905.11960.
  * Bielefeld and Caldwell (2015) J. Bielefeld and R. R. Caldwell, Phys. Rev. D 91, 124004 (2015), eprint 1503.05222.
  * Aghanim et al. (2020) N. Aghanim et al. (Planck), Astron. Astrophys. 641, A6 (2020), eprint 1807.06209.
  * Fields et al. (2020) B. D. Fields, K. A. Olive, T.-H. Yeh, and C. Young, JCAP 03, 010 (2020), [Erratum: JCAP 11, E02 (2020)], eprint 1912.01132.
  * Gertsenshteyn (1962) M. E. Gertsenshteyn, Sov. Phys. JETP 14, 84 (1962).
  * Poznanin (1969) P.-L. Poznanin, Sov. Phys. J. 12, 1296 (1969).
  * Boccaletti et al. (1970) D. Boccaletti, V. Sabbata, P. Fortini, and C. Gualdi, Nuovo Cimento V Serie 70, 129 (1970).
  * Zeldovich (1974) Y. B. Zeldovich, Sov. Phys. JETP 38, 652 (1974).
  * Caldwell et al. (2016) R. R. Caldwell, C. Devulder, and N. A. Maksimova, Phys. Rev. D 94, 063005 (2016), eprint 1604.08939.
  * Gluscevic and Kamionkowski (2010) V. Gluscevic and M. Kamionkowski, Phys. Rev. D 81, 123529 (2010), eprint 1002.1308.
  * Gerbino et al. (2016) M. Gerbino, A. Gruppuso, P. Natoli, M. Shiraishi, and A. Melchiorri, JCAP 07, 044 (2016), eprint 1605.09357.
  * Thorne et al. (2018) B. Thorne, T. Fujita, M. Hazumi, N. Katayama, E. Komatsu, and M. Shiraishi, Phys. Rev. D 97, 043506 (2018), eprint 1707.03240.
  * Thrane and Romano (2013) E. Thrane and J. D. Romano, Phys. Rev. D 88, 124032 (2013), eprint 1310.5300.
  * Smith and Caldwell (2019) T. L. Smith and R. Caldwell, Phys. Rev. D 100, 104055 (2019), eprint 1908.00546.
  * Schmitz (2021) K. Schmitz, JHEP 01, 097 (2021), eprint 2002.04615.
  * Amaro-Seoane et al. (2017) P. Amaro-Seoane et al. (LISA) (2017), eprint 1702.00786.
  * Abbott et al. (2017) B. P. Abbott et al. (LIGO Scientific), Class. Quant. Grav. 34, 044001 (2017), https://dcc.ligo.org/LIGO-P1600143/public, eprint 1607.08697.
  * Hild et al. (2011) S. Hild et al., Class. Quant. Grav. 28, 094013 (2011), eprint 1012.0908.
  * (89) M. Punturo, _ET sensitivities page_ , http://www.et-gw.eu/index.php/etsensitivities.
  * Crowder and Cornish (2005) J. Crowder and N. J. Cornish, Phys. Rev. D 72, 083005 (2005), eprint gr-qc/0506015.
  * Harry et al. (2006) G. M. Harry, P. Fritschel, D. A. Shaddock, W. Folkner, and E. S. Phinney, Class. Quant. Grav. 23, 4887 (2006), [Erratum: Class.Quant.Grav. 23, 7361 (2006)].
  * Maleknejad and Sheikh-Jabbari (2011) A. Maleknejad and M. M. Sheikh-Jabbari, Phys. Rev. D 84, 043515 (2011), eprint 1102.1932.
  * Ma and Bertschinger (1995) C.-P. Ma and E. Bertschinger, Astrophys. J. 455, 7 (1995), eprint astro-ph/9506072.
  * Mukhanov et al. (1992) V. F. Mukhanov, H. A. Feldman, and R. H. Brandenberger, Phys. Rept. 215, 203 (1992).