# Gravitational wave induced baryon acoustic oscillations

C<EMAIL_ADDRESS>Max-Planck-Institut für Kernphysik,
Saupfercheckweg 1, 69117 Heidelberg, Germany,  S.
<EMAIL_ADDRESS>AHEP Group, Institut de Física Corpuscular –
C.S.I.C./Universitat de València, Parc Científic de Paterna.
C/ Catedrático José Beltrán, 2 E-46980 Paterna (Valencia), SPAIN,  M.
<EMAIL_ADDRESS>Max-Planck-Institut für Kernphysik,
Saupfercheckweg 1, 69117 Heidelberg, Germany,  B. M.
<EMAIL_ADDRESS>Zentrum für Astronomie der
Universität Heidelberg, Astronomisches Rechen-Institut, Philosophenweg 12,
69120 Heidelberg, Germany,
M<EMAIL_ADDRESS>Institut für Theoretische Physik,
Heidelberg University, Philosophenweg 16, 69120 Heidelberg, Germany

We study the impact of gravitational waves originating from a first order
phase transition on structure formation. To do so, we perform a second order
perturbation analysis in the $1+3$ covariant framework and derive a wave
equation in which second order, adiabatic density perturbations of the photon-
baryon fluid are sourced by the gravitational wave energy density during
radiation domination and on sub-horizon scales. The scale on which such waves
affect the energy density perturbation spectrum is found to be proportional to
the horizon size at the time of the phase transition times its inverse
duration. Consequently, structure of the size of galaxies and bigger can only
be affected in this way by relatively late phase transitions at $\geq
10^{6}\,\text{s}$. Using cosmic variance as a bound we derive limits on the
strength $\alpha$ and the relative duration $(\beta/H_{*})^{-1}$ of phase
transitions as functions of the time of their occurrence which results in a
new exclusion region for the energy density in gravitational waves today. We
find that the cosmic variance bound forbids only relative long lasting phase
transitions, e.g. $\beta/H_{*}\lesssim 6.8$ for $t_{*}\approx 5\times
10^{11}\,\text{s}$, which exhibit a substantial amount of supercooling
$\alpha>20$ to affect the matter power spectrum.

Keywords: First order phase transitions, gravitational waves, structure
formation, matter power spectrum, 1+3 gauge invariant theory.

##  1 Introduction

With the first ever-measurement of a gravitational wave (GW) signal in 2016
from a black hole binary merger by the LIGO collaboration [1], a new window of
probing the universe has been opened. While this technique probes so far
mostly astrophysical processes, future experiments like the space
interferometer LISA [2] have the potential to explore also cosmological
sources like first order phase transitions (FPT) in the early universe. In
particle physics these phase transitions occur when the dropping temperature
of the universe causes the vacuum expectation value (VEV) of a field to change
discontinuously. If the field is hindered for a while from adapting to the new
VEV by a barrier in its potential then bubbles enclosing the new VEV form,
expand and eventually fill the universe with the new VEV. While such FPTs can
produce GWs via the dynamics of the bubbles like collisions, soundwave
formation and magneto-hydrodynamical effects, second order phase transitions
and cross overs are not expected to produce substantial amounts of GWs,
essentially because they lack the mechanism of vacuum bubble formation. The
latter applies to the standard model of particles physics (SM), well described
by the symmetry group $SU(3)_{QCD}\times SU(2)_{L}\times U(1)_{Y}$. It
undergoes a cross over phase transition during the electro-weak symmetry
breaking $SU(2)_{L}\times U(1)_{Y}\to U(1)_{QED}$ when the Higgs boson
acquires a non-zero VEV [3, 4] and hence no GW signal is expected. The SM has,
however, various problems, motivating for new physics beyond the SM (BSM).
Many alternative models which incorporate new symmetries and particles allow
for FPTs. The observation of GWs has therefore triggered many studies of FPTs
in BSM models [5, 6, 7, 8, 9, 10, 11] where often GWs are expected to be seen
in future GW experiments. For reviews see e.g. [12, 13, 14, 15, 16]. Future GW
experiments can therefore valuably constrain BSM physics.

However, adding a FPT to the history of the universe might also affect other
cosmological processes such as formation of structure. This potential
consequence is studied in this work. In the standard model of cosmology,
linear density perturbations develop from inflation and seed over- and under
densities in the various fluid components of the early cosmological medium.
They propagate through the universe and undergo, depending on their scale,
various changes caused by physical processes like the decoupling of a fluid
component until they eventually form the structure we observe today. FPTs and
the emerging GWs might influence this evolution depending on strength and
duration. GWs are linear tensor perturbations of the metric sourced by an
anisotropic stress distortion in the fluid while density perturbations of the
fluid are scalar perturbations of the metric. At linear order in perturbation
theory, they do not couple, but they can interact at second order and source
second order density perturbations. Hence, we have to perform a second order
expansion in order to capture effects that strong GW events may have.
Typically, phase transitions are expected to occur while the universe is
dominated by radiation and on sub-horizon scales. Consequently, potential
effects on density perturbations are tied to the scale and thus the time of
the transition. We shall work within the $1+3$ covariant approach to gravity
[17, 18, 19, 20, 21, 22] in which an exact, non linear equation of density
perturbations is given.

Imprints of phase transitions in the matter power spectrum666Linear matter
power spectrum and measurements can be found in [23]. have been of interest in
the past, [24, 25]. In contrast to our work these papers focus on the QCD
phase transition during which they predicted a significant drop in the sound
speed. This in turn affects the preexisting linear density perturbations and
induces large peaks in the Harrison-Zel’dovich spectrum. Similar effects could
happen in BSM transitions if the sound speed drops significantly which could
be possible for theories with massive fermions or weakly interacting scalars
but is not expected e.g. in simple scalar extensions of the SM [26].

This work is organized as follows. In Sec. 2 we first give an introduction to
the $1+3$ covariant formulation of gravity and show how to derive the
dynamical equations of the familiar linear density perturbations as well as
for the shear perturbations (GWs). Then we move on to investigate second order
density perturbations and their coupling to GWs in Sec. 3. In Sec. 4 we then
summarize the physics of GWs from FPTs and present our results in Sec. 5.
Subsequently we discuss in which way and under which conditions the GWs from
FPTs do or do not affect the matter power spectrum, but also debate the
limitations of our approach. Finally we conclude and give an outlook in Sec.
6.

##  2 The 1+3-covariant formulation

Quantities in perturbation theory are in general gauge dependent, i.e. they
change under infinitesimal coordinate transformations
$\tilde{x}^{\mu}=x^{\mu}+\epsilon\,\xi^{\mu}$, where $\epsilon\ll 1$ is a
small parameter and $\xi^{\mu}\in\mathbb{R}^{4}$ is some vector field. Under a
linear perturbation an arbitrary tensor field $S$ is split into its zeroth-
order part777Subscripts in parentheses denote the perturbative order.
$S^{(0)}$, also referred to as background value, and its first-order part
$\delta S\equiv\epsilon\,S^{(1)}$, i.e. $S=S^{(0)}+\epsilon\,S^{(1)}$. Under a
gauge transformation the latter perturbative component is not simply mapped to
itself but rather receives an additional term dependent of the gauge vector
$\xi^{\mu}$ according to

$\displaystyle S^{(1)}\rightarrow S^{(1)}+\mathcal{L}_{\xi}S^{(0)},$ (1)

where this additional term is the Lie-derivative $\mathcal{L}_{\xi}$ along
$\xi^{\mu}$ of the background term $S^{(0)}$ [27, 28, p.59]. In order to make
universally valid predictions for a perturbative physical model we need to
introduce gauge invariant quantities.

A tensor field is called gauge invariant to first order if for any vector
field $\xi^{\mu}$ the Lie-derivative vanishes $\mathcal{L}_{\xi}S^{(0)}$.
Based on the gauge transformation rule in Eq. (1) the Stewart & Walker lemma
[29] states that a tensor is gauge invariant888Unless stated otherwise, we
mean by gauge invariant always gauge invariant to linear order. if and only if
it either vanishes in the background, is a constant scalar on the background
or can be written as a sum of products of Kronecker-deltas with constant
coefficients [27].

In the spirit of this lemma, Ellis, Bruni and co-authors developed the so
called $1+3$ covariant formulation of gravity [21, 22] based on earlier papers
by Heckmann and Schücking [17], Raychaudhuri [18], Ehlers [19] and Hawking
[20]. In this section we will follow closely ref. [30]. The advantage of this
approach resides in the simple geometric meaning of the central variables and
their gauge invariance which is due to the fact that they vanish in a
spatially homogeneous background, for example in the background of a
Friedmann-Lemaître-Robertson-Walker (FLRW) metric. These variables are
constructed by decomposing the spacetime into the direction of the four-
velocity of a comoving observer that follows the fluid flow lines $x^{a}$ and
the projection tensor into the instantaneous rest space of $u^{a}$,

$\displaystyle u^{a}=\frac{\text{d}x^{a}}{\text{d}\tau}\quad\text{and}\quad
h_{ab}:=g_{ab}+u_{a}u_{b},$ (2)

with the proper time $\tau$, $u^{a}u_{a}=-1$ and $g_{ab}$ being the metric
tensor with signature $(-+++)$. We follow the convention of the literature and
use Latin indices for four-vectors $a,b,c,\dots=0,1,2,3$ and
$\alpha,\beta,\gamma\dots=1,2,3$ for spacelike three-vectors.

The two tensors are perpendicular projectors

$\displaystyle h_{ab}u^{b}=g_{ab}u^{b}+u_{a}u_{b}u^{b}=u_{a}-u_{a}=0$ (3)

that project a spacetime quantity onto the flow lines or in the orthogonal
direction which enables a unique splitting into irreducible timelike and
spacelike components (establishing the name $1+3$).

Exemplarily the time- and space derivative of a general tensor
$S_{ab\dots}^{\;cd\dots}$ is obtained by projecting the covariant derivative
$\nabla_{a}$:

$\displaystyle\dot{S}_{ab\dots}{}^{cd\dots}:=u^{e}\nabla_{e}S_{ab\dots}{}^{cd\dots}\quad\text{and}\quad{\text{D}}_{e}S_{ab\dots}{}^{cd\dots}:=h_{e}{}^{s}h_{a}{}^{f}h_{q}{}^{c}\cdots\nabla_{s}S_{f\cdots}{}^{q\cdots}.$
(4)

The next step is to describe the kinematics of an observer in this framework
under the influence of gravity and matter represented by the energy momentum
tensor $T_{ab}$. We will set the speed of light $c=1$ and the gravitational
coupling $\kappa:=8\pi G=1$, where $G$ is the gravitational constant.

###  2.1 Kinematic variables

In the $1+3$-covariant approach to gravity, the kinematic quantities that
determine the motion of a test particle are the tracefree shear tensor
$\sigma_{ab}:={\text{D}}_{\langle b}u_{a\rangle}$, the antisymmetric (hence
tracefree) vorticity tensor $\omega_{ab}:={\text{D}}_{[b}u_{a]}$, the volume
expansion scalar $\Theta:={\text{D}}^{a}u_{a}$ and the four-acceleration
$A_{a}=u^{b}\nabla_{b}u_{a}$, which emerge from the irreducible decomposition
of the covariant derivative of the four-velocity

$\displaystyle\nabla_{b}u_{a}=\sigma_{ab}+\omega_{ab}+\frac{1}{3}\Theta
h_{ab}-A_{a}u_{b}.$ (5)

The here used brackets are defined as

$\displaystyle S_{(ab)}:=\frac{1}{2}\left(S_{ab}+S_{ba}\right),\quad
S_{[ab]}:=\frac{1}{2}\left(S_{ab}-S_{ba}\right),$ (6) $\displaystyle
S_{\langle
ab\rangle}:=h_{a}^{\;c}h_{b}^{\;d}S_{cd}-\frac{1}{3}h^{cd}S_{cd}h_{ab},\quad
V_{\langle a\rangle}:=h_{a}^{\;b}V_{b}.$ (7)

Useful identities for these objects are collected in Appendix A. In an FLRW
universe at zeroth-order the shear, the vorticity and the acceleration vanish
and hence the Stewart & Walker Lemma provides gauge invariant quantities [31].
Moreover in a spatially homogeneous model like the FLRW metric the background
value of the scalar $\Theta^{(0)}(t)$ depends solely on time and thus the
spatial derivative ${\text{D}}_{a}\Theta$ is equally gauge-invariant.

###  2.2 Gravity

In general relativity gravitation arises from intrinsic properties of the
spacetime manifold and matter. Einstein’s field equations formulate this
relation by

$\displaystyle R_{ab}-\frac{1}{2}Rg_{ab}=\kappa T_{ab}-\Lambda g_{ab},$ (8)

where $R_{ab}$ and $R$ are the Ricci tensor and scalar, respectively, and
$\Lambda$ is the cosmological constant. The Ricci tensor is the contraction of
the Riemann tensor $R_{abcd}$ which encodes the curvature of the spacetime
manifold. The latter can be split into two parts

$\displaystyle
R_{abcd}=C_{abcd}+\frac{1}{2}(g_{ac}R_{bd}+g_{bd}R_{ac}-g_{bc}R_{ad}-g_{ad}R_{bc})-\frac{1}{6}R(g_{ac}g_{bd}-g_{ad}g_{bc}).$
(9)

While Ricci tensor $R_{ab}$ and scalar $R$ express volume changes due to a
local matter source and hence reflect the local part of the gravitational
field, the Weyl tensor $C_{abcd}$999The Weyl tensor shares all symmetries with
the Riemann tensor $R_{abcd}=R_{cdcb}$, $R_{abcd}=R_{[ab][cd]}$ and
$R_{a[bcd]}=0$ and is, by construction, additionally tracefree. contains
information about the propagating degrees of freedom. Using the four-velocity
vector, $C_{abcd}$ can be decomposed further into the so called electric and
magnetic parts [32, 33]

$\displaystyle E_{ab}=C_{acbd}u^{c}u^{d}\quad\text{and}\quad
H_{ab}=\frac{1}{2}\epsilon_{a}^{cd}C_{cdbe}u^{e},\>\text{respectively}.$ (10)

Both tensors are symmetric, tracefree and gauge invariant due to
$C_{abcd}^{(0)}=0$ in the FLRW background. As we shall see, the propagation of
GWs is mainly governed by the magnetic part $H_{ab}$ while the electric part
$E_{ab}$ is closely related to tidal forces.

Having discussed long range gravitational effects, let us focus on local
gravity which is expressed by the Ricci tensor and the energy-momentum tensor.
For a general fluid the energy-momentum tensor decomposes with respect to the
fundamental timelike velocity field into

$\displaystyle T_{ab}=\rho u_{a}u_{b}+2u_{(a}q_{b)}+ph_{ab}+\pi_{ab},$ (11)

where $\rho:=T^{ab}u_{a}u_{b}$ is the energy density,
$q_{b}:=h_{a}^{\;b}T_{bc}u^{c}$ is the energy current density,
$p:=T_{ab}h^{ab}/3$ is the pressure and $\pi_{ab}:=T_{\langle ab\rangle}$ the
trace-free anisotropic stress. While both the anisotropic stress and the
energy current density again vanish in a FLRW universe and are thus gauge
invariant, the pressure and the energy density depend only on time and hence
their spatial gradients are also gauge invariant to first order.

Rewriting Einstein’s equation as
$R_{ab}=\kappa(T_{ab}-\frac{1}{2}T_{ab})+\Lambda g_{ab}$ leads to three
equations that relate the Ricci tensor and the matter-fields [30],

$\displaystyle R_{ab}u^{a}u^{b}=\kappa{\frac{1}{2}}\,(\rho+3p)-\Lambda,$ (12)
$\displaystyle h_{a}{}^{b}R_{bc}u^{c}=-\kappa q_{a}\quad\text{and}$ (13)
$\displaystyle
h_{a}{}^{c}h_{b}{}^{d}R_{cd}=\kappa{\frac{1}{2}}\,(\rho-p)h_{ab}+\kappa\pi_{ab}+\Lambda
h_{ab}.$ (14)

###  2.3 Equations of motion

As discussed in the previous sections, the $1+3$ covariant approach identifies
gauge invariant components of the energy-momentum tensor, the Riemann tensor
and the four velocity gradient with a clear geometrical and physical meaning.
The equations of motion for these variables are inferred from the Bianchi and
Ricci identities and are accompanied by constraint equations. The equations
quoted in this section have been derived in [34, 35], for details see also the
review [30]. Using Eq. (14) and the Bianchi identities for the Weyl tensor

$\displaystyle\nabla^{d}C_{abcd}=\nabla_{[b}R_{a]c}+\frac{1}{6}g_{c[b}\nabla_{a]}R$
(15)

one finds the non-linear propagation equations of the electric and magnetic
components of the Weyl-tensor,

$\displaystyle\dot{E}_{\langle ab\rangle}$ $\displaystyle=$
$\displaystyle-\Theta
E_{ab}-{\frac{1}{2}}\kappa\,(\rho+p)\sigma_{ab}+\,\text{curl}\,H_{ab}-{\frac{1}{2}}\kappa\,\dot{\pi}_{ab}-{\frac{1}{6}}\kappa\,\Theta\pi_{ab}-{\frac{1}{2}}\kappa\,{\text{D}}_{\langle
a}q_{b\rangle}-\kappa A_{\langle a}q_{b\rangle}$ (16)
$\displaystyle+3\sigma_{\langle a}{}^{c}\left(E_{b\rangle
c}-{\frac{1}{6}}\kappa\,\pi_{b\rangle c}\right)+\varepsilon_{cd\langle
a}\left[2A^{c}H_{b\rangle}{}^{d}-\omega^{c}\left(E_{b\rangle}{}^{d}+{\frac{1}{2}}\kappa\,\pi_{b\rangle}{}^{d}\right)\right],$
$\displaystyle\dot{H}_{\langle ab\rangle}$ $\displaystyle=$
$\displaystyle-\Theta
H_{ab}-\,\text{curl}\,E_{ab}+{\frac{1}{2}}\kappa\,\,\text{curl}\,\pi_{ab}+3\sigma_{\langle
a}{}^{c}H_{b\rangle c}-{\frac{3}{2}}\kappa\,\omega_{\langle a}q_{b\rangle}$
(17) $\displaystyle-\varepsilon_{cd\langle
a}\left(2A^{c}E_{b\rangle}{}^{d}-{\frac{1}{2}}\kappa\,\sigma^{c}{}_{b\rangle}q^{d}+\omega^{c}H_{b\rangle}{}^{d}\right)\,.$

while the spacelike constraints become

$\displaystyle{\text{D}}^{b}E_{ab}=\kappa\left[{\frac{1}{3}}\,{\text{D}}_{a}\rho-{\frac{1}{2}}{{\text{D}}}^{b}\pi_{ab}-{\frac{1}{3}}\,\Theta
q_{a}+{\frac{1}{2}}\,\sigma_{ab}q^{b}\right]-3H_{ab}\omega^{b}+\varepsilon_{abc}\left(\sigma^{b}{}_{d}H^{cd}-{\frac{3}{2}}\kappa\,\omega^{b}q^{c}\right)$
(18)

and

$\displaystyle{\text{D}}^{b}H_{ab}=\kappa(\rho+p)\omega_{a}-{\frac{1}{2}}\kappa\,\,\text{curl}\,q_{a}+3E_{ab}\omega^{b}-{\frac{1}{2}}\kappa\,\pi_{ab}\omega^{b}-\varepsilon_{abc}\sigma^{b}{}_{d}\left(E^{cd}+{\frac{1}{2}}\kappa\,\pi^{cd}\right)\,.$
(19)

Here, the vorticity vector $\omega_{a}:=\epsilon_{abc}\omega^{bc}/2$ has been
introduced together with the projection $\epsilon_{abc}:=\eta_{abcd}u^{d}$ of
the totally antisymmetric tensor $\eta_{abcd}$.

From the Bianchi identity expressing the conservation of energy and momentum
$\nabla^{a}T_{ab}=0$ we find the propagation equations for the energy density
and the energy flux,

$\displaystyle\dot{\rho}=-\Theta(\rho+p)-{\text{D}}^{a}q_{a}-2A^{a}q_{a}-\sigma^{ab}\pi_{ab},$
(20) $\displaystyle\dot{q}_{\langle
a\rangle}=-{\text{D}}_{a}p-(\rho+p)A_{a}-{\frac{4}{3}}\,\Theta
q_{a}-(\sigma_{ab}+\omega_{ab})q^{b}-{{\text{D}}}^{b}\pi_{ab}-\pi_{ab}A^{b}.$
(21)

Finally the Ricci-identities $2\nabla_{[a}\nabla_{b]}u_{c}=R_{abcd}u^{d}$ give
the Raychaudhuri equation for the volume expansion and the propagation
equations for the shear and the vorticity

$\displaystyle\dot{\Theta}=-{\frac{1}{3}}\,\Theta^{2}-{\frac{1}{2}}\kappa\,(\rho+3p)-2(\sigma^{2}-\omega^{2})+{\text{D}}^{a}A_{a}+A_{a}A^{a}+\Lambda,$
(22) $\displaystyle\dot{\sigma}_{\langle
ab\rangle}=-{\frac{2}{3}}\,\Theta\sigma_{ab}-\sigma_{c\langle
a}\sigma^{c}{}_{b\rangle}-\omega_{\langle
a}\omega_{b\rangle}+{\text{D}}_{\langle a}A_{b\rangle}+A_{\langle
a}A_{b\rangle}-E_{ab}+{\frac{1}{2}}\kappa\,\pi_{ab},$ (23)
$\displaystyle\dot{\omega}_{\langle
a\rangle}=-{\frac{2}{3}}\,\Theta\omega_{a}-{\frac{1}{2}}\,\,\text{curl}\,A_{a}+\sigma_{ab}\omega^{b}.$
(24)

These identities also imply the following constraints for the shear, the
vorticity and the magnetic component of the Weyl tensor

$\displaystyle{\text{D}}^{b}\sigma_{ab}={\frac{2}{3}}\,{\text{D}}_{a}\Theta+\,\text{curl}\,\omega_{a}+2\varepsilon_{abc}A^{b}\omega^{c}-\kappa
q_{a},\quad{\text{D}}^{a}\omega_{a}=A_{a}\omega^{a},$ (25) $\displaystyle
H_{ab}=\,\text{curl}\,\sigma_{ab}+{\text{D}}_{\langle
a}\omega_{b\rangle}+2A_{\langle a}\omega_{b\rangle}\,.$ (26)

These equations allow us to find the behavior of density perturbations and GWs
on a FLRW background. For a detailed derivation of these equations see [30].
In table 1 we have summarized the central quantities of the $1+3$ approach,
their interpretation and their gauge properties.

Variable | Symbol |  Perturbative Expansion $S=S^{(0)}+\epsilon\,S^{(1)}$ | First order GI
---|---|---|---
Energy density | $\rho$ | $\rho(t)+\rho(\textbf{x},t)$ | ${\text{D}}_{a}\rho(\textbf{x},t)$
Pressure | $p$ | $p(t)+p(\textbf{x},t)$ | ${\text{D}}_{a}p(\textbf{x},t)$
Anisotropic stress | $\pi$ | $0+\pi(\textbf{x},t)$ | $\pi(\textbf{x},t)$
Energy density current | $q$ | $0+q(\textbf{x},t)$ | $q(\textbf{x},t)$
Volume expansion | $\Theta$ | $\Theta(t)+\Theta(\textbf{x},t)$ | ${\text{D}}_{a}\Theta(\textbf{x},t)$
Shear | $\sigma$ | $0+\sigma(\textbf{x},t)$ | $\sigma(\textbf{x},t)$
Vorticity | $\omega$ | $0+\omega(\textbf{x},t)$ | $\omega(\textbf{x},t)$
Acceleration | $A$ | $0+A(\textbf{x},t)$ | $A(\textbf{x},t)$
Long range grav. field (Weyl tensor) | $C$ | $0+C_{abcd}(\textbf{x},t)$ | $C_{abcd}(\textbf{x},t)$
Table 1: Perturbative expansion of the central quantities in the $1+3$
covariant theory in a FLRW background and their first order gauge invariant
version. In the third column the first term refers to the zeroth-order
component and the second term to $\delta S\equiv\epsilon S^{(1)}$. Since for
most quantities the zeroth-order part is zero and hence $S=\delta S$ we often
neither use the superposed index nor the $\delta$-symbol for their first-order
term.

###  2.4 Linear density perturbations

Spatial inhomogeneities in the matter density are described by the spatial
comoving fractional gradient and the comoving expansion gradient [21]

$\displaystyle\Delta_{a}:=\frac{a}{\rho}{\text{D}}_{a}\rho,$ (27)
$\displaystyle Z_{a}:=a{\text{D}}_{a}\Theta,$ (28)

which are both orthogonal to the fluid flow. In a spatially homogeneous
background they are gauge invariant because $\rho$ and $\Theta$ depend only on
time such that the spatial gradient
${\text{D}}_{a}\rho|_{\text{background}}={\text{D}}_{a}\rho(t)=0$ vanishes in
the background. The time- and space dependent variations of over- and under
densities which are expressed by the orthogonal projected divergence of the
comoving fractional gradient $a{\text{D}}^{a}\Delta_{a}=:\Delta$ are closely
related to the Laplacian of the density contrast $\delta:=\delta\rho/\rho$.
However, besides the usual over- and under densities also distortions
$\Delta_{\langle ab\rangle}$ and vorticity $\Delta_{[ab]}$ can be introduced
by the splitting

$\displaystyle a{\text{D}}_{b}\Delta_{a}=\frac{1}{3}\Delta
h_{ab}+\Delta_{\langle ab\rangle}+\Delta_{[ab]}.$ (29)

Taking into account the equations and constraints from the Bianchi identities,
spatial inhomogeneities evolve according to the full non-linear equations [30]
as

$\displaystyle\dot{\Delta}_{\langle a\rangle}$
$\displaystyle={\frac{p}{\rho}}\,\Theta\Delta_{a}-\left(1+{\frac{p}{\rho}}\right){Z}_{a}+{a\frac{\Theta}{\rho}}\left(\dot{q}_{\langle
a\rangle}+{\frac{4}{3}}\,\Theta
q_{a}\right)-{\frac{a}{\rho}}\,{\text{D}}_{a}{\text{D}}^{b}q_{b}+{a\frac{\Theta}{\rho}}\,{\text{D}}^{b}\pi_{ab}$
$\displaystyle-\left(\sigma^{b}{}_{a}+\omega^{b}{}_{a}\right)\Delta_{b}-{\frac{a}{\rho}}\,{\text{D}}_{a}\left(2A^{b}q_{b}+\sigma^{bc}\pi_{bc}\right)+{a\frac{\Theta}{\rho}}\left(\sigma_{ab}+\omega_{ab}\right)q^{b}+{a\frac{\Theta}{\rho}}\,\pi_{ab}A^{b}$
$\displaystyle+{\frac{1}{\rho}}\left({\text{D}}^{b}q_{b}+2A^{b}q_{b}+\sigma^{bc}\pi_{bc}\right)\left(\Delta_{a}-aA_{a}\right)\,,$
(30)

and

$\displaystyle\dot{Z}_{\langle a\rangle}$
$\displaystyle=-{\frac{2}{3}}\,\Theta
Z_{a}-{\frac{1}{2}}\,\kappa\rho\Delta_{a}-{\frac{3}{2}}\,\kappa
a{\text{D}}_{a}p-a\left[{\frac{1}{3}}\,\Theta^{2}+{\frac{1}{2}}\,\kappa(\rho+3p)-\Lambda\right]A_{a}+a{\text{D}}_{a}{\text{D}}^{b}A_{b}$
$\displaystyle-\left(\sigma^{b}{}_{a}+\omega^{b}{}_{a}\right)Z_{b}-2a{\text{D}}_{a}\left(\sigma^{2}-\omega^{2}\right)+2aA^{b}{\text{D}}_{a}A_{b}$
$\displaystyle-a\left[2\left(\sigma^{2}-\omega^{2}\right)-{\text{D}}^{b}A_{b}-A^{b}A_{b}\right]A_{a}\,.$
(31)

Here $\sigma^{2}:=\frac{1}{2}\sigma_{ab}\sigma^{ab}$ and
$\omega^{2}:=\frac{1}{2}\omega_{ab}\omega^{ab}$.

As the full non-linear equations are too complex to solve we seek to perturb
the equations to first order. Therefore, we have to choose a background model,
the FLRW metric, which reads in spherical coordinates

$\displaystyle\text{d}s^{2}=-\text{d}t^{2}+a^{2}(t)\left[\frac{\text{d}^{2}r}{1-Kr^{2}}+r^{2}\text{d}\Omega(\phi,\theta)\right].$
(32)

The curvature parameter $K$ can be $-1,0,+1$. In this background the volume
expansion is related to the Hubble parameter $H(t):=\dot{a}/a$ by
$\Theta^{(0)}(t)=3H(t)$ and thus Raychaudhuri’s equation, the continuity
equation and the Friedmann equation read

$\displaystyle\dot{H}=-H^{2}-\frac{\kappa}{6}(\rho+3p)+\frac{1}{3}\Lambda,\quad\dot{\rho}=-3H(\rho+p),\,\text{and}$
(33) $\displaystyle
H^{2}=\frac{\kappa}{3}\rho-\frac{K}{a^{2}}+\frac{1}{3}\Lambda.$ (34)

The evolution equation for linear density perturbations in a barotropic
perfect fluid $p=\omega\rho$ in a FLRW universe $\Theta=3H(t)$ with zero
vorticity $\omega_{ab}=0$ is then obtained from these equations setting the
energy current density and the anisotropic stress to zero. This leads to [22]

$\displaystyle\ddot{\Delta}=$
$\displaystyle-2\left(1-3\omega+\frac{3}{2}c_{s}^{2}\right)H\dot{\Delta}$
$\displaystyle+\kappa\left[\left(\frac{1}{2}+4\omega-3c_{s}^{2}-\frac{3}{2}\omega^{2}\right)\rho+(5\omega-3c_{s}^{2})\Lambda-\frac{12(\omega-
c_{s}^{2})K}{a^{2}}\right]\Delta$
$\displaystyle+c_{s}^{2}{\text{D}}^{2}\Delta.$ (35)

Since for first order gauge invariant variables the zeroth order is zero we
omit their perturbative labels and since appearing $\rho$’s, $p$’s and
$\Theta$’s always occur together with a gauge invariant variable they must be
of zeroth order such that we can also omit their superscripts.

It is useful to convert this equation into $k$-space by expanding $\Delta$ in
scalar harmonics $\mathcal{Q}_{k}$ such that
$\Delta=\int_{k}\Delta_{k}\mathcal{Q}_{k}$. The latter ones have the
properties $\dot{{\mathcal{Q}}_{k}}=0$ and
${\text{D}}^{2}{\mathcal{Q}}_{k}=-\frac{k^{2}}{a^{2}}{\mathcal{Q}}_{k}$ (see
appendix B for more information). In a spatially flat spacetime $K=0$ the
scalar harmonics are plane waves. For a radiation dominated, flat universe we
have $\omega=c_{s}^{2}=1/3$, $H=1/(2t)$, $a(t)=a_{0}\sqrt{t/t_{0}}$,
$\kappa\rho\equiv\rho^{(0)}=3/(4t^{2})$ and $K=0$ such that in a comoving
frame the dynamical equation for the $k$-th density perturbation mode yields

$\displaystyle\frac{\text{d}^{2}\Delta_{k}}{\text{d}t^{2}}+\frac{1}{2t}\frac{\text{d}\Delta_{k}}{\text{d}t}-\frac{1}{2t^{2}}\left[1-\frac{1}{6}\left(\frac{k}{a(t)H(t)}\right)^{2}\right]\Delta_{k}=0.$
(36)

which during radiation domination yields an oscillatory solution on sub-
horizon scales and a linearly growing solution on super-horizon scales. During
matter domination all modes grow with $\sim t^{2/3}$.

###  2.5 Linear metric perturbations: Gravitational waves

In the $1+3$ covariant approach long range gravity effects are incorporated by
the Weyl tensor and hence GWs are monitored by means of the transverse and
tracefree components of its electric and magnetic parts. Linearizing the
propagation equations Eq. (16), Eq. (17) and the constraints Eq. (18), Eq.
(19) these equations read [36]

$\displaystyle\dot{E}_{ab}=-\Theta
E_{ab}+\,\text{curl}\,H_{ab}-\frac{1}{2}\kappa\left[(\rho+p)\sigma_{ab}-{\text{D}}_{\langle
a}q_{b\rangle}+\dot{\pi}_{ab}+\frac{1}{3}\Theta\pi_{ab}\right],$ (37)
$\displaystyle\dot{H}_{ab}=-\Theta
H_{ab}-\,\text{curl}\,E_{ab}-\frac{1}{2}\kappa\pi_{ab},$ (38)
$\displaystyle{\text{D}}^{b}E_{ab}=\kappa\left(\frac{1}{3}\Theta
q_{b}+\frac{1}{3}{\text{D}}_{a}\rho+\frac{1}{2}{\text{D}}^{a}\pi_{ab}\right)\quad\text{and}\quad{\text{D}}^{b}H_{ab}=\frac{1}{2}\kappa\left[2(\rho+p)\omega_{a}+\,\text{curl}\,q_{b}\right].$
(39)

Hence in the absence of vorticity, the electric and magnetic parts of the Weyl
tensor that are not sourced by density gradients
(${\text{D}}_{a}\rho={\text{D}}_{a}\Theta=0$) are transverse tensors for a
perfect fluid ($\pi_{ab}=q_{a}=0$) on a FLRW background ($\Theta=3H(t)$) and
due to Eq. (25) this is also true for the shear

$\displaystyle{\text{D}}^{a}E_{ab}=0,\quad{\text{D}}^{a}H_{ab}=0\quad\text{and}\quad{\text{D}}^{a}\sigma_{ab}=0.$
(40)

Using the linearized equations of motion for the shear Eq. (23) and the Weyl
tensors, the latter ones can be eliminated from the discussion, see [37], to
give the propagation equation

$\displaystyle\ddot{\sigma}_{ab}+5H(t)\dot{\sigma}_{ab}+\frac{1}{2}\kappa\rho(1-3\omega)\sigma_{ab}-{\text{D}}^{2}\sigma_{ab}=0,$
(41)

in the absence of curvature $K=0$.

###  2.6 Connection to Bardeen-formalism and Newtonian theory

The standard formalism frequently used for studying structure formation is
based on the approach introduced by Bardeen [38]. In reference [27] Bruni et
al. gave the transformation equations between the $1+3$-formalism presented
here, and the $3+1$-slicing used by Bardeen. For later use and to connect to
the more common formalism of Bardeen let us briefly repeat here the
transformation rules. Primes denote in the following the conformal derivative
$\text{d}\eta=\text{d}t/a$. We introduce the perturbed metric parametrized as

$\displaystyle\text{d}s^{2}=a^{2}(\eta)\left\\{-(1+2A)\text{d}\eta^{2}-2B_{\alpha}\text{d}\eta\text{d}x^{\alpha}+[(1-2D)\delta_{\alpha\beta}+2E_{\alpha\beta}]\text{d}x^{\alpha}\text{d}x^{\beta})\right\\},$
(42) $\displaystyle\text{or}\>g_{ab}=\eta_{ab}+\delta
g_{ab}=a^{2}\left[\eta_{ab}+\begin{pmatrix}-2A&-B_{\alpha}\\\
-B_{\alpha}&-2D\delta_{\alpha\beta}+2E_{\alpha\beta},\end{pmatrix}\right],$

where $D=1/6\cdot\delta g_{a}^{a}$ and
$\delta^{\alpha\beta}E_{\alpha\beta}=0$. The vector B is commonly split into a
curl free, longitudinal part $\nabla\times\textbf{B}^{||}=0$ and a source
free, transverse $\nabla\cdot\textbf{B}^{\bot}=0$ such that
$\textbf{B}=\textbf{B}^{||}+\textbf{B}^{\bot}$. While the first part can be
written in terms of a scalar potential $B$ the latter one originates from a
vector potential. Similarly the tensor $E$ can be split into the components
$E_{\alpha\beta}=E_{\alpha\beta}^{||}+E_{\alpha\beta}^{\bot}+E_{\alpha\beta}^{T}$,
where the first two again can be derived from a scalar and a vector potential
$E$ and E, respectively, and the last one fulfils the transversity conditions
for tensors

$\displaystyle
E_{\alpha\beta}^{||}=(\partial_{\alpha}\partial_{\beta}-\frac{1}{3}\delta_{\alpha\beta}\nabla^{2})E,$
(43) $\displaystyle
E_{\alpha\beta}^{\bot}=-\frac{1}{2}(\partial_{\beta}E_{\alpha}+\partial_{\alpha}E_{\beta})\quad\text{with}\quad\nabla\cdot\textbf{E}=0,$
(44)
$\displaystyle\delta^{\alpha\gamma}\partial_{\gamma}E^{T}_{\alpha\beta}=0\quad\text{and}\quad\delta^{\alpha\beta}E^{T}_{\alpha\beta}=0.$
(45)

The invariant form of the density variation $\delta:=\delta\rho/\rho$ under a
gauge transformation $x^{a}\to x^{a}+\epsilon\xi^{a}$ reads

$\displaystyle\tilde{\delta}=\delta+3\frac{a^{\prime}}{a}(1+\omega)\xi^{0}.$
(46)

With the splitting introduced above and the scalar potential
$\textbf{B}^{||}=:-\nabla B$ the gauge invariant form of the scalar
perturbations in terms of the metric perturbation parameters is

$\displaystyle\tilde{\delta}=\delta-3(1+\omega)\frac{a^{\prime}}{a}(v-B),$
(47)

where the fluid velocity perturbation is $\delta u_{\alpha}=:1/a\cdot
v_{\alpha}$, which also splits like the vector B in
$\textbf{v}=\textbf{v}^{||}+\textbf{v}^{\bot}$ with the potential
$\textbf{v}^{||}=-\nabla v$. A gauge is specified by choosing values for $v$
and $B$. Similarly, the other perturbative quantities can be made gauge
invariant. The projected comoving density gradient $\Delta_{a}$ is

$\displaystyle\mathbf{\Delta}=a\nabla\tilde{\delta}-3a^{\prime}(1+\omega)\textbf{V}_{c},$
(48)

where $\textbf{V}_{c}:=\textbf{v}^{\bot}-\textbf{B}^{\bot}$. Hence its
divergence yields

$\displaystyle\Delta=\nabla^{2}\tilde{\delta},$ (49)

due to $\nabla\cdot\textbf{v}^{\bot}=\nabla\cdot\textbf{B}^{\bot}=0$. Eq. (49)
is the desired connection between the common gauge invariant Bardeen variable
for the density perturbation and the $1+3$ scalar variations.

The shear tensor $\sigma_{\alpha\beta}$, describing GWs in the $1+3$
formalism, is closely related to the transverse and tracefree part of the
tensor perturbation $E^{\bot}_{\alpha\beta}$ which equals the commonly used
$h_{\alpha\beta}$ in the transverse tracefree gauge and is gauge invariant by
itself. The shear tensor expressed in terms of Bardeen parametrized metric
perturbation reads

$\displaystyle\sigma_{\alpha\beta}=a(\nabla_{\alpha\beta}V_{S}+\nabla_{(\alpha}V_{S\>\beta)}+E^{T\prime}_{\alpha\beta}),$
(50)

where
$\Delta_{\alpha\beta}:=\nabla_{\alpha}\nabla_{\beta}-\frac{1}{3}\delta_{\alpha\beta}\nabla^{2}$
and $V_{S}:=v-D^{\prime}$. For the purpose of this work we will only need the
relations between the projected density gradient and the shear with Bardeens
variables given in Eqs. (49) and (50), respectively. Analog expressions for
other quantities can be found in [27].

##  3 Second order density perturbations

This section is devoted to the search for an equation in which density
perturbations are sourced by GWs. As we discussed in section 2.5, to first
order the shear tensor represents the kinematical properties of GWs in the
$1+3$ covariant description. Therefore, we seek a relation between the
orthogonal projected gradient $\Delta_{a}$ and linear perturbations of the
shear tensor $\sigma^{(1)}_{ab}$. This occurs for the first time at second
perturbative order in the density gradient $\Delta_{a}^{(2)}$.

In the following we will deduce the evolution equation for $\Delta_{a}^{(2)}$
with a non-zero linear shear contribution $\sigma^{(1)}_{ab}$ from the full,
non-linear Eqs. (30) and (31).

To do so, we will choose a model for the cosmic fluid which will significantly
simplify the non-linear equations. Then, we will resolve the remaining angular
brackets in the indices and take the orthogonal projected gradient of the
equations in order to obtain scalar equations. While terms that are at least
of third order will be directly neglected during the calculation, the explicit
expansion of the remaining quantities to second order is performed after
obtaining the scalar equation. Finally we set the background cosmology to FLRW
and specialize the result for a radiation dominated fluid. During the
calculation we will set $\kappa=8\pi G=1$ and reintroduce the units at the
end. For our fluid model we impose the requirements

  1. Assumption 1:

At the background level, the matter-energy density is described by a (single
component) perfect fluid.

  2. Assumption 2:

To all orders, we assume a negligible contribution from vorticity
$\omega_{ab}=0=\omega_{a}$, current density $q_{a}=0$ and anisotropy
$\pi_{ab}=0$ in the fluid.

In this model the full non linear equations Eq. (30) and Eq. (31) reduce to

$\displaystyle\dot{\Delta}_{\langle
a\rangle}=\frac{p}{\rho}\Theta\Delta_{a}-\left(1+\frac{p}{\rho}\right)Z_{a}-\sigma_{ab}\Delta^{b},$
(51)

with

$\displaystyle\dot{Z}_{\langle a\rangle}=-\frac{2}{3}\Theta
Z_{a}-\frac{1}{2}\rho\Delta_{a}-\frac{3}{2}a{\text{D}}_{a}p+a\dot{\Theta}A_{a}+a{\text{D}}_{a}{\text{D}}^{b}A_{b}-\sigma_{ab}Z^{b}-2a{\text{D}}_{a}\sigma^{2}+2aA^{b}{\text{D}}_{a}A_{b}\,.$
(52)

Additionally we fix the relation between pressure and energy density by
imposing

  * Assumption 3:

Perfect barotropic fluid: This implies $p=\omega\rho$ with constant $\omega$
and together with Eq. (21) and $q_{a}=\pi_{ab}=0$ the acceleration is to all
orders $A_{a}=\frac{-c_{s}^{2}\Delta_{a}}{a(1+\omega)}$ due to
${\text{D}}_{a}p=\frac{\rho}{a}c_{s}^{2}\Delta_{a}$. The perturbations are
adiabatic. We also assume
${\text{D}}_{a}\omega=\dot{\omega}={\text{D}}_{a}c_{s}^{2}=\dot{c_{s}^{2}}=0$.

Our perturbation procedure will extend to second order and thus we can neglect
terms that are at least of third order in advance. This is the case for the
term
$a\left(2\sigma^{2}-A^{b}A_{b}\right)A_{a}\gtrsim\mathcal{O}(\epsilon^{3})$
since the acceleration $A_{a}$ and the shear $\sigma_{ab}$ are zero at zero
order. Applying the third assumption to Eq. (51) and Eq. (52) our equations
reduce to

$\displaystyle\dot{\Delta}_{\langle a\rangle}=$
$\displaystyle\>\omega\Theta\Delta_{a}-\left(1+\omega\right)Z_{a}-\sigma_{ab}\Delta^{b},$
(53) $\displaystyle\dot{Z}_{\langle a\rangle}=$
$\displaystyle-\frac{2}{3}\Theta
Z_{a}-\left(1+3c_{s}^{2}\right)\frac{\rho}{2}\Delta_{a}-\frac{c_{s}^{2}}{1+\omega}\dot{\Theta}\Delta_{a}-\frac{c_{s}^{2}}{a(1+\omega)}{\text{D}}_{a}\Delta-\sigma_{ab}Z^{b}$
$\displaystyle-2a{\text{D}}_{a}\sigma^{2}+\frac{c_{s}^{4}}{a(1+\omega)^{2}}{\text{D}}_{a}\left(\Delta^{b}\Delta_{b}\right).$
(54)

The next step is to deal with the projected time derivative of the density
perturbation. We expand it by applying the inverse product rule

$\displaystyle\dot{\Delta}_{\langle
a\rangle}:=h_{a}^{\;b}\dot{\Delta}_{b}=\dot{(h_{a}^{\;b}\Delta_{b})}-\dot{h}_{a}^{\;b}\Delta_{b},$
(55)

but since $h_{a}^{\;b}{\text{D}}_{b}={\text{D}}_{a}$ and
$u^{a}\Delta_{a}=u^{a}h_{a}^{\;b}\nabla_{b}=0$ we find for the two terms

$\displaystyle
h_{a}^{\;b}\Delta_{b}:=h_{a}^{\;b}\frac{a}{\rho}{\text{D}}_{b}\rho=\frac{a}{\rho}{\text{D}}_{a}\rho=\Delta_{a},$
(56)
$\displaystyle\dot{h}_{a}^{\;b}\Delta_{b}=(u_{a}A^{b}+u^{b}A_{a})\Delta_{b}=u_{a}A^{b}\Delta_{b}.$
(57)

We thus get

$\displaystyle\dot{\Delta}_{\langle
a\rangle}=\dot{\Delta}_{a}-u_{a}A^{b}\Delta_{b},$ (58)

which reflects the fact that the orthogonally projected time variation of
density inhomogeneities is the same as the complete time derivative of the
density perturbation minus the projection of the density perturbation on the
flow lines. The four acceleration in turn can be expressed by the relation
$A_{a}=-\frac{c_{s}^{2}}{a(1+\omega)}\Delta_{a}$ and we find

$\displaystyle\dot{\Delta}_{\langle
a\rangle}=\dot{\Delta}_{a}+u_{a}\frac{c_{s}^{2}}{a(1+\omega)}\Delta^{b}\Delta_{b}.$
(59)

For the expansion gradient we repeat this calculation and find

$\displaystyle\dot{Z}_{\langle
a\rangle}=h_{a}^{\;b}\dot{Z}_{b}=\dot{\left(h_{a}^{\;b}Z_{b}\right)}-\dot{h}_{a}^{\;b}Z_{b}=\dot{Z}_{a}-u_{a}A^{b}Z_{b},$
(60)

for the same reason as for the density perturbations, $h_{a}^{\;b}Z_{b}=Z_{a}$
and $u^{a}Z_{a}=0$. Plugging these identities into Eqs. (53) and (54) results
in

$\displaystyle\dot{\Delta}_{a}+u_{a}\frac{c_{s}^{2}}{a(1+\omega)}\Delta^{b}\Delta_{b}=$
$\displaystyle\;\omega\Theta\Delta_{a}-\left(1+\omega\right)Z_{a}-\sigma_{ab}\Delta^{b},$
(61)
$\displaystyle\dot{Z}_{a}+u_{a}\frac{c_{s}^{2}}{a(1+\omega)}\Delta^{b}Z_{b}=$
$\displaystyle-\frac{2}{3}\Theta
Z_{a}-\left(1+3c_{s}^{2}\right)\frac{\rho}{2}\Delta_{a}-\frac{c_{s}^{2}}{1+\omega}\dot{\Theta}\Delta_{a}-\frac{c_{s}^{2}}{a(1+\omega)}{\text{D}}_{a}\Delta$
$\displaystyle-\sigma_{ab}Z^{b}-2a{\text{D}}_{a}\sigma^{2}+\frac{c_{s}^{4}}{(1+\omega)^{2}a}{\text{D}}_{a}\left(\Delta^{a}\Delta_{a}\right).$
(62)

###  3.1 Taking the orthogonal projected gradient

We are interested in the density fluctuations described by the comoving
divergence of the density perturbations $\Delta:=a{\text{D}}^{a}\Delta_{a}$.
Hence, we will take the divergence of Eqs. (61) and (62) which will yield a
scalar equation.

#### Comoving fractional density gradient:

Let us start with the divergence of the first term in Eq. (61) which gives
according to [37]

$\displaystyle a{\text{D}}^{a}\dot{\Delta}_{a}$
$\displaystyle=ah^{ab}\nabla_{b}u^{c}\nabla_{c}\Delta_{a}+ah^{ab}u^{c}\nabla_{b}\nabla_{c}\Delta_{a}$
(63) $\displaystyle=\dot{\Delta}+\sigma^{ab}\Delta_{\langle
ab\rangle}-\omega^{ab}\Delta_{[ab]}+\frac{1}{3}a\Theta
A^{a}\Delta_{a}-aA^{a}\dot{\Delta}_{a}-aq^{a}\Delta_{a}+a\left(\sigma^{ab}+\omega^{ab}\right)\Delta_{a}A_{b}.$
(64)

Applying our assumptions to this equation sets $q_{a}=\omega_{ab}=0$. Ignoring
again terms vanishing at second order we find

$\displaystyle a{\text{D}}^{a}\dot{\Delta}_{a}$
$\displaystyle=\dot{\Delta}+\sigma^{ab}\Delta_{\langle
ab\rangle}-\frac{c_{s}^{2}}{3(1+\omega)}\Theta\Delta^{a}\Delta_{a}+\frac{c_{s}^{2}}{1+\omega}\Delta^{a}\dot{\Delta}_{a}+\mathcal{O}(\epsilon^{3})$
$\displaystyle=\dot{\Delta}+\sigma^{ab}\Delta_{\langle
ab\rangle}+\frac{c_{s}^{2}}{1+\omega}\left[\frac{1}{2}\frac{\text{d}}{\text{d}t}-\frac{1}{3}\Theta\right]\Delta^{a}\Delta_{a}.$
(65)

We have also replaced $A_{a}=-\frac{c_{s}^{2}}{a(1+\omega)}\Delta_{a}$. Using
$u^{a}{\text{D}}_{a}=0$ the second term in Eq. (61) becomes

$\displaystyle
a{\text{D}}^{a}u_{a}\frac{c_{s}^{2}}{a(1+\omega)}\Delta^{b}\Delta_{b}=\frac{c_{s}^{2}}{a(1+\omega)}\Delta^{b}\Delta_{b}\,a{\text{D}}^{a}u_{a}=\frac{c_{s}^{2}}{(1+\omega)}\Delta^{b}\Delta_{b}\Theta.$
(66)

The first term on the right hand side of Eq. (61) becomes

$\displaystyle
a{\text{D}}^{a}\omega\Theta\Delta_{a}=\omega(Z^{a}\Delta_{a}+\Theta\Delta),$
(67)

while the second term reads

$\displaystyle-(1+\omega)a{\text{D}}^{a}Z_{a}=-(1+\omega)Z.$ (68)

Moving on to the third term we note that in general the space-like constraint
on the shear is not zero but rather given by Eq. (25). However, in our model
the shear plays the role of GWs. Hence we only consider the transverse
component of the shear, thus ${\text{D}}^{b}\sigma_{ab}=0$ (following [39, 37,
30]). This implies for the last term in Eq. (61)

$\displaystyle a{\text{D}}^{a}(-\sigma_{ab}\Delta^{b})$
$\displaystyle=-(a\underbrace{{\text{D}}^{a}\sigma_{ab}\Delta^{b}}_{=0}+a\,\sigma_{ab}{\text{D}}^{a}\Delta^{b})$
(69) $\displaystyle=-\sigma_{ab}\left(\frac{1}{2}\Delta h^{ba}+\Delta^{\langle
ba\rangle}+\Delta^{[ba]}\right)$ (70)
$\displaystyle=-\sigma_{ab}\Delta^{\langle ab\rangle},$ (71)

where we have made use of the fact that the shear is also tracefree
$\sigma_{ab}h^{ab}=\sigma^{a}{}_{a}=0$ and that the complete contraction of an
antisymmetric with a symmetric tensor vanishes. Altogether, the differential
equation for density fluctuations becomes

$\displaystyle\dot{\Delta}=\omega(Z^{a}\Delta_{a}+\Theta\Delta)-(1+\omega)Z-2\sigma_{ab}\Delta^{\langle
ab\rangle}-\frac{c_{s}^{2}}{1+\omega}\left(\frac{2}{3}\Theta+\frac{1}{2}\frac{\text{d}}{\text{d}t}\right)\Delta^{a}\Delta_{a}.$
(72)

#### Comoving volume expansion:

Now we move on to the differential equation for the volume expansion gradient
$Z_{a}$ in Eq. (62). Starting on the left hand side analogously to
$\dot{\Delta}$ we find for the first term

$\displaystyle a{\text{D}}^{a}\dot{Z}_{a}$
$\displaystyle=\dot{Z}+\sigma^{ab}Z_{\langle ab\rangle}+\frac{1}{3}a\Theta
A^{a}Z_{a}-aA^{a}\dot{Z}_{a}$ (73)
$\displaystyle=\dot{Z}+\sigma^{ab}Z_{\langle
ab\rangle}-\frac{c_{s}^{2}}{3(1+\omega)}\Theta\Delta^{a}Z_{a}+\frac{c_{s}^{2}}{1+\omega}\Delta^{a}\dot{Z}_{a}$
(74)

and for the second term

$\displaystyle
a{\text{D}}^{a}u_{a}\frac{c_{s}^{2}}{a(1+\omega)}\Delta^{b}Z_{b}=\frac{c_{s}^{2}}{1+\omega}\Delta^{b}Z_{b}\Theta.$
(75)

Taking the comoving divergence of the first line on the right hand side of Eq.
(62) gives

$\displaystyle-\frac{2}{3}\left(Z^{a}Z_{a}+\Theta
Z\right)-(1+3c_{s}^{2})\frac{1}{2}(\rho\Delta+\rho\Delta^{a}\Delta_{a})-\frac{c_{s}^{2}}{1+\omega}\left(a\Delta_{a}{\text{D}}^{a}\dot{\Theta}+\dot{\Theta}\Delta\right)-\frac{c_{s}^{2}}{1+\omega}{\text{D}}^{2}\Delta,$
(76)

where we have applied the definition $a{\text{D}}_{a}\Theta=:Z_{a}$ as well as
$a{\text{D}}^{a}Z_{a}=:Z$ and used $a{\text{D}}_{a}\rho=\rho\Delta_{a}$. In
the second line of Eq. (62), we find for the comoving divergences

$\displaystyle-\sigma_{ab}Z^{\langle
ab\rangle},\quad-2a^{2}{\text{D}}^{2}\sigma^{2}\quad\text{and}\quad+\frac{c_{s}^{4}}{(1+\omega)^{2}}{\text{D}}^{2}\left(\Delta^{a}\Delta_{a}\right),$
(77)

where we have used in the first term
$a{\text{D}}_{b}Z_{a}=1/3Zh_{ab}+Z_{\langle ab\rangle}+Z_{[ab]}$. So far we
find for the evolution of the volume expansion gradient

$\displaystyle\dot{Z}=$
$\displaystyle-\frac{2}{3}\frac{c_{s}^{2}}{1+\omega}\Theta\Delta^{a}Z_{a}-\frac{c_{s}^{2}}{1+\omega}\Delta^{a}\dot{Z}_{a}$
$\displaystyle-\frac{2}{3}(Z^{a}Z_{a}+\Theta
Z)-\frac{1}{2}(1+3c_{s}^{2})(\rho\Delta+\rho\Delta^{a}\Delta_{a})$
$\displaystyle-\frac{c_{s}^{2}}{1+\omega}{\text{D}}^{2}\Delta-2\sigma_{ab}Z^{\langle
ab\rangle}-2a^{2}{\text{D}}^{2}\sigma^{2}$
$\displaystyle-\frac{c_{s}^{2}}{1+\omega}(\Delta_{a}a{\text{D}}^{a}\dot{\Theta}+\dot{\Theta}\Delta)$
$\displaystyle+\frac{c_{s}^{4}}{(1+\omega)^{2}}{\text{D}}^{2}\left(\Delta^{a}\Delta_{a}\right).$
(78)

###  3.2 Second order differential equation for second order perturbations

Taking a further time derivative of Eq. (72) leads to an equation of motion
for $\Delta$ which reads

$\displaystyle\ddot{\Delta}=\omega(\dot{Z}_{a}\Delta^{a}+Z^{a}\dot{\Delta}_{a}+\dot{\Theta}\Delta+\Theta\dot{\Delta})-(1+\omega)\dot{Z}$
$\displaystyle-2\frac{\text{d}}{\text{d}t}\left(\sigma_{ab}\Delta^{\langle
ab\rangle}\right)-\frac{c_{s}^{2}}{1+\omega}\left(\frac{2}{3}\dot{\Theta}+\frac{2}{3}\Theta\frac{\text{d}}{\text{d}t}+\frac{1}{2}\frac{\text{d}^{2}}{\text{d}t^{2}}\right)\Delta^{a}\Delta_{a}.$
(79)

This equation still depends on the volume expansion gradient. In what follows
we eliminate this dependence. First, we substitute $\dot{Z}$ from Eq. (78) in
Eq. (79) and find

$\displaystyle\ddot{\Delta}=$ $\displaystyle\frac{2}{3}(1+\omega)\Theta
Z+\omega\Theta\dot{\Delta}+(\omega+c_{s}^{2})\dot{\Theta}\Delta+\frac{(1+\omega)(1+3c_{s}^{2})}{2}\rho\Delta+c_{s}^{2}{\text{D}}^{2}\Delta$
$\displaystyle-2\frac{\text{d}}{\text{d}t}\left(\sigma_{ab}\Delta^{\langle
ab\rangle}\right)+2(1+\omega)\sigma_{ab}Z^{\langle
ab\rangle}+2a^{2}(1+\omega){\text{D}}^{2}\sigma^{2}$
$\displaystyle-\frac{c_{s}^{2}}{1+\omega}\left[\frac{2}{3}\dot{\Theta}-\frac{(1+\omega)^{2}(1+3c_{s}^{2})}{c_{s}^{2}}\frac{\rho}{2}+c_{s}^{2}{\text{D}}^{2}+\frac{2}{3}\Theta\frac{\text{d}}{\text{d}t}+\frac{1}{2}\frac{\text{d}^{2}}{\text{d}t^{2}}\right]\Delta^{a}\Delta_{a}$
$\displaystyle+(\omega+c_{s}^{2})\dot{Z}^{a}\Delta_{a}+(\omega\dot{\Delta}_{a}+\frac{2}{3}c_{s}^{2}\Theta\Delta_{a}+\frac{2}{3}(1+\omega)Z_{a})Z^{a}+c_{s}^{2}(\Delta_{a}a{\text{D}}^{a}\dot{\Theta}).$
(80)

We then replace $\dot{\Theta}$ by the Raychaudhuri equation (22)

$\displaystyle\dot{\Theta}=-\frac{1}{3}\Theta^{2}-\frac{1}{2}(1+3\omega)\rho-2\sigma^{2}-\frac{c_{s}^{2}}{a^{2}(1+\omega)}\Delta+\frac{c_{s}^{4}}{a^{2}(1+\omega)^{2}}\Delta_{a}\Delta^{a}+\Lambda$
(81)

and the scalar version of the volume expansion gradient

$\displaystyle Z=-\frac{1}{1+\omega}\left(\dot{\Delta}-\omega
Z^{a}\Delta_{a}-\omega\Theta\Delta+2\sigma_{ab}\Delta^{\langle
ab\rangle}+\frac{c_{s}^{2}}{1+\omega}\left(\frac{2}{3}\Theta+\frac{1}{2}\frac{\text{d}}{\text{d}t}\right)\Delta^{a}\Delta_{a}\right),$
(82)

which is taken from Eq. (72). The spatial gradient of the Raychaudhuri
equation yields

$\displaystyle
c_{s}^{2}\Delta_{a}a{\text{D}}^{a}\dot{\Theta}=-c_{s}^{2}\frac{2}{3}\Theta\Delta_{a}Z^{a}-\frac{c_{s}^{2}}{2}(1+3\omega)\rho\Delta_{a}\Delta^{a}-\frac{c_{s}^{4}}{(1+\omega)a}\Delta_{a}{\text{D}}^{a}\Delta.$
(83)

Inserting the last three equations into Eq. (80) we find

$\displaystyle\ddot{\Delta}+\left(\frac{2}{3}-\omega\right)$
$\displaystyle\Theta\dot{\Delta}-\left((\omega-
c_{s}^{2})\frac{\Theta^{2}}{3}+(1+2c_{s}^{2}-3\omega^{2})\frac{\rho}{2}+(\omega+c_{s}^{2})\Lambda+c_{s}^{2}{\text{D}}^{2}\right)\Delta$
$\displaystyle=$
$\displaystyle-2\left(\frac{2}{3}\Theta+\frac{\text{d}}{\text{d}t}\right)\left(\sigma_{ab}\Delta^{\langle
ab\rangle}\right)+2(1+\omega)\sigma^{ab}Z_{\langle
ab\rangle}+2a^{2}(1+\omega){\text{D}}^{2}\sigma^{2}$
$\displaystyle-\frac{c_{s}^{2}}{1+\omega}\left[\frac{2}{9}\Theta^{2}-\frac{1}{c_{s}^{2}}\left((1+\omega)^{2}+\left(\frac{8}{3}+4\omega\right)c_{s}^{2}\right)\frac{\rho}{2}+c_{s}^{2}{\text{D}}^{2}+\frac{2}{3}\Lambda+\Theta\frac{\text{d}}{\text{d}t}+\frac{1}{2}\frac{\text{d}^{2}}{\text{d}t^{2}}\right]\Delta^{a}\Delta_{a}$
$\displaystyle+(\omega+c_{s}^{2})\dot{Z}^{a}\Delta_{a}+\left(\omega\dot{\Delta}_{a}+\frac{2}{3}\Theta\omega\Delta_{a}+\frac{2}{3}(1+\omega)Z_{a}\right)Z^{a}$
$\displaystyle-
c_{s}^{2}\frac{\omega+c_{s}^{2}}{(1+\omega)a^{2}}\Delta^{2}-\frac{c_{s}^{4}}{(1+\omega)a}\Delta_{a}{\text{D}}^{a}\Delta.$
(84)

Let us emphasize once more that products of three variables for which
$S^{(0)}=0$ are at least of third order and are thus neglected here.

#### Perturbative expansion:

So far we have neglected terms that are at least of third order in a
perturbative expansion of the dynamical variables. To complete the
perturbative analysis, we expand the remaining variables and truncate the
series at second order

$\displaystyle\Delta\approx\Delta^{(0)}+\Delta^{(1)}+\Delta^{(2)}\equiv\Delta^{(1)}+\Delta^{(2)}\,,$
(85) $\displaystyle Z\approx Z^{(0)}+Z^{(1)}+Z^{(2)}\equiv Z^{(1)}+Z^{(2)}\,,$
(86)
$\displaystyle\sigma\approx\sigma^{(0)}+\sigma^{(1)}+\sigma^{(2)}\equiv\sigma^{(1)}+\sigma^{(2)}\,,$
(87) $\displaystyle\Theta\approx\Theta^{(0)}+\Theta^{(1)}+\Theta^{(2)}\,,$
(88) $\displaystyle\rho\approx\rho^{(0)}+\rho^{(1)}+\rho^{(2)}.$ (89)

Vector and tensor versions of a quantity are expanded in the same manner as
their scalar counterparts given above. In Eq. (84) the shear couples either to
itself, to $\Delta$ or $Z$ whose zeroth order term is zero. Thus in an
expansion to second order only the first order perturbation term of the shear
will survive and thus we can immediately put $\sigma\approx\sigma^{(1)}$. Also
note that $\rho^{(2)}$ and $\Theta^{(2)}$ will not occur in our perturbed
formula because $\rho$ and $\Theta$ always appear in combination with a
linearly gauge independent quantity for which the zeroth order term is zero.

Expanding the variables in Eq. (84) in this manner we get

$\displaystyle\ddot{\Delta}^{(2)}$
$\displaystyle+\left(\frac{2}{3}-\omega\right)\Theta^{(0)}\dot{\Delta}^{(2)}-\left(\frac{1}{3}(\omega-
c_{s}^{2})\Theta^{(0)\,2}+(1+2c_{s}^{2}-3\omega^{2})\frac{\rho^{(0)}}{2}+(\omega+c_{s}^{2})\Lambda+c_{s}^{2}{\text{D}}^{2}\right)\Delta^{(2)}$
$\displaystyle=-\left(\frac{2}{3}-\omega\right)\Theta^{(1)}\dot{\Delta}^{(1)}+\left(\frac{2}{3}(\omega-
c_{s}^{2})\Theta^{(0)}\Theta^{(1)}+(1+2c_{s}^{2}-3\omega^{2})\frac{\rho^{(1)}}{2}\right)\Delta^{(1)}$
$\displaystyle-2\left(\frac{2}{3}\Theta^{(0)}+\frac{\text{d}}{\text{d}t}\right)\left(\sigma_{ab}^{(1)}\Delta^{(1)\,\langle
ab\rangle}\right)+2(1+\omega)\sigma^{(1)\,ab}Z_{\langle
ab\rangle}^{(1)}+2a^{2}(1+\omega){\text{D}}^{2}\sigma^{(1)\,2}$
$\displaystyle-\frac{c_{s}^{2}}{1+\omega}\left[\frac{2}{9}\Theta^{(0)\,2}-\frac{1}{c_{s}^{2}}\left((1+\omega)^{2}+\left(\frac{8}{3}+4\omega\right)c_{s}^{2}\right)\frac{\rho^{(0)}}{2}+c_{s}^{2}{\text{D}}^{2}+\frac{2}{3}\Lambda+\Theta^{(0)}\frac{\text{d}}{\text{d}t}+\frac{1}{2}\frac{\text{d}^{2}}{\text{d}t^{2}}\right]\Delta^{(1)\,a}\Delta_{a}^{(1)}$
$\displaystyle+(\omega+c_{s}^{2})\dot{Z}^{(1)\,a}\Delta_{a}^{(1)}+\left(\omega\dot{\Delta}^{(1)}_{a}+\frac{2}{3}\Theta^{(0)}\omega\Delta_{a}^{(1)}+\frac{2}{3}(1+\omega)Z_{a}^{(1)}\right)Z^{(1)\,a}$
$\displaystyle-
c_{s}^{2}\frac{\omega+c_{s}^{2}}{(1+\omega)a^{2}}\Delta^{(1)\,2}-\frac{c_{s}^{4}}{(1+\omega)a}\Delta_{a}^{(1)}{\text{D}}^{a}\Delta^{(1)}.$
(90)

We have ordered the terms in this equation such that second order quantities
appear on the left hand side of the equation and combinations of first order
variables appear on the right hand side. At this point it is also useful to
replace the contractions with $Z_{a}^{(1)}$ and $\dot{Z}_{a}^{(1)}$. To do so,
we use Eq. (53) which yields

$\displaystyle Z_{a}=\frac{\omega\Theta\Delta_{a}-\dot{\Delta}_{\langle
a\rangle}-\sigma_{ab}\Delta^{b}}{1+\omega},$ (91) $\displaystyle\Rightarrow\,$
$\displaystyle
Z_{a}^{(1)}=\frac{\omega\Theta^{(0)}\Delta_{a}^{(1)}-\dot{\Delta}_{a}^{(1)}}{1+\omega}+\mathcal{O}(\epsilon^{2}).$
(92)

Similarly, from Eqs. (62) and (92) we get for the time derivative of
$Z_{a}^{(1)}$

$\displaystyle\dot{Z}_{a}^{(1)}=$
$\displaystyle\frac{1}{3}\frac{\Theta^{(0)\,2}}{1+\omega}\left(c_{s}^{2}-2\omega\right)\Delta_{a}^{(1)}+\frac{2}{3}\frac{1}{1+\omega}\Theta^{(0)}\dot{\Delta}_{a}^{(1)}-\frac{c_{s}^{2}}{a(1+\omega)}{\text{D}}^{a}\Delta^{(1)}$
$\displaystyle+\left(\frac{c_{s}^{2}}{1+\omega}(1+3\omega)-(1+3c_{s}^{2})\right)\frac{\rho^{(0)}}{2}\Delta_{a}^{(1)}.$
(93)

The respective terms in Eq. (90) become

$\displaystyle(\omega+c_{s}^{2})$
$\displaystyle\dot{Z}^{(1)\,a}\Delta_{a}^{(1)}=$
$\displaystyle\frac{1}{3}\frac{\omega+c_{s}^{2}}{1+\omega}(c_{s}^{2}-2\omega)\Theta^{(0)\,2}\Delta_{a}^{(1)}\Delta^{(1)\,a}+\frac{1}{3}\frac{\omega+c_{s}^{2}}{1+\omega}\Theta^{(0)}\frac{\text{d}}{\text{d}t}(\Delta^{(1)\,a}\Delta_{a}^{(1)})-\frac{c_{s}^{2}(c_{s}^{2}+\omega)}{a(1+\omega)}\Delta_{a}^{(1)}{\text{D}}^{a}\Delta^{(1)}$
$\displaystyle+(\omega+c_{s}^{2})\left(\frac{c_{s}^{2}}{1+\omega}(1+3\omega)-(1+3c_{s}^{2})\right)\frac{\rho^{(0)}}{2}\Delta_{a}^{(1)}\Delta^{(1)\,a},$
(94) and $\displaystyle(\omega\dot{\Delta}^{(1)}_{a}+$
$\displaystyle\frac{2}{3}\Theta^{(0)}\omega\Delta_{a}^{(1)}+\frac{2}{3}(1+\omega)Z_{a}^{(1)})Z^{(1)\,a}=$
$\displaystyle\frac{\omega-2}{2}\frac{\omega}{1+\omega}\Theta^{(0)}\frac{\text{d}}{\text{d}t}(\Delta^{(1)\,a}\Delta_{a}^{(1)})+\frac{\frac{2}{3}-\omega}{1+\omega}\dot{\Delta}^{(1)\,a}\dot{\Delta}_{a}^{(1)}+\frac{4}{3}\frac{\omega^{2}}{1+\omega}\Theta^{(0)\,2}\Delta^{(1)\,a}\Delta_{a}^{(1)}.$
(95)

With these we can eliminate $\dot{Z}_{a}$ and $Z_{a}$ completely and our
perturbed second order differential equation for second order density
perturbations in a perfect fluid with shear reads

$\displaystyle\ddot{\Delta}^{(2)}$
$\displaystyle+\left(\frac{2}{3}-\omega\right)\Theta^{(0)}\dot{\Delta}^{(2)}-\left(\frac{1}{3}(\omega-
c_{s}^{2})\Theta^{(0)\,2}+(1+2c_{s}^{2}-3\omega^{2})\frac{\rho^{(0)}}{2}+(\omega+c_{s}^{2})\Lambda+c_{s}^{2}{\text{D}}^{2}\right)\Delta^{(2)}$
$\displaystyle=-\left(\frac{2}{3}-\omega\right)\Theta^{(1)}\dot{\Delta}^{(1)}+\left(\frac{2}{3}(\omega-
c_{s}^{2})\Theta^{(0)}\Theta^{(1)}+(1+2c_{s}^{2}-3\omega^{2})\frac{\rho^{(1)}}{2}\right)\Delta^{(1)}$
$\displaystyle-2\left(\frac{2}{3}\Theta^{(0)}+\frac{\text{d}}{\text{d}t}\right)\left(\sigma_{ab}^{(1)}\Delta^{(1)\,\langle
ab\rangle}\right)+2(1+\omega)\sigma_{ab}^{(1)}Z^{(1)\,\langle
ab\rangle}+2a^{2}(1+\omega){\text{D}}^{2}\sigma^{(1)\,2}$
$\displaystyle-\frac{c_{s}^{2}}{1+\omega}\left[\left(\frac{2}{3}-\frac{2\omega^{2}+c_{s}^{4}-\omega
c_{s}^{2}}{c_{s}^{2}}\right)\frac{\Theta^{(0)\,2}}{3}+\frac{1}{c_{s}^{2}}\left(-1-\omega-\omega
c_{s}^{2}-\frac{5}{3}c_{s}^{2}+2c_{s}^{4}\right)\frac{\rho^{(0)}}{2}+c_{s}^{2}{\text{D}}^{2}+\frac{2}{3}\Lambda\right.$
$\displaystyle\left.+\left(\frac{2}{3}+\frac{2\omega}{3c_{s}^{2}}-\frac{\omega^{2}}{2c_{s}^{2}}\right)\Theta^{(0)}\frac{\text{d}}{\text{d}t}+\frac{1}{2}\frac{\text{d}^{2}}{\text{d}t^{2}}\right]\Delta^{(1)\,a}\Delta_{a}^{(1)}$
$\displaystyle+\frac{\frac{2}{3}-\omega}{1+\omega}\dot{\Delta}_{a}^{(1)}\dot{\Delta}^{(1)\,a}-\frac{c_{s}^{2}(c_{s}^{2}+\omega)}{(1+\omega)a^{2}}\Delta^{(1)\,2}-\frac{c_{s}^{2}(2c_{s}^{2}+\omega)}{(1+\omega)a}\Delta_{a}^{(1)}{\text{D}}^{a}\Delta^{(1)}.$
(96)

We fix the background cosmology by introducing a further assumption.

  * Assumption 4:

As background we choose an FLRW universe: $\Theta^{(0)}(t)=3H(t)$ and
$\rho^{(0)}=\rho(t)$ is given by the Friedmann and continuity equations

$\displaystyle\dot{H}=-H^{2}-\frac{1}{6}(\rho+3p)+\frac{1}{3}\Lambda,\quad\dot{\rho}=-3H(\rho+p),\,\text{and}$
$\displaystyle H^{2}=\frac{1}{3}\rho-\frac{K}{a^{2}}+\frac{1}{3}\Lambda.$

On this background our equation eventually yields

$\displaystyle\ddot{\Delta}^{(2)}$
$\displaystyle+2H\left(1-\frac{3}{2}\omega\right)\dot{\Delta}^{(2)}$
$\displaystyle-\left[\frac{3}{2}\left(1+2\omega-3\omega^{2}\right)H^{2}-\frac{1}{2}\left(1-2\omega-3\omega^{2}\right)\Lambda+\left(1+2c_{s}^{2}-3\omega^{2}\right)\frac{3K}{2a^{2}}+c_{s}^{2}{\text{D}}^{2}\right]\Delta^{(2)}$
$\displaystyle=$
$\displaystyle-\left(\frac{2}{3}-\omega\right)\Theta^{(1)}\dot{\Delta}^{(1)}+\left(\frac{1}{3}(\omega-
c_{s}^{2})6H\Theta^{(1)}+(1+2c_{s}^{2}-3\omega^{2})\frac{\rho^{(1)}}{2}\right)\Delta^{(1)}$
$\displaystyle-2\left(2H+\frac{\text{d}}{\text{d}t}\right)\sigma_{ab}^{(1)}\Delta^{(1)\,\langle
ab\rangle}+2(1+\omega)\sigma_{ab}^{(1)}Z^{(1)\,\langle
ab\rangle}+2(1+\omega)a^{2}{\text{D}}^{2}\sigma^{(1)\,2}$
$\displaystyle-\frac{c_{s}^{2}}{1+\omega}\left[\left(-1-\omega+\frac{7}{3}\omega
c_{s}^{2}-\frac{5}{3}c_{s}^{2}-4\omega^{2}\right)\frac{\rho^{(0)}}{2c_{s}^{2}}+c_{s}^{2}{\text{D}}^{2}+\left(\omega-2\frac{\omega^{2}}{c_{s}^{2}}-c_{s}^{2}+\frac{2}{3}\right)\Lambda\right.$
$\displaystyle\left.+\left(\frac{2}{3}+\frac{2\omega}{3c_{s}^{2}}-\frac{\omega^{2}}{2c_{s}^{2}}\right)3H\frac{\text{d}}{\text{d}t}+\frac{1}{2}\frac{\text{d}^{2}}{\text{d}t^{2}}\right]\Delta^{(1)\,a}\Delta_{a}^{(1)}$
$\displaystyle+\frac{\frac{2}{3}-\omega}{1+\omega}\dot{\Delta}_{a}^{(1)}\dot{\Delta}^{(1)\,a}-\frac{c_{s}^{2}(c_{s}^{2}+\omega)}{(1+\omega)a^{2}}\Delta^{(1)\,2}-\frac{c_{s}^{2}(2c_{s}^{2}+\omega)}{(1+\omega)a}\Delta_{a}^{(1)}{\text{D}}^{a}\Delta^{(1)}.$
(97)

In this equation the second order density perturbations on the left hand side
are sourced by couplings of first order perturbations on the right hand side.
In particular, we find source terms in which the shear tensor couples to first
order density perturbations but also to itself. The next step is to study our
equation in the two important regimes where either radiation or matter
dominates. For that we reintroduce $\kappa:=8\pi G$.

#### Matter dominated universe:

For a matter dominated universe we have $\omega=c_{s}^{2}=0$. The linearized
evolution equation for the shear tensor and the Gauss-Codazzi equation [37]
determine the projected Ricci tensor by the shear $R_{\langle
ab\rangle}=-3H\sigma_{ab}-\dot{\sigma}_{ab}$. For super-horizon shear modes
${\text{D}}^{2}\sigma^{2}=0$ in a flat universe $K=0$, this identity together
with the relations $\dot{\Delta}_{\langle ab\rangle}=-Z_{\langle ab\rangle}$
and $\dot{\Delta}_{a}=-Z_{a}$ (details see ref. [37]) reduces Eq. (97) to

$\displaystyle\ddot{\Delta}^{(2)}+2H\dot{\Delta}^{(2)}-\frac{1}{2}(3H^{2}-\Lambda)\Delta^{(2)}=$
$\displaystyle\frac{3H^{2}}{2}\kappa\Delta_{a}^{(1)}\Delta^{(1)\,a}+\frac{2}{3}Z_{a}^{(1)}Z^{(1)\,a}+2H\sigma_{ab}^{(1)}\Delta^{(1)\,\langle
ab\rangle}$ $\displaystyle+4\sigma_{ab}^{(1)}Z^{(1)\,\langle
ab\rangle}+2R_{\langle ab\rangle}^{(1)}\Delta^{(1)\,\langle ab\rangle}$
$\displaystyle-\frac{2}{3}\Theta^{(1)}\dot{\Delta}^{(1)}+\frac{1}{2}\rho^{(1)}\Delta^{(1)}.$
(98)

which partially reproduces the result in [37] for $\Lambda=0$. We find two
additional terms
$-2H\Theta^{(1)}\dot{\Delta}^{(1)}+\frac{1}{2}\rho^{(1)}\Delta^{(1)}$ which
were missed by the reduction procedure used in [37].

#### Radiation dominated universe:

Important for this work is the evolution of the perturbations during radiation
domination where we have $\omega=c_{s}^{2}=1/3$. For our purpose it is also
convenient to eliminate $Z_{\langle ab\rangle}$ such that we get an equation
only depending on $\sigma_{ab}$ and $\Delta$ in its various forms. We also
take a flat universe with $K=0$ and $\Lambda=0$. Again by using Eq. (61) the
volume expansion gradients yield to linear order
$(1+\omega)Z_{a}^{(1)}=\omega\Theta\Delta_{a}^{(1)}-\dot{\Delta}^{(1)}$.
Taking the comoving derivative $a{\text{D}}^{b}$ of the previous expression,
using the decomposing rule $a{\text{D}}_{b}Z_{a}=1/3h_{ab}Z+Z_{\langle
ab\rangle}+Z_{[ab]}$ (analogously for $\Delta_{a}$) and the linear rule
$a{\text{D}}_{b}\dot{\Delta}_{a}^{(1)}=a\frac{\text{d}}{\text{d}t}({\text{D}}_{b}\Delta_{a}^{(1)})$
we can estimate

$\displaystyle Z_{\langle
ab\rangle}^{(1)}=\frac{\omega\Theta^{(0)}\Delta^{(1)}h_{ab}+\omega\Theta^{(0)}\Delta_{\langle
ab\rangle}^{(1)}-\dot{\Delta}_{\langle
ab\rangle}^{(1)}}{1+\omega}+\mathcal{O}(\epsilon^{2}).$ (99)

Finally our equation yields in a flat universe without cosmological constant

$\displaystyle\ddot{\Delta}^{(2)}+$ $\displaystyle
H\dot{\Delta}^{(2)}-2H^{2}\Delta^{(2)}-\frac{1}{3}{\text{D}}^{2}\Delta^{(2)}$
$\displaystyle\begin{rcases*}=-2\left(H\sigma_{ab}^{(1)}\Delta^{(1)\,\langle
ab\rangle}+\dot{\sigma}_{ab}^{(1)}\Delta^{(1)\,\langle
ab\rangle}+2\sigma_{ab}^{(1)}\dot{\Delta}^{(1)\,\langle
ab\rangle}\right)+\frac{8}{3}a^{2}{\text{D}}^{2}\sigma^{(1)\,2}\end{rcases*}\text{{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}GW
sources}}$
$\displaystyle\begin{rcases*}\phantom{=}-\frac{1}{3}\Theta^{(1)}\dot{\Delta}^{(1)}+\frac{2}{3}\rho^{(1)}\Delta^{(1)}\\\
\phantom{=}+\left(2H^{2}-\frac{1}{12}{\text{D}}^{2}\right)\Delta_{a}^{(1)}\Delta^{(1)\,a}-\frac{7}{4}H\dot{\Delta}_{a}^{(1)}\Delta^{(1)\,a}-\frac{1}{4}\ddot{\Delta}_{a}^{(1)}\Delta^{(1)\,a}\\\
\phantom{=}-\frac{1}{6}\frac{1}{a^{2}}\Delta^{(1)\,2}-\frac{1}{4}\frac{1}{a}\Delta_{a}^{(1)}{\text{D}}^{a}\Delta^{(1)}.\end{rcases*}\text{{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}Pure
density sources}}$ (100)

##  4 First order phase transitions and GWs

First order phase transitions occur when a configuration does not minimize the
energy anymore. In this process the minimum $\langle\phi\rangle$ (order
parameter) of the temperature dependent potential $V_{T}(\phi)$ evolves from a
symmetric (unordered) phase $\langle\phi\rangle=0$ to an asymmetric (ordered)
phase $\langle\phi\rangle\neq 0$. If the former minimum and the new minimum in
the potential are separated by a potential barrier then the order parameter
does not smoothly roll into the new minimum but rather jumps or tunnels (see
figure 1). The barrier of a FPT prevents the system from continuously relaxing
into a new state with lower energy which results in latent heat being stored
and later released in a shorter time interval. This delayed energy releasable
makes FPTs particularly interesting in comparison to other phase transitions.
The temperature at which the two minima are degenerate is called the critical
temperature $T_{c}$ while the temperature at which the probability per volume
element of reaching the new minimum is unity is called nucleation temperature
$T_{\text{nuc}}$.

In the early universe such phase transitions could have happened and are
realized in many extensions of the SM [5, 6, 7, 8, 9, 10, 11, 13, 14, 15, 16]
where a Lagrangian $\mathcal{L}(\\{\phi_{i}\\})$ with extra fields $\phi_{i}$,
symmetries and couplings is introduced. Here the order parameter is the vacuum
expectation value (VEV) of the field which acquires a non zero mass in the low
temperature phase. This results in a ”spontaneous breaking” of the involved
symmetries, i.e. a non-linear realization of the symmetry in the vacuum at
zero temperature. In a FPT the field does not accept the new VEV everywhere in
space at the same time. Instead bubbles with the new VEV inside nucleate with
some initial sizes and begin to spatially expand into regions formerly
occupied by the high temperature, symmetric VEV. Hereby the bubbles release
the stored energy, the latent heat fraction

$\displaystyle\alpha=\frac{\rho_{\text{vac}}}{\rho_{\text{rad}}(T_{\text{nuc}})},$
(101)

where $\rho_{\text{vac}}\sim|V_{T_{\text{nuc}}}(\langle\phi\rangle)|$, and
release it in the form of motion but also while their surfaces, called walls,
eventually collide. The time needed until the field $\phi$ has acquired the
new VEV everywhere in the universe is denoted by $1/\beta$ (see figure 1
second plot) and is the inverse of the nucleation rate per Hubble volume

$\displaystyle\Gamma=\Gamma_{\text{nuc}}\exp(\beta(t-t_{\text{nuc}})).$ (102)

The rate $\beta$, in turn, is deduced from the $O(3)$-symmetric, effective
action [40]

$\displaystyle
S_{3}=4\pi\int_{0}^{\infty}r^{2}\text{d}r\left(\frac{1}{2}\left(\frac{\text{d}\phi}{\text{d}r}\right)^{2}+V_{T}(\phi)\right),\quad\text{via}\quad\frac{\beta}{H_{\text{nuc}}}=\left.T_{\text{nuc}}\frac{\text{d}S_{3}}{\text{d}T}\right|_{T_{\text{nuc}}},$
(103)

whose minimization determines the tunneling trajectory from the former,
symmetric vacuum to the low temperature, forming vacuum. In this way, via the
form of the potential $V_{T}(\phi)\subset\mathcal{L}(\phi)$, the phase
transition parameters $\alpha$ and $\beta$ are directly connected to the
details of the particle physics model. Whether or not a certain model posseses
a potential barrier and a FPT depends therefore on both the model details and
the choice of coupling parameters.

The formation of bubbles has an important implication. By their expansion and
collision, they induce three different forms of anisotropic stress into the
fluid, which in turn generate GWs.

  1. 1.

Bubble collision: The collision of the forming and expanding bubbles leads to
an anisotropic stress that sources GWs [41, 42, 43],

  2. 2.

Turbulence: The highly ionized plasma can develop magnetic and hydrodynamical
turbulence from the percolation of the fluid induced by the bubble collisions
[44, 45],

  3. 3.

Sound waves: The bulk fluid motion produces pressure waves in the fluid that
source GWs [46, 47, 48],

such that the total abundance of GWs from a FPT reads

$\displaystyle\Omega_{\text{GW}}(k)=\Omega_{\text{BC}}(k)+\Omega_{\text{SW}}(k)+\Omega_{\text{MHD-
turb}}(k).$ (104)

In summary the important parameters are:

  * •

The nucleation temperature $T_{\text{nuc}}$ ($\sim H_{\text{nuc}}$) which sets
the scale of the released energy density in GWs $\rho_{\text{GW}}$.

  * •

The strength of the phase transition is described by the latent heat fraction
$\alpha$ and its duration by the inverse nucleation rate $\beta^{-1}$.

The strength of the GW signal also depends on the bubble wall velocity
$v_{w}$. Its calculation is more involved (see eg. [49, 50, 15, 51]).
Depending on the bubble dynamics not all of the released energy produces GWs.
The fraction of the released energy actually transmitted to the kinetic energy
of the fluid is provided by the efficiency factor $\kappa_{\text{eff}}$. If
the phase transition does not reheat the universe too much one can approximate
the temperature at which the GWs are released as $T_{*}\approx
T_{\text{nuc}}$. On the other hand in models with large $\alpha>1$ there will
be a supercooling phase which separates the nucleation temperature from the
percolation temperature $T_{\text{nuc}}\ll T_{*}$ [15]. The release of GWs
then falls together with reheating of the universe $T_{*}=T_{\text{reh}}$
which turns the universe back into radiation domination. If the reheating is
fast enough we do not need to distinguish between $H(T_{\text{nuc}})$ and
$H_{*}$ [16]. The two cases constitute very different bubble dynamics. In the
former the contribution from bubble walls to the GW signal is only relevant
for so called run-away dynamics in which the bubble walls strongly accelerate
until they reach the speed of light. However, they expand into a still
radiation dominated universe such that parts of the available energy are
absorbed by the plasma and thus the efficiency factor $\kappa_{\text{eff}}$ is
smaller than unity. In the latter case, where the bubbles propagate with the
speed of light $v_{w}=c$ in a vacuum energy dominated universe, the efficiency
factor equals unity and bubble walls are the only contribution to GW
production. The efficiency factor thus reads [15]

$\displaystyle\kappa_{\text{eff}}=1\quad\text{ for}\quad\alpha>1.$ (105)

Note, that for the purpose of this paper only the time and scale of GW
production is important and therefore we will always refer to $t_{*}$ and not
bother about $t_{\text{nuc}}$, also in the case of strong supercooling because
in any case $H(t_{*})\approx H(t_{\text{nuc}})$.

The sufficient set of parameters describing the phase transition is then
$(H_{*},\beta,\alpha,\kappa_{\text{eff}}(\alpha))$.

Figure 1: Left: Qualitative temperature dependence of a typical potential in
particle physics developing a barrier and new global minimum as temperature
drops. Right: Schematic illustration of the temperature dependence of the VEV
(order parameter) in a FPT. Shown are also the two important temperatures
$T_{c}$ and $T_{\text{nuc}}$ at which the minima are degenerated and at which
the nucleation probability reaches one bubble per Hubble volume, respectively.
The temperature at which the transition is completed is denoted by $T_{f}$.

###  4.1 Analytic description of bubble collision

We study the GW energy density originating from bubble collisions as source of
second order density perturbations following the literature, e.g. [41, 52, 43,
49, 53]. Typical central assumptions in these derivations are:

  1. 1.

Thin wall: the bubble walls are infinitesimal thin and all energy is stored on
them,

  2. 2.

Envelope approximation: Already collided walls do not source GWs. Only the
remaining envelope of collided bubbles carries energy and momentum.

  3. 3.

The phase transition performs in less than a Hubble time.

GWs originate from linear tensor perturbations (for consistency with the
literature in this subsection Latin indices are spatial and run from $1$ to
$3$.)

$\displaystyle\text{d}s^{2}=-\text{d}t^{2}+a^{2}(t)(\delta_{ij}+2h_{ij})\text{d}x^{i}\text{d}x^{j},$
(106)

by a tracefree and transverse tensor $h_{ij}(\textbf{x},t)$. The GWs propagate
according to the wave equation [54, 55] and are sourced by the transverse and
tracefree component of the anistropic stress tensor
$\Pi_{ij}^{\bot}(\textbf{x},t)$. In Fourier space ($k$ denotes comoving wave
number) the equation of motion reads

$\displaystyle\ddot{h}_{ij}(\textbf{k},t)+\frac{k^{2}}{a^{2}}h_{ij}(\textbf{k},t)=16\pi
G\,\Pi_{ij}^{\bot}(\textbf{k},t),$ (107)

Solving this equation for a given anisotropic stress tensor allows to derive
the energy density of GWs

$\displaystyle\rho_{\text{GW}}(t):=\frac{\langle\dot{h}_{ij}\dot{h}_{ij}\rangle}{8\pi
G}=\int_{0}^{\infty}\frac{k^{3}}{2\pi}|\dot{h}(k,t)|^{2}\text{d}\ln k.$ (108)

In the case of FPTs, the collision of bubbles produces an anistropic stress
tensor $\Pi_{ij}^{\bot}(\textbf{k},t)$ which drives GWs through Eq. (107).
This leads to an energy density per logarithmic frequency, which is given by
(since the FPT is short we can put $a(t)\approx a_{*}$)

$\displaystyle\Omega^{\text{log}}_{\text{GW}}(k/a_{*},t):=\frac{1}{\rho_{\text{tot}}}\frac{\text{d}\rho_{\text{GW}}}{\text{d}\ln
k}=\kappa_{\text{eff}}^{2}\left(\frac{H_{*}}{\beta}\right)^{2}\left(\frac{\alpha}{1+\alpha}\right)^{2}\Delta(k/a_{*},\beta,t,v_{w}).$
(109)

The challenge for analytical [41, 52, 42, 53] and numerical studies [56, 57,
58] is to find an expression for the dimensionless power spectrum
$\Delta(k/\beta,t,v_{w})$. The essential ingredients are the homogeneous
solution of the wave equation Eq. (107) and the power spectrum of the
anisotropic stress tensor evaluated at different times. Here we simply refer
to the literature and stick to an approximate formula from Caprini et al. [49]
which incorporates the most important features101010In fact the approximation
seems very close to more refined analyses like [53]. The latter reference
explicitly mentions that the underling assumptions listed above are especially
well fulfilled for bubbles expanding into vacuum. This is in particular
important for this work since significant impact will only be generated in
this regime. . Following this reference the dimensionless power spectrum is
well described by

$\displaystyle\Delta(k,\beta,t,v_{w})=\beta^{2}k^{3}\left|\int_{t_{*}}^{t}f(k,\tilde{t})e^{i\tilde{t}k}\text{d}\tilde{t}\right|^{2}\approx
k^{3}\beta^{2}(t-t_{*})^{2}f^{2}(k,(t+t_{*})/2),$ (110)

with the broken rational function

$\displaystyle
f(k,t)^{2}=L(t)^{2}\left(\frac{v_{w}\epsilon}{\beta}\right)\left(\frac{1+(\frac{kL}{3})^{2}}{1+(\frac{kL}{2})^{2}+(\frac{kL}{3})^{6}}\right)$
(111)

and the characteristic length- and time-functions

$\displaystyle L(t)=\frac{v_{w}}{\beta}g(t),$ (112) $\displaystyle
g(t)=4\beta^{2}(t-t_{*})\left(\frac{1}{\beta}-(t-t_{*})\right)\left[\Theta_{\text{Hv}}\left(t-t_{*}\right)\cdot\Theta_{\text{Hv}}\left(\frac{1}{\beta}-t_{*}\right)\right],$
(113)

where $\epsilon$ is a small parameter (taken to $0.01$ in the following),
$v_{w}$ is the bubble expansion speed, $t_{*}$ the starting time of the
release of GWs and $\Theta_{\text{Hv}}(t)$ is the Heaviside step-function. In
Fig. 2 we show the dimensionless power spectrum achieved from these functions.

In terms of the rescaled time $\tau:=H_{*}\cdot t$, rescaled Fourier mode
$\kappa:=ck/(a_{*}H_{*})$ (note that we put $k\to k/a_{*}$) and the relative
phase transition duration $H_{*}/\beta$ we get (also reintroducing the speed
of light $c$)

$\displaystyle k/a_{*}\cdot L(t)$ $\displaystyle=\frac{k}{a_{*}H_{*}}\cdot
H_{*}\cdot\frac{c}{\beta}\left(\frac{v_{w}}{c}\right)g(t)=\kappa\left(\frac{v_{w}}{c}\right)\left(\frac{H_{*}}{\beta}\right)\tilde{g}(\tau)=\kappa
L(\tau),$ (114)

where the time-dependent function $\tilde{g}(\tau)$ and the broken rational
function $f(\kappa,\tau)$ becomes

$\displaystyle\tilde{g}(\tau)=4\left(\frac{\beta}{H_{*}}\right)^{2}(\tau-\tau_{*})\left(\left(\frac{H_{*}}{\beta}\right)-(\tau-\tau_{*})\right)\Theta_{\text{Hv}}(\tau_{*},\tau_{*}+\frac{H_{*}}{\beta}),$
(115) $\displaystyle k^{3}f(k,t)=\kappa^{3}f^{2}(\kappa,\tau)=(\kappa
L(\tau))^{2}\left(\kappa\frac{v_{w}}{c}\frac{H_{*}}{\beta}\epsilon\right)\left(\frac{1+\left(\frac{\kappa
L(\tau)}{3}\right)^{2}}{1+\left(\frac{\kappa
L(\tau)}{2}\right)^{2}+\left(\frac{\kappa L(\tau)}{3}\right)^{6}}\right).$
(116)

Therefore, the rescaled function for the GW abundance is

$\displaystyle\Omega_{\text{GW}}(\kappa,\tau)=\kappa_{\text{eff}}^{2}(\tau-\tau_{*})^{2}\left(\frac{\alpha}{1+\alpha}\right)^{2}\int_{\kappa_{*}}^{\kappa}\text{d}\tilde{\kappa}\,\tilde{\kappa}^{2}f^{2}(\tilde{\kappa},(\tau+\tau_{*})/2),$
(117)

with $k_{*}:=a_{*}H_{*}/c$ the scale of the horizon at $t_{*}$.

The left plot in Fig. 2 shows the spectrum as a function of time for a fixed
mode while the right plot shows the time evolution of the peak
$\Delta(\kappa_{\text{peak}}(\tau),\tau)$. In Fig. 3 we show the approximation
Eq. (110) for different times. The main three features of the power spectrum
of GWs from FPT are

  * •

It peaks around $\sim k_{\text{peak}}=\frac{1.3\,\pi\,\beta}{c\cdot
L(\tau_{(\text{ev}}+\tau_{*})/2)}$ which is $\frac{2\,\pi\,\beta}{c}$ at the
end of the transition,

  * •

For small wave numbers the spectrum grows as $k^{3}$,

  * •

For large wave numbers the spectrum decreases as $k^{-1}$.

Figure 2: Left: The dimensionless power spectrum (with approximated time
integral) of GWs sourced by a FPT evaluated at halftime
$\tau=\tau_{*}+0.5\cdot H_{*}/\beta$. Right: The time evolution of the peak of
the dimensionless power spectrum. The inverse duration in terms of Hubble time
is denoted by $r_{\beta}=\frac{\beta}{H_{*}}$ and
$\bar{\tau}_{\text{ev}}:=\frac{\tau_{\text{ev}}+\tau_{*}}{2}$ is the time mean
while $\Delta\tau_{\text{ev}}:=\tau_{\text{ev}}-\tau_{*}$. Figure 3: The
dimensionless power spectrum of GWs sourced by a FPT (light red and light
green) and the approximation of the time integral in Eq. (110) (bold lines)
close to the start and the end of the FPT. The yellow curve shows the horizon
mode. The notation is the same as in Fig. 2.

##  5 Results

In this section we present the results obtained from Eq. (100) for the
following scenario (see also sketch 4):

During radiation domination a FPT is triggered at time $t_{\text{nuc}}$ and
emits GWs by bubble collision at time $t_{*}$ on sub-horizon scales $k\gg
a_{*}H_{*}$. The GWs manifest themselves as shear perturbations. The
transition completes within less than a Hubble time $t_{*}+1/\beta$, where
$\beta>H(t_{*})$. The shear distortions induce second order density
perturbations via Eq. (100). After sourcing the induced density perturbations
remain imprinted in the spectrum. Hence we need to identify the source terms,
calculate the density perturbations they induce and transfer them to the
matter power spectrum in order to estimate their impact on structure
formation.

Figure 4: Timeline and schematic flow chart of the model described in the
text. The green line depicts the time of the phase transition and the
radiation of GWs. The flow chart on top of the line shows the processes
happening during that phase on sub-horizon scales. We assume that these
processes happen during the whole period of the transition. After the
transition is completed, the induced perturbations affect the matter power
spectrum today. The two angled lines indicate that the time between the end of
the phase transition and today is much longer than the duration of the
transition.

To do so, we have to use the solution of Eq. (100) to derive the transfer
function $T(k)$. In order to induce any changes at all, the FPT must be strong
enough such that the terms not involving the shear tensor are subdominant.
Moreover, as we will see in Eq. (120) the shear is related to the GW abundance
at the transition time via
$|\kappa\sigma^{(1)}|\sim\sqrt{3\Omega_{\text{GW}}(\kappa,\tau)}$ and the
first order perturbations can be estimated as $|\Delta^{(1)}|\sim 10^{-4}$,
see Eq. (148). This implies $|\Delta^{(1)}|^{2}\sim 10^{-8}$ and the
$\sigma_{ab}\Delta^{\langle ab\rangle}$ couplings are of order
$|\kappa\sigma^{(1)}|\times 10^{-4}$. Therefore we find the pure shear term to
be the most interesting and powerful source term if
$|\sigma^{(1)}|>|\Delta^{(1)}|\sim 10^{-4}$. The latter condition is fulfilled
for a relatively wide range of phase transition parameters. For
$\alpha\to\infty$ the duration can be as small as $\beta/H_{*}\approx 100$
until the high-$\kappa$ plateau of $\rho_{\text{GW}}$ reaches a magnitude of
$10^{-4}$.

Thus, in the following we focus on the self coupling of the shear, namely

$\displaystyle\ddot{\Delta}^{(2)}+H\dot{\Delta}^{(2)}-2H^{2}\Delta^{(2)}-\frac{1}{3}{\text{D}}^{2}\Delta^{(2)}=\frac{8}{3}a^{2}{\text{D}}^{2}\sigma^{(1)\,2}.$
(118)

Following references [59, 60, 30] the shear tensor is related to the linear,
tracefree and transverse metric perturbation by
$\sigma_{ab}=a(h_{\alpha\beta})^{\prime}$ and
$\sigma^{ab}=a^{-3}(h_{\alpha\beta})^{\prime}$ (see also Eq. (50)). Recall
that in this work $a,b=0,1,2,3$ and $\alpha,\beta=1,2,3$ and the prime denotes
comoving derivatives. Using the definition of the energy density of GWs

$\displaystyle\rho_{\text{GW}}(\textbf{x},t)=\frac{(h_{\alpha\beta})^{\prime}(\textbf{x},t)(h^{\alpha\beta})^{\prime}(\textbf{x},t)}{2a^{2}\,8\pi
G}=\frac{\dot{h}_{\alpha\beta}(\textbf{x},t)\dot{h}^{\alpha\beta}(\textbf{x},t)}{16\pi
G},$ (119)

we observe that the squared shear
$\sigma^{2}(\textbf{x},t)=1/2\,\sigma_{ab}\sigma^{ab}$ is nothing but

$\displaystyle\rho_{\text{GW}}(\textbf{x},t)=\frac{1}{8\pi
G}\sigma^{2}(\textbf{x},t).$ (120)

We can also transform the divergence of the fractional energy density gradient
into a more familiar variable. To do so, we note, that in a vorticity free and
spatially flat space ($K=0$) the projected derivatives become spatial
Laplacians

$\displaystyle{\text{D}}^{a}{\text{D}}_{a}=\frac{\nabla^{2}}{a^{2}}$ (121)

and thus the divergence of the fractional density gradient can be written as
the Laplacian of some function $\tilde{\delta}$ which depends on the relative
energy density perturbation $\delta:=\delta\rho/\bar{\rho}$

$\displaystyle\Delta(\textbf{x},t)=\nabla^{2}\tilde{\delta}(\textbf{x},t).$
(122)

If $\Delta(\textbf{x},t)\equiv\Delta^{(1)}(\textbf{x},t)$ is a linear
perturbation then $\tilde{\delta}$ is equivalent to Bardeen’s variable for the
relative energy density perturbation
$\tilde{\delta}^{(1)}:=\delta+3H(1+\omega)\xi^{0}$ with the time component
$\xi^{0}$ of an arbitrary gauge transformation $x^{\mu}\to x^{\mu}+\xi^{\mu}$
[27] (see also subsection 2.6). Therefore, using Eq. (120) and Eq. (122) the
$1+3$ covariant variables can be expressed in a more standard manner.

Similarly, in a spatially flat spacetime the harmonic decomposition of the
variables reduces to standard Fourier modes (see [27, 61])

$\displaystyle
f(\textbf{x},t)=\int_{\textbf{k}}f_{\textbf{k}}\,e^{i\textbf{k}\cdot\textbf{x}},$
(123)

where k is the comoving wave vector and x is the comoving space vector.
Applying the Fourier decomposition to our Eq. (118) while the source is active
yields

$\displaystyle-k^{2}\ddot{\tilde{\delta}}_{k}^{(2)}(t)+H(t)(-k^{2})\dot{\tilde{\delta}}_{k}^{(2)}(t)-2H^{2}(t)\left(1-\frac{1}{6}\left(\frac{k}{a(t)H(t)}\right)^{2}\right)(-k^{2})\tilde{\delta}_{k}^{(2)}(t)$
$\displaystyle=-\frac{8}{3}\frac{k^{2}}{a(t)^{2}}a(t)^{2}\,8\pi
G\,\rho_{\text{GW}}(k/a(t)).$ (124)

and the factor $-k^{2}$ can be canceled such that

$\displaystyle\ddot{\tilde{\delta}}_{k}^{(2)}(t)+H(t)\dot{\tilde{\delta}}_{k}^{(2)}(t)-2H^{2}(t)\left(1-\frac{1}{6}\left(\frac{k}{a(t)H(t)}\right)^{2}\right)\tilde{\delta}_{k}^{(2)}(t)$
$\displaystyle=\frac{8}{3}\,8\pi G\,\rho_{\text{GW}}(k/a(t),t).$ (125)

For sub-horizon modes we can neglect the unity on the left hand side of the
equation. Additionally, the right hand side can be formulated in terms of
standard abundance $\Omega_{\text{GW}}$ by replacing
$\rho_{\text{GW}}=\rho_{\text{tot}}\Omega_{\text{GW}}$ and using
$\rho_{\text{tot}}=\frac{3H_{*}^{2}}{8\pi G}$ at the time of the phase
transition. We assume that the generation of GWs coincides with the duration
of the FPT and thus completes within less than a Hubble time. Hence we can
neglect the friction term $H(t)\dot{\tilde{\delta}}_{k}^{(2)}(t)$, approximate
$a(t)\approx a(t_{*})\approx a(t_{*}+1/\beta)$ and $H(t)\approx
H(t_{*})\approx H(t_{*}+1/\beta)$ and use $\rho_{\text{GW}}$ from the previous
subsection. After the phase transition completes the energy density of GWs
simply redshifts as a radiation and we assume that during that time its power
as source is negligible. In total our result reads

$\displaystyle\ddot{\tilde{\delta}}_{k}^{(2)}(t)+\frac{1}{3}\frac{k^{2}c^{2}}{a(t_{*})^{2}}\tilde{\delta}_{k}^{(2)}(t)=8H_{*}^{2}\Omega_{\text{GW}}(k/a(t_{*}),t)\quad\text{for}\>t\in[t_{*},t_{*}+1/\beta],$
(126)
$\displaystyle\ddot{\tilde{\delta}}_{k}^{(2)}(t)+H(t)\dot{\tilde{\delta}}_{k}^{(2)}(t)+\frac{1}{3}\frac{k^{2}c^{2}}{a(t)^{2}}\tilde{\delta}_{k}^{(2)}(t)=8H(t)^{2}\left(\frac{a_{*}}{a(t)}\right)^{4}\Omega_{\text{GW}}(k/a(t),t)\approx
0\quad\text{for}\>t>t_{*}+1/\beta.$ (127)

Note that the right hand side of Eq. (127) decays as $H(t)^{2}/a(t)^{4}$ and
thus can be safely neglected. We have checked this approximation semi-
analytically and found it to be consistent, see Appendix C. Also note that the
choice of gauge is of negligible importance for sub-horizon modes.

#### Solving the equation:

Next, we solve Eq. (126) to find $\tilde{\delta}^{(2)}_{k}(t)$ and estimate
its impact on the matter power spectrum. This requires us to calculate the
energy density in GWs from the fractional, logarithmic energy density
$\Omega^{\text{log}}_{\text{GW}}$ in Eq. (109). Integrating the equation gives

$\displaystyle\Omega_{\text{GW}}(k,t):=\frac{1}{\rho_{\text{tot}}}\rho_{\text{GW}}(k,t)=\kappa_{\text{eff}}^{2}\left(\frac{H_{*}}{\beta}\right)^{2}\left(\frac{\alpha}{1+\alpha}\right)^{2}\int_{k_{\text{min}}}^{k}\Delta(k^{\prime}/\beta,t,v_{w})\text{d}\ln
k^{\prime},$ (128)

where we take for $k_{\text{min}}=\frac{a_{*}H_{*}}{c}$ the inverse size of
the horizon at transition time and $v_{w}=c$ for the bubble wall velocity. The
resulting energy density of GWs as a function of the wave number is shown in
Fig. 5 for the dimensionless time $\tau:=H_{*}t$ and wave number
$\kappa:=\frac{c}{a_{*}H_{*}}k$.

Figure 5: The integrated dimensionless power spectrum of GWs sourced by a FPT.
The spectrum is evaluated in the middle of the phase transition
$\tau_{\text{ev}}$. The yellow line shows the horizon mode and the pink dashed
dotted line indicates the peak wave number of the logarithmic GW abundance. We
chose $\beta/H_{*}=1$ for demonstration reasons.

With this scaling the differential equation becomes

$\displaystyle\tilde{\delta}^{(2)\,\prime\prime}(\kappa,\tau)+\frac{1}{3}\kappa^{2}\tilde{\delta}^{(2)}(\kappa,\tau)=8\cdot\Omega_{\text{GW}}(\kappa,\tau),$
(129)

where primes denote the derivative with respect to unit free time
$\tau:=H_{*}\cdot t$. The numerical solution at the end of the phase
transition $\tilde{\delta}^{(2)}(\kappa,\tau_{*}+\frac{H_{*}}{\beta})$ is
shown in Fig. 6 for initial conditions
$\tilde{\delta}^{(2)}(\kappa,\tau_{*})=\tilde{\delta}^{(2)\prime}(\kappa,\tau_{*})=0$.

Figure 6: Left: In blue the numerical solution of Eq. (129) at the end of the
phase transition and with initial conditions
$\tilde{\delta}^{(2)}(\kappa,t)=\tilde{\delta}^{(2)\,\prime}(\kappa,t)=0$ and
for $\alpha\to\infty$ and $\beta/H_{*}=1$. We also show the source term in
orange at the end of the phase transition
$\tau_{\text{ev}}=\tau_{*}+H_{*}/\beta$. Right: Induced density perturbation
by GWs from FPT. For demonstration purposes we chose
$\alpha\to\infty,\beta=H_{*}$. Shown are solutions for different scales
$a_{*}H_{*}$ in comparison with the horizon at matter-radiation equality.

For an analytical resolution of Eq. (129) in various simplifying limits see
appendix D. How to interpret this equation? From the expansion of the pressure
to second order we see that for adiabatic perturbations and small changes in
the sound speed on sub-horizon scales, that

$\displaystyle p^{(2)}=c_{s}^{2}\rho^{(2)}+\sigma s^{(2)}+\frac{\partial
c_{s}^{2}}{\partial\epsilon}\rho^{(1)\,2}+\frac{\partial c_{s}^{2}}{\partial
s}s^{(1)\,2}+\frac{\partial\sigma}{\partial s}s^{(1)\,2}\approx
c_{s}^{2}\rho^{(2)},$ (130)

with $\sigma:=\left(\partial p/\partial s\right)$ [62]. Therefore, as for
linear perturbations, photon perturbations are characterized by
$c_{s}^{2}=1/3$ and hence Eq. (129) describes the evolution of photon
perturbations $\tilde{\delta}^{(2)}\equiv\tilde{\delta}^{(2)}_{\gamma}$.
Comparing this equation with the wave equation for photon perturbations in the
photon-baryon fluid before photon decoupling [63, 64]

$\displaystyle\ddot{\delta}_{\gamma}+c_{s}^{2}\frac{k^{2}}{a^{2}}\delta_{\gamma}=\frac{4}{3}4\pi
G\left(\rho^{(0)}_{D}\delta_{D}+\rho^{(0)}_{B}\delta_{B}+\rho^{(0)}_{\gamma}\delta_{\gamma}\right),$
(131)
$\displaystyle\delta^{\prime\prime}_{\gamma}+c_{s}^{2}\kappa^{2}\delta_{\gamma}=2\left(\Omega^{(0)}_{D}\delta_{D}+\Omega^{(0)}_{B}\delta_{B}+\Omega^{(0)}_{\gamma}\delta_{\gamma}\right),$
(132)

we notice, that what we found is a very similar system. But in our case their
oscillations are driven by the gravitational wave density instead of matter or
radiation density component. Since at that time baryons are still tightly
coupled to photons they follow almost the same wave equation and thus we
interpret our findings as baryon acoustic oscillations (BAOs) at a second
order perturbative level driven by the GW energy density. As seen in Fig. 6
the oscillations lie on top of a dominant peak. The typical sound horizon of
the oscillations is given by

$\displaystyle
r^{\text{GW}}_{s}:=\int_{t_{*}}^{t_{*}+1/\beta}\frac{\text{d}t}{a_{*}}c_{s}=\frac{1}{\sqrt{3}a_{*}\beta}\frac{2\pi}{\sqrt{3}\,k_{\text{peak}}}.$
(133)

For comparison, the typical sound horizon for standard BAOs and our BAOs is

$\displaystyle
r_{s}=147\,\text{Mpc}\>\cite[cite]{[\@@bibref{}{BOSS:2016wmc,Aghanim:2018eyx}{}{}]},$
(134) $\displaystyle r^{\text{GW}}_{s}=3\,\text{Mpc},$ (135)

respectively. Like for the standard BAOs after photon decoupling the baryons
will transfer this information gravitationally to the dark matter
perturbations and will thus be imprinted in the matter power spectrum.

Let us estimate the time on which a FPT has to occur in order to impact the
matter power spectrum by density fluctuations produced via Eq. (129). The
typical comoving scale on which the GW energy density per logarithmic
frequency $\Omega^{\text{log}}_{\text{GW}}(k,t)$ peaks at the end of the
transition is at

$\displaystyle\frac{k_{\text{peak}}}{a_{*}}\approx\frac{2\pi\beta}{c},$ (136)

with the phase transition duration $1/\beta$. Around this scale, the source
term, the fractional energy density $\Omega_{\text{GW}}$, becomes
approximately constant (see Fig. 5). Hence we can use it as a typical scale
which will also be inherited to the induced density perturbations via Eq.
(129). We rewrite the phase transition duration in terms of the Hubble
parameter $\beta=r_{\beta}H_{*}=r_{\beta}H(t_{*})$, where $r_{\beta}>1$ for
transitions shorter than a Hubble time. Therefore, the comoving wave number
where the density fluctuation spectrum is approximately maximal is

$\displaystyle k_{\text{peak}}=\frac{2\pi
r_{\beta}H_{*}a_{*}}{c},\quad\text{or}$
$\displaystyle\kappa_{\text{peak}}=2\pi r_{\beta}.$ (137)

This is analogous to a primordial density fluctuation which enters the horizon
at $H_{*}a_{*}$, only that in our case we can shift the scale relative to
$H_{*}a_{*}$ by the duration ratio $r_{\beta}$ of the phase transition.

From then on the scale of the density fluctuation is fixed and the time of a
phase transition that impacts the matter power spectrum at its typical scale
must fulfil the condition

$\displaystyle\frac{2\pi r_{\beta}H_{*}a_{*}}{c}\gtrsim k_{\text{eq}}.$ (138)

This condition is met by late phase transitions around

$\displaystyle t:10^{6}\,\text{s}-t_{\text{eq}}\sim
T:(\mathcal{O}(100)-\mathcal{O}(1))\,\text{eV}.$ (139)

We calculate the Hubble rate for these times using

$\displaystyle
H(t)=\frac{\dot{a}}{a}=H_{0}\sqrt{\Omega_{m0}}\frac{\sqrt{a+a_{\text{eq}}}}{a^{2}},$
(140)

where $H_{0}\approx 70\,\text{Mpc}/(\text{km}\,\text{s})\approx 2.27\cdot
10^{-18}\text{s}^{-1}$ denotes the Hubble rate today and
$a_{\text{eq}}=\frac{\Omega_{\text{rad}0}}{\Omega_{m0}}=\frac{8.5\cdot
10^{-5}}{0.3}=2.4\cdot 10^{-4}$ is the scale factor at equality (the Hubble
rate is $H_{\text{eq}}=H(a_{\text{eq}})=9.1\cdot 10^{-14}\frac{1}{\text{s}}$
or $t_{\text{eq}}\approx 70000\,\text{years}$). Integrating this equation
leads to an implicit equation for the scale factor

$\displaystyle t\cdot
H_{0}=\frac{2}{3}\frac{1}{\sqrt{\Omega_{m,0}}}\left[\sqrt{a+a_{\text{eq}}}(a-2a_{\text{eq}})+2a_{\text{eq}}^{3/2}\right].$
(141)

For events sufficiently far enough from equality $a\ll a_{\text{eq}}$ we can
approximate Eq. (141) and get as limiting equation for the scale factor
$a(t)=\sqrt{3\cdot H_{0}\sqrt{\Omega_{\text{rad},0}}\cdot t}$.

#### Impact on matter power spectrum:

Due to the production of extra deviations $\tilde{\delta}^{(2)}$ from the
energy density by the phase transition the primordial modes around $k_{*}=2\pi
a_{*}H_{*}/c$ experience a modification compared to their standard evolution.
The change is captured by the transfer function

$\displaystyle
T^{2}_{\tilde{\delta}^{(2)}}(k):=1+\left(\frac{\tilde{\delta}^{(2)}(k)}{\delta_{*}^{(1)}(k)}\right)^{2},$
(142)

where $\delta_{*}^{(1)}(k)$ are the primordial perturbations inside the
horizon at $t_{*}$. Then the altered matter power spectrum with the amplitude
at matter radiation equality compared to the spectrum today
$\mathcal{P}_{0}(k)$ is

$\displaystyle\tilde{\mathcal{P}}_{\text{eq}}(k)=T_{\tilde{\delta}^{(2)}}^{2}(k)\mathcal{P}_{0}(k)D_{+}^{2}(a_{\text{eq}}),$
(143)

with the approximate linear growth function

$\displaystyle
D_{+}(a)\approx\frac{5}{2}\frac{a\Omega_{m0}}{\Omega_{m0}^{3/4}-\Omega_{\Lambda}+(1+\Omega_{m0}/2)(1+\Omega_{\Lambda}/70)},$
(144)

which is $D_{+}(a_{\text{eq}})\approx 2.5\cdot 10^{-4}$ around equality. The
linear matter power spectrum linearly extrapolated to today is given by the
fitting formula [67]

$\displaystyle\mathcal{P}_{0}(k)=A_{0}\,k\cdot\frac{\ln(1+c_{1}q)}{c_{1}q}\cdot\left(1+(c_{2}q)+(c_{3}q)^{2}+(c_{4}q)^{3}+(c_{5}q)^{4}\right)^{-\frac{1}{4}},$
(145)

with

$\displaystyle
q:=\frac{k}{\Omega_{m0}h\cdot\mathrm{exp}\left(-\Omega_{\text{baryon}0}-\sqrt{2h}\cdot\frac{\Omega_{\text{baryon}0}}{\Omega_{m0}}\right)}\approx
0.073\frac{k}{k_{\text{eq}}},$ (146)

and $c_{1}=2.34,c_{2}=3.89,c_{3}=16.1,c_{4}=5.46$ and $c_{5}=6.71$. The
reduced Hubble parameter is set to $h=0.7$ and the abundance of baryons today
is $\Omega_{\text{baryon}0}=0.05$ and $\Omega_{m0}=0.3$ [66]. The amplitude
$A_{0}$ of $\mathcal{P}_{0}(k)$ is calibrated such that the variance becomes
[66]

$\displaystyle
0.8^{2}=\sigma^{2}_{8}=\int_{0}^{\infty}\text{d}\tilde{k}\frac{\tilde{k}^{2}}{2\pi^{2}}\mathcal{P}_{0}(\tilde{k})\times\left(\frac{3j_{1}(\tilde{k}R)}{(\tilde{k}R)}\right)^{2}.$
(147)

Here $j_{1}(x):=\sin(x)/x^{2}-\cos(x)/x$ and $R=8\,\text{Mpc}/h$.

In order to derive the transfer function $T_{\tilde{\delta}^{(2)}}(k)$ we
estimate the amplitude of a typical density perturbation $\delta_{*}^{(1)}(k)$
for sub-horizon modes as standard deviation from $\mathcal{P}$

$\displaystyle\delta_{*}^{(1)}(k)\cong
D_{+}(a_{\text{eq}})\sqrt{\int_{0}^{\infty}\text{d}\tilde{k}\frac{\tilde{k}^{2}}{2\pi^{2}}\mathcal{P}_{0}(\tilde{k})W_{k}^{2}(\tilde{k})}\quad\text{for}\quad
k_{*}\leq k,$ (148)

where
$W_{k}(\tilde{k})=\frac{3}{(\tilde{k}/k)^{3}}(\sin(\tilde{k}/k)-\tilde{k}/k\cos(\tilde{k}/k))$
is called window function. The restriction to modes with $k_{*}<k$ is
necessary since only modes that have entered the horizon at the time of the
phase transition are relevant. In Fig. 7 we show Eq. (148) in terms of the
dimensionless wave number $\kappa$, $\delta_{*}^{(1)}(\kappa\cdot k_{*})$ for
$10^{0}\leq\kappa$, at different transition times.

Figure 7: The estimated first order density fluctuations from the matter power
spectrum in Eq. (148) as a function of dimensionless wave number for different
times. Here $k_{*}=\frac{a*H_{*}}{c}$. The $\kappa$-axis and the Hubble
horizon shown in yellow are given in terms of the respective time for each of
the curves.

We use the estimated primordial density fluctuations to define the transfer
function Eq. (142) which is show in Fig. 9 for some example cases together
with the modified matter power spectrum $\tilde{\mathcal{P}}$ at equality.
Note that $\delta^{(2)}\ll 1$ is fulfilled at all times, see for example Fig.
7. The whole procedure is schematically summarized in Fig. 8.

Figure 8: Schematic summary of the deduction of the matter power spectrum with
the effect from gravitational wave induced second order density fluctuations.

As seen in Figs. 9, 10 and 11 the GWs produced by the FPT imprint a peak on
the matter power spectrum around the comoving scale $k_{\text{peak}}$. The
transfer functions decrease rapidly with smaller phase transition duration
$r_{\beta}$ and also with smaller strength $\alpha$. This behaviour is
expected from the prefactors of the GW energy density in Eq. (128).

Hence, the height of the peak is determined by the parameters $t_{*}$,
$\alpha$ and $r_{\beta}$. We can put limits on them by requiring that the
height of the peak should not exceed the bound set by the cosmic variance of
the linear matter power spectrum. The latter is defined via [68]

$\displaystyle\sigma(k):=\sqrt{\text{cov}(\mathcal{P}_{0}(k),\mathcal{P}_{0}(k))}=\mathcal{P}_{0}(k)\sqrt{\left(\frac{2}{N}+\frac{1}{n}\right)},$
(149)

which holds for Gaussian random fields. $N$ denotes the number of modes and
$n$ is related to the so called band-averaged trispectrum which can be
estimated to $1/n\approx 0.0079^{2}\,(\text{Gpc}/h)^{3}/V$ [69]. We estimate
the number of modes as $2/N=(2\pi)^{2}/(V\cdot k^{2}\Delta k)$ with $V\approx
1(\text{Gpc}/h)^{3}$ and $\Delta k\sim
0.02\cdot\log{\left(k\,\text{Mpc}\right)}\,(\text{Mpc}/h)^{-1}$ (typical
distance between galaxies) which reproduces approximately the cosmic variance
found in [69].

The modified matter power spectrum $\tilde{\mathcal{{P}}}$ should not exceed
this bound, i.e

$\displaystyle\tilde{\mathcal{{P}}}(k)\ll\sigma(k),\quad\forall k.$ (150)

In Fig. 12 we show a parameter scan in the $\alpha$-$\beta$-plane for
different phase transition times $t_{*}$. The red shaded regions are excluded
by the cosmic variance bound, while values in the white region are consistent
with it. We observe that only very long $r_{\beta}<5-6.8$ and strong
$\alpha>1$ phase transitions can be ruled out.

The earlier the phase transition takes place the less is it constrained. FPTs
with such extreme parameter values have been proposed in the past. Long
lasting transitions are realized for example in SUSY [70] and models with a
lot of supercooling are for example Randall-Sundrum, composite Higgs models
[51] and models with an almost conformal symmetry in general [15].

,

Figure 9: The impact of a late phase transition for different starting times
$t_{*}$ on the linear matter power spectrum. We chose for each time an inverse
duration of $\beta=3H_{*}$ and latent heat $\alpha=3$, respectively.

Figure 10: The impact of a late phase transition for different strength
$\alpha$ on the linear matter power spectrum. We chose for each $\alpha$ an
inverse duration of $\beta=3H_{*}$ and fixed the transition time to
$t_{*}=5\times 10^{9}\,\text{s}$, respectively.

Figure 11: The impact of a late phase transition for different duration ratios
$r_{\beta}$ on the linear matter power spectrum. We chose for each $r_{\beta}$
a strength $\alpha=3$ and fixed the transition time to $t_{*}=5\times
10^{9}\,\text{s}$, respectively.

In Fig. 13 we convert contour line values into a bound on the GW signal today
in the standard frequency - GW abundance plane. The logarithmic GW abundance
today due to bubble collisions is [43]

$\displaystyle h^{2}\Omega^{\text{log}}_{\text{GW}\,0}(f)=1.67\times
10^{-5}\cdot
r_{\beta}^{-2}\cdot\left(\frac{\kappa_{\text{eff}}(\alpha)\alpha}{1+\alpha}\right)^{2}\left(\frac{100}{g_{*}}\right)^{\frac{1}{3}}\left(\frac{0.11\,v_{w}^{3}}{0.42+v_{w}^{2}}\right)\frac{3.8(ff_{\text{peak}})^{2.8}}{1+2.8(f/f_{\text{peak}})^{3.8}},$
(151)

with the peak frequency $f_{\text{peak}}$ today

$\displaystyle f_{\text{peak}}=16.5\times
10^{-6}\,\text{Hz}\frac{0.62}{v_{w}^{2}-0.1v_{w}+1.8}r_{\beta}\left(\frac{T_{*}}{100}\right)\left(\frac{g_{*}}{100}\right)^{\frac{1}{6}}.$
(152)

The time of the phase transition can be converted into the temperature of the
plasma by

$\displaystyle T_{*}=\frac{30}{\pi^{2}}\frac{3}{4g_{*}^{3}}\left(\frac{1}{8\pi
Gt_{*}^{2}}\right)^{\frac{1}{4}}$ (153)

and the number of relativistic degrees of freedom $g_{*}$ after the QCD phase
transition (and hence for a late phase transition) is $3.36$. Note, that for
BSM models this value differs depending on the field content and their
properties. The frequency window is set to $f\approx 1.5\times
10^{-16}\,\text{Hz}$ as lower bound and $f\approx 1.5\times
10^{-14}\,\text{Hz}$ as upper bound which approximately corresponds to the
right, big $k$ slope of the matter power spectrum. For smaller frequencies our
assumption of a radiation dominated universe becomes very weak.

For comparison, we show also the projected bounds for LISA [2], the timing
pulsar arrays [71] NANOGrav [72, 73], PPTA [74], EPTA [75] and CMB [76, 77,
78].

Following reference [79, 80, 16] the GW energy density can be limited by the
effective number of neutrino species $N_{\nu}$ via

$\displaystyle h^{2}\Omega_{\text{GW}}(f)\leq 5.6\cdot 10^{-6}\Delta N_{\nu},$
(154)

where $\Delta N_{\nu}$ denotes the deviation from the SM value $N_{\nu}=3$.
BBN constrains this number to $\Delta N_{\nu}\leq 0.2$ [81] giving the bound
on the allowed amount of GW before BBN shown in Fig. 13. The indirect bound
for the CMB is taken from [76].

Figure 12: Parameter scan of $\alpha$ and $r_{\beta}$ for different transition
times $t_{*}$. The red colored regions are excluded by cosmic variance and
thus a late FPT with these parameter combination would change structure
formation too strongly. Figure 13: The GW abundance today as a function of
today frequency. Shown are the projected bound by LISA [2] (dashed-dot blue),
the bound set by PTAs [72, 73, 74, 75] (solid blue) and the indirect bounds
from CMB and BBN [76, 77, 78] (hatched areas). As the result of our
calculation, the red colored regions are excluded by structure formation.
However, many assumptions went into our calculation so the bound should be
taken with care.

##  6 Conclusion

Let us summarize the results. In this work we have studied the possible impact
of a FPT on small scale structure via the production of GWs in the radiation
dominated epoch. A linear relation between the energy density of GWs
$\rho_{\text{GW}}(k)$ and adiabatic density perturbations has been found by
expanding the full non-linear Eqs. (30) and (31) to second order in the $1+3$
covariant formulation. In this formalism the spacetime is decomposed into the
direction of fluid flow and its orthogonal hyper-surface. Then, a set of gauge
invariants to first order with clear geometrical interpretations can be
constructed.

When only considering parameters for which the GW energy density surpasses the
other source terms during the transition, the adiabatic density perturbations
follow a wave equation which is driven by the GW energy density. In this case
our equation describes photon acoustic oscillations induced by the GW energy
density. Since the photons are still coupled to the baryons at such times, the
baryons undergo the same oscillations which manifest themselves eventually in
the matter power spectrum.

Since phase transitions are typically taking place within the Hubble horizon
$H_{*}$ at the time of the transition the scale on which the perturbations are
affected is bounded by the horizon size $k_{*}=a_{*}H_{*}/c$. However, we
found that the scale that is maximally impacted equals the scale where the GW
energy density per logarithmic frequency has a maximum $k_{*}=2\pi
r_{\beta}\,a_{*}H_{*}/c$. This implies that the linear matter power spectrum,
if at all, can only be affected on the length scales of galaxies and above if
the transition occurred at very late times $\geq 10^{6}\,\text{s}$ but still
within the radiation dominated regime. Late phase transitions and their impact
on structure formation (also due to gravitational waves) have been discussed
in the past, for example in the matter dominated era [82, 83, 84] (in the
literature the phrase late time phase transition is sometimes used for
transitions after equality or photon decoupling). Specific particle models
have been discussed in [85] and a model with a very late phase transition
including a dark energy component is presented in [86].

The maximally allowed duration $\beta^{-1}$ and strength $\alpha$ of a phase
transition is bounded by cosmic variance and depends on the time of the
transition. We find that this bound constrains these parameters only very
weakly, excluding transitions that last longer than $\gtrsim 1/6.8$ Hubble
times in the case the transition is close to equality and $\gtrsim 1/5$ in the
case the transition takes place on galaxy scales. From the parameter set
$t_{*},\beta$ and $\alpha$ we derived the GW abundance per logarithmic
frequency today and translated the bounds from structure formation into an
exclusion region in Fig. 13.

Our results are based on the following assumptions. First of all we looked at
adiabatic perturbations only. We simplified our calculation further by
neglecting anisotropy and vorticity effects as well as current density
effects. In principle the anisotropic stress could be also have effects on the
matter power spectrum directly. As a next step it would be reasonable to study
the possible effects of the anisotropic stress on the density perturbations in
more details. For example, its scalar part (corresponding to the quadruple
term in the momentum distribution caused by the bubble collision in the fluid)
could constitute a difference in the Bardeen potentials $\Psi-\Phi\sim\Pi$
analogous to neutrino and photon anisotropies and in this way even affect
linear perturbations. The effect of an extra anisotropic stress on the CMB and
on curvature perturbations has been discussed in [87] also using the $1+3$
covariant formalism. Additionally, one could consider effects of the
anisotropic stress on a second perturbative level. A non zero and transverse
anisotropic stress tensor can appear in the non-linear Eqs. (30), (31) and the
conservation laws Eqs. (20), (21) coupled to the acceleration $A_{a}$ and the
shear $\sigma_{ab}$. The acceleration $A_{a}^{(2)}$ is thus not parallel to
the density gradient any more which will make the calculation much more
complex when including the anisotropic stress.

In our derivation we assumed the equation of state parameter $\omega$ and the
sound speed $c_{s}^{2}$ are constant in time and space and also that
$\frac{\delta\rho}{\delta p}=\frac{\delta^{2}\rho}{\delta^{2}p}=c_{s}^{2}$. In
general these parameters could depend on space and time. However, on the one
hand the decline of $\omega$ close to equality is very gentle and on the other
hand the change within a Hubble time is expected to be negligible. As closer
we get to matter-radiation equality $\omega$ departs more and more from being
$1/3$. A rough estimation gives $\omega\approx 0.27$ at $t\approx
10^{11}\,\text{s}$.

In this work we found that only strong GWs sourced by phase transitions with a
lot of supercooling can have effects on structure. In this regime the bubble
dynamics is fixed to bubbles expanding into vacuum and hence the only source
of GWs are bubble collisions.

Note also that our study is limited to phase transitions on sub-horizon scales
which complete within a Hubble time $\beta>H_{*}$. Our results are very close
to this boundary and hence effects of the Hubble friction terms in the wave
equation for the GWs and the density perturbations might suppress the
amplitudes even further, shrinking the constrained region in parameter space.

In future work we will look at more direct consequences of the phase
transition on structure formation. One idea is to take up on the work done
Schmid et al. [24, 25] and study the effect on linear perturbations by
changing sound speed. As mentioned, in [26] it was shown that the sound speed
does not change a lot in particle models with many scalar fields, but could
depart from $1/\sqrt{3}$ in fermion rich models. Another possibility is to
study the direct impact of the anisotropic stress on linear perturbations, as
mentioned before through the difference in the Bardeen potentials
$\Psi-\Phi\sim\Pi$. Also, the huge amount of supercooling $\alpha\gg 1$ in our
calculation turns the background cosmology from radiation dominated to vacuum
energy dominated such that the equation of motion of the linear density
perturbations changes which could also lead to direct effect on the matter
power spectrum.

## Acknowledgements

The authors thank Ruth Durrer for advice in the initial stages of this project
and for helpful comments on the final version of the paper. CD would like to
acknowledge insightful discussions with Andreas Trautner. S.C.C. would like to
thank the Max-Planck-Institut für Kernphysik in Heidelberg for hospitality
during his visit, where this work was initiated. The work of S.C.C. is
supported by the Spanish grants SEV-2014-0398, FPA2017-85216-P (AEI/FEDER,
UE), PROMETEO/2018/165 (Generalitat Valenciana) and BES-2016-076643. This work
is supported by the Deutsche Forschungsgemeinschaft (DFG, German Research
Foundation) under Germany’s Excellence Strategy EXC 2181/1 - 390900948 (the
Heidelberg STRUCTURES Excellence Cluster).

## APPENDIX

##  Appendix A Important identities in 1+3 covariant theory

The orthogonal projected gradient and the time derivative of the orthogonal
projection operator $h_{ab}$ meet the relations [27, 30]

$\displaystyle{\text{D}}_{a}h_{bc}=0,$ (A.1)
$\displaystyle{\text{D}}^{a}h_{ab}=u_{b}\Theta,$ (A.2)
$\displaystyle\dot{h}_{ab}=u_{b}A_{a}+u_{a}A_{b}.$ (A.3)

Calculating the projected gradient of the four velocity gives

$\displaystyle{\text{D}}_{b}u_{a}=\sigma_{ab}+\omega_{ab}+\frac{1}{3}\Theta
h_{ab}.$ (A.4)

Also important are the commutation laws for the derivatives which we simply
repeat from reference [30]. For a scalar $f$, a vector $v_{a}$ and a tensor
$S_{ab}$ we have for the spatial derivative

$\displaystyle{\text{D}}_{[a}{\text{D}}_{b]}f=-\omega_{ab}\dot{f},$ (A.5)
$\displaystyle{\text{D}}_{[a}{\text{D}}_{b]}v_{c}=-\omega_{ab}\dot{v}_{\langle
a\rangle}+\frac{1}{2}\mathcal{R}_{dcba}v^{d},$ (A.6)
$\displaystyle{\text{D}}_{[a}{\text{D}}_{b]}S_{cd}=-\omega_{ab}h_{c}{}^{e}h_{d}{}^{f}\dot{S}_{ef}+\frac{1}{2}(\mathcal{R}_{ecba}S^{e}{}_{d}+\mathcal{R}_{edba}S_{c}{}^{e}),$
(A.7)

where $\mathcal{R}_{abcd}$ is the Riemann tensor in the local rest space of
the observer.

Similarly the time derivative and the space derivative do not commute in
general

$\displaystyle{\text{D}}_{a}\dot{f}-h_{a}{}^{b}\dot{({\text{D}}_{b}f)}=-\dot{f}A_{a}+\frac{1}{3}\Theta{\text{D}}_{a}f+{\text{D}}_{b}f\left(\sigma^{b}{}_{a}+\omega^{b}{}_{a}\right).$
(A.8)

##  Appendix B Harmonic Decomposition

It is convenient to expand all scalars, vectors and tensors in harmonic
functions. No matter if scalar harmonic $\mathcal{Q}_{k}$, vector harmonic
$\mathcal{Q}_{k,a}$ or tensor harmonic $\mathcal{Q}_{k,ab}$, their defining
property is to be an eigenfunction of the orthogonal projected Laplace
operator (Laplace-Beltrami equation)

$\displaystyle{\text{D}}^{2}\mathcal{Q}_{k,\\{\hphantom{a},a,ab\\}}=-\frac{k^{2}}{a^{2}}\mathcal{Q}_{k,\\{\hphantom{a},a,ab\\}},$
(B.1)

with eigenvalue $-k^{2}/a^{2}$. In case of a flat space $K=0$ the orthogonal
projected Laplace operator ${\text{D}}^{2}$ reduces to the usual Laplace
operator $\nabla^{2}/a^{2}$ such that the harmonic functions $\mathcal{Q}$ are
Fourier transforms [28, 88]. Scalar, vector and tensor modes thus transform
like

$\displaystyle
f(\textbf{x},t)=\int\text{d}\textbf{k}f_{\textbf{k}}(t)\mathcal{Q}_{\textbf{k}},$
$\displaystyle
V_{a}^{\bot}(\textbf{x},t)=\int\text{d}\textbf{k}\sum_{m=-1,1}V_{\textbf{k}}^{\bot\,[m]}\mathcal{Q}_{\textbf{k},a}^{[m]},$
(B.2) $\displaystyle
S_{ab}^{T}(\textbf{x},t)=\int\text{d}\textbf{k}\sum_{m=-2,2}S_{\textbf{k}}^{T\,[m]}\mathcal{Q}_{\textbf{k},ab}^{[m]},$
(B.3)

with

$\displaystyle\mathcal{Q}_{\textbf{k}}=\exp{(i\textbf{k}\cdot\textbf{x})},$
$\displaystyle\mathcal{Q}_{\textbf{k},a}^{[\pm
1]}=\frac{-i}{\sqrt{2}}(\textbf{e}_{1}\pm
i\textbf{e}_{2})_{a}\exp{(i\textbf{k}\cdot\textbf{x})},$
$\displaystyle\mathcal{Q}_{\textbf{k},ab}^{[\pm
2]}=-\sqrt{\frac{3}{8}}(\textbf{e}_{1}\pm
i\textbf{e}_{2})_{a}(\textbf{e}_{1}\pm
i\textbf{e}_{2})_{b}\exp{(i\textbf{k}\cdot\textbf{x})},$

where $\textbf{e}_{1}$ and $\textbf{e}_{2}$ are orthonormal basis vectors. The
scalar part of a vector and the scalar- and vector part of a tensor expand
like

$\displaystyle
V_{a}^{||}=\int\text{d}\textbf{k}\left(-i\frac{k_{a}}{k}\right)V_{\textbf{k}}\mathcal{Q}_{\textbf{k}},$
(B.4) $\displaystyle
S_{ab}^{||}=\int\text{d}\textbf{k}\left(-\frac{k_{a}k_{b}}{k^{2}}+\frac{1}{3}h_{ab}\right)S_{\textbf{k}}\mathcal{Q}_{\textbf{k}},$
(B.5) $\displaystyle
S_{ab}^{\bot}=\int\text{d}\textbf{k}\sum_{m=-1,+1}\left(-\frac{i}{2k}\right)(k_{a}S_{\textbf{k}}^{\bot\,[m]}\mathcal{Q}_{\textbf{k},b}^{[m]}+k_{b}S_{\textbf{k}}^{\bot\,[m]}\mathcal{Q}_{\textbf{k},a}^{[m]}),$
(B.6)

respectively.

##  Appendix C Decay of the source after the FoPT

We will now show that the GW source decays sufficiently fast after the PT such
that we can take the right hand side of Eq. (127) to zero. We first perform
the adimensional change of variables

$\displaystyle\tilde{\delta}^{\prime\prime}(\kappa,\tau)+\frac{a[\tau]^{\prime}}{a[\tau]}\tilde{\delta}^{\prime}(\kappa,\tau)+\frac{1}{3}\frac{a_{*}^{2}}{a[\tau]^{2}}\kappa^{2}\tilde{\delta}(\kappa,\tau)=8\frac{H[\tau]^{2}}{H_{*}^{2}}\frac{a_{*}^{4}}{a[\tau]^{4}}\cdot\Omega_{\text{GW}}(\kappa,\tau).$
(C.1)

Note that
$\Omega_{\text{GW}}(\kappa,\tau>\tau_{f})=\Omega_{\text{GW}}(\kappa,\tau_{f}):=\Omega_{\text{GW}}(\kappa)$
is a constant in time, where $\tau_{f}=\tau_{*}+1/r_{\beta}$ is the time when
the PT ends. We can now solve the homogeneous part of Eq. (C.1) for a general
$a[\tau]$:

$\displaystyle\tilde{\delta}_{h}(\kappa,\tau>\tau_{f})=C_{\kappa}\cos\left[\frac{a_{*}\kappa}{\sqrt{3}}\int^{\tau}_{\tau_{f}}\frac{1}{a[\tau^{\prime}]}d\tau^{\prime}\right]+D_{\kappa}\sin\left[\frac{a_{*}\kappa}{\sqrt{3}}\int^{\tau}_{\tau_{f}}\frac{1}{a[\tau^{\prime}]}d\tau^{\prime}\right],$
(C.2)

which in the radiation dominated universe simplifies to
$a(\tau)=a_{*}\sqrt{\frac{\tau}{\tau_{*}}}$ and $H(\tau)=H_{*}\frac{1}{2\tau}$

$\displaystyle\tilde{\delta}_{h}(\kappa,\tau>\tau_{f})=C_{\kappa}\cos\left[\frac{2\tau_{*}\kappa}{\sqrt{3}}\left(\sqrt{\tau}-\sqrt{\tau_{f}}\right)\right]+D_{\kappa}\sin\left[\frac{2\tau_{*}\kappa}{\sqrt{3}}\left(\sqrt{\tau}-\sqrt{\tau_{f}}\right)\right].$
(C.3)

This is the solution sourced solely by the GW energy during the PT. Turning
now to the solution sourced by the GW after the PT, using the variation of
parameters method we find

$\displaystyle\tilde{\delta}(\kappa,\tau>\tau_{f})=\tilde{\delta}_{h}(\kappa,\tau>\tau_{f})+\tilde{\delta}_{nh}(\kappa,\tau>\tau_{f}),$
(C.4)
$\displaystyle\tilde{\delta}_{nh}(\kappa,\tau>\tau_{f})=C_{\Omega}(\kappa,\tau)\cos\left[\frac{2\sqrt{\tau_{*}}\kappa}{\sqrt{3}}\left(\sqrt{\tau}-\sqrt{\tau_{f}}\right)\right]+D_{\Omega}(\kappa,\tau)\sin\left[\frac{2\sqrt{\tau_{*}}\kappa}{\sqrt{3}}\left(\sqrt{\tau}-\sqrt{\tau_{f}}\right)\right],$
(C.5) $\displaystyle
C_{\Omega}(\kappa,\tau)=-\int_{\tau_{f}}^{\tau}\frac{\sqrt{3}8\frac{H[\tau^{\prime}]^{2}}{H_{*}^{2}}\frac{a_{*}^{4}}{a[\tau^{\prime}]^{4}}\cdot\Omega_{\text{GW}}(\kappa)\sqrt{\frac{\tau^{\prime}}{\tau_{*}}}\sin\left(\frac{2\kappa\sqrt{\tau_{*}}\left(\sqrt{\tau^{\prime}}-\sqrt{\tau_{f}}\right)}{\sqrt{3}}\right)}{\kappa}\,d\tau^{\prime},$
(C.6) $\displaystyle
D_{\Omega}(\kappa,\tau)=\int_{\tau_{f}}^{\tau}\frac{\sqrt{3}8\frac{H[\tau^{\prime}]^{2}}{H_{*}^{2}}\frac{a_{*}^{4}}{a[\tau^{\prime}]^{4}}\cdot\Omega_{\text{GW}}(\kappa)\sqrt{\frac{\tau^{\prime}}{\tau_{*}}}\cos\left(\frac{2\kappa\sqrt{\tau_{*}}\left(\sqrt{\tau^{\prime}}-\sqrt{\tau_{f}}\right)}{\sqrt{3}}\right)}{\kappa}\,d\tau^{\prime}.$
(C.7)

For demonstration purposes we will now recombine the trigonometric functions
into a single one with a phase by using the identity

$C\cos x+D\sin x=A\sin(x+\phi),$ (C.8)

where the new amplitude is given by $A=\sqrt{C^{2}+D^{2}}$ and the relative
phase $\phi=\arctan\frac{C}{D}$. Applying this identity into Eqs. (C.3) and
(C.5) we obtain

$\displaystyle\tilde{\delta}_{h}(\kappa,\tau>\tau_{f})=A_{\kappa}\sin\left[\frac{2\tau_{*}\kappa}{\sqrt{3}}\left(\sqrt{\tau}-\sqrt{\tau_{f}}\right)+\phi_{\kappa}\right],$
(C.9)
$\displaystyle\tilde{\delta}_{nh}(\kappa,\tau>\tau_{f})=A_{\Omega}(\kappa,\tau)\sin\left[\frac{2\sqrt{\tau_{*}}\kappa}{\sqrt{3}}\left(\sqrt{\tau}-\sqrt{\tau_{f}}\right)+\phi_{\Omega}\right].$
(C.10)

We now compare the relative sizes of the amplitudes $A_{\kappa}$, the
amplitude of the solution sourced by the gravitational wave energy at $t_{f}$,
and $A_{\Omega}(\kappa,\tau_{\text{eq}})$, the amplitude of the solution
sourced by the gravitational wave energy after $t_{f}$ at $t_{\text{eq}}$. As
can be seen in Fig. (14), the homogeneous solution dominates and thus taking
the right hand side of Eq. (127) to zero is a sound approximation.

Figure 14: Ratio between the amplitudes in Eqs. (C.9) and (C.10) evaluated at
$t_{\text{eq}}$ for a benchamark strong PT with $\alpha\rightarrow\infty$,
$r_{\beta}=1$ and $t_{*}=10^{10}\,\text{s}$. Note that the ratio is always
bigger than 20, meaning that the solution sourced by the decaying GW source
after $t_{f}$ can be safely neglected in Eq. (127).

##  Appendix D Analytical solution of the GW sourced wave equation in the
small and high wave number limit

We start from the source given by Eq. (110). By performing the variable
tranformations $t=H_{*}\tau$ and $k=\frac{H_{*}a_{*}}{c}\kappa$ and by
slightly abusing the notation, we obtain

$\displaystyle
L(\tau)=\frac{4v\beta}{H_{*}^{2}}(\tau-\tau_{*})(\tau_{f}-\tau),$ (D.1)
$\displaystyle
f^{2}(\kappa,\tau)=L(\tau)^{2}\frac{v\epsilon}{\beta}\frac{1+\left(\frac{\kappa\frac{H_{*}a_{*}}{c}}{3}L(\tau)\right)^{2}}{1+\left(\frac{\kappa\frac{H_{*}a_{*}}{c}}{3}L(\tau)\right)^{2}+\left(\frac{\kappa\frac{H_{*}a_{*}}{c}}{3}L(\tau)\right)^{6}},$
(D.2)
$\displaystyle\Delta\left(\frac{\kappa}{a_{*}},\tau\right)=\left(\frac{\kappa
H_{*}}{c}\right)^{3}\left(\frac{\beta}{H_{*}}\right)^{2}(\tau-\tau_{*})^{2}f^{2}\left(\frac{\kappa}{a_{*}},\frac{\tau+\tau_{*}}{2}\right),$
(D.3)
$\displaystyle\Omega_{\text{GW}}^{\text{log}}\left(\frac{\kappa^{\prime}}{a_{*}},\tau\right)=\kappa_{\text{eff}}^{2}\left(\frac{H_{*}}{\beta}\right)^{2}\frac{\alpha^{2}}{(1+\alpha)^{2}}\Delta\left(\frac{\kappa}{a_{*}},\tau\right),$
(D.4)
$\displaystyle\Omega_{\text{GW}}\left(\frac{\kappa}{a_{*}},\tau\right)=\int_{\kappa_{\text{min}}}^{\kappa}\Omega_{\text{GW}}^{\text{log}}\left(\frac{\kappa^{\prime}}{a_{*}},\tau\right)d\log\kappa^{\prime}.$
(D.5)

Then, we want to expand the quantity
$\Omega_{\text{GW}}\left(\frac{\kappa}{a_{*}},\tau\right)$ around $\kappa=0$
and $\kappa\rightarrow\infty$. The expansion around $\kappa=0$ can simply be
done by first expanding $\Omega_{\text{GW}}^{\text{log}}$ and then
integrating:

$\displaystyle\Omega_{\text{GW}}^{\text{low}}=C_{\text{GW}}^{\text{low}}(\kappa-\kappa_{\text{min}})^{3}(\tau_{*}-\tau)^{4}(\tau+\tau_{*}-2\tau_{f})^{2}.$
(D.6)

With
$C_{\text{GW}}^{\text{low}}=\frac{1}{3}\frac{v^{3}}{c^{3}}\frac{\beta}{H_{*}}\frac{\alpha^{2}}{(1+\alpha)^{2}}\epsilon\kappa_{\text{eff}}$.
Note, however, that the same cannot be applied for the high $\kappa$ limit,
since the integration of $\Omega_{\text{GW}}^{\text{low}}$ will necessarily
run over small values of $\kappa^{\prime}$. In order to solve this issue, we
first obtain the value of $\kappa_{\text{cross}}$ defined as

$\displaystyle\left(\Omega_{\text{GW}}^{\text{log}}(\kappa_{\text{cross}}/a_{*},\tau_{f})\right)^{\text{low}}=\left(\Omega_{\text{GW}}^{\text{log}}(\kappa_{\text{cross}}/a_{*},\tau_{f})\right)^{\text{high}}\rightarrow\kappa_{\text{cross}}=3\frac{H_{*}}{\beta}\frac{c}{v}\frac{1}{(\tau_{f}-\tau_{*})^{2}},$
(D.7)

which is the value of $\kappa$ for which the low and high $\kappa$ limits of
the GW energy density intersect at $\tau=\tau_{f}$. Then for the high $\kappa$
approximation of $\Omega_{\text{GW}}$ we can write

$\displaystyle\Omega_{\text{GW}}^{\text{high}}$
$\displaystyle=\int_{\kappa_{\text{min}}}^{\kappa_{\text{cross}}}\left(\Omega_{\text{GW}}^{\text{log}}(\kappa/a_{*},\tau_{f})\right)^{\text{low}}d\log\kappa^{\prime}+\int_{\kappa_{\text{cross}}}^{\kappa}\left(\Omega_{\text{GW}}^{\text{log}}(\kappa/a_{*},\tau_{f})\right)^{\text{high}}d\log\kappa^{\prime}$
$\displaystyle=C_{\text{GW}}^{\text{low}}(\kappa_{\text{cross}}-\kappa_{\text{min}})^{3}(\tau_{*}-\tau)^{4}(\tau+\tau_{*}-2\tau_{f})^{2}+C_{\text{GW}}^{\text{high}}\frac{\kappa-\kappa_{\text{min}}}{\kappa\kappa_{\text{min}}}\frac{1}{(\tau+\tau_{*}-2\tau_{f})^{2}}$
(D.8)

and
$C_{\text{GW}}^{\text{high}}=81\frac{c}{v}\frac{H_{*}^{3}}{\beta^{3}}\frac{\alpha^{2}}{(1+\alpha)^{2}}\epsilon\kappa_{\text{eff}}$.

In these two regimes the differential equation Eq. (129) can be solved
analytically. Note, that the equation is only valid during the FPT, i.e.
$\tau_{*}<\tau<\tau_{f}$. The homogeneus part of the equation is given by
$\tilde{\delta}^{\prime\prime}(\kappa,\tau)+\frac{1}{3}\kappa^{2}\tilde{\delta}(\kappa,\tau)=0$
with a trivial solution:
$\tilde{\delta}_{h}(\kappa,\tau)=A_{\kappa}\cos\left(\frac{\kappa}{\sqrt{3}}\tau\right)+B_{\kappa}\sin\left(\frac{\kappa}{\sqrt{3}}\tau\right)$,
where $A_{\kappa}$ and $B_{\kappa}$ are given by the initial conditions. We
can then use the variation of parameters method to obtain the solution to the
non-homogeneous solution which is given by

$\displaystyle\tilde{\delta}(\kappa,\tau)=\tilde{\delta}_{h}(\kappa,\tau)+$
$\displaystyle\frac{8\sqrt{3}}{\kappa}\left[\sin\left(\frac{\kappa}{\sqrt{3}}\tau\right)\int_{0}^{\tau}\cos\left(\frac{\kappa}{\sqrt{3}}\tilde{\tau}\right)\Omega_{\text{GW}}(\tilde{\tau},\kappa)d\tilde{\tau}-\cos\left(\frac{\kappa}{\sqrt{3}}\tau\right)\int_{0}^{\tau}\sin\left(\frac{\kappa}{\sqrt{3}}\tilde{\tau}\right)\Omega_{\text{GW}}(\tilde{\tau},\kappa)d\tilde{\tau}\right].$
(D.9)

Since $\Omega_{\text{GW}}(\tau,\kappa)$ is $0$ for times before the phase
transition the non-homogeneus part is $0$ when $\tau<\tau_{*}$. The lower
limit of the integral becomes $\tau_{*}$ if $\tau>\tau_{*}$. If we assume that
the source decays quickly after $\tau_{f}$ the upper limit of the integral
becomes $\tau_{f}$ if $\tau>t_{f}$. Therefore, the expression evaluated at
$\tau>\tau_{f}$, i.e. right after the end of the phase transition, becomes

$\displaystyle\tilde{\delta}(\kappa,\tau>\tau_{f})=\,\tilde{\delta}_{h}(\kappa,\tau)+$
$\displaystyle\frac{8\sqrt{3}}{\kappa}\left[\sin\left(\frac{\kappa}{\sqrt{3}}\tau\right)\int_{\tau_{*}}^{\tau_{f}}\cos\left(\frac{\kappa}{\sqrt{3}}\tilde{\tau}\right)\Omega_{\text{GW}}(\tilde{\tau},\kappa)d\tilde{\tau}-\cos\left(\frac{\kappa}{\sqrt{3}}\tau\right)\int_{\tau_{*}}^{\tau_{f}}\sin\left(\frac{\kappa}{\sqrt{3}}\tilde{\tau}\right)\Omega_{\text{GW}}(\tilde{\tau},\kappa)d\tilde{\tau}\right].$
(D.10)

We can see that before $\tau_{*}$, $\tilde{\delta}$ behaves as an harmonic
oscillating function with constants $A_{k}$ and $B_{k}$ given by an initial
value. During the phase transition, the time dependence of $\tilde{\delta}$
will be a complicated function. However, for times $\tau>\tau_{f}$ the
integrals become constants in time (they still depend on $\kappa$) and
therefore $\tilde{\delta}$ is again an harmonic oscillation with modified
amplitudes for each $\kappa$.

We can now solve these integrals in the limits for low and high $\kappa$ by
substituting $\Omega_{\text{GW}}$ in Eq. (D) by Eqs. (D.6) and (D.8). The
solution for $\tilde{\delta}^{(2)\,\text{low}}$ becomes

$\displaystyle\tilde{\delta}^{\text{low}}(\kappa,\tau>\tau_{f})$
$\displaystyle=\tilde{\delta}_{h}(\kappa,\tau)+\tilde{\delta}_{GW}(\kappa,\tau)$
$\displaystyle=\left(A_{\kappa}+A_{GW\kappa}^{\text{low}}\right)\cos\left(\frac{\kappa\tau}{\sqrt{3}}\right)+\left(B_{\kappa}+B_{GW\kappa}^{\text{low}}\right)\sin\left(\frac{\kappa\tau}{\sqrt{3}}\right),$
(D.11)

where $A_{GW\kappa}^{\text{low}}$ and $B_{GW\kappa}^{\text{low}}$ are the
’modified’ amplitudes due to the effect of the FPT in the low k limit. They
are given by

$\displaystyle
A_{GW\kappa}^{\text{low}}=-\frac{24C_{\text{GW}}^{\text{low}}\left(\kappa^{3}-\kappa_{\text{min}}^{3}\right)}{\beta^{6}\kappa^{8}}$
$\displaystyle\left[2\sqrt{3}\beta
H_{*}\kappa\left(2160\beta^{4}\sin\left(\frac{\kappa\tau_{*}}{\sqrt{3}}\right)+\left(1080\beta^{4}+36\beta^{2}H_{*}^{2}\kappa^{2}+H_{*}^{4}\kappa^{4}\right)\sin\left(\frac{\kappa\tau_{f}}{\sqrt{3}}\right)\right)+\right.$
$\displaystyle\left.\left(19440\beta^{6}+216\beta^{4}H_{*}^{2}\kappa^{2}-6\beta^{2}H_{*}^{4}\kappa^{4}-H_{*}^{6}\kappa^{6}\right)\cos\left(\frac{\kappa\tau_{f}}{\sqrt{3}}\right)+432\beta^{4}\left(2H_{*}^{2}\kappa^{2}-45\beta^{2}\right)\cos\left(\frac{\kappa\tau_{*}}{\sqrt{3}}\right)\right],$
(D.12) $\displaystyle
B_{GW\kappa^{\text{low}}}=\frac{24C_{\text{GW}}^{\text{low}}\left(\kappa^{3}-\kappa_{\text{min}}^{3}\right)}{\beta^{6}\kappa^{8}}$
$\displaystyle\left[2\sqrt{3}\beta
H_{*}\kappa\left(1080\beta^{4}+36\beta^{2}H_{*}^{2}\kappa^{2}+H_{*}^{4}\kappa^{4}\right)\cos\left(\frac{\kappa\tau_{f}}{\sqrt{3}}\right)+4320\sqrt{3}\beta^{5}H_{*}\kappa\cos\left(\frac{\kappa\tau_{*}}{\sqrt{3}}\right)+\right.$
$\displaystyle\left.\left(-19440\beta^{6}-216\beta^{4}H_{*}^{2}\kappa^{2}+6\beta^{2}H_{*}^{4}\kappa^{4}+H_{*}^{6}\kappa^{6}\right)\sin\left(\frac{\kappa\tau_{f}}{\sqrt{3}}\right)+432\beta^{4}\left(45\beta^{2}-2H_{*}^{2}\kappa^{2}\right)\sin\left(\frac{\kappa\tau_{*}}{\sqrt{3}}\right)\right].$
(D.13)

Analogously, for high $\kappa$ we have

$\displaystyle\tilde{\delta}^{\text{high}}(\kappa,\tau>\tau_{f})$
$\displaystyle=\tilde{\delta}_{h}(\kappa,\tau)+\tilde{\delta}_{GW}^{\text{high}}(\kappa,\tau)$
$\displaystyle=\left(A_{\kappa}+A_{GW\kappa}^{\text{high}}\right)\cos\left(\frac{\kappa\tau}{\sqrt{3}}\right)+\left(B_{\kappa}+B_{GW\kappa}^{\text{high}}\right)\sin\left(\frac{\kappa\tau}{\sqrt{3}}\right),$
(D.14)

with the constants $A_{GW\kappa}^{\text{high}}$ and
$B_{GW\kappa}^{\text{high}}$ given by

$\displaystyle A_{GW}^{\text{high}}=$ $\displaystyle
A_{GW\kappa_{\text{cross}}}^{\text{low}}-\frac{4C_{\text{GW}}^{\text{high}}(\kappa-\kappa_{\text{cross}})}{3\sqrt{3}H_{*}\kappa^{2}\kappa_{\text{cross}}}\left[6\beta\sin\left(\frac{\kappa\tau_{f}}{\sqrt{3}}\right)-3\beta\sin\left(\frac{\kappa\tau_{*}}{\sqrt{3}}\right)+\right.$
$\displaystyle\left.2\sqrt{3}\kappa\left(H_{*}\left(\text{Ci}\left(\frac{H_{*}\kappa}{\sqrt{3}\beta}\right)-\text{Ci}\left(\frac{2H_{*}\kappa}{\sqrt{3}\beta}\right)\right)\cos\left(\frac{\kappa(2\tau_{f}-\tau_{*})}{\sqrt{3}}\right)+\right.\right.$
$\displaystyle\left.\left.H_{*}\left(\text{Si}\left(\frac{H_{*}\kappa}{\sqrt{3}\beta}\right)-\text{Si}\left(\frac{2H_{*}\kappa}{\sqrt{3}\beta}\right)\right)\sin\left(\frac{\kappa(2\tau_{f}-\tau_{*})}{\sqrt{3}}\right)\right)\right],$
(D.15) $\displaystyle B_{GW}^{\text{high}}=$ $\displaystyle
B_{GW\kappa_{\text{cross}}}^{\text{low}}-\frac{4C_{\text{GW}}^{\text{high}}(\kappa-\kappa_{\text{cross}})}{3\sqrt{3}H_{*}\kappa^{2}\kappa_{\text{cross}}}\left[6\beta\cos\left(\frac{\kappa\tau_{f}}{\sqrt{3}}\right)-3\beta\cos\left(\frac{\kappa\tau_{*}}{\sqrt{3}}\right)+\right.$
$\displaystyle\left.2\sqrt{3}\kappa\left(H_{*}\left(\text{Ci}\left(\frac{2H_{*}\kappa}{\sqrt{3}\beta}\right)-\text{Ci}\left(\frac{H_{*}\kappa}{\sqrt{3}\beta}\right)\right)\sin\left(\frac{\kappa(2\tau_{f}-\tau_{*})}{\sqrt{3}}\right)+\right.\right.$
$\displaystyle\left.\left.H_{*}\left(\text{Si}\left(\frac{H_{*}\kappa}{\sqrt{3}\beta}\right)-\text{Si}\left(\frac{2H_{*}\kappa}{\sqrt{3}\beta}\right)\right)\cos\left(\frac{\kappa(2\tau_{f}-\tau_{*})}{\sqrt{3}}\right)\right)\right].$
(D.16)

Here the functions $\text{Ci}(x)$ and $\text{Si}(x)$ are the CosIntegral and
SinIntegral functions, respectively, defined by
$\text{Ci}(x)=\int_{0}^{x}\frac{\cos t}{t}dt$ and
$\text{Si}(x)=\int_{0}^{x}\frac{\sin t}{t}dt$.

In order to compare these results with the numerical solution, we evaluate
Eqs. (D.11) and (D.14) in $\tau=\tau_{f}$. As a a benchmark scenario, we
choose extreme values for the phase transition parameters: $\beta=H_{*}$,
$\alpha\rightarrow\infty$, $v=c$, $k_{\text{eff}}=1$, $\epsilon=0.01$, which
yields $\kappa_{\text{cross}}=3$. In Fig. 15 we show the analytic results in
the two $\kappa$-regimes compared to the numerical solution for zero initial
conditions.

Figure 15: Comparison between the analytic solution for the density
perturbations at second order $\tilde{\delta}^{(2)}$ evaluated at
$\tau=\tau_{f}$ in the low- and high $\kappa$ limit with the full numerical
solution.

## References

  * [1] B.P. Abbott et al.  Binary Black Hole Mergers in the first Advanced LIGO Observing Run.  Phys. Rev. X, 6(4):041015, 2016.  [Erratum: Phys.Rev.X 8, 039903 (2018)].
  * [2] Pau Amaro-Seoane et al.  Laser interferometer space antenna, 2017.
  * [3] K. Kajantie, M. Laine, K. Rummukainen, and Mikhail E. Shaposhnikov.  The Electroweak phase transition: A Nonperturbative analysis.  Nucl. Phys. B, 466:189–258, 1996.
  * [4] F. Csikor, Z. Fodor, and J. Heitger.  Endpoint of the hot electroweak phase transition.  Phys. Rev. Lett., 82:21–24, 1999.
  * [5] Stefano Profumo, Michael J. Ramsey-Musolf, and Gabe Shaughnessy.  Singlet Higgs phenomenology and the electroweak phase transition.  JHEP, 08:010, 2007.
  * [6] Joerg Jaeckel, Valentin V. Khoze, and Michael Spannowsky.  Hearing the signal of dark sectors with gravitational wave detectors.  Phys. Rev. D, 94(10):103519, 2016.
  * [7] Pedro Schwaller.  Gravitational Waves from a Dark Phase Transition.  Phys. Rev. Lett., 115(18):181101, 2015.
  * [8] Mitsuru Kakizaki, Shinya Kanemura, and Toshinori Matsui.  Gravitational waves as a probe of extended scalar sectors with the first order electroweak phase transition.  Phys. Rev. D, 92(11):115007, 2015.
  * [9] Mikael Chala, Germano Nardini, and Ivan Sobolev.  Unified explanation for dark matter and electroweak baryogenesis with direct detection and gravitational wave signatures.  Phys. Rev. D, 94(5):055006, 2016.
  * [10] Andrea Addazi, Antonino Marcianò, António P. Morais, Roman Pasechnik, Rahul Srivastava, and José W.F. Valle.  Gravitational footprints of massive neutrinos and lepton number breaking.  Phys. Lett. B, 807:135577, 2020.
  * [11] Benoit Laurent, James M. Cline, Avi Friedlander, Dong-Ming He, Kimmo Kainulainen, and David Tucker-Smith.  Baryogenesis and gravity waves from a UV-completed electroweak phase transition.  2 2021.