††journal: iScience

A topic review on probing primordial black hole dark matter with scalar
induced gravitational waves

Chen<EMAIL_ADDRESS>Qing-Guo
<EMAIL_ADDRESS>
1School of Physical Sciences, University of Chinese Academy of Sciences, No.
19A Yuquan Road, Beijing 100049, China

2CAS Key Laboratory of Theoretical Physics, Institute of Theoretical Physics,
Chinese Academy of Sciences, Beijing 100190, China

3Center for Gravitation and Cosmology, College of Physical Science and
Technology, Yangzhou University, Yangzhou 225009, China

## abstract

Primordial black holes (PBHs) might form from the collapse of over-densed
regions generated by large scalar curvature perturbations in the radiation
dominated era. Despite decades of various independent observations, the nature
of dark matter (DM) remains highly puzzling. Recently, PBH DM have aroused
interest since they provide an attracting explanation to the merger events of
binary black holes discovered by LIGO/VIRGO and may play an important role on
DM. During the formation of PBH, gravitational waves will be sourced by linear
scalar perturbations at second-order, known as the scalar-induced
gravitational waves (SIGWs), which provides a new way to hunt for PBH DM. This
topic review mainly focus on the physics about SIGWs accompanying the
formation of PBH DM.

###### Contents

  1. 1 Introduction
  2. 2 Formation of Primordial Black Holes
    1. 2.1 Primordial non-Gaussianities
    2. 2.2 The non-linear effects between density contrast and curvature perturbation
    3. 2.3 Window Function and Power spectrum
    4. 2.4 The formation criterion
  3. 3 Scalar-induced gravitational waves
    1. 3.1 Scalar Induced Gravitational Waves in a General Cosmological Background
    2. 3.2 Higher-Order Corrections to Scalar Induced Gravitational Waves
    3. 3.3 Scalar Induced Gravitational Waves within Non-Gaussianities
    4. 3.4 Gauge issue of Scalar Induced Gravitational Waves
  4. 4 Searching for Primordial Black Hole Dark Matter Using Scalar Induced Gravitational Waves
  5. 5 Summary and Outlook

##  1 Introduction

Primordial black holes (PBHs) have aroused interest recently, not only because
they can represent the dark matter (DM) in our Universe, but also can explain
the binary black hole mergers events Sasaki  _et al._ (2016); Chen and Huang
(2018); Raidal  _et al._ (2019); De Luca  _et al._ (2020a); Hall  _et al._
(2020); Bhagwat  _et al._ (2021); Hütsi  _et al._ (2020) discovered by LIGO
Abbott  _et al._ (2016) if $\sim 10^{-3}$ DM is in the form of PBHs. So far,
there is no evidence for PBHs. Although various observations have constrained
the fraction of DM in the form of PBHs Tisserand  _et al._ (2007); Carr  _et
al._ (2010); Barnacka  _et al._ (2012); Griest  _et al._ (2013); Graham  _et
al._ (2015); Brandt (2016); Chen  _et al._ (2016); Wang  _et al._ (2018);
Gaggero  _et al._ (2017); Ali-Haimoud and Kamionkowski (2017); Aloni  _et al._
(2017); Horowitz (2016); Niikura  _et al._ (2019a); Zumalacarregui and Seljak
(2018); Nakama  _et al._ (2018); Abbott  _et al._ (2018); Magee  _et al._
(2018); Chen  _et al._ (2019a); Niikura  _et al._ (2019b); Chen and Huang
(2019); Abbott  _et al._ (2019); Wang  _et al._ (2019a), $f_{\mathrm{pbh}}$,
there still exist an open window in the mass range of
$[10^{-16},10^{-14}]\cup[10^{-13},10^{-12}]M_{\odot}$ where PBHs are possible
to present all the DM in our Universe.

PBH is an old conception and it can date back to 1974 when Hawking and Carr
proposed that black holes can be generated due to the collapse of over-densed
regions in the early universe Carr and Hawking (1974); Carr (1975). The
formation of PBHs is a threshold process. Once scalar perturbations exceed a
critical value, they would generate an over-densed region which would
immediately undergo gravitational collapse to form a single PBH when the
comoving size of such region is of the order of the horizon size. The exact
calculation of the PBH mass function, $\beta$, which describes the mass
fraction of the Universe contained within PBHs at the formation time is still
a debating and complicated question by today.

Among all the constraints on PBH DM, the scalar-induced gravitational waves
(SIGWs) provide a quite stringent constraint which can be several orders of
magnitude better than the other constraints Chen  _et al._ (2019b) in a
certain mass range of PBHs. During radiation dominant (RD) epoch, scalar
perturbations will alter the quadrupolar moment of the radiation and thus emit
GWs at second-order Tomita (1967); Matarrese  _et al._ (1993, 1994, 1998); Noh
and Hwang (2004); Carbone and Matarrese (2005); Nakamura (2007). Therefore,
SIGWs were inevitably generated during the formation of PBHs, providing a
powerful tool to hunt for PBH DM. Moreover, PBHs are generated by large
perturbations at small scales much larger than those on CMB scales, the
second-order GWs induced by the enhanced perturbations sourced by the linear
perturbations may exceed the first-order tensor inflationary modes Saito and
Yokoyama (2009). See more relevant studies for SIGWs in Ananda  _et al._
(2007); Baumann  _et al._ (2007); Saito and Yokoyama (2009); Arroja  _et al._
(2009); Assadullahi and Wands (2010); Bugaev and Klimai (2010a, b); Saito and
Yokoyama (2010); Bugaev and Klimai (2011); Alabidi  _et al._ (2013); Nakama
and Suyama (2016); Nakama  _et al._ (2017); Inomata  _et al._ (2017a);
Orlofsky  _et al._ (2017); Garcia-Bellido  _et al._ (2017); Sasaki  _et al._
(2018); Espinosa  _et al._ (2018); Kohri and Terada (2018a); Cai  _et al._
(2019a); Bartolo  _et al._ (2019a, b); Unal (2019); Byrnes  _et al._ (2019);
Inomata and Nakama (2019); Clesse  _et al._ (2018); Cai  _et al._ (2019b);
Inomata  _et al._ (2019a, b); Cai  _et al._ (2019c); Yuan  _et al._ (2019a);
Cai  _et al._ (2019d); Lu  _et al._ (2019); Yuan  _et al._ (2019b); Tomikawa
and Kobayashi (2019); De Luca  _et al._ (2019a); Yuan  _et al._ (2019c);
Inomata  _et al._ (2020a, b, c); Yuan and Huang (2020); Papanikolaou  _et al._
(2020); Zhang  _et al._ (2020a); Kapadia  _et al._ (2020a); Zhang  _et al._
(2020b); Domènech  _et al._ (2020a); Dalianis and Kouvaris (2020); Atal and
Domènech (2021). The SIGWs from PBHs are first calculated by Saito and
Yokoyama in Saito and Yokoyama (2009) where they evaluate the energy density
of SIGWs from monochromatic PBHs. They found that SIGWs from the current PBH
DM in our Universe could be detected by pulsar timing arrays and space-based
GW detectors. After the detection of GWs, intriguing studies emerged in this
field in the recent years and there are hundreds of studies concerning SIGWs
from PBHs so far.

Since there are a lot of reviews on PBHs in literature, e.g. some recent
reviews given in Carr  _et al._ (2020); Green and Kavanagh (2020); Carr and
Kuhnel (2020), we mainly focus on SIGWs inevitably generated during the
formation of PBHs in this paper. This paper will be organized as follows. In
Sec. II, we give a brief introduction to the formation of PBHs. The physics
about the SIGWs will be reviewed in Sec. III and then we discuss how to use
the SIGWs to probe PBHs in Sec. IV. Finally, summary and outlook are given in
Sec. V.

##  2 Formation of Primordial Black Holes

In this section, we will introduce the formation of PBHs and take a brief
review on calculating the mass function of PBHs, $\beta$. PBHs are generated
from the collapse of all the matter inside the Hubble volume. Therefore, there
exists a one-to-one correspondence between the mass of PBHs and the comoving
frequency $f_{*}$, namely Carr and Hawking (1974); Carr (1975)

${m_{\mathrm{pbh}}^{*}}\approx 2.3\times
10^{18}M_{\odot}\left(\frac{3.91}{g_{*}^{\mathrm{form}}}\right)^{1/6}\left(\frac{H_{0}}{f_{*}}\right)^{2},$
(1)

where $g_{*}^{\mathrm{form}}$ is the corresponding degrees of freedom and
$H_{0}$ is the Hubble constant by today. Moreover, the one-to-one
correspondence can be transferred to another useful form, namely

${m_{\mathrm{pbh}}^{*}}\approx 2\times
10^{5}M_{\odot}\left({t\over{1s}}\right).$ (2)

The fraction of PBHs in all the dark matter,
$f_{\mathrm{pbh}}\equiv\Omega_{\mathrm{PBH}}/\Omega_{\mathrm{DM}}$ can be
estimated by Nakama  _et al._ (2017)

$\displaystyle f_{\mathrm{pbh}}\simeq 2.5\times
10^{8}\beta\left(\frac{g_{*}^{\mathrm{form}}}{10.75}\right)^{-\frac{1}{4}}\left(\frac{m_{\mathrm{pbh}}}{M_{\odot}}\right)^{-\frac{1}{2}}.$
(3)

On comoving slices, the relation between the primordial comoving curvature
perturbation, $\zeta(k)$, and the density contrast, $\Delta(k)$, at linear
order is given by

$\Delta(k)=\frac{2(1+w)}{5+3w}\left({k\over aH}\right)^{2}\zeta(k),$ (4)

where $w$ is the equation of state and $H$ is the Hubble parameter. The
comoving curvature perturbation $\zeta$ is related to the metric perturbations
by

$\zeta\equiv\psi-\mathcal{H}(v+B),$ (5)

where our notations for the scalar perturbations $\psi$, $v$ and $B$ are
introduced in Eq. (3) and Eq. (3) below. Moreover, the comoving curvature
perturbation is related to the Bardeen potential (see Eq. (23) below) by

$\displaystyle\zeta=\Psi-\frac{2}{3(1+w)}(\mathcal{H}^{-1}\Psi^{\prime}+\Phi)$
(6)

For adiabatic perturbations, $\zeta$ stays constant on superhorizon scales.
Then we have $\Phi=-3(1+w)/(5+3w)\zeta$ by assuming the absence of
anisotropies and $\Phi=-2/3\zeta$ during RD. The density contrast smoothed
over a scale, $R$, is calculated as

$\Delta(\bm{x},R)=\int\mathrm{d}^{3}x^{\prime}W(|\bm{x}-\bm{x^{\prime}}|,R)\Delta(\bm{x^{\prime}}),\qquad\Delta(\bm{k},R)=W(k,R)\Delta(k),$
(7)

with $W$ to be the window function chosen to smooth the density contrast. The
variance of $\Delta(\bm{k},R)$ is given by

$\left\langle\Delta^{2}\right\rangle=\int_{0}^{\infty}\frac{\mathrm{d}k}{k}W^{2}(k,R){4(1+w)^{2}\over(5+3w)^{2}}(kR)^{4}P_{\zeta}(k),$
(8)

where we define the dimensionless curvature power spectrum as
$\left\langle\zeta(k)\zeta(k^{\prime})\right\rangle\equiv\frac{2\pi^{2}}{k^{3}}\delta(\bm{k}+\bm{k^{\prime}})P_{\zeta}(k)$.
Usually, the amplitude of perturbations is assumed to obey Gaussian
distribution. Once the perturbation satisfies the formation criterion, then
PBHs are generated. Therefore, the formation of PBHs can be regarded as the
statistics of peaks of a three-dimensional Gaussian random field, also known
as peak theory Bardeen  _et al._ (1986). PBHs might be generated by very large
perturbations with amplitudes $\mathcal{O}(0.01-0.1)$ if PBHs constitute most
of the DM in our Universe. As a result, the primordial power spectrum that
generate PBHs will have a “bump” on scales much smaller than the CMB scales
Ivanov  _et al._ (1994); Garcia-Bellido  _et al._ (1996); Ivanov (1998);
Yokoyama (1997); Kawasaki  _et al._ (2006); Hertzberg and Yamada (2018);
Inomata  _et al._ (2018, 2017b); Kohri and Terada (2018b). As a benchmark
example, we will assume that the power spectrum is “spicky” at a particular
scale, $k_{*}$, to calculate $\beta$. In this case, the production of PBHs is
almost monochromatic. Otherwise, for a broad power spectrum, one has take into
account the so-called “cloud-in-cloud” problem where a single PBH is swallowed
by the formation of a bigger PBH. Moreover, there are perturbations with
different frequency and different amplitudes which may blur the calculation.

For spicky spectrum, using peak theory, the number density of peaks can be
approximated by (see e.g., Eq. (4.14) in Bardeen  _et al._ (1986))

$n_{\mathrm{pk}}(\nu_{c})\simeq\frac{\left(\left\langle
k^{2}\right\rangle/3\right)^{3/2}}{(2\pi)^{2}}(\nu_{c}^{2}-1)\mathrm{e}^{-\nu_{c}^{2}/2},$
(9)

where the dimensionless threshold is defined as
$\nu_{c}\equiv\Delta_{c}/\sqrt{\left\langle\Delta^{2}\right\rangle}$.
$\left\langle k^{2}\right\rangle$ is defined as

$\left\langle
k^{2}\right\rangle=\frac{1}{\Delta^{2}}\int_{0}^{\infty}\frac{\mathrm{d}k}{k}k^{2}W^{2}(k,R)P_{\zeta}(k)$
(10)

Moreover, $\beta$ is related to $n_{\mathrm{pk}}$ by
$\beta=n_{\mathrm{pk}}(\nu_{c})(2\pi)^{3/2}R^{3}$. Another commonly used
method to calculate $\beta$ is the Press-Schechter formalism, where $\beta$ is
evaluated by simply integrating the probability density function (PDF) beyond
the threshold value

$\beta=\int_{\nu_{c}}^{+\infty}\frac{\mathrm{d}\nu}{\sqrt{2\pi}}e^{-\nu^{2}/2}={1\over
2}\mathrm{erfc}\left({\nu_{c}\over\sqrt{2}}\right).$ (11)

The Press-Schechter formalism only considers the amplitude of the perturbation
but neglecting the higher derivatives. The comparison between Press-Schechter
formalism and peak theory can be found in Young  _et al._ (2014) and the
result showed close agreement differing by a factor of $10$ for large
$\nu_{c}$.

The formation of PBHs depend on many aspects and there are some potential
problems in calculating the mass function of PBHs. We will discuss the
relevant aspects next.

###  2.1 Primordial non-Gaussianities

Formation of PBHs takes place at the tail of the PDF of the perturbations.
Hence any non-Gaussianities that alter the PDF could significantly change the
formation probability of PBHs. The impacts of non-Gaussianities on the
formation of PBHs have been discussed long ago Bullock and Primack (1997);
Ivanov (1998); Pina Avelino (2005); Hidalgo (2007); Klimai and Bugaev (2012).
Besides considering a certain inflation model, a commonly used non-Gaussian
model is the local-type non-Gaussianities, where the perturbation is expanded
by the Gaussian part such that (up to cubic order)

$\zeta=f(\zeta_{g})=\zeta_{g}+F_{\mathrm{NL}}\left(\zeta_{g}^{2}-\left\langle\zeta_{g}^{2}\right\rangle\right)+G_{\mathrm{NL}}\zeta_{g}^{3}.$
(12)

Here $\zeta_{g}$ is the Gaussian part whose PDF is Gaussian. The non-Gaussian
parameter $F_{\mathrm{NL}}$ would skew the PDF and $G_{\mathrm{NL}}$ will
change the kurtosis of the PDF. Since the Press-Schechter formalism simply
integrating the PDF above the threshold value, it is convenient in estimating
the mass function $\beta$ in the presence of non-Gaussianities. Let
${P_{\rm{NG}}}$ to be the PDF of $\zeta$ and ${P_{\rm{G}}}$ to be the PDF of
the Gaussian part, $\zeta_{g}$. Then ${P_{\rm{NG}}}$ can be obtained by
changing the variables, $\zeta_{g}=f^{-1}(\zeta)$, such that

$\displaystyle\beta=\int_{\zeta_{c}}^{\infty}{P_{\rm{NG}}}(\zeta)\mathrm{d}\zeta=\int_{f_{i}^{-1}(\zeta)>\zeta_{c}}\sum_{i=1}^{n}{\mathrm{d}f_{i}^{-1}(\zeta)\over\mathrm{d}\zeta}P_{\mathrm{G}}[f_{i}^{-1}(\zeta)]\mathrm{d}\zeta.$
(13)

This is equivalent to integrate the Gaussian PDF in the region,
$f_{i}^{-1}(\zeta)>\zeta_{c}$. Here, the lower index $i$ indicates the $i$-th
solution of the total $n$ real solutions. In Byrnes  _et al._ (2012), Byrnes
et al. used the above method to calculated $\beta$ within local-type non-
Gaussianities up to cubic order. After that, the formation of PBHs within non-
Gaussianities was studied up to fifth-order Young and Byrnes (2013) and their
results showed a highly sensitive relation between $\beta$ and non-Gaussian
parameters ($F_{\mathrm{NL}}$, $G_{\mathrm{NL}}$, etc.).

Except for changing the variables in Eq. (13), Franciolini  _et al._ (2018)
adopted another method, the path-integral formulation, to calculate $\beta$ in
non-Gaussian regions. The authors expressed $\beta$ by the sum of the N-point
correlation function. Readers interested in the path-integral formulation can
refer to their work and the references therein.

However, the path integral formulation might be impractical to calculate
$\beta$. Riccardi et al. argued that the non-Gaussianities would affect
cumulants at any order Riccardi  _et al._ (2021). Therefore, to get the exact
result for $\beta$, one has to sum over all the N-point correlation function,
which is impractical. They also proposed a semi-analytical expression to
estimate the mass function of PBHs, see Eq. (8) in Riccardi  _et al._ (2021).

###  2.2 The non-linear effects between density contrast and curvature
perturbation

The standard procedure to calculate the density contrast is based on Eq. (4).
However, this equation is just the linear relation between the density
contrast and the comoving curvature perturbation. In the comoving slicing, the
non-linear relation in the long-wavelength approximation is given by Harada
_et al._ (2015); Yoo  _et al._ (2018); Musco (2019)

$\displaystyle\Delta(r,t)=-\frac{4(1+w)}{5+3w}\left({1\over
aH}\right)^{2}\mathrm{e}^{-5\zeta(r)/2}\nabla^{2}\mathrm{e}^{\zeta(r)/2}.$
(14)

Due to the linear relation, the PDF of density contrast will no longer obey
Gaussian distribution even if the PDF of curvature perturbations is Gaussian.
This is an unavoidably generated non-Gaussianities that make the PBHs
inevitably form in non-Gaussian regions. Recently, these “intrinsic non-
Gaussianities” generated by the non-linear effects were studied in Yoo  _et
al._ (2018); Kawasaki and Nakatsuka (2019); De Luca  _et al._ (2019b); Young
_et al._ (2019) using either peak theory or Press-Schechter formalism. They
found that the “intrinsic non-Gaussianities would slightly suppress the PBH
formation. More precisely, to produce the same abundance of PBH in our
Universe, the power spectrum of $\zeta$ need to be amplified by a factor of
$\sim\mathcal{O}(2)$ if using the linear relation, Eq. (4), to calculate the
abundance.

###  2.3 Window Function and Power spectrum

From Eq. (8), it is clear that the formation of PBHs depends on the choice of
window functions and the power spectrum of the primordial scalar
perturbations. The power spectrum can be given by the inflation model, while
the choice of the window function is a coarse-graining procedure. Despite
several commonly used window functions in literature, there is no physical
interpretation on which window function should be used. The choice of window
function will lead to uncertainties in calculating the mass function of PBHs,
see Ando  _et al._ (2018a); Young (2019); Tokeshi  _et al._ (2020)

The abundance of PBHs introduced in Eq. (9) and Eq. (11) can be applied to a
narrow power spectrum, where the PBHs are about monochromatic and one can
assume that the PBHs are formed at the same time. However, for a broad power
spectrum, PBHs of different masses are expected to be formed at a different
time, and the formation process will be rather complicated.

For the broad case, the Press-Schechter formalism failed, and one has to adopt
the peak theory to evaluate the results. In Moradinezhad Dizgah  _et al._
(2019); De Luca  _et al._ (2020b), the authors studied the formation of PBHs
for a tilted broad power spectrum, namely

$P_{\zeta}=A(k/k_{s})^{n_{p}}\Theta(k_{s}-k)\Theta(k-k_{l}),$ (15)

where $\Theta$ denotes the Heaviside theta function. For a broad and flat
spectrum, $n_{p}=0$, the corresponding mass distribution of PBHs is dominated
by a single $m_{\rm{pbh}}$. They also found that the mass function has a
power-law tail, scaling as $m_{\rm{pbh}}^{-3/2}$.

Recently, Yoo et al. Yoo  _et al._ (2021) proposed a modified procedure to
calculate PBH abundance in which they consider the intrinsic non-Gaussian
effects within an arbitrary power spectrum in the absence of primordial non-
Gaussianities. The results for a narrow spectrum showed no window function
dependence, while the results for a broad power spectrum depends largely on
the choice of window functions. They concluded that the top-hat window
function in Fourier space would be the best choice since it minimizes the
required property in theoretical PBH estimation.

###  2.4 The formation criterion

Studies on the critical value of PBH formation can be traced back since 1975
when Carr first estimated the value of $\Delta_{c}$ by considering simplified
Jeans length in Newtonian gravity, $\Delta_{c}=c_{s}^{2}$ Carr (1975), where
$c_{s}=1/\sqrt{3}$ is the sound speed during RD. However, the exact value of
$\Delta_{c}$ is still a debating question since it depends on the initial
density profile Polnarev and Musco (2007), the equation of state and the sound
speed Jedamzik (1997); Byrnes  _et al._ (2018), primordial non-Gaussianities
Kehagias  _et al._ (2019), the primordial scalar power spectrum Germani and
Musco (2019) and even the window function Young (2019).

Despite some spherically numeric simulations Jedamzik and Niemeyer (1999);
Shibata and Sasaki (1999); Musco  _et al._ (2005); Hawke and Stewart (2002);
Polnarev and Musco (2007); Musco  _et al._ (2009); Musco and Miller (2013);
Nakama  _et al._ (2014); Harada and Jhingan (2016); Harada  _et al._ (2015);
Musco (2019); Escrivà  _et al._ (2020a) on the formation of PBHs, Harada et
al. proposed an analytical expression,
$\Delta_{c}=[3(1+w)/(5+3w)]\sin^{2}[\pi\sqrt{w}/(1+3w)]\simeq 0.41$ during RD
Harada  _et al._ (2013). Escriva et al. also proposed an approximate
expression for $\Delta_{c}$ for $w\in[1/3,1]$ Escrivà  _et al._ (2020b). The
analytical $\Delta_{c}$ in a general cosmological background is studied by
Escrivà  _et al._ (2020b) where they consider the equation of state to be
$w\in(0,1]$.

Recently, Kehagias  _et al._ (2019) estimated the effects of “intrinsic non-
Gaussianities” generated by the non linear relation, Eq. (14), on
$\Delta_{c}$. They found that the relative change of $\Delta_{c}$ is at the
percent level, which might not be significant due to other uncertainties in
estimating the formation of PBHs. More recently, a detailed study based on
numeric simulation was made by Musco et al. Musco  _et al._ (2020) where they
consider the intrinsic non-Gaussian effects and some commonly used power
spectrum that generate the PBHs. They also gave an analytical approach to
estimate $\Delta_{c}$ for all possible shapes of the power spectrum, see Eq.
(19) in Musco  _et al._ (2020).

Apart from using $\Delta_{c}$, Shibata and Sasaki proposed that the compaction
function, which equals to half of the volume average of the density per
turbations in the long-wavelength limit, can be used to described the
formation of PBHs Shibata and Sasaki (1999). A general definition for the
compaction function was introduced in Harada  _et al._ (2015) such that
$\mathcal{C}=2\delta M/R$. Here $R$ represents the areal radius and $\delta M$
stands for the mass difference between the Misner-Sharp mass insides a sphere
of radius $r$ and the mass inside a sphere with the same radius in the FLRW
universe.

Using the maximum value of the compaction function, $\mathcal{C}_{\max}$, as a
formation criterion, there would be less dependence on the shape of the
density perturbation Shibata and Sasaki (1999) and this is later confirmed
numerically by Harada and Jhingan (2016). Harada and Jhingan (2016) also found
that there would be less dependence on the lapse function if using
$\mathcal{C}_{\max}$ instead of $\Delta_{c}$. Recently, it was shown in
Escrivà  _et al._ (2020a) that the threshold value of compaction function,
$\mathcal{C}_{\max,c}$, is only sensitive to the curvature at the maximum and
$\mathcal{C}_{\max,c}$ is, to some approximations, a universal quantity to
describe the PBH formation.

##  3 Scalar-induced gravitational waves

In this section, we will introduce the calculations of SIGWs during RD. Let’s
begin from the most generic perturbed metric, which contains scalar modes and
tensors modes, namely

$\displaystyle g_{00}$ $\displaystyle=$ $\displaystyle-1-2\phi,$
$\displaystyle g_{0i}$ $\displaystyle=$ $\displaystyle a\partial_{i}B,$
$\displaystyle g_{ij}$ $\displaystyle=$ $\displaystyle
a^{2}\delta_{ij}+a^{2}\left({1\over
2}h_{ij}-2\delta_{ij}\psi+2\partial_{i}\partial_{j}E\right),$ (16)

where $\phi$, $B$, $\psi$ and $E$ are linear scalar perturbations while the
vector perturbations are not considered here. For vector-induced GWs and
tensor-induced GWs, we refer the readers to Gong (2019). Here, $h_{ij}$
denotes the second-order tensor mode. We do not consider the first-order
tensor mode, since the second-order effects will far exceed the linear order
effects during the formation of PBHs Saito and Yokoyama (2009). During RD, the
cosmological background is described by perturbed perfect fluid, which reads

$\displaystyle T_{00}$ $\displaystyle=$
$\displaystyle\rho+2\rho\phi+\delta\rho,$ $\displaystyle T_{0i}$
$\displaystyle=$
$\displaystyle-\rho\partial_{i}B-P\partial_{i}v-\rho\partial_{i}v,$
$\displaystyle T_{ij}$ $\displaystyle=$ $\displaystyle(P+\delta
P)\delta_{ij}-2P\psi\delta_{ij}+2P\partial_{i}\partial_{j}E,$ (17)

where $v$ is the velocity potential, $P$ and $\rho$ are the pressure and the
energy density of the background while $\delta P$ and $\delta\rho$ are the
corresponding first-order perturbations. Here we neglect the anisotropies
caused by neutrinos and photon Saga  _et al._ (2015).

As we shall see, besides the GWs, there is only one physical degree of freedom
in our calculation. Consider a first-order change in the coordinate such that

$\tilde{\eta}=\eta+T,~{}\tilde{x}^{i}=x^{i}+\partial^{i}L,$ (18)

and the scalar modes will transform as Malik and Wands (1998)

$\displaystyle\tilde{\phi}$ $\displaystyle=$
$\displaystyle\phi+\mathcal{H}T+T^{\prime},$ (19) $\displaystyle\tilde{\psi}$
$\displaystyle=$ $\displaystyle\psi-\mathcal{H}T,$ (20)
$\displaystyle\tilde{B}$ $\displaystyle=$ $\displaystyle B-T+L^{\prime},$ (21)
$\displaystyle\tilde{E}$ $\displaystyle=$ $\displaystyle E+L,$ (22)

where a prime denotes a derivative with respect to the conformal time $\eta$.
The most generic perturbed metric has four scalar modes. However, one can fix
$T$ to determine the time-slicing and choose $L$ to decide the spatial
coordinates on the hyper-surfaces. This would reduce two degrees of freedom.
Furthermore, the Einstein equation will reduce one more degree of freedom and
give the equation of motion for the last degree of freedom.

The simplest gauge to calculate the SIGWs is the Newton gauge, also known
widely as orthogonal zero-shear gauge, longitudinal gauge. Newton gauge
demands that $\tilde{B}=\tilde{E}=0$. This can be done by choosing $L=-E$ and
$T=B-E^{\prime}$. Then, the remaining two modes are just the Bardeen potential

$\displaystyle\Phi\equiv\phi+\mathcal{H}\sigma+\sigma^{\prime},\quad\Psi\equiv\psi-\mathcal{H}\sigma,$
(23)

where $\sigma\equiv E^{\prime}-B$ is the shear potential. The advantage of the
Newton gauge is that the degrees of freedom are completely fixed, and there
are no more gauge modes. The only modes left are just the Bardeen potential
$\phi=\Phi$ and $\psi=\Psi$. Write down the first-order Einstein equation in
terms of the Bardeen potential, one finds that $\Phi=\Psi$ and the equation of
motion for the only scalar mode left is

$\Phi^{\prime\prime}+4\mathcal{H}\Phi-{1\over 3}\nabla^{2}\Phi=0.$ (24)

Keeping the decay modes of Eq. (24), the solution of $\Phi$ in Fourier space
is given by

$\displaystyle\Phi(\bm{k},\eta)\equiv\Phi_{k}T_{\Phi}(k\eta)=\Phi_{k}{9\over
x^{2}}\left({\sin(x/\sqrt{3})\over x/\sqrt{3}}-\cos(x/\sqrt{3})\right),$ (25)

where we label the primordial value as $\Phi_{k}$ and $T_{\Phi}$ is the
normalized transfer function such that $T_{\Phi}(0)=1$. We also introduce the
dimensionless variable $x\equiv k\eta$. The value of $\Phi_{k}$ is given by
certain inflation models and we will treat it as free parameters in our
discussion. Moreover, the first-order GWs is given by

$h_{ij}^{\prime\prime(1)}+2\mathcal{H}h_{ij}^{\prime(1)}-\nabla^{2}h_{ij}^{(1)}=0,$
(26)

which is a source-free equation. The linear GWs will decay from the primordial
value as $\eta^{-1}$ during RD. During the formation of PBHs, the scalar power
spectrum are enhanced to $\mathcal{O}(0.01)$ if PBH represents the main part
of the DM. In this case, the second-order SIGWs would be stronger than the
linear GWs, and hence we will not consider $h_{ij}^{(1)}$. Readers interested
in second-order GWs induced by $h_{ij}^{(1)}$ can refer to Gong (2019).

Up to second-order, the evolution of the GWs can be written as

$h_{ij}^{\prime\prime}+2\mathcal{H}h_{ij}^{\prime}-\nabla^{2}h_{ij}=-4\mathcal{T}_{ij}^{\ell
m}S_{\ell m},$ (27)

where $\mathcal{T}_{ij}^{\ell
m}=e_{ij}^{(+)}(\bm{k})e^{(+)lm}(\bm{k})+e_{ij}^{(\times)}(\bm{k})e^{(\times)lm}(\bm{k})$
is the transverse and traceless projection operator. The polarization tensors
are defined as $(e_{i}e_{j}-\bar{e}_{i}\bar{e}_{j})/\sqrt{2}$ and
$(e_{i}\bar{e}_{j}+\bar{e}_{i}e_{j})/\sqrt{2}$ for $+$ and $\times$ modes.
Here we choose $e=(1,0,0)$, $\bar{e}=(0,1,0)$ and $\bm{k}=(0,0,k)$.

The source term in eq. (27) reads

$\displaystyle
S_{ij}=3\Phi\partial_{i}\partial_{j}\Phi-{2\over\mathcal{H}}\partial_{i}\Phi^{\prime}\partial_{j}\Phi-{1\over\mathcal{H}^{2}}\partial_{i}\Phi^{\prime}\partial_{j}\Phi^{\prime},$
(28)

where we have expressed the density perturbation $\delta\rho$, the velocity
potential $v$ and the pressure perturbation $\delta P$ in terms of the scalar
modes $\Phi$. Overall, the perfect fluid acts like an intermediary, and the
GWs seem to be generated by the scalar modes. This is why SIGWs gain its name.
Eq. (27) can be solved by Green’s function in Fourier space, namely

$h(\eta,\bm{k})={1\over
a(\eta)}\int_{0}^{\eta}g_{k}(\eta;\eta^{\prime})a(\eta^{\prime})S(\eta^{\prime},\bm{k})\mathrm{d}\eta^{\prime}.$
(29)

We define $S(\eta,\bm{k})\equiv-4e^{ij}(\bm{k})\mathcal{S}(\eta,k)$ and
$\mathcal{S}(\eta,k)$ is the source term transferred to Fourier space. Here,
we define the Fourier transform of $h_{ij}$ to be

$h_{ij}(\eta,\bm{x})=\int{\mathrm{d}^{3}k\over(2\pi)^{3}}\left[h^{(+)}(\eta,\bm{k})e_{ij}^{(+)}(\bm{k})+h^{(\times)}(\eta,\bm{k})e_{ij}^{(\times)}(\bm{k})\right]e^{i\bm{k}\cdot\bm{x}}.$
(30)

In the following part, we write $h(\eta,\bm{k})$ to denote either the plus
mode or the cross mode. In our convention, $S(\eta,\bm{k})$ takes the form

$\displaystyle
S(\eta,\bm{k})=-4\int\frac{\mathrm{d}^{3}p}{(2\pi)^{3/2}}\Big{(}e^{ij}p_{i}p_{j}\Big{)}\Phi_{p}\Phi_{|\bm{p}-\bm{k}|}F(|\bm{p}|,|\bm{k}-\bm{p}|,\eta).$
(31)

The transfer function is defined as

$\displaystyle{F}(u,v,x)=3T_{\phi}(ux)T_{\phi}(vx)+uxT_{\phi}^{\prime}(ux)T_{\phi}(vx)+vxT_{\phi}^{\prime}(vx)T_{\phi}(vx)+uvx^{2}T_{\phi}^{\prime}(ux)T_{\phi}^{\prime}(vx).$
(32)

For convenience, we introduce the dimensionless variable $u\equiv p/k$,
$v\equiv|\bm{p}-\bm{k}|/k$ and $x\equiv k\eta$. Unless otherwise being stated,
the prime with $T^{\prime}(y)$ denotes the derivative with respect to $y$,
other than the conformal time. The Green’s function in Eq. (29) takes the form
$g_{k}(\eta;\eta^{\prime})={1\over k}\sin(k\eta-k\eta^{\prime})$ during RD. A
more important quantity in observation is the density parameter of the
stochastic GW background defined as the energy of GWs per logarithm frequency
normalized by the critical energy $\rho_{c}(\eta)$

$\Omega_{\mathrm{GW}}(k,\eta)\equiv{1\over\rho_{c}}{\mathrm{d}\rho_{\mathrm{GW}}\over\mathrm{d}\ln
k}={k^{3}\over
48\pi^{2}}\left(k\over\mathcal{H}\right)^{2}\overline{\left\langle|h(\eta,{\bm{k}})|^{2}\right\rangle},$
(33)

where the overline denotes the oscillating average. To calculate the
stochastic GW background of the SIGWs, we have to know the primordial scalar
power spectrum, which is given by specific inflation models. In this paper, we
do not consider certain inflation model, instead, we will parameterize the
power spectrum by some commonly used function in literature. The density
parameter is evaluated as (see e.g. Kohri and Terada (2018a))

$\displaystyle\Omega_{\mathrm{GW}}(k,\eta)\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!$
$\displaystyle=\left(\frac{k}{\mathcal{H}}\right)^{2}{k^{3}\over
48\pi^{2}{a(\eta)}^{2}}\int\mathrm{d}\tilde{\eta_{1}}\mathrm{d}\tilde{\eta_{2}}{{g_{k}(\eta;\tilde{\eta_{1}})g_{k^{\prime}}(\eta;\tilde{\eta_{2}})}}a({\tilde{\eta_{1}}})a(\tilde{\eta_{2}}){\left\langle{S}(\tilde{\eta_{1}},\bm{k}){S}(\tilde{\eta_{2}},\bm{k^{\prime}})\right\rangle}.$
(34) $\displaystyle=\left(\frac{k}{\mathcal{H}}\right)^{2}{4\pi^{2}k^{3}\over
3a(\eta)^{2}}\int\frac{\mathrm{d}^{3}p}{(2\pi)^{3}}\int\mathrm{d}\tilde{\eta_{1}}\mathrm{d}\tilde{\eta_{2}}{a(\tilde{\eta_{1}})a(\tilde{\eta_{2}})}{g_{k}(\eta;\tilde{\eta_{1}})g_{k}(\eta;\tilde{\eta_{2}})}\Big{(}e^{ij}p_{i}p_{j}\Big{)}^{2}{1\over
p^{3}|\bm{k}-\bm{p}|^{3}}$ $\displaystyle\times
P_{\Phi}(k)P_{\Phi}(|\bm{k}-\bm{p}|)\Big{[}{F}(|\bm{p}|,|\bm{k}-\bm{p}|,\tilde{\eta_{1}}){F}(|\bm{p}|,|\bm{k}-\bm{p}|,\tilde{\eta_{2}})+{F}(|\bm{p}|,|\bm{k}-\bm{p}|,\tilde{\eta_{1}}){F}(|\bm{k}-\bm{p}|,|\bm{p}|,\tilde{\eta_{2}})\Big{]}$
$\displaystyle=\left(\frac{k}{\mathcal{H}}\right)^{2}{8\pi^{2}k\over
3}\int\frac{\mathrm{d}^{3}p}{(2\pi)^{3}}\left(\int_{0}^{\eta}\mathrm{d}\tilde{\eta_{1}}{a(\tilde{\eta_{1}})\over
a(\eta)}kg_{k}(\eta;\tilde{\eta_{1}})\tilde{F}(|\bm{p}|,|\bm{k}-\bm{p}|,\tilde{\eta_{1}})\right)^{2}\Big{(}e^{ij}p_{i}p_{j}\Big{)}^{2}{1\over
p^{3}|\bm{k}-\bm{p}|^{3}}$ $\displaystyle\qquad\times
P_{\Phi}(k)P_{\Phi}(|\bm{k}-\bm{p}|)$ $\displaystyle\simeq{1\over
6}\int_{0}^{\infty}\mathrm{d}u\int_{|1-u|}^{1+u}\mathrm{d}v~{}{v^{2}\over
u^{2}}\Big{[}1-\left({1+v^{2}-u^{2}\over
2v}\right)^{2}\Big{]}^{2}P_{\Phi}(uk)P_{\Phi}(vk)\overline{I^{2}(u,v,x)},$

where in the last step of Eq. (34), we define the power spectrum as

$\displaystyle\left\langle\Phi(\bm{k})\Phi(\bm{k^{\prime}})\right\rangle\equiv{2\pi^{2}\over
k^{3}}P_{\Phi}(k)\delta(\bm{k}+\bm{k^{\prime}}),$ (35)

and the kernel function to be

$\displaystyle
I(u,v,x)\equiv\int_{0}^{x}\mathrm{d}\tilde{x}~{}\tilde{x}\sin(x-\tilde{x})\tilde{F}(u,v,x),$
(36)

where $\tilde{F}(u,v,x)\equiv(F(u,v,x)+F(v,u,x))/2$ is the symmetric part of
the transfer function. We also used that $\mathcal{H}=\eta^{-1}$ and
$a(\tilde{\eta})/a(\eta)=\tilde{\eta}/\eta$ during RD. We sum over the two
polarization modes in Eq. (34) and the $e^{ij}$ in the expression refers to
either the “$+$” mode or the “$\times$” mode. In the last step of Eq. (34),
the overline denotes the oscillating average, namely $\sin^{2}x\to{1\over 2}$,
$\cos^{2}x\to{1\over 2}$ and $\sin x\cos x\to 0$. The $\Omega_{\mathrm{GW}}$
has a time dependence through $I(u,v,x)$ and, as we will demonstrate below,
the function $I(u,v,x)$ will converge to a finite value at late time.

The density parameter given by Eq. (34) is independent of the position,
rendering the SIGWs are isotropic GWs. This is due to the statistical
homogeneity of the FLRW metric. However, the GW detectors would detect GWs
from different directions, resulting in angular anisotropies of the SIGW
signal. This was studied in Bartolo  _et al._ (2020). In the absence of
primordial non-Gaussianities, the authors found that the anisotropies are
negligible by today due to the propagation effects.

After the horizon entry, the energy of GWs redshift as radiation,
$\rho_{\rm{GW}}\propto a^{-4}$. Therefore, the current GW density parameter,
$\Omega_{\mathrm{GW,0}}$ would be

$\Omega_{\mathrm{GW,0}}(\eta_{0},f)=\Omega_{r}\Omega_{\mathrm{GW}}(\eta_{c},f),$
(37)

where $\Omega_{r}$ is the energy density fraction of radiation by today and we
neglect the effects of relativistic degrees of freedom. Here,
$\Omega_{\mathrm{GW}}(\eta_{c},f)$ is the late time value when
$\Omega_{\mathrm{GW}}(\eta,f)$ becomes a constant and can be computed by
taking $I(u,v,x\to\infty)$ (Here, a constant means that
$\Omega_{\mathrm{GW}}(\eta,f)$ is independent of time). This result would be
more clear if one considers the evolution of tensor modes during both MD and
RD. See e.g., Kohri and Terada (2018a) where the author studied analytically
the SIGWs during MD, RD and RD-to-MD transition. The analytical expression for
$\overline{I^{2}(u,v,x)}$ was derived in Espinosa  _et al._ (2018); Kohri and
Terada (2018a). The indefinite integral of Eq. (36) is given by

$\displaystyle I(u,v,x)\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!$
$\displaystyle=-{27(u^{2}+v^{2}-3)^{2}\over
16u^{3}v^{3}}\Bigg{\\{}\Bigg{[}\mathrm{Si}\left(\left(1-{(u+v)\over\sqrt{3}}\right)x\right)+\mathrm{Si}\left(\left(1+{(u+v)\over\sqrt{3}}\right)x\right)-\mathrm{Si}\left(\left(1+{(u-v)\over\sqrt{3}}\right)x\right)$
(38)
$\displaystyle-\mathrm{Si}\left(\left(1-{(u-v)\over\sqrt{3}}\right)x\right)\Bigg{]}\cos
x+\Bigg{[}\mathrm{Ci}\left(\left(1+{(u-v)\over\sqrt{3}}\right)x\right)+\mathrm{Ci}\left(\left(1-{(u-v)\over\sqrt{3}}\right)x\right)$
$\displaystyle-\mathrm{Ci}\left(\Bigg{|}1-{(u+v)\over\sqrt{3}}\Bigg{|}x\right)-\mathrm{Ci}\left(\left(1+{(u+v)\over\sqrt{3}}\right)x\right)\Bigg{]}\sin
x+{1\over u^{2}+v^{2}-3}$
$\displaystyle\times\left(4uv+(u^{2}+v^{2}-3)\ln|1-\frac{4uv}{(u+v)^{2}-3}\Bigg{|}\right)\sin
x+{12\over x^{2}(u^{2}+v^{2}-3)^{2}}$
$\displaystyle\times\Bigg{[}2u\cos{ux\over\sqrt{3}}\left(vx\cos{vx\over\sqrt{3}}-\sqrt{3}\sin{vx\over\sqrt{3}}\right)+2u\left(-vx+\sqrt{3}\sin{vx\over\sqrt{3}}\right)$
$\displaystyle+\sin{ux\over\sqrt{3}}\left(2\sqrt{3}v-2\sqrt{3}v\cos{vx\over\sqrt{3}}+(u^{2}+v^{2}-3)x\sin{vx\over\sqrt{3}}\right)\Bigg{]}\Bigg{\\}}.$

Then the late time expression, $\overline{I^{2}(u,v,x\to\infty)}$, takes the
form

$\displaystyle\overline{I^{2}}=\frac{729(u^{2}+v^{2}-3)^{2}}{512u^{6}v^{6}}\Bigg{\\{}\Big{(}-4uv+(u^{2}+v^{2}-3)\ln|{3-(u+v)^{2}\over
3-(u-v)^{2}}\Big{|}\Big{)}^{2}+\pi^{2}\left(u^{2}+v^{2}-3\right)^{2}\Theta(u+v-\sqrt{3})\Bigg{\\}}.$
(39)

Below we use the comoving curvature perturbation $\zeta$, which has the
relation $\zeta=(3/2)\Phi$ to calculate the SIGWs. A widely used model for the
power spectrum is the infinite narrow spectrum

$P_{\zeta}(k)=Ak_{*}\delta(k-k_{*}),$ (40)

with a dimensionless amplitude $A$ peaked at $k_{*}$. This power spectrum
corresponds to a monochromatic PBH formation and is convenient to study the
properties of SIGWs from PBHs. For this spectrum, the density parameter has an
analytical form, namely

$\Omega_{\mathrm{GW}}(k)={3A^{2}\Omega_{\rm{r}}\over
64}\tilde{k}^{2}\left(1-{\tilde{k}^{2}\over
4}\right)^{2}(2-3\tilde{k}^{2})^{2}\Theta(2-\tilde{k})\Bigg{(}(2-3\tilde{k}^{2})^{2}\pi^{2}\Theta(2-\sqrt{3}\tilde{k})+\left(-4+(2-3\tilde{k}^{2})\ln|1-{4\over
3\tilde{k}^{2}}\Big{|}\right)^{2}\Bigg{)},$ (41)

where we define $\tilde{k}\equiv k/k_{*}$. Another example is a log-normal
power spectrum parameterized by

$P_{\zeta}(k)={A\over\sqrt{2\pi\sigma_{*}^{2}}}\exp\left(-\frac{\ln\tilde{k}^{2}}{2\sigma_{*}^{2}}\right),$
(42)

where $\sigma_{*}$ is a dimensionless parameter which denotes the width of the
power spectrum. The SIGWs generated by this power spectrum has been carefully
studied in Pi and Sasaki (2020).

Finally, we consider a box spectrum described by

$P_{\zeta}(k)=\frac{A}{\ln\left({k_{\mathrm{max}}\over
k_{\mathrm{min}}}\right)}\Theta(\tilde{k}-k_{\mathrm{min}})\Theta(k_{\mathrm{max}}-\tilde{k}).$
(43)

In the above case, we have normalized the power spectrum so that $\int
P_{\zeta}(k)~{}\mathrm{d}\ln k=A=\left\langle\zeta^{2}\right\rangle$ is the
variance of the perturbations. The $\Omega_{\mathrm{GW}}(k)$ for these three
power spectrum is given in Fig. 1.

As shown in the figure, on small scales, the SIGWs will have a cut-off
wavelength. The cut-off is due to the momentum conservation such that two
modes $\bf{k_{1}}$ and $\bf{k_{2}}$ generate a mode with
$\bf{k}=\bf{k_{1}}+\bf{k_{2}}$. Therefore, for a power spectrum defined on
$[k_{\min},k_{\max}]$, the cut-off wavelength of the corresponding SIGWs will
be $2k_{\max}$. Since a constant gravitational potential will not change the
distribution of matter, it will not induce GWs. After the perturbation re-
enters the horizon at $k\simeq k_{*}$, they will induce GWs, and the shape of
$\Omega_{\mathrm{GW}}$ is dominated by the power spectrum. For instance, a
power spectrum with a peak located at $k_{*}$ will induce GWs with a peak at
$2k_{*}/\sqrt{3}$ where $1/\sqrt{3}$ comes from the sound of speed during RD
while a box spectrum will induce a smooth signal. Finally, on large scales,
the scaling of SIGWs shows a log-dependent behavior,
$\Omega_{\mathrm{GW}}\propto f^{2}\ln\left({4f_{*}^{2}\over 3f^{2}}\right)$ or
$\Omega_{\mathrm{GW}}\propto f^{3}\ln\left({4f_{*}^{2}\over 3f^{2}}\right)$
where $f_{*}$ is a pivot scale of the power spectrum. On large scales, the
power spectrum has decreased to negligible value, and the signal is dominated
by the evolution of scalar perturbations. A strict proof of the log-dependent
scaling is given in Yuan  _et al._ (2019b) for a general spectrum. This
special scaling is very important in distinguishing the signal of SIGWs from
other stochastic GW background, and it could be smoking guns in detecting
SIGWs. In Yuan  _et al._ (2019b), Yuan et al. estimated the distinguishability
of LISA. They consider a fiducial case, $\Omega_{\mathrm{GW}}^{fid}(k)$,
generated by a log-normal power spectrum with $\sigma_{*}=0.5$ and the other
one is described by

$\Omega_{\mathrm{GW}}^{m}(k)=\left\\{\begin{aligned}
&\Omega_{\mathrm{GW}}^{fid}(0.1k_{*}){\left({k\over
0.1k_{*}}\right)^{3}}~{},&&\ \hbox{for}\ k<0.1k_{*},\\\
&\Omega_{\mathrm{GW}}^{fid}(k)~{},&&\ \hbox{for}\ k\geq 0.1k_{*}.\\\
\end{aligned}\right.$ (44)

These two cases have different scalings in the infrared region, as shown in
Fig. 2. A statistic quantity, $\delta\chi^{2}$, which characterize the
discrepancy of two models, takes the form Kuroyanagi  _et al._ (2018)

$\displaystyle\delta\chi^{2}\simeq
T\int_{0}^{\infty}\mathrm{d}f\left(\frac{\Omega_{\mathrm{GW}}^{fid}-\Omega_{\mathrm{GW}}^{m}}{\Omega_{\mathrm{GW}}^{m}+\Omega_{n}}\right)^{2},$
(45)

where $T$ is the observation time and $\Omega_{n}$ is the noise density
parameter of the detector. Scanning the mass of PBHs, the result of
$\delta\chi^{2}$ is Fig. 3. It is shown that LISA can well distinguish the
scaling of two models beyond $5\sigma$ in a wide mass range, especially for
those PBHs of
$m_{\rm{pbh}}\in[10^{-16},10^{-14}]\cup[10^{-13},10^{-12}]M_{\odot}$ which
could still represent all the DM in our Universe. To summarize, the log-
dependent scaling, $\Omega_{\mathrm{GW}}\propto f^{2}\ln\left({4f_{*}^{2}\over
3f^{2}}\right)$ or $\Omega_{\mathrm{GW}}\propto f^{3}\ln\left({4f_{*}^{2}\over
3f^{2}}\right)$ can be smoking guns in detecting SIGWs.

Since there are no observational results concerning the primordial scalar
power spectrum on small scales, the power spectrum is usually parameterized by
some given function in literature. A more realistic way is to consider the
primordial power spectrum from a given inflation model and then calculate the
corresponding SIGWs. See reflevant works Kawasaki  _et al._ (2013); Choudhury
and Mazumdar (2014); Inomata  _et al._ (2017a); Di and Gong (2018); Ando  _et
al._ (2018b, c); Xu  _et al._ (2020); Özsoy and Tasinato (2020); Gao and Yang
(2019); Lin  _et al._ (2020); Ballesteros  _et al._ (2020); Liu  _et al._
(2020); Braglia  _et al._ (2020a); Fu  _et al._ (2020); Dalianis and Kritos
(2021); Yi  _et al._ (2020a); Dalianis (2020); Aldabergenov  _et al._ (2021);
Ragavendra  _et al._ (2020a); Bhaumik and Jain (2020); Zhou  _et al._ (2020);
Ragavendra  _et al._ (2020b); Yi  _et al._ (2020b); Braglia  _et al._ (2020b);
Gao  _et al._ (2020); Gao (2021).

###  3.1 Scalar Induced Gravitational Waves in a General Cosmological
Background

For a constant equation of state and let $c_{s}^{2}=w$, the semi-analytical
solutions of SIGWs were obtained by Domènech in Domènech (2019) for $0<w\leq
1$. Then, the results for constant $w$ were extended to $w<0$ in Domènech  _et
al._ (2020b). For the most general situation where $w$ and $c_{s}^{2}$ may
varied, there are no semi-analytical solutions to SIGWs, and a numeric method
should be adopted. For instance, during the QCD phase transition, $w$ and
$c_{s}^{2}$ may become slightly smaller than $1/3$ and the equation
$c_{s}^{2}=w$ no longer holds. The SIGWs generated during the QCD epoch were
numerically calculated in Abe  _et al._ (2020). The derivation within a
general $w(\eta)$ and $c_{s}(\eta)$ differs from that in RD in the following
aspects.

Firstly, the equation of motion for linear scalar perturbations become

$\Phi^{\prime\prime}+3\mathcal{H}(\eta)(1+c_{s}^{2})\Phi^{\prime}+3\mathcal{H}^{2}(c_{s}^{2}-w)\Phi-
c_{s}^{2}\nabla^{2}\Phi=0,$ (46)

where the conformal Hubble parameter now become
$\mathcal{H}(\eta)={2\over(1+3w)\eta}$. This equation should be numerically
solved with the initial condition to be $\Phi_{\bm{k}}(\eta\to
0)=\Phi_{k}T_{\Phi}(k\eta\to 0)=\Phi_{k}$.

Secondly, the equation of motion for the second-order tensor modes takes the
same form as Eq. (27) but the source term changes to

$\displaystyle S_{ij}=2\Phi\partial_{i}\partial_{j}\Phi-{4\over
3(1+w)}\left(\partial_{i}\Phi+{\partial_{i}\Phi^{\prime}\over\mathcal{H}(\eta)}\right)\left(\partial_{j}\Phi+{\partial_{j}\Phi^{\prime}\over\mathcal{H}(\eta)}\right).$
(47)

The Green’s function and the transfer function methods can still be applied to
solve $h_{ij}$. Since the form of the equation of motion for $h_{ij}$ is
unchanged, the solution to $h_{ij}$ takes the same form as Eq. (29) except for
the Green’s function, which now satisfies the equation

$\displaystyle\left(\partial_{\eta}^{2}+k^{2}-{1-3w(\eta)\over
2}\mathcal{H}^{2}\right)g_{k}(\eta,\tilde{\eta})=\delta(\eta-\tilde{\eta}),$
(48)

The solution of $g_{k}(\eta,\tilde{\eta})$ can be obtained by

$g_{k}(\eta,\tilde{\eta})=\frac{u(\eta)v(\tilde{\eta})-u(\tilde{\eta})v(\eta)}{u^{\prime}(\tilde{\eta})v(\tilde{\eta})-u(\tilde{\eta})v^{\prime}(\tilde{\eta})},$
(49)

with $u(\eta)$ and $v(\eta)$ to be the two independent homogeneous solutions.
Move on to the density parameter, the final expression for
$\Omega_{\mathrm{GW}}$ takes the same form as the last line of Eq. (34) but
the kernel function should now be evaluated as

$I(u,v,x)=k^{2}\int_{0}^{\eta}\mathrm{d}\tilde{\eta}{a(\tilde{\eta})\over
a(\eta)}g_{k}(\eta,\tilde{\eta})\tilde{F}(|\bm{p}|,|\bm{k}-\bm{p}|,\tilde{\eta}).$
(50)

After calculating the oscillating average of $I(u,v,x)$, the density
parameter, $\Omega_{\mathrm{GW}}$, can be obtained. Readers interested in
numerically solutions may refer to Abe  _et al._ (2020) for details where the
authors introduce some technique to reduce the computational complexity. The
$\Omega_{\mathrm{GW}}$ of SIGWs generated during the QCD phase transition is
shown in Fig. 4 in Abe  _et al._ (2020) for a delta-power spectrum.

###  3.2 Higher-Order Corrections to Scalar Induced Gravitational Waves

SIGWs are generated at second-order soured by linear scalar perturbations,
which lead to $h_{ij}\sim\Phi^{2}$. Given that the scalar perturbations are
enhanced to $\sim\mathcal{O}(0.01-0.1)$ during the formation of PBHs, the
higher-order corrections to SIGWs are expected to be significant. The higher-
order corrections to SIGWs were first calculated by Yuan et al. where they
proposed a semi-analytical method to evaluate the SIGWs generated by a
$\delta$-spectrum during RD Yuan  _et al._ (2019a). Yuan et al. computed the
tensor modes sourced by the quadratic, cubic and biquadratic terms of the
linear perturbations during RD. Here we give the source term in Newton gauge
up to fourth-order in general cosmological background. Our convention for the
higher-order tensor modes are (up to fourth-order)

$\delta g_{ij}=a^{2}\left((1-2\Phi)\delta_{ij}+{1\over
2}(h_{ij}^{(2)}+h_{ij}^{(3)}+h_{ij}^{(4)})\right),$ (51)

and the perturbed perfect fluid up to third-order in Newton gauge is given by

$\displaystyle T_{00}^{(2)}\\!\\!\\!\\!\\!$ $\displaystyle=$
$\displaystyle\\!\\!\\!\\!\\!{1\over
2}\delta\rho^{(2)}+2\delta\rho^{(1)}\Phi^{(1)}+\rho\Phi^{(2)}+(P+\rho)v_{k}^{(1)}v^{k(1)}$
$\displaystyle T_{0i}^{(2)}\\!\\!\\!\\!\\!$ $\displaystyle=$
$\displaystyle\\!\\!\\!\\!\\!{1\over 2}\Big{(}-\rho
B_{i}^{(2)}-2\delta\rho^{(1)}B_{i}^{(1)}-2(\delta
P^{(1)}+\delta\rho^{(1)})v_{i}^{(1)}-(P+\rho)\left(4h_{ik}^{(1)}v^{k(1)}+2\left(\Phi^{(1)}-2\Psi^{(1)}\right)v_{i}^{(1)}+v_{i}^{(2)}\right)\Big{)}$
$\displaystyle T_{ij}^{(2)}\\!\\!\\!\\!\\!$ $\displaystyle=$
$\displaystyle\\!\\!\\!\\!\\!Ph_{ij}^{(2)}+2h_{ij}^{(1)}\delta
P^{(1)}+(P+\rho)\left(B_{i}^{(1)}+v_{i}^{(1)}\right)\left(B_{j}^{(1)}+v_{j}^{(1)}\right)+{1\over
2}\delta_{ij}\left(\delta P^{(2)}-4\delta
P^{(1)}\Psi^{(1)}-2P\Psi^{(2)}\right)$ $\displaystyle
T_{00}^{(3)}\\!\\!\\!\\!\\!$ $\displaystyle=$
$\displaystyle\\!\\!\\!\\!\\!\left(\delta
P^{(1)}+\delta\rho^{(1)}\right)v_{k}^{(1)}v^{k(1)}+(P+\rho)v^{k(1)}v_{k}^{(2)}+2(P+\rho)h_{lm}^{(1)}v^{l(1)}v^{m(1)}+{1\over
6}\delta\rho^{(3)}+{1\over 3}\rho\Phi^{(3)}$
$\displaystyle+2(P+\rho)\Phi^{(1)}v_{k}^{(1)}v^{k(1)}+\delta\rho^{(2)}\Phi^{(1)}+\delta\rho^{(1)}\Phi^{(2)}-2(P+\rho)\Psi^{(1)}v_{k}^{(1)}v^{k(1)}$
$\displaystyle T_{0i}^{(3)}\\!\\!\\!\\!\\!$ $\displaystyle=$
$\displaystyle\\!\\!\\!\\!\\!{1\over
6}\Bigg{[}(P+\rho)\Big{[}-6B_{i}^{(1)}B^{k(1)}v_{k}^{(1)}-6\left(h_{ik}^{(2)}v^{k(1)}+h_{ik}^{(1)}v^{k(2)}\right)+v_{i}^{(3)}-3\Big{(}(2B_{i}^{(1)}+v_{i}^{(1)})v_{k}^{(1)}v^{k(1)}$
$\displaystyle+4h_{ik}^{(1)}v^{k(1)}\Phi^{(1)}+v_{i}^{(2)}(\Phi^{(1)}-2\Psi^{(1)})-v_{i}^{(1)}\left((\Phi^{(1)})^{2}-\Phi^{(2)}+4\Phi^{(1)}\Psi^{(1)}+2\Psi^{(2)}\right)\Big{)}\Big{]}-3\left(\delta
P^{(1)}+\delta\rho^{(1)}\right)$
$\displaystyle\times\left(4h_{ik}^{(1)}v^{k(1)}+v_{i}^{(2)}+2v_{i}^{(1)}(\Phi^{(1)}-2\Psi^{(1)})\right)-3\left(\delta
P^{(2)}+\delta\rho^{(2)}\right)v_{i}^{(1)}-3\left(B_{i}^{(2)}\delta\rho^{(1)}+B_{i}^{(1)}\delta\rho^{(2)}\right)-\rho
B_{i}^{(3)}\Bigg{]}$ $\displaystyle T_{ij}^{(3)}\\!\\!\\!\\!\\!$
$\displaystyle=$
$\displaystyle\\!\\!\\!\\!\\!\frac{1}{6}\Bigg{(}2h_{ij}^{(3)}P+6\left(h_{ij}^{(2)}\delta
P^{(1)}+h_{ij}^{(1)}\delta P^{(2)}\right)+6(\delta
P^{(1)}+\delta\rho^{(1)})v_{i}^{(1)}v_{j}^{(1)}+3(P+\rho)\Big{[}v_{i}^{(1)}B_{j}^{(2)}+4h_{jk}^{(1)}v^{k(1)}v_{i}^{(1)}$
(52)
$\displaystyle+v_{i}^{(1)}v_{j}^{(2)}+v_{j}^{(1)}\left(B_{i}^{(2)}+4h_{ik}^{(1)}v^{k(1)}+v_{i}^{(2)}-8v_{i}^{(1)}\Psi^{(1)}\right)\Big{]}+3B_{j}^{(1)}\Big{[}2(B_{i}^{(1)}+v_{i}^{(1)})(\delta
P^{(1)}+\delta\rho^{(1)})$
$\displaystyle+(P+\rho)\times\left(B_{i}^{(2)}+4h_{ik}^{(1)}v^{k(1)}+v_{i}^{(2)}-4B_{i}^{(1)}\Phi^{(1)}-2v_{i}^{(1)}(\Phi^{(1)}+2\Psi^{(1)})\right)\Big{]}+3B_{i}^{(1)}\Big{[}$
$\displaystyle 2v_{j}^{(1)}(\delta
P^{(1)}+\delta\rho^{(1)})+(P+\rho)\Big{(}B_{j}^{(2)}+4h_{jk}^{(1)}v^{k(1)}+v_{j}^{(2)}-2v_{j}^{(1)}(\Phi^{(1)}+2\Psi^{(1)})\Big{)}\Big{]}$
$\displaystyle+\delta_{ij}\left[\delta P^{(3)}-2P\Psi^{(3)}-6(\delta
P^{(1)}\Psi^{(2)}+\delta P^{(2)}\Psi^{(1)})\right]\Bigg{)},$

where $v_{i}\equiv v_{i}^{\perp}+\partial_{i}v$ and $B_{i}$ are the velocity
perturbation and vector perturbation in Newton Gauge respectively.

The equation of motion for each $h_{ij}^{(n)}$ takes the same form as Eq. (27)
and the third-order source term reads

$\begin{split}S_{ij}^{(3)}=-\frac{4}{9\mathcal{H}^{4}\left(1+w\right)^{2}}&\Big{\\{}\left(\mathcal{H}\partial_{i}\Phi+\partial_{i}\Phi^{\prime}\right)\left(\mathcal{H}\partial_{j}\Phi+\partial_{j}\Phi^{\prime}\right)\Big{[}-2\left(1+c_{s}^{2}\right)\partial^{2}\Phi+6\mathcal{H}^{2}\left(c_{s}^{2}-w\right)\Phi\\\
&+3\mathcal{H}\left(3+2c_{s}^{2}+w\right)\Phi^{\prime}\Big{]}+3\mathcal{H}\left(1+w\right)\left(4\mathcal{H}\Phi-\Phi^{\prime}\right)\partial_{i}\Phi^{\prime}\partial_{j}\Phi^{\prime}\\\
&-3\mathcal{H}^{3}\left(1+w\right)\left(2\mathcal{H}\left(5+3w\right)\Phi-\Phi^{\prime}\right)\partial_{i}\Phi\partial_{j}\Phi\Big{\\}},\end{split}$
(53)

and the fourth-order source term is given by

$\begin{split}S_{ij}^{(4)}=&16\Phi^{3}\partial_{i}\partial_{j}\Phi-\frac{4}{27\mathcal{H}^{6}\left(1+w\right)^{3}}\Big{\\{}\left(\mathcal{H}\partial_{i}\Phi+\partial_{i}\Phi^{\prime}\right)\left(\mathcal{H}\partial_{j}\Phi+\partial_{j}\Phi^{\prime}\right)\Big{[}4\left(1+c_{s}^{2}\right)^{2}\left(\partial^{2}\Phi\right)^{2}\\\
&\qquad\quad-\mathcal{H}^{2}\left(9+5c_{s}^{2}+9w+9wc_{s}^{2}\right)\partial_{k}\phi\partial^{k}\Phi+4c_{s}^{2}\left(2\mathcal{H}\partial^{k}\Phi+\partial^{k}\Phi^{\prime}\right)\partial_{k}\Phi^{\prime}\\\
&\qquad\quad-6\mathcal{H}\left(1+c_{s}^{2}\right)\left(2\mathcal{H}\left(3+2c_{s}^{2}+w\right)\Phi+\left(5+4c_{s}^{2}+w\right)\Phi^{\prime}\right)\partial^{2}\Phi\\\
&\qquad\quad+18\mathcal{H}^{3}\left(7+4c_{s}^{4}+7c_{s}^{2}-wc_{s}^{2}+7w+4w^{2}\right)\Phi\Phi^{\prime}\\\
&\qquad\quad+36\mathcal{H}^{4}\left(1+c_{s}^{4}-2wc_{s}^{2}+2w+2w^{2}\right)\Phi^{2}+9\mathcal{H}^{2}\left(1+c_{s}^{2}\right)\left(5+4c_{s}^{2}+w\right)\left(\Phi^{\prime}\right)^{2}\Big{]}\\\
&+6\mathcal{H}\left(1+w\right)\partial_{i}\Phi^{\prime}\partial_{j}\Phi^{\prime}\Big{[}\left(1+c_{s}^{2}\right)\left(12\mathcal{H}^{3}\Phi^{2}+\Phi^{\prime}\partial^{2}\Phi-4\mathcal{H}\Phi\partial^{2}\Phi-3\mathcal{H}\left(\Phi^{\prime}\right)^{2}\right)\\\
&+3\mathcal{H}^{2}\left(-1+3c_{s}^{2}-4w\right)\Phi\Phi^{\prime}\Big{]}-3\mathcal{H}^{3}\left(1+w\right)\partial_{i}\Phi\partial_{j}\Phi\Big{[}\left(1+c_{s}^{2}\right)\left(2\Phi^{\prime}\partial^{2}\Phi-8\mathcal{H}\Phi\partial^{2}\Phi+18\mathcal{H}^{2}\Phi\Phi^{\prime}\right)\\\
&\qquad\quad-3\mathcal{H}\left(3+2c_{s}^{2}+w\right)\left(\Phi^{\prime}\right)^{2}+6\mathcal{H}^{3}\left(11+4c_{s}^{2}+22w+15w^{2}\right)\Phi^{2}\Big{]}\Big{\\}}.\end{split}$
(54)

Setting $w=c_{s}^{2}=1/3$ during RD, they return to the source terms used in
Yuan  _et al._ (2019a). In our convention, the density parameter for high-
order corrections takes the same form as Eq. (33). After defining the Fourier
transform of $S_{ij}$ to be

$\mathcal{S}({\bm{k},\eta})=-4e^{ij}S_{ij}(\bm{k})=\mathcal{S}^{(2)}({\bm{k},\eta})+\mathcal{S}^{(3)}({\bm{k},\eta})+\mathcal{S}^{(4)}({\bm{k},\eta}),$
(55)

then we can write the source terms as

$\displaystyle\mathcal{S}^{(2)}({\bm{k},\eta})$ $\displaystyle=$
$\displaystyle
4\int\frac{\mathrm{d}^{3}p}{(2\pi)^{3/2}}\mathbf{e}_{\bm{k}}(\bm{p},\bm{p})F^{(2)}(p,|\bm{k}-\bm{p}|,\eta)\Phi_{p}\Phi_{|\bm{k-p}|},$
(56) $\displaystyle\mathcal{S}^{(3)}({\bm{k},\eta})$ $\displaystyle=$
$\displaystyle
4\int\frac{\mathrm{d}^{3}p\mathrm{d}^{3}q}{(2\pi)^{3}}\mathbf{e}_{\bm{k}}(\bm{p},\bm{q})F^{(3)}(p,q,|\bm{k}-\bm{p}-\bm{q}|,\eta)\Phi_{p}\Phi_{q}\Phi_{|\bm{k-p-q}|},$
(57) $\displaystyle\mathcal{S}^{(4)}({\bm{k},\eta})$ $\displaystyle=$
$\displaystyle
4\int\frac{\mathrm{d}^{3}p\mathrm{d}^{3}q\mathrm{d}^{3}l}{(2\pi)^{9/2}}\Big{[}\mathbf{e}_{\bm{k}}(\bm{l},\bm{l})F^{(4)}_{1}(p,q,l,|\bm{k}-\bm{p}-\bm{q}-\bm{l}|,\eta)$
(59)
$\displaystyle\qquad\qquad\qquad+\mathbf{e}_{\bm{k}}(\bm{p},\bm{q})F^{(4)}_{2}(p,q,l,|\bm{k}-\bm{p}-\bm{q}-\bm{l}|,\eta)\Big{]}\Phi_{p}\Phi_{q}\Phi_{l}\Phi_{|\bm{k-p-
q-l}|},\qquad$

where we have defined $\mathbf{e}_{\bm{k}}(\bm{p},\bm{q})\equiv
e^{ij}(\bm{k})p_{i}q_{j}$. The contribution for the cross mode is omitted here
since the plus mode, and cross mode have the same energy density and we can
sum over the polarization mode at the final step. In Eq. (56), the time
evolution of $\Phi$ is absorbed in the transfer function, namely

$\begin{split}F^{(2)}(\bm{q}_{1},\bm{q}_{2},\eta)&=\frac{1}{\mathcal{H}^{2}}\left(3\mathcal{H}^{2}T_{q_{1}}T_{q_{2}}+\mathcal{H}(T_{q_{1}}T^{\prime}_{q_{2}}+T^{\prime}_{q_{1}}T_{q_{2}})+T^{\prime}_{q_{1}}T^{\prime}_{q_{2}}\right),\\\
F^{(3)}(\bm{q}_{1},\bm{q}_{2},\bm{q}_{3},\eta)&=\frac{1}{3\mathcal{H}^{4}}\Big{\\{}2\mathcal{H}^{2}T_{q_{1}}T_{q_{2}}\left[\left(18\mathcal{H}^{2}-q_{3}^{2}\right)T_{q_{3}}-6\mathcal{H}T_{q_{3}}^{\prime}\right]-2T_{q_{1}}^{\prime}T_{q_{2}}^{\prime}\left[\left(6\mathcal{H}^{2}+q_{3}^{2}\right)T_{q_{3}}+3\mathcal{H}T_{q_{3}}^{\prime}\right]\\\
&\quad-\mathcal{H}\left(T_{q_{1}}T_{q_{2}}^{\prime}+T_{q_{1}}^{\prime}T_{q_{2}}\right)\left(2q_{3}^{2}T_{q_{3}}+9\mathcal{H}T_{q_{3}}^{\prime}\right)\Big{\\}},\\\
F^{(4)}_{1}(\bm{q}_{1},\bm{q}_{2},\bm{q}_{3},\bm{q}_{4},\eta)&=16\,T_{q_{1}}T_{q_{2}}T_{q_{3}}T_{q_{4}},\\\
F^{(4)}_{2}(\bm{q}_{1},\bm{q}_{2},\bm{q}_{3},\bm{q}_{4},\eta)&=\frac{1}{36\mathcal{H}^{6}}\Big{\\{}\mathcal{H}^{2}T_{q_{1}}T_{q_{2}}\Big{[}-108\mathcal{H}^{2}T^{\prime}_{q_{3}}T^{\prime}_{q_{4}}+48\mathcal{H}^{2}\left(24\mathcal{H}^{2}+q_{3}^{2}+q_{4}^{2}\right)T_{q_{3}}T_{q_{4}}\\\
&\qquad\qquad\qquad+12\mathcal{H}\left(9\mathcal{H}^{2}-q_{3}^{2}\right)T_{q_{3}}T^{\prime}_{q_{4}}+12\mathcal{H}\left(9\mathcal{H}^{2}-q_{4}^{2}\right)T^{\prime}_{q_{3}}T_{q_{4}}\Big{]}\\\
&\qquad\quad+T^{\prime}_{q_{1}}T^{\prime}_{q_{2}}\Big{[}72\mathcal{H}^{2}T^{\prime}_{q_{3}}T^{\prime}_{q_{4}}-48\mathcal{H}^{2}\left(6\mathcal{H}^{2}+q_{3}^{2}+q_{4}^{2}\right)T_{q_{3}}T_{q_{4}}\\\
&\qquad\qquad\qquad+12\mathcal{H}\left(3\mathcal{H}^{2}+q_{3}^{2}\right)T_{q_{3}}T^{\prime}_{q_{4}}+12\mathcal{H}\left(3\mathcal{H}^{2}+q_{4}^{2}\right)T^{\prime}_{q_{3}}T_{q_{4}}\Big{]}\\\
&\qquad\quad-\left(\mathcal{H}T_{q_{1}}+T^{\prime}_{q_{1}}\right)\left(\mathcal{H}T_{q_{2}}+T^{\prime}_{q_{2}}\right)\Big{[}12\mathcal{H}\left(21\mathcal{H}^{2}+5q_{3}^{2}\right)T_{q_{3}}T^{\prime}_{q_{4}}\\\
&\qquad\qquad\qquad+12\mathcal{H}\left(21\mathcal{H}^{2}+5q_{4}^{2}\right)T^{\prime}_{q_{3}}T_{q_{4}}+180\mathcal{H}^{2}T^{\prime}_{q_{3}}T^{\prime}_{q_{4}}\\\
&\qquad\qquad\qquad+8\left(18\mathcal{H}^{4}+9\mathcal{H}^{2}q_{3}^{2}+9\mathcal{H}^{2}q_{4}^{2}+2q_{3}^{2}q_{4}^{2}\right)T_{q_{3}}T_{q_{4}}\\\
&\qquad\qquad\qquad+3\bm{q}_{3}\cdot\bm{q}_{4}\left(11\mathcal{H}^{2}T_{q_{3}}T_{q_{4}}-\mathcal{H}T_{q_{3}}T^{\prime}_{q_{4}}-\mathcal{H}T^{\prime}_{q_{3}}T_{q_{4}}-T^{\prime}_{q_{3}}T^{\prime}_{q_{4}}\right)\Big{]}\Big{\\}},\end{split}$
(60)

where we set $w=c_{s}^{2}=1/3$. The notation $T_{k}$ is short for the time
evolution function, $T_{\Phi}(k\eta)$. Applying the Green’s function method,
the solution to $h_{ij}^{(n)}$ takes the same form as Eq. (29). The next-order
corrections to $\Omega_{\mathrm{GW}}$ comes from two parts. The first part is
the coupling of $S_{ij}^{(3)}$ itself,

$\begin{split}\Omega_{\mathrm{GW}}^{I}(k,\eta)=&\frac{k^{3}}{384\pi^{2}}\left({k\over\mathcal{H}}\right)^{2}\int\mathrm{d}\tilde{\eta}_{1}\frac{a(\tilde{\eta}_{1})}{a(\eta)}g_{k}(\eta;\tilde{\eta}_{1})\int\tilde{\eta}_{2}\frac{a(\tilde{\eta}_{2})}{a(\eta)}g_{k}(\eta,\tilde{\eta}_{2})\int\mathrm{d}^{3}p\mathrm{d}^{3}q\frac{P_{\Phi}(|\bm{p}-\bm{q}|)}{|\bm{p}-\bm{q}|^{3}}\frac{P_{\Phi}(q)}{q^{3}}\frac{P_{\Phi}(|\bm{k}-\bm{p}|)}{|\bm{k}-\bm{p}|^{3}}\\\
&\quad\times\Big{[}2\bm{e}_{\bm{k}}(\bm{p}-\bm{q},\bm{q})^{2}F^{(3)}(q,|\bm{p}-\bm{q}|,|\bm{k}-\bm{p}|,\tilde{\eta}_{1})F^{(3)}(q,|\bm{p}-\bm{q}|,|\bm{k}-\bm{p}|,\tilde{\eta}_{2})\\\
&\qquad+4\bm{e}_{\bm{k}}(\bm{p},\bm{p}-\bm{q})\bm{e}_{\bm{k}}(\bm{p},\bm{q})F^{(3)}(|\bm{k}-\bm{p}|,|\bm{p}-\bm{q}|,q,\tilde{\eta}_{1})F^{(3)}(q,|\bm{k}-\bm{p}|,|\bm{p}-\bm{q}|,\tilde{\eta}_{2})\Big{]},\end{split}$
(61)

while the second part is the coupling of $S_{ij}^{(2)}$ and $S_{ij}^{(4)}$

$\begin{split}\Omega_{\mathrm{GW}}^{II}(k,\eta)=&\frac{k^{3}}{384\pi^{2}}\left({k\over\mathcal{H}}\right)^{2}\int\mathrm{d}\tilde{\eta}_{1}\frac{a(\tilde{\eta}_{1})}{a(\eta)}g_{k}(\eta;\tilde{\eta}_{1})\int\tilde{\eta}_{2}\frac{a(\tilde{\eta}_{2})}{a(\eta)}g_{k}(\eta,\tilde{\eta}_{2})\int\mathrm{d}^{3}p\mathrm{d}^{3}q\frac{P_{\Phi}(p)}{p^{3}}\frac{P_{\Phi}(q)}{q^{3}}\frac{P_{\Phi}(|\bm{k}-\bm{p}|)}{|\bm{k}-\bm{p}|^{3}}\\\
&\times\mathbf{e}_{\bm{k}}(\bm{p},\bm{p})F^{(2)}(\bm{p},\bm{k}-\bm{p},\tilde{\eta}_{1})\Big{\\{}6\left[\mathbf{e}_{\bm{k}}(\bm{p},\bm{p})+\mathbf{e}_{\bm{k}}(\bm{q},\bm{q})\right]F_{1}^{(4)}(\bm{p},\bm{k}-\bm{p},\bm{q},\bm{q},\tilde{\eta}_{2})\\\
&-2\mathbf{e}_{\bm{k}}(\bm{p},\bm{p})F_{2}^{(4)}(\bm{p},\bm{k}-\bm{p},\bm{q},-\bm{q},\tilde{\eta}_{2})-8\mathbf{e}_{\bm{k}}(\bm{p},\bm{q})F_{2}^{(4)}(\bm{p},\bm{q},\bm{q},\bm{k}-\bm{p},\tilde{\eta}_{2})\\\
&-2\mathbf{e}_{\bm{k}}(\bm{q},\bm{q})F_{2}^{(4)}(\bm{q},\bm{q},\bm{p},\bm{k}-\bm{p},\tilde{\eta}_{2})\Big{\\}}.\end{split}$
(62)

The next-order corrections of $\Omega_{\mathrm{GW}}$ is the sum of Eq. (61)
and Eq. (62). For a $\delta$-spectrum, these two expressions can be further
simplified after integrate over the $\delta$ function (see Yuan  _et al._
(2019a) for results). After considering the higher-order corrections, the deep
valley at $k=\sqrt{2/3}k_{*}$ generated at second-order will be smoothed to a
finite value. Moreover, the cut-off frequency will be extended from $2k_{*}$
to $3k_{*}$. See Fig. 4.

The above calculation only focus on the GWs induced by the linear-order scalar
perturbations. However, the higher-order scalar, vector and tensor
perturbations will contribute to the higher-order SIGWs. Recently, the SIGWs
induced by the second-order perturbations were analytically studied in Zhou
_et al._ (2021). But the complete calculation of the next order correction to
SIGWs has not been explored yet.

###  3.3 Scalar Induced Gravitational Waves within Non-Gaussianities

GWs generated by non-Gaussian scalar perturbations were first estimated by
Nakama Nakama  _et al._ (2017) where they found that the amplitude of SIGWs
could be suppressed by non-Gaussianities to several orders of magnitude. After
Espinosa  _et al._ (2018); Kohri and Terada (2018a) proposed the semi-
analytical method to calculate SIGWs, Cai  _et al._ (2019a); Unal (2019)
calculate the non-Gaussian effects on SIGWs by considering a local-type non-
Gaussianities up to second-order (or $F_{\mathrm{NL}}$-order). They both
argued that SIGWs within non-Gaussianities have observable signatures which
can be used to probe primordial non-Gaussianities. However, in a more recent
study Yuan and Huang (2020), Yuan et al. revisited the non-Gaussian effects
and extended the calculation to third-order (or $G_{\mathrm{NL}}$-order non-
Gaussianities). They argue that all the non-Gaussian effects are degenerate
with the power spectrum. Hence it is impossible to read any information about
non-Gaussianities only through the signal of SIGWs.

To see this more clear, let us go back to the first line of Eq. (34), where
one needs to compute the two-point correlator of the source terms. This will
lead to the four-point correlator of $\Phi$, namely
$\left\langle\Phi_{\bm{p}}\Phi_{\bm{k-p}}\Phi_{\bm{q}}\Phi_{\bm{k^{\prime}-q}}\right\rangle$.
If $\Phi$ obeys Gaussian distribution, then the four-point correlator can be
simplified to three two-point correlators through Wick’s theorem. However, the
general formula for the four-point correlator is

$\displaystyle\left\langle\Phi_{\bm{p}}\Phi_{\bm{k-p}}\Phi_{\bm{q}}\Phi_{\bm{k^{\prime}-q}}\right\rangle$
$\displaystyle=$
$\displaystyle\left\langle\Phi_{\bm{p}}\Phi_{\bm{k-p}}\Phi_{\bm{q}}\Phi_{\bm{k^{\prime}-q}}\right\rangle_{c}+\left\langle\Phi_{\bm{p}}\Phi_{\bm{k-p}}\right\rangle\left\langle\Phi_{\bm{q}}\Phi_{\bm{k-q}}\right\rangle+\left\langle\Phi_{\bm{p}}\Phi_{\bm{q}}\right\rangle\left\langle\Phi_{\bm{k-p}}\Phi_{\bm{k^{\prime}-q}}\right\rangle$
(63) $\displaystyle+$
$\displaystyle\left\langle\Phi_{\bm{p}}\Phi_{\bm{k^{\prime}-q}}\right\rangle\left\langle\Phi_{\bm{k-p}}\Phi_{\bm{q}}\right\rangle$

where we assume that $\left\langle\Phi\right\rangle=0$ and
$\left\langle\Phi_{\bm{p}}\Phi_{\bm{k-p}}\Phi_{\bm{q}}\Phi_{\bm{k^{\prime}-q}}\right\rangle_{c}$
is the connected four-point correlation function (4PCF) (or called the fourth
cumulant in statistics) which is related to the trispectrum,
$\mathcal{T}_{\Phi}$, by

$\left\langle\Phi_{\bm{p}}\Phi_{\bm{k-p}}\Phi_{\bm{q}}\Phi_{\bm{k^{\prime}-q}}\right\rangle_{c}=(2\pi)^{3}\delta(\bm{k}+\bm{k^{\prime}})\mathcal{T}_{\Phi}(\bm{p},\bm{q},\bm{k},\bm{k^{\prime}}).$
(64)

The connected 4PCF will vanish if $\Phi$ is a Gaussian variable. Otherwise it
will contribute to be SIGWs Unal (2019). In Unal (2019); Cai  _et al._
(2019a), the authors consider a local-type non-Gaussianities up to
$F_{\mathrm{NL}}$-order and calculate the SIGWs in the absence of the
connected 4PCF. After that, the $G_{\mathrm{NL}}$-order is studied by Yuan and
Huang (2020) still in the absence of the connected 4PCF. Recently, Atal and
Domènech (2021) considered the leading order of the connected 4PCF and Adshead
_et al._ (2021) included the complete contribution of the connected 4PCF.
Although non-Gausiannities can alter the waveform of the SIGW signal, one can
not read any information about non-Gaussianities only through the signal of
$\Omega_{\mathrm{GW}}(f)$. This is because the non-Gaussian effects are
absorbed into the total 4PCF and there will be a degeneracy between the non-
Gaussian effects and the total 4PCF. Therefore, one can not tell from the
total 4PCF whether there is a evidence for non-Gaussianities (or for the
connected 4PCF).

###  3.4 Gauge issue of Scalar Induced Gravitational Waves

Although the tensor mode are gauge invariant at first-order, they fail to
remain gauge invariant at second-order (see, e.g., Noh and Hwang (2003)).
Therefore, a natural question arises: Is the energy density of SIGWs gauge
dependent or not?

For a long time, Newton gauge has been a commonly used gauge to compute SIGWs
since the degrees of freedom are fixed completely in this gauge, and there is
no residual gauge freedom. In Newton gauge, the scalar mode $\phi$ is the just
Bardeen potential $\Phi$ and all the degrees of freedom are fixed, which make
Newton gauge the most mathematically convenient gauge for evaluating SIGWs.
SIGWs in other gauges were first calculated numerically by Hwang et al. Hwang
_et al._ (2017) where they investigated $\Omega_{\mathrm{GW}}$ in uniform
expansion gauge, comoving gauge and uniform curvature gauge and they found
that $\Omega_{\mathrm{GW}}$ is gauge dependent. After that, the semi-
analytical method was developed, and Gong calculated the SIGWs in comoving
gauge during MD Gong (2019). Gong’s result showed that $h_{ij}$ increases as
$\eta^{2}$ in comoving gauge. Another study by Tomikawa and Kobayashi appeared
at the same time Tomikawa and Kobayashi (2019), where they investigate the
SIGWs in comoving gauge and uniform curvature gauge. They found that the
results in comoving gauge increase with time for $w\geq 0$ while the result in
uniform curvature gauge is identical with that in Newton gauge for $w>0$.
Moreover, they found that the results are all different in Newton gauge,
comoving gauge, and uniform curvature gauge if $w=0$.

Due to the gauge dependence of $\Omega_{\mathrm{GW}}$, it is natural to ask
which gauge is relevant in interpreting the observations. This question is
discussed by De Luca et al. De Luca  _et al._ (2019a). They calculated the
SIGWs in synchronous gauge since they argued that the sensitivity curves of
LISA are given in that gauge. By neglecting the $E$ mode but keeping its
derivatives in the source term , De Luca et al. found that the SIGWs in
synchronous gauge are identical to that of Newton gauge. Their results are
later confirmed by Yuan  _et al._ (2019c); Inomata and Terada (2019). However,
the calculation of SIGWs in synchronous gauge is a a tough nut since there is
residual gauge freedom in $E$ mode. The complete calculation in synchronous
gauge and which gauge is relevant for observations are still open questions.

In synchronous gauge, $\delta g_{00}=\delta g_{0i}=0$, corresponding to
$\tilde{\phi}=\tilde{B}=0$ in Eq. (19). This leads to Lu  _et al._ (2020)

$\displaystyle T^{S}(\eta)$ $\displaystyle=$ $\displaystyle-{1\over
a}\left(\int_{0}^{\eta}a(\tilde{\eta})\phi(\tilde{\eta})\mathrm{d}\tilde{\eta}-\mathcal{C}_{1}(\bm{x})\right)$
(65) $\displaystyle L^{S}(\eta)$ $\displaystyle=$
$\displaystyle\int_{0}^{\eta}[T^{S}(\tilde{\eta})-B(\tilde{\eta})]\mathrm{d}\tilde{\eta}+\mathcal{C}_{2}(\bm{x}),$
(66)

where $\mathcal{C}_{1}$ and $\mathcal{C}_{2}$ are two arbitrary spatial
functions. The presence of $\mathcal{C}_{1}$ and $\mathcal{C}_{2}$ come from
the residual gauge freedom in synchronous gauge. To determine the time slicing
and the spatial coordinate on the hypersurface, one has to fix
$\mathcal{C}_{1}$ and $\mathcal{C}_{2}$. To see the impacts of residual gauge
freedom more clearly, we derive the equation of motion for the scalar modes
during RD, namely

$\displaystyle 2\mathcal{H}E^{\prime}+E^{\prime\prime}+\psi$ $\displaystyle=$
$\displaystyle 0,$ (67) $\displaystyle
6\psi^{\prime\prime}+2\mathcal{H}\left(9\psi^{\prime}-4\partial^{2}E\right)-3\partial^{2}E^{\prime\prime}-5\partial^{2}\psi$
$\displaystyle=$ $\displaystyle 0.$ (68)

After some algebra, one can get the equation of motion for $E$, (see also Lu
_et al._ (2020))

$x^{3}T_{E}^{\prime\prime\prime\prime}(x)+5x^{2}T_{E}^{\prime\prime\prime}(x)+\left(2+{x^{2}\over
3}\right)xT_{E}^{\prime\prime}(x)-\left(2-{x^{2}\over
3}\right)T_{E}^{\prime}(x)=0,$ (69)

where the transfer function for $E$ and $\psi$ are defined as

$\displaystyle k^{2}E(\bm{k})$ $\displaystyle\equiv$
$\displaystyle\Phi_{k}T_{E}(k\eta)$ (70) $\displaystyle\psi(\bm{k})$
$\displaystyle\equiv$ $\displaystyle\Phi_{k}T_{\psi}(k\eta).$ (71)

In our convention, the boundary condition for the transfer function should be
$T_{\psi}(x\to 0)=3/2$ in order to match the relation between curvature
perturbation and the Bardeen potential. The general solution for Eq. (69) is

$T_{E}(x)=\mathcal{C}_{3}+\mathcal{C}_{4}\left(\mathrm{Ci}(z)-{\sin(z)\over
z}\right)+\mathcal{C}_{5}\ln(z)+\mathcal{C}_{6}\left(\mathrm{Si}(z)+{\cos(z)\over
z}\right),$ (72)

where $z\equiv x/\sqrt{3}$ and $\mathcal{C}_{i},(i=3,4,5,6)$ are integrate
constants. One can check that under the residual gauge transformation,

$\displaystyle T$ $\displaystyle=$ $\displaystyle-{\mathcal{C}_{7}\over x}$
(73) $\displaystyle L$ $\displaystyle=$ $\displaystyle-\mathcal{C}_{7}\ln
x+\mathcal{C}_{8},$ (74)

the gauge condition $\tilde{\phi}=\tilde{B}=0$ still satisfy. Combining Eq.
(72) and Eq. (74), we see that $\mathcal{C}_{8}$ can be absorbed by
$\mathcal{C}_{3}$ while $\mathcal{C}_{7}$ corresponds to $\mathcal{C}_{5}$. In
other words, the constant $\mathcal{C}_{3}$ and the logarithm term are two
independent pure gauge modes. One can arbitrarily choose $\mathcal{C}_{3}$ and
$\mathcal{C}_{5}$ to work in a specific hypersurface. After selecting the
hypersurface, the solution can be well determined after performing the only
physical condition, $T_{\psi}(x\to 0)=3/2$. Here comes a big problem in
synchronous gauge. The solution to $E$ mode will definitely diverge either at
the initial time or at a late time. For instance, let’s assume that $T_{E}(x)$
have a finite value at $x=0$. Notice that the Taylor expansion of cosine
integral function at $x=0$ is
$\mathrm{Ci}(x)=\gamma+\ln(x)+\mathcal{O}(x^{2})$, where $\gamma$ is the Euler
gamma constant. Therefore, one must choose $\mathcal{C}_{5}=-\mathcal{C}_{4}$
so that the logarithm term can be cancelled by the $\mathrm{Ci}(x)$ term. By
applying the boundary condition $T_{\psi}(0)=3/2$, one gets
$\mathcal{C}_{6}=0$ and $\mathcal{C}_{4}=-\mathcal{C}_{5}=9$ and the finite
value at $x=0$ is given by $T_{E}(0)=\mathcal{C}_{3}+9\gamma$. However, under
this situation, we see that the $T_{E}(x)$ at a late time will diverge
according to the logarithm term in Eq. (72). The consequence of the divergence
is that the $E$ mode will induce the GWs continuously during RD, resulting in
the divergence of $\Omega_{\mathrm{GW}}$ for SIGWs. On the other hand, if we
let the gauge mode $\mathcal{C}_{5}$ vanish, the presence of the
$\mathrm{Ci}(x)$ term will result in the divergence at the initial time. A
physical interpretation is that if the observer requires a converged initial
condition, the geodesics of the observer using as a reference may collapse or
cross each other. So one has to change the residual gauge freedom to adjust
the reference system. This makes the calculation of SIGWs in synchronous gauge
very complicated.

In De Luca  _et al._ (2019a); Yuan  _et al._ (2019c); Inomata and Terada
(2019), the authors calculated the SIGWs in the absence of $E$ mode but
considering its derivatives where the gauge freedom are completely fixed. In
this situation, the energy density of SIGWs agrees with that in Newton gauge.
Moreover, it was shown in Lu  _et al._ (2020) that if throwing away the pure
gauge modes, the result returns to that in Newton gauge. These studies
indicate that the large divergence of SIGWs in synchronous gauge comes from
the residual gauge freedom and the synchronous gauge seems to be ill-defined
for calculating the SIGWs. Readers interested in the calculation of different
gauge can refer to Lu  _et al._ (2020); Ali  _et al._ (2020); Inomata (2020).

Since the standard procedure to calculate the energy density of $h_{ij}$
brings divergence, one way to tackle the problem is to find a new gauge
invariant quantity for the GWs. One can construct infinite gauge invariant
variables for scalar and tensor perturbations. However, the point is to find a
quantity that can interpret the physical world and the measurements. This idea
has been tried by Zhang, Wang and Zhu in Chang  _et al._ (2020a, b). By
applying the technique of Lie derivative, they constructed gauge invariant
second-order GWs in synchronous gauge, which is related to the measurement,
and their results coincide with that computed in Newton gauge. They also used
this idea to evaluate the energy density of SIGWs in uniform density gauge,
which is supposed to diverge as $\eta^{6}$. They found the new gauge invariant
second-order GWs converged and were identical with the one in Newton gauge
Chang  _et al._ (2020c). In a more recent study Domènech and Sasaki (2020), it
is shown that the current $\Omega_{\mathrm{GW}}$ is well defined in most of
gauges if one takes the sub-horizon limit, namely $\eta\to\infty$. The authors
found that $\Omega_{\mathrm{GW}}$ will be the same as that in Newton gauge if
the trace part of the metric perturbation at sub-horizon limit is the order of
or smaller than the scalar perturbation $\Phi$ in Newton gauge (see Eq. (3.2)
or Eq. (3.5) in Domènech and Sasaki (2020)). They also pointed out that, for
particular choices of the residual gauge freedom, synchronous gauge can give
the same result in Newton gauge.

Although the synchronous gauge is relevant for observations, Newton gauge is
still the most popular gauge since it not only gives the same result as
synchronous gauge but also is the simplest gauge to perform the calculation.

##  4 Searching for Primordial Black Hole Dark Matter Using Scalar Induced
Gravitational Waves

SIGWs provide the most efficient way so far to search for PBH DM. This is due
to the fact that PBHs might from the peak of scalar perturbations, which is
extremely sensitive to the amplitude of the power spectrum, $A$. Let’s
consider a monochromatic power spectrum for example, described by Eq. (40).
The variance of the perturbations generated by this spectrum is $\sigma=\int
P_{\zeta}(k)\mathrm{d}\ln k=A$. Using the Press-Schechter formalism, Eq. (11),
the mass function can be approximated as

$\beta={1\over 2}\mathrm{erfc}\left({\nu_{c}\over\sqrt{2}}\right)=\sqrt{1\over
2\pi}{\mathrm{e}^{-\nu_{c}^{2}/2}\over\nu_{c}}\propto\sqrt{A}~{}\mathrm{e}^{-{\Delta_{c}^{2}\over
2A}},$ (75)

which is exponentially dependent on $A$. On the other hand, the inevitably
generated SIGWs satisfy $h_{ij}\propto A$. Therefore, a small change in $A$
would change $f_{\mathrm{pbh}}$ by several orders of magnitude.

To further quantify the power of SIGWs, we shall estimate the expected signal-
to-noise ratio (SNR) by the GW detector, which is evaluated as Allen and
Romano (1999); Thrane and Romano (2013),

$\rho^{2}\\!=T\\!\int\\!\mathrm{d}f\frac{\Gamma(f)^{2}S_{h}(f)^{2}}{\left[{1\over
25}+\Gamma(f)^{2}\right]S_{h}(f)^{2}+P_{n}(f)^{2}+{2\over
5}S_{h}(f)P_{n}(f)},$ (76)

The strain power spectral density is given by
$S_{h}(f)=3H_{0}^{2}\Omega_{\mathrm{GW},0}/(2\pi^{2}f^{3})$. $\Gamma(f)$ is
the overlap function and $P_{n}(f)$ is the noise power spectral density. For
LISA, $\Gamma(f)=R(f)$ where $R(f)$ is the signal transfer function of LISA
and its expression should be computed numerically Larson  _et al._ (2000) but
it can be well fit by Robson  _et al._ (2019)

$R(f)={3\over 10(1+0.6(f/f_{\star})^{2})},$ (77)

where $f_{\star}=c/(2\pi L)=19.09~{}\mathrm{mHz}$ is the transfer frequency
and $L=2.5~{}\mathrm{Gm}$ for the current LISA design. On the other hand,
$P_{n}(f)$ can be approximated by Robson  _et al._ (2019)

$P_{n}(f)={P_{\mathrm{oms}}\over
L^{2}}+2(1+\cos^{2}(f/f_{\star})){P_{\mathrm{acc}}\over{(2\pi f)^{4}L^{2}}}.$
(78)

The optical metrology noise spectrum and the acceleration noise is
parameterized by Robson  _et al._ (2019)

$\displaystyle P_{\mathrm{oms}}(f)$ $\displaystyle=$ $\displaystyle(1.5\times
10^{-11}\mathrm{m})^{2}\left(1+\left({2\mathrm{mHz}\over
f}\right)^{4}\right)\mathrm{Hz}^{-1}$ (79) $\displaystyle P_{\mathrm{acc}}(f)$
$\displaystyle=$ $\displaystyle(3\times
10^{-15}~{}\mathrm{m~{}s^{-2}})^{2}\left(1+\left({0.4\mathrm{mHz}\over
f}\right)^{2}\right)\left(1+\left({f\over
8\mathrm{mHz}}\right)^{4}\right)\mathrm{Hz}^{-1}.$ (80)

For pulsar timing arrays (PTA), the expressions are
$\Gamma(f)=R(f)=1/(12\pi^{2}f^{2})$ and $P_{n}(f)=2\Delta t\sigma^{2}$ where
$\Delta t$ is the observation time and $\sigma$ is the root-mean-square timing
noise. For PTAs, each of the millisecond pulsars can be regarded as a single
GW detector, and one has to average the spatial contribution such that Allen
and Romano (1999); Siemens  _et al._ (2013)

$\rho^{2}=2T\sum_{I,J}^{M}\zeta_{IJ}^{2}\left[\int\\!\\!\mathrm{d}f\frac{R(f)^{2}S_{h}(f)^{2}}{\left[{1\over
25}+R(f)^{2}\right]S_{h}(f)^{2}+P_{n}(f)^{2}+{2\over
5}S_{h}(f)P_{n}(f)}\right],$ (81)

where $\zeta_{IJ}$ is the normalized Hellings and Downs coefficient for
pulsars $I$ and $J$. Assuming the pulsars are distributed homogeneously on the
sky, then $\zeta_{IJ}$ takes the form

$\zeta_{IJ}={N(N-1)}\frac{4.74}{20\times 19},$ (82)

where $N$ is the number of pulsars. Combing the above equations, we can
estimate the expected SNR from monochromatic PBHs and scan the mass of PBHs.
The result is shown in Fig. 5. PBHs lighter than $\lesssim 10^{-18}M_{\odot}$
have already evaporated through Hawking radiation by today and is not shown in
the figure.

If the PTAs and LISA fail to detect the SIGWs from PBHs, an upper limit can be
placed for $f_{\mathrm{pbh}}$. In Wang  _et al._ (2019b), Wang et al. found
the constraints on $f_{\mathrm{pbh}}$ can reach $10^{-13}$ for
$m_{\rm{pbh}}\in[10^{-8},1]M_{\odot}$. Chen et al. Chen  _et al._ (2019b)
searched the SIGW signals in the NANOGrav 11-yr data set for monochromatic
PBHs. They scan the mass of PBHs and dit not find the SIGW signal. Therefore
they placed an upper limit on $f_{\mathrm{pbh}}$, see Fig. 6. As shown in the
figure, the constraints from SIGW are several orders of magnitude better than
the other constraints. Recently, NANOGrav has reported strong evidence for a
common-spectrum process modelled as power-law in the 12.5-year dataset.
Various models have been put forward to explain the signal assumed to be
gravitational waves. SIGWs generated by a broad and flat power spectrum seems
to give one possible explanation De Luca  _et al._ (2021); Sugiyama  _et al._
(2021) among the various cosmological and astrophysical models. See also
Vaskonen and Veermäe (2021); Kohri and Terada (2021); Domènech and Pi (2020);
Inomata  _et al._ (2020d). However, further analysis indicates that the signal
is preferred to be scalar transverse modes in the general metric theory rather
than the gravitational waves predicted in general relativity Chen  _et al._
(2021).

Another study was made by Kapadia et al. Kapadia  _et al._ (2020b) where they
explored the SIGW signal in the LIGO data. Kapadia et al. reported a null
detection of SIGW signal and thus placing a severe constraint on
$f_{\mathrm{pbh}}$ to be less than a few parts in million. Their work provides
another independent test in addition to the $\gamma$-ray burst Carr  _et al._
(2010) in this mass window.

Up to now there are various independent cosmological observations constraining
the fraction of PBHs in Fig. 7. For $m_{\rm{pbh}}\lesssim 10^{-18}M_{\odot}$,
PBHs are severely constrained through the null detections of extra-galactic
Gamma-ray background (EGB) from PBH evaporation Carr  _et al._ (2010). On the
other hand, PBHs heavier than $\sim 10^{3}M_{\odot}$ are tightly constrained
by the CMB observations through the accreting PBHs Ali-Haimoud and
Kamionkowski (2017); Aloni  _et al._ (2017); Horowitz (2016); Chen  _et al._
(2016). In the mass range $m_{\rm{pbh}}\in[10^{-18},10^{3}]M_{\odot}$,
$f_{\mathrm{pbh}}$ are constrained to no more than a few in thousand by
numerous astrophysical observations. For example, microlensing events such as
Subaru/HSC Niikura  _et al._ (2019a), Kepler Griest  _et al._ (2013), OGLE
Niikura  _et al._ (2019b) and EROS/MACHO Tisserand  _et al._ (2007) and the
dynamical heating of ultra-faint dwarf galaxies Brandt (2016). On the other
hand, GWs physics also placed broad constraints on $f_{\mathrm{pbh}}$, such as
the merger events from binary PBHs in the subsolar mass range Abbott  _et al._
(2018); Magee  _et al._ (2018); Chen and Huang (2019); Abbott  _et al._ (2019)
and the null detections of stochastic GW background (SGWB) from binary PBHs
with LIGO Wang  _et al._ (2018); Chen and Huang (2019). Although the existence
of white dwarfs (WDs) in our local galaxy have also placed constraints on
$f_{\mathrm{pbh}}$ Graham  _et al._ (2015). However, this constraint has been
challenged by a recent paper Montero-Camacho  _et al._ (2019), and we label
the constraint from WDs with the dashed line in Fig. 7. Despite the above
constraints, PBHs in a substantial window
$\sim[10^{-16},10^{-14}]\cup[10^{-13},10^{-12}]M_{\odot}$ are still allowed to
account for all of the DM.

Although SIGWs place the currently severest constraint on $f_{\mathrm{pbh}}$,
it is based on the standard PBH formation model in which PBHs are formed from
the critical collapse of overdensity. However, other mechanism of the PBH
formation has been proposed (see e.g., Cotner  _et al._ (2019)) and the
constraints from SIGWs are not valid for those PBHs.

##  5 Summary and Outlook

PBHs can represent the DM in our Universe and explain the merger events
detected by LIGO/VIRGO if $f_{\mathrm{pbh}}\sim 10^{-3}$ in that mass range.
We focus on the standard PBH formation model, where the PBHs are formed
through the collapse of the overdensity in the very early Universe. We have
review several aspects of PBH formation, and the SIGWs produced from PBHs
studied over the past decades.

First of all, We review the calculation of the PBH formation. The formation of
PBHs depends on the power spectrum of the scalar perturbations, the window
function, the primordial non-Gaussianities, intrinsic non-Gaussianities and
the critical value beyond which PBH can form.

Secondly, we consider the SIGWs inevitably generated during the formation of
PBHs. We review the relevant works concerning SIGWs over the past decades.
SIGW is a solid prediction of general relativity. But due to the unknown of
primordial scalar power spectrum on small scales, the waveform of the SIGW
background is model dependent which make the detection of SIGWs very
difficult. Fortunately, the logarithm scaling in the infrared region is a
model-independent feature for SIGWs, making SIGWs unique from other SGWB. We
also review the impacts from the cosmological background, higher-order
corrections to SIGWs, primordial non-Gaussianities. Finally, we discuss the
gauge issues of SIGWs. On the other hand, the presence of primordial non-
Gaussianities is an essential aspect since it could dramatic change the
waveform and the amplitude of the SIGWs. For local-type non-Gaussianities, as
long as monochromatic PBH in the allowing mass window represents all the DM,
then LISA and PTAs should detect the signals irrespective of
$F_{\mathrm{NL}}$. However, up to $G_{\mathrm{NL}}$-order, this conclusion no
longer holds. The $G_{\mathrm{NL}}$-order could further suppress the SIGWs
than $F_{\mathrm{NL}}$-order and thus avoiding the detection of LISA. Another
important impact of non-Gaussianities is that the power spectrum and the non-
Gaussianities are degenerate. As a result, one cannot tell the presence of
non-Gaussianities only through the SIGW signal unless considering another
independent observation.

So far, SIGW has placed the most severe constraint on $f_{\mathrm{pbh}}$ Chen
_et al._ (2019b). However, this results from monochromatic PBHs in the absence
of primordial non-Gaussianities and neglecting the QCD phase transition.
During the QCD epoch, the equation of state and the sound speed decrease
slightly from $1/3$, thus changing the waveform of SIGWs. This will affect the
PBHs within the PTA frequency band. Searching the SIGW signal for realistic
models such as considering the QCD phase transition, a specific inflation
model, and its loop correction to the power spectrum is necessary in this
field.

Over the past decade, SIGWs from PBH DM have been widely studied. SIGW
provides a promising tool to verify or falsify the PBH DM hypothesis. We are
expected to witness intriguing progress in the coming decades.

Acknowledgments. Cosmological perturbations are derived using the xPand Pitrou
_et al._ (2013) package. We would like to thank Zu-Cheng Chen, Sai Wang, Misao
Sasaki and Guillem Domènech for useful discussion. We also acknowledge the use
of HPC Cluster of ITP-CAS. This work is supported by the National Key Research
and Development Program of China Grant No.2020YFC2201502, the grants from NSFC
(grant No. 11975019, 11690021, 11991052, 12047503), the Key Research Program
of the Chinese Academy of Sciences (Grant NO. XDPB15), Key Research Program of
Frontier Sciences, CAS, Grant NO. ZDBS-LY-7009, and the science research
grants from the China Manned Space Project with NO. CMS-CSST-2021-B01.

Author contribution. Q.G.H. organizes the whole project. C.Y. performs the
main calculations and write the paper.

Declaration of interests. The authors declare no competing interests.

Figure 1:  SIGWs generated by different power spectrum. The parameters for the
log-normal spectrum is $\sigma_{*}=0.5$. For the box spectrum, we set
$k_{\mathrm{min}}=0.1$ and $k_{\mathrm{max}}=3$.  Figure 2:  The two models
described by Eq. (44). Both lines correspond to the SIGWs induced by a log-
normal power spectrum with $\sigma_{*}=0.5$. The solid line is the full
numeric result, while the dashed line has a $k^{3}$ scaling in the infrared
region, $k<0.1k_{*}$.  Figure 3:  Taken from Fig. 3 in Yuan  _et al._ (2019b).
The relation between $\delta\chi^{2}$ and the peak mass of PBHs
($m_{\rm{pbh}}^{*}$ corresponds to $k_{*}$) generated by the log-normal power
spectrum given in Eq. (42) with $\sigma_{*}=0.5$. The amplitude of ${A}$ is
fixed by assuming PBHs represent $10^{-3}$ of DM. The $5\sigma$ dashed line
corresponds to $\delta\chi^{2}=28.74$.  Figure 4:  Taken from Fig. 1 in Yuan
_et al._ (2019a). The energy density parameter of SIGWs generated by a
$\delta$-spectrum up to third-order. The sensitivity curves for LISA, IPTA,
FAST and SKA are shown. The dashed lines denote the second-order SIGWs, and
the solid black lines include the third-order correction.  Figure 5:  The
expected SNR by LISA/IPTA/FAST/IPTA in detecting SIGWs generated by
monochromatic PBHs. The dotted line corresponds to $\mathrm{SNR}=5$. We assume
the observation time to be $T=4$yr for LISA and $T=30$yr for PTAs. The value
of $\Delta t$, $N$ and $\sigma$ for current PTAs can be found in Kuroda  _et
al._ (2015) (see Table. 5).  Figure 6:  Taken from Fig.2 in Chen  _et al._
(2019b). The constraints on $f_{\mathrm{pbh}}$ through the null detection of
SIGWs by NANOGrav 11-yr observations. The dotted line corresponds to
$10^{-6}$. The colored region are the results excluded by OGLE microlensing
Niikura  _et al._ (2019b), EROS microlensing Tisserand  _et al._ (2007) and
the SGWB from PBH binaries Chen and Huang (2019).  Figure 7:  An overview on
the current observational constraints on $f_{\mathrm{pbh}}$. Constraints are
taken from EGB Carr  _et al._ (2010), WDs Graham  _et al._ (2015), HSC Niikura
_et al._ (2019a), EROS Tisserand  _et al._ (2007), OGLE Niikura  _et al._
(2019b), SGWB from binary PBHs Wang  _et al._ (2018); Chen and Huang (2019),
UFD Brandt (2016) and SIGWs using NANOGrav 11-yr data set Chen  _et al._
(2019b)

## References

  * Sasaki  _et al._ (2016) M. Sasaki, T. Suyama, T. Tanaka, and S. Yokoyama, Phys. Rev. Lett. 117, 061101 (2016), [erratum: Phys. Rev. Lett.121,no.5,059901(2018)], arXiv:1603.08338 [astro-ph.CO] .
  * Chen and Huang (2018) Z.-C. Chen and Q.-G. Huang, Astrophys. J. 864, 61 (2018), arXiv:1801.10327 [astro-ph.CO] .
  * Raidal  _et al._ (2019) M. Raidal, C. Spethmann, V. Vaskonen, and H. Veermäe, JCAP 1902, 018 (2019), arXiv:1812.01930 [astro-ph.CO] .
  * De Luca  _et al._ (2020a) V. De Luca, G. Franciolini, P. Pani, and A. Riotto, JCAP 06, 044 (2020a), arXiv:2005.05641 [astro-ph.CO] .
  * Hall  _et al._ (2020) A. Hall, A. D. Gow, and C. T. Byrnes, Phys. Rev. D 102, 123524 (2020), arXiv:2008.13704 [astro-ph.CO] .
  * Bhagwat  _et al._ (2021) S. Bhagwat, V. De Luca, G. Franciolini, P. Pani, and A. Riotto, JCAP 01, 037 (2021), arXiv:2008.12320 [astro-ph.CO] .
  * Hütsi  _et al._ (2020) G. Hütsi, M. Raidal, V. Vaskonen, and H. Veermäe, (2020), arXiv:2012.02786 [astro-ph.CO] .
  * Abbott  _et al._ (2016) B. P. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. Lett. 116, 061102 (2016), arXiv:1602.03837 [gr-qc] .
  * Tisserand  _et al._ (2007) P. Tisserand _et al._ (EROS-2), Astron. Astrophys. 469, 387 (2007), arXiv:astro-ph/0607207 [astro-ph] .
  * Carr  _et al._ (2010) B. J. Carr, K. Kohri, Y. Sendouda, and J. Yokoyama, Phys. Rev. D81, 104019 (2010), arXiv:0912.5297 [astro-ph.CO] .
  * Barnacka  _et al._ (2012) A. Barnacka, J. F. Glicenstein, and R. Moderski, Phys. Rev. D86, 043001 (2012), arXiv:1204.2056 [astro-ph.CO] .
  * Griest  _et al._ (2013) K. Griest, A. M. Cieplak, and M. J. Lehner, Phys. Rev. Lett. 111, 181302 (2013).
  * Graham  _et al._ (2015) P. W. Graham, S. Rajendran, and J. Varela, Phys. Rev. D92, 063007 (2015), arXiv:1505.04444 [hep-ph] .
  * Brandt (2016) T. D. Brandt, Astrophys. J. 824, L31 (2016), arXiv:1605.03665 [astro-ph.GA] .
  * Chen  _et al._ (2016) L. Chen, Q.-G. Huang, and K. Wang, JCAP 1612, 044 (2016), arXiv:1608.02174 [astro-ph.CO] .
  * Wang  _et al._ (2018) S. Wang, Y.-F. Wang, Q.-G. Huang, and T. G. F. Li, Phys. Rev. Lett. 120, 191102 (2018), arXiv:1610.08725 [astro-ph.CO] .
  * Gaggero  _et al._ (2017) D. Gaggero, G. Bertone, F. Calore, R. M. T. Connors, M. Lovell, S. Markoff, and E. Storm, Phys. Rev. Lett. 118, 241101 (2017), arXiv:1612.00457 [astro-ph.HE] .
  * Ali-Haimoud and Kamionkowski (2017) Y. Ali-Haimoud and M. Kamionkowski, Phys. Rev. D95, 043534 (2017), arXiv:1612.05644 [astro-ph.CO] .
  * Aloni  _et al._ (2017) D. Aloni, K. Blum, and R. Flauger, JCAP 1705, 017 (2017), arXiv:1612.06811 [astro-ph.CO] .
  * Horowitz (2016) B. Horowitz, (2016), arXiv:1612.07264 [astro-ph.CO] .
  * Niikura  _et al._ (2019a) H. Niikura _et al._ , Nat. Astron. 3, 524 (2019a), arXiv:1701.02151 [astro-ph.CO] .
  * Zumalacarregui and Seljak (2018) M. Zumalacarregui and U. Seljak, Phys. Rev. Lett. 121, 141101 (2018), arXiv:1712.02240 [astro-ph.CO] .
  * Nakama  _et al._ (2018) T. Nakama, T. Suyama, K. Kohri, and N. Hiroshima, Phys. Rev. D 97, 023539 (2018), arXiv:1712.08820 [astro-ph.CO] .
  * Abbott  _et al._ (2018) B. P. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. Lett. 121, 231103 (2018), arXiv:1808.04771 [astro-ph.CO] .
  * Magee  _et al._ (2018) R. Magee, A.-S. Deutsch, P. McClincy, C. Hanna, C. Horst, D. Meacher, C. Messick, S. Shandera, and M. Wade, Phys. Rev. D98, 103024 (2018), arXiv:1808.04772 [astro-ph.IM] .
  * Chen  _et al._ (2019a) Z.-C. Chen, F. Huang, and Q.-G. Huang, Astrophys. J. 871, 97 (2019a), arXiv:1809.10360 [gr-qc] .
  * Niikura  _et al._ (2019b) H. Niikura, M. Takada, S. Yokoyama, T. Sumi, and S. Masaki, Phys. Rev. D99, 083503 (2019b), arXiv:1901.07120 [astro-ph.CO] .
  * Chen and Huang (2019) Z.-C. Chen and Q.-G. Huang, (2019), arXiv:1904.02396 [astro-ph.CO] .
  * Abbott  _et al._ (2019) B. P. Abbott _et al._ (LIGO Scientific, Virgo), (2019), arXiv:1904.08976 [astro-ph.CO] .
  * Wang  _et al._ (2019a) Y.-F. Wang, Q.-G. Huang, T. G. F. Li, and S. Liao, (2019a), arXiv:1910.07397 [astro-ph.CO] .
  * Carr and Hawking (1974) B. J. Carr and S. W. Hawking, Mon. Not. Roy. Astron. Soc. 168, 399 (1974).
  * Carr (1975) B. J. Carr, Astrophys. J. 201, 1 (1975).
  * Chen  _et al._ (2019b) Z.-C. Chen, C. Yuan, and Q.-G. Huang, (2019b), arXiv:1910.12239 [astro-ph.CO] .
  * Tomita (1967) K. Tomita, Progress of Theoretical Physics 37, 831 (1967).
  * Matarrese  _et al._ (1993) S. Matarrese, O. Pantano, and D. Saez, Phys. Rev. D47, 1311 (1993).
  * Matarrese  _et al._ (1994) S. Matarrese, O. Pantano, and D. Saez, Phys. Rev. Lett. 72, 320 (1994), arXiv:astro-ph/9310036 [astro-ph] .
  * Matarrese  _et al._ (1998) S. Matarrese, S. Mollerach, and M. Bruni, Phys. Rev. D58, 043504 (1998), arXiv:astro-ph/9707278 [astro-ph] .
  * Noh and Hwang (2004) H. Noh and J.-c. Hwang, Phys. Rev. D69, 104011 (2004).
  * Carbone and Matarrese (2005) C. Carbone and S. Matarrese, Phys. Rev. D71, 043508 (2005), arXiv:astro-ph/0407611 [astro-ph] .
  * Nakamura (2007) K. Nakamura, Prog. Theor. Phys. 117, 17 (2007), arXiv:gr-qc/0605108 [gr-qc] .
  * Saito and Yokoyama (2009) R. Saito and J. Yokoyama, Phys. Rev. Lett. 102, 161101 (2009), [Erratum: Phys. Rev. Lett.107,069901(2011)], arXiv:0812.4339 [astro-ph] .
  * Ananda  _et al._ (2007) K. N. Ananda, C. Clarkson, and D. Wands, Phys. Rev. D75, 123518 (2007), arXiv:gr-qc/0612013 [gr-qc] .
  * Baumann  _et al._ (2007) D. Baumann, P. J. Steinhardt, K. Takahashi, and K. Ichiki, Phys. Rev. D76, 084019 (2007), arXiv:hep-th/0703290 [hep-th] .
  * Arroja  _et al._ (2009) F. Arroja, H. Assadullahi, K. Koyama, and D. Wands, Phys. Rev. D 80, 123526 (2009), arXiv:0907.3618 [astro-ph.CO] .
  * Assadullahi and Wands (2010) H. Assadullahi and D. Wands, Phys. Rev. D81, 023527 (2010), arXiv:0907.4073 [astro-ph.CO] .
  * Bugaev and Klimai (2010a) E. V. Bugaev and P. A. Klimai, JETP Lett. 91, 1 (2010a), arXiv:0911.0611 [astro-ph.CO] .
  * Bugaev and Klimai (2010b) E. Bugaev and P. Klimai, Phys. Rev. D81, 023517 (2010b), arXiv:0908.0664 [astro-ph.CO] .
  * Saito and Yokoyama (2010) R. Saito and J. Yokoyama, Prog. Theor. Phys. 123, 867 (2010), [Erratum: Prog. Theor. Phys.126,351(2011)], arXiv:0912.5317 [astro-ph.CO] .
  * Bugaev and Klimai (2011) E. Bugaev and P. Klimai, Phys. Rev. D83, 083521 (2011), arXiv:1012.4697 [astro-ph.CO] .
  * Alabidi  _et al._ (2013) L. Alabidi, K. Kohri, M. Sasaki, and Y. Sendouda, JCAP 05, 033 (2013), arXiv:1303.4519 [astro-ph.CO] .
  * Nakama and Suyama (2016) T. Nakama and T. Suyama, Phys. Rev. D94, 043507 (2016), arXiv:1605.04482 [gr-qc] .
  * Nakama  _et al._ (2017) T. Nakama, J. Silk, and M. Kamionkowski, Phys. Rev. D95, 043511 (2017), arXiv:1612.06264 [astro-ph.CO] .
  * Inomata  _et al._ (2017a) K. Inomata, M. Kawasaki, K. Mukaida, Y. Tada, and T. T. Yanagida, Phys. Rev. D 95, 123510 (2017a), arXiv:1611.06130 [astro-ph.CO] .
  * Orlofsky  _et al._ (2017) N. Orlofsky, A. Pierce, and J. D. Wells, Phys. Rev. D 95, 063518 (2017), arXiv:1612.05279 [astro-ph.CO] .
  * Garcia-Bellido  _et al._ (2017) J. Garcia-Bellido, M. Peloso, and C. Unal, JCAP 09, 013 (2017), arXiv:1707.02441 [astro-ph.CO] .
  * Sasaki  _et al._ (2018) M. Sasaki, T. Suyama, T. Tanaka, and S. Yokoyama, Class. Quant. Grav. 35, 063001 (2018), arXiv:1801.05235 [astro-ph.CO] .
  * Espinosa  _et al._ (2018) J. R. Espinosa, D. Racco, and A. Riotto, JCAP 1809, 012 (2018), arXiv:1804.07732 [hep-ph] .
  * Kohri and Terada (2018a) K. Kohri and T. Terada, Phys. Rev. D97, 123532 (2018a), arXiv:1804.08577 [gr-qc] .
  * Cai  _et al._ (2019a) R.-g. Cai, S. Pi, and M. Sasaki, Phys. Rev. Lett. 122, 201101 (2019a), arXiv:1810.11000 [astro-ph.CO] .
  * Bartolo  _et al._ (2019a) N. Bartolo, V. De Luca, G. Franciolini, A. Lewis, M. Peloso, and A. Riotto, Phys. Rev. Lett. 122, 211301 (2019a), arXiv:1810.12218 [astro-ph.CO] .
  * Bartolo  _et al._ (2019b) N. Bartolo, V. De Luca, G. Franciolini, M. Peloso, D. Racco, and A. Riotto, Phys. Rev. D99, 103521 (2019b), arXiv:1810.12224 [astro-ph.CO] .
  * Unal (2019) C. Unal, Phys. Rev. D99, 041301 (2019), arXiv:1811.09151 [astro-ph.CO] .
  * Byrnes  _et al._ (2019) C. T. Byrnes, P. S. Cole, and S. P. Patil, JCAP 06, 028 (2019), arXiv:1811.11158 [astro-ph.CO] .
  * Inomata and Nakama (2019) K. Inomata and T. Nakama, Phys. Rev. D99, 043511 (2019), arXiv:1812.00674 [astro-ph.CO] .
  * Clesse  _et al._ (2018) S. Clesse, J. García-Bellido, and S. Orani, (2018), arXiv:1812.11011 [astro-ph.CO] .
  * Cai  _et al._ (2019b) R.-G. Cai, S. Pi, S.-J. Wang, and X.-Y. Yang, JCAP 1905, 013 (2019b), arXiv:1901.10152 [astro-ph.CO] .
  * Inomata  _et al._ (2019a) K. Inomata, K. Kohri, T. Nakama, and T. Terada, (2019a), arXiv:1904.12878 [astro-ph.CO] .
  * Inomata  _et al._ (2019b) K. Inomata, K. Kohri, T. Nakama, and T. Terada, Phys. Rev. D100, 043532 (2019b), arXiv:1904.12879 [astro-ph.CO] .
  * Cai  _et al._ (2019c) R.-G. Cai, S. Pi, S.-J. Wang, and X.-Y. Yang, (2019c), arXiv:1907.06372 [astro-ph.CO] .
  * Yuan  _et al._ (2019a) C. Yuan, Z.-C. Chen, and Q.-G. Huang, Phys. Rev. D100, 081301 (2019a), arXiv:1906.11549 [astro-ph.CO] .
  * Cai  _et al._ (2019d) R.-G. Cai, S. Pi, and M. Sasaki, (2019d), arXiv:1909.13728 [astro-ph.CO] .
  * Lu  _et al._ (2019) Y. Lu, Y. Gong, Z. Yi, and F. Zhang, JCAP 12, 031 (2019), arXiv:1907.11896 [gr-qc] .
  * Yuan  _et al._ (2019b) C. Yuan, Z.-C. Chen, and Q.-G. Huang, (2019b), arXiv:1910.09099 [astro-ph.CO] .
  * Tomikawa and Kobayashi (2019) K. Tomikawa and T. Kobayashi, (2019), arXiv:1910.01880 [gr-qc] .
  * De Luca  _et al._ (2019a) V. De Luca, G. Franciolini, A. Kehagias, and A. Riotto, (2019a), arXiv:1911.09689 [gr-qc] .
  * Yuan  _et al._ (2019c) C. Yuan, Z.-C. Chen, and Q.-G. Huang, (2019c), arXiv:1912.00885 [astro-ph.CO] .
  * Inomata  _et al._ (2020a) K. Inomata, K. Kohri, T. Nakama, and T. Terada, J. Phys. Conf. Ser. 1468, 012002 (2020a).
  * Inomata  _et al._ (2020b) K. Inomata, K. Kohri, T. Nakama, and T. Terada, J. Phys. Conf. Ser. 1468, 012001 (2020b).
  * Inomata  _et al._ (2020c) K. Inomata, M. Kawasaki, K. Mukaida, T. Terada, and T. T. Yanagida, Phys. Rev. D 101, 123533 (2020c), arXiv:2003.10455 [astro-ph.CO] .
  * Yuan and Huang (2020) C. Yuan and Q.-G. Huang, (2020), arXiv:2007.10686 [astro-ph.CO] .
  * Papanikolaou  _et al._ (2020) T. Papanikolaou, V. Vennin, and D. Langlois, (2020), arXiv:2010.11573 [astro-ph.CO] .
  * Zhang  _et al._ (2020a) F. Zhang, A. Ali, Y. Gong, J. Lin, and Y. Lu, (2020a), arXiv:2008.12961 [gr-qc] .
  * Kapadia  _et al._ (2020a) S. J. Kapadia, K. L. Pandey, T. Suyama, S. Kandhasamy, and P. Ajith, (2020a), arXiv:2009.05514 [gr-qc] .
  * Zhang  _et al._ (2020b) F. Zhang, Y. Gong, J. Lin, Y. Lu, and Z. Yi, (2020b), arXiv:2012.06960 [astro-ph.CO] .
  * Domènech  _et al._ (2020a) G. Domènech, C. Lin, and M. Sasaki, (2020a), arXiv:2012.08151 [gr-qc] .
  * Dalianis and Kouvaris (2020) I. Dalianis and C. Kouvaris, (2020), arXiv:2012.09255 [astro-ph.CO] .
  * Atal and Domènech (2021) V. Atal and G. Domènech, (2021), arXiv:2103.01056 [astro-ph.CO] .
  * Carr  _et al._ (2020) B. Carr, K. Kohri, Y. Sendouda, and J. Yokoyama, (2020), arXiv:2002.12778 [astro-ph.CO] .
  * Green and Kavanagh (2020) A. M. Green and B. J. Kavanagh, (2020), arXiv:2007.10722 [astro-ph.CO] .
  * Carr and Kuhnel (2020) B. Carr and F. Kuhnel, Ann. Rev. Nucl. Part. Sci. 70, 355 (2020), arXiv:2006.02838 [astro-ph.CO] .
  * Bardeen  _et al._ (1986) J. M. Bardeen, J. R. Bond, N. Kaiser, and A. S. Szalay, Astrophys. J. 304, 15 (1986).
  * Ivanov  _et al._ (1994) P. Ivanov, P. Naselsky, and I. Novikov, Phys. Rev. D50, 7173 (1994).
  * Garcia-Bellido  _et al._ (1996) J. Garcia-Bellido, A. D. Linde, and D. Wands, Phys. Rev. D54, 6040 (1996), arXiv:astro-ph/9605094 [astro-ph] .
  * Ivanov (1998) P. Ivanov, Phys. Rev. D57, 7145 (1998), arXiv:astro-ph/9708224 [astro-ph] .
  * Yokoyama (1997) J. Yokoyama, Astron. Astrophys. 318, 673 (1997), arXiv:astro-ph/9509027 [astro-ph] .
  * Kawasaki  _et al._ (2006) M. Kawasaki, T. Takayama, M. Yamaguchi, and J. Yokoyama, Phys. Rev. D74, 043525 (2006), arXiv:hep-ph/0605271 [hep-ph] .
  * Hertzberg and Yamada (2018) M. P. Hertzberg and M. Yamada, Phys. Rev. D 97, 083509 (2018), arXiv:1712.09750 [astro-ph.CO] .
  * Inomata  _et al._ (2018) K. Inomata, M. Kawasaki, K. Mukaida, and T. T. Yanagida, Phys. Rev. D 97, 043514 (2018), arXiv:1711.06129 [astro-ph.CO] .
  * Inomata  _et al._ (2017b) K. Inomata, M. Kawasaki, K. Mukaida, Y. Tada, and T. T. Yanagida, Phys. Rev. D 96, 043504 (2017b), arXiv:1701.02544 [astro-ph.CO] .
  * Kohri and Terada (2018b) K. Kohri and T. Terada, Class. Quant. Grav. 35, 235017 (2018b), arXiv:1802.06785 [astro-ph.CO] .
  * Young  _et al._ (2014) S. Young, C. T. Byrnes, and M. Sasaki, JCAP 1407, 045 (2014), arXiv:1405.7023 [gr-qc] .
  * Bullock and Primack (1997) J. S. Bullock and J. R. Primack, Phys. Rev. D55, 7423 (1997), arXiv:astro-ph/9611106 [astro-ph] .
  * Pina Avelino (2005) P. Pina Avelino, Phys. Rev. D 72, 124004 (2005), arXiv:astro-ph/0510052 .
  * Hidalgo (2007) J. C. Hidalgo, (2007), arXiv:0708.3875 [astro-ph] .
  * Klimai and Bugaev (2012) P. A. Klimai and E. V. Bugaev, in  _17th International Seminar on High Energy Physics_ , Vol. 2 (INR, Moscow, Russia, 2012) pp. 163–174, arXiv:1210.3262 [astro-ph.CO] .
  * Byrnes  _et al._ (2012) C. T. Byrnes, E. J. Copeland, and A. M. Green, Phys. Rev. D86, 043512 (2012), arXiv:1206.4188 [astro-ph.CO] .
  * Young and Byrnes (2013) S. Young and C. T. Byrnes, JCAP 08, 052 (2013), arXiv:1307.4995 [astro-ph.CO] .
  * Franciolini  _et al._ (2018) G. Franciolini, A. Kehagias, S. Matarrese, and A. Riotto, JCAP 03, 016 (2018), arXiv:1801.09415 [astro-ph.CO] .
  * Riccardi  _et al._ (2021) F. Riccardi, M. Taoso, and A. Urbano, (2021), arXiv:2102.04084 [astro-ph.CO] .
  * Harada  _et al._ (2015) T. Harada, C.-M. Yoo, T. Nakama, and Y. Koga, Phys. Rev. D 91, 084057 (2015), arXiv:1503.03934 [gr-qc] .
  * Yoo  _et al._ (2018) C.-M. Yoo, T. Harada, J. Garriga, and K. Kohri, PTEP 2018, 123E01 (2018), arXiv:1805.03946 [astro-ph.CO] .
  * Musco (2019) I. Musco, Phys. Rev. D 100, 123524 (2019), arXiv:1809.02127 [gr-qc] .
  * Kawasaki and Nakatsuka (2019) M. Kawasaki and H. Nakatsuka, Phys. Rev. D 99, 123501 (2019), arXiv:1903.02994 [astro-ph.CO] .
  * De Luca  _et al._ (2019b) V. De Luca, G. Franciolini, A. Kehagias, M. Peloso, A. Riotto, and C. Ünal, JCAP 07, 048 (2019b), arXiv:1904.00970 [astro-ph.CO] .
  * Young  _et al._ (2019) S. Young, I. Musco, and C. T. Byrnes, JCAP 11, 012 (2019), arXiv:1904.00984 [astro-ph.CO] .
  * Ando  _et al._ (2018a) K. Ando, K. Inomata, and M. Kawasaki, Phys. Rev. D 97, 103528 (2018a), arXiv:1802.06393 [astro-ph.CO] .
  * Young (2019) S. Young, Int. J. Mod. Phys. D 29, 2030002 (2019), arXiv:1905.01230 [astro-ph.CO] .
  * Tokeshi  _et al._ (2020) K. Tokeshi, K. Inomata, and J. Yokoyama, JCAP 12, 038 (2020), arXiv:2005.07153 [astro-ph.CO] .
  * Moradinezhad Dizgah  _et al._ (2019) A. Moradinezhad Dizgah, G. Franciolini, and A. Riotto, JCAP 11, 001 (2019), arXiv:1906.08978 [astro-ph.CO] .
  * De Luca  _et al._ (2020b) V. De Luca, G. Franciolini, and A. Riotto, Phys. Lett. B 807, 135550 (2020b), arXiv:2001.04371 [astro-ph.CO] .
  * Yoo  _et al._ (2021) C.-M. Yoo, T. Harada, S. Hirano, and K. Kohri, PTEP 2021, 013E02 (2021), arXiv:2008.02425 [astro-ph.CO] .
  * Polnarev and Musco (2007) A. G. Polnarev and I. Musco, Class. Quant. Grav. 24, 1405 (2007), arXiv:gr-qc/0605122 .
  * Jedamzik (1997) K. Jedamzik, Phys. Rev. D 55, 5871 (1997), arXiv:astro-ph/9605152 .
  * Byrnes  _et al._ (2018) C. T. Byrnes, M. Hindmarsh, S. Young, and M. R. S. Hawkins, JCAP 08, 041 (2018), arXiv:1801.06138 [astro-ph.CO] .
  * Kehagias  _et al._ (2019) A. Kehagias, I. Musco, and A. Riotto, JCAP 12, 029 (2019), arXiv:1906.07135 [astro-ph.CO] .
  * Germani and Musco (2019) C. Germani and I. Musco, Phys. Rev. Lett. 122, 141302 (2019), arXiv:1805.04087 [astro-ph.CO] .
  * Jedamzik and Niemeyer (1999) K. Jedamzik and J. C. Niemeyer, Phys. Rev. D 59, 124014 (1999), arXiv:astro-ph/9901293 .
  * Shibata and Sasaki (1999) M. Shibata and M. Sasaki, Phys. Rev. D 60, 084002 (1999), arXiv:gr-qc/9905064 .
  * Musco  _et al._ (2005) I. Musco, J. C. Miller, and L. Rezzolla, Class. Quant. Grav. 22, 1405 (2005), arXiv:gr-qc/0412063 .
  * Hawke and Stewart (2002) I. Hawke and J. M. Stewart, Class. Quant. Grav. 19, 3687 (2002).
  * Musco  _et al._ (2009) I. Musco, J. C. Miller, and A. G. Polnarev, Class. Quant. Grav. 26, 235001 (2009), arXiv:0811.1452 [gr-qc] .
  * Musco and Miller (2013) I. Musco and J. C. Miller, Class. Quant. Grav. 30, 145009 (2013), arXiv:1201.2379 [gr-qc] .
  * Nakama  _et al._ (2014) T. Nakama, T. Harada, A. Polnarev, and J. Yokoyama, JCAP 01, 037 (2014), arXiv:1310.3007 [gr-qc] .
  * Harada and Jhingan (2016) T. Harada and S. Jhingan, PTEP 2016, 093E04 (2016), arXiv:1512.08639 [gr-qc] .
  * Escrivà  _et al._ (2020a) A. Escrivà, C. Germani, and R. K. Sheth, Phys. Rev. D 101, 044022 (2020a), arXiv:1907.13311 [gr-qc] .
  * Harada  _et al._ (2013) T. Harada, C.-M. Yoo, and K. Kohri, Phys. Rev. D88, 084051 (2013), [Erratum: Phys. Rev.D89,no.2,029903(2014)], arXiv:1309.4201 [astro-ph.CO] .
  * Escrivà  _et al._ (2020b) A. Escrivà, C. Germani, and R. K. Sheth, (2020b), arXiv:2007.05564 [gr-qc] .
  * Musco  _et al._ (2020) I. Musco, V. De Luca, G. Franciolini, and A. Riotto, (2020), arXiv:2011.03014 [astro-ph.CO] .
  * Gong (2019) J.-O. Gong, (2019), arXiv:1909.12708 [gr-qc] .
  * Saga  _et al._ (2015) S. Saga, K. Ichiki, and N. Sugiyama, Phys. Rev. D 91, 024030 (2015), arXiv:1412.1081 [astro-ph.CO] .
  * Malik and Wands (1998) K. A. Malik and D. Wands, (1998), arXiv:gr-qc/9804046 [gr-qc] .
  * Bartolo  _et al._ (2020) N. Bartolo, D. Bertacca, V. De Luca, G. Franciolini, S. Matarrese, M. Peloso, A. Ricciardone, A. Riotto, and G. Tasinato, JCAP 02, 028 (2020), arXiv:1909.12619 [astro-ph.CO] .
  * Pi and Sasaki (2020) S. Pi and M. Sasaki, JCAP 09, 037 (2020), arXiv:2005.12306 [gr-qc] .
  * Kuroyanagi  _et al._ (2018) S. Kuroyanagi, T. Chiba, and T. Takahashi, JCAP 11, 038 (2018), arXiv:1807.00786 [astro-ph.CO] .
  * Kawasaki  _et al._ (2013) M. Kawasaki, N. Kitajima, and S. Yokoyama, JCAP 08, 042 (2013), arXiv:1305.4464 [astro-ph.CO] .
  * Choudhury and Mazumdar (2014) S. Choudhury and A. Mazumdar, Phys. Lett. B 733, 270 (2014), arXiv:1307.5119 [astro-ph.CO] .
  * Di and Gong (2018) H. Di and Y. Gong, JCAP 07, 007 (2018), arXiv:1707.09578 [astro-ph.CO] .
  * Ando  _et al._ (2018b) K. Ando, K. Inomata, M. Kawasaki, K. Mukaida, and T. T. Yanagida, Phys. Rev. D97, 123512 (2018b), arXiv:1711.08956 [astro-ph.CO] .
  * Ando  _et al._ (2018c) K. Ando, M. Kawasaki, and H. Nakatsuka, Phys. Rev. D 98, 083508 (2018c), arXiv:1805.07757 [astro-ph.CO] .
  * Xu  _et al._ (2020) W.-T. Xu, J. Liu, T.-J. Gao, and Z.-K. Guo, Phys. Rev. D 101, 023505 (2020), arXiv:1907.05213 [astro-ph.CO] .
  * Özsoy and Tasinato (2020) O. Özsoy and G. Tasinato, JCAP 04, 048 (2020), arXiv:1912.01061 [astro-ph.CO] .
  * Gao and Yang (2019) T.-J. Gao and X.-Y. Yang, Int. J. Mod. Phys. A 34, 1950213 (2019).
  * Lin  _et al._ (2020) J. Lin, Q. Gao, Y. Gong, Y. Lu, C. Zhang, and F. Zhang, Phys. Rev. D 101, 103515 (2020), arXiv:2001.05909 [gr-qc] .
  * Ballesteros  _et al._ (2020) G. Ballesteros, J. Rey, M. Taoso, and A. Urbano, JCAP 07, 025 (2020), arXiv:2001.08220 [astro-ph.CO] .
  * Liu  _et al._ (2020) J. Liu, Z.-K. Guo, and R.-G. Cai, Phys. Rev. D 101, 083535 (2020), arXiv:2003.02075 [astro-ph.CO] .
  * Braglia  _et al._ (2020a) M. Braglia, D. K. Hazra, F. Finelli, G. F. Smoot, L. Sriramkumar, and A. A. Starobinsky, JCAP 08, 001 (2020a), arXiv:2005.02895 [astro-ph.CO] .
  * Fu  _et al._ (2020) C. Fu, P. Wu, and H. Yu, Phys. Rev. D 102, 043527 (2020), arXiv:2006.03768 [astro-ph.CO] .
  * Dalianis and Kritos (2021) I. Dalianis and K. Kritos, Phys. Rev. D 103, 023505 (2021), arXiv:2007.07915 [astro-ph.CO] .
  * Yi  _et al._ (2020a) Z. Yi, Y. Gong, B. Wang, and Z.-h. Zhu, (2020a), arXiv:2007.09957 [gr-qc] .
  * Dalianis (2020) I. Dalianis, PoS CORFU2019, 096 (2020).
  * Aldabergenov  _et al._ (2021) Y. Aldabergenov, A. Addazi, and S. V. Ketov, Phys. Lett. B 814, 136069 (2021), arXiv:2008.10476 [hep-th] .
  * Ragavendra  _et al._ (2020a) H. V. Ragavendra, P. Saha, L. Sriramkumar, and J. Silk, (2020a), arXiv:2008.12202 [astro-ph.CO] .
  * Bhaumik and Jain (2020) N. Bhaumik and R. K. Jain, (2020), arXiv:2009.10424 [astro-ph.CO] .
  * Zhou  _et al._ (2020) Z. Zhou, J. Jiang, Y.-F. Cai, M. Sasaki, and S. Pi, Phys. Rev. D 102, 103527 (2020), arXiv:2010.03537 [astro-ph.CO] .
  * Ragavendra  _et al._ (2020b) H. V. Ragavendra, L. Sriramkumar, and J. Silk, (2020b), arXiv:2011.09938 [astro-ph.CO] .
  * Yi  _et al._ (2020b) Z. Yi, Q. Gao, Y. Gong, and Z.-h. Zhu, (2020b), arXiv:2011.10606 [astro-ph.CO] .
  * Braglia  _et al._ (2020b) M. Braglia, X. Chen, and D. K. Hazra, (2020b), arXiv:2012.05821 [astro-ph.CO] .
  * Gao  _et al._ (2020) Q. Gao, Y. Gong, and Z. Yi, (2020), arXiv:2012.03856 [gr-qc] .
  * Gao (2021) Q. Gao, (2021), arXiv:2102.07369 [gr-qc] .
  * Domènech (2019) G. Domènech, (2019), arXiv:1912.05583 [gr-qc] .
  * Domènech  _et al._ (2020b) G. Domènech, S. Pi, and M. Sasaki, JCAP 08, 017 (2020b), arXiv:2005.12314 [gr-qc] .
  * Abe  _et al._ (2020) K. T. Abe, Y. Tada, and I. Ueda, (2020), arXiv:2010.06193 [astro-ph.CO] .
  * Zhou  _et al._ (2021) J.-Z. Zhou, X. Zhang, Q.-H. Zhu, and Z. Chang, (2021), arXiv:2106.01641 [astro-ph.CO] .
  * Adshead  _et al._ (2021) P. Adshead, K. D. Lozanov, and Z. J. Weiner, (2021), arXiv:2105.01659 [astro-ph.CO] .
  * Noh and Hwang (2003) H. Noh and J.-c. Hwang, (2003), arXiv:astro-ph/0305123 [astro-ph] .
  * Hwang  _et al._ (2017) J.-C. Hwang, D. Jeong, and H. Noh, Astrophys. J. 842, 46 (2017), arXiv:1704.03500 [astro-ph.CO] .
  * Inomata and Terada (2019) K. Inomata and T. Terada, (2019), arXiv:1912.00785 [gr-qc] .
  * Lu  _et al._ (2020) Y. Lu, A. Ali, Y. Gong, J. Lin, and F. Zhang, (2020), arXiv:2006.03450 [gr-qc] .
  * Ali  _et al._ (2020) A. Ali, Y. Gong, and Y. Lu, (2020), arXiv:2009.11081 [gr-qc] .
  * Inomata (2020) K. Inomata, (2020), arXiv:2008.12300 [gr-qc] .
  * Chang  _et al._ (2020a) Z. Chang, S. Wang, and Q.-H. Zhu, (2020a), arXiv:2009.11025 [astro-ph.CO] .
  * Chang  _et al._ (2020b) Z. Chang, S. Wang, and Q.-H. Zhu, (2020b), arXiv:2009.11994 [gr-qc] .
  * Chang  _et al._ (2020c) Z. Chang, S. Wang, and Q.-H. Zhu, (2020c), arXiv:2010.01487 [gr-qc] .
  * Domènech and Sasaki (2020) G. Domènech and M. Sasaki, (2020), arXiv:2012.14016 [gr-qc] .
  * Allen and Romano (1999) B. Allen and J. D. Romano, Phys. Rev. D 59, 102001 (1999), arXiv:gr-qc/9710117 .
  * Thrane and Romano (2013) E. Thrane and J. D. Romano, Phys. Rev. D88, 124032 (2013), arXiv:1310.5300 [astro-ph.IM] .
  * Larson  _et al._ (2000) S. L. Larson, W. A. Hiscock, and R. W. Hellings, Phys. Rev. D 62, 062001 (2000), arXiv:gr-qc/9909080 .
  * Robson  _et al._ (2019) T. Robson, N. J. Cornish, and C. Liu, Class. Quant. Grav. 36, 105011 (2019), arXiv:1803.01944 [astro-ph.HE] .
  * Siemens  _et al._ (2013) X. Siemens, J. Ellis, F. Jenet, and J. D. Romano, Class. Quant. Grav. 30, 224015 (2013), arXiv:1305.3196 [astro-ph.IM] .
  * Kuroda  _et al._ (2015) K. Kuroda, W.-T. Ni, and W.-P. Pan, Int. J. Mod. Phys. D24, 1530031 (2015), arXiv:1511.00231 [gr-qc] .
  * Wang  _et al._ (2019b) S. Wang, T. Terada, and K. Kohri, Phys. Rev. D99, 103531 (2019b), arXiv:1903.05924 [astro-ph.CO] .
  * De Luca  _et al._ (2021) V. De Luca, G. Franciolini, and A. Riotto, Phys. Rev. Lett. 126, 041303 (2021), arXiv:2009.08268 [astro-ph.CO] .
  * Sugiyama  _et al._ (2021) S. Sugiyama, V. Takhistov, E. Vitagliano, A. Kusenko, M. Sasaki, and M. Takada, Phys. Lett. B 814, 136097 (2021), arXiv:2010.02189 [astro-ph.CO] .
  * Vaskonen and Veermäe (2021) V. Vaskonen and H. Veermäe, Phys. Rev. Lett. 126, 051303 (2021), arXiv:2009.07832 [astro-ph.CO] .
  * Kohri and Terada (2021) K. Kohri and T. Terada, Phys. Lett. B 813, 136040 (2021), arXiv:2009.11853 [astro-ph.CO] .
  * Domènech and Pi (2020) G. Domènech and S. Pi, (2020), arXiv:2010.03976 [astro-ph.CO] .
  * Inomata  _et al._ (2020d) K. Inomata, M. Kawasaki, K. Mukaida, and T. T. Yanagida, (2020d), arXiv:2011.01270 [astro-ph.CO] .
  * Chen  _et al._ (2021) Z.-C. Chen, C. Yuan, and Q.-G. Huang, (2021), arXiv:2101.06869 [astro-ph.CO] .
  * Kapadia  _et al._ (2020b) S. J. Kapadia, K. L. Pandey, T. Suyama, and P. Ajith, Phys. Rev. D 101, 123535 (2020b), arXiv:2005.05693 [astro-ph.CO] .
  * Montero-Camacho  _et al._ (2019) P. Montero-Camacho, X. Fang, G. Vasquez, M. Silva, and C. M. Hirata, (2019), arXiv:1906.05950 [astro-ph.CO] .
  * Cotner  _et al._ (2019) E. Cotner, A. Kusenko, M. Sasaki, and V. Takhistov, JCAP 10, 077 (2019), arXiv:1907.10613 [astro-ph.CO] .
  * Pitrou  _et al._ (2013) C. Pitrou, X. Roy, and O. Umeh, Class. Quant. Grav. 30, 165002 (2013), arXiv:1302.6174 [astro-ph.CO] .