USTC-ICTS/PCFT-21-46

# Gravitational Waves from Fully General Relativistic Oscillon Preheating

Xiao-Xiao Kou<EMAIL_ADDRESS>Interdisciplinary Center for
Theoretical Study, University of Science and Technology of China, Hefei, Anhui
230026, China  Peng Huanwu Center for Fundamental Theory, Hefei, Anhui 230026,
China  James B. Mertens<EMAIL_ADDRESS>Department of Physics and
McDonnell Center for the Space Sciences, Washington University, St. Louis, MO
63130, USA  Chi Tian<EMAIL_ADDRESS>School of Physics and Optoelectronics
Engineering, Anhui University, Hefei, Anhui 230601, China  Department of
Physics and McDonnell Center for the Space Sciences, Washington University,
St. Louis, MO 63130, USA  Shuang-Yong Zhou<EMAIL_ADDRESS>Interdisciplinary Center for Theoretical Study, University of Science and
Technology of China, Hefei, Anhui 230026, China  Peng Huanwu Center for
Fundamental Theory, Hefei, Anhui 230026, China

###### Abstract

As long-lived quasi-solitons from the fragmentation of a scalar condensate,
oscillons may dominate the preheating era after inflation. During this period,
stochastic gravitational waves can also be generated. We quantify the
gravitational wave production in this period with simulations accounting for
full general relativity to capture all possible non-perturbative effects. We
compute the gravitational wave spectra across a range of choices of the
oscillon preheating models and compare our results to a conventional
perturbative approach on an FLRW background. We clarify the gauge ambiguities
in computing induced gravitational waves from scenarios where dense non-
perturbative objects such as oscillons are being formed. In particular, we
find that the synchronous gauge tends to contain large artificial enhancements
in the gravitational wave spectrum due to gauge modes if gravity plays an
important role in the formation of the oscillons, while other gauge choices,
such as the radiation gauge or a suitably chosen “1+log” gauge, can
efficiently reduce the contributions of gauge modes. The full general
relativistic simulations indicate that gravitational wave spectra obtained
from the perturbative approach on the FLRW background are fairly accurate,
except when oscillon formation induces strong gravitational effects, for which
case there can be an order unity enhancement.

##  I Introduction

During the preheating period after inflation, when the inflaton field still
oscillates near the minimum of its potential, strong parametric resonance can
trigger the exponential growth of fluctuations and quickly push the universe
into a highly nonlinear state. In such a phase, fragmentation of the inflaton
condensate can generate copious non-perturbative structures, including
oscillons Bogolyubsky (1976); Gleiser (1994); Honda and Choptuik (2002); Amin
_et al._ (2012, 2010); McDonald (2002); Broadhead and McDonald (2005);
Copeland  _et al._ (2014); Xie  _et al._ (2021). For a large class of
inflationary models, oscillons may account for the majority of the energy
density during the preheating phase Amin  _et al._ (2012), which is referred
to as the oscillon preheating scenario.

During the oscillon preheating phase, ample relics, such as primordial black
holes or stochastic gravitational waves (GWs) can be produced Zhou  _et al._
(2013); Amin  _et al._ (2018); Liu  _et al._ (2018); Antusch  _et al._ (2018a,
b, 2017); Zhou (2015); Sang and Huang (2019); Cotner  _et al._ (2018, 2019);
Kou  _et al._ (2021); Nazari  _et al._ (2021); Hiramatsu  _et al._ (2021).
While the potential for observable signatures of these relics makes the
studies of such phenomena extremely fascinating, the high degree of non-
linearity during this phase makes analytical studies very challenging. Thus,
numerical approaches that perturbatively solve the Einstein field equations
around an FLRW background have been developed and used to quantify the
generation of oscillons and their relics Gleiser (1994); Copeland  _et al._
(1995); Honda and Choptuik (2002); Copeland  _et al._ (2002); Farhi  _et al._
(2005); Gleiser (2007); Graham and Stamatopoulos (2006); Hindmarsh and Salmi
(2006); Fodor  _et al._ (2006); Saffin and Tranberg (2007); Graham (2007);
Farhi  _et al._ (2008); Hindmarsh and Salmi (2008); Amin  _et al._ (2012);
Zhou  _et al._ (2013); Amin  _et al._ (2018); Liu  _et al._ (2018); Antusch
_et al._ (2018a, b, 2017); Zhou (2015); Sang and Huang (2019); Cotner  _et
al._ (2018, 2019); Kou  _et al._ (2021); Nazari  _et al._ (2021); Hiramatsu
_et al._ (2021). More recently, non-perturbative approaches using techniques
from numerical relativity have been applied to preheating scenarios in order
to capture all possible nonlinear physics Bastero-Gil  _et al._ (2008, 2010a);
Giblin and Tishue (2019); Kou  _et al._ (2021); Nazari  _et al._ (2021). It
has been shown that solving the full Einstein equations will result in an
enhancement in the oscillon abundance in models where oscillons are formed
with significant help from gravitational attraction, and can even lead to the
formation of primordial black holes Kou  _et al._ (2021); Nazari  _et al._
(2021).

Recently, an increasing amount of attention has been devoted towards better
understanding ambiguities in computing the GW power spectra in different
gauges. Those ambiguities originate from the fact that the energy content of
GWs is, when working beyond linear order, gauge dependent. This problem was
first pointed out in Matarrese  _et al._ (1998), and has more recently been
discussed in Hwang  _et al._ (2017), followed by extensive studies Domènech
and Sasaki (2018); Gong (2019); Tomikawa and Kobayashi (2020); De Luca  _et
al._ (2020); Inomata and Terada (2020); Yuan  _et al._ (2020); Lu  _et al._
(2020); Ali  _et al._ (2021); Domènech and Sasaki (2021); Gurian  _et al._
(2021). However, the discussions so far focus on tensor modes induced by
second order scalar perturbations around a homogeneous cosmological
background. It remains an open question how gauge choice can effect the
evaluation of the GW spectrum in the early universe where large
inhomogeneities and non-linearities are present in the form of localized and
dense structures such as oscillons.

Figure 1: From left to right: oscillon formation in simulations with Model I,
Model II and Model III respectively, which correspond to different choices of
potential parameters $\alpha$ and $\beta$ (see Sec. IV). Contours of a
constant field value are shown in a spatial volume at time $t=187m^{-1}$.
Simulations were run using the “1+log” gauge choice with $\eta=0.01$.

In this paper, we analyze the GWs generated during oscillon preheating with
numerical simulations utilizing the full power of numerical relativity. We
will see that a perturbative approach on an FLRW background turns out to be a
relatively good approximation for the class of models we consider, except when
gravitational effects strongly impact oscillon formation, which may lead to
enhancements up to about a factor of a few. We will discuss the gauge
ambiguities in extracting the GW spectra from the simulations with full
general relativity (GR), which is important in the models when oscillons are
formed in strong gravity. Our results suggest that the “1+log” or the
radiation gauge is a more appropriate choice than synchronous gauge when
computing GWs beyond linear order, as large artificial enhancements may arise
due to gauge modes in synchronous gauge.

The paper is organized as follows. In Section II, we will introduce the
oscillon preheating model we consider and the methodology we use to
numerically solve the system and compute the GWs generated during oscillon
preheating. In Section III, we introduce the gauge conditions we employ with
our numerical GR simulations. The main results and discussions are presented
in Section IV, where we show the spectra of the GWs computed in different
gauges and compare them with the ones calculated with the perturbative method
on an FLRW background. We present a semi-quantitative analysis to explain such
discrepancies between different gauges in Section V. We summarize and conclude
in Section VI.

##  II Oscillon Preheating and Models

As localized quasi-stationary field configurations, oscillons can emerge due
to field dynamics in models containing a scalar field with a shallow
potential. We will focus on a representative class of minimally coupled models
given by the action

$\displaystyle S=\int\mathrm{d}^{4}x\sqrt{-g}\left(\frac{M_{\rm
pl}^{2}}{2}R-\frac{1}{2}\partial_{\mu}\phi\partial^{\mu}\phi-V(\phi)\right),$
(1)

with a class of potentials motivated by inflationary model constructions in
string/M theory Silverstein and Westphal (2008); McAllister  _et al._ (2010);
Dong  _et al._ (2011),

$\displaystyle
V(\phi)=\frac{m^{2}M^{2}}{2\alpha}\left[\left(1+\frac{\phi^{2}}{M^{2}}\right)^{\alpha}-1\right].$
(2)

To facilitate our analysis of the parameter space, we parametrize the
potential with the exponent $\alpha$ and the ratio $\beta\equiv M_{\rm pl}/M$,
where $M_{\rm pl}$ is the reduced Planck mass. The range of parameters
$0<\alpha<1$ and $0<\beta<100$ is considered. Within this range, the potential
is flatter than the quadratic mass term at large $\phi$ and approximates to
the quadratic term at smaller $\phi$. This type of potential supports quasi-
stable oscillon configurations, as the interacting potential,
$V(\phi)-\frac{1}{2}m^{2}\phi^{2}$, is negative and thus attractive. In other
words, “particles” under the influence of this kind of potential prefer to
condense to localized objects rather than propagate to infinity, which allows
long-lived oscillons to form from quite generic initial conditions. A
visualization of oscillon formation during the preheating period is shown in
Fig. 1.

The value of scalar field mass $m$ can be estimated by matching to the power
spectrum of curvature perturbations from most recent CMB observations Aghanim
_et al._ (2020) through the relation

$\displaystyle A_{s}=\frac{(4\alpha{\cal
N})^{1+\alpha}}{96\pi^{2}\alpha^{3}}\left(\frac{m}{M_{\rm
pl}}\right)^{2}\\!\\!\left(\frac{M}{M_{\rm pl}}\right)^{2-2\alpha}\\!\\!\simeq
2.1\\!\times\\!10^{-9},$ (3)

where ${\cal N}\simeq 50$ is chosen to be the fiducial number of e-foldings
before inflation ends.

Let us first review how to extract the stochastic GW power spectrum in an FLRW
background. Oscillations of the inflaton field obey the Klein-Gordon equation.
When coupled to a homogeneous and isotropic cosmological background, the
Klein-Gordon equation can be simplified to

$\displaystyle\ddot{\phi}+3H\dot{\phi}-\frac{\nabla^{2}\phi}{a^{2}}+\frac{{\rm
d}V(\phi)}{\rm{d}\phi}=0,$ (4)

where $H$ is the Hubble parameter and $a$ is the scale factor.

The GWs generated around the FLRW background can be quantified by analysing
the transverse-traceless part of the metric perturbations, the tensor
perturbation $h_{ij}^{TT}$, sourced by matter fields, that is, by solving the
linearized tensor perturbation equations

$\displaystyle\ddot{h}_{ij}^{TT}$
$\displaystyle+3H\dot{h}_{ij}^{TT}-\frac{1}{a^{2}}\nabla^{2}h_{ij}^{TT}=16\pi
GS_{ij}^{TT},$ (5) $\displaystyle S_{ij}^{TT}$
$\displaystyle=\frac{1}{a^{2}}\left(\partial_{i}\phi\partial_{j}\phi\right)^{TT}.$
(6)

Numerically calculating $S_{ij}^{TT}(t,\mathbf{x})$ at every step in the
simulation is unnecessarily expensive. Rather, due to the linearity of this
equation, one can simply evolve $h_{ij}$, which satisfies

$\displaystyle\ddot{h}_{ij}+3H\dot{h}_{ij}-\frac{1}{a^{2}}\nabla^{2}h_{ij}=16\pi
GS_{ij}$ (7)

and, when needed, project to the physical transverse-traceless part
$h_{ij}^{TT}(t,\mathbf{x})$ with:

$\displaystyle
h_{ij}^{TT}(t,\mathbf{x})=\int\frac{\mathrm{d}^{3}k}{(2\pi)^{3}}e^{i\mathbf{k}\cdot\mathbf{x}}\Lambda_{ij,lm}(\hat{\mathbf{k}})h_{lm}(t,\mathbf{k})$
(8)

where $\Lambda_{ij,lm}(\hat{\mathbf{k}})$ is a projection operator defined as

$\displaystyle\Lambda_{ij,lm}(\hat{\mathbf{k}})\equiv
P_{il}(\hat{\mathbf{k}})P_{jm}(\hat{\mathbf{k}})-\frac{1}{2}P_{ij}(\hat{\mathbf{k}})P_{lm}(\hat{\mathbf{k}})$
(9)

with $\hat{\mathbf{k}}=\mathbf{k}/k$ and
$P_{ij}(\hat{\mathbf{k}})=\delta_{ij}-\hat{k}_{i}\hat{k}_{j}$. To guarantee
tracelessness and transversality in the projection, we use the projector based
on nearest-neighbour spatial derivatives suggested by Figueroa  _et al._
(2011) after trying different projectors (see Appendix B for details of tests
for different projectors). The energy density of GWs can be extracted by
ensemble averaging of several wavelengths

$\displaystyle\rho_{GW}$ $\displaystyle=\frac{1}{32\pi
G}\left\langle\dot{h}_{ij}^{TT}\dot{h}^{ijTT}\right\rangle.$ (10)

The GW energy density spectrum per logarithmic interval is defined as

$\displaystyle\rho_{\mathrm{GW}}(t)$
$\displaystyle\equiv\int\frac{d\rho_{\mathrm{GW}}}{d\log k}d\log k,$ (11)
$\displaystyle\frac{d\rho_{\mathrm{GW}}}{d\log k}$
$\displaystyle=\frac{k^{3}}{(4\pi)^{3}GV}\int\frac{d\Omega_{k}}{4\pi}\dot{h}_{ij}^{TT}(\mathbf{k},t)\dot{h}_{ij}^{*TT}(\mathbf{k},t),$
(12)

where $\Omega_{k}$ represents the solid angle in $\bm{k}$-space, and the
volume $V=L^{3}$ in the lattice implementation. It is customary to normalize
the GW energy density spectrum with the critical density $\rho_{c}\equiv
3H_{0}^{2}/8\pi G$, leading to

$\displaystyle\Omega_{\mathrm{GW}}(k,t)\equiv\frac{1}{\rho_{c}}\frac{d\rho_{\mathrm{GW}}}{d\log
k}$ (13)

This is the conventional perturbative approach on an FLRW background, which
involves relatively simple differential equations and where the non-linear
dynamics of oscillons and their GW generation can be well resolved with fairly
low computational costs. However, in this approach, gravitational backreaction
is neglected, which a priori may not be a valid approximation for simulations
when gravitational effects are strong during the oscillon formation. Indeed,
it was found in Kou  _et al._ (2021) that, with the full GR approach, the
abundance of oscillons in preheating is enhanced for those circumstances.

We shall investigate whether solving full Einstein equations might introduce
corrections to the relic GW spectra. To this end, we perform the 3+1
decomposition of the Einstein equations, which decomposes the metric into the
form

$\displaystyle g_{\mu\nu}=\begin{bmatrix}-N^{2}+N_{k}N^{k}&N_{j}\\\
N_{i}&\gamma_{ij}\end{bmatrix},$ (14)

where $N$ and $N^{i}$ are the lapse and shift respectively and $\gamma_{ij}$
is the 3D metric on the spatial hypersurfaces foliated by $N$ and $N^{i}$. The
extrinsic curvature tensor $K_{ij}$ describing the embedding of the spatial
hypersurfaces is defined as the Lie derivative of $\gamma_{ij}$ with respect
to the normal vector of the hypersurfaces, namely

$\displaystyle K_{ij}=\mathcal{L}_{n}\gamma_{ij}.$ (15)

We also cast the Klein-Gordon equation coupled to GR into the hyperbolic form

$\displaystyle\dot{\phi}$ $\displaystyle=N^{i}\partial_{i}\phi-N\Pi$ (16)
$\displaystyle\dot{\Pi}$
$\displaystyle=N^{k}\partial_{k}\Pi-N\gamma^{ij}\partial_{i}\psi_{j}+N\gamma^{ij}\Gamma^{k}_{ij}\psi_{k}+NK\Pi$
$\displaystyle-\gamma^{ij}\psi_{i}\partial_{j}N+N{\rm d}V(\phi)/{\rm d}\phi$
(17) $\displaystyle\dot{\psi}_{i}$
$\displaystyle=N^{j}\partial_{j}\psi_{i}+\psi_{j}\partial_{i}N^{j}-N\partial_{i}\Pi-\Pi\partial_{i}N,$
(18)

where $\psi_{i}$ is the auxiliary field $\psi_{i}\equiv\partial_{i}\phi$, $K$
is the trace of the extrinsic curvature and $\Pi$ is the canonical momentum of
the scalar field

$\displaystyle\Pi\equiv-\left(\frac{1}{N}\dot{\phi}-\frac{1}{N}N^{k}\partial_{k}\phi\right).$
(19)

We will from now on refer to this scheme as the GR scheme, and the scheme
relying on Eq. (4) as the FLRW scheme.

To compute the energy of GWs in full GR and compare with the FLRW scheme, we
need a background expanding metric. A nature choice is to average the spatial
hypersurface to get the background, then subtract the expanding background
from the full metric Isaacson (1968a, b); Bastero-Gil  _et al._ (2010b); Ota
_et al._ (2021). This allows us to define

$\displaystyle\rho_{GW}$ $\displaystyle=\frac{1}{32\pi
G}\left\langle{\hat{\gamma}}_{ij;0}^{TT}{\hat{\gamma}}^{ijTT}_{;0}\right\rangle,$
(20)

where ;0 means covariant time derivative compatible with background metric.
And $\hat{\gamma}_{ij}$ is defined as

$\hat{\gamma}_{ij}\equiv\chi^{-2}\tilde{\gamma}_{ij}-a^{2}(t)\delta_{ij},$
(21)

where $\chi$ is the conformal factor and $\tilde{\gamma}_{ij}$ is the
conformal metric defined by decomposing the spatial metric $\gamma_{ij}$ as
$\tilde{\gamma}_{ij}\equiv\chi^{2}\gamma_{ij}$, with the constraint ${\rm
det}\tilde{\gamma}_{ij}=1$. When the metric is conformally flat
($\tilde{\gamma}_{ij}=\delta_{ij}$) and $\chi$ has no spatial dependence,
$\chi^{-1}$ is equivalent to the scale factor $a$. However, the inhomogeneous
distribution of the scalar field will introduce small spatial dependence to
$\chi$ through the Einstein equations. Its average still has the time
evolution that approximately corresponds to the scale factor in the FLRW
universe, which allows us to use $\langle\chi^{-1}\rangle$ to mimic the scale
factor in Eq. (21). Note that $\hat{\gamma}_{ij;0}^{TT}$ essentially reduces
to the $2K_{ab}^{\mathrm{TT}}$ of Ota  _et al._ (2021) when the fluctuations
in the field $\chi$ are very small, and these small fluctuations can be
initially fixed by solving the constraints equations. Also note that non-
perturbatively defining a GW energy that is gauge independent and can be
implemented in a grid simulation scheme is highly nontrivial. Indeed, the
$\rho_{GW}$ defined above, as we will see later, is gauge dependent. As we
will explain later, we will need to choose gauge conditions judiciously so as
to minimize the gauge effects from the redundant gauge modes. In particular,
despite the background being in synchronous gauge, we will see later that the
synchronous gauge may not be a good gauge choice for the perturbations on top
of the background.

To solve the metric components together with the scalar field they couple to,
we employ the grid-based numerical relativity code CosmoGRaPH Mertens  _et
al._ (2016), which makes use of the popular Baumgarte–Shapiro–Shibata–Nakamura
(BSSN) formalism Nakamura  _et al._ (1987); Shibata and Nakamura (1995);
Baumgarte and Shapiro (1998) and integrates the Adaptive Mesh Refinement (AMR)
framework into its spatial grid scheme. However, for simulations performed in
this work, only uniform grid configuration is employed.

Different from the FLRW scheme, solving the equations in the GR scheme
requires obtaining exact solutions to the constraint equations in the initial
condition setup. We will seek solutions that mostly resemble the FLRW metric
to facilitate our comparisons between the two schemes. To this end, after
splitting the extrinsic curvature into a traceless part and a trace part
$K_{ij}\equiv A_{ij}+\frac{1}{3}\gamma_{ij}K$, we choose the initial
conditions such that $\tilde{\gamma}_{ij}=\delta_{ij}$, $A_{ij}=0$. Since
$K=-3H$, the initial homogeneous $K$ is set such that the corresponding Hubble
parameter satisfies the Friedmann equation. The initial configuration of the
scalar field is chosen to have a standard spectrum of the initial vacuum
fluctuations for $\phi=M$ and $\dot{\phi}=0$, for which the momentum
constraint is trivialized, and the Hamitonian constraint equation can be
written as

$\nabla^{2}\Psi+\pi\left(\partial_{i}\phi\partial_{i}\phi\right)\Psi+\left(\pi
V(\phi)-K^{2}/12\right)\Psi^{5}=0,$ (22)

where $\Psi\equiv\chi^{-1/2}$. We then solve this non-linear Hamiltonian
constraint equation and get the conformal factor $\chi$ by employing the
multigrid constraint solver integrated within CosmoGRaPH. We use a periodic
box with a spatial size $L=50\,m^{-1}$. To ensure numerical convergence of the
elliptic constraint solver, a cut-off at wavenumber $k=12\,(2\pi/L)$ is
implemented.

Figure 2: Time evolution of GW power spectra in the FLRW scheme, with darker
colours representing later times. Model I to III are presented from left to
right.

##  III Gauge choices

Even though the perturbative definition of GW energy Eq. (10) is a gauge
invariant quantity to linear order (when the box size is larger than the
wavelengths of GWs, see Misner  _et al._ (1973)), it becomes gauge dependent
when higher order perturbations are included. This gauge dependence is by
construction absent from the FLRW scheme, as the Klein-Gordon equation is only
coupled to the background cosmology (Eq. (4)). We will consider multiple gauge
choices in the GR scheme to assess how large the gauge dependence is for our
simulations. By comparing different gauges, we look for gauge choices that are
insensitive to the gauge modes from higher order perturbations, which will
allow us to establish discrepancies between the GR and the FLRW scheme.

The simplest gauge condition we consider here is the synchronous gauge, or
geodesic slicing, which simply sets the lapse $N=1$ and a vanishing $N^{i}$.
Indeed, our background is in the synchronous gauge. Despite its simplicity,
the synchronous gauge tends to focus geodesics and create coordinate
singularities very easily. In addition, it is well known that the synchronous
gauge is not a complete gauge fixing, even to leading order. As we shall see
shortly, we find this gauge generally an unsuitable choice for extracting GW
spectra from oscillon preheating.

Another widely used set of gauge conditions in numerical relativity is the
combination of the “1+log” and the “Gamma-driver” shift conditions:

$\displaystyle\partial_{t}N$ $\displaystyle=$ $\displaystyle-2\eta
N\left(K-\langle K\rangle\right)+N^{i}\partial_{i}N,$ (23)
$\displaystyle\partial_{t}N^{i}$ $\displaystyle=$ $\displaystyle
B^{i},\;\;\;\;\partial_{t}B^{i}=\frac{3}{4}\partial_{t}\tilde{\Gamma}^{i}-B^{i},$
(24)

where $\eta$ is a constant to be chosen, $B^{i}$ is only an auxiliary field
and
$\tilde{\Gamma}^{i}\\!\equiv\\!\tilde{\gamma}^{jk}\tilde{\Gamma}_{jk}^{i}$,
with $\tilde{\Gamma}_{jk}^{i}$ being the Christoffel symbols of
$\tilde{\gamma}_{ij}$. We will refer to this gauge condition as the “1+log”
gauge for simplicity. This gauge condition has been proven to have a powerful
singularity-avoiding property, and thus is widely used in solving spacetime
evolutions containing black holes.

We also examine a less-used gauge choice, called the radiation gauge Chen and
Zhu (2011): $g^{ij}\Gamma_{ij}^{\rho}=0$. This gauge condition might look like
the conventional harmonic gauge $g^{\mu\nu}\Gamma_{\mu\nu}^{\rho}=0$ in
appearance, but it is rather different in nature. In fact, it can be proven to
be equivalent to time evolution of the lapse according to:

$\displaystyle\partial_{t}N=-\frac{2\tau N^{2}}{1+N^{2}}\left(K-\langle
K\rangle\right),$ (25)

with a vanishing shift $N^{i}$, where $\tau$ is chosen to be 0.05 in this
paper. Note that this gauge condition has the late-time asymptotic solution
that $K\rightarrow\langle K\rangle$ and $\partial_{t}N=0$, which corresponds
to a homogeneous expansion rate.

##  IV Quantifying the GW spectra

Figure 3: Time evolution (from top to bottom) of GWs spectra for $3$ sets of
Models with different gauge choices. $(a)$ Model I – strong resonance: The top
panel is taken at $t\approx 6m^{-1}$; the middle panel, at $t\approx
94m^{-1}$, is when the simulation in the synchronous gauge breaks down; the
bottom panel, at $t\approx 112m^{-1}$, shows the GW spectra in the final stage
of preheating. $(b)$ Model II – median resonance: The top panel is taken at
$t\approx 8m^{-1}$; the middle panel, at $t\approx 156m^{-1}$ is when the
simulation in the synchronous gauge breaks down; the bottom panel, at
$t\approx 165m^{-1}$, shows the spectrum in the final stage of preheating.
$(c)$ Model III – weak resonance and strong self-gravity: The top panel is
taken at $t\approx 15m^{-1}$; the middle panel, at $t\approx 140m^{-1}$ is
when the simulation in the synchronous gauge breaks down; the bottom panel, at
$t\approx 172m^{-1}$, shows the spectrum in the final stage of preheating. The
insets in the bottom panels depict the ratio between the GW spectrum in the
corresponding gauge condition and that of the FLRW simulation. For Model I, we
run the simulation based on a $192^{3}$ grid, but for Model II and III, we
utilize $256^{3}$ simulations to capture potential self-gravity effects.

We select three representative sets of model parameters to show the influence
of the gauge choice and possible enhancements of GW power spectra compared
with the conventional FLRW based scheme.

Model I | $\alpha=0.18,\beta=50$
---|---
Model II | $\alpha=0.50,\beta=50$
Model III | $\alpha=0.18,\beta=22$

Model I has very strong parametric resonance and generates a large number of
oscillons in the preheating period, Model II corresponds to the case with
slightly weaker parametric resonance, and Model III corresponds to the case
with very weak parametric resonance from the scalar self-interaction but with
strong self-gravity during the oscillon formation 111Some authors may refer to
an oscillon with strong self-gravity as an oscillaton.. Note that even though
the scalar parametric resonance in Model III is weak, strong gravitational
effects can boost the generations of oscillons (see Kou  _et al._ (2021) for
more detail).

For all the 3 sets of parameters of the potential, the inflaton condensate can
efficiently fragment into oscillons; see Fig. 1 for the illustrations of
oscillons for the $3$ models in the GR scheme. During the fragmentation,
significant stochastic GW backgrounds are generated, and their power spectra
keep growing until the oscillons are stabilized. We compute the power spectra
of GWs $\rho_{\rm GW}(k)$ after oscillons form and are stabilized, and the
growth of the GWs spectra for all $3$ models is shown in Fig. 2. Note that at
the earlier stage during the evolution, the high frequency ends of the spectra
have approximated $k$ dependence $\Omega_{\rm GW}\propto k^{3}$, which result
in straight lines in the lower right corner in this logarithmic plot. This is
due to the numerical artifact of applying the cutoff at higher $k$ — the
cutoff results in a small and flat spectrum of $h_{ij}^{\rm TT}$ at higher
$k$, which translates into the relation $\Omega_{\rm GW}\propto k^{3}$
approximately according to Eq. (12).

We depict the evolution of the GW spectra with different gauges for all $3$
models in Fig. 3. The GWs from Model I, which gives rise to strong parametric
resonance and formation of a large number of oscillons in preheating, is
plotted in the left column for 3 time instances in Fig. 3. We can see that,
for the “1+log” gauges and the radiation gauge, the results generally coincide
with those of the FLRW scheme until the wavelengths are significantly smaller
than the box size. The enhancements seen for these small wavelengths are
mostly numerical artifacts due to the lack of spatial resolution for these
modes (see Appendix A for details of the convergence tests). These results are
consistent with the GW spectra computed in Zhou (2015), which uses a pseudo-
spectral method and evolves the metric perturbatively on an FLRW background.
The multi-peak patterns in the spectra are intimately linked to the fact that
oscillons evolve to a spherically symmetric shape after emerging from the
inflaton condensate fragmentation Zhou (2015). We also note that the choice of
the synchronous gauge stands out as an exception, where an order of magnitude
of enhancements are observed before the simulation breaks down at $t\approx
95m^{-1}$ (see the middle panel in the left column in Fig. 3). This suggests
that there are huge redundant gauge modes in evaluating the GW production in
oscillon preheating with the synchronous gauge, meaning that the synchronous
gauge is a rather inappropriate choice for calculating the GW production in
full numerical relativity, at least for oscillon preheating.

For Model II, as shown in the middle column in Fig. 3, similarly as in the
Model I: the GW spectra in “1+log” gauges and the radiation gauge are also
consistent with the GW spectra from the FLRW scheme, within the wavelengths
that are well resolved numerically. Synchronous gauge breaks down at around
$156\,m^{-1}$, where the GW power spectrum is again significantly enhanced due
to large artificial gauge modes.

As suggested by Kou  _et al._ (2021), Model III is subjected to strong
gravitational backreaction, which can only be adequately resolved by the full
GR scheme. As shown in the right column in Fig. 3, apart from similar large
gauge modes in synchronous gauge (red line in the middle panel), the GW
spectra computed in “1+log” gauges are also subjected to some gauge dependence
for the wavelengths that are well resolved in our simulations. We find that
the gauge with $\eta=0.1$ is a better choice compared to that with $\eta=0.01$
and $\eta=0.5$, as it contains relatively fewer gauge modes. The reason for
this is that: when $\eta$ approaches 0, the 1+log gauge choice approaches the
synchronous gauge, which has been shown to have a rather large gauge
redundancy; when $\eta$ is large, on the other hand, it is overdamping the
lapse function. An overdamped lapse function will make the high frequency part
within $\tilde{\gamma}_{ij}$ less stable compared with the choice of
$\eta=0.1$, and these unstable high frequency gauge modes will enter
$\hat{\gamma}_{ij}^{TT}$ within its evolution, as can be seen in the high
frequency peaks for $\eta=0.5$ in the last panel of Fig. 3.

Also, the GW spectra from “1+log” gauge with $\eta=0.1$ is consistent with
radiation gauge with $\tau=0.05$ (cf. Eq. (25)), which suggests that the GW
spectra in these gauges might be taken as a good estimate of the “physical”
spectra in this model. With this interpretation, a small enhancement of a
factor of $3$ in the medium wavelengths can identified.

Similar to the gauge ambiguities of the GW energy density we have shown in the
full GR simulations, the gauge dependence of GWs in higher order cosmological
perturbation theory has been widely discussed in the context of GWs induced by
higher order scalar perturbations Matarrese  _et al._ (1998); Hwang  _et al._
(2017); Tomikawa and Kobayashi (2020); Gong (2019); Lu  _et al._ (2020) (see
Domènech (2021) for a recent comprehensive review). In these scenarios,
synchronous gauge is often found to contain more gauge modes and thus is less
reliable in evaluating quantities related to GWs. Similar gauge issues are
also identified here in the oscillon preheating model we study. We show that
when oscillons are in a stably oscillating mode, some gauges do contain fewer
gauge ambiguities, and are thus more suitable for performing computations and
comparing with observations. On the other hand, as pointed out recently
Inomata and Terada (2020); Yuan  _et al._ (2020); Domènech and Sasaki (2021),
the differences between different gauge choices will eventually decrease and
die down when the sources that induce those GWs decay at a later stage and
gravity becomes weaker. As our observations are always at much later times,
these gauge ambiguities thus may not present any true difficulties in theory,
if we can follow the evolution for a sufficiently long time. For our current
case, this means the time when oscillons decays, which is technically
challenging, as oscillons in our model are very stable. We leave this line of
investigation for future work.

##  V Gauge dependence of GW energy density

In Fig. 3, we can observe at least one order of magnitude enhancement of GW
power spectra in the synchronous gauge before the simulation breaks down,
compared with properly chosen “1+log” gauges and radiation gauge. In this
section, we will show how such big gauge dependence may arise by performing a
semi-quantitative analysis of higher order perturbations in the GWs content.

We want to compare the GW energy density (cf. Eq. (10)) between the two
different gauges, i.e. between two different coordinate systems. Consider
small coordinate transformation
$x^{\alpha}\rightarrow\tilde{x}^{\alpha}=x^{\alpha}+\epsilon\xi^{\alpha}$
which induces gauge transformation

$\displaystyle\tilde{h}_{ij}\simeq
h_{ij}-\partial_{j}\xi_{i}-\partial_{i}\xi_{j}-h_{ik}\partial_{j}\xi^{k}-h_{jk}\partial_{i}\xi^{k}+O(\epsilon^{3}).$
(26)

where we have neglected terms of order $\epsilon^{3}$ and terms that, by
integrating by part, lead to divergences when evaluating the GW energy
density. In the $\tilde{x}^{\alpha}$ coordinates, the terms
$\partial_{j}\xi_{i}+\partial_{i}\xi_{j}$ do not contribute to the GW energy
density, thanks to the transverse-traceless projection on $h_{ij}$, which
means that the GW spectrum is gauge invariant to leading order. However, in
the full GR scheme, higher order terms like $h_{ik}\partial_{j}\xi^{k}$ do
give rise to nonvanishing contributions to the GW energy density in the
$\tilde{x}^{\alpha}$ coordinates, which leads to gauge dependence in the GW
energy density for higher order perturbations. Because the GW energy density
is represented schematically as $\rho_{GW}\sim\partial h\partial h$, the
contributions from the higher order terms in the GW energy content
representing the gauge modes can be estimated by

$\displaystyle\rho_{GW,\rm 2}$ $\displaystyle\sim\partial
h\partial(h\partial\xi)+O(h^{4})$
$\displaystyle\sim(\partial^{2}h)h\left(\partial\xi\right)\sim Th\delta h$
(27)

where we have used the leading order equation of motion of $h_{ij}$ to replace
$\partial^{2}h$ with the matter stress-energy tensor $T_{ij}$, and also
approximated $\partial\xi$ with $h-\tilde{h}$. As a simple estimate of the GW
energy density from the gauge modes, we will simply approximate $\rho_{GW,\rm
g}$ with $\rho^{S}_{GW,\rm 2}=(\sum T)(\sum h)(\sum\delta h)$, where $\sum T$
stands for $\sum_{ij}T_{ij}/6$ and so on. We compare this estimation with the
difference in the GW energy density between the synchronous gauge and the
radiation gauge in Fig. 4, where the green arrows representing the second
order estimations for gauge modes approximately match the difference between
the two gauges. It is also worth noting that Eq. (V) implies that the gauge
modes decay when the source $T$ is suppressed.

Figure 4: Time evolution of GW energy densities for Model I with the
synchronous and the radiation gauge. The inset depicts the estimate of the
second order gauge dependence $\rho_{GW,\rm 2}\sim{T}{h}\delta{h}$ as well as
the difference of the GW energy density between the synchronous and the
radiation gauge. Both energy densities are extracted from the full GR
simulation of Model I and normalized to the critical density at the time when
the simulation with the synchronous gauge breaks down.

##  VI Summary

In this paper, we have studied the production of stochastic GWs during the
oscillon preheating scenario by numerically solving the full Einstein
equations under different gauge conditions. The results show that in the
models where gravitational backreation is negligible during oscillon
preheating, the GWs extracted from full GR simulations coincide with the
results given by the conventional simulations in an FLRW background, except
for the spectrum computed in the synchronous gauge. In the synchronous gauge,
the computed GW power spectrum is typically 1-2 orders of magnitude higher
than that of the FLRW simulation or GR simulation with the “1+log” or the
radiation gauge, which means that significant redundant gauge effects are
present in this gauge, making it an unsuitable gauge condition to be used in
oscillon preheating. On the other hand, for models where gravitational
backreaction is significant, the stochastic GW power spectra are enhanced by
up to a factor of $3$ in the full GR simulations under proper gauge choices.

The GWs generated by the particular single-field inflationary models
considered in this paper have a frequency band beyond the current and next
generation of GW detectors. However, once the role of the inflaton, which
should correctly generate the power spectrum of the observed cosmic microwave
background, is released for the scalar field that generate oscillons, an
oscillon preheating-like scenario can occur at lower energy scales, and the
frequencies of the generalized GW background can be in the bands detectable in
the GW detectors. Potential detection of GW signals for such models can be
used to test general relativity and probe new physics at high energy scales
and in the early universe.

Our study suggests that it is important to choose appropriate gauge conditions
to evaluate GW power spectra in full GR preheating simulations, and the
synchronous gauge appears to be a very unreliable gauge condition to be used
to extract GW spectra in full GR simulations, at least when dense objects such
as oscillons are being formed. Incidentally, Bastero-Gil  _et al._ (2010a)
found an order of magnitude enhancement for the stochastic GW spectrum in a
hybrid inflation model, using the full GR simulations with the synchronous
gauge. Of course, in that model no oscillons or other non-perturbative objects
are forming during the preheating. Nevertheless, it would be interesting to
revisit it and exclude any possible gauge redundancies, which are left for
future work.

## Acknowledgements

We are grateful to Li-ming Cao, Glenn Starkman and John T. Giblin, Jr for
helpful discussions. SYZ acknowledges support from the starting grants from
University of Science and Technology of China under grant No. KY2030000089 and
GG2030040375, and is also supported by National Natural Science Foundation of
China under grant No. 11947301, 12075233 and 12047502, and supported by the
Fundamental Research Funds for the Central Universities under grant No.
WK2030000036. JBM is partly supported by a Department of Energy grant DE-
SC0017987 to the particle astrophysics theory group at WUSTL

##  Appendix A Convergence Tests

In this work, we compute GWs spectra in a periodic box with a spatial size
$L=50/m$. For Model I, we run the simulations using a $192^{3}$ grid, since
self-gravity effects are not strong when the parametric resonance is
significant, and the resolution of $192^{3}$ is enough to capture essential
physics. For Models II and III, we utilize $256^{3}$ simulations to capture
potential self-gravity effects, as fully general relativistic corrections are
often important when the parametric resonance is weak Kou  _et al._ (2021).
The validity of these choices of resolutions is also examined by the
convergence tests.

Here we show the convergence tests for the most demanding case, Model III,
which is the model that needs the most simulation steps and converges the
slowest. The convergence tests for the FLRW, the “1+log” GR gauge, and
radiation GR gauge are shown in Figs 5, 6 and 7 respectively. We See that
generally our simulations are sufficiently convergent for the lower part of
the GW spectrum with the $192^{3}$ resolution. The conservation of the
Hamiltonian constraint is also kept well under control throughout the whole
duration of the simulations, and roughly a 3rd order of convergence is
achieved; see Fig. 8.

Figure 5: GW power spectra in the FLRW scheme under different resolutions.
Figure 6: GW power spectra in the “1+log” gauge with $\eta=0.1$ under
different resolutions. Figure 7: GW power spectra in the radiation gauge under
different resolutions. Figure 8: L2 norms of the Hamiltonian constraint under
different resolutions in the radiation gauge. The model parameters are chosen
as $\alpha=0.18$ and $\beta=22$.

##  Appendix B GWs in the lattice

We now discuss the approach to obtain the GW spectra in the lattice applied to
both the FLRW and GR scheme. In order to derive a discrete expression for Eq.
(12) in a lattice of volume $V=L^{3}$, we need to first specify our convention
for the discrete Fourier transform:

$\displaystyle f(\mathbf{n})$
$\displaystyle=\frac{1}{N^{3}}\sum_{\tilde{n}}e^{-\frac{2\pi
i}{N}\tilde{\mathbf{n}}\mathbf{n}}\tilde{f}(\tilde{\mathbf{n}})$ (28)
$\displaystyle\tilde{f}(\tilde{\mathbf{n}})$
$\displaystyle=\sum_{n}e^{+\frac{2\pi
i}{N}\tilde{\mathbf{n}}\mathbf{n}}f(\mathbf{n})$ (29)

where $\mathbf{n}=\left(n_{1},n_{2},n_{3}\right)$ and
$\tilde{\mathbf{n}}=\left(\tilde{n}_{1},\tilde{n}_{2},\tilde{n}_{3}\right)$
label the lattice sites in configuration space and in the momentum space
respectively. It was shown in Figueroa  _et al._ (2011)

$\displaystyle\left(\frac{d\rho_{GW}}{d\log
k}\right)(\tilde{\mathbf{n}})\equiv\frac{dx^{6}k^{3}(|\tilde{\mathbf{n}}|)}{(4\pi)^{3}GL^{3}}\left\langle\dot{h}_{ij}(|\tilde{\mathbf{n}}|,t)\dot{h}_{ij}^{*}(|\tilde{\mathbf{n}}|,t)\right\rangle_{R(\tilde{\mathbf{n}})}$
(30)

where
$\left\langle\dot{h}_{ij}(|\tilde{\mathbf{n}}|,t)\dot{h}_{ij}^{*}(|\tilde{\mathbf{n}}|,t)\right\rangle_{R(\tilde{\mathbf{n}})}$
is an average over configurations with lattice momenta
$\tilde{\mathbf{n}}^{\prime}\in[|\tilde{\mathbf{n}}|,|\tilde{\mathbf{n}}|+\delta\tilde{n}]$.
As mentioned in the main text, the TT part of tensor perturbation $h_{ij}$ is
obtained in Fourier space via the projector
$\Lambda_{ij,lm}(\hat{\mathbf{k}})$ given in Eq. (9). In the lattice, we want
to construct a projector $\Lambda_{ij,lm}$ that satisfies the ”transversality”
and ”tracelessness” conditions, which may come in different forms Figueroa
_et al._ (2011). Nevertheless, for our case, the differences caused by
choosing different projectors are very small, less than the percentage level
prior to the Nyquist frequency band, as we show in Fig. 9. Thus, we choose the
more stable projectors

$\displaystyle\Lambda_{ij,lm}^{(L)}(\tilde{\mathbf{n}})$ $\displaystyle\equiv
P_{il}^{(L)}(\tilde{\mathbf{n}})P_{jm}^{(L)}(\tilde{\mathbf{n}})-\frac{1}{2}P_{ij}^{(L)}(\tilde{\mathbf{n}})P_{lm}^{(L)}(\tilde{\mathbf{n}})$
(31) $\displaystyle P_{ij}^{(L)}(\tilde{\mathbf{n}})$
$\displaystyle=\delta_{ij}-\frac{k_{i}^{(L)}k_{j}^{(L)}}{\left|k^{(L)}\right|^{2}},~{}~{}k_{i}^{(L)}=2\frac{\sin\left(\pi\tilde{n}_{i}/N\right)}{dx},$

which are based on the nearest-neighbour spatial derivatives Figueroa  _et
al._ (2011) and fit well with our finite difference scheme in our code.

Figure 9: GW spectra in the FLRW scheme for Model III at $t\approx 172m^{-1}$
for three projectors: projector based on the central finite difference,
conventional projector based on k and projector based on the symmetric
nearest-neighbor derivatives.

## References

  * Bogolyubsky (1976) I. L. Bogolyubsky, JETP Lett. 24, 535 (1976).
  * Gleiser (1994) M. Gleiser, Phys. Rev. D 49, 2978 (1994), arXiv:hep-ph/9308279 .
  * Honda and Choptuik (2002) E. P. Honda and M. W. Choptuik, Phys. Rev. D 65, 084037 (2002), arXiv:hep-ph/0110065 .
  * Amin  _et al._ (2012) M. A. Amin, R. Easther, H. Finkel, R. Flauger, and M. P. Hertzberg, Phys. Rev. Lett. 108, 241302 (2012), arXiv:1106.3335 [astro-ph.CO] .
  * Amin  _et al._ (2010) M. A. Amin, R. Easther, and H. Finkel, JCAP 12, 001 (2010), arXiv:1009.2505 [astro-ph.CO] .
  * McDonald (2002) J. McDonald, Phys. Rev. D 66, 043525 (2002), arXiv:hep-ph/0105235 .
  * Broadhead and McDonald (2005) M. Broadhead and J. McDonald, Phys. Rev. D 72, 043519 (2005), arXiv:hep-ph/0503081 .
  * Copeland  _et al._ (2014) E. J. Copeland, P. M. Saffin, and S.-Y. Zhou, Phys. Rev. Lett. 113, 231603 (2014), arXiv:1409.3232 [hep-th] .
  * Xie  _et al._ (2021) Q.-X. Xie, P. M. Saffin, and S.-Y. Zhou, JHEP 07, 062 (2021), arXiv:2101.06988 [hep-th] .
  * Zhou  _et al._ (2013) S.-Y. Zhou, E. J. Copeland, R. Easther, H. Finkel, Z.-G. Mou, and P. M. Saffin, JHEP 10, 026 (2013), arXiv:1304.6094 [astro-ph.CO] .
  * Amin  _et al._ (2018) M. A. Amin, J. Braden, E. J. Copeland, J. T. Giblin, C. Solorio, Z. J. Weiner, and S.-Y. Zhou, Phys. Rev. D 98, 024040 (2018), arXiv:1803.08047 [astro-ph.CO] .
  * Liu  _et al._ (2018) J. Liu, Z.-K. Guo, R.-G. Cai, and G. Shiu, Phys. Rev. Lett. 120, 031301 (2018), arXiv:1707.09841 [astro-ph.CO] .
  * Antusch  _et al._ (2018a) S. Antusch, F. Cefala, S. Krippendorf, F. Muia, S. Orani, and F. Quevedo, JHEP 01, 083 (2018a), arXiv:1708.08922 [hep-th] .
  * Antusch  _et al._ (2018b) S. Antusch, F. Cefala, and S. Orani, JCAP 03, 032 (2018b), arXiv:1712.03231 [astro-ph.CO] .
  * Antusch  _et al._ (2017) S. Antusch, F. Cefala, and S. Orani, Phys. Rev. Lett. 118, 011303 (2017), [Erratum: Phys.Rev.Lett. 120, 219901 (2018)], arXiv:1607.01314 [astro-ph.CO] .
  * Zhou (2015) S.-Y. Zhou, JCAP 06, 033 (2015), arXiv:1501.01217 [astro-ph.CO] .
  * Sang and Huang (2019) Y. Sang and Q.-G. Huang, Phys. Rev. D 100, 063516 (2019), arXiv:1905.00371 [astro-ph.CO] .
  * Cotner  _et al._ (2018) E. Cotner, A. Kusenko, and V. Takhistov, Phys. Rev. D 98, 083513 (2018), arXiv:1801.03321 [astro-ph.CO] .
  * Cotner  _et al._ (2019) E. Cotner, A. Kusenko, M. Sasaki, and V. Takhistov, JCAP 10, 077 (2019), arXiv:1907.10613 [astro-ph.CO] .
  * Kou  _et al._ (2021) X.-X. Kou, C. Tian, and S.-Y. Zhou, Class. Quant. Grav. 38, 045005 (2021), arXiv:1912.09658 [gr-qc] .
  * Nazari  _et al._ (2021) Z. Nazari, M. Cicoli, K. Clough, and F. Muia, JCAP 05, 027 (2021), arXiv:2010.05933 [gr-qc] .
  * Hiramatsu  _et al._ (2021) T. Hiramatsu, E. I. Sfakianakis, and M. Yamaguchi, JHEP 03, 021 (2021), arXiv:2011.12201 [hep-ph] .
  * Copeland  _et al._ (1995) E. J. Copeland, M. Gleiser, and H. R. Muller, Phys. Rev. D 52, 1920 (1995), arXiv:hep-ph/9503217 .
  * Copeland  _et al._ (2002) E. J. Copeland, S. Pascoli, and A. Rajantie, Phys. Rev. D 65, 103517 (2002), arXiv:hep-ph/0202031 .
  * Farhi  _et al._ (2005) E. Farhi, N. Graham, V. Khemani, R. Markov, and R. Rosales, Phys. Rev. D 72, 101701 (2005), arXiv:hep-th/0505273 .
  * Gleiser (2007) M. Gleiser, Int. J. Mod. Phys. D 16, 219 (2007), arXiv:hep-th/0602187 .
  * Graham and Stamatopoulos (2006) N. Graham and N. Stamatopoulos, Phys. Lett. B 639, 541 (2006), arXiv:hep-th/0604134 .
  * Hindmarsh and Salmi (2006) M. Hindmarsh and P. Salmi, Phys. Rev. D 74, 105005 (2006), arXiv:hep-th/0606016 .
  * Fodor  _et al._ (2006) G. Fodor, P. Forgacs, P. Grandclement, and I. Racz, Phys. Rev. D 74, 124003 (2006), arXiv:hep-th/0609023 .
  * Saffin and Tranberg (2007) P. M. Saffin and A. Tranberg, JHEP 01, 030 (2007), arXiv:hep-th/0610191 .
  * Graham (2007) N. Graham, Phys. Rev. Lett. 98, 101801 (2007), [Erratum: Phys.Rev.Lett. 98, 189904 (2007)], arXiv:hep-th/0610267 .
  * Farhi  _et al._ (2008) E. Farhi, N. Graham, A. H. Guth, N. Iqbal, R. R. Rosales, and N. Stamatopoulos, Phys. Rev. D 77, 085019 (2008), arXiv:0712.3034 [hep-th] .
  * Hindmarsh and Salmi (2008) M. Hindmarsh and P. Salmi, Phys. Rev. D 77, 105025 (2008), arXiv:0712.0614 [hep-th] .
  * Bastero-Gil  _et al._ (2008) M. Bastero-Gil, M. Tristram, J. F. Macias-Perez, and D. Santos, Phys. Rev. D 77, 023520 (2008), arXiv:0709.3510 [astro-ph] .
  * Bastero-Gil  _et al._ (2010a) M. Bastero-Gil, J. Macias-Perez, and D. Santos, Phys. Rev. Lett. 105, 081301 (2010a), arXiv:1005.4054 [astro-ph.CO] .
  * Giblin and Tishue (2019) J. T. Giblin and A. J. Tishue, Phys. Rev. D 100, 063543 (2019), arXiv:1907.10601 [gr-qc] .
  * Matarrese  _et al._ (1998) S. Matarrese, S. Mollerach, and M. Bruni, Phys. Rev. D 58, 043504 (1998), arXiv:astro-ph/9707278 .
  * Hwang  _et al._ (2017) J.-C. Hwang, D. Jeong, and H. Noh, Astrophys. J. 842, 46 (2017), arXiv:1704.03500 [astro-ph.CO] .
  * Domènech and Sasaki (2018) G. Domènech and M. Sasaki, Phys. Rev. D 97, 023521 (2018), arXiv:1709.09804 [gr-qc] .
  * Gong (2019) J.-O. Gong, (2019), arXiv:1909.12708 [gr-qc] .
  * Tomikawa and Kobayashi (2020) K. Tomikawa and T. Kobayashi, Phys. Rev. D 101, 083529 (2020), arXiv:1910.01880 [gr-qc] .
  * De Luca  _et al._ (2020) V. De Luca, G. Franciolini, A. Kehagias, and A. Riotto, JCAP 03, 014 (2020), arXiv:1911.09689 [gr-qc] .
  * Inomata and Terada (2020) K. Inomata and T. Terada, Phys. Rev. D 101, 023523 (2020), arXiv:1912.00785 [gr-qc] .
  * Yuan  _et al._ (2020) C. Yuan, Z.-C. Chen, and Q.-G. Huang, Phys. Rev. D 101, 063018 (2020), arXiv:1912.00885 [astro-ph.CO] .
  * Lu  _et al._ (2020) Y. Lu, A. Ali, Y. Gong, J. Lin, and F. Zhang, Phys. Rev. D 102, 083503 (2020), arXiv:2006.03450 [gr-qc] .
  * Ali  _et al._ (2021) A. Ali, Y. Gong, and Y. Lu, Phys. Rev. D 103, 043516 (2021), arXiv:2009.11081 [gr-qc] .
  * Domènech and Sasaki (2021) G. Domènech and M. Sasaki, Phys. Rev. D 103, 063531 (2021), arXiv:2012.14016 [gr-qc] .
  * Gurian  _et al._ (2021) J. Gurian, D. Jeong, J.-c. Hwang, and H. Noh, Phys. Rev. D 104, 083534 (2021), arXiv:2104.03330 [astro-ph.CO] .
  * Silverstein and Westphal (2008) E. Silverstein and A. Westphal, Phys. Rev. D 78, 106003 (2008), arXiv:0803.3085 [hep-th] .
  * McAllister  _et al._ (2010) L. McAllister, E. Silverstein, and A. Westphal, Phys. Rev. D 82, 046003 (2010), arXiv:0808.0706 [hep-th] .
  * Dong  _et al._ (2011) X. Dong, B. Horn, E. Silverstein, and A. Westphal, Phys. Rev. D 84, 026011 (2011), arXiv:1011.4521 [hep-th] .
  * 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] .
  * Figueroa  _et al._ (2011) D. G. Figueroa, J. Garcia-Bellido, and A. Rajantie, JCAP 11, 015 (2011), arXiv:1110.0337 [astro-ph.CO] .
  * Isaacson (1968a) R. A. Isaacson, Phys. Rev. 166, 1263 (1968a).
  * Isaacson (1968b) R. A. Isaacson, Phys. Rev. 166, 1272 (1968b).
  * Bastero-Gil  _et al._ (2010b) M. Bastero-Gil, J. Macias-Perez, and D. Santos, Phys. Rev. Lett. 105, 081301 (2010b), arXiv:1005.4054 [astro-ph.CO] .
  * Ota  _et al._ (2021) A. Ota, H. J. Macpherson, and W. R. Coulton, (2021), arXiv:2111.09163 [gr-qc] .
  * Mertens  _et al._ (2016) J. B. Mertens, J. T. Giblin, and G. D. Starkman, Phys. Rev. D 93, 124059 (2016), arXiv:1511.01106 [gr-qc] .
  * Nakamura  _et al._ (1987) T. Nakamura, K. Oohara, and Y. Kojima, Prog. Theor. Phys. Suppl. 90, 1 (1987).
  * Shibata and Nakamura (1995) M. Shibata and T. Nakamura, Phys. Rev. D 52, 5428 (1995).
  * Baumgarte and Shapiro (1998) T. W. Baumgarte and S. L. Shapiro, Phys. Rev. D 59, 024007 (1998), arXiv:gr-qc/9810065 .
  * Misner  _et al._ (1973) C. Misner, K. Thorne, K. Thorne, J. Wheeler, W. Freeman, and Company,  _Gravitation_, Gravitation (W. H. Freeman, 1973).
  * Chen and Zhu (2011) X.-S. Chen and B.-C. Zhu, Phys. Rev. D 83, 061501 (2011), arXiv:1006.3927 [gr-qc] .
  * Domènech (2021) G. Domènech, Universe 7, 398 (2021), arXiv:2109.01398 [gr-qc] .