# Induced gravitational waves: the effect of first order tensor perturbations

Raphaël Picard  and Karim A. Malik

###### Abstract

Scalar induced gravitational waves contribute to the cosmological
gravitational wave background. They can be related to the primordial density
power spectrum produced towards the end of inflation and therefore are a
convenient new tool to constrain models of inflation. These waves are sourced
by terms quadratic in perturbations and hence appear at second order in
cosmological perturbation theory. While the focus of research so far was on
purely scalar source terms we also study the effect of including first order
tensor perturbations as an additional source. This gives rise to two
additional source terms: a term quadratic in the tensor perturbations and a
cross term involving mixed scalar and tensor perturbations. We present full
analytical expressions for the spectral density of these new source terms and
discuss their general behaviour. To illustrate the generation mechanism we
study two toy models containing a peak on small scales. For these models we
show that the scalar-tensor contribution becomes non-negligible compared to
the scalar-scalar contribution on smaller scales. We also consider
implications for future gravitational wave surveys.

##  1 Introduction

The first direct detection of gravitational waves by the LIGO/Virgo
collaboration [1] resulting from the merger of two black holes and the
subsequent detections, reignited cosmologists interests in using gravitational
waves as a novel probe into the early universe. Recently, thanks to cosmic
timing arrays, the NANOGrav collaboration [2], the EPTA/InPTA [3, 4, 5], the
CPTA [6] and the PPTA [7, 8] published data that provides evidence for a low
frequency stochastic gravitational wave background (SGWB). Although the source
of the SGWB is as yet undetermined, this evidently is a milestone for
researchers who have predicted a SGWB of cosmological origin. The NANOGrav
collaboration also fitted cosmological sources [9] such as domain walls,
cosmic strings and cosmological phase transitions to their data release, but
stress that their results should not (yet) be taken as evidence for new
physics.

On the theory side, inflation allows for two different and independent ways of
production of gravitational waves, leading to different observational
signatures. At linear order, gravitational waves produced during inflation
have been studied for some time now and have become a wells established tool
to discriminate between models of inflation [10]. They can leave a distinct
imprint in the Cosmic Microwave Background (CMB) in form of B-mode
polarisation. For example, the Planck collaboration has extracted an upper
limit on the tensor-to-scalar ratio $r<0.065$, at $k=0.002\text{
Mpc}^{-1}$[12]. Going beyond linear order in perturbation theory allows us to
study scalar induced gravitational waves (SIGWs), which are gravitational
waves that are sourced by terms quadratic in linear order scalar quantities.
SIGWs were first studied in the $1960^{\prime}s$ by K. Tomita [11] and later
on by Refs. [13] and [14], amongst others. Unlike first order gravitational
waves, they have a non-negligible source term which is related to the
primordial scalar power spectrum. These fluctuations couple to produce
gravitational waves. The Planck collaboration determined the scalar power
spectrum and gives the value of the spectral index $n_{s}=0.9626\pm 0.0057$ at
$k=0.05\text{ Mpc}^{-1}$[12].

In order for SIGWs to be detectable there needs to be some enhancement of the
scalar power spectrum (see e.g. Ref. [15] for a review). Both the parameters
$r$ and $n_{s}$ are constrained at CMB scales, and a priori there is no reason
that these values hold for scales that exit the horizon in the last stages of
inflation. There are multiple ways to enhance the primordial scalar power
spectrum: for example in single field inflation via a period of ultra-slow-
roll [18, 20, 21, 19] or in multi-field models of inflation [22, 23, 24]. The
basic idea, is that the power spectrum of SIGWs ($\mathcal{P}_{SS}$) is
related to the linear power spectrum, $\mathcal{P}_{\zeta}$, as:
$\mathcal{P}_{ss}\propto\mathcal{P}_{\zeta}^{2}$, so that if there is an
enhancement/peak in $\mathcal{P}_{\zeta}$ there will also be one in
$\mathcal{P}_{ss}$. Depending on the scale at which this enhancement happens,
the resulting background of gravitational waves can be in principle observable
by future detectors such as the Laser Interferometer Space Antenna (LISA)
[28], DECi-hertz Interferometer Gravitational wave Observatory (DECIGO) [29]
and the Einstein Telescope (ET) [30] to name a few.

Most recent research has been centered around how scalar fluctuations couple
to produce SIGWs, however, there is also a possibility that a scalar
fluctuation couples with a tensor mode or that two tensor modes couple to
induce more general: second order gravitational waves (SOGWs). Moreover, if
there is also a peak (on smaller scales) of primordial gravitational waves,
then we could expect to see its imprint also at second order. Models of
inflation that produce such a peak typically have a spectator axion coupled to
a $SU(2)$ gauge field [25, 26, 27]. This was looked at in Ref. [17] for
monochromatic power spectra and in Ref. [33] analytical expressions for the
kernels are derived. Recently, Ref. [34] studied scalar-tensor induced
gravitational waves in a radiation dominated universe with a possibility of
having an initial parity violation and Ref. [35] looked at the case where the
primordial power spectrum follows a sudden-broken power-law.

In this paper we look at the production of SOGWs in a homogeneous and
isotropic universe sourced by first order scalar and tensor modes. We will
derive the power spectrum of SOGWs $\mathcal{P}_{h^{(2)}}(\eta,k)$ defined as

$\langle
h^{(2)}_{\lambda}(\eta,\mathbf{k})h^{(2)}_{\lambda^{\prime}}(\eta,\mathbf{k^{\prime}})\rangle=\delta^{(3)}(\mathbf{k}+\mathbf{k^{\prime}})\delta^{\lambda\lambda^{\prime}}\frac{2\pi^{2}}{k^{3}}\mathcal{P}_{h^{(2)}}(\eta,k)\text{
,}$ (1.1)

where $\lambda,\lambda^{\prime}=+,\times$ are the polarization states of the
GWs. We will use the fact that the spectrum can be decomposed into the sum of
the independent contributions of the scalar-scalar, tensor-tensor
($\mathcal{P}_{tt}$) and scalar-tensor induced waves ($\mathcal{P}_{st}$)
power spectra, such that:

$\mathcal{P}_{h^{(2)}}(\eta,k)=\mathcal{P}_{ss}(\eta,k)+\mathcal{P}_{tt}(\eta,k)+\mathcal{P}_{st}(\eta,k)\text{
.}$ (1.2)

Assuming there is a peaked power spectra for both first the order power
spectrum of scalars and tensors, we will parameterise them by four parameters:
$\mathcal{A}_{\zeta}$ (the amplitude of the scalar power spectrum),
$k_{\zeta}$ (the scale at which the scalar power spectrum peaks),
$\mathcal{A}_{h}$ (the amplitude of the tensor power spectrum) and $k_{h}$
(the scale at which the tensor power spectrum peaks). Given these parameters
we will show that there is non-negligible contribution from considering
tensors as a source term. Our results include the density of SOGWs and
expressions for their corresponding kernels.

We begin in Section 2 by investigating what terms source gravitational waves
and solve the equation of motion for these induced tensor modes. Section 3
shows how to obtain the power spectrum of SOGWs, followed by Section 4 where
we derive expressions for the spectral density and plot the latter for peaked
power spectras. We use the ‘mostly plus’ metric, see e.g. Ref. [39], and work
in natural units, i.e. $c=1$, and define $M_{Pl}^{2}=(8\pi G)^{-1}$.

##  2 The second order gravitational wave equation

###  2.1 Extracting gravitational wave source terms

We use the conformal Newtonian gauge in a flat FLRW background, see e.g. Ref.
[39] for a review of cosmological perturbation theory. We do not include
vectors at any order since they are diluted away during inflation. The line
element is given by [36]

$ds^{2}=a^{2}(\eta)\left[-(1+2\Phi+\Phi^{(2)})d\eta^{2}+\left((1-2\Psi-\Psi^{(2)})\delta_{ij}+2\overline{h}_{ij}+\overline{h}^{(2)}_{ij}\right)dx^{i}dx^{j}\right]\text{
,}$ (2.1)

where $\eta$ is conformal time, $a(\eta)$ the conformal scale factor and
$\delta_{ij}$ is the Kronecker delta. Metric perturbations
$\Phi^{(n)},\Psi^{(n)}$ are respectively the lapse and curvature perturbations
and $\overline{h}_{ij}^{(n)}$ are tensorial perturbations, which are
transverse-traceless (TT):
$\partial^{i}\overline{h}^{(n)}_{ij}=\delta^{ij}\overline{h}^{(n)}_{ij}=0$.
Additionally, to reduce clutter, we have taken the convention to only show the
order of perturbation for second order variables and get rid of the overline
which indicates when a tensorial quantity is TT.

We are interested in the second order equation of motion for tensors,
therefore only TT tensorial quantities will contribute. We can extract TT
quantities by applying the TT projection operator $\Lambda^{ij}_{ab}$ to the
Einstein equations at second order. This operator can be defined in the
following way, see Ref. [16],

$\Lambda^{ij}_{ab}=\left(\delta^{i}_{a}-\frac{\partial^{i}\partial_{a}}{\nabla^{2}}\right)\left(\delta^{j}_{b}-\frac{\partial^{j}\partial_{b}}{\nabla^{2}}\right)-\frac{1}{2}\left(\delta_{ab}-\frac{\partial_{a}\partial_{b}}{\nabla^{2}}\right)\left(\delta^{ij}-\frac{\partial^{i}\partial^{j}}{\nabla^{2}}\right)\text{
,}$ (2.2)

where the inverse Laplacian defined by its action on the Laplacian of a tensor
field $X$ as: $\nabla^{-2}(\nabla^{2}X)=X$. The second order Einstein
equations can be found in App. A.3. Hence, in real space we find

$\Lambda^{ij}_{ab}G_{ij}^{(2)}=\frac{1}{M_{Pl}^{2}}\Lambda^{ij}_{ab}T_{ij}^{(2)}\text{
.}$ (2.3)

We study SOGWs sourced by modes produced during inflation and re-entering
during the radiation dominated (RD) epoch, and describe the energy content of
the universe as an adiabatic perfect fluid with energy momentum-tensor

$T_{\mu\nu}=(\rho+P)u_{\mu}u_{\nu}+Pg_{\mu\nu}\text{ .}$ (2.4)

where $\rho$ and $P$ are, respectively, the energy density and pressure, and
$u^{\mu}$ the 4-velocity. Throughout this paper we neglect anisotropic stress,
as it has been shown to have a negligible effect on the spectrum of SIGWs
[14]. Then Eq. (2.3) becomes:

$h_{ab}^{\prime\prime(2)}+2\mathcal{H}h_{ab}^{\prime(2)}-\nabla^{2}h_{ab}^{(2)}=S_{ab}\text{
,}$ (2.5)

where $S_{ab}$ is the source term. We separate the source term into quadratic
pure first order scalar contributions (ss), quadratic pure first order tensor
contributions (tt) and a scalar-tensor cross term (st), such that:

$S_{ab}=S_{ab}^{ss}+S_{ab}^{tt}+S_{ab}^{st}\text{ .}$ (2.6)

These terms are given by:

$\displaystyle S^{ss}_{ab}$
$\displaystyle=\frac{8}{3(1+w)}\left[(\partial_{a}\Psi+\frac{\partial_{a}\Psi^{\prime}}{\mathcal{H}})(\partial_{b}\Psi+\frac{\partial_{b}\Psi^{\prime}}{\mathcal{H}})\right]+4\partial_{a}\Psi\partial_{b}\Psi\text{
,}$ (2.7a)
$\displaystyle\begin{split}S^{tt}_{ab}&=-4h^{cd}\partial_{c}\partial_{d}h_{ab}+4\partial_{d}h_{bc}\partial^{c}h_{a}^{d}-4\partial_{d}h_{bc}\partial^{d}h^{c}_{a}+8h^{dc}\partial_{a}\partial_{c}h_{bd}\\\
&+4h^{c\prime}_{a}h_{bc}^{\prime}+2\partial_{a}h^{cd}\partial_{b}h_{cd}\text{
,}\end{split}$ (2.7b)
$\displaystyle\begin{split}S^{st}_{ab}&=8\Psi\nabla^{2}h_{ab}+8\partial_{c}h_{ab}\partial^{c}\Psi+4h_{ab}(\mathcal{H}(1+3c_{s}^{2})\Psi^{\prime}+(1-c_{s}^{2})\nabla^{2}\Psi)\text{
.}\end{split}$ (2.7c)

Note that these expressions have been simplified using first order equations
of motion, which can be found in App. A.2. Now that we have extracted the
source terms we can move on to solving Eq. (2.5), which we will do in the next
section.

###  2.2 Solving the GW equation

It is convenient to solve this equation in Fourier space, see App. B for
conventions. For each polarization $\lambda=+,\times$, Eq. (2.5) turns into

$h_{\lambda}^{(2)\prime\prime}(\eta,\mathbf{k})+2\mathcal{H}h_{\lambda}^{(2)\prime}(\eta,\mathbf{k})+k^{2}h_{\lambda}^{(2)}(\eta,\mathbf{k})=4\mathcal{S}_{\lambda}(\eta,\mathbf{k})\text{
,}$ (2.8)

with
$S_{\lambda}(\eta,\mathbf{k})=S_{\lambda}^{ss}(\eta,\mathbf{k})+S_{\lambda}^{tt}(\eta,\mathbf{k})+S_{\lambda}^{st}(\eta,\mathbf{k})$.
To solve the above, we use Green’s method:

$h^{(2)}_{\lambda}(\eta,\mathbf{k})=\frac{4}{a(\eta)}\int^{\eta}d\overline{\eta}G_{\mathbf{k}}(\eta,\overline{\eta})a(\overline{\eta})S_{\lambda}(\overline{\eta},\mathbf{k})\text{
,}$ (2.9)

where $G_{\mathbf{k}}(\eta,\overline{\eta})$ solves:

$G_{\mathbf{k}}^{\prime\prime}(\eta,\overline{\eta})+2\mathcal{H}G_{\mathbf{k}}^{\prime}(\eta,\overline{\eta})+k^{2}G_{\mathbf{k}}(\eta,\overline{\eta})=\delta(\eta-\overline{\eta})\text{
.}$ (2.10)

In our work we look at modes which re-enter the horizon and interact during
RD, hence $\mathcal{H}=1/\eta$, which leads to:

$G_{\mathbf{k}}(\eta,\overline{\eta})=\frac{1}{k}\sin{(k\eta-k\overline{\eta})}\text{
.}$ (2.11)

The solution of Eq. (2.9) depends on the evolution of the first order scalar
and tensor modes. When these modes are super-Hubble their amplitudes are
constant (i.e. when $k\eta<1$) and set during inflation. This implies we can
split their respective equations of motion as a combination of its primordial
value and a transfer function. For the first order scalar modes we define

$\Psi(\eta,\mathbf{k})=\Psi_{\mathbf{k}}\frac{3}{(c_{s}x)^{2}}\left(\frac{\sin{c_{s}x}}{c_{s}x}-\cos{c_{s}x}\right)\equiv\frac{2}{3}\mathbf{\zeta}_{\mathbf{q}}T_{\Psi}(c_{s}x)\text{
,}$ (2.12)

where $\zeta$ is the comoving curvature perturbation conserved when the mode
is super-Hubble. For the first order tensor modes:

$h^{\lambda}(\eta,\mathbf{q})=h^{\lambda}_{\mathbf{k}}\frac{\sin{x}}{x}\equiv
h^{\lambda}_{\mathbf{k}}T_{h}(x)\text{ ,}$ (2.13)

where we have defined $x=k\eta$ in order to solve Eq. (A.6) and Eq. (A.8) in
Fourier space. Both transfer functions are symmetric
$\mathbf{k}\rightarrow-\mathbf{k}$.

Finally, splitting the contributions to Eq. (2.9) into pure scalar, tensor and
mixed we arrive at using Eqs. (2.7), (2.12) and (2.13):

$\displaystyle h_{\lambda,ss}^{(2)}$
$\displaystyle=\frac{16}{9}\int\frac{d^{3}\mathbf{p}}{(2\pi)^{\frac{3}{2}}}q^{\lambda,ss}(\mathbf{k},\mathbf{p})\zeta_{\mathbf{p}}\zeta_{\mathbf{k}-\mathbf{p}}\int^{\eta}_{0}d\overline{\eta}\frac{a(\overline{\eta})}{a(\eta)}G_{\mathbf{k}}(\eta,\overline{\eta})f^{ss}(c_{s}\overline{\eta}p,c_{s}\overline{\eta}|\mathbf{k}-\mathbf{p}|)\text{
,}$ (2.14a)
$\displaystyle\begin{split}h_{\lambda,tt}^{(2)}&=4\int\frac{d^{3}\mathbf{p}}{(2\pi)^{\frac{3}{2}}}q_{1}^{\lambda,tt}(\mathbf{k},\mathbf{p})h_{\mathbf{p}}^{\lambda_{1}}h_{\mathbf{k}-\mathbf{p}}^{\lambda_{2}}\int^{\eta}_{0}d\overline{\eta}\frac{a(\overline{\eta})}{a(\eta)}G_{\mathbf{k}}(\eta,\overline{\eta})f_{1}^{tt}(\overline{\eta}p,\overline{\eta}|\mathbf{k}-\mathbf{p}|)\\\
&+4\int\frac{d^{3}\mathbf{p}}{(2\pi)^{\frac{3}{2}}}q_{2}^{\lambda,tt}(\mathbf{k},\mathbf{p})h_{\mathbf{p}}^{\lambda_{1}}h_{\mathbf{k}-\mathbf{p}}^{\lambda_{2}}\int^{\eta}_{0}d\overline{\eta}\frac{a(\overline{\eta})}{a(\eta)}G_{\mathbf{k}}(\eta,\overline{\eta})f_{2}^{tt}(\overline{\eta}p,\overline{\eta}|\mathbf{k}-\mathbf{p}|)\text{
,}\end{split}$ (2.14b)
$\displaystyle\begin{split}h_{\lambda,st}^{(2)}&=\frac{8}{3}\int\frac{d^{3}\mathbf{p}}{(2\pi)^{\frac{3}{2}}}q^{\lambda,st}(\mathbf{k},\mathbf{p})h_{\mathbf{p}}^{\lambda_{3}}\zeta_{\mathbf{k}-\mathbf{p}}\int^{\eta}_{0}d\overline{\eta}\frac{a(\overline{\eta})}{a(\eta)}G_{\mathbf{k}}(\eta,\overline{\eta})f^{st}(\overline{\eta}p,c_{s}\overline{\eta}|\mathbf{k}-\mathbf{p}|)\text{
.}\end{split}$ (2.14c)

Note that there is an implicit sum over
$\lambda_{1},\lambda_{2},\lambda_{3}=+,\times$ (i.e. for each first order
gravitational wave source terms there is a sum over the two polarizations).
Each contribution (scalar-scalar, tensor-tensor and scalar-tensor) has been
split into three separate parts: First, the functions
$q^{\lambda}(\mathbf{k},\mathbf{p})$ encoding the polarization, shown below:

$\displaystyle q^{\lambda,ss}(k,p)$
$\displaystyle=q^{ab}_{\lambda}(\mathbf{k})p_{a}p_{b}\text{ ,}$ (2.15a)
$\displaystyle\begin{split}q^{\lambda,tt}_{1}(k,p)&=q^{ab}_{\lambda}(\mathbf{k})\bigg{(}q_{ab}^{\lambda_{2}}(\mathbf{p})q^{cd}_{\lambda_{1}}(\mathbf{k}-\mathbf{p})p_{c}p_{d}-q_{bc}^{\lambda_{1}}(\mathbf{k}-\mathbf{p})q^{d}_{a,\lambda_{2}}(\mathbf{p})(k-p)_{d}p^{c}\\\
&+q_{bc}^{\lambda_{1}}(\mathbf{k}-\mathbf{p})q^{c}_{a,\lambda_{2}}(\mathbf{p})(k-p)_{d}p^{d}-2q^{dc}_{\lambda_{1}}(\mathbf{k}-\mathbf{p})q^{\lambda_{2}}_{bd}(\mathbf{p})p_{a}p_{c}\\\
&-\frac{1}{2}q^{cd}_{\lambda_{1}}(\mathbf{k}-\mathbf{p})q^{\lambda_{2}}_{cd}(\mathbf{p})(k-p)_{a}p_{b}\bigg{)}\text{
,}\end{split}$ (2.15b) $\displaystyle q^{\lambda,tt}_{2}(k,p)$
$\displaystyle=q^{ab}_{\lambda}(\mathbf{k})q_{a}^{c,\lambda_{1}}(\mathbf{k}-\mathbf{p})q^{\lambda_{2}}_{bc}(\mathbf{p})\text{
.}$ (2.15c) $\displaystyle q^{\lambda,st}(k,p)$
$\displaystyle=q^{ab}_{\lambda}(\mathbf{k})q^{\lambda_{3}}_{ab}(\mathbf{p})\text{
.}$ (2.15d)

The second part is the primordial values of the perturbations $\mathbf{\zeta}$
and $h^{\lambda}$. The third part of Eq. (2.14) is an integral over the
retarded time $\overline{\eta}$ of the Green’s and transfer functions which
describes how modes behave once they re-enter the horizon:

$\displaystyle\begin{split}f^{ss}(c_{s}\eta
p,c_{s}\eta|\mathbf{k}-\mathbf{p}|)&=\frac{2}{3(1+w)}\bigg{[}\bigg{(}(T_{\Psi}(c_{s}\eta|\mathbf{k}-\mathbf{p}|)+\frac{T^{\prime}_{\Psi}(c_{s}\eta|\mathbf{k}-\mathbf{p}|)}{\mathcal{H}}\bigg{)}\bigg{(}T_{\Psi}(c_{s}\eta
p)\\\ &+\frac{T^{\prime}_{\Psi}(c_{s}\eta
p)}{\mathcal{H}}\bigg{)}\bigg{]}+T_{\Psi}(c_{s}\eta|\mathbf{k}-\mathbf{p}|)T_{\Psi}(c_{s}\eta
p)\text{ ,}\end{split}$ (2.16a) $\displaystyle f^{tt}_{1}(\eta
p,\eta|\mathbf{k}-\mathbf{p}|)$
$\displaystyle=T_{h}(\eta|\mathbf{k}-\mathbf{p}|)T_{h}(\eta p)\text{ ,}$
(2.16b) $\displaystyle f^{tt}_{2}(\eta p,\eta|\mathbf{k}-\mathbf{p}|)$
$\displaystyle=T^{\prime}_{h}(\eta|\mathbf{k}-\mathbf{p}|)T^{\prime}_{h}(\eta
p)\text{ ,}$ (2.16c) $\displaystyle\begin{split}f^{st}(\eta
p,c_{s}\eta|\mathbf{k}-\mathbf{p}|)&=-2p^{2}T_{h}(\eta
p)T_{\Psi}(c_{s}\eta|\mathbf{k}-\mathbf{p}|)-2(k-p)pT_{h}(\eta
p)T_{\Psi}(c_{s}\eta|\mathbf{k}-\mathbf{p}|)\\\
&+\mathcal{H}(1+3c_{s}^{2})T_{h}(\eta
p)T_{\Psi}^{\prime}(c_{s}\eta|\mathbf{k}-\mathbf{p}|)-(k-p)^{2}(1-c_{s}^{2})T_{h}(\eta
p)\\\ &\times T_{\Psi}(c_{s}\eta|\mathbf{k}-\mathbf{p}|)\text{ .}\end{split}$
(2.16d)

Note that the source terms in Fourier space are given by a convolution of the
first order perturbations, hence there is a freedom in relabelling (i.e. a
symmetry between) $\mathbf{p}\leftrightarrow\mathbf{k}-\mathbf{p}$. This
symmetry will be used to simplify calculations when considering the two-point
correlation function of $h_{\lambda}^{(2)}$, in the following section.

##  3 Extracting the power spectrum

To extract the power spectrum we substitute Eq. (2.9) into the left-hand side
of Eq. (1.1). By factorizing out each terms primordial value in Eq. (2.14), it
becomes more apparent what contributes to the final power spectrum:

$\displaystyle h^{\lambda(2)}_{ss}(\eta,\mathbf{k})$
$\displaystyle=\frac{16}{9k^{2}}\int\frac{d^{3}\mathbf{p}}{(2\pi)^{\frac{3}{2}}}q^{\lambda,ss}(k,p)I^{ss}\zeta_{\mathbf{p}}\zeta_{\mathbf{k}-\mathbf{p}}\text{
,}$ (3.1a) $\displaystyle h^{\lambda(2)}_{tt}(\eta,\mathbf{k})$
$\displaystyle=\frac{4}{k^{3}\eta}\int\frac{d^{3}\mathbf{p}}{(2\pi)^{\frac{3}{2}}}\left(q_{1}^{\lambda,tt}(k,p)\mathcal{I}_{tt,1}+q_{2}^{\lambda,tt}(k,p)\mathcal{I}_{tt,2}\
\right)h_{\mathbf{p}}^{\lambda_{1}}h_{\mathbf{k}-\mathbf{p}}^{\lambda_{2}}\text{
,}$ (3.1b) $\displaystyle h^{\lambda(2)}_{st}(\eta,\mathbf{k})$
$\displaystyle=\frac{8}{3k^{3}\eta}\int\frac{d^{3}\mathbf{p}}{(2\pi)^{\frac{3}{2}}}q^{\lambda,st}(k,p)\mathcal{I}_{st}h_{\mathbf{p}}^{\lambda_{3}}\zeta_{\mathbf{k}-\mathbf{p}}\text{
.}$ (3.1c)

The tensor-tensor contribution has two terms due to the nature of the
polarization functions Eq. (2.15b) and Eq. (2.15c), which means we can’t
factor out a common factor. In the above we introduced the kernels,
$I\equiv(x,p,|\mathbf{k}-\mathbf{p}|)$, defined as:

$\displaystyle\mathcal{I}_{ss}(x,p,|\mathbf{k}-\mathbf{p}|)$
$\displaystyle=\int^{x}_{0}d\overline{x}f^{ss}(c_{s}\overline{\eta}p,c_{s}\overline{\eta}|\mathbf{k}-\mathbf{p}|)\frac{\overline{x}}{x}\sin(x-\overline{x})\text{
,}$ (3.2a) $\displaystyle\mathcal{I}_{tt,1}(x,p,|\mathbf{k}-\mathbf{p}|)$
$\displaystyle=\int^{x}_{0}d\overline{x}f_{1}^{tt}(\overline{\eta}p,\overline{\eta}|\mathbf{k}-\mathbf{p}|)\overline{x}\sin(x-\overline{x})\text{
,}$ (3.2b) $\displaystyle\mathcal{I}_{tt,2}(x,p,|\mathbf{k}-\mathbf{p}|)$
$\displaystyle=\int^{x}_{0}d\overline{x}f_{2}^{tt}(\overline{\eta}p,\overline{\eta}|\mathbf{k}-\mathbf{p}|)\overline{x}\sin(x-\overline{x})\text{
,}$ (3.2c) $\displaystyle\mathcal{I}_{st}(x,p,|\mathbf{k}-\mathbf{p}|)$
$\displaystyle=\int^{x}_{0}d\overline{x}f^{st}(\overline{\eta}p,c_{s}\overline{\eta}|\mathbf{k}-\mathbf{p}|)\overline{x}\sin(x-\overline{x})\text{
,}$ (3.2d)

where $\overline{x}=k\overline{\eta}$. We note that the pure scalar Kernel has
been defined in this particular way in order to use the results from Ref.
[38].

Using Wick’s theorem and assuming total gaussianity of the linear order
perturbations, the two point function reduces to the sum of the three
individual contributions111 We do not include contributions from third order
terms, such as $\langle
h^{\lambda}(\eta,\mathbf{k})h^{\lambda^{\prime}(3)}(\eta,\mathbf{k^{\prime}})\rangle$.
We thank G. Domènech for pointing this out. (scalar-scalar, tensor-tensor and
scalar-tensor):

$\begin{split}\langle
h^{\lambda(2)}(\eta,\mathbf{k})h^{\lambda^{\prime}(2)}(\eta,\mathbf{k^{\prime}})\rangle=&\langle
h^{\lambda(2)}_{ss}(\eta,\mathbf{k})h^{\lambda^{\prime}(2)}_{ss}(\eta,\mathbf{k^{\prime}})\rangle+\langle
h^{\lambda(2)}_{tt}(\eta,\mathbf{k})h^{\lambda^{\prime}(2)}_{tt}(\eta,\mathbf{k^{\prime}})\rangle\\\
+&\langle
h^{\lambda(2)}_{st}(\eta,\mathbf{k})h^{\lambda^{\prime}(2)}_{st}(\eta,\mathbf{k^{\prime}})\rangle\text{
.}\end{split}$ (3.3)

Here we have used the fact that the two point function
$\langle\mathbf{\zeta_{\mathbf{p}}}h_{\mathbf{p^{\prime}}}\rangle$ vanishes
for any momenta ${\mathbf{p}}$, ${\mathbf{p^{\prime}}}$. This can be seen as a
consequence of the SVT theorem: everything decouples at linear order.

As stated in the introduction, the total power spectrum
$\mathcal{P}_{h^{(2)}}$ can therefore be expressed as the sum of the three
separate contributions $\mathcal{P}_{ss}$, $\mathcal{P}_{tt}$ and
$\mathcal{P}_{st}$, as

$\displaystyle\langle
h^{\lambda(2)}(\eta,\mathbf{k})h^{\lambda^{\prime}(2)}(\eta,\mathbf{k^{\prime}})\rangle$
$\displaystyle=$
$\displaystyle\delta^{(3)}(\mathbf{k}+\mathbf{k^{\prime}})\delta^{\lambda\lambda^{\prime}}\frac{2\pi^{2}}{k^{3}}\\{\mathcal{P}_{ss}+\mathcal{P}_{tt}+\mathcal{P}_{st}\\}\text{
.}$ (3.4)

It is useful to introduce the variables:

$v=\frac{p}{k},\quad u=\frac{|\mathbf{k}-\mathbf{p}|}{k}\text{ ,}$ (3.5)

so that the integral over the wave-number
$d^{3}\mathbf{p}=dpd\theta_{p}d\phi_{p}p^{2}\sin{\theta_{p}}$ can be written
as:

$\int
d^{3}\mathbf{p}=k^{3}\int_{0}^{\infty}dv\int_{|v-1|}^{1+v}du(uv)\int_{0}^{2\pi}d\phi_{p}\text{
.}$ (3.6)

The reader can refer to App. D.1 to see how the functions presented in Eq.
(2.16) are expressed in the new coordinate system.

In the following, we calculate each contribution to the power spectrum of
SOGWs and assume no parity violation in the primordial tensor power spectrum.

###  3.1 Scalar-scalar power spectrum

In this section we present the results for the pure scalar contribution to the
final power spectrum. This has also been done in Refs.[15, 38, 37], which the
reader can refer to for more details. We need to consider: $\langle
h^{\lambda(2)}_{ss}(\eta,\mathbf{k})h^{\lambda^{\prime}(2)}_{ss}(\eta,\mathbf{k^{\prime}})\rangle$,
which is given by:

$\displaystyle\langle
h^{\lambda(2)}_{ss}(\eta,\mathbf{k})h^{\lambda^{\prime}(2)}_{ss}(\eta,\mathbf{k^{\prime}})\rangle$
$\displaystyle=$
$\displaystyle\left(\frac{16}{9}\right)^{2}\frac{1}{k^{2}k^{\prime
2}}\int\frac{d^{3}\mathbf{p}}{(2\pi)^{\frac{3}{2}}}\int\frac{d^{3}\mathbf{p^{\prime}}}{(2\pi)^{\frac{3}{2}}}q^{\lambda,ss}(k,p)q^{\lambda^{\prime},ss}(k^{\prime},p^{\prime})$
(3.7) $\displaystyle\times
I^{ss}(x,v,u)I^{ss}(x^{\prime},u^{\prime},v^{\prime})\langle\zeta_{\mathbf{p}}\zeta_{\mathbf{k-p}}\zeta_{\mathbf{p^{\prime}}}\zeta_{\mathbf{k^{\prime}-p^{\prime}}}\rangle\text{
.}$

The scalar four point function can be broken down into two-point correlations
(Wick’s theorem):

$\displaystyle\langle\zeta_{\mathbf{p}}\zeta_{\mathbf{k-p}}\zeta_{\mathbf{p^{\prime}}}\zeta_{\mathbf{k^{\prime}-p^{\prime}}}\rangle$
$\displaystyle=$
$\displaystyle\langle\zeta_{\mathbf{p}}\zeta_{\mathbf{p^{\prime}}}\rangle\langle\zeta_{\mathbf{k-p}}\zeta_{\mathbf{k^{\prime}-p^{\prime}}}\rangle+\langle\zeta_{\mathbf{p}}\zeta_{\mathbf{k^{\prime}-p^{\prime}}}\rangle\langle\zeta_{\mathbf{k-p}}\zeta_{\mathbf{p^{\prime}}}\rangle$
(3.8) $\displaystyle=$
$\displaystyle\delta(\mathbf{k}+\mathbf{k^{\prime}})\left(\delta(\mathbf{p}+\mathbf{p^{\prime}})+\delta(\mathbf{p}+\mathbf{k^{\prime}}-\mathbf{p^{\prime}})\right)$
$\displaystyle\times$
$\displaystyle\frac{2\pi^{2}}{p^{3}}\frac{2\pi^{2}}{|\mathbf{k}-\mathbf{p}|^{3}}\mathcal{P}_{\zeta}(p)\mathcal{P}_{\zeta}(|\mathbf{k}-\mathbf{p}|)\text{
.}$

Inserting the above into Eq. (3.7) we find:

$\displaystyle\langle
h^{\lambda(2)}_{ss}(\eta,\mathbf{k})h^{\lambda^{\prime}(2)}_{ss}(\eta,\mathbf{k^{\prime}})\rangle$
$\displaystyle=$
$\displaystyle\delta^{(3)}(\mathbf{k}+\mathbf{k^{\prime}})\delta^{\lambda\lambda^{\prime}}\frac{2\pi^{2}}{k^{3}}\int_{0}^{\infty}du\int_{|1-u|}^{u+1}dv\frac{64}{81}\mathcal{P}_{h}(kv)\mathcal{P}_{\zeta}(ku)$
(3.9) $\displaystyle\times$
$\displaystyle\frac{v^{2}}{u^{2}}\left(1-\frac{(1v^{2}-u^{2})^{2}}{4v^{2}}\right)^{2}\mathcal{I}_{ss}^{2}(x,v,u)$

and we can extract the power spectrum:

$\mathcal{P}_{ss}=\frac{64}{81}\int_{0}^{\infty}dv\int_{|v-1|}^{1+v}du\frac{v^{2}}{u^{2}}\bigg{(}1-\frac{(1v^{2}-u^{2})^{2}}{4v^{2}}\bigg{)}^{2}\mathcal{I}_{ss}^{2}(x,v,u)\mathcal{P}_{\zeta}(ku)\mathcal{P}_{\zeta}(kv)$
(3.10)

in agreement with the literature.

###  3.2 Tensor-tensor power spectrum

In this section we are interested in extracting the pure tensor power spectrum
contribution, corresponding to $\langle
h^{\lambda(2)}_{tt}(\eta,\mathbf{k})h^{\lambda^{\prime}(2)}_{tt}(\eta,\mathbf{k^{\prime}})\rangle$
in Eq. (3.3):

$\displaystyle\langle
h^{\lambda(2)}_{tt}(\eta,\mathbf{k})h^{\lambda^{\prime}(2)}_{tt}(\eta,\mathbf{k^{\prime}})\rangle$
$\displaystyle=$
$\displaystyle\left(\frac{4}{\eta}\right)^{2}\frac{1}{k^{3}k^{\prime
3}}\int\frac{d^{3}\mathbf{p}}{(2\pi)^{\frac{3}{2}}}\int\frac{d^{3}\mathbf{p^{\prime}}}{(2\pi)^{\frac{3}{2}}}\bigg{(}q_{1}^{\lambda,tt}(k,p)\mathcal{I}_{tt,1}(x,v,u)$
(3.11)
$\displaystyle+q_{2}^{\lambda,tt}(k,p)\mathcal{I}_{tt,2}(x,v,u)\bigg{)}\bigg{(}q_{1}^{\lambda^{\prime},tt}(k^{\prime},p^{\prime})\mathcal{I}_{tt,1}(x^{\prime},u^{\prime},v^{\prime})$
$\displaystyle+q_{2}^{\lambda^{\prime},tt}(k^{\prime},p^{\prime})\mathcal{I}_{tt,2}(x^{\prime},u^{\prime},v^{\prime})\bigg{)}\langle
h^{\lambda_{1}}_{\mathbf{p}}h^{\lambda_{2}}_{\mathbf{k-p}}h^{\lambda_{1}^{\prime}}_{\mathbf{p^{\prime}}}h^{\lambda_{2}^{\prime}}_{\mathbf{p^{\prime}-k^{\prime}}}\rangle\text{
.}$

The tensor four-point correlation function is given by:

$\displaystyle\langle
h^{\lambda_{1}}_{\mathbf{p}}h^{\lambda_{2}}_{\mathbf{k-p}}h^{\lambda_{1}^{\prime}}_{\mathbf{p^{\prime}}}h^{\lambda_{2}^{\prime}}_{\mathbf{p^{\prime}-k^{\prime}}}\rangle$
$\displaystyle=$ $\displaystyle\langle
h^{\lambda_{1}}_{\mathbf{p}}h^{\lambda_{1}^{\prime}}_{\mathbf{p^{\prime}}}\rangle\langle
h^{\lambda_{2}}_{\mathbf{k-p}}h^{\lambda_{2}^{\prime}}_{\mathbf{p^{\prime}-k^{\prime}}}\rangle+\langle
h^{\lambda_{1}}_{\mathbf{p}}h^{\lambda_{2}^{\prime}}_{\mathbf{p^{\prime}-k^{\prime}}}\rangle\langle
h^{\lambda_{2}}_{\mathbf{k-p}}h^{\lambda_{1}^{\prime}}_{\mathbf{p^{\prime}}}\rangle$
(3.12) $\displaystyle=$
$\displaystyle\delta(\mathbf{k}+\mathbf{k^{\prime}})\left(\delta(\mathbf{p}+\mathbf{p^{\prime}})\delta^{\lambda_{1}\lambda_{1}^{\prime}}\delta^{\lambda_{2}\lambda_{2}^{\prime}}+\delta(\mathbf{p}+\mathbf{k^{\prime}}-\mathbf{p^{\prime}})\delta^{\lambda_{1}\lambda_{2}^{\prime}}\delta^{\lambda_{2}\lambda_{1}^{\prime}}\right)$
$\displaystyle\times$
$\displaystyle\frac{2\pi^{2}}{p^{3}}\frac{2\pi^{2}}{|\mathbf{k}-\mathbf{p}|^{3}}\mathcal{P}_{h}(p)\mathcal{P}_{h}(|\mathbf{k}-\mathbf{p}|)\text{
.}$

Here $\delta^{\lambda_{1}\lambda_{1}^{\prime}}$ is the Kronecker delta which
ensures the correct polarizations are selected. In App. C.1 we show that the
integral over $\mathbf{p}$ is invariant under
$\mathbf{p}\rightarrow\mathbf{k}-\mathbf{p}$. This implies that when carrying
out the integral over $\mathbf{p}^{\prime}$ the contribution from integrating
over $\delta(\mathbf{p}+\mathbf{p^{\prime}})$ or
$\delta(\mathbf{p}+\mathbf{k^{\prime}}-\mathbf{p^{\prime}})$ is identical. By
remembering that the Kronecker deltas multiplying these Dirac delta functions
are ultimately switchable labels, we can carry out the integral over
$\mathbf{p^{\prime}}$ in Eq. (3.11) over one of the Dirac deltas and multiply
the final result by two. Switching to coordinates $v$ and $u$ and integrating
over the azimuthal angle leads to:

$\displaystyle\langle
h^{\lambda(2)}_{tt}(\eta,\mathbf{k})h^{\lambda^{\prime}(2)}_{tt}(\eta,\mathbf{k^{\prime}})\rangle=\delta^{(3)}(\mathbf{k}+\mathbf{k^{\prime}})\delta^{\lambda\lambda^{\prime}}\frac{2\pi^{2}}{k^{3}}\int_{0}^{\infty}dv\int_{|v-1|}^{1+v}du\frac{1}{x^{2}}\mathcal{P}_{h}(ku)\mathcal{P}_{h}(kv)$
$\displaystyle\bigg{[}\left(\frac{4v^{2}}{u^{2}}\left(1-\frac{(1+v^{2}-u^{2})^{2}}{4v^{2}}\right)^{2}\mathcal{I}_{tt,1}^{2}(x,v,u)\right)+\left(\frac{4}{k^{4}u^{2}v^{2}}\left(1-\frac{(1+v^{2}-u^{2})^{2}}{4v^{2}}\right)^{2}\mathcal{I}_{tt,2}^{2}(x,v,u)\right)$
$\displaystyle-\left(\frac{8}{k^{2}u^{2}}\left(1-\frac{(1+v^{2}-u^{2})^{2}}{4v^{2}}\right)^{2}\mathcal{I}_{tt,1}(x,v,u)\mathcal{I}_{tt,2}(x,v,u)\right)\bigg{]}\text{
,}$ (3.13)

where the expressions that appear in front of the kernels result from integral
of the polarization functions, as shown in App. C. Similarly to Ref. [37] we
split our Kernels as:

$\displaystyle\begin{split}\mathcal{I}_{tt,1}^{2}(x,v,u)&=\cos^{2}(x)[\mathcal{I}^{tt}_{c,1}(x,v,u)]^{2}+\sin^{2}(x)[\mathcal{I}^{tt}_{s,1}(x,v,u)]^{2}\\\
&+\sin(2x)\mathcal{I}^{tt}_{c,1}(x,v,u)\mathcal{I}^{tt}_{s,1}(x,v,u)\text{
,}\end{split}$ (3.14a)
$\displaystyle\begin{split}\mathcal{I}_{tt,2}^{2}(x,v,u)&=\cos^{2}(x)[\mathcal{I}^{tt}_{c,2}(x,v,u)]^{2}+\sin^{2}(x)[\mathcal{I}^{tt}_{s,2}(x,v,u)]^{2}\\\
&+\sin(2x)\mathcal{I}^{tt}_{c,2}(x,v,u)\mathcal{I}^{tt}_{s,2}(x,v,u)\text{
,}\end{split}$ (3.14b)
$\displaystyle\begin{split}\mathcal{I}_{tt,1}(x,v,u)\mathcal{I}_{tt,2}(x,v,u)&=\cos^{2}(x)\mathcal{I}^{tt}_{c,1}(x,v,u)\mathcal{I}^{tt}_{c,2}(x,v,u)+\sin^{2}(x)\mathcal{I}^{tt}_{s,1}(x,v,u)\mathcal{I}^{tt}_{s,2}(x,v,u)\\\
&+\cos{x}\sin{x}\mathcal{I}^{tt}_{c,1}(x,v,u)\mathcal{I}^{tt}_{s,2}(x,v,u)+\cos{x}\sin{x}\mathcal{I}^{tt}_{s,1}(x,v,u)\mathcal{I}^{tt}_{c,2}(x,v,u)\end{split}$
(3.14c)

where

$\displaystyle\mathcal{I}_{c,1}^{tt}(x,v,u)=\int_{0}^{x}d\overline{x}\sin(-\overline{x})\overline{x}f^{tt}_{1}(\overline{x},v,u)\text{
,}$ (3.15a)
$\displaystyle\mathcal{I}_{s,1}^{tt}(x,v,u)=\int_{0}^{x}d\overline{x}\cos(\overline{x})\overline{x}f^{tt}_{1}(\overline{x},v,u)\text{
,}$ (3.15b)
$\displaystyle\mathcal{I}_{c,2}^{tt}(x,v,u)=\int_{0}^{x}d\overline{x}\sin(-\overline{x})\overline{x}f^{tt}_{2}(\overline{x},v,u)\text{
,}$ (3.15c)
$\displaystyle\mathcal{I}_{s,2}^{tt}(x,v,u)=\int_{0}^{x}d\overline{x}\cos(\overline{x})\overline{x}f^{tt}_{2}(\overline{x},v,u)\text{
.}$ (3.15d)

By comparing Eq. (3.4) to Eq. (3.2) it is straightforward to see what
contributes to the final power spectrum:

$\displaystyle\mathcal{P}_{tt}=$
$\displaystyle\frac{1}{x^{2}}\int_{0}^{\infty}dv\int_{|v-1|}^{1+v}du\left[1-\frac{(1+v^{2}-u^{2})^{2}}{4v^{2}}\right]^{2}\mathcal{P}_{h}(ku)\mathcal{P}_{h}(kv)$
(3.16)
$\displaystyle\times\bigg{(}\frac{4v^{2}}{u^{2}}\mathcal{I}_{tt,1}^{2}(x,v,u)-\frac{8}{k^{2}u^{2}}\mathcal{I}_{tt,1}(x,v,u)\mathcal{I}_{tt,2}(x,v,u)+\frac{4}{k^{4}u^{2}v^{2}}\mathcal{I}_{tt,2}^{2}(x,v,u)\bigg{)}\text{.}$

###  3.3 Scalar-tensor power spectrum

In this section we show how to extract the mixed-scalar power spectrum
contribution, corresponding to $\langle
h^{\lambda(2)}_{st}(\eta,\mathbf{k})h^{\lambda^{\prime}(2)}_{st}(\eta,\mathbf{k^{\prime}})\rangle$
in Eq. (3.3):

$\displaystyle\langle
h^{\lambda(2)}_{st}(\eta,\mathbf{k})h^{\lambda^{\prime}(2)}_{st}(\eta,\mathbf{k^{\prime}})\rangle$
$\displaystyle=$
$\displaystyle\left(\frac{8}{3\eta}\right)^{2}\frac{1}{k^{3}k^{\prime
3}}\int\frac{d^{3}\mathbf{p}}{(2\pi)^{\frac{3}{2}}}\int\frac{d^{3}\mathbf{p^{\prime}}}{(2\pi)^{\frac{3}{2}}}q^{\lambda,st}(k,p)q^{\lambda^{\prime},st}(k^{\prime},p^{\prime})$
(3.17)
$\displaystyle\times\mathcal{I}_{st}(x,v,u)\mathcal{I}_{st}(x^{\prime},u^{\prime},v^{\prime})\langle
h^{\lambda_{3}}_{\mathbf{p}}\zeta_{\mathbf{k-p}}h^{\lambda_{3}^{\prime}}_{\mathbf{p^{\prime}}}\zeta_{\mathbf{k^{\prime}-p^{\prime}}}\rangle\text{
.}$

The four-point correlation function is given by:

$\displaystyle\langle
h^{\lambda_{3}}_{\mathbf{p}}\zeta_{\mathbf{k-p}}h^{\lambda_{3}^{\prime}}_{\mathbf{p^{\prime}}}\zeta_{\mathbf{k^{\prime}-p^{\prime}}}\rangle$
$\displaystyle=$ $\displaystyle\langle
h^{\lambda_{3}}_{\mathbf{p}}h^{\lambda_{3}^{\prime}}_{\mathbf{p^{\prime}}}\rangle\langle\zeta_{\mathbf{k-p}}\zeta\mathbf{{}_{k^{\prime}-p^{\prime}}}\rangle$
(3.18) $\displaystyle=$
$\displaystyle\delta(\mathbf{k}+\mathbf{k^{\prime}})\left(\delta(\mathbf{p}+\mathbf{p^{\prime}})\delta^{\lambda_{3}\lambda_{3}^{\prime}}\right)\cdot\frac{2\pi^{2}}{p^{3}}\frac{2\pi^{2}}{|\mathbf{k}-\mathbf{p}|^{3}}\mathcal{P}_{h}(p)\mathcal{P}_{\zeta}(|\mathbf{k}-\mathbf{p}|)\text{
.}$

Following the same steps outlined in the previous section:

$\displaystyle\langle
h^{\lambda(2)}_{st}(\eta,\mathbf{k})h^{\lambda^{\prime}(2)}_{st}(\eta,\mathbf{k^{\prime}})\rangle$
$\displaystyle=$
$\displaystyle\delta^{(3)}(\mathbf{k}+\mathbf{k^{\prime}})\delta^{\lambda\lambda^{\prime}}\frac{2\pi^{2}}{k^{3}}\int_{0}^{\infty}du\int_{|1-u|}^{u+1}dv\frac{16}{9x^{2}}\mathcal{P}_{h}(kv)\mathcal{P}_{\zeta}(ku)$
(3.19) $\displaystyle\times$
$\displaystyle\frac{1}{k^{4}u^{2}v^{2}}M_{st}(v,u)\mathcal{I}_{st}^{2}(x,v,u)\text{
.}$

The function $M_{st}(v,u)$ comes from the integration in the azimuthal
direction and is defined in Eq. C.14. Similarly, the kernels are split up as:

$\begin{split}\mathcal{I}_{st}^{2}(x,v,u)&=\cos^{2}(x)[\mathcal{I}^{st}_{c}(x,v,u)]^{2}+\sin^{2}(x)[\mathcal{I}^{st}_{s}(x,v,u)]^{2}\\\
&+\sin(2x)\mathcal{I}^{st}_{c}(x,v,u)\mathcal{I}^{st}_{s}(x,v,u)\text{
,}\end{split}$ (3.20)

where

$\displaystyle\mathcal{I}_{c}^{st}(x,v,u)=\int_{0}^{x}d\overline{x}\sin(-\overline{x})\overline{x}f^{st}(\overline{x},v,u)\text{
,}$ (3.21a)
$\displaystyle\mathcal{I}_{s}^{st}(x,v,u)=\int_{0}^{x}d\overline{x}\cos(\overline{x})\overline{x}f^{st}(\overline{x},v,u)\text{
.}$ (3.21b)

Hence we have:

$\displaystyle\mathcal{P}_{st}=$
$\displaystyle\frac{1}{x^{2}}\frac{16}{9}\int_{0}^{\infty}dv\int_{|v-1|}^{1+v}du\mathcal{P}_{h}(kv)\mathcal{P}_{\zeta}(ku)M_{st}(v,u)\frac{1}{k^{4}u^{2}v^{2}}\mathcal{I}_{st}^{2}(x,v,u)\text{
.}$ (3.22)

##  4 Results: peaked power spectra

In this section we present our results: expressions for the current spectral
density of SOGWs. Furthermore, we show the resulting spectral density for
peaked power spectra sources, comment on the behaviour of the individual
contributions as well as the total contribution. We also make some comments on
detectability.

###  4.1 Current energy density of gravitational waves

The observable cosmologist use is the spectral density [46], related to the
power spectrum by

$\Omega_{GW}(\eta,k)=\frac{1}{6}\left(\frac{k}{\mathcal{H(\eta)}}\right)^{2}\overline{\mathcal{P}_{h^{(2)}}(\eta,k)}\text{
,}$ (4.1)

where the overline here denotes an oscillation average over a few periods.
Since these waves are generated in an early radiation dominated universe, we
have to take into account how their density will change as the universe
evolves. The density of gravitational waves scales like radiation and
therefore it is possible to relate the spectral density at the time of
creation, $\Omega_{GW}(k)$, to the present the spectral density,
$\Omega_{GW}(\eta_{0},k)$, which takes into account how the universe becomes
matter dominated after the matter-radiation equality. For a detailed
explanation the reader can refer to Ref. [15]:

$\Omega_{GW}(\eta_{0},k)h^{2}=1.62\times
10^{-5}\times\Omega_{GW}(k)=1.62\times
10^{-5}\times\frac{1}{6}\frac{k^{2}}{\mathcal{H}(\eta)^{2}}\overline{\mathcal{P}_{h^{(2)}}(\eta,k)}\text{
,}$ (4.2)

where $h$ is the dimensionless reduced Hubble constant to take into account
the uncertainty in its value. Using $\mathcal{H}(\eta)=1/\eta$ we define:

$\displaystyle\Omega_{GW}(\eta_{0},k)h^{2}=1.62\times
10^{-5}\times\left\\{\Omega_{ss}(k)+\Omega_{tt}(k)+\Omega_{st}(k)\right\\}\text{
,}$ (4.3)

where $\Omega_{ss}(k)$, $\Omega_{tt}(k)$ and $\Omega_{st}(k)$ are respectively
the spectral density from the induced scalar-scalar, tensor-tensor and scalar-
tensor at the time of creation. Below we present the corresponding expressions
and their kernels which are the main results of our paper. To see the details
behind our kernel computations see App. D.2. We note that the spectral density
depends on the time-averaged power spectrum. The time dependence in our
expressions is contained in the kernels, via the variable $x=k\eta$, hence the
following results are used to derive the final expressions:

$\overline{\sin x}=\overline{\cos
x}=0,\quad\overline{\sin^{2}x}=\overline{\cos^{2}x}=\frac{1}{2}\text{ .}$

####  4.1.1 Spectral density of the scalar-scalar contribution

Inserting the power spectrum Eq. (3.10) into Eq. (4.1), we arrive to an
expression for the spectral density from the scalar-scalar contribution

$\Omega_{ss}(k)=\frac{32}{243}\int_{0}^{\infty}dv\int_{|v-1|}^{1+v}du\bigg{(}\frac{4v^{2}-(1+v^{2}-u^{2})^{2}}{4uv}\bigg{)}^{2}\overline{\mathcal{I}_{ss}^{2}(v,u)}\mathcal{P}_{\zeta}(ku)\mathcal{P}_{\zeta}(kv)\text{
,}$ (4.4)

with

$\displaystyle\overline{\mathcal{I}_{ss}^{2}(v,u)}$ $\displaystyle=$
$\displaystyle
2^{2}\frac{1}{2}\left(\frac{3(u^{2}+v^{2}-3)}{4u^{3}v^{3}}\right)^{2}\bigg{[}\left(-4uv+(u^{2}+v^{2}-3)\log\bigg{|}\frac{3-(u+v)^{2}}{3-(u-v)^{2}}\bigg{|}\right)^{2}$
(4.5)
$\displaystyle+\pi^{2}(u^{2}+v^{2}-3)^{2}\Theta(v+u-\sqrt{3})\bigg{]}\text{
.}$

The expression for the Kernel has been taking from Ref. [38]. The factor of
$2^{2}$ infront comes from the difference in decomposing the metric in our two
papers.

####  4.1.2 Spectral density of the tensor-tensor contribution

Similarly, the tensor-tensor spectral density is given by:

$\displaystyle\Omega_{tt}(k)$ $\displaystyle=$
$\displaystyle\frac{1}{6}\int_{0}^{\infty}dv\int_{|v-1|}^{1+v}du\left[1-\frac{(1+v^{2}-u^{2})^{2}}{4v^{2}}\right]^{2}\mathcal{P}_{h}(ku)\mathcal{P}_{h}(kv)$
(4.6)
$\displaystyle\times\bigg{(}\frac{4v^{2}}{u^{2}}\overline{\mathcal{I}_{tt,1}^{2}(v,u)}-\frac{8}{k^{2}u^{2}}\overline{\mathcal{I}_{tt,1}(v,u)\mathcal{I}_{tt,2}(v,u)}+\frac{4}{k^{4}u^{2}v^{2}}\overline{\mathcal{I}_{tt,2}^{2}(v,u)}\bigg{)}\text{
,}$

The time-averaged kernels are shown below:

$\displaystyle\overline{\mathcal{I}_{tt,1}^{2}(v,u)}$
$\displaystyle=\frac{1}{32u^{2}v^{2}}\left(\pi^{2}\Theta(u+v-1)+\log\bigg{|}\frac{1-(u+v)^{2}}{1-(u-v)^{2}}\bigg{|}^{2}\right)\text{
,}$ (4.7a) $\displaystyle\overline{\mathcal{I}_{tt,2}^{2}(v,u)}$
$\displaystyle=\frac{k^{4}}{128u^{2}v^{2}}\left(\pi^{2}(u^{2}+v^{2}-1)^{2}\Theta(u+v-1)+\left(\log\bigg{|}\frac{1-(u+v)^{2}}{1-(u-v)^{2}}\bigg{|}-4vu\right)^{2}\right)\text{
,}$ (4.7b)
$\displaystyle\overline{\mathcal{I}_{tt,1}(v,u)\mathcal{I}_{tt,2}(v,u)}$
$\displaystyle=\frac{k^{2}}{64u^{2}v^{2}}\bigg{(}\pi^{2}(u^{2}+v^{2}-1)\Theta(u+v-1)+\log\bigg{|}\frac{1-(u+v)^{2}}{1-(u-v)^{2}}\bigg{|}\bigg{(}-4vu$
$\displaystyle+\log\bigg{|}\frac{1-(u+v)^{2}}{1-(u-v)^{2}}\bigg{|}\bigg{)}\bigg{)}\text{
.}$ (4.7c)

####  4.1.3 Spectral density of the scalar-tensor contribution

Finally, the spectral density of scalar-tensor induced gravitational waves
reads:

$\displaystyle\Omega_{st}(k)=\frac{32}{243}\int_{0}^{\infty}dv\int_{|v-1|}^{1+v}du\mathcal{P}_{h}(kv)\mathcal{P}_{\zeta}(ku)\frac{1}{k^{4}u^{2}v^{2}}M_{st}(v,u)\overline{\mathcal{I}_{st}^{2}(v,u)}\text{
,}$ (4.8)

with:

$\displaystyle\overline{\mathcal{I}_{st}^{2}(v,u)}$ $\displaystyle=$
$\displaystyle\frac{9k^{4}}{2048v^{2}u^{6}c_{s}^{6}}\bigg{(}\pi^{2}(i^{st}_{-}(v,u))^{2}\Theta(v+uc_{s}-1)+\bigg{(}i_{+}^{st}(v,u)\log\bigg{|}\frac{c_{s}u-v-1}{c_{s}u+v+1}\bigg{|}$
(4.9)
$\displaystyle\times(c_{s}u-v-1)(c_{s}u+v+1)+i_{-}^{st}(v,u)\log\bigg{|}\frac{c_{s}u-v+1}{c_{s}u+v-1}\bigg{|}(c_{s}u-v+1)(c_{s}u+v-1)$
$\displaystyle-$ $\displaystyle
20c_{s}uv\bigg{[}-\frac{3}{5}+\frac{v^{2}}{5}(7-c_{s}^{2})+\frac{8uv}{5}+u^{2}\left(c_{s}^{4}-\frac{7c_{s}^{2}}{15}+\frac{4}{5}\right)-\frac{9c_{s}^{2}}{5}\bigg{]}\bigg{)}^{2}\bigg{)}\text{
.}$

The function $i^{st}_{\pm}(v,u)$ is defined as:

$\displaystyle i^{st}_{\pm}(v,u)$ $\displaystyle=$
$\displaystyle-1+4u^{2}-3c_{s}^{2}(u^{2}+(v\pm
1)^{2}-c_{s}^{2}u^{2})+8uv+v(\mp 2+7v)\text{ .}$ (4.10)

###  4.2 Comments on potential divergences and resonances

The behaviour of the integrands (4.4), (4.6) and (4.8) in the IR and UV limits
contains information of any possible divergences that may appear. The IR limit
(long wave-length) corresponds to taking $k\rightarrow 0$, then $v\rightarrow
1/k$ and $u\rightarrow 1/k$, and all three contributions are well behaved in
this limit. In fact, in the IR limit the scalar-scalar induced waves have a
log-dependent slope [31] and it was shown that this is not the case for the
scalar-tensor induced waves due to scalar and tensor fluctuations not having
the same propagation speed [34].

In the UV limit (short wave-length), as $k\rightarrow\infty$, we have
$v\rightarrow 0$ and $u\rightarrow 1$ for the pure quadratic contributions.
The scalar-scalar limit is finite, whilst the tensor-tensor limit goes as
$v^{-2}$. For the cross scalar-tensor term, taking the UV limit breaks the
symmetry we have between $v$ and $u$, there are now two distinct limits:
$v\rightarrow 0$ and $u\rightarrow 1$, corresponding to the large wavelength
for tensors, and $v\rightarrow 1$ and $u\rightarrow 0$ corresponding to the
large wavelength for scalars. In the former case the integrand diverges as
$v^{-4}$ whilst in the latter as $u^{-6}$. Hence we see that if the first
order power spectra are given by a power law, below a certain order, there
could be an unphysical contribution to the final spectral density. This was
first pointed out for the scalar-tensor case in Ref. [34].

Furthermore, the scalar-scalar, tensor-tensor and scalar-tensor integrands
contain logarithmic terms that could lead to possible divergences, provided
$v$ and $u$ take particular values. By considering ‘momentum conservation’ we
can extract a range allowed values for $v$ and $u$. The wave vectors in
Fourier space of the three modes form a triangle, with sides: $\mathbf{k}$,
$\mathbf{p}$ and $\mathbf{k}-\mathbf{p}$, where the magnitude of the vectors
satisfy the triangle inequality, which in terms of the variables $v$ and $u$
is given by: $|u-v|<1<u+v$. Hence, we see that the tensor-tensor spectral
density will not enter a logarithmic resonance since the value of $v$ and $u$
needed do not satisfy the triangle inequality. On the other hand, the scalar-
scalar contribution can enter into a logarithmic resonance when $v+u=\sqrt{3}$
since this corresponds to a possible configuration of the wave-vectors.
Similarly, there is a possible logarithmic divergence for the scalar-tensor
contribution when $v=1\pm c_{s}u$, but this is prevented from the terms
multiplying the logarithmic terms. However, in the limit $v\rightarrow 1\pm
c_{s}u$ the scalar-tensor integrand tends to $u^{-6}$ which causes a
divergences in the large scalar wavelength limit.

For the reasons outlined above we will exclusively input peaked power spectra,
which we do in the following sections, by first looking at the result induced
waves from a monochromatic power spectra and secondly lognormal power spectra.

###  4.3 Monochromatic power spectra

It is informative to study the spectrum of SOGWs that results from having a
sharp peak in the input power spectrum, modelled by a Dirac delta function.
Although ‘not physical’ (i.e. there is no such thing as a infinitely sharp
peak in nature), the resulting spectral density tells us a lot about the
properties and shape of resulting spectra in general.

The monochromatic power spectrum is given by

$\mathcal{P}_{\zeta,h}(k)=\mathcal{A}_{\zeta,h}\delta\left(\log\frac{k}{k_{\zeta,h}}\right)\text{
,}$ (4.11)

with $\mathcal{A}_{\zeta,h}$ the amplitude of the power spectrum and
$k_{\zeta,h}$ the location of the scalar and tensor peak in the primordial
power spectrum, respectively. Hence we get

$\mathcal{P}_{\zeta,h}(kv)=\mathcal{A}_{\zeta,h}\delta\left(\log\frac{kv}{k_{\zeta,h}}\right)=\mathcal{A}_{\zeta,h}\frac{k_{\zeta,h}}{k}\delta(v-\frac{k_{\zeta,h}}{k})\text{
.}$ (4.12)

The same holds for $\mathcal{P}_{\zeta,h}(ku)$, if we swap $v\leftrightarrow
u$. The way we have expressed the power spectrum makes it straightforward to
use the sifting property of the Dirac delta function, we just have to take the
limits of integration of Eqs. (4.4), (4.6) and (4.8) into account

$\displaystyle\begin{cases}|1-v|<u<1+v\\\
0<v<\infty\end{cases}=\begin{cases}|k-k_{h}|<k_{\zeta}<k+k_{h}\\\
0<k_{h}/k<\infty\end{cases}=|k_{\zeta}-k_{h}|<k<k_{\zeta}+k_{h}\text{.}$
(4.13)

Below we study two different scenarios: first, the case where the primordial
power spectrum for scalars and tensors peaks at the same scale and then the
case where the two peaks are at different scales.

####  4.3.1 Peak at the same scale

In this section, the primordial power spectrum of scalar and tensor
perturbations peaks at the same scale so that $k_{\zeta}=k_{h}\equiv k_{p}$.
The inequality (4.13) reduces to: $0<k<2k_{p}$, which can be encoded by a
Heaviside step function. For convenience, we define $\tilde{k}=k/k_{p}$.

The analytic expression for scalar-scalar part Eq. (4.4) is

$\displaystyle\Omega_{ss}(k)$ $\displaystyle=$
$\displaystyle\frac{32}{243}\mathcal{A}_{\zeta}^{2}\tilde{k}^{-2}(\frac{\tilde{k}^{2}-4}{4})^{2}\overline{I^{ss}(v=\tilde{k}^{-1},u=\tilde{k}^{-1})^{2}}\Theta(2k_{p}-k)\text{
.}$ (4.14)

Figure 1: The spectral density $\Omega_{ss}(k)$ of second order gravitational
waves induced by scalar-scalar modes, normalised by $\mathcal{A}_{\zeta}$,
against a range of $k$ around some pivot scale $k_{p}$. There is a logarithmic
resonance that occurs when $\tilde{k}=2/\sqrt{3}$, represented by the black-
dotted line, as expected from looking at the kernel (4.5). After the peak, the
spectral density decreases rapidly before hitting the cut-off at $k=2k_{p}$.

The resulting behaviour of Eq. (4.14) is shown in Fig. 1. There is a sharp
peak at $\tilde{k}=2/\sqrt{3}$, or when $k=2c_{s}k_{p}$, which corresponds to
$u+v=\sqrt{3}$ as expected following our discussion on resonances in the
scalar-scalar SOGWs in Sect. 4.2. Furthermore, due to momentum conservation,
there is the expected sharp cutoff at $k=2k_{p}$.

For the Dirac peaked spectra the tensor-tensor part Eq. (4.6)

$\displaystyle\Omega_{tt}(k)$ $\displaystyle=$
$\displaystyle\frac{2}{3}\mathcal{A}_{h}^{2}(\frac{\tilde{k}^{2}-4}{4})^{2}\bigg{(}\frac{4}{\tilde{k}^{2}}\overline{(\mathcal{I}_{tt,1}(v=\tilde{k}^{-1},u=\tilde{k}^{-1})^{2}}+\frac{4\tilde{k}^{2}}{k^{4}}\overline{(\mathcal{I}_{tt,2}(v=\tilde{k}^{-1},u=\tilde{k}^{-1})^{2}}$
(4.15)
$\displaystyle-\frac{8}{k^{2}}\overline{\mathcal{I}_{tt,1}(v=\tilde{k}^{-1},u=\tilde{k}^{-1})}\times\overline{\mathcal{I}_{tt,2}(v=\tilde{k}^{-1},u=\tilde{k}^{-1})}\bigg{)}\Theta(2k_{p}-k)\text{
.}$

Figure 2: The spectral density $\Omega_{tt}(k)$ of second order gravitational
waves induced by tensor-tensor modes, normalised by $\mathcal{A}_{h}$, against
a range of $k$ around some pivot scale $k_{p}$. There is no logarithmic
resonance since the magnitude of the wave-vectors needed for the resonance,
$u+v=1$, are not allowed due to momentum conservation. The orders of magnitude
spanned by the tensor-tensor contribution peaks just above $10^{-1}$, an
amplitude much smaller than the scalar-scalar induced waves.

Equation (4.15) is plotted in Fig 2. We note that, as expected, for tensor-
tensor induced gravitational waves there is no resonant peak.

The scalar-tensor or cross-term part of Eq. (4.8) is given by

$\displaystyle\Omega_{st}(k)$ $\displaystyle=$
$\displaystyle\frac{8}{27}\mathcal{A}_{\zeta}\mathcal{A}_{h}\frac{\tilde{k}^{2}}{k^{4}}M_{st}(v=\tilde{k}^{-1},u=\tilde{k}^{-1})\overline{\mathcal{I}_{st}(k,v=\tilde{k}^{-1},u=\tilde{k}^{-1})^{2}}\Theta(2k_{p}-k)\text{.}$
(4.16)

Figure 3: The spectral density of second order gravitational waves induced by
scalar-tensor modes, where the primordial power spectra peak at the same scale
$k_{p}$. The values needed to have a logarithmic resonance are allowed by
momentum conservation but due to the structure of (4.9), the resonance does
not occur. The amplitude spanned by the scalar-tensor induced waves is
comparable to the scalar-scalar case.

Equation (4.16) is shown in Fig. 3. As expected, there is no resonant
behavior. Furthermore, we have compared our results with the results in Ref.
[34]222See their Eqs. (3.12) and (4.2) for comparison. and find that our
spectral densities have an overall similar behaviour.

Figure 4: The total spectral density, where we have set
$\mathcal{A}_{h}=0.1\mathcal{A}_{\zeta}$. The tensor-tensor contribution (blue
line) is subdominant compared to the the scalar-scalar (red line) and scalar-
tensor (green line) contributions. This is not only due to the fact that the
amplitude of the primordial power spectrum for tensors is smaller than the one
for scalars, but mainly due to the structure of the tensor-tensor source
terms. Secondly, when considering the total contribution to the induced waves
(black line) it becomes clear that for high values of $k$ the scalar-tensor
contribution becomes the dominant mode of production for induced gravitational
waves in this case.

Finally, we can compare all three contributions to the total spectral density,
shown in Fig. 4. We see that the tensor-tensor contribution is several orders
of magnitude smaller than the other two types of sources. This is partly due
to the fact that the amplitude of the primordial tensor power spectrum is
smaller than the primordial scalar power spectrum but this is not solely the
reason. Expanding the expressions for the spectral density in $k$ for the
scalar-scalar contribution (4.14) and the scalar-tensor contribution (4.16),
we see that the highest power of $k$ is $k^{14}$ whilst for the tensor-tensor
contribution (4.15) it is $k^{10}$. This implies that the nature of the
tensor-tensor source terms inherently makes them a smaller contribution as
compared to the other modes of production. The second aspect of Fig. 4 that is
important is the behaviour of the total spectral density for high values of
$k$. The scalar-scalar contribution drops around four orders of magnitude
between the resonant peak and the cut-off at $k=2k_{p}$, whilst the scalar-
tensor contribution does not drop off and stays roughly constant until the
cut-off. This shows that the mixed contribution is in fact important and must
be taken into account, if there is a peak in the primordial scalar and tensor
power spectrum.

We find that the total induced gravitational wave spectrum is in overall
agreement with the results in Ref. [17] (see the blue curve in their Fig. 2).
However we found that the amplitude of the spectral density drops off more
steeply after the resonant peak on small scales.

####  4.3.2 Peaks at a different scale

In this section we show how the second order scalar-tensor power spectrum
changes when the primordial power spectra for scalars and tensors peak at
different scales. The expression for the spectral density becomes

$\displaystyle\Omega_{st}(k)$ $\displaystyle=$
$\displaystyle\frac{2}{27}\mathcal{A}_{\zeta}\mathcal{A}_{h}\frac{1}{k^{2}k_{\zeta}k_{h}}M_{st}(v=k_{h}/k,u=k_{\zeta}/k)\overline{\mathcal{I}_{st}(k,v=k_{h}/k,u=k_{\zeta}/k)^{2}}$
(4.17)
$\displaystyle\times\Theta(k-|k_{\zeta}-k_{h}|)\Theta(-k+k_{\zeta}+k_{h})\text{
,}$

where the two Heaviside step functions encode the integral limits in Eq.
(4.13). We look at two different informative cases: first when the tensor
power spectrum peaks at a scale higher than the resonance from the scalar-
scalar contribution ($k_{h}\approx c_{s}k_{\zeta}$), and then at a scale
smaller than the resonance ($k_{h}\approx 4c_{s}k_{\zeta}$). The results are
presented in Fig. 5. We see that when $k_{h}\approx c_{s}k_{\zeta}$ the
scalar-tensor contribution is smaller than when the power spectra peak at the
same scale, in fact we chose this value for $k_{h}$ because this value
corresponds to an approximate threshold: when $k_{h}\leq c_{s}k_{\zeta}$ the
signal will always be smaller than the case when $k_{h}=k_{\zeta}$. Secondly,
we chose the peak scale $k_{h}\approx 4c_{s}k_{\zeta}$ to show that we can get
higher amplitude signals than the previous setup. The short ‘burst’ like
signal is due to the allowed $k$ values encoded in the Heaviside step
functions in Eq. (4.17). Indeed, the closer $k_{h}$ is $k_{\zeta}$ the more
the range of $k$ increases. In Fig. 6 we show all three contributions to the
spectral density. The sharp feature in this signal is due to the fact that the
input is a Dirac delta power spectrum.

Figure 5: The spectral density of gravitational waves induced by scalar-tensor
modes for a range of $k$ values. Here we present three different scenarios:
when the primordial power spectra peak at the same scale (solid line), when
the primordial tensor peaks on larger scales than the scalar-scalar resonant
peak (dotted line), and when the primordial tensor peaks on smaller scales
than the scalar-scalar resonant peak (dashed line). The ‘impulse’ like
behaviour of the dashed line can be understood from momentum conservation
encoded in the Heaviside step functions. Figure 6: Full spectral density,
where we have set $\mathcal{A}_{h}=0.1\mathcal{A}_{\zeta}$, with the scalar-
tensor contribution (dashed green line) this time coming from a primordial
power spectrum which peaks at a smaller scale than the scalar-scalar
resonance.

###  4.4 Lognormal Peak

As discussed earlier, a Dirac delta peak is not realistic, and hence to
discuss the plausibility of detecting these induced waves we should move to a
power spectrum which has some finite width to its peak. This can be achieved
with a Lognormal power spectrum,

$\mathcal{P}_{\zeta,h}=\frac{\mathcal{A}_{\zeta,h}}{\sqrt{2\pi}\sigma}\exp\left({-\frac{\log^{2}(k/k_{\zeta,h})}{2\sigma^{2}}}\right)\,,$
(4.18)

with $\sigma$ controlling the width of the peak and normalised such that
$\int^{\infty}_{-\infty}d\log k\mathcal{P}_{\zeta,h}=\mathcal{A}_{\zeta,h}$.
In the following we take $\sigma=0.1$, so that the general behaviour of the
spectral densities discussed in the previous section (Dirac delta input power
spectra), holds in this scenario. For a detailed discussion of scalar-scalar
induced waves from a lognormal peak, see Ref. [41]. The scale at which we
centre the scalar power spectrum ($k_{\zeta}$) is chosen so that the physical
frequency is centred in the LISA band: $k=f\times(2\pi a_{0})$, with
$f_{\zeta}=3.4$ mHz. We take $\mathcal{A}_{\zeta}\approx 2.1\times 10^{-2}$,
which corresponds to an enhancement of $\mathcal{O}(10^{7})$ of the Planck
2018 data [44], motivated by the enhancement needed to the power spectrum in
order to produce primordial black holes [21]. All spectral densities in the
following are found by numerically integrating the expressions in Eqs. (4.4),
(4.6) and (4.8). In Fig. 7 we present the results of our toy model for
different scenarios. We consider the cases where the power spectrum of scalars
and tensors peak at identical and different scales for different values of
$\mathcal{A}_{h}$ relative to $\mathcal{A}_{\zeta}$. A feature that becomes
apparent is that when $\mathcal{A}_{h}=\mathcal{A}_{\zeta}$ and
$k_{h}=k_{\zeta}$ (top-left plot in Fig. 7), the scalar-tensor contribution
becomes dominant. This could be due to the fact that the scalar-scalar
resonance is less sharp and hence has lower amplitude than in the Dirac delta
case. When we decrease the amplitude of the tensor power spectrum by a factor
of ten (bottom-left plot in Fig. 7) the resonant behaviour becomes apparent in
the total signal and the scalar-tensor signal contributes narrowly on smaller
scales. For the case when $k_{h}=4c_{s}k_{\zeta}$ (top and bottom right plots
in Fig. 7) we can still distinguish the resonant peak and also have an extra
contribution coming from the scalar-tensor terms, in the form of a peak.

Figure 7: Plots of the spectral density against frequency. Each graph
corresponds to a different setup in the value of $\mathcal{A}_{\zeta}$,
$\mathcal{A}_{h}$, $k_{\zeta}$ and $k_{h}$. We have also included the noise
power spectral densities curves [32] for LISA, ET [42] and DECIGO [29], as
well as bounds on primordial gravitational waves from BBN [43]. Furthermore,
the dashed dark blue line is the primordial tensor power spectrum. We see that
the scalar-tensor contribution is non-negligible and in some cases dominant.
The sharp signals on larger scales and below the threshold for detection are
artifacts of numerical integration.

##  5 Discussion

In higher order cosmological perturbation theory lower order terms can couple
to source the higher order terms. Hence second order gravitational waves will
be sourced by coupled first order scalar, vector, and tensor terms. Until
recently, research has focused on source terms which are quadratic in first
order scalar fluctuations, so-called SIGWs (see e.g. Ref. [15]). In this
paper, we investigated what happens if we include first order tensor
fluctuations to contribute as a possible source. Two new sources appear: terms
quadratic in tensor fluctuations and terms that couple a scalar and a tensor
fluctuation. We have shown that the tensor-tensor contribution is subdominant
due to its inherit structure whilst the scalar-tensor contribution is non
negligible, particularly on small scales, in agreement with recent studies
[17, 34, 35]. Analytical results for the spectral density and kernels for the
scalar-scalar, tensor-tensor and scalar-tensor contribution can be found in
Sec. 4.1.1, 4.1.2 and Sec. 4.1.3 respectively.

In a $\Lambda$CDM universe, SOGWs sourced during the radiation dominated epoch
of universe have a distinctive signal compared to the SIGWs. Due to the
behaviour of source terms including tensors in the UV limit, we studied the
spectral density of waves induced by peaked spectra. Particularly, we showed
if there is peak in the first order tensor power spectrum, then the overall
spectral density will be different to waves induced by quadratic scalar
fluctuations. First of all, in the case of the primordial power spectra having
the same amplitude, the scalar-tensor contribution becomes dominant in the
power spectrum and if the primordial spectra peak at the same scale the
distinctive resonant behaviour of SIGWs is not apparent anymore. We also
showed that when the primordial tensor power spectrum peaks at smaller scale
then the scalar one, the spectral density gets an additional contribution
which shows up as an additional peak (see Fig. 7). Since the power spectrum
for first order scalars and tensors is directly linked to inflation, our
results can be used to constrain certain models which have a peak in the
primordial tensor power spectrum.

Work on second order gravitational waves, sourced by first order scalar and
tensor modes has so far focused mainly on toy model input power spectra. A
natural next step is therefore to use more realistic input power spectra, both
for the scalar and the tensor perturbations. In most cases this will no longer
allow the use of analytical techniques and will require the use of numerical
methods. Luckily there are already numerical packages to calculate the power
spectrum for any single or multi-field model available, such as PyTransport
[47]. In a future extension of the work presented here we therefore plan to
combine our results with PyTransport, which will allow us to further constrain
the parameter space of models of inflation using gravitational waves.

##  Appendix A Evolution equations

Below we give the background, first order and second order relevant Einstein
equations. We have used the $xPand$ package [45] for our computations. All
calculations carried out below are done in the Newtonian gauge.

###  A.1 Background equations

The energy-momentum tensor is given in Eq. (2.4). At $0^{th}$ order
(background), the two Friedmann equations are

$\mathcal{H}^{2}=\frac{a^{2}}{3M_{Pl}^{2}}\rho,\quad\mathcal{H}^{\prime}=-\frac{a^{2}\rho}{6M_{Pl}^{2}}(1+3w),\quad\text{with}\quad
w=\frac{P}{\rho}\text{ ,}$ (A.1)

where $w$ is the equation of state parameter. There is also a conservation
equation

$\rho^{\prime}+3\mathcal{H}\left(\rho+P\right)=0\text{ .}$ (A.2)

From the above we infer that during RD ($w=1/3$):

$a(\eta)\propto\eta,\quad\mathcal{H}=\frac{1}{\eta}\text{ .}$ (A.3)

###  A.2 First order equations

The source terms in Eq. (2.7) have been simplified using first order relations
we obtain from the Einstein equations, below we present only the relevant
ones. The first order density perturbation and the scalar part of the velocity
perturbation can be related to metric variables, from the $0-0$ and by taking
the divergence of the $0-i$ first order Einstein equations respectively, we
arrive to

$\displaystyle\delta\rho$ $\displaystyle=$
$\displaystyle-2\rho\Phi+\frac{M_{Pl}^{2}}{a^{2}}\left(-6\mathcal{H}\Psi^{\prime}+2\nabla^{2}\Psi\right)\text{
,}$ (A.4) $\displaystyle v$ $\displaystyle=$
$\displaystyle-\frac{2}{3(1+w)\mathcal{H}}\left(\Phi+\frac{\Psi^{\prime}}{\mathcal{H}}\right)\text{
.}$ (A.5)

Furthermore, in the case of vanishing anisotropic stress, the symmetric trace-
free part of the Einstein equations for a perfect fluid implies $\Phi=\Psi$.
This simplifies the calculations significantly, namely we can get an evolution
equation for $\Psi$ (the trace of the spatial part of the Einstein equations)
in Fourier space:

$\Psi^{\prime\prime}+3\mathcal{H}\left(1+c_{s}^{2}\right)\Psi^{\prime}+c_{s}^{2}k^{2}\Psi=0\text{
.}$ (A.6)

The solution to this equation is shown in Eq. (2.12). We have introduced the
speed of sound:

$c_{s}^{2}=\frac{\delta P}{\delta\rho}\text{ ,}$ (A.7)

for adiabatic perturbations. We further note that when the equation of state
is constant, then $c_{s}^{2}=w$, as is the case in our work.

Additionally, we are interested in the first order tensor evolution equation,
using the TT projector, defined in Eq. (2.2), on the first order spatial
Einstein equations we arrive at

$h_{ij}^{\prime\prime}+2\mathcal{H}h_{ij}^{{}^{\prime}}-\nabla^{2}h_{ij}=0\text{
.}$ (A.8)

The solution to the first order tensor mode equation is shown in Eq. (2.13).

###  A.3 Second order equations

The spatial part of the second order Einstein tensor is given by:

$\begin{split}G_{ij}^{(2)}&=h_{ij}^{\prime\prime(2)}+2\mathcal{H}h_{ij}^{\prime(2)}-h_{ij}^{(2)}(2\mathcal{H}^{2}+4\mathcal{H}^{\prime})-\nabla^{2}h_{ij}^{(2)}+\delta_{ij}(\Phi^{(2)}+\Psi^{(2)})(2\mathcal{H}^{2}+4\mathcal{H}^{\prime})\\\
&+2\delta_{ij}\Phi^{\prime(2)}\mathcal{H}+4\delta_{ij}\Psi^{\prime(2)}\mathcal{H}+2\delta_{ij}\Psi^{\prime\prime(2)}+\delta_{ij}\nabla^{2}\Phi^{(2)}-\delta_{ij}\nabla^{2}\Psi^{(2)}\\\
&-\partial_{i}\partial_{j}\Phi^{(2)}+\partial_{i}\partial_{j}\Psi^{(2)}-4h_{i}^{m\prime}h_{jm}^{\prime}+\delta_{ij}\left(3h_{mk}^{\prime}h^{mk\prime}+4h^{mk}h_{mk}^{\prime\prime}+8\mathcal{H}h^{mk}h_{mk}^{\prime}\right)\\\
&-\Phi\left(4h_{ij}^{\prime\prime}+8\mathcal{H}h_{ij}^{\prime}-8\mathcal{H}^{2}h_{ij}-16\mathcal{H}^{\prime}h_{ij}\right)-\delta_{ij}\left(8\mathcal{H}^{2}+16\mathcal{H}^{\prime}\right)\Phi^{(2)}-16\delta_{ij}\mathcal{H}\Phi\Phi^{\prime}\\\
&+\Phi^{\prime}\left(8\mathcal{H}h_{ij}-2h_{ij}^{\prime}\right)+\delta_{ij}\left(16\mathcal{H}^{\prime}-8\mathcal{H}^{2}\right)\Phi\Psi-8\delta_{ij}\mathcal{H}\Psi\Phi^{\prime}+2h_{ij}^{\prime}\Psi^{\prime}+24\mathcal{H}h_{ij}\Psi^{\prime}\\\
&-16\delta_{ij}\mathcal{H}\Phi\Psi^{\prime}-4\delta_{ij}\Phi^{\prime}\Psi^{\prime}+2\delta_{ij}\Psi^{\prime
2}+12h_{ij}\Psi^{\prime\prime}-8\delta_{ij}\Phi\Psi^{\prime\prime}-4\Psi\nabla^{2}h_{ij}+4h_{i}^{m}\partial_{m}\partial_{j}\Psi\\\
&+\left(4h_{ij}-4\delta_{ij}\Phi\right)\nabla^{2}\Phi-\left(8h_{ij}+4\delta_{ij}\Psi\right)\nabla^{2}\Psi+4h_{j}^{m}\partial_{m}\partial_{i}\Psi-2\partial_{k}h_{ij}\partial^{k}\Phi\\\
&-2\delta_{ij}\partial_{k}\Phi\partial^{k}\Phi-6\partial_{k}h_{ij}\partial^{k}\Psi-4\delta_{ij}\partial_{k}\Psi\partial^{k}\Psi+4h^{mk}\partial_{m}\partial_{k}h_{ij}-4h^{mk}\delta_{ij}\partial_{m}\partial_{k}\Phi\\\
&-4\partial_{m}h_{jk}\partial^{k}h_{i}^{m}+4\partial_{k}h_{jm}\partial^{k}h_{i}^{m}-4h^{mk}\delta_{ij}\nabla^{2}h_{mk}+2\delta_{ij}\partial_{k}h_{md}\partial^{d}h^{mk}\\\
&-3\delta_{ij}\partial_{d}h_{mk}\partial^{d}h^{mk}+2\partial^{m}\Phi\partial_{i}h_{jm}+2\partial^{m}\Psi\partial_{i}h_{jm}-4h^{mk}\partial_{i}\partial_{k}h_{jm}+2\partial_{i}h^{mk}\partial_{j}h_{mk}\\\
&+2\partial^{m}\Psi\partial_{j}h_{im}+2\partial^{m}\Phi\partial_{j}h_{im}+2\partial_{i}\Phi\partial_{j}\Phi-2\partial_{i}\Psi\partial_{j}\Phi-2\partial_{i}\Phi\partial_{j}\Psi+6\partial_{i}\Psi\partial_{j}\Psi\\\
&-4h^{mk}\partial_{j}\partial_{k}h_{im}+4h^{mk}\partial_{j}\partial_{i}h_{mk}+4\Phi\partial
j\partial_{i}\Phi+4\Psi\partial_{j}\partial_{i}\Psi\end{split}\text{ ,}$ (A.9)

where we have not used the background and first order equations for
simplification yet. Similarly, the energy momentum tensor is given by:

$\displaystyle T_{ij}^{(2)}$ $\displaystyle=$ $\displaystyle
a^{2}\bigg{(}\delta_{ij}\left[\frac{1}{2}P^{(2)}-2c_{s}^{2}\delta\rho\Psi-w\rho\Psi^{(2)}\right]+2c_{s}^{2}\delta\rho
h_{ij}+(1+w)\rho\partial_{i}v\partial_{j}v$ (A.10) $\displaystyle+$
$\displaystyle w\rho h_{ij}^{(2)}\bigg{)}\text{ .}$

##  Appendix B Fourier space conventions

In this Appendix we present our choice of conventions in Fourier space in
which we carry out our work.

A scalar quantity, $S(\mathbf{x},\eta)$, can be expressed as a Fourier
integral

$S(\mathbf{x},\eta)=\frac{1}{2\pi^{\frac{3}{2}}}\int
d^{3}\mathbf{k}S(\mathbf{k},\eta)e^{i\mathbf{k}\cdot\mathbf{x}}\text{ .}$
(B.1)

Similarly, a symmetric transverse-traceless tensor, $T_{ab}(\mathbf{x},\eta)$,
in real space can be expressed as a Fourier integral

$T_{ab}(\mathbf{x},\eta)=\frac{1}{2\pi^{\frac{3}{2}}}\int
d^{3}\mathbf{k}\left\\{T(\mathbf{k},\eta)q_{ab}(\mathbf{k})+\overline{T}(\mathbf{k},\eta)\overline{q}_{ab}(\mathbf{k})\right\\}e^{i\mathbf{k}\cdot\mathbf{x}}\text{
,}$ (B.2)

where we note that the two polarizations decouple. The tensor basis is given
in terms of two time-independent polarization tensors $q_{ab}(\mathbf{k})$ and
$\overline{q}_{ab}(\mathbf{k})$ which, respectively, are the ‘plus’ and
‘cross’ polarizations of a gravitational wave. These tensors can be expressed
in terms of vectors that are orthogonal to the direction of propagation of the
gravitational wave, $\mathbf{k}$:

$\begin{split}q_{ab}^{+}(\mathbf{k})=\frac{1}{\sqrt{2}}\left(e_{a}(\mathbf{k})e_{b}(\mathbf{k})-\overline{e}_{a}(\mathbf{k})\overline{e}_{b}(\mathbf{k})\right)\\\
q_{ab}^{\times}(\mathbf{k})=\frac{1}{\sqrt{2}}\left(e_{a}(\mathbf{k})\overline{e}_{b}(\mathbf{k})+\overline{e}_{a}(\mathbf{k})e_{b}(\mathbf{k})\right)\end{split}$
(B.3)

The basis vectors obey:

$\displaystyle e_{a}e^{a}=\overline{e}_{a}\overline{e}^{a}=1,\quad
e_{a}\overline{e}^{a}=0,\quad e_{a}k^{a}=\overline{e}_{a}k^{a}=0\text{ ,}$
(B.4)

and hence the symmetric tensor basis satisfies ($\lambda=+,\times$):

$q_{ab}^{\lambda}(\mathbf{k})\delta^{ab}=0,\quad
q_{ab}^{\lambda}(\mathbf{k})k^{a}=0,\quad
q_{ab}^{\lambda}(\mathbf{k})q^{ab,\lambda^{\prime}}(\mathbf{k})=\delta^{\lambda\lambda^{\prime}}\text{
.}$ (B.5)

Additionally $q_{ij}^{\lambda}(\mathbf{k})$ can be taken to be real, so that:
$q_{ij}^{\lambda}(\mathbf{k})=q_{ij}^{\lambda}(\mathbf{-k})$.

It is useful to parameterize the vectors and basis tensors in spherical
coordinates, since FLRW spacetime is spherically symmetric. The wave-vector
$\mathbf{k}$ is expressed as:

$\mathbf{k}=k\left(\sin{\theta_{k}}\cos{\phi_{k}},\sin{\theta_{k}}\sin{\phi_{k}},\cos{\theta_{k}}\right)\text{
,}$ (B.6)

in the $(k,\theta,\phi)$ coordinate system, with $\phi\in[0;2\pi[$ and
$\theta\in[0;\pi]$. We can construct two orthonormal vectors in the subspace
perpendicular to $\mathbf{k}$:

$\begin{split}&\mathbf{e}(\mathbf{k})=\left(\cos{\theta_{k}}\cos{\phi_{k}},\cos{\theta_{k}}\sin{\phi_{k}},-\sin{\theta_{k}}\right)\text{
,}\\\
&\overline{\mathbf{e}}(\mathbf{k})=\left(-\sin{\phi_{k}},\cos{\phi_{k}},0\right)\text{
.}\end{split}$ (B.7)

Furthermore, due to isotropy we have the freedom to align $\mathbf{k}$ with
the z-$axis$, therefore $\theta_{k}=\phi_{k}=0$. Similarly, for the the vector
$\mathbf{p}$:

$\mathbf{p}=p\left(\sin{\theta_{p}}\cos{\phi_{p}},\sin{\theta_{p}}\sin{\phi_{p}},\cos{\theta_{p}}\right)\text{
,}$ (B.8)

from which we get the important result:

$|\mathbf{k}-\mathbf{p}|^{2}=k^{2}+p^{2}-2kp\cos{\theta}_{p}\text{ .}$ (B.9)

##  Appendix C Polarization functions

In this Appendix we show that the convolution integral is invariant under the
switching of the two variables of the convolution and provide the results for
the polarization function integrals in the azimuthal direction.

###  C.1 Convolution symmetry in the tensor-tensor sector

As stated previously, due to he nature of the second order computation (i.e. a
convolution) there is a symmetry between $\mathbf{p}$ and
$\mathbf{k}-\mathbf{p}$ that can be exploited to simplify computations when
extracting the power spectrum for the scalar-scalar and tensor-tensor induced
waves. This has been shown in Ref. [40] and for the case of scalar-scalar
induced waves in Ref. [37]. Here we extend this reasoning for convoluted terms
quadratic in tensor perturbations.

When we substitute Eq. (3.12) into Eq. (3.11), we use the fact that
integrating over $\delta(\mathbf{p}+\mathbf{p^{\prime}})$ or
$\delta(\mathbf{p}+\mathbf{k^{\prime}}-\mathbf{p^{\prime}})$ is the same since
there is a symmetry originating from the convolution. Consider a term that
appears in Eq. (2.14), for example:

$\int
d^{3}\mathbf{p}q^{ab}_{\lambda}(\mathbf{k})q^{m,\lambda_{1}}_{a}(\mathbf{k}-\mathbf{p})q^{\lambda_{2}}_{bm}(\mathbf{p})f(p)f(k-p)$
(C.1)

We want to show that the integral is unchanged by switching $\mathbf{p}$ and
$\mathbf{k}-\mathbf{p}$. Let us consider the change of variable
$\mathbf{\tilde{p}}=\mathbf{k}-\mathbf{p}$, then the above changes to:

$\int
d^{3}\mathbf{\tilde{p}}q^{ab}_{\lambda}(\mathbf{k})q^{m,\lambda_{1}}_{a}(\mathbf{\tilde{p}})q^{\lambda_{2}}_{bm}(\mathbf{k}-\mathbf{\tilde{p}})f(k-\tilde{p})f(\tilde{p})$
(C.2)

We have the freedom to relabel $a\leftrightarrow b$ and using the fact that
$q^{ab}(\mathbf{k})$ is symmetric we have:

$\displaystyle\int
d^{3}\mathbf{\tilde{p}}q^{ab}_{\lambda}(\mathbf{k})q^{m,\lambda_{1}}_{b}(\mathbf{\tilde{p}})q^{\lambda_{2}}_{am}(\mathbf{k}-\mathbf{\tilde{p}})f(k-\tilde{p})f(\tilde{p})$
(C.3)

Finally, remembering that there is an implicit sum over $\lambda_{1}$ and
$\lambda_{2}$ (i.e. they are just labels), we can relabel them, so that:

$\int
d^{3}\mathbf{\tilde{p}}q^{ab}_{\lambda}(\mathbf{k})q^{m,\lambda_{1}}_{a}(\mathbf{k}-\mathbf{\tilde{p}})q^{\lambda_{2}}_{bm}(\mathbf{\tilde{p}})f(\tilde{p})f(k-\tilde{p})$
(C.4)

This reasoning can be extended to all the tensor-tensor terms that appear in
Eq. (2.15).

###  C.2 Integral evaluations of the polarization functions

After switching to spherical coordinates, the polarization functions in Eq.
(2.15) are evaluated to be:

$q^{\lambda,ss}(k,p)=\begin{cases}\frac{1}{\sqrt{2}}p^{2}\sin^{2}{\theta}\cos{2\phi}\quad&\lambda=+\\\
\frac{1}{\sqrt{2}}p^{2}\sin^{2}{\theta}\sin{2\phi}\quad&\lambda=\times\end{cases}$
(C.5)
$q^{\lambda,tt}_{1}(k,p)=\begin{cases}\frac{1}{\sqrt{2}}p^{2}\sin^{2}{\theta}\cos{2\phi}\quad&\lambda=+\\\
\frac{1}{\sqrt{2}}p^{2}\sin^{2}{\theta}\sin{2\phi}\quad&\lambda=\times\end{cases}$
(C.6)
$q^{\lambda,tt}_{2}(k,p)=\begin{cases}-\frac{1}{\sqrt{2}}\sin^{2}{\theta}\cos{2\phi}\quad&\lambda=+\\\
-\frac{1}{\sqrt{2}}\sin^{2}{\theta}\sin{2\phi}\quad&\lambda=\times\end{cases}$
(C.7) $\displaystyle
q^{\lambda,st}(k,p)=\begin{cases}\frac{1}{2}\left((1+\cos{\theta_{p}}^{2})(\cos{\phi_{p}}^{2}-\sin{\phi_{p}}^{2})-4\cos{\theta_{p}}\cos{\phi_{p}}\sin{\phi_{p}}\right)\quad&\lambda=+\\\
\cos{\theta_{p}}\cos{2\phi_{p}}+\cos{\phi_{p}}\sin{\phi_{p}}(1+\cos{\theta_{p}}^{2})\quad&\lambda=\times\end{cases}$
(C.8)

The expressions that contribute to the power spectrum are listed below (once
integrated over the azimuthal angle). For the pure scalar part:

$\displaystyle\int_{0}^{2\pi}q^{\lambda,ss}q^{\lambda^{\prime},ss}$
$\displaystyle=$
$\displaystyle\frac{1}{2}p^{4}\pi\sin^{4}\theta_{p}\delta^{\lambda\lambda^{\prime}}$
(C.9) $\displaystyle=$ $\displaystyle\frac{1}{2}k^{4}\pi
v^{4}\left(1-\frac{(1+v^{2}-u^{2})^{2}}{4v^{2}}\right)^{2}\delta^{\lambda\lambda^{\prime}}\text{
,}$

where in the last equality we used:

$\sin^{2}{\theta_{p}}=1-\frac{(1+v^{2}-u^{2})^{2}}{4v^{2}}=1-\cos^{2}{\theta_{p}}\text{
.}$ (C.10)

For the pure tensor part, the resulting expressions are:

$\displaystyle\int_{0}^{2\pi}q^{\lambda,tt}_{1}q^{\lambda^{\prime},tt}_{1}$
$\displaystyle=$ $\displaystyle\frac{1}{2}k^{4}\pi
v^{4}\left(1-\frac{(1+v^{2}-u^{2})^{2}}{4v^{2}}\right)^{2}\delta^{\lambda\lambda^{\prime}}\text{
,}$ (C.11)
$\displaystyle\int_{0}^{2\pi}q^{\lambda,tt}_{1}q^{\lambda^{\prime},tt}_{2}$
$\displaystyle=$ $\displaystyle-\frac{1}{2}\pi
k^{2}v^{2}\left(1-\frac{\left(v^{2}-u^{2}+1\right)^{2}}{4v^{2}}\right)^{2}\delta^{\lambda\lambda^{\prime}}\text{
,}$ (C.12)
$\displaystyle\int_{0}^{2\pi}q^{\lambda,tt}_{2}q^{\lambda^{\prime},tt}_{2}$
$\displaystyle=$
$\displaystyle\frac{1}{2}\pi\left(1-\frac{\left(v^{2}-u^{2}+1\right)^{2}}{4v^{2}}\right)^{2}\delta^{\lambda\lambda^{\prime}}\text{
.}$ (C.13)

Finally the cross term contributes as:

$\displaystyle\int_{0}^{2\pi}q^{\lambda,st}q^{\lambda^{\prime},st}$
$\displaystyle=$
$\displaystyle\frac{\pi}{64v^{4}}\bigg{(}(u^{2}-1)^{4}-4v^{2}(u^{2}-7)(u^{2}-1)^{2}+2v^{4}(3u^{4}-30u^{2}+35)+v^{8}$
(C.14) $\displaystyle-$ $\displaystyle
4v^{6}(u^{2}-7)\bigg{)}\delta^{\lambda\lambda^{\prime}}$ $\displaystyle\equiv$
$\displaystyle\pi M_{st}(v,u)\delta^{\lambda\lambda^{\prime}}\text{ .}$

##  Appendix D Evaluation of kernels

###  D.1 Definition of the functions

In Eq. (2.16) we defined several functions. Here they are:

$\begin{split}f^{ss}(\eta,k,|\mathbf{k}-\mathbf{p}|)&=\frac{2}{3(1+w)}\left[\left((T_{\Psi}(\eta|\mathbf{k}-\mathbf{p}|)+\frac{T^{\prime}_{\Psi}(\eta|\mathbf{k}-\mathbf{p}|)}{\mathcal{H}}\right)\left(T_{\Psi}(\eta
p)+\frac{T^{\prime}_{\Psi}(\eta p)}{\mathcal{H}}\right)\right]\\\
&+T_{\Psi}(\eta|\mathbf{k}-\mathbf{p}|)T_{\Psi}(\eta p)\text{ ,}\end{split}$
(D.1)

Tensor terms:

$\displaystyle f^{tt}_{1}(\eta,k,|\mathbf{k}-\mathbf{p}|)$
$\displaystyle=T_{E}(\eta|\mathbf{k}-\mathbf{p}|)T_{h}(\eta p)\text{ ,}$
(D.2a) $\displaystyle f^{tt}_{2}(\eta,k,|\mathbf{k}-\mathbf{p}|)$
$\displaystyle=T^{\prime}_{E}(\eta|\mathbf{k}-\mathbf{p}|)T^{\prime}_{E}(\eta
p)\text{ ,}$ (D.2b)

Cross term:

$\displaystyle f^{st}(\eta p,c_{s}\eta|\mathbf{k}-\mathbf{p}|)$
$\displaystyle=$ $\displaystyle-2p^{2}T_{h}(\eta
p)T_{\Psi}(c_{s}\eta|\mathbf{k}-\mathbf{p}|)-2(k-p)pT_{h}(\eta
p)T_{\Psi}(c_{s}\eta|\mathbf{k}-\mathbf{p}|)$ (D.3) $\displaystyle+$
$\displaystyle\mathcal{H}(1+3c_{s}^{2})T_{h}(\eta
p)T_{\Psi}^{\prime}(c_{s}\eta|\mathbf{k}-\mathbf{p}|)-(k-p)^{2}(1-c_{s}^{2})T_{h}(\eta
p)$ $\displaystyle\times$ $\displaystyle
T_{\Psi}(c_{s}\eta|\mathbf{k}-\mathbf{p}|)\text{ .}$

We can write everything in terms of the variables $v$ and $u$ in Eq. (3.5).
The scalar-scalar terms:

$\begin{split}f^{ss}(x,u,v,c_{s})&=\frac{2}{3(1+w)}\left[\left((T_{\Psi}(c_{s}xv)+x\frac{\partial
T_{\Psi}(c_{s}xv)}{\partial x}\right)\left(T_{\Psi}(c_{s}xu)+x\frac{\partial
T_{\Psi}(c_{s}xu)}{\partial x}\right)\right]\\\
&+T_{\Psi}(c_{s}xv)T_{\Psi}(c_{s}xu)\text{ ,}\end{split}$ (D.4)

Tensor-tensor terms:

$\displaystyle f^{tt}_{1}(x,v,u)$ $\displaystyle=T_{h}(xv)T_{h}(xu)\text{ ,}$
(D.5a) $\displaystyle f^{tt}_{2}(x,v,u)$ $\displaystyle=k^{2}\frac{\partial
T_{h}(xv)}{\partial x}\frac{\partial T_{h}(xu)}{\partial x}\text{ ,}$ (D.5b)

and finally the cross term:

$\displaystyle f^{st}(x,u,v,c_{s})$ $\displaystyle=$
$\displaystyle-2v^{2}k^{2}T_{h}(xv)T_{\Psi}(c_{s}xu)-2vuk^{2}T_{h}(xv)T_{\Psi}(c_{s}xu)$
$\displaystyle+$
$\displaystyle(1+3c_{s}^{2})\frac{k^{2}}{x}T_{h}(xv)\frac{\partial
T_{\Psi}(c_{s}ux)}{\partial
x}-u^{2}k^{2}(1-c_{s}^{2})T_{h}(xv)T_{\Psi}(c_{s}xu)\text{ .}$

###  D.2 Kernel results

All the Kernels presented in Eq. (3.15) and Eq. (3.21) can be evaluated
analytically and our presented below. These are trigonometric integrals, so it
is useful to recall the following definitions:

$Si(z)=\int_{0}^{z}\frac{\sin t}{t}dt\qquad Ci(z)=-\int_{z}^{\infty}\frac{\cos
t}{t}dt$

and their finite limits: $\lim_{x\to\infty}Si(z)=\pm\frac{\pi}{2}$ and
$\lim_{x\to\infty}Ci(z)=0$.

First we give the result from [38] for the scalar-scalar Kernel,
$c_{s}=1/\sqrt{3}$:

$\displaystyle I^{ss}(v,u,x\rightarrow\infty)$ $\displaystyle=$
$\displaystyle\frac{3(u^{2}+v^{2}-3)}{4u^{3}v^{3}x}\bigg{[}\sin{x}\left(-4uv+(u^{2}+v^{2}-3)\log\bigg{|}\frac{3-(u+v)^{2}}{3-(u-v)^{2}}\bigg{|}\right)$
(D.7)
$\displaystyle-\pi\cos{x}(u^{2}+v^{2}-3)\Theta(v+u-\sqrt{3})\bigg{]}\text{ .}$

The evaluation of the tensor-tensor kernels are presented below. For
concreteness, we show the steps for $\mathcal{I}_{c,1}^{tt}(x,v,u)$ and
$\mathcal{I}_{s,1}^{tt}(x,v,u)$ in Eq. (3.15):

$\displaystyle\mathcal{I}_{c,1}^{tt}(x,v,u)$ $\displaystyle=$
$\displaystyle-\frac{1}{4uv}\bigg{[}\text{Si}((u-v+1)x)+\text{Si}((-u+v+1)x)+\text{Si}((u+v-1)x)$
(D.8) $\displaystyle-\text{Si}((u+v+1)x)\bigg{]}\text{ ,}$
$\displaystyle\mathcal{I}_{s,1}^{tt}(x,v,u)$ $\displaystyle=$
$\displaystyle\frac{1}{4uv}\bigg{[}\text{Ci}(|u-v-1|x)+\text{Ci}(|u-v+1|x)-\text{Ci}(|u+v+1|)x)$
(D.9)
$\displaystyle-\text{Ci}(|u+v-1|x)-\log\bigg{|}\frac{(u-v-1)(u-v+1)}{(u+v+1)(u+v-1)}\bigg{|}\bigg{]}\text{
.}$

Since we are interested in modes that are deep inside the horizon we take the
limit where $k\eta>>1$, or equivalently sending $x\rightarrow\infty$.

$\displaystyle\mathcal{I}_{c,1}^{tt}(x\rightarrow\infty,v,u)=-\frac{\pi}{4uv}\Theta(u+v-1)\text{
,}$ (D.10)
$\displaystyle\mathcal{I}_{s,1}^{tt}(x\rightarrow\infty,v,u)=\frac{1}{4uv}\log\bigg{|}\frac{1-(u+v)^{2}}{1-(u-v)^{2}}\bigg{|}\text{
.}$ (D.11)

Similarly for $\mathcal{I}_{c,2}^{tt}(x,v,u)$ and
$\mathcal{I}_{s,2}^{tt}(x,v,u)$:

$\displaystyle\mathcal{I}_{c,2}^{tt}(x,v,u)$ $\displaystyle=$
$\displaystyle\frac{k^{2}}{2uvx^{2}}[\sin{x}\sin{ux}\sin{vx}]+\frac{4k^{2}}{8uvx}[-v\sin{x}\sin{ux}\cos{xx}$
(D.12)
$\displaystyle+\cos{x}\sin{ux}\sin{vx}-u\sin{x}\cos{ux}\sin{vx}]-\frac{k^{2}\left(u^{2}+v^{2}-1\right)}{8uv}\bigg{[}$
$\displaystyle\text{Si}((u-v+1)x)+\text{Si}((-u+v+1)x)-\text{Si}((u+v+1)x)$
$\displaystyle-\text{Si}(-ux-vx+x)\bigg{]}\text{ ,}$
$\displaystyle\mathcal{I}_{c,2}^{tt}(x\rightarrow\infty,v,u)$ $\displaystyle=$
$\displaystyle-\frac{\pi
k^{2}\left(u^{2}+v^{2}-1\right)}{8uv}\Theta(u+v-1)\text{ .}$ (D.13)

and

$\displaystyle\mathcal{I}_{s,2}^{tt}(x,v,u)$ $\displaystyle=$
$\displaystyle\frac{k^{2}}{2uvx^{2}}(\cos{x}\sin{ux}\sin{vx})+\frac{k^{2}}{2uvx}(\sin{x}\sin{ux}\sin{vx}+v\cos{x}\sin{ux}\cos{vx}$
$\displaystyle+u\cos{x}\cos{ux}\sin{vx})-\frac{k^{2}}{8uv}\left(u^{2}+v^{2}-1\right)\bigg{(}\text{Ci}(|u+v-1|x)+\text{Ci}(|u+v+1|x)$
$\displaystyle-\text{Ci}(|u-v-1|x)-\text{Ci}(|u-v+1|)x+\log\bigg{|}\frac{(u-v-1)(u-v+1)}{(u+v-1)(u+v+1)}\bigg{|}\bigg{)}$
, (D.14) $\displaystyle\mathcal{I}_{s,2}^{tt}(x\rightarrow\infty,v,u)$
$\displaystyle=$
$\displaystyle\frac{k^{2}}{8uv}\bigg{(}(u^{2}+v^{2}-1)\log\bigg{|}\frac{1-(u+v)^{2}}{1-(u-v)^{2}}\bigg{|}-4vu\bigg{)}\text{
.}$ (D.15)

Below we present the results for the scalar-tensor kernels, which is the
evaluation of Eq. (3.21). These expressions are much more involved than the
tensor-tensor kernels, therefore we display only the part that will contribute
to the final expression, i.e. the limit where $x\rightarrow\infty$.
$\mathcal{I}_{c}^{st}(x\rightarrow\infty,v,u)$ reads:

$\displaystyle\mathcal{I}_{c}^{st}(x\rightarrow\infty,v,u)$ $\displaystyle=$
$\displaystyle\frac{3\pi
k^{2}}{32u^{3}vc_{s}^{3}}(1+c_{s}u-v)(-1+c_{s}u+v)i^{st}_{-}(v,u)\Theta(v+uc_{s}-1)\text{,}$
(D.16)

where for compactness we defined:

$\displaystyle i^{st}_{\pm}(v,u)$ $\displaystyle=$
$\displaystyle-1+4u^{2}-3c_{s}^{2}(u^{2}+(v\pm
1)^{2}-c_{s}^{2}u^{2})+8uv+v(\mp 2+7v)\text{ .}$ (D.17)

and $\mathcal{I}_{s}^{st}(x\rightarrow\infty,v,u)$:

$\displaystyle\mathcal{I}_{s}^{st}(x\rightarrow\infty,v,u)$ $\displaystyle=$
$\displaystyle\frac{3k^{2}}{32c_{s}^{3}u^{3}v}\bigg{(}i_{+}^{st}(v,u)\log\bigg{|}\frac{c_{s}u-v-1}{c_{s}u+v+1}\bigg{|}(c_{s}u-v-1)(c_{s}u+v+1)$
(D.18) $\displaystyle+$ $\displaystyle
i_{-}^{st}(v,u)\log\bigg{|}\frac{c_{s}u-v+1}{c_{s}u+v-1}\bigg{|}(c_{s}u-v+1)(c_{s}u+v-1)-\frac{4}{9}c_{s}uv\bigg{[}-9$
$\displaystyle+$ $\displaystyle
12u^{2}+15c^{4}u^{2}+3v(8u+7v)-c_{s}^{2}(7u^{2}+3(9+v^{2}))\bigg{]}\bigg{)}\text{
.}$

## Acknowledgments

The authors thank Chris Addis, Chris Clarkson, Charles Dalang, Laura Iaconni
and David Mulryne for useful discussions. The authors are also grateful to
Guillem Domènech for useful discussions and sharing an early draft of their
work. RP is funded by STFC grant ST/P000592/1 and KAM is supported in part by
STFC grant ST/T000341/1.

##### Note added.

Whilst in the final stages of this project, Ref. [34] and Ref. [35] were
published on the arXiv.

## References

  * [1] B. P. Abbott et al. [LIGO Scientific and Virgo], “Observation of Gravitational Waves from a Binary Black Hole Merger,” Phys. Rev. Lett. 116 (2016) no.6, 061102 doi:10.1103/PhysRevLett.116.061102 [arXiv:1602.03837 [gr-qc]].
  * [2] G. Agazie et al. [NANOGrav], “The NANOGrav 15 yr Data Set: Evidence for a Gravitational-wave Background,” Astrophys. J. Lett. 951 (2023) no.1, L8 doi:10.3847/2041-8213/acdac6 [arXiv:2306.16213 [astro-ph.HE]].
  * [3] J. Antoniadis, S. Babak, A. S. Bak Nielsen, C. G. Bassa, A. Berthereau, M. Bonetti, E. Bortolas, P. R. Brook, M. Burgay and R. N. Caballero, et al. “The second data release from the European Pulsar Timing Array I. The dataset and timing analysis,” doi:10.1051/0004-6361/202346841 [arXiv:2306.16224 [astro-ph.HE]].
  * [4] J. Antoniadis, P. Arumugam, S. Arumugam, S. Babak, M. Bagchi, A. S. B. Nielsen, C. G. Bassa, A. Bathula, A. Berthereau and M. Bonetti, et al. “The second data release from the European Pulsar Timing Array III. Search for gravitational wave signals,” [arXiv:2306.16214 [astro-ph.HE]].
  * [5] J. Antoniadis, P. Arumugam, S. Arumugam, P. Auclair, S. Babak, M. Bagchi, A. S. Bak Nielsen, E. Barausse, C. G. Bassa and A. Bathula, et al. “The second data release from the European Pulsar Timing Array: V. Implications for massive black holes, dark matter and the early Universe,” [arXiv:2306.16227 [astro-ph.CO]].
  * [6] H. Xu, S. Chen, Y. Guo, J. Jiang, B. Wang, J. Xu, Z. Xue, R. N. Caballero, J. Yuan and Y. Xu, et al. “Searching for the Nano-Hertz Stochastic Gravitational Wave Background with the Chinese Pulsar Timing Array Data Release I,” Res. Astron. Astrophys. 23 (2023) no.7, 075024 doi:10.1088/1674-4527/acdfa5 [arXiv:2306.16216 [astro-ph.HE]].
  * [7] D. J. Reardon, A. Zic, R. M. Shannon, G. B. Hobbs, M. Bailes, V. Di Marco, A. Kapur, A. F. Rogers, E. Thrane and J. Askew, et al. “Search for an Isotropic Gravitational-wave Background with the Parkes Pulsar Timing Array,” Astrophys. J. Lett. 951 (2023) no.1, L6 doi:10.3847/2041-8213/acdd02 [arXiv:2306.16215 [astro-ph.HE]].
  * [8] A. Zic, D. J. Reardon, A. Kapur, G. Hobbs, R. Mandow, M. Curyło, R. M. Shannon, J. Askew, M. Bailes and N. D. R. Bhat, et al. “The Parkes Pulsar Timing Array Third Data Release,” [arXiv:2306.16230 [astro-ph.HE]].
  * [9] A. Afzal et al. [NANOGrav], “The NANOGrav 15 yr Data Set: Search for Signals from New Physics,” Astrophys. J. Lett. 951 (2023) no.1, L11 doi:10.3847/2041-8213/acdc91 [arXiv:2306.16219 [astro-ph.HE]].
  * [10] A. R. Liddle and D. H. Lyth, “Cosmological inflation and large scale structure,” doi:10.1017/CBO9781139175180
  * [11] K. Tomita, “Non-Linear Theory of Gravitational Instability in the Expanding Universe”, Progress of Theoretical Physics 37 (1967) 831 [https://academic.oup.com/ptp/article-pdf/37/5/831/5234391/37-5-831.pdf]
  * [12] N. Aghanim et al. [Planck], “Planck 2018 results. VI. Cosmological parameters,” Astron. Astrophys. 641 (2020), A6 [erratum: Astron. Astrophys. 652 (2021), C4] doi:10.1051/0004-6361/201833910 [arXiv:1807.06209 [astro-ph.CO]].
  * [13] K. N. Ananda, C. Clarkson and D. Wands, “The Cosmological gravitational wave background from primordial density perturbations,” Phys. Rev. D 75 (2007), 123518 doi:10.1103/PhysRevD.75.123518 [arXiv:gr-qc/0612013 [gr-qc]].
  * [14] D. Baumann, P. J. Steinhardt, K. Takahashi and K. Ichiki, “Gravitational Wave Spectrum Induced by Primordial Scalar Perturbations,” Phys. Rev. D 76 (2007), 084019 doi:10.1103/PhysRevD.76.084019 [arXiv:hep-th/0703290 [hep-th]].
  * [15] G. Domènech, “Scalar Induced Gravitational Waves Review,” Universe 7 (2021) no.11, 398 doi:10.3390/universe7110398 [arXiv:2109.01398 [gr-qc]].
  * [16] G. Domènech, C. Lin and M. Sasaki, “Gravitational wave constraints on the primordial black hole dominated early universe,” JCAP 04 (2021), 062 [erratum: JCAP 11 (2021), E01] doi:10.1088/1475-7516/2021/11/E01 [arXiv:2012.08151 [gr-qc]].
  * [17] Z. Chang, X. Zhang and J. Z. Zhou, “Gravitational waves from primordial scalar and tensor perturbations,” [arXiv:2209.07693 [astro-ph.CO]].
  * [18] W. H. Kinney, “A Hamilton-Jacobi approach to nonslow roll inflation,” Phys. Rev. D 56 (1997), 2002-2009 doi:10.1103/PhysRevD.56.2002 [arXiv:hep-ph/9702427 [hep-ph]].
  * [19] J. Garcia-Bellido and E. Ruiz Morales, “Primordial black holes from single field models of inflation,” Phys. Dark Univ. 18 (2017), 47-54 doi:10.1016/j.dark.2017.09.007 [arXiv:1702.03901 [astro-ph.CO]].
  * [20] L. Iacconi, H. Assadullahi, M. Fasiello and D. Wands, “Revisiting small-scale fluctuations in $\alpha$-attractor models of inflation,” JCAP 06 (2022) no.06, 007 doi:10.1088/1475-7516/2022/06/007 [arXiv:2112.05092 [astro-ph.CO]].
  * [21] H. Motohashi and W. Hu, “Primordial Black Holes and Slow-Roll Violation,” Phys. Rev. D 96 (2017) no.6, 063503 doi:10.1103/PhysRevD.96.063503 [arXiv:1706.06784 [astro-ph.CO]].
  * [22] J. Fumagalli, S. Renaux-Petel, J. W. Ronayne and L. T. Witkowski, “Turning in the landscape: A new mechanism for generating primordial black holes,” Phys. Lett. B 841 (2023), 137921 doi:10.1016/j.physletb.2023.137921 [arXiv:2004.08369 [hep-th]].
  * [23] G. A. Palma, S. Sypsas and C. Zenteno, “Seeding primordial black holes in multifield inflation,” Phys. Rev. Lett. 125 (2020) no.12, 121301 doi:10.1103/PhysRevLett.125.121301 [arXiv:2004.06106 [astro-ph.CO]].
  * [24] M. Braglia, D. K. Hazra, F. Finelli, G. F. Smoot, L. Sriramkumar and A. A. Starobinsky, “Generating PBHs and small-scale GWs in two-field models of inflation,” JCAP 08 (2020), 001 doi:10.1088/1475-7516/2020/08/001 [arXiv:2005.02895 [astro-ph.CO]].
  * [25] N. Barnaby, E. Pajer and M. Peloso, “Gauge Field Production in Axion Inflation: Consequences for Monodromy, non-Gaussianity in the CMB, and Gravitational Waves at Interferometers,” Phys. Rev. D 85 (2012), 023525 doi:10.1103/PhysRevD.85.023525 [arXiv:1110.3327 [astro-ph.CO]].
  * [26] B. Thorne, T. Fujita, M. Hazumi, N. Katayama, E. Komatsu and M. Shiraishi, “Finding the chiral gravitational wave background of an axion-SU(2) inflationary model using CMB observations and laser interferometers,” Phys. Rev. D 97 (2018) no.4, 043506 doi:10.1103/PhysRevD.97.043506 [arXiv:1707.03240 [astro-ph.CO]].
  * [27] E. Dimastrogiovanni, M. Fasiello and T. Fujita, “Primordial Gravitational Waves from Axion-Gauge Fields Dynamics,” JCAP 01 (2017), 019 doi:10.1088/1475-7516/2017/01/019 [arXiv:1608.04216 [astro-ph.CO]].
  * [28] N. Bartolo, C. Caprini, V. Domcke, D. G. Figueroa, J. Garcia-Bellido, M. C. Guzzetti, M. Liguori, S. Matarrese, M. Peloso and A. Petiteau, et al. “Science with the space-based interferometer LISA. IV: Probing inflation with gravitational waves,” JCAP 12 (2016), 026 doi:10.1088/1475-7516/2016/12/026 [arXiv:1610.06481 [astro-ph.CO]].
  * [29] S. Kawamura, M. Ando, N. Seto, S. Sato, M. Musha, I. Kawano, J. Yokoyama, T. Tanaka, K. Ioka and T. Akutsu, et al. “Current status of space gravitational wave antenna DECIGO and B-DECIGO,” PTEP 2021 (2021) no.5, 05A105 doi:10.1093/ptep/ptab019 [arXiv:2006.13545 [gr-qc]].
  * [30] S. Hild, M. Abernathy, F. Acernese, P. Amaro-Seoane, N. Andersson, K. Arun, F. Barone, B. Barr, M. Barsuglia and M. Beker, et al. “Sensitivity Studies for Third-Generation Gravitational Wave Observatories,” Class. Quant. Grav. 28 (2011), 094013 doi:10.1088/0264-9381/28/9/094013 [arXiv:1012.0908 [gr-qc]].
  * [31] C. Yuan, Z. C. Chen and Q. G. Huang, “Log-dependent slope of scalar induced gravitational waves in the infrared regions,” Phys. Rev. D 101 (2020) no.4, 043019 doi:10.1103/PhysRevD.101.043019 [arXiv:1910.09099 [astro-ph.CO]].
  * [32] C. J. Moore, R. H. Cole and C. P. L. Berry, “Gravitational-wave sensitivity curves,” Class. Quant. Grav. 32 (2015) no.1, 015014 doi:10.1088/0264-9381/32/1/015014 [arXiv:1408.0740 [gr-qc]].
  * [33] J. O. Gong, ’‘Analytic Integral Solutions for Induced Gravitational Waves,” Astrophys. J. 925 (2022) no.1, 102 doi:10.3847/1538-4357/ac3a6c [arXiv:1909.12708 [gr-qc]].
  * [34] P. Bari, N. Bartolo, G. Domènech and S. Matarrese, “Gravitational waves induced by scalar-tensor mixing,” [arXiv:2307.05404 [astro-ph.CO]].
  * [35] Y. H. Yu and S. Wang, “Primordial Gravitational Waves Assisted by Cosmological Scalar Perturbations,” [arXiv:2303.03897 [astro-ph.CO]].
  * [36] P. Peter and J. P. Uzan, “Primordial Cosmology,” Oxford University Press, 2013, ISBN 978-0-19-966515-0, 978-0-19-920991-0
  * [37] J. R. Espinosa, D. Racco and A. Riotto, “A Cosmological Signature of the SM Higgs Instability: Gravitational Waves,” JCAP 09 (2018), 012 doi:10.1088/1475-7516/2018/09/012 [arXiv:1804.07732 [hep-ph]].
  * [38] K. Kohri and T. Terada, “Semianalytic calculation of gravitational wave spectrum nonlinearly induced from primordial curvature perturbations,” Phys. Rev. D 97 (2018) no.12, 123532 doi:10.1103/PhysRevD.97.123532 [arXiv:1804.08577 [gr-qc]].
  * [39] K. A. Malik and D. Wands, “Cosmological perturbations,” Phys. Rept. 475 (2009), 1-51 doi:10.1016/j.physrep.2009.03.001 [arXiv:0809.4944 [astro-ph]].
  * [40] K. A. Malik, “A not so short note on the Klein-Gordon equation at second order,” JCAP 03 (2007), 004 doi:10.1088/1475-7516/2007/03/004 [arXiv:astro-ph/0610864 [astro-ph]].
  * [41] S. Pi and M. Sasaki, “Gravitational Waves Induced by Scalar Perturbations with a Lognormal Peak,” JCAP 09 (2020), 037 doi:10.1088/1475-7516/2020/09/037 [arXiv:2005.12306 [gr-qc]].
  * [42] B. S. Sathyaprakash and B. F. Schutz, “Physics, Astrophysics and Cosmology with Gravitational Waves,” Living Rev. Rel. 12 (2009), 2 doi:10.12942/lrr-2009-2 [arXiv:0903.0338 [gr-qc]].
  * [43] P. Auclair et al. [LISA Cosmology Working Group], “Cosmology with the Laser Interferometer Space Antenna,” Living Rev. Rel. 26 (2023) no.1, 5 doi:10.1007/s41114-023-00045-2 [arXiv:2204.05434 [astro-ph.CO]].
  * [44] Y. Akrami et al. [Planck], “Planck 2018 results. X. Constraints on inflation,” Astron. Astrophys. 641 (2020), A10 doi:10.1051/0004-6361/201833887 [arXiv:1807.06211 [astro-ph.CO]].
  * [45] C. Pitrou, X. Roy and O. Umeh, “xPand: An algorithm for perturbing homogeneous cosmologies,” Class. Quant. Grav. 30 (2013), 165002 doi:10.1088/0264-9381/30/16/165002 [arXiv:1302.6174 [astro-ph.CO]].
  * [46] C. Caprini and D. G. Figueroa, “Cosmological Backgrounds of Gravitational Waves,” Class. Quant. Grav. 35 (2018) no.16, 163001 doi:10.1088/1361-6382/aac608 [arXiv:1801.04268 [astro-ph.CO]]..
  * [47] D. J. Mulryne and J. W. Ronayne, “PyTransport: A Python package for the calculation of inflationary correlation functions,” J. Open Source Softw. 3 (2018) no.23, 494 doi:10.21105/joss.00494 [arXiv:1609.00381 [astro-ph.CO]].