# Non-Gaussianity and Secondary Gravitational Waves from Primordial Black
Holes Production in $\alpha$-attractor Inflation

Kazem<EMAIL_ADDRESS>Zeinab
<EMAIL_ADDRESS>Saeid<EMAIL_ADDRESS>and Kayoomars<EMAIL_ADDRESS>1School of Physics, Institute for
Research in Fundamental Sciences (IPM), P.O. Box 19395-5531, Tehran, Iran
2Department of Physics, University of Kurdistan, Pasdaran Street, P.O. Box
66177-15175, Sanandaj, Iran

###### Abstract

We study the non-Gaussianity and secondary Gravitational Waves (GWs) in the
process of the Primordial Black Holes (PBHs) production from inflation. In our
work, we focus on the $\alpha$-attractor inflation model in which a tiny bump
in the inflaton potential enhances the amplitude of the curvature
perturbations at some scales and consequently leads to the PBHs production
with different mass scales. We implement the computational code BINGO which
calculates the non-Gaussianity parameter in different triangle configurations.
Our examination implies that in this setup, the non-Gaussianity gets amplified
significantly in the equilateral shape around the scales in which the power
spectrum of the scalar perturbations undergoes a sharp declination. The
imprints of these non-Gaussianities can be probed in the scales corresponding
to the BBN and $\mu$-distortion events, or in smaller scales, and detection of
such signatures in the future observations may confirm the idea of our model
for the generation of PBHs or rule it out. Moreover, we investigate the
secondary GWs in this framework and show that in our model, the peak of the
present fractional energy density is obtained as $\Omega_{\rm GW0}\sim
10^{-8}$ at different frequencies which depends on the model parameters. These
results lie well within the sensitivity region of some GWs detectors at some
frequencies, and therefore the observational compatibility of our model can be
evaluated by the forthcoming data from these detectors. We further provide
some estimations for the tilts of the induced GWs spectrum in the different
intervals of frequency, and demonstrate that the spectrum obeys the power-law
relation $\Omega_{\rm GW0}\sim f^{n}$ in those frequency bands.

Inflation, Primordial black holes, Non-Gaussianity, Secondary gravitational
waves, Dark matter, Early universe

###### pacs:

98.80.Cq, 04.50.Kd, 05.70.Fh

##  I Introduction

The idea that a considerable part of the Dark Matter (DM) can be constituted
from Primordial Black Holes (PBHs) has absorbed much interests in the recent
years Bird:2016dcv ; Sasaki:2016jop ; Blinnikov:2016bxu ; Carr:2016drx ;
Clesse:2016vqa ; Carr:2017jsz ; mroz2017no ; Garcia-Bellido:2017mdw ;
Germani:2017bcs ; Clesse:2017bsw ; Ezquiaga:2017fvi ; Gong:2017qlj ;
ballesteros2018primordial ; Zumalacarregui:2017qqd ; Ballesteros:2018wlw ;
ashoorioon2021eft ; Ashoorioon:2020hln ; Cai:2018tuh ; Kamenshchik:2018sig ;
Niikura:2019kqi ; Chen:2019zza ; mishra2020primordial ; Ballesteros:2017fsr ;
Bhaumik:2019tvl ; germani2017primordial ; garcia2017primordial ;
Motohashi:2017kbs ; dalianis2019primordial ; Mahbub:2019uhl ; Kawai:2021edk ;
Motohashi:2019rhu ; Liu:2020oqe ; Fu:2019ttf ; dalianis2020generalized ;
Dalianis:2019vit ; Teimoori:2021thk ; Lin:2020goi ; Yi:2020cut ; Gao:2020tsa ;
Gao:2021vxb ; Solbi:2021rse ; Solbi:2021wbo ; Teimoori:2021pte ;
Heydari:2021gea ; Heydari:2021qsr . The first detection of Gravitational Wave
(GW), GW150914, from merging of binary black holes with mass $\sim
30M_{\odot}$ where $M_{\odot}$ is the solar mass, by the LIGO-Virgo
Collaboration Abbott:2016blz ; Abbott:2016nmj , has attracted more attention
to the physics of PBHs and increased the possibility that the origin of some
part of DM or even all of it may be PBHs seeded in the early universe.
Besides, the failure of direct detection of particle DM amplifies this idea.
PBHs can also be regarded as a DM candidate at other mass scales. Due to the
Hawking radiation Hawking:1975vcx , light PBHs are expected to be evaporated
by now, and this fact together with the implications of the recent
observations gould1992femtolensing ; dalcanton1994observational ;
alcock2001macho ; nemiroff2001limits ; wilkinson2001limits ;
tisserand2007limits ; carr2010new ; griest2013new ; Jacobs:2014yca ; Ali-
Haimoud:2016mbv ; Wang:2016ana ; Sato-Polito:2019hws ; Laha:2020ivk  put tight
constraints on the range of mass distribution of dark matter PBHs. PBHs with
earth mass ($\sim\mathcal{O}\left(10^{-5}\right)M_{\astrosun}$) can comprise
of order $\mathcal{O}\left(10^{-2}\right)$ of DM observed today. Moreover,
PBHs can constitute a substantial portion of DM in the two asteroid mass
ranges $10^{-16}-10^{-14}M_{\astrosun}$ and $10^{-13}-10^{-11}M_{\astrosun}$
Niikura:2017zjd ; Katz:2018zrn ; Laha:2019ssq ; Dasgupta:2019cae ; Carr:2009jm
; Graham:2015apa .

The origin of PBHs formation may return to the primordial perturbations
generated in the inflationary stage of the Universe. The relevant scales to
those perturbations exit the Hubble horizon at some epochs during inflation
and reenter the horizon later in the radiation-dominated era. The
gravitational collapse of those perturbations after their reentry gives rise
to PBHs formation if they have sufficiently large amplitude. This idea was
first proposed by Zel’dovich and Novikov in 1966, and then by Hawking and Carr
in the early 1970s zel1967hypothesis ; Hawking:1971ei ; Carr:1974nx . However,
providing such a large amplitude usually raises some difficulties, because CMB
observations constrain the amplitude of the scalar perturbations to be of
order ${\cal
P}_{s}(k_{*})=\mathcal{P}_{s*}\sim\mathcal{O}\left(10^{-9}\right)$ at the
pivot scale $k_{*}=0.05\,\mathrm{Mpc}^{-1}$ akrami2020planck , while in order
to have PBHs, their amplitude should grow seven orders of magnitude and reach
$\mathcal{P}_{s}\sim\mathcal{O}\left(10^{-2}\right)$ at some smaller scales.

So far, many suggestions have been proposed to overcome this problem and
provide a significant growth in the curvature power spectrum giving rise
ultimately to PBHs formation in some scales during inflation. These
suggestions are mainly characterized by different considerations for the
statistics of the curvature perturbations, because the statistical feature of
the primordial perturbations can affect on the abundance of PBHs Saito:2008em
; Byrnes:2012yx ; Young:2013oia ; Tada:2015noa ; Young:2015kda ; Young:2015cyn
; Garcia-Bellido:2017aan ; Franciolini:2018vbk ; ezquiaga2018quantum ;
Passaglia:2018ixg ; Atal:2018neu ; Belotsky:2018wph ; DeLuca:2019qsy ;
Yoo:2019pma ; Ezquiaga:2019ftu ; Serpico:2020ehh ; Braglia:2020taf ;
Figueroa:2020jkf ; Davies:2021loj ; Zhang:2021rqs . The PBHs formation due to
the drastic amplification of the Gaussian power spectrum in the framework of
the single-field inflation has been well investigated in Cai:2018tuh ;
Kamenshchik:2018sig ; Niikura:2019kqi ; Chen:2019zza ; mishra2020primordial ;
Ballesteros:2017fsr ; Bhaumik:2019tvl ; dalianis2020generalized ;
germani2017primordial ; garcia2017primordial ; Motohashi:2017kbs ;
dalianis2019primordial ; Mahbub:2019uhl ; Kawai:2021edk ; Motohashi:2019rhu ;
Liu:2020oqe ; Fu:2019ttf ; dalianis2020generalized ; Dalianis:2019vit ;
Teimoori:2021thk ; Lin:2020goi ; Yi:2020cut ; Gao:2020tsa ; Gao:2021vxb ;
Solbi:2021rse ; Solbi:2021wbo ; Teimoori:2021pte ; Heydari:2021gea ;
Heydari:2021qsr . In the literature, most studies focus on the scenarios in
which inflation takes place in a multi-phase process, and the slow-roll
condition is violated at least in one phase so that the inflaton field
experiences a so-called ultra slow-roll (or friction-dominated) phase that
leads to efficient growth in the primordial power spectrum, see e.g.
germani2017primordial ; garcia2017primordial ; Motohashi:2017kbs ;
Bhaumik:2019tvl ; Ballesteros:2017fsr ; dalianis2019primordial ;
Mahbub:2019uhl ; Kawai:2021edk ; Motohashi:2019rhu ; Liu:2020oqe ; Fu:2019ttf
; dalianis2020generalized ; Dalianis:2019vit ; Teimoori:2021thk ; Lin:2020goi
; Yi:2020cut ; Gao:2020tsa ; Gao:2021vxb ; Solbi:2021rse ; Solbi:2021wbo ;
Teimoori:2021pte ; Heydari:2021gea . For instance, in nonminimal derivative
coupling inflation, the mechanism of the gravitationally enhanced friction
leads to an ultra slow-roll inflationary phase which in turn provides a
remarkable enhancement of the power spectrum amplitude that is large enough
for PBHs production Fu:2019ttf ; dalianis2020generalized ; Dalianis:2019vit ;
Teimoori:2021thk . The origin of PBHs may come back to the existence of an
inflection-point in the inflaton potential germani2017primordial ;
garcia2017primordial ; Motohashi:2017kbs ; Bhaumik:2019tvl ;
Ballesteros:2017fsr ; dalianis2019primordial ; Mahbub:2019uhl . This idea has
already been investigated in the framework of $\alpha$-attractor inflation
dalianis2019primordial ; Mahbub:2019uhl . The $\alpha$-attractors belong to
the class of supergravity and superconformal inflationary theories which have
been proposed based on the conformal symmetry in the Jordan frame, and have a
universal attractor behavior in the Einstein frame Kallosh:2013hoa ;
Kallosh:2013yoa ; Kallosh:2014rga ; Kallosh:2014laa ; Kallosh:2015lwa ;
Linde:2015uga ; Galante:2014ifa ; Carrasco:2015rva ; Carrasco:2015pla ;
Roest:2015qya ; scalisi2015cosmological ; Eshaghi:2016kne . The free parameter
$\alpha$ is inversely proportional to the curvature of the inflaton Kähler
manifold Kallosh:2013yoa ; Galante:2014ifa . The $\alpha$-attractor models, in
the large $e$-fold number $N$ and small $\alpha$ limit, lead to the same
inflationary predictions including the scalar spectral index and tensor-to-
scalar ratio as $n_{s}=1-2/N$ and $r=12\alpha/N^{2}$, respectively, which for
$\alpha=1$, these inflationary observables coincide exactly with the
predictions of the Starobinsky model Starobinsky:1980te  as well as the Higgs
inflation scenario Bezrukov:2007ep .

PBHs also can be seeded in the multiple-field inflationary frameworks Garcia-
Bellido:1996mdl ; Sasaki:2006kq ; Kawasaki:2006zv ; Kawaguchi:2007fz ;
Frampton:2010sw ; Lyth:2011kj ; Bugaev:2011wy ; Kohri:2012yw ; Kawasaki:2012wr
; Linde:2012bt ; Bugaev:2013vba ; Clesse:2015wea ; Garcia-Bellido:2016dkw ;
Kawasaki:2016pql ; Inomata:2016rbd ; Domcke:2017fix ; Braglia:2020eai ;
Ahmed:2021ucx . In these settings, the sufficiently large peaks required for
the PBHs formation can be supplied both with curvature (adiabatic)
perturbations Polarski:1992dq ; Starobinsky:1994mh  and isocurvature ones
Starobinsky:1994mh ; Kofman:1986wm ; Kofman:1988xg ; Polarski:1994rz ;
Starobinsky:2001xq . In particular, in the context of the more traditional
hybrid inflation, the authors of Clesse:2015wea  have succeeded within a two-
field scenario to prepare a large peak in the scalar power spectrum. It should
be noted that the multi-field inflationary models are capable to produce
strong non-Gaussianities in the curvature power spectrum, and the curvaton
scenario is a well-motivated example for these class of models Young:2013oia ;
Sasaki:2006kq .

Most of the scenarios of PBHs formation from inflation are based on the
assumption of Gaussian statistics for the curvature perturbations. However,
the spectrum of the curvature fluctuations may follow from the non-Gaussian
statistics Saito:2008em ; Byrnes:2012yx ; Young:2013oia ; Tada:2015noa ;
Young:2015kda ; Young:2015cyn ; Garcia-Bellido:2017aan ; Franciolini:2018vbk ;
ezquiaga2018quantum ; Passaglia:2018ixg ; Atal:2018neu ; Belotsky:2018wph ;
DeLuca:2019qsy ; Yoo:2019pma ; Ezquiaga:2019ftu ; Serpico:2020ehh ;
Braglia:2020taf ; Figueroa:2020jkf ; Davies:2021loj ; Zhang:2021rqs .
Recently, the authors of Garcia-Bellido:2017aan  considered the effect of
primordial non-Gaussianity on the PBHs abundance and concluded that the non-
Gaussian spectrum can provide a similar abundance for PBHs from a lower
amplitude of the scalar perturbations relative to the Gaussian spectrum. It is
crucial to point out that due to the impact of non-Gaussianity on the PBHs
abundance, the formation of PBHs can be regarded as an instrument to evaluate
the amplitude and non-Gaussianity of the curvature perturbations at the small
scales.

The production of PBHs due to enhancement of the power spectrum of the
curvature perturbations during inflation may be accompanied by induction of
secondary GWs Matarrese:1997ay ; Mollerach:2003nq ; Ananda:2006af ;
Baumann:2007zm ; Saito:2008jc ; Saito:2009jt ; Bugaev:2009zh ; Bugaev:2010bb ;
Alabidi:2012ex ; Nakama:2016gzw ; Peirone:2017vcq ; Cheng:2018yyr .
Consequently, if we observe PBHs or secondary GWs, then it is possible to use
them to constrain the amplification of the curvature perturbations power
spectrum during inflation which provides valuable information about the
physics of the early universe for us. Various ideas are suggested so far to
generate secondary GWs by virtue of the production of PBHs from inflation. One
idea is application of a phenomenological delta function giving rise to the
required growth in the primordial power spectrum Saito:2008jc ;
Bartolo:2018evs ; Cai:2018dig . This goal can also be achieved appropriately
by using a broken power-law Lu:2019sti  or Gaussian power spectrum
Namba:2015gja ; Garcia-Bellido:2017aan ; Lu:2019sti . The running-mass model
Stewart:1996ey ; Drees:2011hb ; Datta:2019euh , the axion-curvaton model
Kasuya:2009up ; Kawasaki:2012wr , oscillons after inflation Easther:2006vd ;
Antusch:2016con ; Liu:2017hua , and parametric resonance due to the double
inflation Kawasaki:2006zv ; Kawasaki:2016pql  are another mechanisms leading
to the generation of PBHs and secondary GWs in the context of inflationary
cosmology. Note that, like the abundance of PBHs, the spectrum of induced GWs
depends on both the amplitude and statistics of the scalar quantum
fluctuations Garcia-Bellido:2017aan . However, in contrast to the PBHs
abundance which is very sensitive to the statistics of the curvature
perturbations, the amplitude of the secondary GWs is not very sensitive to
this feature and it is hardly affected by the non-Gaussianity of the
primordial curvature perturbations Zhang:2021vak .

It is possible to produce PBHs only with a trivial modification in the
potential shape, like the approach of mishra2020primordial  in which the
source of DM originates from a tiny bump in the potential of the
$\alpha$-attractor inflation. In this scenario, the bump behaves like a speed-
breaker and slows down the scalar field evolution. It leads to a sharp
increase in the amplitude of the curvature perturbations that is required for
PBHs production.

In the present work, we are interested to calculate the primordial non-
Gaussianity on the different scales and also the secondary gravitational waves
within the setup of the $\alpha$-attractor inflation model in the presence of
a tiny bump in its potential. The main reason for us to concentrate this model
originates from the fact that this model firstly has well-based theoretical
motivation from string theory Kallosh:2013hoa ; Kallosh:2013yoa . Secondly,
this model is in excellent agreement with the latest CMB data released from
the Planck 2018 akrami2020planck  measurements. The prediction of this model
for the scalar spectral index $n_{s}$ and the tensor-to-scalar ratio $r$, is
located satisfactorily within the 68% CL marginalized joint region of the
Planck 2018 data akrami2020planck , as it has been shown in
mishra2020primordial . Most of the inflationary models predicts rather large
values for $r$ so that they are ruled out explicitly by the current
observations. The $\alpha$-attractor scenario gives a small value for this
observable, and therefore it remains in the class of acceptable inflationary
models in light of the current data. This model is further capable to provide
smaller values for $r$, and so if the future data prefer less value for this
observable, the $\alpha$-attractor model can still be remain compatible with
the observations.

Although the production of PBHs from $\alpha$-attractor inflation has already
been investigated in mishra2020primordial , there exist some aspects of this
model that has not been regarded so far. Moreover, although the authors of
mishra2020primordial  have computed the curvature power spectrum
$\mathcal{P}_{s}$ of this model, they did not compare their findings with the
current upper bounds on this spectrum. At the CMB scales, this spectrum is
tightly constrained by the Planck 2018 observations akrami2020planck . At the
smaller scales, there are still a lot of constraints on $\mathcal{P}_{s}$,
such as the bounds from CMB spectral distortions Kohri:2014lza ; Chluba:2012we
, Acoustic Reheating (AR) Inomata:2016uip , European Pulsar Timing Array
(EPTA) observations lentati2015european , Big Bang Nucleosynthesis (BBN)
Kohri:2018awv , and advanced Laser Interferometer Gravitational wave
Observatory (aLIGO) O2 LIGOScientific:2017zlf ; LIGOScientific:2017vwq . In
addition, there are several sensitivity regions for the curvature perturbation
spectrum from different detectors like Square Kilometre Array (SKA)
Moore:2014lga ; Janssen:2014dka , the Big Bang Observer (BBO) Yagi:2011wg ;
phinney2003big , The DECi-hertz Interferometer Gravitational wave Observatory
(DECIGO) Seto:2001qf ; Yagi:2011wg , Einstein Telescope (ET) Punturo:2010zz ,
Laser Interferometer Space Antenna (LISA) LISA:2017pwj , and aLIGO design
KAGRA:2013rdx , that can be used to check viability of PBHs scenarios when
their data are released. Especially, the sensitivity of ET on the scales of
$10^{14}-10^{19}\mathrm{Mpc}^{-1}$ is high enough to evaluate the amplitude of
$\sim 10^{-5}$ for the spectrum of the curvature perturbations. Therefore, it
will be much useful if the result of our model for the scalar power spectrum
is compared with these sensitivity regions.

Besides, it will be useful if the results of the PBHs abundances
$f_{\rm{PBH}}$ in different mass scales are compared with the current
constraints on $f_{\rm{PBH}}$ such as the constraints of Extra-Galactic gamma-
ray (EG$\gamma$) Carr:2009jm , White Dwarf explosion (WD) Graham:2015apa ,
microlensing events with Subaru Hyper Suprime-Cam (HSC) Niikura:2017zjd ,
Kepler milli/microlensing (Kepler) Griest:2013esa , the Earth Resources
Observation System(EROS)/Massive Astrophysical Compact Halo Object (MACHO)
microlensing (EROS) EROS-2:2006ryy , the Optical Gravitational Lensing
Experiment (OGLE) Niikura:2019kqi , dynamical heating of Ultra-Faint Dwarf
galaxies (UFD) Brandt:2016aco , and accretion constraints by CMB Ali-
Haimoud:2016mbv ; aloni2017cosmic ; Poulin:2017bwe . The results of the model
also can be appraised in light of the aLIGO O1 & O2 LIGOScientific:2016dsl ;
LIGOScientific:2017zlf ; LIGOScientific:2017vwq , aLIGO design KAGRA:2013rdx ,
ET and Cosmic Explorer (CE) LIGOScientific:2016wof  detectable limits.
Accordingly, in the present work, we examine the possibility of PBHs
production with masses around ${\cal O}(10)M_{\odot}$, ${\cal
O}(10^{-5})M_{\odot}$, and ${\cal O}(10^{-13})M_{\odot}$ in the framework of
the $\alpha$-attractor inflation with a tiny bump in the potential, in light
of the current constraints on $\mathcal{P}_{s}$ and $f_{\rm{PBH}}$. We show
that for the model satisfy the current constraints on the curvature power
spectrum and the PBHs abundance, the model parameters should be somewhat
different from those used in mishra2020primordial .

Furthermore, we study the subjects of non-Gaussianity and secondary GWs in the
process of PBHs formation in the framework of $\alpha$-attractor Inflation.
These concepts are very important issues which have absorbed much attentions
recently (see, e.g. Saito:2008em ; Byrnes:2012yx ; Young:2013oia ;
Tada:2015noa ; Young:2015kda ; Young:2015cyn ; Garcia-Bellido:2017aan ;
Franciolini:2018vbk ; ezquiaga2018quantum ; Passaglia:2018ixg ; Atal:2018neu ;
Belotsky:2018wph ; DeLuca:2019qsy ; Yoo:2019pma ; Ezquiaga:2019ftu ;
Serpico:2020ehh ; Braglia:2020taf ; Figueroa:2020jkf ; Davies:2021loj ;
Zhang:2021rqs ). The non-Gaussianity parameter is an important observable that
can provide valuable information about the inflationary era of the early
Universe. We also investigate the secondary GWs in this framework. We
calculate the present fractional energy density of these GWs and compare our
findings with the sensitivity regions of different GWs detectors.
Additionally, we estimate the slope of the secondary GWs spectrum in different
frequency ranges. This feature of the secondary GWs spectrum also may be
assessed by the forthcoming data, and then can be used to discriminate between
the different mechanisms of PBHs production from inflation. In the recent
years, several GWs detectors have been designed, and their data can be used to
scrutinize the PBHs seeded from the primordial perturbations. Among these
detectors, one can refer to aLIGO O2 & O3 runs LIGOScientific:2017zlf ;
LIGOScientific:2017vwq ; KAGRA:2021kbb , the HLV (LIGO-Hanford, LIGOLivingston
and Virgo), the mid-scale upgrade of aLIGO (A+ detectors), ET Punturo:2010zz ,
LISA LISA:2017pwj , EPTA lentati2015european , North American Nanohertz
Observatory for Gravitational Waves (NANOGrav) NANOGrav:2020gpb ,
International Pulsar Timing Array (IPTA), and the Square Kilometer Array (SKA)
Moore:2014lga ; Janssen:2014dka ; Schmitz:2020syl . Our results in this paper
for the abundance of PBHs and the induced GWs spectrum can be produced by our
numerical codes which are available
online555https://github.com/krezazadeh/PBHs_alpha_attractor.

To compute the non-Gaussianity parameter in various cosmological scales, we
employ the numerical code BINGO which is developed by Hazra et al.
Hazra:2012yn  and is publicly available666https://github.com/dkhaz/bingo. This
code computes the inflationary non-Gaussianities in different triangle
configurations. On the scales far from the bump, the slow-roll regime is
completely valid, and hence the non-Gaussianity is small as expected for the
single-field standard inflation model Maldacena:2002vr ; Acquaviva:2002ud ;
Rigopoulos:2004ba . But on the scales around the bump, the slow-roll condition
may be disturbed, and therefore the non-Gaussianity may take large values. The
values of non-Gaussianity on these scales may be extensively comparable in
front of the values evaluated on the CMB scales ($k_{*}\sim
0.05\,\mathrm{Mpc}^{-1}$), which are strongly constrained by the observational
data of Planck satellite Planck:2019kim . However, the prediction of large
non-Gaussianity on the scales beyond the CMB scales may be checked with more
precise measurements in the future. In this way, one can discriminate between
the inflationary models proposed to explain PBHs production. Besides, such
measurement provides a useful insight for us to understand whether the PBHs
can be regarded as the origin of DM or not.

The rest of this paper is structured as follows: first, in Sec. II, we present
a brief review of the base formulas for the calculation of the mass and
abundance of PBHs. In Sec. III, we review briefly single field inflationary
model. In Sec. IV, we discuss the process of PBH production in the framework
of the $\alpha$-attractor scenario with a tiny bump in the potential. Then, in
Sec. V, we calculate the non-Gaussianity parameter for three different cases
of the model which have different values for the parameters. We subsequently
investigate the secondary GWs in our model in Sec VI. Finally, in Sec. VII, we
summarize our concluding remarks.

##  II Abundance of Primordial Black Holes

If the amplitude of the scalar modes of perturbations produced during the
phase of inflation is large enough, they may gravitationally collapse to form
PBHs after horizon re-entry in the radiation dominated era. The relation
between the mass of formed PBHs, $M(k)$ and the horizon mass, $M_{\rm hor}$ is
Gong:2017qlj

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

where $\gamma$ is the efficiency factor and $g_{*}$ is the effective degrees
of freedom for the energy density. These parameters usually assumed to be
$\gamma\simeq 0.2$ and $g_{*}=107.5$ carr1975primordial .

We assume that the curvature perturbations obey Gaussian statistics. In this
way the production rate of PBHs with mass $M(k)$ has the following form
Young:2014ana ; ozsoy2018mechanisms ; Tada:2019amh

$\beta(M)=\sqrt{\frac{2}{\pi}}\frac{\sigma(M)}{\delta_{c}}\exp{\left(-\frac{\delta_{c}^{2}}{2\sigma^{2}(M)}\right)},$
(2)

where $\delta_{c}$ is the threshold of the density perturbation for the PBHs
formation, and we assume to be $\delta_{c}=0.4$ Musco:2012au ; Harada:2013epa
; Escriva:2019phb . $\sigma^{2}(k)$ is the coarse-grained variance of the
density contrast smoothed on a scale $k$ and is given by Young:2014ana ;
ozsoy2018mechanisms

$\sigma^{2}(k)=\int{\frac{dq}{q}}\,W^{2}(q/k)\frac{16}{81}(q/k)^{4}{\cal
P}_{s}(q),$ (3)

where ${\cal P}_{s}$ is the power spectrum of the curvature perturbations and
$W(x)$ is the window function. In this work we choose Gaussian window function
as $W(x)=\exp{(-x^{2}/2)}$ that is very popular in the literature. Other
choices for window functions have been studied in Young:2019osy ; Ando:2018qdb
.

The current energy fraction of PBHs with a mass $M(k)$ over the total DM is
given by Gong:2017qlj ; Carr:2016drx

$f_{\rm{PBH}}(M)\equiv\frac{\Omega_{\rm{PBH}}}{\Omega_{\rm{DM}}}=\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_{\rm{DM}}h^{2}}\right)\left(\frac{M}{M_{\odot}}\right)^{-1/2},$
(4)

where $\Omega_{\rm{DM}}$ is the current density parameter of DM and we take
its value as $\Omega_{\rm{DM}}h^{2}\simeq 0.12$ from the Planck 2018 results
akrami2020planck .

##  III Single Field Inflationary Model

The standard single field inflationary paradigm is described by the action
Baumann:2009ds

$S=\int{\rm d}^{4}x\sqrt{-g}\left[\frac{1}{2}R+X-V(\phi)\right],$ (5)

where $g$ is the determinant of the metric $g_{{\mu}{\nu}}$, $R$ is the Ricci
scalar,
$X\equiv-\frac{1}{2}g^{\mu\nu}\partial_{\mu}{\phi}\partial_{\nu}{\phi}$ is the
kinetic term, and $V(\phi)$ is the scalar field potential. Throughout this
paper, we set the reduced Planck mass equal to unity, i.e., $M_{\rm{pl}}=(8\pi
G)^{-1/2}=1$.

For the flat Friedmann-Robertson-Walker (FRW) metric $g_{{\mu}{\nu}}={\rm
diag}\Big{(}-1,a^{2}(t),a^{2}(t),a^{2}(t)\Big{)}$, we have
$X=\dot{\phi}^{2}/2$ and the background equations turn into Baumann:2009ds

$\displaystyle 3H^{2}-\frac{1}{2}\dot{\phi}^{2}-V(\phi)=0,$ (6) $\displaystyle
2\dot{H}+3H^{2}+\frac{1}{2}\dot{\phi}^{2}-V(\phi)=0,$ (7)
$\displaystyle\ddot{\phi}+3H\dot{\phi}+V_{,\phi}=0,$ (8)

where $H\equiv\dot{a}/a$ is the Hubble parameter and the overdot represents
the derivative with respect to the cosmic time $t$.

In the context of standard single field inflation, the slow-roll parameters
are defined as Baumann:2009ds

$\varepsilon_{H}\equiv-\frac{\dot{H}}{H^{2}},\hskip
14.22636pt\eta_{H}\equiv-\frac{\ddot{\phi}}{H\,\dot{\phi}}.$ (9)

The first relation in Eq. (9) indicates that to have inflation, the condition
$\varepsilon<1$ should be satisfied. Note that, during the slow-roll regime,
two parameters $\varepsilon_{H}$ and $\eta_{H}$ should be much smaller than
unity. Using Eqs. (6) and (7), one can easily find the slow-roll parameter
$\varepsilon_{H}$ as

$\varepsilon_{H}=\frac{\dot{\phi}^{2}}{2H^{2}}.$ (10)

In the framework of standard inflation, the power spectrum of the scalar
perturbation ${\cal P}_{s}$, is in the following form Baumann:2009ds

${\cal P}_{s}=\frac{H^{2}}{8\pi^{2}}\frac{1}{\varepsilon_{H}}\Big{|}_{k=aH},$
(11)

where ${\cal P}_{s}$ is computed at the time of sound horizon exit, i.e.,
$k=aH$ in which $k$ denotes the comoving wavenumber. The observational value
of the amplitude of scalar perturbations at the CMB pivot scale
$k_{*}=0.05\,{\rm Mpc}^{{\rm-1}}$ reported by the Planck team is ${\cal
P}_{s}(k_{*})\simeq 2.1\times 10^{-9}$ akrami2020planck .

Under the slow-roll approximation, the field equations (6) and (8) reduce to

$\displaystyle 3H^{2}\simeq V(\phi),$ (12) $\displaystyle
3H\dot{\phi}+V_{,\phi}\simeq 0.$ (13)

Using Eqs. (12) and (13), the slow-roll parameter $\varepsilon_{H}$ in (10)
can be written as

$\varepsilon_{H}\simeq\varepsilon_{V}\equiv\frac{1}{2}\left(\frac{V_{,\phi}}{V}\right)^{2}.$
(14)

With the help of the first Friedmann Eq. (12) and Eq. (14), the power spectrum
(11) in the slow-roll limit reduces to

${\cal P}_{s}\simeq\frac{V^{3}}{12\pi^{2}V_{,\phi}^{2}}.$ (15)

In the slow-roll approximation, the Hubble parameter $H$ and its variation
change much slower than the scale factor $a$ of the Universe Garriga:1999vw .
Therefore, using the relation $k=aH$, we have $d\ln k\approx Hdt$. Applying
this approximation and the definition $n_{s}-1\equiv d\ln{\cal P}_{s}/d\ln k$,
one can easily find $n_{s}-1\simeq\dot{{\cal P}_{s}}/(H{\cal P}_{s})$. Using
this relation and Eqs. (12), (13), (15), and the definition (14), we can find
the scalar spectral index $n_{s}$ as

$n_{s}-1\simeq 2\eta_{V}-6\varepsilon_{V},$ (16)

where

$\eta_{V}\equiv\frac{V_{,\phi\phi}}{V}=\eta_{H}-\varepsilon_{H}.$ (17)

The observational value of the scalar spectral index is $n_{s}=0.9627\pm
0.0060$ (68% CL, Planck 2018 TT+lowE) akrami2020planck . Using the
approximation $d\ln k\approx Hdt$ and Eq. (13), the running of the scalar
spectral index as

$\frac{d{n_{s}}}{d\ln
k}\simeq-\left(\frac{V_{,\phi}}{V(\phi)}\right)n_{s,\phi}.$ (18)

The observational constraint on the running of the scalar spectral index is
about $d{n_{s}}/d\ln k=-{\rm{0}}{\rm{.0078}}\pm 0.0082$ (68% CL, Planck 2018
TT+lowE) akrami2020planck .

The amount of inflation is usually described by the $e$-fold number $N$ which
is defined as

$N\equiv\ln\left(\frac{a_{e}}{a}\right),$ (19)

where $a_{e}$ is the scale factor at the end of inflation. The definition 19
is equivalent to

$dN=-H\,dt.$ (20)

It is important to note that the number of $e$-folds between the horizon exit
and the end of inflation should be about $50-60$ Dodelson:2003vq ;
Liddle:2003as .

The tensor power spectrum in the framework of the standard canonical
inflationary model is given by

${\cal P}_{t}=\frac{H^{2}}{2\pi^{2}}\Big{|}_{k=aH}.$ (21)

From Eqs. (11) and (21), we get the tensor-to-scalar ratio in the slow-roll
regime as

$r=16\varepsilon_{H}\simeq 16\varepsilon_{V}.$ (22)

The Planck 2018 data determined an upper bound on the tensor-to-scalar ratio
as $r<0.0654$ (68% CL, Planck 2018 TT+lowE) akrami2020planck .

##  IV $\alpha$-attractor inflation with a tiny bump in the potential

Including a tiny bump, the general form of the potential can be written as the
following form mishra2020primordial

$V(\phi)=V_{b}(\phi)\left[1+\varepsilon(\phi)\right],$ (23)

where $V_{b}$ is the base potential of the model and here it is taken in the
form of $\alpha$-attractor potential Kallosh:2013hoa ; Kallosh:2013yoa

$V_{b}(\phi)=V_{0}\tanh^{2n}\left(\frac{\phi}{\sqrt{6\alpha}}\right).$ (24)

This is the original inflationary potential responsible for providing
primordial quantum fluctuations compatible with the CMB constraints on the
inflationary observables including the scalar spectral index $n_{s}$ and
tensor-to-scalar ratio $r$. Following mishra2020primordial , we set two free
parameters $\alpha$ and $n$ to be $\alpha=1$ and $n=1$. The function
$\varepsilon$ in (23) specifies the shape of the bump, and it is considered as
mishra2020primordial

$\varepsilon(\phi)=A\cosh^{-2}\left(\frac{\phi-\phi_{0}}{\sigma}\right),$ (25)

where $A$, $\sigma$, and $\phi_{0}$ characterize the height, width, and
position of the peak, respectively. As discussed in mishra2020primordial , the
presence of such a modification in the potential function, leads to
sufficiently large enhancement in the curvature power spectrum required for
the PBHs formation.

In mishra2020primordial , the authors have considered three cases for the
model that have different values for these parameters, and estimated the PBHs
abundance for each case. In the present paper, we take the amounts of these
parameters as listed in Table 1. The value of $V_{0}$ is fixed by imposing the
CMB normalization at the pivot scale $k_{*}=0.05\,{\rm Mpc}^{{\rm-1}}$
corresponding to $N_{*}=60$. It is worth to mentioning that the values for
Case 1 and Case 3 are very close to those adopted in mishra2020primordial ,
but we have changed the values slightly here in order to the cases satisfy the
existing constraints on the curvature power spectrum as well as the bounds on
the PBHs abundances at different cosmological scales. In addition, we have
considered another parameter set as Case 2 whose parameters are completely
different with those adopted in mishra2020primordial . We will show this case
provides properly the current observational bounds too. We further compute the
inflationary observables $n_{s}$, $r$, and $dn_{s}/d\ln k$ at the
$k_{*}=0.05\,{\rm Mpc}^{{\rm-1}}$, and also quantities relevant for producing
PBHs of the considered three cases. The results are presented in Table 2.

Table 1: The successful parameter sets for the potential (23) that can produce
PBHs. The value of $V_{0}$ is fixed by imposing the CMB normalization at the
pivot scale $k_{*}=0.05\,{\rm Mpc}^{{\rm-1}}$ corresponding to $N_{*}=60$.

$\qquad\\#\qquad$ | $\qquad\phi_{0}\qquad$ | $\qquad\qquad A\qquad\qquad$ | $\qquad\qquad\sigma\qquad\qquad$ | $\qquad V_{0}\qquad$
---|---|---|---|---
Case 1 | $4.850001$ | $2.044880\times 10^{-3}$ | $2.524999\times 10^{-2}$ | $1.448\times 10^{-10}$
Case 2 | $5.301327$ | $1.080552\times 10^{-3}$ | $1.940550\times 10^{-2}$ | $1.495\times 10^{-10}$
Case 3 | $5.580000$ | $6.402075\times 10^{-4}$ | $1.428997\times 10^{-2}$ | $1.420\times 10^{-10}$

The evolution of the scalar field $\phi$ as a function of the $e$-fold number
$N$, for Case 1 (solid line), Case 2 (dashed line), and Case 3 (dash-dotted
line) are plotted in Fig. 1 by solving the full equations of motion (6)-(8)
numerically. We set the initial conditions by using Eqs. (12) and (13) at
$N_{*}=60$. In this figure we see that a plateau-like region appears at
$\phi=\phi_{0}$ corresponding to $20\lesssim N\lesssim 29$ for Case 1,
$26\lesssim N\lesssim 39$ for Case 2, and $37\lesssim N\lesssim 46$ for Case
3. During this region, the inflaton undergoes a period of ultra slow-roll
inflation in which its evolution is slowing down by increasing friction which
leads to enhance the curvature power spectrum by several orders of magnitude
that is required to PBHs formation.

Table 2: Results of the scalar spectral index $n_{s}$, the tensor-to-scalar
ratio $r$, the running of the scalar spectral index $dn_{s}/d\ln k$, the mass
of the corresponding PBHs $M_{\rm{PBH}}^{\rm{peak}}$, and the PBH fractional
abundance $f_{\rm{PBH}}$ for the three cases of Table 1. The inflationary
observables $n_{s}$, $r$, and $dn_{s}/d\ln k$ are found at the CMB scale
corresponding to $N_{*}=60$.

# |  $n_{s}$ |  $r$ |  $dn_{s}/d\ln k$ |  $M_{\rm{PBH}}^{\rm{peak}}/M_{\odot}$ |  $f_{\rm{PBH}}^{\rm{peak}}$
---|---|---|---|---|---
Case 1 |  $0.96031$ |  $0.00465$ |  $-0.000783$ |  $1.66\times 10^{-13}$ |  $0.9550$
Case 2 |  $0.96003$ |  $0.00471$ |  $-0.000794$ |  $1.28\times 10^{-5}$ |  $0.0217$
Case 3 |  $0.96028$ |  $0.00465$ |  $-0.000784$ |  $20.98$ |  $0.0070$

Figure 1: Evolution of the scalar field $\phi$ as a function of the $e$-fold
number $N$ for Case 1 (solid line), Case 2 (dashed line), and Case 3 (dash-
dotted line).

In Fig. 2 we plot the evolution of the slow-roll parameters $\varepsilon_{H}$
(left panel) and $\eta_{H}$ (right panel) as functions of $N$ for the three
cases of Table 1. From the left panel of this figure we see that
$\varepsilon_{H}$ remains below unity up until the end of inflation
$(\varepsilon_{H}<1)$ but the right panel shows that the slow-roll
approximation $\eta_{H}\ll 1$ is violated on the scales of the bump. Of
course, it is worth mentioning that for Case 2 the first slow-roll parameter
doesn’t reach unity at the end of inflation, i.e., $N=0$ and inflation ends
after a finite duration. It continues $3.307$ $e$-folds due to strong slow
down of inflaton in the ultra slow-roll phase.

A quick look on Fig. 2 demonstrates that at the time of sound horizon exit
corresponding to $N_{*}=60$, the slow-roll approximation remains valid.
Therefore, we can compute the values of scalar spectral index $n_{s}$, the
tensor-to-scalar ratio $r$, and the running of the scalar spectral index
$dn_{s}/d\ln k$ with the help of Eqs. (16), (18), (22), and also using Eqs.
(14), (17) for the three cases of Table 1. The numerical results are presented
in Table 2. We see that the scalar spectral index, the tensor-to-scalar ratio,
and the running of the scalar spectral index for all three cases are
$n_{s}\approx 0.96$, $r\approx 0.0047$, and $dn_{s}/d\ln k\approx-0.0008$,
respectively. These values are in agreement with the $68\%$ CL constraints of
Planck 2018 TT+lowE data akrami2020planck .

Figure 2: Evolution of the first slow-roll parameter $\varepsilon_{H}$ (left)
and the second slow-roll parameter $\eta_{H}$ (right) versus the $e$-fold
number $N$ for Case 1 (solid line), Case 2 (dashed line), and Case 3 (dash-
dotted line).

As discussed above, the slow-roll approximation in this model is violated at
the scales of the bump. Therefore, in order to examine the evolution of the
curvature perturbation, it is essential to solve numerically the Mukhanov-
Sasaki equation which is in the following form

$v_{k}^{\prime\prime}+\left(k^{2}-\frac{z^{\prime\prime}}{z}\right)v_{k}=0,$
(26)

where the prime indicates the derivative with respect to the conformal time
$\tau$. Also, $v$ is the canonically normalized variable for the curvature
perturbation $\mathcal{R}$, and it is defined as $v\equiv z\mathcal{R}$, with
$z=a\dot{\phi}/H$. To solve the differential equation (26), we use the Bunch-
Daveis vacuum initial condition,

$v_{k}\to\frac{1}{\sqrt{2k}}e^{-ik\tau},$ (27)

which is quite valid for the scales deep inside the horizon. The solution of
the Mukhanov-Sasaki equation (26) is usually used to calculate the power
spectrum of the scalar perturbation which is given by

$\mathcal{P}_{s}=\frac{k^{3}}{2\pi^{2}}\frac{\left|v_{k}\right|^{2}}{z^{2}}.$
(28)

We evaluate the scalar power spectrum (28) for the three cases of Table 1, and
plot their scalar power spectrum versus comoving wavenumber $k$ in Fig. 3. The
results are shown by the solid (Case 1), dashed (Case 2), and dash-dotted
(Case 3) black lines. Note that through this paper, we assume that the pivot
scale $k_{*}\sim 0.05\,\mathrm{Mpc}^{-1}$ exits the Hubble horizon at
$N_{*}=60$ $e$-folds before the end of inflation. As we see in Fig. 3, the
scalar power spectra of Cases 1, 2, and 3 represent their peaks at the scales
$\mathcal{O}(10^{12})$, $\mathcal{O}(10^{8})$, and $\mathcal{O}(10^{5})$
$\mathrm{Mpc}^{-1}$, respectively. In these scales the amplitude of the power
spectra get amplified considerably and reach the magnitudes of order
$10^{-2}$, while on the CMB scales ($k_{*}\sim 0.05\,\mathrm{Mpc}^{-1}$), they
approach the value $\mathcal{P}_{s*}\sim\mathcal{O}\left(10^{-9}\right)$ which
is compatible with the Planck 2018 measurements akrami2020planck . The
numerical results for the peak of the scalar power spectrum are presented in
Table 3.

In Fig. 3, the existing and expected limits on the curvature power spectrum at
the small scales are also specified. The colored shaded regions are excluded
by the current observations, whereas the colored dashed lines show the
expected limits from future experiments Inomata:2018epa . The dark green,
brown, pink, blue, light green and purple shaded regions indicate the
constraints from the Cosmic Microwave Background/Large-Scale Structure
(CMB/LSS) observations hunt2015search , CMB spectral distortions Kohri:2014lza
; Chluba:2012we , Acoustic Reheating (AR) Inomata:2016uip , European Pulsar
Timing Array (EPTA) observations lentati2015european , Big Bang
Nucleosynthesis (BBN) Kohri:2018awv , and advanced Laser Interferometer
Gravitational wave Observatory (aLIGO) O2 LIGOScientific:2017zlf ;
LIGOScientific:2017vwq , respectively. The sensitivity curves of the Square
Kilometre Array (SKA) Moore:2014lga ; Janssen:2014dka , the Big Bang Observer
(BBO) Yagi:2011wg ; phinney2003big , The DECi-hertz Interferometer
Gravitational wave Observatory (DECIGO) Seto:2001qf ; Yagi:2011wg , Einstein
Telescope (ET) Punturo:2010zz , Laser Interferometer Space Antenna (LISA)
LISA:2017pwj , and aLIGO design KAGRA:2013rdx  are also shown by colored
dashed lines. We see in the figure that the power spectra of all the three
cases satisfy properly the existing constraints, and hence our setup can
successfully explain the enhancement of the scalar perturbations required for
the PBHs formation. From the figure, we further see that the predictions of
our setup are located inside the sensitivity region of SKA, LISA, BBO and
DECIGO detectors. Therefore, the viability of our model can be assessed in
light of the future observations.

In Fig. 3, we also see that, before the enhancement, the scalar power spectrum
shows a sudden declination in its shape, and this is a special feature coming
from the existence of the tiny bump in the potential of the $\alpha$-attractor
model. We will discuss in the next section that this feature is the one that
sources a considerable enhancement for the equilateral non-Gaussianity in this
model.

Figure 3: The curvature power spectra founded by solving numerically the
Mukhanov-Sasaki equation versus the comoving wavenumber $k$ for Case 1 (solid
black line), Case 2 (dashed black line), and Case 3 (dash-dotted black line).
The colored shaded regions are excluded by the current observations and the
colored dashed lines show the expected limits from the future experiments
Inomata:2018epa . Table 3: The comoving wavenumber and amplitude of the peak
of the scalar power spectrum for Case 1, 2, and 3.

# |  $k^{\rm{peak}}\,\left[\mathrm{Mpc}^{-1}\right]$ | $\mathcal{P}_{s}^{\rm{peak}}$
---|---|---
Case 1 | $4.833\times 10^{12}$ | $3.079\times 10^{-2}$
Case 2 | $5.254\times 10^{8}$ | $3.518\times 10^{-2}$
Case 3 | $3.886\times 10^{5}$ | $4.740\times 10^{-2}$

Using the exact power spectrum obtained by solving the Mukhanov-Sasaki
equation numerically and with the help of Eqs. (1)-(4), we find the mass of
the formed PBHs and their abundance for parameter sets 1, 2, and 3. The
numerical results are tabulated in Table 2 which are demonstrated that the
model can generate PBHs in a wide range of mass scales. As we see, in Case 1,
our model predicts PBHs with the mass $M\simeq 1.66\times 10^{-13}M_{\odot}$
and PBH abundance $f_{\rm{PBH}}\simeq 0.9550$ which means that the produced
PBHs in this class constitute around $0.96\%$ of DM content of the present
Universe. It is an interesting result and implies that these PBHs can be taken
as a suitable candidate for DM.

In Fig. 4, we plot the abundance of formed PBHs, $f_{\rm{PBH}}$ versus the
mass for Case 1 (solid black line), Case 2 (dashed black line) and Case 3
(dash-dotted black line). In this figure, the current constraints on
$f_{\rm{PBH}}$ are also shown Chen:2019irf . The black, purple, magenta, and
orange curves indicate the results of targeted search from aLIGO O1 & O2
LIGOScientific:2016dsl ; LIGOScientific:2017zlf ; LIGOScientific:2017vwq ,
aLIGO design KAGRA:2013rdx , ET and Cosmic Explorer (CE)
LIGOScientific:2016wof , respectively Chen:2019irf . The dashed lines
illustrate the detectable limits of $f_{\rm{PBH}}$ from these detectors. The
gray vertical line shows that the constraints from the targeted search are
only valid for the sub-solar mass PBHs Chen:2019irf . The red curve is the
updated upper bound of $f_{\rm{PBH}}$ constrained by the non-detection of
Stochastic Gravitational Wave Background (SGWB) from both aLIGO O1 & O2
searches Chen:2019irf . The other shaded regions indicate the current
observational constraints on the abundance of PBHs: Extra-Galactic gamma-ray
(EG$\gamma$) Carr:2009jm , White Dwarf explosion (WD) Graham:2015apa ,
microlensing events with Subaru Hyper Suprime-Cam (HSC) Niikura:2017zjd ,
Kepler milli/microlensing (Kepler) Griest:2013esa , the Earth Resources
Observation System(EROS)/Massive Astrophysical Compact Halo Object (MACHO)
microlensing (EROS) EROS-2:2006ryy , the Optical Gravitational Lensing
Experiment (OGLE) Niikura:2019kqi , dynamical heating of Ultra-Faint Dwarf
galaxies (UFD)Brandt:2016aco , and accretion constraints by CMB Ali-
Haimoud:2016mbv ; aloni2017cosmic ; Poulin:2017bwe . From the figure, we see
that the predictions of our model for the three cases are in agreement with
the existing constraints. Moreover, the figure shows that the obtained results
for Case 3 can be located within the detectable limits of $f_{\rm{PBH}}$ by
the targeted search from aLIGO O1 & O2, aLIGO design, ET and CE.

Figure 4: The fractional abundance of PBHs as a function of PBH mass for Case
1 (solid black line), Case 2 (dashed black line), and Case 3 (dash-dotted
black line). The black, purple, magenta, and orange curves indicate the
results of targeted search from aLIGO O1 & O2, aLIGO design, ET Punturo:2010zz
, and CE LIGOScientific:2016wof , respectively Chen:2019irf . The dashed lines
illustrate the detectable limits of $f_{\rm{PBH}}$ from these detectors. The
gray vertical line shows that the constraints from the targeted search are
only valid for the sub-solar mass PBHs Chen:2019irf . The red curve is the
updated upper bound of $f_{\rm{PBH}}$ constrained by the non-detection of SGWB
from both aLIGO O1 & O2 searches. The other shaded regions indicate the
current observational constraints on the abundance of PBHs: Extra-Galactic
gamma-ray (EG$\gamma$) Carr:2009jm , White Dwarf explosion (WD) Graham:2015apa
, microlensing events with Subaru Hyper Suprime-Cam (HSC) Niikura:2017zjd ,
Kepler milli/microlensing (Kepler) Griest:2013esa , the Earth Resources
Observation System/Massive Astrophysical Compact Halo Object (EROS/MACHO)
microlensing (EROS) EROS-2:2006ryy , the Optical Gravitational Lensing
Experiment (OGLE) Niikura:2019kqi , dynamical heating of Ultra-Faint Dwarf
galaxies (UFD)Brandt:2016aco , and accretion constraints by CMB Ali-
Haimoud:2016mbv ; aloni2017cosmic ; Poulin:2017bwe .

##  V Calculation of non-Gaussianity

Planck 2018 collaboration Planck:2019kim  has analyzed the cosmic microwave
background temperature and $E$-mode polarization maps to obtain constraints on
primordial non-Gaussianity. Their combined temperature and polarization
analysis produces the following final results for the local, equilateral, and
orthogonal bispectrum amplitudes: $f_{\mathrm{NL}}^{\mathrm{local}}=-0.9\pm
5.1$, $f_{\mathrm{NL}}^{\mathrm{equil}}=-26\pm 47$, and
$f_{\mathrm{NL}}^{\mathrm{ortho}}=-38\pm 24$ (68% CL) Planck:2019kim . Here,
we are interested to calculate the non-Gaussianity parameter for the
$\alpha$-attractor inflation model with a tiny bump in its potential, and
compare our results at the CMB scales with the Planck 2018 constraints
Planck:2019kim . To do so, we use the publicly available computational code
BINGO Hazra:2012yn  which is capable to compute the non-Gaussianity parameter
in different triangle configurations.

First, we use the code to calculate the equilateral non-Gaussianity
$f_{\mathrm{NL}}^{\mathrm{equil}}$ for the three cases regarded in this paper.
The results of these cases are presented in Fig. 5. We see in the figure that
at the CMB scales, the equilateral non-Gaussianity for the three cases take
the values which are of order the slow-roll parameters,
$f_{\mathrm{NL}}^{\mathrm{equil}}\sim\mathcal{O}\left(10^{-2}\right)$, and
therefore are in excellent agreement with the 68% CL constraint of Planck 2018
data Planck:2019kim . At the smaller scales, however, as we see in the figure,
the equilateral non-Gaussianity in each case represents a peak in its diagram
which is much larger than unity. For all the three cases the peak of
equilateral non-Gaussianity is of order
$f_{\mathrm{NL}}^{\mathrm{equil}}\sim\mathcal{O}\left(10^{3}\right)$, but the
comoving wavenumber of the peak is different for each case. For Cases 1, 2,
and 3, the comoving wavenumber of the peak is equal to $\mathcal{O}(10^{11})$,
$\mathcal{O}(10^{10})$, and $\mathcal{O}(10^{4})$ $\mathrm{Mpc}^{-1}$,
respectively. Our numerical findings are reported in Table 4. If we compare
these values with the wavenumbers of the peaks in Table 3, we see that the
wavenumbers at which the equilateral non-Gaussianity peaks are close to the
one at which the scalar power spectrum peaks. More precisely, the peaks of the
equilateral non-Gaussianity appear at the scales where the power spectrum of
the scalar perturbations undergoes a sharp declination in its shape before it
enhances and reaches its maximum value. Our conclusion about the large
equilateral non-Gaussianities in the production of PBHs from the models with a
sudden declination is in agreement with the result of Davies:2021loj  for an
ultra-slow-roll inflation model, and with the findings of Zhang:2021rqs  in
the setup of Gauss-Bonnet inflation. The existence of such a sharp declination
appears also in some other models like constant-roll Motohashi:2019rhu  and
ultra constant-roll Liu:2020oqe  inflations. We expect that in these models
the equilateral non-Gaussianity presents a substantial enhancement too, that
can be examined in future investigations. But there are some models, such as
the non-minimal derivative Fu:2019ttf ; dalianis2020generalized ;
Dalianis:2019vit ; Teimoori:2021thk  and $k/G$ Lin:2020goi ; Yi:2020cut ;
Gao:2020tsa ; Gao:2021vxb ; Solbi:2021rse ; Solbi:2021wbo ; Teimoori:2021pte ;
Heydari:2021gea ; Heydari:2021qsr  inflationary scenarios, in which the
variation of the scalar power spectrum is rather mild over all scales, and so
we expect that these models do not show a large enhancement in the equilateral
non-Gaussianity and it remains of order
$f_{\mathrm{NL}}^{\mathrm{equil}}\lesssim 1$ in all scales. This point also
may be verified by further examinations in the future.

Figure 5: The equilateral non-Gaussianity parameter
$f_{\mathrm{NL}}^{\mathrm{equil}}$ for the three Cases 1, 2, and 3. Table 4:
The comoving wavenumber and amount of non-Gaussianity parameter for the peak
of the equilateral non-Gaussianity in Case 1, 2, and 3.

# |  $k\,\left[\mathrm{Mpc}^{-1}\right]$ | $f_{\mathrm{NL}}^{\mathrm{equil}}$
---|---|---
Case 1 | $2.190\times 10^{10}$ | $3577.85$
Case 2 | $2.330\times 10^{6}$ | $2022.86$
Case 3 | $2.999\times 10^{3}$ | $1864.40$

In the next step, we turn to compute the non-Gaussianity in the general form
of scalene triangle configuration in which the triangle sides $k_{1}$,
$k_{2}$, and $k_{3}$ can take different values that satisfy the triangular
inequalities. The parameter $k_{1}$ in these figures is fixed to the comoving
wavenumber at which the equilateral non-Gaussianity peaks and its numerical
values are presented in Table 4. The results of Cases 1, 2, and 3 in this
configuration are presented in Figs. 6, 7, and 8, respectively. The top-left
point in these figures is related to the squeezed configuration ($k_{1}\approx
k_{2}\gg k_{3}$), and the signal of the local non-Gaussianity
$f_{\mathrm{NL}}^{\mathrm{local}}$ is maximized in this configuration
Babich:2004gb ; Bartolo:2004if ; Baumann:2009ds . In other words, if we see
any enhancement at that point in comparison to the other points of the
diagram, then the local non-Gaussianity is the most important form in that
model. The squeezed configuration is the dominant mode in the multi-field
inflationary models Bartolo:2001cw ; Sasaki:2008uc , the curvaton scenario
Linde:1996gt ; Lyth:2002my , the models with inhomogeneous reheating
Dvali:2003em ; Kofman:2003nx , and the New Ekpyrotic frameworks
Creminelli:2007aq ; Koyama:2007if ; Buchbinder:2007at .

The top-right corner in these diagrams corresponds to the equilateral shape
($k_{1}=k_{2}=k_{3}$) Babich:2004gb ; Bartolo:2004if ; Baumann:2009ds , and
this kind of non-Gaussianity arises in the models with higher-derivative
interactions and non-trivial sound speeds Alishahiha:2004eh ; Chen:2006nt ;
Li:2008qc ; Tolley:2009fg ; Amani:2018ueu ; Rasouli:2018kvy . The bottom
middle-point in the diagrams are related to the folded triangle form
($k_{1}=2k_{2}=2k_{3}$), and a signal at this point for a model implies that
the orthogonal non-Gaussianity $f_{\mathrm{NL}}^{\mathrm{ortho}}$ is the
dominant kind in that model Baumann:2009ds . This type of non-Gaussianity
finds importance in the inflation models with non-standard initial states
Chen:2006nt ; Holman:2007na . Additionally, there are the intermediate cases,
namely the elongated ($k_{1}=k_{2}+k_{3}$) and isosceles ($k_{1}>k_{2}=k_{3}$)
that their signatures are pronounced respectively along with the left and
right sides of the triangle in the figures. In our setting for PBHs
production, since the bispectrum enhances at the top-right point, as we see in
Figs. 6, 7, and 8, the equilateral non-Gaussianity is the dominant kind.

From Fig. 5, or also from Figs. 6, 7, and 8, we see that the value of the
equilateral non-Gaussianity in our model can reach very large values. This
enhancement appears at the scales related to BBN Jeong:2014gna ;
Inomata:2016uip , $\mu$-distortion Fixsen:1996nj , or smaller scales. Around
the CMB scales, the non-Gaussianity is very small and so it is in very good
agreement with the Planck 2018 constraints Planck:2019kim . Therefore, we
conclude that in order for our model can explain the formation of PBHs, then
we can observe large equilateral non-Gaussianities in some special scales. In
other words, the prediction of large equilateral non-Gaussianities in those
scales is one important new feature of our model. This is a key feature to
check the viability of the model in light of future observations. If such non-
Gaussianities are observed at the BBN Jeong:2014gna ; Inomata:2016uip  and
$\mu$-distortion Fixsen:1996nj  scales or smaller scales in future
experiments, then our model will be confirmed or will be constrained.
Conversely, if we cannot detect such a signal in the forthcoming data, the
model will be ruled out. In addition, detection or non-detection of these non-
Gaussianities provides us a useful insight to find out whether the PBHs can
source DM or not.

Figure 6: The non-Gaussianity in the general form of scalene triangle
configuration for Case 1. Figure 7: The non-Gaussianity in the general form of
scalene triangle configuration for Case 2. Figure 8: The non-Gaussianity in
the general form of scalene triangle configuration for Case 3.

At the end of this section, it is important to note that in our work, we
estimated the abundance of PBHs by assuming the Gaussian statistics for the
curvature perturbations. But, as we saw above, our computations imply that
this assumption is ruled out around the peak of the scalar power spectrum in
this model. On the other hand, the abundance of PBHs is very sensitive to the
statistics of the curvature perturbations. Therefore, the estimation of the
PBHs abundance in the framework of this model should be reconsidered by using
the formulation based on the non-Gaussian statistics, and this is out of the
scope of the present paper, and we leave it for future studies.

##  VI Secondary gravitational waves

The enhancement of the power spectrum of the curvature perturbations, in
addition to the production significant amount of PBHs at special small scales
during inflation, may also lead to the generation of secondary GWs. These GWs
have gathered considerable attention among the inflationary contexts recently
because they can be tested through the data of several earth-based or space-
based GW observatories. The detection of GWs can provide crucial information
to check the validity of inflationary models. In this section, we examine the
secondary GWs in the $\alpha$-attractor inflation scenario with a tiny bump in
the potential energy of the inflaton and then compare its predictions with the
sensitivity regions of some designed GWs detectors.

Although the spectrum of the secondary GWs depends on the statistics of the
primordial density fluctuations, but unlike the PBHs abundance, the
sensitivity of the spectrum in this case is not very efficient and it is
hardly affected Zhang:2021vak . Therefore, it will be useful here to use the
formalism which is valid for the Gaussian perturbations to evaluate the
induced GWs spectrum to provide some estimations for height and frequency of
the spectrum peak in our model. According to this formalism, the fractional
energy density of the induced GWs in the radiation dominated era is expressed
as Kohri:2018awv ; Lu:2019sti

$\Omega_{GW}(k,\eta)=\frac{1}{6}\left(\frac{k}{aH}\right)^{2}\int_{0}^{\infty}dv\int_{|1-v|}^{|1+v|}du\left(\frac{4v^{2}-\left(1-u^{2}+v^{2}\right)^{2}}{4uv}\right)^{2}\overline{I_{RD}^{2}(u,v,x)}\mathcal{P}_{s}(ku)\mathcal{P}_{s}(kv),$
(29)

where $\eta$ indicates the conformal time, and the time average of the source
terms is given by

$\displaystyle\overline{I_{RD}^{2}(u,v,x\to\infty)}=$
$\displaystyle\frac{1}{2x^{2}}\Bigg{[}\left(\frac{3\pi\left(u^{2}+v^{2}-3\right)^{2}\Theta\left(u+v-\sqrt{3}\right)}{4u^{3}v^{3}}+\frac{T_{c}(u,v,1)}{9}\right)^{2}$
$\displaystyle+\left(\frac{\tilde{T}_{s}(u,v,1)}{9}\right)^{2}\Bigg{]}.$ (30)

The functions $T_{c}$ and $T_{s}$ in the above equation are defined in the
following forms

$\displaystyle T_{c}=$
$\displaystyle-\frac{27}{8u^{3}v^{3}x^{4}}\Bigg{\\{}-48uvx^{2}\cos\left(\frac{ux}{\sqrt{3}}\right)\cos\left(\frac{vx}{\sqrt{3}}\right)\left(3\sin(x)+x\cos(x)\right)+$
$\displaystyle
48\sqrt{3}x^{2}\cos(x)\left(v\sin\left(\frac{ux}{\sqrt{3}}\right)\cos\left(\frac{vx}{\sqrt{3}}\right)+u\cos\left(\frac{ux}{\sqrt{3}}\right)\sin\left(\frac{vx}{\sqrt{3}}\right)\right)+$
$\displaystyle
8\sqrt{3}x\sin(x)\Bigg{[}v\left(18-x^{2}\left(u^{2}-v^{2}+3\right)\right)\sin\left(\frac{ux}{\sqrt{3}}\right)\cos\left(\frac{vx}{\sqrt{3}}\right)+$
$\displaystyle
u\left(18-x^{2}\left(-u^{2}+v^{2}+3\right)\right)\cos\left(\frac{ux}{\sqrt{3}}\right)\sin\left(\frac{vx}{\sqrt{3}}\right)\Bigg{]}+$
$\displaystyle
24x\cos(x)\left(x^{2}\left(-u^{2}-v^{2}+3\right)-6\right)\sin\left(\frac{ux}{\sqrt{3}}\right)\sin\left(\frac{vx}{\sqrt{3}}\right)+$
$\displaystyle
24\sin(x)\left(x^{2}\left(u^{2}+v^{2}+3\right)-18\right)\sin\left(\frac{ux}{\sqrt{3}}\right)\sin\left(\frac{vx}{\sqrt{3}}\right)\Bigg{\\}}$
$\displaystyle-\frac{\left(27\left(u^{2}+v^{2}-3\right)^{2}\right)}{4u^{3}v^{3}}\Bigg{\\{}\text{Si}\left[\left(\frac{u-v}{\sqrt{3}}+1\right)x\right]-\text{Si}\left[\left(\frac{u+v}{\sqrt{3}}+1\right)x\right]$
$\displaystyle+\text{Si}\left[\left(1-\frac{u-v}{\sqrt{3}}\right)x\right]-\text{Si}\left[\left(1-\frac{u+v}{\sqrt{3}}\right)x\right]\Bigg{\\}},$
(31)

$\displaystyle T_{s}=$
$\displaystyle\frac{27}{8u^{3}v^{3}x^{4}}\Bigg{\\{}48uvx^{2}\cos\left(\frac{ux}{\sqrt{3}}\right)\cos\left(\frac{vx}{\sqrt{3}}\right)\left(x\sin(x)-3\cos(x)\right)-$
$\displaystyle
48\sqrt{3}x^{2}\sin(x)\left(v\sin\left(\frac{ux}{\sqrt{3}}\right)\cos\left(\frac{vx}{\sqrt{3}}\right)+u\cos\left(\frac{ux}{\sqrt{3}}\right)\sin\left(\frac{vx}{\sqrt{3}}\right)\right)+$
$\displaystyle
8\sqrt{3}x\cos(x)\Bigg{[}v\left(18-x^{2}\left(u^{2}-v^{2}+3\right)\right)\sin\left(\frac{ux}{\sqrt{3}}\right)\cos\left(\frac{vx}{\sqrt{3}}\right)+$
$\displaystyle
u\left(18-x^{2}\left(-u^{2}+v^{2}+3\right)\right)\cos\left(\frac{ux}{\sqrt{3}}\right)\sin\left(\frac{vx}{\sqrt{3}}\right)\Bigg{]}+$
$\displaystyle
24x\sin(x)\left(6-x^{2}\left(-u^{2}-v^{2}+3\right)\right)\sin\left(\frac{ux}{\sqrt{3}}\right)\sin\left(\frac{vx}{\sqrt{3}}\right)+$
$\displaystyle
24\cos(x)\left(x^{2}\left(u^{2}+v^{2}+3\right)-18\right)\sin\left(\frac{ux}{\sqrt{3}}\right)\sin\left(\frac{vx}{\sqrt{3}}\right)\Bigg{\\}}-\frac{27\left(u^{2}+v^{2}-3\right)}{u^{2}v^{2}}+$
$\displaystyle\frac{\left(27\left(u^{2}+v^{2}-3\right)^{2}\right)}{4u^{3}v^{3}}\Bigg{\\{}-\text{Ci}\left[\left|1-\frac{u+v}{\sqrt{3}}\right|x\right]+\ln\left|\frac{3-(u+v)^{2}}{3-(u-v)^{2}}\right|+$
$\displaystyle\text{Ci}\left[\left(\frac{u-v}{\sqrt{3}}+1\right)x\right]-\text{Ci}\left[\left(\frac{u+v}{\sqrt{3}}+1\right)x\right]+\text{Ci}\left[\left(1-\frac{u-v}{\sqrt{3}}\right)x\right]\Bigg{\\}}.$
(32)

In these equations, the sine-integral $\text{Si}(x)$ and cosine-integral
$\text{Ci}(x)$ functions are defined respectively as follows

$\text{Si}(x)=\int_{0}^{x}\frac{\sin(y)}{y}dy,\qquad\text{Ci}(x)=-\int_{x}^{\infty}\frac{\cos(y)}{y}dy.$
(33)

The relation between the current energy densities of the included GWs and the
corresponding values well after their horizon reentry in the radiation
domination era, is given by the following equation

$\Omega_{GW}\left(k,\eta_{0}\right)=\Omega_{GW}(k,\eta)\frac{\Omega_{r0}}{\Omega_{r}(\eta)},$
(34)

where $\Omega_{r}$ denotes the fractional energy density of radiation, and the
subscript $0$ refers to the present epoch. Here, we take
$\Omega_{r0}h^{2}\simeq 4.2\times 10^{-5}$ Cai:2019bmk ; Fu:2019vqc ;
Fu:2020lob . It is worth mentioning that the conformal time $\eta\gg\eta_{k}$
in Eq. (34) is considered to be earlier than the matter-radiation equality,
and of course late enough so that $\Omega_{GW}(k,\eta)$ be converged to a
constant value during the radiation dominated era.

In Fig. 9, the fractional energy density of the secondary GWs is plotted
versus frequency for the three sets of parameters that are listed in Table 1,
by using Eq. (34). The sensitivity regions of the GWs ground-based
interferometers including aLIGO O2 (dotted orange) and O3 (solid orange) runs
LIGOScientific:2017zlf ; LIGOScientific:2017vwq ; KAGRA:2021kbb , projections
for the HLV (LIGO-Hanford, LIGOLivingston and Virgo) network at design
sensitivity (dashed purple), the mid-scale upgrade of aLIGO (A+ detectors,
solid purple) and the ET (light green) Punturo:2010zz  are also marginalized
in the figure Bagui:2021dqi . Moreover, the sensitivity of the future space-
based interferometer LISA is shown (best experimental design (solid black) and
worst experimental design (dashed black)) LISA:2017pwj . The Pulsar Timing
Arrays (PTA’s) considered here are the EPTA (solid pink) lentati2015european ,
North American Nanohertz Observatory for Gravitational Waves (NANOGrav) (solid
light blue) NANOGrav:2020gpb , International Pulsar Timing Array (IPTA) (solid
light pink), and the Square Kilometer Array (SKA) (dashed red) Moore:2014lga ;
Janssen:2014dka ; Schmitz:2020syl . Fig. 9 clearly illustrates that the peak
amplitude of the spectra is of order $10^{-8}$ for all three cases, but the
peaks occur in different frequencies. For Case 1, the peak appears at the
critical frequency $f_{c}\sim{\cal O}(10^{-3})\,\mathrm{Hz}$, and the result
is located within the joint region of LISA (best and worst) detector. The peak
of Case 2 appears at the critical frequency $f_{c}\sim{\cal
O}(10^{-7})\,\mathrm{Hz}$, and the result can be located within the
sensitivity regions of LISA (best), and SKA. For Case 3, the peak takes place
at the frequency $f_{c}\sim{\cal O}(10^{-10})\,\mathrm{Hz}$, and it can lie
within the sensitivity regions of SKA, IPTA, and EPTA. The spectrum of Case 3
can also be located within the region of the detected NANOGrav signal, and
therefore the induced GWs in our PBHs formation scenario are capable to
explain this observation. We present the exact values of the critical
frequencies and peak heights for the three cases in Table 5. As we see, the
predictions of our model in three cases, lie inside the sensitivity
marginalized joints regions of some designed GWs detectors. It is a remarkable
consequence of our scenario that may be regarded in the future to check the
viability of the model in light of the observational data.

Here, we also estimate the tilt of the GWs spectrum at different ranges of
frequency. Recently, it has been argued that the power spectrum of
$\Omega_{\rm GW0}$ can be expressed in terms of frequency as the power-law
form $\Omega_{\rm GW0}\sim f^{n}$, where $n$ is a constant parameter
Xu:2019bdp ; Fu:2019vqc ; Kuroyanagi:2018csn . In the setup of
$\alpha$-attractor inflation with a tiny bump, we compute the parameter $n$
for the frequency bands $f\ll f_{c}$, $f<f_{c}$, and $f>f_{c}$. We denote the
results of these regions by $n_{f\ll f_{c}}$, $n_{f<f_{c}}$, and
$n_{f>f_{c}}$, respectively. Our numerical findings are presented in Table 5.
The results obtained in the infrared regime $f\ll f_{c}$ verify properly the
analytical relation $\Omega_{\rm GW0}\sim f^{3-2/\ln\left(f_{c}/f\right)}$
derived by Yuan:2019wwo ; Cai:2019cdl .

Figure 9: The present fractional energy density of the secondary GWs in terms
of frequency. The thick solid, thick dashed, and thick dash-dotted black plots
correspond to Cases 1, 2, and 3, respectively. The sensitivity regions of the
GWs ground-based interferometers including aLIGO O2 (dotted orange) and O3
(solid orange) runs LIGOScientific:2017zlf ; LIGOScientific:2017vwq ;
KAGRA:2021kbb , projections for the HLV (LIGO-Hanford, LIGOLivingston and
Virgo) network at design sensitivity (dashed purple), the mid-scale upgrade of
aLIGO (A+ detectors, solid purple) and the ET (light green) Punturo:2010zz
are also demonstrated in the figure Bagui:2021dqi . Moreover, the sensitivity
of the future space-based interferometer LISA is shown (best experimental
design (solid black) and worst experimental design (dashed black))
LISA:2017pwj . The Pulsar Timing Arrays (PTA’s) considered here are the EPTA
(solid pink) lentati2015european , North American Nanohertz Observatory for
Gravitational Waves (NANOGrav) (solid light blue) NANOGrav:2020gpb ,
International Pulsar Timing Array (IPTA) (solid light pink), and the Square
Kilometer Array (SKA) (dashed red) Moore:2014lga ; Janssen:2014dka ;
Schmitz:2020syl . Table 5: The frequencies and heights of the peak of the
spectrum of the current fractional energy density of the secondary GWs for
Cases 1, 2, and 3.

# |  $f_{c}$ |  $\Omega_{\rm GW0}\left(f_{c}\right)$ |  $n_{f\ll f_{c}}$ |  $n_{f<f_{c}}$ |  $n_{f>f_{c}}$
---|---|---|---|---|---
Case 1 | $8.163\times 10^{-3}$ | $2.081\times 10^{-8}$ | $3.010$ | $1.584$ | $-2.841$
Case 2 | $8.807\times 10^{-7}$ | $2.954\times 10^{-8}$ | $3.021$ | $1.582$ | $-1.702$
Case 3 | $7.389\times 10^{-10}$ | $4.832\times 10^{-8}$ | $3.016$ | $1.559$ | $-2.736$

##  VII Conclusions

We studied the process of PBHs formation in the setup of $\alpha$-attractor
inflation with a tiny bump in its potential mishra2020primordial . For this
purpose, we focused on three cases which the values of these cases are listed
in Table 1. We evaluated the scalar power spectrum in terms of the comoving
wavenumber by solving numerically the Mukhanov-Sasaki equation for these cases
and showed that in each case, the power spectrum exhibits a peak on a special
scale that is beyond the CMB scales. The diagrams of the scalar power spectrum
are shown in Fig. 3. The height of the peak in all cases is of order
$10^{-2}$, which is seven orders of magnitudes larger than the one observed in
the CMB scales ($\mathcal{P}_{s*}\approx 2.1\times 10^{-9}$) akrami2020planck
. This enhancement is large enough for the formation of PBHs. Therefore, we
confirmed that the inclusion of a tiny bump in the $\alpha$-attractor
potential can explain successfully the PBHs production from inflation. The
abundance of PBHs for each case in this setting has also been estimated, and
it has been shown that in this model the PBHs can be formed in a wide range of
mass scales. The diagrams of the abundance of the produced PBHs in terms of
the mass are plotted in Fig. 4. We concluded that the formed PBHs in Case 1,
can constitute a significant contribution to the DM observed today. We
compared the results of our model for the curvature power spectrum with the
current observational bounds on this observable, and found that all the three
cases satisfy the constraints properly. The results of these cases are also
compatible with the observational bounds on abundance of PBHs.

Subsequently, we turned to calculate the inflationary non-Gaussianity for the
given three cases. For this purpose, we applied the computational code BINGO
Hazra:2012yn  which calculates the primordial non-Gaussianity in different
triangle configurations. At the CMB scales, the results of the models for
different types of non-Gaussianity are of the order of the inflationary slow-
roll parameter,
$f_{\mathrm{NL}}^{\mathrm{equil}}\sim\mathcal{O}\left(10^{-2}\right)$, and
therefore they are in good agreement with the Planck 2018 constraints
Planck:2019kim . However, away from the CMB scales the equilateral non-
Gaussianity in this model gets amplified and can obtain values much larger
than unity which are of order
$f_{\mathrm{NL}}^{\mathrm{equil}}\sim\mathcal{O}\left(10^{3}\right)$. The
peaks of the equilateral non-Gaussianity for the considered cases appear at
the scales which are related to the BBN Jeong:2014gna ; Inomata:2016uip ,
$\mu$-distortion Fixsen:1996nj  events or smaller scales. We also found out
that the equilateral non-Gaussianity for each case peaks at a scale close to
the one for the peak of the power spectrum of the scalar perturbations. In
fact, the peaks of the equilateral non-Gaussianity take place at the scales
where the scalar power spectrum reveals a sudden declination in its diagram
versus comoving wavenumber before it gets amplified and reaches its maximum
value. This result is in agreement with the findings of Davies:2021loj  for an
ultra-slow-roll inflation model, and with the findings of Zhang:2021rqs  in
the setup of Gauss-Bonnet inflation. The existence of such a sharp declination
appears also in some other models like constant-roll Motohashi:2019rhu  and
ultra constant-roll Liu:2020oqe  inflation models. Due to this property, we
expect that these models result in a substantial enhancement for the
equilateral non-Gaussianity too, and this prediction can be checked by further
investigations. Contrariwise, there are some models, such as the non-minimal
derivative Fu:2019ttf ; Dalianis:2019vit ; Teimoori:2021thk  and $k/G$
Lin:2020goi ; Yi:2020cut ; Gao:2020tsa ; Gao:2021vxb ; Solbi:2021rse ;
Solbi:2021wbo ; Teimoori:2021pte ; Heydari:2021gea ; Heydari:2021qsr
inflationary scenarios, where the variation of the scalar power spectrum is
rather mild in all scales, and thus we expect that these models do not show a
large enhancement in the equilateral non-Gaussianity and it remains of order
$f_{\mathrm{NL}}^{\mathrm{equil}}\lesssim 1$ in all scales. This point also
may be verified by future examinations.

We then used the BINGO code to compute the non-Gaussianity in the general form
of a scalene triangle. The results of this configuration for the considered
three cases are presented in two-dimensional diagrams of Figs. 6, 7, and 8.
From these diagrams, we concluded that in the $\alpha$-attractor inflation
with a tiny bump in the potential, the non-Gaussianity peaks in the
equilateral triangle configuration, and therefore the equilateral non-
Gaussianity is the most prominent type of non-Gaussianity in this scenario.
Our investigation implies that the $\alpha$-attractor model of inflation with
a tiny bump in its potential predicts significant non-Gaussianities at scales
beyond the CMB scales which include the BBN Jeong:2014gna ; Inomata:2016uip
and $\mu$-distortion Fixsen:1996nj  scales or smaller scales. If such non-
Gaussianities are observed in the future measurements, then our model will be
confirmed or will be constrained. Moreover, observation or non-observation of
such non-Gaussianities provides a powerful insight for us to understand
whether the PBHs from inflation can be regarded as the origin of DM or not. We
further conclude that study of non-Gaussianity can discriminate between the
models proposed so far for PBHs production through inflation.

It is important to note that in our investigation, we estimated the abundance
of PBHs for the $\alpha$-attractor scenario with a tiny bump by assuming
Gaussian statistics for the curvature perturbations. However, as we discussed
so far, this scenario manifests strong non-Gaussianities around the peak of
the scalar power spectrum, and on the other hand, the abundance of PBHs is
very sensitive to the statistical properties of the primordial curvature
fluctuations. Therefore, the estimation of the PBHs abundance in this scenario
should be reconsidered by applying the formalism which is valid for the non-
Gaussian perturbations. We leave this issue for further investigations in the
future.

We further studied the secondary GWs induced in the process of PBHs production
in our setup. Although the spectrum of the included GWs also depends on the
statistics of the scalar perturbations, but in this case, unlike the PBHs
abundance, the sensitivity is not high such that the results of the Gaussian
and non-Gaussian treatments are almost close to each other with a good
approximation Zhang:2021vak . To study the secondary GWs in our scenario, we
followed the same procedure which is based on the Gaussian perturbations to
provide some assessments for the height and frequency of the spectrum peak for
the different cases of our model. In particular, we computed the present
fractional energy density ($\Omega_{\rm GW0}$) for the three sets of
parameters listed in Table 1. We showed that for all the three cases, the
spectra of $\Omega_{\rm GW0}$ represent their peaks at different frequencies,
while the height of the peaks is of order $10^{-8}$ for all the cases. For
Case 1, the frequency of the peak is $8.163\times 10^{-3}\mathrm{Hz}$, and it
is located inside the sensitivity region of LISA (best and worst) detector.
For cases 2 and 3, the frequencies of the peaks appeared at $8.807\times
10^{-7}\mathrm{Hz}$ and $7.389\times 10^{-10}\,\mathrm{Hz}$, respectively. The
peak of the spectrum of $\Omega_{\rm GW0}$ for case 2 can lie within the
sensitivity regions of LISA (best), and SKA, and for case 3 within the
sensitivity regions of SKA, IPTA, and EPTA. It is interesting that the
spectrum of Case 3 can enter the observational joint region of the NANOGrav
signal, and therefore the induced GWs in our PBHs formation scenario are
capable to explain this signal. Since the results of our scenario can be
located inside the detectable limits of GWs detectors, it seems that we can
hope to test our model in light of the future observational data.

Additionally, we presented some estimations for the tilt of the secondary GWs
spectrum in our scenario for different frequency bands. Our results imply that
the power spectrum of induced GWs behaves almost as the power-law function
$\Omega_{\rm GW0}\sim f^{n}$ in various intervals of frequency. For each case
of our model, we presented the values of the constant parameter $n$ for
different frequency domains and showed that the results for the infrared
regime $f\ll f_{c}$ are in good agreement with the analytical expression
$\Omega_{\rm GW0}\sim f^{3-2/\ln\left(f_{c}/f\right)}$ obtained by
Yuan:2019wwo ; Cai:2019cdl .

## References

  * (1) S. Bird et al.,  Phys. Rev. Lett. 116, 201301 (2016).
  * (2) M. Sasaki, T. Suyama, T. Tanaka, and S. Yokoyama,  Phys. Rev. Lett. 117, 061101 (2016),  [Erratum: Phys.Rev.Lett. 121, 059901 (2018)].
  * (3) S. Blinnikov, A. Dolgov, N. Porayko, and K. Postnov,  JCAP 11, 036 (2016).
  * (4) B. Carr, F. Kuhnel, and M. Sandstad,  Phys. Rev. D 94, 083504 (2016).
  * (5) S. Clesse and J. García-Bellido,  Phys. Dark Univ. 15, 142.
  * (6) B. Carr, M. Raidal, T. Tenkanen, V. Vaskonen, and H. Veermäe,  Phys. Rev. D 96, 023514 (2017).
  * (7) P. Mróz et al.,  Nature 548, 183 (2017).
  * (8) J. Garcia-Bellido and E. Ruiz Morales,  Phys. Dark Univ. 18, 47 (2017).
  * (9) C. Germani and T. Prokopec,  Phys. Dark Univ. 18, 6 (2017).
  * (10) S. Clesse and J. García-Bellido,  Phys. Dark Univ. 22, 137 (2018).
  * (11) J. M. Ezquiaga, J. Garcia-Bellido, and E. Ruiz Morales,  Phys. Lett. B 776, 345 (2018).
  * (12) H. Di and Y. Gong,  JCAP 07, 007 (2018).
  * (13) G. Ballesteros and M. Taoso,  Physical Review D 97, 023501 (2018).
  * (14) M. Zumalacarregui and U. Seljak,  Phys. Rev. Lett. 121, 141101 (2018).
  * (15) G. Ballesteros, J. Beltran Jimenez, and M. Pieroni,  JCAP 06, 016 (2019).
  * (16) A. Ashoorioon, A. Rostami, and J. T. Firouzjaee,  Journal of High Energy Physics 2021, 1 (2021).
  * (17) A. Ashoorioon, A. Rostami, and J. T. Firouzjaee,  Phys. Rev. D 103, 123512 (2021).
  * (18) Y.-F. Cai, X. Tong, D.-G. Wang, and S.-F. Yan,  Phys. Rev. Lett. 121, 081306 (2018).
  * (19) A. Y. Kamenshchik, A. Tronconi, T. Vardanyan, and G. Venturi,  Phys. Lett. B 791, 201 (2019).
  * (20) H. Niikura, M. Takada, S. Yokoyama, T. Sumi, and S. Masaki,  Phys. Rev. D 99, 083503 (2019).
  * (21) C. Chen and Y.-F. Cai,  JCAP 10, 068 (2019).
  * (22) S. S. Mishra and V. Sahni,  Journal of Cosmology and Astroparticle Physics 2020, 007 (2020).
  * (23) G. Ballesteros and M. Taoso,  Phys. Rev. D 97, 023501 (2018).
  * (24) N. Bhaumik and R. K. Jain,  JCAP 01, 037 (2020).
  * (25) C. Germani and T. Prokopec,  Physics of the dark universe 18, 6 (2017).
  * (26) J. Garcia-Bellido and E. R. Morales,  Physics of the dark universe 18, 47 (2017).
  * (27) H. Motohashi and W. Hu,  Phys. Rev. D 96, 063503 (2017).
  * (28) I. Dalianis, A. Kehagias, and G. Tringas,  Journal of Cosmology and Astroparticle Physics 2019, 037 (2019).
  * (29) R. Mahbub,  Phys. Rev. D 101, 023533 (2020).
  * (30) S. Kawai and J. Kim,  Phys. Rev. D 104, 083545 (2021).
  * (31) H. Motohashi, S. Mukohyama, and M. Oliosi,  JCAP 03, 002 (2020).
  * (32) J. Liu, Z.-K. Guo, and R.-G. Cai,  Phys. Rev. D 101, 083535 (2020).
  * (33) C. Fu, P. Wu, and H. Yu,  Phys. Rev. D 100, 063532 (2019).
  * (34) I. Dalianis, S. Karydas, and E. Papantonopoulos,  Journal of Cosmology and Astroparticle Physics 2020, 040 (2020).
  * (35) I. Dalianis, S. Karydas, and E. Papantonopoulos,  JCAP 06, 040 (2020).
  * (36) Z. Teimoori, K. Rezazadeh, and K. Karami,  Astrophys. J. 915, 118 (2021).
  * (37) J. Lin et al.,  Phys. Rev. D 101, 103515 (2020).
  * (38) Z. Yi, Q. Gao, Y. Gong, and Z.-h. Zhu,  Phys. Rev. D 103, 063534 (2021).
  * (39) Q. Gao, Y. Gong, and Z. Yi,  Nuclear Physics B 969, 115480 (2021).
  * (40) Q. Gao,  Science China Physics, Mechanics & Astronomy 64 (2021).
  * (41) M. Solbi and K. Karami,  Eur. Phys. J. C 81, 884 (2021).
  * (42) M. Solbi and K. Karami,  JCAP 08, 056 (2021).
  * (43) Z. Teimoori, K. Rezazadeh, M. A. Rasheed, and K. Karami,  Journal of Cosmology and Astroparticle Physics 10 (2021).
  * (44) S. Heydari and K. Karami,  Eur. Phys. J. C 82, 83 (2022).
  * (45) S. Heydari and K. Karami,  JCAP 03, 033 (2022).
  * (46) B. Abbott et al.,  Phys. Rev. Lett. 116, 061102 (2016).
  * (47) B. P. Abbott et al.,  Phys. Rev. Lett. 116, 241103 (2016).
  * (48) S. W. Hawking,  Commun. Math. Phys. 43, 199 (1975),  [Erratum: Commun.Math.Phys. 46, 206 (1976)].
  * (49) A. Gould,  The Astrophysical Journal 386, L5 (1992).
  * (50) J. J. Dalcanton, C. R. Canizares, A. Granados, C. C. Steidel, and J. T. Stocke,  The Astrophysical Journal 424, 550 (1994).
  * (51) C. Alcock et al.,  The Astrophysical Journal Letters 550, L169 (2001).
  * (52) R. J. Nemiroff, G. F. Marani, J. P. Norris, and J. T. Bonnell,  Physical review letters 86, 580 (2001).
  * (53) P. Wilkinson et al.,  Physical Review Letters 86, 584 (2001).
  * (54) P. Tisserand et al.,  Astronomy & Astrophysics 469, 387 (2007).
  * (55) B. Carr, K. Kohri, Y. Sendouda, and J. Yokoyama,  Physical Review D 81, 104019 (2010).
  * (56) K. Griest, A. M. Cieplak, and M. J. Lehner,  Physical review letters 111, 181302 (2013).
  * (57) D. M. Jacobs, G. D. Starkman, and B. W. Lynn,  Mon. Not. Roy. Astron. Soc. 450, 3418 (2015).
  * (58) Y. Ali-Haïmoud and M. Kamionkowski,  Phys. Rev. D 95, 043534 (2017).
  * (59) S. Wang, Y.-F. Wang, Q.-G. Huang, and T. G. F. Li,  Phys. Rev. Lett. 120, 191102 (2018).
  * (60) G. Sato-Polito, E. D. Kovetz, and M. Kamionkowski,  Phys. Rev. D 100, 063521 (2019).
  * (61) R. Laha, J. B. Muñoz, and T. R. Slatyer,  Phys. Rev. D 101, 123514 (2020).
  * (62) H. Niikura et al.,  Nature Astron. 3, 524 (2019).
  * (63) A. Katz, J. Kopp, S. Sibiryakov, and W. Xue,  JCAP 12, 005 (2018).
  * (64) R. Laha,  Phys. Rev. Lett. 123, 251101 (2019).
  * (65) B. Dasgupta, R. Laha, and A. Ray,  Phys. Rev. Lett. 125, 101101 (2020).
  * (66) B. Carr, K. Kohri, Y. Sendouda, and J. Yokoyama,  Phys. Rev. D 81, 104019 (2010).
  * (67) P. W. Graham, S. Rajendran, and J. Varela,  Phys. Rev. D 92, 063007 (2015).
  * (68) Y. B. Zel’dovich and I. Novikov,  Soviet Astronomy 10, 602 (1967).
  * (69) S. Hawking,  Mon. Not. Roy. Astron. Soc. 152, 75 (1971).
  * (70) B. J. Carr and S. Hawking,  Mon. Not. Roy. Astron. Soc. 168, 399 (1974).
  * (71) Y. Akrami et al.,  Astronomy & Astrophysics 641, A10 (2020).
  * (72) R. Saito, J. Yokoyama, and R. Nagata,  JCAP 06, 024 (2008).
  * (73) C. T. Byrnes, E. J. Copeland, and A. M. Green,  Phys. Rev. D 86, 043512 (2012).
  * (74) S. Young and C. T. Byrnes,  JCAP 08, 052 (2013).
  * (75) Y. Tada and S. Yokoyama,  Phys. Rev. D 91, 123534 (2015).
  * (76) S. Young and C. T. Byrnes,  JCAP 04, 034 (2015).
  * (77) S. Young, D. Regan, and C. T. Byrnes,  JCAP 02, 029 (2016).
  * (78) J. Garcia-Bellido, M. Peloso, and C. Unal,  JCAP 09, 013 (2017).
  * (79) G. Franciolini, A. Kehagias, S. Matarrese, and A. Riotto,  JCAP 03, 016 (2018).
  * (80) J. M. Ezquiaga and J. García-Bellido,  Journal of Cosmology and Astroparticle Physics 2018, 018 (2018).
  * (81) S. Passaglia, W. Hu, and H. Motohashi,  Phys. Rev. D 99, 043536 (2019).
  * (82) V. Atal and C. Germani,  Phys. Dark Univ. 24, 100275 (2019).
  * (83) K. M. Belotsky et al.,  Eur. Phys. J. C 79, 246 (2019).
  * (84) V. De Luca et al.,  JCAP 07, 048 (2019).
  * (85) C.-M. Yoo, J.-O. Gong, and S. Yokoyama,  JCAP 09, 033 (2019).
  * (86) J. M. Ezquiaga, J. García-Bellido, and V. Vennin,  JCAP 03, 029 (2020).
  * (87) P. D. Serpico, V. Poulin, D. Inman, and K. Kohri,  Phys. Rev. Res. 2, 023204 (2020).
  * (88) M. Braglia, X. Chen, and D. K. Hazra,  JCAP 03, 005 (2021).
  * (89) D. G. Figueroa, S. Raatikainen, S. Rasanen, and E. Tomberg,  Phys. Rev. Lett. 127, 101302 (2021).
  * (90) M. W. Davies, P. Carrilho, and D. J. Mulryne,  JCAP 06, 019 (2022).
  * (91) F. Zhang,  Phys. Rev. D 105, 063539 (2022).
  * (92) R. Kallosh and A. Linde,  JCAP 07, 002 (2013).
  * (93) R. Kallosh, A. Linde, and D. Roest,  JHEP 11, 198 (2013).
  * (94) R. Kallosh, A. Linde, and D. Roest,  JHEP 08, 052 (2014).
  * (95) R. Kallosh, A. Linde, and D. Roest,  JHEP 09, 062 (2014).
  * (96) R. Kallosh and A. Linde,  Phys. Rev. D 91, 083528 (2015).
  * (97) A. Linde,  JCAP 05, 003 (2015).
  * (98) M. Galante, R. Kallosh, A. Linde, and D. Roest,  Phys. Rev. Lett. 114, 141302 (2015).
  * (99) J. J. M. Carrasco, R. Kallosh, and A. Linde,  Phys. Rev. D 92, 063519 (2015).
  * (100) J. J. M. Carrasco, R. Kallosh, and A. Linde,  JHEP 10, 147 (2015).
  * (101) D. Roest and M. Scalisi,  Phys. Rev. D 92, 043525 (2015).
  * (102) M. Scalisi,  Journal of High Energy Physics 2015, 1 (2015).
  * (103) M. Eshaghi, M. Zarei, N. Riazi, and A. Kiasatpour,  Phys. Rev. D 93, 123517 (2016).
  * (104) A. A. Starobinsky,  Adv. Ser. Astrophys. Cosmol. 3, 130 (1987).
  * (105) F. L. Bezrukov and M. Shaposhnikov,  Phys. Lett. B 659, 703 (2008).
  * (106) J. Garcia-Bellido, A. D. Linde, and D. Wands,  Phys. Rev. D 54, 6040 (1996).
  * (107) M. Sasaki, J. Valiviita, and D. Wands,  Phys. Rev. D 74, 103003 (2006).
  * (108) M. Kawasaki, T. Takayama, M. Yamaguchi, and J. Yokoyama,  PhRvD 74, 043525 (2006).
  * (109) T. Kawaguchi, M. Kawasaki, T. Takayama, M. Yamaguchi, and J. Yokoyama,  Mon. Not. Roy. Astron. Soc. 388, 1426 (2008).
  * (110) P. H. Frampton, M. Kawasaki, F. Takahashi, and T. T. Yanagida,  JCAP 04, 023 (2010).
  * (111) D. H. Lyth,  arXiv:1107.1681 (2011).
  * (112) E. Bugaev and P. Klimai,  Phys. Rev. D 85, 103504 (2012).
  * (113) K. Kohri, C.-M. Lin, and T. Matsuda,  Phys. Rev. D 87, 103527 (2013).
  * (114) M. Kawasaki, N. Kitajima, and T. T. Yanagida,  PhRvD 87, 063519 (2013).
  * (115) A. Linde, S. Mooij, and E. Pajer,  Phys. Rev. D 87, 103506 (2013).
  * (116) E. V. Bugaev and P. A. Klimai,  Int. J. Mod. Phys. D 22, 1350034 (2013).
  * (117) S. Clesse and J. García-Bellido,  Phys. Rev. D 92, 023524 (2015).
  * (118) J. Garcia-Bellido, M. Peloso, and C. Unal,  JCAP 12, 031 (2016).
  * (119) M. Kawasaki, A. Kusenko, Y. Tada, and T. T. Yanagida,  PhRvD 94, 083523 (2016).
  * (120) K. Inomata, M. Kawasaki, K. Mukaida, Y. Tada, and T. T. Yanagida,  PhRvD 95, 123510 (2017).
  * (121) V. Domcke, F. Muia, M. Pieroni, and L. T. Witkowski,  JCAP 07, 048 (2017).
  * (122) M. Braglia et al.,  JCAP 08, 001 (2020).
  * (123) W. Ahmed, M. Junaid, and U. Zubair,  arXiv:2109.14838 (2021).
  * (124) D. Polarski and A. A. Starobinsky,  Nucl. Phys. B 385, 623 (1992).
  * (125) A. A. Starobinsky and J. Yokoyama,  Density fluctuations in Brans-Dicke inflation,  in 4th Workshop on General Relativity and Gravitation, page 381, 1994.
  * (126) L. A. Kofman and A. D. Linde,  Nucl. Phys. B 282, 555 (1987).
  * (127) L. A. Kofman and D. Y. Pogosian,  Phys. Lett. B 214, 508 (1988).
  * (128) D. Polarski and A. A. Starobinsky,  Phys. Rev. D 50, 6123 (1994).
  * (129) A. A. Starobinsky, S. Tsujikawa, and J. Yokoyama,  Nucl. Phys. B 610, 383 (2001).
  * (130) S. Matarrese, S. Mollerach, and M. Bruni,  PhRvD 58, 043504 (1998).
  * (131) S. Mollerach, D. Harari, and S. Matarrese,  PhRvD 69, 063002 (2004).
  * (132) K. N. Ananda, C. Clarkson, and D. Wands,  PhRvD 75, 123518 (2007).
  * (133) D. Baumann, P. J. Steinhardt, K. Takahashi, and K. Ichiki,  PhRvD 76, 084019 (2007).
  * (134) R. Saito and J. Yokoyama,  PhRvL 102, 161101 (2009),  [Erratum: PhRvL 107, 069901 (2011)].
  * (135) R. Saito and J. Yokoyama,  PThPh 123, 867 (2010),  [Erratum: PThPh 126, 351–352 (2011)].
  * (136) E. Bugaev and P. Klimai,  PhRvD 81, 023517 (2010).
  * (137) E. Bugaev and P. Klimai,  PhRvD 83, 083521 (2011).
  * (138) L. Alabidi, K. Kohri, M. Sasaki, and Y. Sendouda,  JCAP 09, 017 (2012).
  * (139) T. Nakama, J. Silk, and M. Kamionkowski,  PhRvD 95, 043511 (2017).
  * (140) S. Peirone, N. Frusciante, B. Hu, M. Raveri, and A. Silvestri,  PhRvD 97, 063518 (2018).
  * (141) S.-L. Cheng, W. Lee, and K.-W. Ng,  JCAP 07, 001 (2018).
  * (142) N. Bartolo et al.,  Phys. Rev. Lett. 122, 211301 (2019).
  * (143) R.-g. Cai, S. Pi, and M. Sasaki,  Phys. Rev. Lett. 122, 201101 (2019).
  * (144) Y. Lu, Y. Gong, Z. Yi, and F. Zhang,  JCAP 12, 031 (2019).
  * (145) R. Namba, M. Peloso, M. Shiraishi, L. Sorbo, and C. Unal,  JCAP 01, 041 (2016).
  * (146) E. D. Stewart,  PhLB 391, 34 (1997).
  * (147) M. Drees and E. Erfani,  JCAP 04, 005 (2011).
  * (148) S. Datta and S. Bose,  PhRvD 99, 084001 (2019).
  * (149) S. Kasuya and M. Kawasaki,  PhRvD 80, 023516 (2009).
  * (150) R. Easther, J. Giblin, John T., and E. A. Lim,  PhRvL 99, 221301 (2007).
  * (151) S. Antusch, F. Cefala, and S. Orani,  PhRvL 118, 011303 (2017),  [Erratum: PhRvL 120, 219901 (2018)].
  * (152) J. Liu, Z.-K. Guo, R.-G. Cai, and G. Shiu,  PhRvL 120, 031301 (2018).
  * (153) F. Zhang, J. Lin, and Y. Lu,  Phys. Rev. D 104, 063515 (2021).
  * (154) K. Kohri, T. Nakama, and T. Suyama,  Phys. Rev. D 90, 083514 (2014).
  * (155) J. Chluba, A. L. Erickcek, and I. Ben-Dayan,  Astrophys. J. 758, 76 (2012).
  * (156) K. Inomata, M. Kawasaki, and Y. Tada,  Phys. Rev. D 94, 043527 (2016).
  * (157) L. Lentati et al.,  Monthly Notices of the Royal Astronomical Society 453, 2576 (2015).
  * (158) K. Kohri and T. Terada,  Phys. Rev. D 97, 123532 (2018).
  * (159) B. P. Abbott et al.,  Phys. Rev. Lett. 120, 091101 (2018).
  * (160) B. P. Abbott et al.,  Phys. Rev. Lett. 119, 161101 (2017).
  * (161) C. J. Moore, R. H. Cole, and C. P. L. Berry,  Class. Quant. Grav. 32, 015014 (2015).
  * (162) G. Janssen et al.,  PoS AASKA14, 037 (2015).
  * (163) K. Yagi and N. Seto,  Phys. Rev. D 83, 044011 (2011),  [Erratum: Phys.Rev.D 95, 109901 (2017)].
  * (164) E. Phinney et al.,  Big bang observer mission concept study, 2003.
  * (165) N. Seto, S. Kawamura, and T. Nakamura,  Phys. Rev. Lett. 87, 221103 (2001).
  * (166) M. Punturo et al.,  Class. Quant. Grav. 27, 194002 (2010).
  * (167) P. Amaro-Seoane et al.,  arXiv:1702.00786 (2017).
  * (168) B. P. Abbott et al.,  Living Rev. Rel. 21, 3 (2018).
  * (169) K. Griest, A. M. Cieplak, and M. J. Lehner,  Phys. Rev. Lett. 111, 181302 (2013).
  * (170) P. Tisserand et al.,  Astron. Astrophys. 469, 387 (2007).
  * (171) T. D. Brandt,  Astrophys. J. Lett. 824, L31 (2016).
  * (172) D. Aloni, K. Blum, and R. Flauger,  Journal of Cosmology and Astroparticle Physics 2017, 017 (2017).
  * (173) V. Poulin, P. D. Serpico, F. Calore, S. Clesse, and K. Kohri,  Phys. Rev. D 96, 083524 (2017).
  * (174) B. P. Abbott et al.,  Phys. Rev. X 6, 041015 (2016),  [Erratum: Phys.Rev.X 8, 039903 (2018)].
  * (175) B. P. Abbott et al.,  Class. Quant. Grav. 34, 044001 (2017).
  * (176) R. Abbott et al.,  Phys. Rev. D 104, 022004 (2021).
  * (177) M. F. Alam et al.,  Astrophys. J. Suppl. 252, 4 (2021).
  * (178) K. Schmitz,  JHEP 01, 097 (2021).
  * (179) D. K. Hazra, L. Sriramkumar, and J. Martin,  JCAP 05, 026 (2013).
  * (180) J. M. Maldacena,  JHEP 05, 013 (2003).
  * (181) V. Acquaviva, N. Bartolo, S. Matarrese, and A. Riotto,  Nucl. Phys. B 667, 119 (2003).
  * (182) G. I. Rigopoulos, E. P. S. Shellard, and B. J. W. van Tent,  Phys. Rev. D 72, 083507 (2005).
  * (183) Y. Akrami et al.,  Astron. Astrophys. 641, A9 (2020).
  * (184) B. J. Carr,  Astrophys. J. 201, 1 (1975).
  * (185) S. Young, C. T. Byrnes, and M. Sasaki,  JCAP 07, 045 (2014).
  * (186) O. Özsoy, S. Parameswaran, G. Tasinato, and I. Zavala,  Journal of Cosmology and Astroparticle Physics 2018, 005 (2018).
  * (187) Y. Tada and S. Yokoyama,  Phys. Rev. D 100, 023537 (2019).
  * (188) I. Musco and J. C. Miller,  Class. Quant. Grav. 30, 145009 (2013).
  * (189) T. Harada, C.-M. Yoo, and K. Kohri,  Phys. Rev. D 88, 084051 (2013),  [Erratum: Phys.Rev.D 89, 029903 (2014)].
  * (190) A. Escrivà, C. Germani, and R. K. Sheth,  Phys. Rev. D 101, 044022 (2020).
  * (191) S. Young,  Int. J. Mod. Phys. D 29, 2030002 (2019).
  * (192) K. Ando, K. Inomata, and M. Kawasaki,  Phys. Rev. D 97, 103528 (2018).
  * (193) D. Baumann,  Inflation,  in Theoretical Advanced Study Institute in Elementary Particle Physics: Physics of the Large and the Small, pages 523–686, 2011.
  * (194) J. Garriga and V. F. Mukhanov,  Phys. Lett. B 458, 219 (1999).
  * (195) S. Dodelson and L. Hui,  Phys. Rev. Lett. 91, 131301 (2003).
  * (196) A. R. Liddle and S. M. Leach,  Phys. Rev. D 68, 103503 (2003).
  * (197) K. Inomata and T. Nakama,  Phys. Rev. D 99, 043511 (2019).
  * (198) P. Hunt and S. Sarkar,  Journal of Cosmology and Astroparticle Physics 2015, 052 (2015).
  * (199) Z.-C. Chen and Q.-G. Huang,  JCAP 08, 039 (2020).
  * (200) D. Babich, P. Creminelli, and M. Zaldarriaga,  JCAP 08, 009 (2004).
  * (201) N. Bartolo, E. Komatsu, S. Matarrese, and A. Riotto,  Phys. Rept. 402, 103 (2004).
  * (202) N. Bartolo, S. Matarrese, and A. Riotto,  Phys. Rev. D 65, 103505 (2002).
  * (203) M. Sasaki,  Prog. Theor. Phys. 120, 159 (2008).
  * (204) A. D. Linde and V. F. Mukhanov,  Phys. Rev. D 56, R535 (1997).
  * (205) D. H. Lyth, C. Ungarelli, and D. Wands,  Phys. Rev. D 67, 023503 (2003).
  * (206) G. Dvali, A. Gruzinov, and M. Zaldarriaga,  Phys. Rev. D 69, 023505 (2004).
  * (207) L. Kofman,  arXiv:astro-ph/0303614 (2003).
  * (208) P. Creminelli and L. Senatore,  JCAP 11, 010 (2007).
  * (209) K. Koyama, S. Mizuno, F. Vernizzi, and D. Wands,  JCAP 11, 024 (2007).
  * (210) E. I. Buchbinder, J. Khoury, and B. A. Ovrut,  Phys. Rev. Lett. 100, 171302 (2008).
  * (211) M. Alishahiha, E. Silverstein, and D. Tong,  Phys. Rev. D 70, 123505 (2004).
  * (212) X. Chen, M.-x. Huang, S. Kachru, and G. Shiu,  JCAP 01, 002 (2007).
  * (213) M. Li, T. Wang, and Y. Wang,  JCAP 03, 028 (2008).
  * (214) A. J. Tolley and M. Wyman,  Phys. Rev. D 81, 043502 (2010).
  * (215) R. Amani, K. Rezazadeh, A. Abdolmaleki, and K. Karami,  Astrophys. J. 853, 188 (2018).
  * (216) S. Rasouli, K. Rezazadeh, A. Abdolmaleki, and K. Karami,  Eur. Phys. J. C 79, 79 (2019).
  * (217) R. Holman and A. J. Tolley,  JCAP 05, 001 (2008).
  * (218) D. Jeong, J. Pradler, J. Chluba, and M. Kamionkowski,  Phys. Rev. Lett. 113, 061301 (2014).
  * (219) D. Fixsen et al.,  Astrophys. J. 473, 576 (1996).
  * (220) R.-G. Cai, Z.-K. Guo, J. Liu, L. Liu, and X.-Y. Yang,  JCAP 06, 013 (2020).
  * (221) C. Fu, P. Wu, and H. Yu,  Phys. Rev. D 101, 023529 (2020).
  * (222) C. Fu, P. Wu, and H. Yu,  PhRvD 102, 043527 (2020).
  * (223) E. Bagui and S. Clesse,  arXiv:2110.07487 (2021).
  * (224) W.-T. Xu, J. Liu, T.-J. Gao, and Z.-K. Guo,  Phys. Rev. D 101, 023505 (2020).
  * (225) S. Kuroyanagi, T. Chiba, and T. Takahashi,  JCAP 11, 038 (2018).
  * (226) C. Yuan, Z.-C. Chen, and Q.-G. Huang,  Phys. Rev. D 101, 043019 (2020).
  * (227) R.-G. Cai, S. Pi, and M. Sasaki,  Phys. Rev. D 102, 083528 (2020).