# Primordial black holes and scalar-induced gravitational waves from the
polynomial attractor model

Zhongkai Wang<EMAIL_ADDRESS>School of Physics, Huazhong University of
Science and Technology, Wuhan, Hubei 430074, China  Shengqing Gao
<EMAIL_ADDRESS>School of Physics, Huazhong University of Science and
Technology, Wuhan, Hubei 430074, China  Yungui Gong  Corresponding author.
<EMAIL_ADDRESS>School of Physics, Huazhong University of Science and
Technology, Wuhan, Hubei 430074, China  Department of Physics, School of
Physical Science and Technology, Ningbo University, Ningbo, Zhejiang 315211,
China  Yue Wang<EMAIL_ADDRESS>School of Physics, Huazhong University
of Science and Technology, Wuhan, Hubei 430074, China

###### Abstract

Primordial black holes (PBHs) generated in the early universe are considered
as one of the candidates for dark matter. To produce PBHs with sufficient
abundance, the primordial scalar power spectrum needs to be enhanced to the
order of 0.01. Considering the third-order polynomial potential with
polynomial $\alpha$ attractors, we show that PBHs with the mass about
$10^{17}$g can be produced while satisfying the constraints from the cosmic
microwave background observations at the 2$\sigma$ confidence level. The mass
of PBHs produced in the polynomial $\alpha$ attractors can be much bigger than
that in the exponential $\alpha$ attractors. By adding a negative power-law
term to the polynomials, abundant PBHs with different masses and the accompany
scalar induced gravitational waves (SIGWs) with different peak frequency are
easily generated. The PBHs with masses around $10^{-15}-10^{-12}$ $M_{\odot}$
can account for almost all dark matter. The SIGWs generated in the nanohertz
band can explain the recent detection of stochastic gravitational wave
background by the pulsar timing array observations. The non-Gaussianity of the
primordial curvature perturbations in the squeezed and equilateral limits are
calculated numerically. We find that the non-Gaussianity correction greatly
enhances the PBH abundance which makes the production of PBHs much easier, but
the effect of non-Gaussianity on the generation of SIGWs is negligible.

##  I Introduction

The direct detection of gravitational waves (GWs) by the Laser Interferometer
Gravitational-Wave Observatory(LIGO) and Virgo collaborations opened a new
window to test gravity and understand the Universe [1, 2, 3, 4, 5, 6]. In
particular, some of the binary black holes observed may consist of primordial
black holes (PBHs) [7, 8, 9, 10, 11, 12], and the detections of stochastic GW
background by Pulsar Timing Array (PTA) also hint the existence of PBHs [13,
14, 15, 16, 17, 18, 19, 20, 21]. Furthermore, PBHs can be regarded as one of
the candidates for dark matter (DM) and their abundance and mass ranges are
tightly constrained by observations [22, 23, 24, 25, 26, 27].

If overdense inhomogeneities seeded from the primordial curvature
perturbations at horizon reentry exceed the threshold value $\delta_{c}$
during radiation domination, then PBHs could form in the overdense regions
through gravitational collapse [28, 29, 30]. Accompanied by the generation of
PBHs, the large curvature perturbation also produces scalar-induced GWs
(SIGWs) that contribute to the stochastic gravitational-wave background (SGWB)
[31, 32, 33, 34, 35, 36, 37, 38, 39, 40]. To generate significant abundance of
PBHs, the amplitude of the power spectrum $\mathcal{P}_{\zeta}$ of the
primordial curvature perturbation $\zeta$ should be
$A_{s}\sim\mathcal{O}(0.01)$ at small scales [41, 42, 43, 44, 45, 46, 47, 48,
49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63]. Compared with the
constraints on the amplitude of the power spectrum at the cosmic microwave
background (CMB) scales $A_{s}\approx 2.1\times 10^{-9}$ [64], the power
spectrum at small scales has to be enhanced by at least seven orders of
magnitude.

To amplify the power spectrum at small scales, a straightforward way is to
construct a flat plateau, i.e., an inflection point in the canonical inflation
potential [41, 65, 66, 43, 45, 44, 67, 68, 42, 69, 47, 48, 70, 50, 54, 71, 72,
73], leading to the so-called ultra-slow-roll (USR) [65, 66, 41] stage during
inflation. From the fact that $m$-order polynomial has $m-2$ inflection
points, the potential with $m=3$ can have one inflection point in the inflaton
potential. Thus, there are some attempts to combine the exponential $\alpha$
attractor models [74, 75, 76, 77, 78, 79, 80, 81, 82, 83] with the third-order
polynomial
$V(\psi)=V_{0}\left[1+c_{1}f(\psi)+c_{2}f^{2}(\psi)+c_{3}f^{3}(\psi)\right]^{2}$
to amplify the power spectrum and generate abundant PBHs [70, 84, 85, 86].
However, the masses of PBHs produced from these models are too small. Adding a
negative power-law term $f^{-2}(\psi)$ can solve this problem [86]. Apart from
the exponential $\alpha$ attractors, there are also polynomial $\alpha$
attractors [87, 88, 89, 90, 83]. In this paper, we replace the exponential
$\alpha$ attractors with the polynomial $\alpha$ attractors in the the third-
order polynomial potential with and without the additional negative power-law
term.

Near the inflection point, the slow-roll (SR) conditions are violated [67] and
there may exist a non-negligible non-Gaussianity. The abundance of PBHs and
the energy density of accompanying SIGWs can be affected by large non-
Gaussianities of the primordial curvature perturbations [91, 92, 93, 94, 95,
96, 97, 98, 99, 100]. The impacts of non-Gaussianity of the primordial
curvature perturbation $\zeta$ on the formation of PBHs were discussed in the
literature [101, 102, 103]. There are some methods to calculate the PBH
fractional energy density $\beta$ in the presence of non-Gaussianities, such
as changing the variables in the Gaussian probability density function [104]
and adopting the path-integral formulation [92, 91]. In this paper, we use the
latter method to calculate $\beta$.

This paper is organized as follows. In Sec. II, we show the models in detail
and present the power spectra of the models. In Sec. III, we discuss the PBH
abundance produced from the models. In Sec. IV, we compute the generation of
SIGWs during radiation domination. We also compare the SIGWs produced in the
models with the recent detection of nanohertz stochastic GW background by the
North American Nanohertz Observatory for Gravitational Waves (NANOGrav) [105]
and the European Pulsar Timing Array (EPTA) [106]. In Sec. V, we calculate the
non-Gaussianity of one model in the squeezed and equilateral limits and
discuss the effects of non-Gaussianity on the production of PBHs and SIGWs. We
conclude the paper in Sec. VI.

##  II The Models

For the enhancement of the primordial scalar power spectrum at small scales
with $\alpha$ attractors, third-order polynomial potentials with an inflection
point

$V(\psi)=V_{0}\left[1+c_{1}f(\psi)+c_{2}f^{2}(\psi)+c_{3}f^{3}(\psi)\right]^{2}$
(1)

were usually discussed [70, 84, 85, 86]. For T-models,
$f=f_{t}(\psi)=\tanh{(\psi/\sqrt{6\alpha})}$ [70, 84]; for E-models,
$f=f_{e}(\psi)=\exp{(-\sqrt{2/3\alpha}\,\psi)}$ [85], this is also belongs to
exponential $\alpha$ attractors [83]. The masses of PBHs produced by T-models
with the potential (1) are smaller than $10^{8}$ g [84], and the PBHs produced
by E-models with the potential (1) are asteroid-sized black holes [85]. To
obtain PBHs with larger masses, one must decrease the spectral tilt $n_{s}$
[86]. Because of the observational constraints by Planck, $n_{s}=0.9649\pm
0.0042$ (68% C.L.) [64], the spectral tilt $n_{s}$ cannot decrease too much.
To overcome this problem, a negative power-law term $f^{-2}(\psi)$ was added
to the potential for E-models with extreme fine-tuning of the parameters [86].
On the other hand, the polynomial potential with inverse powers

$V(\psi)=V_{0}\left(1-\frac{\mu^{k}}{\psi^{k}}+\cdots\right)$ (2)

gives the polynomial $\alpha$ attractors [83]

$n_{s}=1-\frac{2}{N}\frac{k+1}{k+2},\quad\mathit{r}=\frac{2\left(2\mathit{\mu}^{k}k\right)^{\frac{2}{k+2}}}{((k/2+1)\mathit{N})^{2-\frac{2}{k+2}}}.$
(3)

The factor $(k+1)/(k+2)$ in Eq. (3) helps to increase the value of $n_{s}$, so
it is possible to get larger mass PBHs with the polynomial $\alpha$ attractors
than T-models and E-models. In this paper, we use the polynomial $\alpha$
attractors [83] to discuss the enhancement of the primordial scalar power
spectrum at small scales from the polynomial potential (1) with the addition
of one negative power-law term,

$V(\psi)=V_{0}\left(1+c_{4}f_{n}^{-1}(\psi)+c_{1}f_{n}(\psi)+c_{2}f_{n}^{2}(\psi)+c_{3}f_{n}^{3}(\psi)\right)^{2},$
(4)

where the function $f_{n}(\psi)=\psi^{-n}$ and we consider three cases with
$n=1,2,3$.

For the potential (1), the inflection point is

$f(\psi^{\text{IP}})=-\frac{c_{2}}{3c_{3}},$ (5)

the derivative of the potential $V(f(\psi))$ at the inflection point is

$V^{\prime}(f(\psi^{\text{IP}}))=c_{1}-\frac{c_{2}^{2}}{3c_{3}}.$ (6)

To avoid a big e-folds $N$, the derivative of the potential at the inflection
point $\psi^{\text{IP}}$ should be close to zero, but not exactly zero, so we
take

$c_{1}=\frac{c_{2}^{2}}{3c_{3}}+\delta,\qquad\delta\ll 1.$ (7)

To compare the enhancement of the power spectrum between exponential $\alpha$
attractors and polynomial $\alpha$ attractors, we first discuss the difference
between the results for the potential (1) with $f=f_{e}$ and $f=f_{n}$, then
we will discuss the power spectrum of the polynomial $\alpha$ attractors (4).

###  II.1 Exponential and polynomial $\alpha$ attractors

There are three parameters in the polynomial potential (1). The condition (6)
for the inflection point reduces one parameter, so there are two free
parameters left. At large $\psi$, the potential (1) is dominated by the linear
term $c_{1}f(\psi)$ and the results depend on the parameter $c_{1}$ only
because of the attractors. The value of $V_{0}$ is determined by the amplitude
of the scalar power spectrum.

To show how to determine the model parameters, we take the model
$f(\psi)=f_{2}=\psi^{-2}$ as an example. Even though the attractor (3) is
independent of $c_{1}$, but the value of the scalar field $\psi_{*}$ at the
pivot scale $k_{*}=0.05$ Mpc-1 is
$\psi_{*}\approx(-16c_{1}N_{\text{sr}})^{1/4}$, so $\psi_{*}$ is determined by
$c_{1}$ and the e-folds $N_{\text{sr}}$ for the slow-roll period. Because of
the attractors (3) and the CMB constraints on $n_{s}$ and $r$, $N_{\text{sr}}$
cannot be too small. It cannot be too large either in order to enhance the
power spectrum at small scales. For the convenience of discussion on the
choices of parameters, we take $N_{\text{sr}}=40$. Note that we don’t fix
$N_{\text{sr}}$ in this paper. For fixed $N_{\text{sr}}$, $\psi_{*}$ is
determined by $c_{1}$. The value of $\psi_{*}$ is larger if $|c_{1}|$ is
larger. Once we take a value of $c_{1}$, then $c_{2}$ and $\psi^{\text{IP}}$
can be determined from $c_{1}$ and $c_{3}$. Combining Eqs (5) and (7), we get
$\psi^{\text{IP}}=(3c_{3}/c_{1})^{1/4}=2(-3c_{3}N_{\text{sr}})^{1/4}/\psi_{*}$,
so $\psi^{\text{IP}}$ is smaller if $|c_{1}|$ is larger. Now we fix
$c_{3}=-1.2$ and $N_{\text{sr}}=40$ to get $\psi_{*}$ by choosing
$c_{1}=-1.7$, $-1.9$, $-2.1$, respectively, then we fine-tune the value of
$\delta$, and numerically solve the background and perturbed cosmological
equations to obtain the scalar power spectrum $\mathcal{P}_{\zeta}$ of the
primordial curvature perturbation $\zeta$, the spectral tilt $n_{s}$ and the
tensor-to-scalar ratio $r$ at the pivot scale $k_{*}$. The results are shown
in Table 1. In Table 1, we also show the total number of e-folds $N$, the peak
scale $k^{\text{peak}}$ and the peak value $\mathcal{P}_{\zeta}^{\text{peak}}$
of the power spectrum at the peak scale $k^{\text{peak}}$. From Table 1, we
see that as $|c_{1}|$ becomes larger, $\psi_{*}$ and $r$ are bigger, but $N$
and $k^{\text{peak}}$ become smaller. For bigger $|c_{1}|$, the power spectrum
can be enhanced to the order $0.01$ with smaller number of e-folds, so the
scale at the end of inflation is also smaller, this leads to smaller peak
scale $k^{\text{peak}}$ even though $\psi^{\text{IP}}$ is smaller. Therefore,
to get massive PBHs, we should take bigger value of $|c_{1}|$. However, if
$|c_{1}|$ is too large, then inflation ends before the inflection point and
enhancement cannot happen. On the other hand, if $|c_{1}|$ is too small, then
inflaton stays at the inflection point for much longer time leading to a large
number of $N$. Since $N$ becomes smaller for larger $|c_{1}|$, we can increase
$N_{\text{sr}}$ to get bigger $n_{s}$, as shown in the last row of Table 1.
Even though we get bigger $n_{s}$ and smaller $r$, but the peak scale
$k^{\text{peak}}$ becomes larger and the mass of generated PBHs becomes
smaller, so we still use smaller $N_{\text{sr}}$. Now we fix $c_{1}=-2.1$ and
$\psi_{*}=6.055$, then vary the value of $c_{3}$, the results are shown in
Table 2. The results tell us that $N$, $n_{s}$ and $k^{\text{peak}}$ are
smaller with smaller $|c_{3}|$. Again, if we increase $\psi_{*}$ to get larger
$n_{s}$ for the same $c_{3}$, $k^{\text{peak}}$ becomes bigger as shown in the
last row of Table 2. The above discussion tells us that we need to choose
larger $|c_{1}|$ and smaller $|c_{3}|$ to get smaller $k^{\text{peak}}$.

$c_{1}$ | $\psi_{*}$ | $\delta$ | $N$ | $n_{s}$ | $r$ | $\mathcal{P}_{\zeta}^{\text{peak}}$ | $k^{\text{peak}}$
---|---|---|---|---|---|---|---
-1.7 | 5.743 | $4.51\times 10^{-4}$ | 66.2 | 0.9648 | 0.0091 | 0.011 | 5.9$\times 10^{18}$
-1.9 | 5.905 | $7.5882\times 10^{-3}$ | 50.5 | 0.9639 | 0.0099 | 0.01 | 1.2$\times 10^{18}$
-2.1 | 6.055 | $2.32285\times 10^{-2}$ | 46.7 | 0.9632 | 0.0107 | 0.012 | 4.4$\times 10^{17}$
-2.1 | 6.123 | $2.32285\times 10^{-2}$ | 48.6 | 0.9648 | 0.0100 | 0.011 | 3.1$\times 10^{18}$
Table 1: The results for the model $f_{2}$ with $c_{3}=-1.2$ and different
values of $c_{1}$. The unit of $k^{\text{peak}}$ is Mpc-1.

To make comparison for the models $f_{n}$ with differnt $n$, we fix
$n_{s}=0.9565$ which is consistent with CMB observations at the $2\sigma$
confidence level. Following the process discussed above, we can find the
appropriate parameters $c_{1}$ and $c_{3}$ for each model $f_{n}$, then we
calculate the power spectra and the results are shown in Table 3. Even though
we take the lower limit of the $2\sigma$ confidence level of the CMB results,
$k^{\text{peak}}$ can only be as small as $10^{14}\text{ Mpc}^{-1}\sim
10^{15}\text{ Mpc}^{-1}$ for the models $f_{n}$, and $k^{\text{peak}}$ is
about $10^{18}\text{ Mpc}^{-1}$ for the model $f_{e}$. Furthermore, for the
polynomial attractors $f_{n}$, $N$ is smaller than $40$. If we increase
$\psi_{*}$, $N$ and $n_{\text{s}}$ can be bigger, but then $k^{\text{peak}}$
becomes bigger too. On the other hand, the results show that $k^{\text{peak}}$
in the model $f_{n}$ can be smaller than that in the model $f_{e}$, so we only
discuss the polynomial attractors in the following sections.

$c_{3}$ | $\delta\ (10^{-2})$ | $N$ | $n_{s}$ | $r$ | $\mathcal{P}_{\zeta}^{\text{peak}}$ | $k^{\text{peak}}$
---|---|---|---|---|---|---
-1.2 | $2.32285$ | 46.7 | 0.9632 | 0.0107 | 0.012 | 4.4$\times 10^{17}$
-1.0 | $3.9718$ | 45.0 | 0.9628 | 0.0108 | 0.013 | 2.7$\times 10^{17}$
-0.8 | $6.48853$ | 43.5 | 0.9624 | 0.0109 | 0.011 | 1.3$\times 10^{17}$
-0.8 | $6.48853$ | 43.9 | 0.9628 | 0.0108 | 0.011 | 2.0$\times 10^{17}$
Table 2: The results of the model $f_{2}$ with $c_{1}=-2.1$ and different values of $c_{3}$. The unit of $k^{\text{peak}}$ is Mpc-1. Models | $V_{0}$ | $c_{1}$ | $c_{3}$ | $\delta$ | $\psi_{*}$ | $N$ | $r$ | $\mathcal{P}_{\zeta}^{\text{peak}}$ | $k^{\text{peak}}$
---|---|---|---|---|---|---|---|---|---
$f_{1}$ | 3.3$\times 10^{-9}$ | -2.9 | -1.16 | $2.84193\times 10^{-2}$ | 8.787 | 39.4 | 0.0503 | 0.011 | 2.8$\times 10^{14}$
$f_{2}$ | 5.5$\times 10^{-10}$ | -2.5 | -1.16 | $8.48823\times 10^{-2}$ | 6.113 | 37.8 | 0.0147 | 0.01 | 4.3$\times 10^{14}$
$f_{3}$ | 1.8$\times 10^{-10}$ | -2.2 | -1.17 | $1.284273\times 10^{-1}$ | 4.729 | 38.8 | 0.0053 | 0.014 | 1.6$\times 10^{15}$
$f_{e}$ | 1.25$\times 10^{-10}$ | -2.6 | -1.15 | $3.74916\times 10^{-2}$ | 5.461 | 51.1 | 0.0037 | 0.010 | 7.2$\times 10^{18}$
Table 3: The parameters and results of different models. The unit of $k^{\text{peak}}$ is Mpc-1. Models | $V_{0}$ | $\delta$ | $c_{1}$ | $c_{4}$ | $\psi_{*}$ | $N$ | $n_{s}$ | $r$ | $\mathcal{P}_{\zeta}^{\text{peak}}$ | $k^{\text{peak}}$
---|---|---|---|---|---|---|---|---|---|---
A1 | 8.33$\times 10^{-10}$ | 1.1007$\times 10^{-3}$ | -1.07 | 3.2$\times 10^{-4}$ | 3.75 | 63.2 | 0.9658 | 0.025 | 0.0338 | 3.0$\times 10^{8}$
A2 | 3.53$\times 10^{-10}$ | 6.3843$\times 10^{-4}$ | -1.1 | 1.3$\times 10^{-4}$ | 4.0 | 54.8 | 0.9652 | 0.011 | 0.029 | 1.7$\times 10^{12}$
A3 | 2.75$\times 10^{-10}$ | 7.574$\times 10^{-4}$ | -1.12 | 9.5$\times 10^{-5}$ | 4.08 | 52.4 | 0.9655 | 0.008 | 0.0269 | 3.6$\times 10^{13}$
Table 4: The parameters and results of modified polynomial $\alpha$
attractors. The unit of $k^{\text{peak}}$ is Mpc-1.

###  II.2 Modified polynomial $\alpha$ attractors

As discussed above, it is very difficult to produce bigger PBHs with the
third-order polynomial potentials (1), so we consider the modified polynomial
potential (4). From Table 3, we see that the tensor-to-scalar ratio $r$ in the
model $f_{1}$ is $0.05$ and $r=0.015$ in the model $f_{2}$, so we only
consider the polynomial attractors $f_{3}$. We choose three sets of parameters
for the model $f_{3}$ and denote them as A1, A2 and A3, respectively.

We take $c_{3}=-1.2$ for the models A1, A2 and A3, and the values of the other
parameters are respectively shown in Table 4. Note that the $c_{4}$ term will
affect the inflection point (5) and the relation (7), but its effection is
very small near the inflection point, so we still use Eq. (7) to determine
$c_{2}$ and absorb the effect of $c_{4}$ term into $\delta$. The results for
the number of e-folds $N$, the scalar spectral tilt $n_{\text{s}}$, the
tensor-to-scalar ratio $r$, the peak value of the power spectra
$\mathcal{P}^{\text{peak}}_{\zeta}$ and the peak scale $k^{\text{peak}}$ are
also shown in Table 4. The scalar power spectra are shown in Fig. 1. From Fig.
1, we see that the power spectrum for model A1 is broad and the power spectra
for models A2 and A3 are narrow. From Table 4, we see that the peak scale
$k^{\text{peak}}$ can be much smaller while satisfying the CMB constraints. In
the model A1, we even get $k^{\text{peak}}=3\times 10^{8}$ Mpc-1. As the peak
scale $k^{\text{peak}}$ becomes smaller, the tensor-to-scalar ratio $r$
becomes larger, and $r$ reaches $0.025$ in the model A1, so the peak scale
$k^{\text{peak}}$ cannot be too much smaller while satisfying the CMB
constraints.

Figure 1: The results of the power spectra for models A1, A2 and A3. The black
solid line is for model A1, the black dotted line is for model A2, and the
black dashed line is for model A3. The light green shaded region is excluded
by the CMB observations [64]. The light purple, cyan, and gray regions show
the constraints from the PTA observations [107], the effect on the ratio
between neutron and proton during the big bang nucleosynthesis (BBN) [108] and
$\mu$-distortion of CMB [109], respectively.

##  III PBH formation

The overdense region generated by large primordial curvature perturbations at
small scales may form PBHs through gravitational collapse after the horizon
reentry during radiation domination. The mass $M$ of generated PBHs is of the
same order as the horizon mass $M_{H}=(2GH)^{-1}$, $M=\gamma M_{H}$, here we
choose $\gamma=0.2$ [28]. Given a scale $k$, the PBH mass $M$ is [26, 58]

$M(k)=3.68\left(\frac{\gamma}{0.2}\right)\left(\frac{g_{*}}{10.75}\right)^{-1/6}\left(\frac{k}{10^{6}\text{
Mpc}^{-1}}\right)^{-2}M_{\odot},$ (8)

where $M_{\odot}$ is the solar mass, $g_{*}$ is the effective degrees of
freedom at the PBH formation, we don’t distinguish the effective degrees of
freedom in the entropy and energy density. The current fractional energy
density of PBHs with mass $M$ to DM is [41]

$f_{\text{pbh}}(M)=\frac{\beta(M)}{3.94\times
10^{-9}}\left(\frac{\gamma}{0.2}\right)^{1/2}\left(\frac{g_{*}}{10.75}\right)^{-1/4}\left(\frac{0.12}{\Omega_{\text{DM}}h^{2}}\right)\left(\frac{M}{M_{\odot}}\right)^{-1/2},$
(9)

where $\Omega_{\text{DM}}$ is the current energy density of DM, the Hubble
constant $H_{0}=100h$ km/s/Mpc, the PBH fractional energy density $\beta(M)$
at the formation for Gaussian curvature perturbation $\zeta$ is [110]

$\beta^{G}(M)\approx\sqrt{\frac{2}{\pi}}\frac{\sigma_{R}(M)}{\delta_{c}}\mathrm{exp}\left(-\frac{\delta_{c}^{2}}{2\sigma^{2}_{R}(M)}\right),$
(10)

$\delta_{c}$ is the critical density perturbation for the PBH formation, the
mass variance $\sigma_{R}(k)$ on the comoving scale $k=aH=R^{-1}$ is [110]

$\sigma^{2}_{R}(k)=\bigg{(}\frac{4}{9}\bigg{)}^{2}\int\frac{dq}{q}W^{2}(q/k)(q/k)^{4}\mathcal{P}_{\zeta}(q),$
(11)

and the Gaussian window function $W(x)=\exp(-x^{2}/2)$. In this paper, we
choose the parameters as $\delta_{c}=0.4$ [111, 112, 113, 114, 114, 115],
$g_{*}=106.75$ and $\Omega_{\text{DM}}h^{2}=0.12$. Using the power spectra
given in Fig. 1, we numerically compute the abundance of PBHs with the above
formulae and the results for the models A1, A2, and A3 are shown in Fig. 2.
The results show that all the produced PBHs satisfy the observational
constraints. The peak values of the mass and the abundance of PBH produced at
the peak scale $k^{\text{peak}}$ are shown in Table 5. As seen from Table 5,
in the models A1, A2 and A3, the PBH masses are around $10^{-5}M_{\odot}$,
$10^{-12}M_{\odot}$ and $10^{-15}M_{\odot}$, respectively. PBHs with masses
around $10^{-12}M_{\odot}$ and $10^{-15}M_{\odot}$ can account for most of DM.

Models | $M^{\text{peak}}/M_{\odot}$ | $f_{\text{pbh}}$ |  $f_{\text{c}}$/Hz
---|---|---|---
A1 | 2.51$\times 10^{-5}$ | 0.04 | 5.56$\times 10^{-7}$
A2 | 9.99$\times 10^{-13}$ | 0.765 | 3.13$\times 10^{-3}$
A3 | 2.29$\times 10^{-15}$ | 0.669 | 6.44$\times 10^{-2}$
Table 5: The results of the abundance and peak mass of PBHs, and the peak
frequency of SIGWs for the models A1, A2 and A3. Figure 2: The abundance of
PBHs produced by models A1, A2 and A3. The black solid line is for model A1,
the black dotted line is for model A2, and the black dashed line is for model
A3. The shaded regions are the observational constraints on PBH abundance: the
light purple region from extragalactic gamma-rays by PBH evaporation
(EG$\gamma$) [116], the green region from the galactic center 511 keV gamma-
ray line (INTEGRAL) [117], the yellow region from the white dwarf explosion
(WD) [118], the cyan region from microlensing events with Subaru HSC [119],
the deep blue region from the Kepler satellite [120], the gray region from the
EROS/MACHO [121], and the khaki region from accretion constraints by CMB
[122].

##  IV Scalar induced gravitational waves

Accompanied by the production of PBHs, the large primordial curvature
perturbations at small scales can source the tensor perturbation at the second
order after the horizon reentry to generate SIGWs. In this section, we
calculate the energy density of SIGWs generated in models A1, A2, and A3.

Figure 3: The SIGWs generated from models A1, A2 and A3. The black solid line
is for model A1, the black dotted line is for model A2, and the black dashed
line is for model A3. The green dashed curve denotes the SKA limit [123], The
purple dashed curve denotes the TianQin limit [124], the orange dotted curve
denotes the LISA limit [125], and the cyan dot-dashed curve denotes the Taiji
limit [126]. The blue and gray violins represent the NANOGrav 15-year and EPTA
DR2 data set, respectively [105, 106].

The second order tensor perturbation in Fourier space is

$h_{ij}(\bm{x},\eta)=\int\frac{d^{3}\bm{k}}{(2\pi)^{3/2}}e^{i\bm{k}\cdot\bm{x}}[h_{\bm{k}}(\eta)e_{ij}(\bm{k})+\bar{h}_{\bm{k}}(\eta)\bar{e}_{ij}(\bm{k})],$
(12)

where $e_{ij}(\bm{k})$ and $\bar{e}_{ij}(\bm{k})$ are the plus and cross
polarization tensors. The equation for the Fourier components of the tensor
perturbation with either polarization ($h_{\bm{k}}(\eta)$ or
$\bar{h}_{\bm{k}}(\eta)$) is

$h_{\bm{k}}^{\prime\prime}+2\mathcal{H}h_{\bm{k}}^{\prime}+k^{2}h_{\bm{k}}=4S_{\bm{k}},$
(13)

with the source term

$S_{\bm{k}}=\int\frac{d^{3}\tilde{\bm{k}}}{(2\pi)^{3/2}}e_{lm}(\bm{k})\tilde{k}^{l}\tilde{k}^{m}\left\\{2\Phi_{\tilde{\bm{k}}}\Phi_{\bm{k}-\tilde{\bm{k}}}+\frac{4}{3(1+w)}\left(\frac{1}{\mathcal{H}}\Phi^{\prime}_{\tilde{\bm{k}}}+\Phi_{\tilde{\bm{k}}}\right)\left(\frac{1}{\mathcal{H}}\Phi^{\prime}_{\bm{k}-\tilde{\bm{k}}}+\Phi_{\bm{k}-\tilde{\bm{k}}}\right)\right\\},$
(14)

where the equation of state $w=p/\rho=1/3$ and
$\mathcal{H}=a^{\prime}/a=1/\eta$ during the radiation domination, the prime
denotes the derivative with respect to the conformal time $\eta$, the Bardeen
potential $\Phi_{\bm{k}}$ can be expressed as the primordial value
$\phi_{\bm{k}}$ times the transfer function $T(k\eta)$

$\Phi_{\bm{k}}=T_{\phi}(k\eta)\phi_{\bm{k}},$ (15)

the transfer function during the radiation domination is

$T(x)=\frac{9}{x^{2}}\left(\frac{\mathrm{sin}\left(x/\sqrt{3}\right)}{x/\sqrt{3}}-\mathrm{cos}\left(x/\sqrt{3}\right)\right),$
(16)

and $x=k\eta$. The primordial value $\phi_{k}$ is related to the primordial
curvature perturbation as

$\left\langle\phi_{\bm{k}}\phi_{\bm{q}}\right\rangle=\frac{2\pi^{2}}{k^{3}}\left(\frac{3+3\omega}{5+3\omega}\right)^{2}\delta^{(3)}(\bm{k}+\bm{q})\mathcal{P}_{\zeta}(k).$
(17)

The power spectrum of the SIGWs is

$\left\langle
h_{\bm{k}}(\eta)h_{\bm{q}}(\eta)\right\rangle=\frac{2\pi^{2}}{k^{3}}\delta^{(3)}(\bm{k}+\bm{q})\mathcal{P}_{h}(k,\eta).$
(18)

Using the Green function method, the solution to Eq. (13) is

$\begin{split}h_{\bm{k}}(\eta)&=\frac{1}{a(\eta)}\int
d\tilde{\eta}G_{k}(\eta,\tilde{\eta})a(\tilde{\eta})S_{\bm{k}}(\tilde{\eta}),\\\
G_{k}(\eta,\tilde{\eta})&=\frac{1}{k}\mathrm{sin}\left[k(\eta-\tilde{\eta})\right].\end{split}$
(19)

Combining Eqs. (14)–(19), we obtain the power spectrum of SIGWs [34, 35, 52]

$\mathcal{P}_{h}(k,\eta)=4\int_{0}^{\infty}dv\int_{|1-v|}^{1+v}du\left[\frac{4v^{2}-(1-u^{2}+v^{2})^{2}}{4uv}\right]^{2}I^{2}_{\mathrm{RD}}(u,v,x\rightarrow\infty)\mathcal{P}_{\zeta}(kv)\mathcal{P}_{\zeta}(ku),$
(20)

where $u=|\bm{k}-\tilde{\bm{k}}|/k$, $v=|\tilde{\bm{k}}|/k$. The time average
of $I^{2}_{\text{RD}}(u,v,x\rightarrow\infty)$ is

$\begin{split}\overline{I^{2}_{RD}(u,v,x\rightarrow\infty)}=\frac{1}{2}\left(\frac{3(u^{2}+v^{2}-3)}{4u^{3}v^{3}x}\right)^{2}&\bigg{[}\left(-4uv+(u^{2}+v^{2}-3)\mathrm{log}\left|\frac{3-(u+v)^{2}}{3-(u-v)^{2}}\right|\right)^{2}\\\
&+\pi^{2}(u^{2}+v^{2}-3)^{2}\Theta(u+v-\sqrt{3})\bigg{]}.\end{split}$ (21)

With these results, we are ready to calculate the fractional energy density of
SIGWs [35]

$\Omega_{\mathrm{GW}}(k)=\frac{1}{24}\left(\frac{k}{\mathcal{H}}\right)^{2}\overline{\mathcal{P}_{h}(k,\eta)}.$
(22)

Since GWs behave like radiation, the current fractional energy density of
SIGWs is

$\Omega_{\text{GW}0}(k)=\Omega_{\text{GW}}(k)\frac{\Omega_{r0}}{\Omega_{r}(\eta)},$
(23)

where $\Omega_{r}(\eta)$ is the fraction energy density of radiation and
$\Omega_{r0}$ is its current value. Using the power spectra given in Fig. 1,
we calculate $\Omega_{\text{GW}0}$ for the models A1, A2 and A3 and the
results are shown in Fig. 3. The value of the peak frequency $f_{\text{c}}$
where $\Omega_{\text{GW}0}$ takes the maximum value is shown in Table 5. For
the model A1, even though the peak frequency of the SIGWs is around $10^{-7}$
Hz, the SIGWs have a broad spectrum and they can be detected in both the
nanohertz and millihertz frequency bands. The infrared parts of the SIGWs in
the model A1 is consistent with the stochastic GW background detected by
NANOGrav and EPTA. For the models A2 and A3, the spectra of the SIGWs are
relatively narrow. The peak frequency of SIGWs for model A2 is around
$10^{-3}$ Hz and the peak frequency of SIGWs for model A3 is around $10^{-2}$
Hz. The SIGWs for models A2 and A3 can be tested by future space-borne GW
detectors such as TianQin, Taiji and LISA.

##  V Primordial Non-Gaussianity

In this section, we calculate the non-Gaussianity of the primordial curvature
perturbation $\hat{\zeta}$ by taking the model A1 as an example. From the
definition of the bispectrum $B_{\zeta}$ [127, 128]

$\langle\hat{\zeta}_{k_{1}}\hat{\zeta}_{k_{2}}\hat{\zeta}_{k_{3}}\rangle=(2\pi)^{3}\delta^{(3)}(\bm{k}_{1}+\bm{k}_{2}+\bm{k}_{3})B_{\zeta}(k_{1},k_{2},k_{3}),$
(24)

we get the non-Gaussianity parameter $f_{\text{NL}}$ [127]

$f_{\text{NL}}(k_{1},k_{2},k_{3})=\frac{5}{6}\frac{B_{\zeta}(k_{1},k_{2},k_{3})}{P_{\zeta}(k_{1})P_{\zeta}(k_{2})+P_{\zeta}(k_{1})P_{\zeta}(k_{3})+P_{\zeta}(k_{2})P_{\zeta}(k_{3})},$
(25)

where $P_{\zeta}(k)=2\pi^{2}\mathcal{P}_{\zeta}(k)/k^{3}$, and the explicit
form of the bispectrum $B_{\zeta}(k_{1},k_{2},k_{3})$ can be found in [129,
97]. In the squeezed limit $k_{3}\to 0$, there is a consistency relationship
between the non-Gaussianity parameter $f_{\text{NL}}$ and the scalar spectral
tilt $n_{s}$ [130]

$\lim_{k_{3}\to 0}f_{\text{NL}}(k_{1},k_{2},k_{3})=\frac{5}{12}(1-n_{s}).$
(26)

Although the consistency relation was originally derived in the canonical
single-field inflation with slow-roll conditions, it was then proved to be
true for any inflation in which the inflaton is the only dynamical field
during inflation [131].

For the model A1, we numerically calculate $f_{\text{NL}}$ in the squeezed and
the equilateral limits and the results are shown in Fig. 4. As shown in the
upper panel of Fig. 4, the consistency relation holds in the squeezed limit.
In this model, the value of $f_{\text{NL}}$ is small when the power spectrum
reaches the peak and is pretty large at scales around $k\sim 10^{6}$ $\text{
Mpc}^{-1}$. The large values of $f_{\text{NL}}$ are due to the plunge of the
power spectrum by several orders of magnitude. The change of the scalar
spectral tilt near the peak is small, so in the squeezed limit $f_{\text{NL}}$
is almost constant when the power spectrum climbs up to the peak and rolls
down the peak.

Figure 4: The primordial scalar power spectrum $\mathcal{P}_{\zeta}$ and the
non-Gaussianity parameter $f_{\text{NL}}$ for model A1. In the upper panel,
the purple dashed line represents the power spectrum $\mathcal{P}_{\zeta}$,
the blue solid line represents $-\frac{12}{5}f_{\text{NL}}$ in the squeezed
limit with $k_{1}=k_{2}=10^{6}k_{3}=k$, and the orange dashed line represents
the scalar spectral tilt $n_{s}-1$. The lower panel shows
$\frac{12}{5}f_{\text{NL}}$ in the equilateral limit with
$k_{1}=k_{2}=k_{3}=k$.

Considering the non-Gaussianity of the primordial curvature perturbation
$\zeta$, the PBH fractional energy density becomes [92, 94, 95, 96]

$\beta=e^{\Delta_{3}}\beta^{G},$ (27)

where for almost monochromatic mass, the third cumulant $\Delta_{3}$ is
approximated as [98]

$\Delta_{3}\approx\frac{23\,\delta^{3}_{c}}{\mathcal{P}_{\zeta}(k^{\text{peak}})}f_{\text{NL}}(k^{\text{peak}},k^{\text{peak}},k^{\text{peak}}).$
(28)

Since
$23\,\delta_{c}^{3}f_{\text{NL}}(k^{\text{peak}},k^{\text{peak}},k^{\text{peak}})\sim\mathcal{O}(1)$
and $\mathcal{P}_{\zeta}(k^{\text{peak}})\sim\mathcal{O}(0.01)$, so the third
cumulant $\Delta_{3}\gg 1$ and the non-Gaussianity correction greatly enhances
the PBH abundance, this makes the production of PBHs become much easier.

In terms of the non-Gaussianity parameter $f_{\text{NL}}$, the primordial
curvature perturbation can be parameterized as

$\zeta(\textbf{x})=\zeta^{G}(\textbf{x})+\frac{3}{5}f_{\text{NL}}(\zeta^{G}(\textbf{x})^{2}-\langle\zeta^{G}(\textbf{x})^{2}\rangle),$
(29)

where $\zeta^{G}(\textbf{x})$ is the Gaussian part of the curvature
perturbation. The power spectrum of the curvature perturbation is

$\mathcal{P}_{\zeta}(k)=\mathcal{P}^{\text{G}}_{\zeta}(k)+\mathcal{P}^{\text{NG}}_{\zeta}(k),$
(30)

where the non-Gaussianity correction to the power spectrum is [93, 98]

$\mathcal{P}_{\zeta}^{\text{NG}}(k)=\bigg{(}\frac{3}{5}\bigg{)}^{2}\frac{k^{3}}{2\pi}f^{2}_{\text{NL}}\int
d^{3}\bm{p}\frac{\mathcal{P}_{\zeta}^{\text{G}}(p)}{p^{3}}\frac{\mathcal{P}^{\text{G}}_{\zeta}(|\bm{k}-\bm{p}|)}{|\bm{k}-\bm{p}|^{3}}.$
(31)

If the non-Gaussianity correction $\mathcal{P}_{\zeta}^{\text{NG}}(k)$ is
negligible, then the effect of non-Gaussianity on the generation of SIGWs can
be neglected. In our models, the non-Gaussianity parameter $f_{\text{NL}}$ at
the peak scale is small, so the contribution from the non-Gaussian curvature
perturbations to SIGWs is negligible.

##  VI Conclusion

Taking advantage of the existence of one inflection point in the third-order
polynomial potential and inflationary attractors, it was shown that the
primordial scalar power spectrum in the models with the combination of the
third-order polynomial and exponential $\alpha$ attractors can satisfy the CMB
constraints at large scales and get enhancement at small scales to produce
abundant PBHs and generate detectable SIGWs. For the exponential $\alpha$
attractors, the peak scale of the power spectrum $k^{\text{peak}}$ is about
$10^{18}\text{ Mpc}^{-1}$, and the PBHs produced from these models can reach
asteroid size only. To produce PBHs with bigger mass, we combine the third-
order polynomial and the polynomial attractor models to amplify the power
spectrum at small scales. For the polynomial attractors, $k^{\text{peak}}$ can
be as small as $10^{14}\text{ Mpc}^{-1}\sim 10^{15}\text{ Mpc}^{-1}$, and PBHs
with the mass about $10^{17}$g can be produced. However, the mass of PBHs
produced in these models is still too small. By adding a negative power-law
term to the polynomials, abundant PBHs with different masses and the accompany
SIGWs with different peak frequency are generated. For the three models
considered in this paper, the peak scales $k^{\text{peak}}$ are $3\times
10^{8}\text{ Mpc}^{-1}$, $1.7\times 10^{12}\text{ Mpc}^{-1}$, and $3.6\times
10^{13}\text{ Mpc}^{-1}$, respectively; the peak mass of PBHs are $2.51\times
10^{-5}M_{\odot}$, $9.99\times 10^{-13}M_{\odot}$ and $2.29\times
10^{-15}M_{\odot}$, respectively; the peak frequency of SIGWs are $5.56\times
10^{-7}$ Hz, $3.13\times 10^{-3}$ Hz and $6.44\times 10^{-2}$ Hz,
respectively. The PBHs with masses around $10^{-15}-10^{-12}$ $M_{\odot}$ can
account for almost all dark matter. The power spectrum of primordial curvature
perturbations for the model A1 is broad and the power spectra for the models
A2 and A3 are narrow, so the SIGWs produced have both broad and narrow
spectra. The infrared parts of the SIGWs in the model A1 can explain the
stochastic GW background detected by NANOGrav and EPTA. The non-Gaussianity of
the primordial curvature perturbations in the squeezed and equilateral limits
are calculated. We find that the non-Gaussianity correction greatly enhances
the PBH abundance which makes the production of PBHs much easier, but the
effect of non-Gaussianity on the generation of SIGWs is negligible.

apsrev4-2.bst 2019-01-14 (MD) hand-edited version of apsrev4-1.bst Control:
key (0) Control: author (72) initials jnrlst Control: editor formatted (1)
identically to author Control: production of article title (-1) disabled
Control: page (0) single Control: year (1) truncated Control: production of
eprint (0) enabled

## References

  * Abbott et al. [2016a] B. P. Abbott et al. (LIGO Scientific and Virgo Collaborations), Observation of Gravitational Waves from a Binary Black Hole Merger, Phys. Rev. Lett. 116, 061102 (2016a).
  * Abbott et al. [2016b] B. P. Abbott et al. (LIGO Scientific and Virgo Collaborations), GW150914: The Advanced LIGO Detectors in the Era of First Discoveries, Phys. Rev. Lett. 116, 131103 (2016b).
  * Abbott et al. [2019] B. P. Abbott et al. (LIGO Scientific and Virgo Collaborations), GWTC-1: A Gravitational-Wave Transient Catalog of Compact Binary Mergers Observed by LIGO and Virgo during the First and Second Observing Runs, Phys. Rev. X 9, 031040 (2019).
  * Abbott et al. [2021] R. Abbott et al. (LIGO Scientific and Virgo Collaborations), GWTC-2: Compact Binary Coalescences Observed by LIGO and Virgo During the First Half of the Third Observing Run, Phys. Rev. X 11, 021053 (2021).
  * Abbott et al. [2024] R. Abbott et al. (LIGO Scientific and Virgo Collaborations), GWTC-2.1: Deep extended catalog of compact binary coalescences observed by LIGO and Virgo during the first half of the third observing run, Phys. Rev. D 109, 022001 (2024).
  * Abbott et al. [2023] R. Abbott et al. (KAGRA, VIRGO and LIGO Scientific Collaborations), GWTC-3: Compact Binary Coalescences Observed by LIGO and Virgo during the Second Part of the Third Observing Run, Phys. Rev. X 13, 041039 (2023).
  * Bird et al. [2016] S. Bird, I. Cholis, J. B. Muñoz, Y. Ali-Haïmoud, M. Kamionkowski, E. D. Kovetz, A. Raccanelli, and A. G. Riess, Did LIGO detect dark matter?, Phys. Rev. Lett. 116, 201301 (2016).
  * Sasaki et al. [2016] M. Sasaki, T. Suyama, T. Tanaka, and S. Yokoyama, Primordial Black Hole Scenario for the Gravitational-Wave Event GW150914, Phys. Rev. Lett. 117, 061101 (2016), [Erratum: Phys.Rev.Lett. 121, 059901 (2018)].
  * Inomata et al. [2017] K. Inomata, M. Kawasaki, K. Mukaida, Y. Tada, and T. T. Yanagida, Inflationary primordial black holes for the LIGO gravitational wave events and pulsar timing array experiments, Phys. Rev. D 95, 123510 (2017).
  * De Luca et al. [2021a] V. De Luca, V. Desjacques, G. Franciolini, P. Pani, and A. Riotto, GW190521 Mass Gap Event and the Primordial Black Hole Scenario, Phys. Rev. Lett. 126, 051101 (2021a).
  * De Luca et al. [2021b] V. De Luca, G. Franciolini, P. Pani, and A. Riotto, Bayesian Evidence for Both Astrophysical and Primordial Black Holes: Mapping the GWTC-2 Catalog to Third-Generation Detectors, J. Cosmol. Astropart. Phys. 05 (2021) 003.
  * Franciolini et al. [2022] G. Franciolini, V. Baibhav, V. De Luca, K. K. Y. Ng, K. W. K. Wong, E. Berti, P. Pani, A. Riotto, and S. Vitale, Searching for a subpopulation of primordial black holes in LIGO-Virgo gravitational-wave data, Phys. Rev. D 105, 083526 (2022).
  * De Luca et al. [2021c] V. De Luca, G. Franciolini, and A. Riotto, NANOGrav Data Hints at Primordial Black Holes as Dark Matter, Phys. Rev. Lett. 126, 041303 (2021c).
  * Vaskonen and Veermäe [2021] V. Vaskonen and H. Veermäe, Did NANOGrav see a signal from primordial black hole formation?, Phys. Rev. Lett. 126, 051303 (2021).
  * Kohri and Terada [2021] K. Kohri and T. Terada, Solar-Mass Primordial Black Holes Explain NANOGrav Hint of Gravitational Waves, Phys. Lett. B 813, 136040 (2021).
  * Domènech and Pi [2022] G. Domènech and S. Pi, NANOGrav hints on planet-mass primordial black holes, Sci. China Phys. Mech. Astron. 65, 230411 (2022).
  * Atal et al. [2021] V. Atal, A. Sanglas, and N. Triantafyllou, NANOGrav signal as mergers of Stupendously Large Primordial Black Holes, J. Cosmol. Astropart. Phys. 06 (2021) 022.
  * Yi and Zhu [2022] Z. Yi and Z.-H. Zhu, NANOGrav signal and LIGO-Virgo primordial black holes from the Higgs field, J. Cosmol. Astropart. Phys. 05 (2022) 046.
  * Afzal et al. [2023] A. Afzal et al. (NANOGrav Collaboration), The NANOGrav 15 yr Data Set: Search for Signals from New Physics, Astrophys. J. Lett. 951, L11 (2023).
  * Antoniadis et al. [2023a] J. Antoniadis et al. (EPTA Collaboration), The second data release from the European Pulsar Timing Array: V. Implications for massive black holes, dark matter and the early Universe, arXiv:2306.16227 .
  * Yi et al. [2023] Z. Yi, Q. Gao, Y. Gong, Y. Wang, and F. Zhang, Scalar induced gravitational waves in light of Pulsar Timing Array data, Sci. China Phys. Mech. Astron. 66, 120404 (2023).
  * Ivanov et al. [1994] P. Ivanov, P. Naselsky, and I. Novikov, Inflation and primordial black holes as dark matter, Phys. Rev. D 50, 7173 (1994).
  * Khlopov et al. [2005] M. Y. Khlopov, S. G. Rubin, and A. S. Sakharov, Primordial structure of massive black hole clusters, Astropart. Phys. 23, 265 (2005).
  * Frampton et al. [2010] P. H. Frampton, M. Kawasaki, F. Takahashi, and T. T. Yanagida, Primordial Black Holes as All Dark Matter, J. Cosmol. Astropart. Phys. 04 (2010) 023.
  * Belotsky et al. [2014] K. M. Belotsky, A. D. Dmitriev, E. A. Esipova, V. A. Gani, A. V. Grobov, M. Y. Khlopov, A. A. Kirillov, S. G. Rubin, and I. V. Svadkovsky, Signatures of primordial black hole dark matter, Mod. Phys. Lett. A 29, 1440005 (2014).
  * Carr et al. [2016] B. Carr, F. Kuhnel, and M. Sandstad, Primordial Black Holes as Dark Matter, Phys. Rev. D 94, 083504 (2016).
  * Carr and Kuhnel [2020] B. Carr and F. Kuhnel, Primordial Black Holes as Dark Matter: Recent Developments, Annu. Rev. Nucl. Part. Sci. 70, 355 (2020).
  * Carr [1975] B. J. Carr, The Primordial black hole mass spectrum, Astrophys. J. 201, 1 (1975).
  * Hawking [1971] S. Hawking, Gravitationally collapsed objects of very low mass, Mon. Not. R. Astron. Soc. 152, 75 (1971).
  * Carr and Hawking [1974] B. J. Carr and S. W. Hawking, Black holes in the early Universe, Mon. Not. R. Astron. Soc. 168, 399 (1974).
  * Ananda et al. [2007] K. N. Ananda, C. Clarkson, and D. Wands, The Cosmological gravitational wave background from primordial density perturbations, Phys. Rev. D 75, 123518 (2007).
  * Baumann et al. [2007] D. Baumann, P. J. Steinhardt, K. Takahashi, and K. Ichiki, Gravitational Wave Spectrum Induced by Primordial Scalar Perturbations, Phys. Rev. D 76, 084019 (2007).
  * Saito and Yokoyama [2009] R. Saito and J. Yokoyama, Gravitational wave background as a probe of the primordial black hole abundance, Phys. Rev. Lett. 102, 161101 (2009), [Erratum: Phys.Rev.Lett. 107, 069901 (2011)].
  * Kohri and Terada [2018] K. Kohri and T. Terada, Semianalytic calculation of gravitational wave spectrum nonlinearly induced from primordial curvature perturbations, Phys. Rev. D 97, 123532 (2018).
  * Espinosa et al. [2018a] J. R. Espinosa, D. Racco, and A. Riotto, A Cosmological Signature of the SM Higgs Instability: Gravitational Waves, J. Cosmol. Astropart. Phys. 09 (2018) 012.
  * Cai et al. [2019a] R.-G. Cai, S. Pi, S.-J. Wang, and X.-Y. Yang, Pulsar Timing Array Constraints on the Induced Gravitational Waves, J. Cosmol. Astropart. Phys. 10 (2019) 059.
  * Inomata et al. [2019] K. Inomata, K. Kohri, T. Nakama, and T. Terada, Enhancement of Gravitational Waves Induced by Scalar Perturbations due to a Sudden Transition from an Early Matter Era to the Radiation Era, Phys. Rev. D 100, 043532 (2019), [Erratum: Phys.Rev.D 108, 049901 (2023)].
  * Domènech et al. [2020] G. Domènech, S. Pi, and M. Sasaki, Induced gravitational waves as a probe of thermal history of the universe, J. Cosmol. Astropart. Phys. 08 (2020) 017.
  * Braglia et al. [2021] M. Braglia, X. Chen, and D. K. Hazra, Probing Primordial Features with the Stochastic Gravitational Wave Background, J. Cosmol. Astropart. Phys. 03 (2021) 005.
  * Yi [2023] Z. Yi, Primordial black holes and scalar-induced gravitational waves from the generalized Brans-Dicke theory, J. Cosmol. Astropart. Phys. 03 (2023) 048.
  * Di and Gong [2018] H. Di and Y. Gong, Primordial black holes and second order gravitational waves from ultra-slow-roll inflation, J. Cosmol. Astropart. Phys. 07 (2018) 007.
  * Espinosa et al. [2018b] J. R. Espinosa, D. Racco, and A. Riotto, Cosmological Signature of the Standard Model Higgs Vacuum Instability: Primordial Black Holes as Dark Matter, Phys. Rev. Lett. 120, 121301 (2018b).
  * Garcia-Bellido and Ruiz Morales [2017] J. Garcia-Bellido and E. Ruiz Morales, Primordial black holes from single field models of inflation, Phys. Dark Univ. 18, 47 (2017).
  * Germani and Prokopec [2017] C. Germani and T. Prokopec, On primordial black holes from an inflection point, Phys. Dark Univ. 18, 6 (2017).
  * Ezquiaga et al. [2018] J. M. Ezquiaga, J. Garcia-Bellido, and E. Ruiz Morales, Primordial Black Hole production in Critical Higgs Inflation, Phys. Lett. B 776, 345 (2018).
  * Ballesteros and Taoso [2018] G. Ballesteros and M. Taoso, Primordial black hole dark matter from single field inflation, Phys. Rev. D 97, 023501 (2018).
  * Sasaki et al. [2018] M. Sasaki, T. Suyama, T. Tanaka, and S. Yokoyama, Primordial black holes—perspectives in gravitational wave astronomy, Classical Quantum Gravity 35, 063001 (2018).
  * Gao and Guo [2018] T.-J. Gao and Z.-K. Guo, Primordial Black Hole Production in Inflationary Models of Supergravity with a Single Chiral Superfield, Phys. Rev. D 98, 063526 (2018).
  * Dalianis et al. [2020] I. Dalianis, S. Karydas, and E. Papantonopoulos, Generalized Non-Minimal Derivative Coupling: Application to Inflation and Primordial Black Hole Production, J. Cosmol. Astropart. Phys. 06 (2020) 040.
  * Ballesteros et al. [2019] G. Ballesteros, J. Beltran Jimenez, and M. Pieroni, Black hole formation from a general quadratic action for inflationary primordial fluctuations, J. Cosmol. Astropart. Phys. 06 (2019) 016.
  * Sato-Polito et al. [2019] G. Sato-Polito, E. D. Kovetz, and M. Kamionkowski, Constraints on the primordial curvature power spectrum from primordial black holes, Phys. Rev. D 100, 063521 (2019).
  * Lu et al. [2019] Y. Lu, Y. Gong, Z. Yi, and F. Zhang, Constraints on primordial curvature perturbations from primordial black hole dark matter and secondary gravitational waves, J. Cosmol. Astropart. Phys. 12 (2019) 031.
  * Ballesteros et al. [2020] G. Ballesteros, J. Rey, M. Taoso, and A. Urbano, Primordial black holes as dark matter and gravitational waves from single-field polynomial inflation, J. Cosmol. Astropart. Phys. 07 (2020) 025.
  * Passaglia et al. [2019] S. Passaglia, W. Hu, and H. Motohashi, Primordial black holes and local non-Gaussianity in canonical inflation, Phys. Rev. D 99, 043536 (2019).
  * Fu et al. [2019] C. Fu, P. Wu, and H. Yu, Primordial Black Holes from Inflation with Nonminimal Derivative Coupling, Phys. Rev. D 100, 063532 (2019).
  * Lin et al. [2020] J. Lin, Q. Gao, Y. Gong, Y. Lu, C. Zhang, and F. Zhang, Primordial black holes and secondary gravitational waves from $k$ and $G$ inflation, Phys. Rev. D 101, 103515 (2020).
  * Yi et al. [2021a] Z. Yi, Y. Gong, B. Wang, and Z.-h. Zhu, Primordial black holes and secondary gravitational waves from the Higgs field, Phys. Rev. D 103, 063535 (2021a).
  * Yi et al. [2021b] Z. Yi, Q. Gao, Y. Gong, and Z.-h. Zhu, Primordial black holes and scalar-induced secondary gravitational waves from inflationary models with a noncanonical kinetic term, Phys. Rev. D 103, 063534 (2021b).
  * Ragavendra et al. [2021] H. V. Ragavendra, P. Saha, L. Sriramkumar, and J. Silk, Primordial black holes and secondary gravitational waves from ultraslow roll and punctuated inflation, Phys. Rev. D 103, 083510 (2021).
  * Kawaguchi and Tsujikawa [2023] R. Kawaguchi and S. Tsujikawa, Primordial black holes from Higgs inflation with a Gauss-Bonnet coupling, Phys. Rev. D 107, 063508 (2023).
  * Aldabergenov and Ketov [2023] Y. Aldabergenov and S. V. Ketov, Primordial Black Holes from Volkov–Akulov–Starobinsky Supergravity, Fortsch. Phys. 71, 2300039 (2023).
  * Cheung et al. [2023] K. Cheung, C. J. Ouseph, and P.-Y. Tseng, NANOGrav Signal and PBH from the Modified Higgs Inflation, arXiv:2307.08046 .
  * Hooper et al. [2023] D. Hooper, A. Ireland, G. Krnjaic, and A. Stebbins, Supermassive Primordial Black Holes From Inflation, arXiv:2308.00756 .
  * Akrami et al. [2020] Y. Akrami et al. (Planck Collaboration), Planck 2018 results. X. Constraints on inflation, Astron. Astrophys. 641, A10 (2020).
  * Martin et al. [2013] J. Martin, H. Motohashi, and T. Suyama, Ultra Slow-Roll Inflation and the non-Gaussianity Consistency Relation, Phys. Rev. D 87, 023514 (2013).
  * Motohashi et al. [2015] H. Motohashi, A. A. Starobinsky, and J. Yokoyama, Inflation with a constant rate of roll, J. Cosmol. Astropart. Phys. 09 (2015) 018.
  * Motohashi and Hu [2017] H. Motohashi and W. Hu, Primordial Black Holes and Slow-Roll Violation, Phys. Rev. D 96, 063503 (2017).
  * Bezrukov et al. [2018] F. Bezrukov, M. Pauly, and J. Rubio, On the robustness of the primordial power spectrum in renormalized Higgs inflation, J. Cosmol. Astropart. Phys. 02 (2018) 040.
  * Cicoli et al. [2018] M. Cicoli, V. A. Diaz, and F. G. Pedro, Primordial Black Holes from String Inflation, J. Cosmol. Astropart. Phys. 06 (2018) 034.
  * Dalianis et al. [2019] I. Dalianis, A. Kehagias, and G. Tringas, Primordial black holes from $\alpha$-attractors, J. Cosmol. Astropart. Phys. 01 (2019) 037.
  * Passaglia et al. [2020] S. Passaglia, W. Hu, and H. Motohashi, Primordial black holes as dark matter through Higgs field criticality, Phys. Rev. D 101, 123523 (2020).
  * Bhaumik and Jain [2020] N. Bhaumik and R. K. Jain, Primordial black holes dark matter from inflection point models of inflation and the effects of reheating, J. Cosmol. Astropart. Phys. 01 (2020) 037.
  * Xu et al. [2020] W.-T. Xu, J. Liu, T.-J. Gao, and Z.-K. Guo, Gravitational waves from double-inflection-point inflation, Phys. Rev. D 101, 023505 (2020).
  * Starobinsky [1980] A. A. Starobinsky, A New Type of Isotropic Cosmological Models Without Singularity, Phys. Lett. B 91, 99 (1980).
  * Salopek et al. [1989] D. S. Salopek, J. R. Bond, and J. M. Bardeen, Designing Density Fluctuation Spectra in Inflation, Phys. Rev. D 40, 1753 (1989).
  * Kallosh and Linde [2013a] R. Kallosh and A. Linde, Universality Class in Conformal Inflation, J. Cosmol. Astropart. Phys. 07 (2013) 002.
  * Kallosh and Linde [2013b] R. Kallosh and A. Linde, Superconformal generalizations of the Starobinsky model, J. Cosmol. Astropart. Phys. 06 (2013) 028.
  * Kallosh and Linde [2013c] R. Kallosh and A. Linde, Non-minimal Inflationary Attractors, J. Cosmol. Astropart. Phys. 10 (2013) 033.
  * Kallosh and Linde [2013d] R. Kallosh and A. Linde, Multi-field Conformal Cosmological Attractors, J. Cosmol. Astropart. Phys. 12 (2013) 006.
  * Kallosh et al. [2014] R. Kallosh, A. Linde, and D. Roest, Universal Attractor for Inflation at Strong Coupling, Phys. Rev. Lett. 112, 011303 (2014).
  * Kallosh et al. [2013] R. Kallosh, A. Linde, and D. Roest, Superconformal Inflationary $\alpha$-Attractors, J. High Energ. Phys. 11 (2013) 198.
  * Galante et al. [2015] M. Galante, R. Kallosh, A. Linde, and D. Roest, Unity of Cosmological Inflation Attractors, Phys. Rev. Lett. 114, 141302 (2015).
  * Kallosh and Linde [2022] R. Kallosh and A. Linde, Polynomial $\alpha$-attractors, J. Cosmol. Astropart. Phys. 04 (2022) 017.
  * Iacconi et al. [2022] L. Iacconi, H. Assadullahi, M. Fasiello, and D. Wands, Revisiting small-scale fluctuations in $\alpha$-attractor models of inflation, J. Cosmol. Astropart. Phys. 06 (2022) 007.
  * Frolovsky et al. [2022] D. Frolovsky, S. V. Ketov, and S. Saburov, E-models of inflation and primordial black holes, Front. in Phys. 10, 1005333 (2022).
  * Frolovsky and Ketov [2023] D. Frolovsky and S. V. Ketov, Production of Primordial Black Holes in Improved E-Models of Inflation, Universe 9, 294 (2023).
  * Dvali and Tye [1999] G. R. Dvali and S. H. H. Tye, Brane inflation, Phys. Lett. B 450, 72 (1999).
  * Burgess et al. [2001] C. P. Burgess, M. Majumdar, D. Nolte, F. Quevedo, G. Rajesh, and R.-J. Zhang, The Inflationary brane anti-brane universe, J. High Energ. Phys. 07 (2001) 047.
  * Lorenz et al. [2008] L. Lorenz, J. Martin, and C. Ringeval, Brane inflation and the WMAP data: A Bayesian analysis, J. Cosmol. Astropart. Phys. 04 (2008) 001.
  * Kallosh et al. [2019] R. Kallosh, A. Linde, and Y. Yamada, Planck 2018 and Brane Inflation Revisited, J. High Energ. Phys. 01 (2019) 008.
  * Saito et al. [2008] R. Saito, J. Yokoyama, and R. Nagata, Single-field inflation, anomalous enhancement of superhorizon fluctuations, and non-Gaussianity in primordial black hole formation, J. Cosmol. Astropart. Phys. 06 (2008) 024.
  * Franciolini et al. [2018] G. Franciolini, A. Kehagias, S. Matarrese, and A. Riotto, Primordial Black Holes from Inflation and non-Gaussianity, J. Cosmol. Astropart. Phys. 03 (2018) 016.
  * Cai et al. [2019b] R.-G. Cai, S. Pi, and M. Sasaki, Gravitational Waves Induced by non-Gaussian Scalar Perturbations, Phys. Rev. Lett. 122, 201101 (2019b).
  * Kehagias et al. [2019] A. Kehagias, I. Musco, and A. Riotto, Non-Gaussian Formation of Primordial Black Holes: Effects on the Threshold, J. Cosmol. Astropart. Phys. 12 (2019) 029.
  * Atal and Germani [2019] V. Atal and C. Germani, The role of non-gaussianities in Primordial Black Hole formation, Phys. Dark Univ. 24, 100275 (2019).
  * Riccardi et al. [2021] F. Riccardi, M. Taoso, and A. Urbano, Solving peak theory in the presence of local non-gaussianities, J. Cosmol. Astropart. Phys. 08 (2021) 060.
  * Zhang et al. [2021a] F. Zhang, Y. Gong, J. Lin, Y. Lu, and Z. Yi, Primordial non-Gaussianity from G-inflation, J. Cosmol. Astropart. Phys. 04 (2021) 045.
  * Zhang et al. [2021b] F. Zhang, J. Lin, and Y. Lu, Double-peaked inflation model: Scalar induced gravitational waves and primordial-black-hole suppression from primordial non-Gaussianity, Phys. Rev. D 104, 063515 (2021b), [Erratum: Phys.Rev.D 104, 129902 (2021)].
  * Gao et al. [2021] Q. Gao, Y. Gong, and Z. Yi, Primordial black holes and secondary gravitational waves from natural inflation, Nucl. Phys. B 969, 115480 (2021).
  * Gao [2021] Q. Gao, Primordial black holes and secondary gravitational waves from chaotic inflation, Sci. China Phys. Mech. Astron. 64, 280411 (2021).
  * Bullock and Primack [1997] J. S. Bullock and J. R. Primack, NonGaussian fluctuations and primordial black holes from inflation, Phys. Rev. D 55, 7423 (1997).
  * Ivanov [1998] P. Ivanov, Nonlinear metric perturbations and production of primordial black holes, Phys. Rev. D 57, 7145 (1998).
  * Yokoyama [1998] J. Yokoyama, Cosmological constraints on primordial black holes produced in the near critical gravitational collapse, Phys. Rev. D 58, 107502 (1998).
  * Byrnes et al. [2012] C. T. Byrnes, E. J. Copeland, and A. M. Green, Primordial black holes as a tool for constraining non-Gaussianity, Phys. Rev. D 86, 043512 (2012).
  * Agazie et al. [2023] G. Agazie et al. (NANOGrav Collaboration), The NANOGrav 15 yr Data Set: Detector Characterization and Noise Budget, Astrophys. J. Lett. 951, L10 (2023).
  * Antoniadis et al. [2023b] J. Antoniadis et al. (EPTA and InPTA Collaborations), The second data release from the European Pulsar Timing Array - III. Search for gravitational wave signals, Astron. Astrophys. 678, A50 (2023b).
  * Inomata and Nakama [2019] K. Inomata and T. Nakama, Gravitational waves induced by scalar perturbations as probes of the small-scale primordial spectrum, Phys. Rev. D 99, 043511 (2019).
  * Inomata et al. [2016] K. Inomata, M. Kawasaki, and Y. Tada, Revisiting constraints on small scale perturbations from big-bang nucleosynthesis, Phys. Rev. D 94, 043527 (2016).
  * Fixsen et al. [1996] D. J. Fixsen, E. S. Cheng, J. M. Gales, J. C. Mather, R. A. Shafer, and E. L. Wright, The Cosmic Microwave Background spectrum from the full COBE FIRAS data set, Astrophys. J. 473, 576 (1996).
  * Young et al. [2014] S. Young, C. T. Byrnes, and M. Sasaki, Calculating the mass fraction of primordial black holes, J. Cosmol. Astropart. Phys. 07 (2014) 045.
  * Musco and Miller [2013] I. Musco and J. C. Miller, Primordial black hole formation in the early universe: critical behaviour and self-similarity, Classical Quantum Gravity 30, 145009 (2013).
  * Harada et al. [2013] T. Harada, C.-M. Yoo, and K. Kohri, Threshold of primordial black hole formation, Phys. Rev. D 88, 084051 (2013), [Erratum: Phys.Rev.D 89, 029903 (2014)].
  * Tada and Yokoyama [2019] Y. Tada and S. Yokoyama, Primordial black hole tower: Dark matter, earth-mass, and LIGO black holes, Phys. Rev. D 100, 023537 (2019).
  * Escrivà et al. [2020] A. Escrivà, C. Germani, and R. K. Sheth, Universal threshold for primordial black hole formation, Phys. Rev. D 101, 044022 (2020).
  * Yoo et al. [2020] C.-M. Yoo, T. Harada, and H. Okawa, Threshold of Primordial Black Hole Formation in Nonspherical Collapse, Phys. Rev. D 102, 043526 (2020), [Erratum: Phys.Rev.D 107, 049901 (2023)].
  * Carr et al. [2010] B. J. Carr, K. Kohri, Y. Sendouda, and J. Yokoyama, New cosmological constraints on primordial black holes, Phys. Rev. D 81, 104019 (2010).
  * Dasgupta et al. [2020] B. Dasgupta, R. Laha, and A. Ray, Neutrino and positron constraints on spinning primordial black hole dark matter, Phys. Rev. Lett. 125, 101101 (2020).
  * Graham et al. [2015] P. W. Graham, S. Rajendran, and J. Varela, Dark Matter Triggers of Supernovae, Phys. Rev. D 92, 063007 (2015).
  * Niikura et al. [2019] H. Niikura et al., Microlensing constraints on primordial black holes with Subaru/HSC Andromeda observations, Nat. Astron. 3, 524 (2019).
  * Griest et al. [2013] K. Griest, A. M. Cieplak, and M. J. Lehner, New Limits on Primordial Black Hole Dark Matter from an Analysis of Kepler Source Microlensing Data, Phys. Rev. Lett. 111, 181302 (2013).
  * Tisserand et al. [2007] P. Tisserand et al. (EROS-2 Collaboration), Limits on the Macho Content of the Galactic Halo from the EROS-2 Survey of the Magellanic Clouds, Astron. Astrophys. 469, 387 (2007).
  * Poulin et al. [2017] V. Poulin, P. D. Serpico, F. Calore, S. Clesse, and K. Kohri, CMB bounds on disk-accreting massive primordial black holes, Phys. Rev. D 96, 083524 (2017).
  * Moore et al. [2015] C. J. Moore, R. H. Cole, and C. P. L. Berry, Gravitational-wave sensitivity curves, Classical Quantum Gravity 32, 015014 (2015).
  * Luo et al. [2016] J. Luo et al. (TianQin Collaboration), TianQin: a space-borne gravitational wave detector, Classical Quantum Gravity 33, 035010 (2016).
  * Amaro-Seoane et al. [2017] P. Amaro-Seoane et al. (LISA Collaboration), Laser Interferometer Space Antenna, arXiv:1702.00786 .
  * Hu and Wu [2017] W.-R. Hu and Y.-L. Wu, The Taiji Program in Space for gravitational wave physics and the nature of gravity, Natl. Sci. Rev. 4, 685 (2017).
  * Byrnes et al. [2010] C. T. Byrnes, M. Gerstenlauer, S. Nurmi, G. Tasinato, and D. Wands, Scale-dependent non-Gaussianity probes inflationary physics, J. Cosmol. Astropart. Phys. 10 (2010) 004.
  * Ade et al. [2016] P. A. R. Ade et al. (Planck Collaboration), Planck 2015 results. XVII. Constraints on primordial non-Gaussianity, Astron. Astrophys. 594, A17 (2016).
  * Hazra et al. [2013] D. K. Hazra, L. Sriramkumar, and J. Martin, BINGO: A code for the efficient computation of the scalar bi-spectrum, J. Cosmol. Astropart. Phys. 05 (2013) 026.
  * Maldacena [2003] J. M. Maldacena, Non-Gaussian features of primordial fluctuations in single field inflationary models, J. High Energ. Phys. 05 (2003) 013.
  * Creminelli and Zaldarriaga [2004] P. Creminelli and M. Zaldarriaga, Single field consistency relation for the 3-point function, J. Cosmol. Astropart. Phys. 10 (2004) 006.