{'fraction_non_alphanumeric': 0.04223031387945239, 'fraction_numerical': 0.019853943326756038, 'mean_word_length': 4.236472092032382, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A topic review on probing primordial black hole dark matter with scalar induced gravitational waves

Summary 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 focuses on the physics about SIGWs accompanying the formation of PBH DM.


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 ∼ 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) 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, β, 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  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. (1993Matarrese et al. ( , 1994Matarrese et al. ( , 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) 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 ; 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.

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, β. 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) where g 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 * pbh ≈ 2 × 10 5 M t 1s . (2) The fraction of PBHs in all the dark matter, f pbh ≡ Ω PBH /Ω DM can be estimated by Nakama et al. (2017) f pbh 2.5 × 10 8 β g form *

10.75
On comoving slices, the relation between the primordial comoving curvature perturbation, ζ(k), and the density contrast, ∆(k), at linear order is given by where w is the equation of state and H is the Hubble parameter. The comoving curvature perturbation ζ is related to the metric perturbations by where our notations for the scalar perturbations ψ, v and B are introduced in Eq. (16) and Eq. (17) below. Moreover, the comoving curvature perturbation is related to the Bardeen potential (see Eq. (23) below) by For adiabatic perturbations, ζ stays constant on superhorizon scales. Then we have Φ = −3(1 + w)/(5 + 3w)ζ by assuming the absence of anisotropies and Φ = −2/3ζ during RD. The density contrast smoothed over a scale, R, is calculated as with W to be the window function chosen to smooth the density contrast. The variance of ∆(k, R) is given by where we define the dimensionless curvature power spectrum as ζ(k)ζ(k ) ≡ 2π 2 k 3 δ(k+k )P ζ (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 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. (2018Inomata et al. ( , 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 β. 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)) where the dimensionless threshold is defined as ν c ≡ ∆ c / ∆ 2 . k 2 is defined as Moreover, β is related to n pk by β = n pk (ν c )(2π) 3/2 R 3 . Another commonly used method to calculate β is the Press-Schechter formalism, where β is evaluated by simply integrating the probability density function (PDF) beyond the threshold value 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 ν 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.

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) Here ζ g is the Gaussian part whose PDF is Gaussian. The non-Gaussian parameter F NL would skew the PDF and G 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 β in the presence of non-Gaussianities. Let P NG to be the PDF of ζ and P G to be the PDF of the Gaussian part, ζ g . Then P NG can be obtained by changing the variables, ζ g = f −1 (ζ), such that This is equivalent to integrate the Gaussian PDF in the region, f −1 i (ζ) > ζ 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 β 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 β and non-Gaussian parameters (F NL , G NL , etc.).
Except for changing the variables in Eq. (13), Franciolini et al. (2018) adopted another method, the path-integral formulation, to calculate β in non-Gaussian regions. The authors expressed β 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 β. 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 β, 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).

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 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. 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 ζ need to be amplified by a factor of ∼ O(2) if using the linear relation, Eq. (4), to calculate the abundance.

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 where Θ 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 pbh . They also found that the mass function has a power-law tail, scaling as m −3/2 pbh . 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.

The formation criterion
Studies on the critical value of PBH formation can be traced back since 1975 when Carr first estimated the value of ∆ c by considering simplified Jeans length in Newtonian gravity, ∆ c = c 2 s Carr (1975), where c s = 1/ √ 3 is the sound speed during RD. However, the exact value of ∆ c is still a debating question since it depends on the initial density profile Polnarev and Musco (2007) (2020b). The analytical ∆ c in a general cosmological background is studied by Escrivà et al. (2020b) where they consider the equation of state to be w ∈ (0, 1]. Recently, Kehagias et al. (2019) estimated the effects of "intrinsic non-Gaussianities" generated by the non linear relation, Eq. (14), on ∆ c . They found that the relative change of ∆ 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 ∆ c for all possible shapes of the power spectrum, see Eq. (19) in Musco et al. (2020).
Apart from using ∆ 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 C = 2δM/R. Here R represents the areal radius and δ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, 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 C max instead of ∆ c . Recently, it was shown in Escrivà et al. (2020a) that the threshold value of compaction function, C max,c , is only sensitive to the curvature at the maximum and C max,c is, to some approximations, a universal quantity to describe the PBH formation.

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 where φ, B, ψ 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 where v is the velocity potential, P and ρ are the pressure and the energy density of the background while δP and δρ 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 and the scalar modes will transform as Malik and Wands (1998) where a prime denotes a derivative with respect to the conformal time η. 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 thatB =Ẽ = 0. This can be done by choosing L = −E and T = B − E . Then, the remaining two modes are just the Bardeen potential where σ ≡ E − 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 φ = Φ and ψ = Ψ. Write down the first-order Einstein equation in terms of the Bardeen potential, one finds that Φ = Ψ and the equation of motion for the only scalar mode left is Keeping the decay modes of Eq. (24), the solution of Φ in Fourier space is given by where we label the primordial value as Φ k and T Φ is the normalized transfer function such that T Φ (0) = 1. We also introduce the dimensionless variable x ≡ kη. The value of Φ 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 which is a source-free equation. The linear GWs will decay from the primordial value as η −1 during RD. During the formation of PBHs, the scalar power spectrum are enhanced to 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 (1)

ij . Readers interested in second-order GWs induced by h
(1) ij can refer to Gong (2019).
The source term in eq. (27) reads where we have expressed the density perturbation δρ, the velocity potential v and the pressure perturbation δP in terms of the scalar modes Φ. 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 We define S(η, k) ≡ −4e ij (k)S(η, k) and S(η, k) is the source term transferred to Fourier space. Here, we define the Fourier transform of h ij to be In the following part, we write h(η, k) to denote either the plus mode or the cross mode. In our convention, S(η, k) takes the form The transfer function is defined as For convenience, we introduce the dimensionless variable u ≡ p/k, v ≡ |p − k|/k and x ≡ kη. Unless otherwise being stated, the prime with T (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 (η; η ) = 1 k sin(kη − kη ) 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 ρ c (η) 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)) where in the last step of Eq. (34), we define the power spectrum as and the kernel function to be )/2 is the symmetric part of the transfer function. We also used that H = η −1 and a(η)/a(η) =η/η 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 "×" mode. In the last step of Eq. (34), the overline denotes the oscillating average, namely sin 2 x → 1 2 , cos 2 x → 1 2 and sin x cos x → 0. The Ω 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, ρ GW ∝ a −4 . Therefore, the current GW density parameter, Ω GW,0 would be where Ω r is the energy density fraction of radiation by today and we neglect the effects of relativistic degrees of freedom. Here, Ω GW (η c , f ) is the late time value when Ω GW (η, f ) becomes a constant and can be computed by taking I(u, v, x → ∞) (Here, a constant means that Ω GW (η, 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 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 Then the late time expression, I 2 (u, v, x → ∞), takes the form Below we use the comoving curvature perturbation ζ, which has the relation ζ = (3/2)Φ to calculate the SIGWs. A widely used model for the power spectrum is the infinite narrow spectrum 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 where we definek ≡ k/k * . Another example is a log-normal power spectrum parameterized by where σ * is a dimensionless parameter which denotes the width of the power spectrum. The SIGWs generated by this power spectrum has been carefully studied in . Finally, we consider a box spectrum described by In the above case, we have normalized the power spectrum so that P ζ (k) d ln k = A = ζ 2 is the variance of the perturbations. The Ω 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 k 1 and k 2 generate a mode with k = k 1 + 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 k * , they will induce GWs, and the shape of Ω 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 * / √ 3 where 1/ √ 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, Ω GW ∝ f 2 ln 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. They consider a fiducial case, Ω f id GW (k), generated by a log-normal power spectrum with σ * = 0.5 and the other one is described by These two cases have different scalings in the infrared region, as shown in Fig. 2. A statistic quantity, δχ 2 , which characterize the discrepancy of two models, takes the form Kuroyanagi et al. (2018) where T is the observation time and Ω n is the noise density parameter of the detector. Scanning the mass of PBHs, the result of δχ 2 is Fig. 3. It is shown that LISA can well distinguish the scaling of two models beyond 5σ in a wide mass range, especially for those PBHs of m pbh ∈ [10 −16 , 10 −14 ] ∪ [10 −13 , 10 −12 ]M which could still represent all the DM in our Universe. To summarize, the log-dependent scaling, 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

Scalar Induced Gravitational Waves in a General Cosmological Background
For a constant equation of state and let c 2 s = w, the semi-analytical solutions of SIGWs were obtained by Domènech in Domènech (2019) for 0 < w ≤ 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 2 s 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 2 s may become slightly smaller than 1/3 and the equation c 2 s = 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(η) and c s (η) differs from that in RD in the following aspects.
Firstly, the equation of motion for linear scalar perturbations become where the conformal Hubble parameter now become H(η) = 2 (1+3w)η . This equation should be numerically solved with the initial condition to be Φ Secondly, the equation of motion for the second-order tensor modes takes the same form as Eq. (27) but the source term changes to 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 The solution of g k (η,η) can be obtained by with u(η) and v(η) to be the two independent homogeneous solutions. Move on to the density parameter, the final expression for Ω GW takes the same form as the last line of Eq. (34) but the kernel function should now be evaluated as After calculating the oscillating average of I(u, v, x), the density parameter, Ω 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 Ω GW of SIGWs generated during the QCD phase transition is shown in Fig. 4 in Abe et al. (2020) for a delta-power spectrum.

Higher-Order Corrections to Scalar Induced Gravitational Waves
SIGWs are generated at second-order soured by linear scalar perturbations, which lead to h ij ∼ Φ 2 . Given that the scalar perturbations are enhanced to ∼ 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 δ-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) and the perturbed perfect fluid up to third-order in Newton gauge is given by (1) j (δP (1) + δρ (1) ) + (P + ρ) B (2) j + 4h (1) where v i ≡ v ⊥ i + ∂ i v and B i are the velocity perturbation and vector perturbation in Newton Gauge respectively.
The equation of motion for each h (n) ij takes the same form as Eq. (27) and the third-order source term reads and the fourth-order source term is given by Setting w = c 2 s = 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 then we can write the source terms as S (4) (k, η) = 4 d 3 pd 3 qd 3 l (2π) 9/2 e k (l, l)F +e k (p, q)F (4) mode at the final step. In Eq. (56), the time evolution of Φ is absorbed in the transfer function, namely where we set w = c 2 s = 1/3. The notation T k is short for the time evolution function, T Φ (kη). Applying the Green's function method, the solution to h (n) ij takes the same form as Eq. (29). The next-order corrections to Ω GW comes from two parts. The first part is the coupling of S while the second part is the coupling of S (2) ij and S The next-order corrections of Ω GW is the sum of Eq. (61) and Eq. (62). For a δ-spectrum, these two expressions can be further simplified after integrate over the δ function (see Yuan et al. (2019a) for results). After considering the higher-order corrections, the deep valley at k = 2/3k * 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.

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 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 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 Φ, namely Φ p Φ k−p Φ q Φ k −q . If Φ 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 where we assume that Φ = 0 and Φ p Φ k−p Φ q Φ k −q c is the connected four-point correlation function (4PCF) (or called the fourth cumulant in statistics) which is related to the trispectrum, T Φ , by The connected 4PCF will vanish if Φ 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 NL -order and calculate the SIGWs in the absence of the connected 4PCF. After that, the G 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 Ω 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).

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 φ is the just Bardeen potential Φ 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 Ω GW in uniform expansion gauge, comoving gauge and uniform curvature gauge and they found that Ω 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 η 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 ≥ 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 Ω GW , it is natural to ask which gauge is relevant in interpreting the observations. This question is discussed by De . 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, δg 00 = δg 0i = 0, corresponding toφ =B = 0 in Eq. (19). This leads to where C 1 and C 2 are two arbitrary spatial functions. The presence of C 1 and 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 C 1 and 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 After some algebra, one can get the equation of motion for E, (see also ) where the transfer function for E and ψ are defined as In our convention, the boundary condition for the transfer function should be T ψ (x → 0) = 3/2 in order to match the relation between curvature perturbation and the Bardeen potential. The general solution for Eq. (69) is where z ≡ x/ √ 3 and C i , (i = 3, 4, 5, 6) are integrate constants. One can check that under the residual gauge transformation, the gauge conditionφ =B = 0 still satisfy. Combining Eq. (72) and Eq. (74), we see that C 8 can be absorbed by C 3 while C 7 corresponds to C 5 . In other words, the constant C 3 and the logarithm term are two independent pure gauge modes. One can arbitrarily choose C 3 and 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 ψ (x → 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 Ci(x) = γ + ln(x) + O(x 2 ), where γ is the Euler gamma constant. Therefore, one must choose C 5 = −C 4 so that the logarithm term can be cancelled by the Ci(x) term. By applying the boundary condition T ψ (0) = 3/2, one gets C 6 = 0 and C 4 = −C 5 = 9 and the finite value at x = 0 is given by T E (0) = C 3 + 9γ. 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 Ω GW for SIGWs. On the other hand, if we let the gauge mode C 5 vanish, the presence of the 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  (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  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 ; 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 η 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 Ω GW is well defined in most of gauges if one takes the sub-horizon limit, namely η → ∞. The authors found that Ω 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 Φ 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.

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 σ = P ζ (k)d ln k = A. Using the Press-Schechter formalism, Eq. (11), the mass function can be approximated as which is exponentially dependent on A. On the other hand, the inevitably generated SIGWs satisfy h ij ∝ A. Therefore, a small change in A would change f 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), The strain power spectral density is given by S h (f ) = 3H 2 0 Ω GW,0 /(2π 2 f 3 ). Γ(f ) is the overlap function and P n (f ) is the noise power spectral density. For LISA, Γ(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) where f = c/(2πL) = 19.09 mHz is the transfer frequency and L = 2.5 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 oms L 2 + 2(1 + cos 2 (f /f )) The optical metrology noise spectrum and the acceleration noise is parameterized by Robson et al.
P oms (f ) = (1.5 × 10 −11 m) 2 1 + 2mHz f 4 Hz −1 (79) For pulsar timing arrays (PTA), the expressions are Γ(f ) = R(f ) = 1/(12π 2 f 2 ) and P n (f ) = 2∆tσ 2 where ∆t is the observation time and σ 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) where ζ IJ is the normalized Hellings and Downs coefficient for pulsars I and J. Assuming the pulsars are distributed homogeneously on the sky, then ζ IJ takes the form 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 10 −18 M 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 pbh . In , Wang et al. found the constraints on f pbh can reach 10 −13 for m pbh ∈ [10 −8 , 1]M .  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 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 pbh to be less than a few parts in million. Their work provides another independent test in addition to the γ-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 Although SIGWs place the currently severest constraint on f 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.

Summary and Outlook
PBHs can represent the DM in our Universe and explain the merger events detected by LIGO/VIRGO if f pbh ∼ 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 NL . However, up to G NL -order, this conclusion no longer holds. The G NL -order could further suppress the SIGWs than F 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 pbh . 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.   Table. 5).   (2019), UFD Brandt (2016) and SIGWs using NANOGrav 11-yr data set