# Primordial black holes formation in the inflationary model with field-
dependent kinetic term for quartic and natural potentials

Milad<EMAIL_ADDRESS>and Kayoomars<EMAIL_ADDRESS>Department of Physics, University of Kurdistan, Pasdaran Street, P.O. Box
66177-15175, Sanandaj, Iran

###### Abstract

Within the framework of inflationary model with field-dependent kinetic term
for quartic and natural potentials, we investigate generation of the
primordial black holes (PBHs) and induced gravitational waves (GWs). In this
setup, we consider a kinetic function as
$G(\phi)=g_{I}(\phi)\big{(}1+g_{II}(\phi)\big{)}$ and show that in the
presence of first term $g_{I}(\phi)$ both quartic and natural potentials, in
contrast to the standard model of inflation, can be consistent, with the 68%
CL of Planck observations. Besides, the second term $g_{II}(\phi)$ can cause a
significant enhancement in the primordial curvature perturbations at the small
scales which results the PBHs formation. For the both potentials, we obtain an
enhancement in the scalar power spectrum at the scales $k\sim 10^{12}~{}\rm
Mpc^{-1}$, $10^{8}~{}\rm Mpc^{-1}$, and $10^{5}~{}\rm Mpc^{-1}$, which causes
PBHs production in mass scales around $10^{-13}M_{\odot}$, $10^{-5}M_{\odot}$,
and $10M_{\odot}$, respectively. Observational constraints confirm that PBHs
with a mass scale of $10^{-13}M_{\odot}$ can constitute the total of dark
matter in the universe. Furthermore, we estimate the energy density parameter
of induced GWs which can be examined by the observation. Also we conclude that
it can be parametrized as a power-law function $\Omega_{\rm
GW}\sim(f/f_{c})^{n}$, where the power index equals $n=3-2/\ln(f_{c}/f)$ in
the infrared limit $f\ll f_{c}$.

##  I Introduction

In the early universe, primordial curvature perturbations can generate
overdense regions. The gravitational collapse of the overdense areas may form
primordial black holes (PBHs) after the horizon re-entry Hawking:1971 ;
Carr:1974 . In contrast to the stellar black holes (BHs), PBHs have a broad
mass range, and so they can explain the binary BH detected by the LIGO and
Virgo collaborations Abbott:2016-a ; Abbott:2016-b ; Abbott:2017-a ;
Abbott:2017-b ; Abbott:2017-c . The PBHs formed in the early universe can be
considered as a candidate for dark matter (DM) Ivanov:1994 ; Khlopov1:2005 ;
Frampton:2010 ; Belotsky:2014 ; Clesse:2015 ; Carr:2016 ; Inomata:2017 .
Sasaki et al. Sasaki:2016  showed that PBHs with a merger rate of around
$12-213~{}{\rm Gpc^{-3}yr^{-1}}$ and a mass scale of about $10M_{\odot}$ can
make up roughly ${\mathcal{O}(1)}\%$ of the total DM. Also, this type of PBHs
is well consistent with the LIGO and Virgo observations Abbott:2017-a .

Recently, the ultra-short timescale microlensing events were detected in the
OGLE data and they prepare an allowed region for PBHs formation. The favored
area of OGLE data shows that the abundance of PBHs with a mass scale of
${\mathcal{O}(10^{-5})}M_{\odot}$ can reach ${\mathcal{O}(10^{-2})}$ OGLE .
The observations tightly constrain the wide mass range of PBHs as DM, but
there is no constraint on the mass scale from
${\mathcal{O}(10^{-13})}M_{\odot}$ to ${\mathcal{O}(10^{-11})}M_{\odot}$
Ali:2017 ; WD ; HSC ; EGG ; femto ; kepler ; EROS ; CMB-a ; CMB-b ; Katz:2018
; Montero:2019 . Therefore, PBHs in this mass range can demonstrate all the DM
in the universe.

The production of PBHs by the collapse of the overdense regions demands an
amplitude of the primordial curvature perturbation $A_{s}$ in order of
$\sim{\mathcal{O}(10^{-2})}$ at small scales sato:2019 . The cosmic microwave
background (CMB) anisotropy measurements show that $A_{s}=2.1\times 10^{-9}$
at the pivot scale $k_{*}=0.05{\rm Mpc}^{-1}$ akrami:2018 . It means that
$A_{s}$ must increase nearly seven orders of magnitude at small scales for PBH
formation. A large abundance of PBHs cannot be generated in slow-roll
inflationary models Motohashi:2017 ; Passaglia:2019 . Therefore, the violation
of the slow-roll conditions is needed for PBHs formation. Recent studies have
suggested various scenarios for PBHs production Cai:2018 ; Ballesteros:2019 ;
Ballesteros:2020a ; Ballesteros:2020b ; Kamenshchik:2019 ; Inomata:2018 ;
Ezquiaga:2018 ; Germani:2017 ; Di:2018 ; Ballesteros:2018 ; Dalianis:2019 ;
chen:2019 ; Ozsoy:2018 ; Atal:2019 ; mishra:2020 ; fu:2019 ; lin:2020 ;
Khlopov:2010 ; Belotsky1:2014 ; Belotsky:2019 ; Braglia:2020 ; Braglia2:2020 ;
shiPi:2018 ; Fumagalli:2020a ; Sypsas:2020 ; Dalianis:2020 . For instance, a
parametric resonance, which is due to the oscillating sound speed square, can
enhance the primordial curvature perturbations Cai:2018 ; chen:2019 . Also, in
the single-field model with a non-canonical kinetic term, the curvature
perturbations may increase if the sound speed becomes zero Ballesteros:2019 ;
Kamenshchik:2019 . Additionally, in the multi-field inflationary models,
turning trajectories may lead to the power spectrum enhance at small scales
and consequently the PBHs can be produced Fumagalli:2020a ; Braglia:2020 ;
Sypsas:2020 . One of the most common approach for PBHs formation is using the
inflationary models which have an inflection point. A single field model with
an inflection point can lead to the violation of the slow-roll regime
Germani:2017 ; Di:2018 ; Ezquiaga:2018 ; Dalianis:2019 . Around the inflection
point, the inflaton velocity decreases significantly, and the amplitude of the
primordial curvature perturbation enhance firmly fu:2019 ; lin:2020 . For
increasing the primordial curvature perturbation, we require fine-tuning of
the model parameters. It should be mentioned that the total number of $e$-fold
should remain between 50 and 60 Passaglia:2019 ; Sasaki:2018 . Moreover, the
models must be consistent with the Planck observations at the CMB scale
akrami:2018 .

Furthermore, the PBHs formation is followed by the generation of the induced
gravitational waves (GWs), when the primordial curvature perturbations enhance
significantly Matarrese:1998 ; Mollerach:2004 ; Saito:2009 ; Garcia:2017 ;
Cai:2019-a ; Cai:2019-b ; Cai:2019-c ; Bartolo:2019-a ; Bartolo:2019-b ;
Wang:2019 ; Fumagalli:2020b ; Domenech:2020a ; Domenech:2020b ; Hajkarim:2019
; Kohri:2018 ; Xu:2020 . In other words, after the horizon re-entry, the
collapse of the overdense regions can generate large metric perturbations
besides PBHs. In the second-order, the scalar and tensor perturbations may be
coupled to each other. The scalar metric perturbations, through the second-
order effect, can generate the stochastic GW background Cai:2019-a ;
Cai:2019-b ; Cai:2019-c ; Bartolo:2019-a ; Bartolo:2019-b ; Wang:2019 ;
Fumagalli:2020b . Thus, induced GW signal detection indicates a novel approach
to search for PBHs.

The Galileon inflation model is one of the most popular models in the
inflation context, which is placed in the subset of Horndeski’s theory
Horndeski:1974 ; kobayashi:2010 ; Burrage:2010 ; Tumurtushaa:2019 ;
teimoori:2018 . In this scenario, the Galileon field acts like a scalar field
which is responsible for inflation. In the Minkowski spacetime, the action of
this specific scalar field is invariant under the Galilean symmetry
$\partial_{\mu}\phi\rightarrow\partial_{\mu}\phi+b_{\mu}$ Nicolis:2009 ;
Deffayet:2009a ; Deffayet:2009b . Also, in the Galileon inflation models, the
CMB anomalies and the decaying of the CMB power spectrum at the largest scales
can be explained when inflaton undergoes the ultra slow-roll phase Hirano:2016
. In addition, in lin:2020  it was shown that in the inflationary model driven
by a suitable Galileon term $G(\phi)$, the scalar power spectrum can be
enhanced and PBHs are produced. In lin:2020 , the authors studied the
possibility of PBH formation with the Galileon term
$G(\phi)=-\frac{1}{2}\frac{d}{\sqrt{\big{(}\frac{\phi-\phi_{c}}{c}\big{)}^{2}+1}}$
for the power-law potential $V(\phi)=\lambda\phi^{p}$, where $p=2/5$.

In this paper, our main goal is to investigate the possibility of PBH
formation in the framework of inflation with field-dependent kinetic term for
quartic and natural potentials. The structure of the paper is as follows. We
review the inflation model with field-dependent kinetic term. in Sec. II. The
mechanism of PBH formation is explained in Sec. III. In Sec. IV, we study the
consequences of reheating to verify whether the primordial curvatures re-enter
the horizon during reheating or after. In Sec. V, we estimate the abundance of
PBHs. The induced GWs are investigated in Sec. VI. Finally, Sec. VII is
dedicated to our conclusions.

##  II Inflation with field-dependent kinetic term

The action of our model is given by lin:2020

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

where $X\equiv-\frac{1}{2}g^{\mu\nu}\phi_{,\mu}\phi_{,\nu}$. In addition, $g$
and $R$ are the determinant of the metric $g_{\mu\nu}$ and Ricci scalar,
respectively. Also, $G$ are considered as general functions of the scalar
field $\phi$.

From the action (1), for a spatially flat Friedmann-Robertson-Walker (FRW)
universe the Friedmann equations are obtained as follows lin:2020 ;
kobayashi:2010 ; Ohashi:2012

$\displaystyle
3H^{2}=\frac{1}{2}\dot{\phi}^{2}+V(\phi)-\dot{\phi}^{2}G(\phi),$ (2)
$\displaystyle
2\dot{H}+3H^{2}+\frac{1}{2}\dot{\phi}^{2}-V(\phi)-\dot{\phi}^{2}G(\phi)=0,$
(3)

where we take $M_{\rm{pl}}=1/\sqrt{8\pi G}=1$.

Using the action (1), the equation of motion governing the scalar field $\phi$
reads

$\displaystyle\ddot{\phi}+3H\dot{\phi}+\frac{V_{,\phi}-\dot{\phi}^{2}G_{,\phi}}{1-2G(\phi)}=0,$
(4)

where $,_{\phi}\equiv{\rm d}/{\rm d}\phi$, and the dot describe the derivative
with respect to the cosmic time.

The quadratic action for curvature perturbation ${\cal R}$ at the first order
approximation is given by kobayashi:2010

$\displaystyle S^{(2)}=\frac{1}{2}\int{\rm d}\tau{\rm
d}^{3}x\tilde{z}^{2}(1-2G_{,\phi})\left[({\cal
R}_{\phi}^{\prime})^{2}-(\vec{\nabla}{\cal R}_{\phi})^{2}\right],$ (5)

where $\tilde{z}=\frac{a\dot{\phi}}{H},$ and the prime indicates the
derivative with respect to the conformal time $\tau$ kobayashi:2010 .

In the Fourier space, the Mukhanov-Sasaki (MS) equation can be calculated by
varying the action (5) with the respect to the curvature perturbation ${\cal
R}$ as follows

$\displaystyle
u^{\prime\prime}_{k}+\left(k^{2}-\frac{z^{\prime\prime}}{z}\right)u_{k}=0,$
(6)

where $z=(1-2G)^{1/2}\tilde{z}$, and $u_{k}=z{\cal R}_{\phi,k}$ lin:2020 ;
kobayashi:2010 . Consequently, the scalar power spectrum can be obtained as

${\cal P}_{\cal R}(k)=(2\pi^{2})^{-1}k^{3}|u_{k}/z|^{2}.$ (7)

The slow-roll parameters here are defined as

$\varepsilon_{1}\equiv-\frac{\dot{H}}{H^{2}},\
\varepsilon_{2}\equiv-\frac{\ddot{\phi}}{H\dot{\phi}},\
\varepsilon_{3}\equiv\frac{G_{,\phi}\dot{\phi}^{2}}{V_{\phi}},$ (8)

where $|\varepsilon_{i}|\ll 1$ for $i=1,2,3$. In the slow-roll conditions, the
background Eqs. (2) and (4) can be turned to

$\displaystyle 3H^{2}\simeq V,$ (9) $\displaystyle
3H\dot{\phi}(1-2G)+V_{\phi}\simeq 0.$ (10)

Under the slow-roll approximation, the power spectrum of curvature
perturbations is given by lin:2020

$P_{\cal
R}=\frac{H^{2}}{8\pi^{2}\varepsilon_{1}}\simeq\frac{V^{3}}{12\pi^{2}V_{\phi}^{2}}(1-2G),\\\
$ (11)

and, the scalar spectral index can be estimated as

$n_{s}-1=\frac{1}{1-2G}\left(2\eta_{V}-6\varepsilon_{V}+\frac{2G_{\phi}}{1-2G}\sqrt{2\varepsilon_{V}}\right),$
(12)

where $\varepsilon_{V}\equiv\frac{1}{2}(V^{\prime}/V)^{2}$ and $\eta_{V}\equiv
V^{\prime\prime}/V$. The recent value of the scalar spectral index measured by
the Planck satellite is $n_{s}=0.9627\pm 0.0060$ (68$\%$ CL, Planck 2018
TT+lowE) akrami:2018 .

The tensor power spectrum is given by lin:2020

$P_{T}=\frac{H^{2}}{2\pi^{2}},$ (13)

and the tensor-to-scalar ratio can be written as

$r\equiv\frac{P_{T}}{P_{\cal R}}=\frac{16X(1-2G)}{H^{2}}.$ (14)

There is an upper bound on the tensor-to-scalar ratio provided by the Planck
observation as $r<0.0654$ (68$\%$ CL, Planck 2018 TT+lowE) akrami:2018 . In
the following section, we solve numerically the background equations (3)-(4)
to obtain evolution of the both Hubble parameter and the scalar field. Then,
with the help of numerical solution of the MS equation (6), we estimate the
exact value of the scalar power spectrum. Note that in our numerical
calculations, we use the slow roll solutions as initial conditions.

##  III PBH formation mechanism

A proper kinetic term can amplify the curvature perturbations at small scales.
Also, the proposed function must cause the model to be consistent with the
Planck measurements at the pivot scale $k_{*}=0.05{\rm Mpc}^{-1}$. To this
aim, we suggest the kinetic function $G(\phi)$ to be parameterized as follows

$G(\phi)=g_{I}(\phi)\big{(}1+g_{II}(\phi)\big{)},$ (15)

where

$g_{I}(\phi)=-\left(\frac{\phi}{M}\right)^{\alpha},$ (16)
$g_{II}(\phi)=\frac{1}{2}\frac{d}{\sqrt{\big{(}\frac{\phi-\phi_{c}}{c}\big{)}^{2}+1}}\,.$
(17)

In the above, $g_{I}(\phi)$ is the base kinetic term, and we utilize it to
sure that the model can satisfy the observational constraints on $n_{s}$ and
$r$. Also $M$ is a constant with the mass dimension and for
$M\rightarrow\infty$, our model recovers the standard inflation. Note that the
term $g_{II}(\phi)$ is responsible for generating the peak in the scalar power
spectrum at $\phi=\phi_{c}$. Also, the value of $g_{II}(\phi)$ is vanishing at
distances far from $\phi=\phi_{c}$. Here, the parameters $\alpha$ and $d$ are
dimensionless constants. Also $\phi_{c}$ and $c$ have the dimension of mass.
Additionally, $d$ and $c$ control the height and width of the scalar power
spectrum at peak position, respectively.

Note that Lin et al. lin:2020  has already been studied a similar model like
(15) in which the selected kinetic term has only one part as
$G(\phi)=-\frac{1}{2}\frac{d}{\sqrt{\big{(}\frac{\phi-\phi_{c}}{c}\big{)}^{2}+1}}$.
Also, they considered the power-law potential $V(\phi)=\lambda\phi^{p}$, where
$p=2/5$.

###  III.1 Quartic inflation with field-dependent kinetic term

Here, we are interested in investigating the possibility of PBH formation in
our setup with quartic potential. The inflationary quartic potential has the
following form

$V(\phi)=\frac{\lambda}{4}\phi^{4},$ (18)

where $\lambda$ is a constant parameter. At the horizon exit, we set $N_{*}=0$
and take $\phi_{*}\simeq 0.765$ at the pivot scale, which is estimated by
slow-roll approximation. Note that in the framework of standard inflation,
i.e. the Einstein gravity, the quartic potential (18) is completely rolled out
by the Planck 2018 data akrami:2018 . This motivates us to examine the quartic
potential in this scenario to check viability of the model in light of the
Planck observations.

In the case of quartic inflation, the model is described by six free
parameters ($\alpha$, $\lambda$, $d$, $c$, $\phi_{c}$, $M$). By setting
suitable values for $M$ and $\alpha$ one can get the observational parameters
$n_{s}$ and $r$ compatible with the Planck 2018 data. Also, the parameter
$\lambda$ is fixed by the amplitude of scalar power spectrum at the CMB scale
$k_{*}=0.05~{}{\rm Mpc}^{-1}$ akrami:2018 . Thus, $d$, $c$, and $\phi_{c}$ are
the remaining free parameters of our model, and they can affect the PBH
production.

Table 1: The parameter sets for PBHs production in quartic inflation model. The parameters $\phi_{c}$ and $c$ are in units of $M_{\rm{pl}}=1$. Here, $M=0.46$ and $\alpha=15$. Sets |  $\phi_{c}$ |  $d$ | $c$
---|---|---|---
$\rm{Case~{}}A_{\rm H}$ |  $0.707$ |  $2.659\times 10^{8}$ |  $2.3\times 10^{-11}$
$\rm{Case~{}}B_{\rm H}$ |  $0.739$ |  $1.415\times 10^{8}$ |  $2.1\times 10^{-11}$
$\rm{Case~{}}C_{\rm H}$ |  $0.753$ |  $1.265\times 10^{8}$ |  $1.8\times 10^{-11}$
Table 2: The values of $n_{s}$, $r$, $k_{\text{peak}}$, ${\cal P}_{\cal R}^{\text{peak}}$, $M_{\text{PBH}}^{\text{peak}}$ and $f_{\text{PBH}}^{\text{peak}}$ for the cases listed in Table 1. Sets |  $n_{s}$ | $r$ |  $k_{\text{peak}}/\text{\rm Mpc}^{-1}$ |  ${\cal P}_{\cal R}^{\text{peak}}$ |  $M_{\text{PBH}}^{\text{peak}}/M_{\odot}$ |  $f_{\text{PBH}}^{\text{peak}}$
---|---|---|---|---|---|---
$\rm{Case~{}}A_{\rm H}$ |  $0.963$ |  $0.057$ |  $2.29\times 10^{12}$ |  $0.035$ |  $4.77\times 10^{-13}$ | $0.91$
$\rm{Case~{}}B_{\rm H}$ |  $0.965$ |  $0.057$ |  $4.52\times 10^{8}$ |  $0.042$ |  $1.82\times 10^{-5}$ | $0.031$
$\rm{Case~{}}C_{\rm H}$ |  $0.972$ |  $0.056$ |  $2.61\times 10^{5}$ |  $0.052$ |  $36.9$ | $0.0017$

The slow-roll regime fails on the small scales where the peak of the power
spectrum appears. It means that, for calculating the exact value of the scalar
power spectrum ${\cal P}_{\cal R}(k)=(2\pi^{2})^{-1}k^{3}|u_{k}/z|^{2}$, we
need to solve the MS equation (6), numerically. We find three sets of
parameters, which are listed in Table 1. Here, we set $M=0.46$ and $\alpha=15$
to keep our model predictions consistent with the Planck measurements for
$n_{s}$ and $r$ akrami:2018 . As shown in Table 2, in the cases $A_{H}$ and
$B_{H}$ the predictions of the model for $n_{s}$ and $r$ take place inside the
68$\%$ CL region of the Planck 2018 TT+lowE data akrami:2018 . In the case
$C_{H}$, the results of $n_{s}$ and $r$, respectively, are compatible with the
95$\%$ and 68$\%$ CL of the Planck 2018 data akrami:2018 . Note that these
results are in contrast with the result of quartic potential in the standard
model of inflation, in which the prediction of the model is completely rolled
out in light of the Planck observations akrami:2018 . The parameter $d$ should
be at least in order of ${\cal O}(10^{8})$ to the power spectrum peak enhances
seven orders of magnitude at peak position. For fine tuning the parameter $c$
which controls the width of the scalar power spectrum at peak position, we
have two restrictions. In one hand, the total number of $e$-folds should
remain between 50 and 60. On the other hand, the quantity $g_{II}(\phi)$
should be negligible away from the peak to the usual slow-roll is guaranteed.
Using these limitations, we set the parameter $c$ as in order of ${\cal
O}(10^{-11})$. The values of the power spectrum and the corresponding PBHs
abundance for these sets are shown in Table 2. The evolution of the scalar
field $\phi$ versus the $e$-fold number $N$ is depicted in Fig. 1(a) for the
parameter set $A_{\rm H}$. The flat region in Fig. 1(a) is because of the
decreasing inflaton velocity in this area where the inflaton undergoes the
ultra slow-roll (USR) phase. The value of $\varepsilon_{1}$ is reduced
critically in the USR stage, as shown in the Fig. 1(b). The intense reduction
can give a substantial accretion in the scalar power spectrum. In Fig. 1(c),
we see that during the USR period, the slow-roll condition
$|\varepsilon_{2}|\ll 1$ is violated.

(a) (b) (c)

Figure 1: Evolutions of (a) the scalar filed $\phi$, (b) the first slow-roll
parameter $\varepsilon_{1}$ and (c) the second slow-roll parameter
$\varepsilon_{2}$ versus the $e$-fold number $N$ for the quartic inflation
model. The dashed line shows the violation of the slow-roll condition. The
auxiliary parameters are given by the parameter set of case $A_{\rm H}$ in
Table 1 and we take $N_{*}=0$ at the horizon exit.

Figure 2 presents the scalar power spectrum for the three parameter sets
tabulated in Table 1. As displayed in Fig. 2, the scalar power spectrum grows
seven orders from $\sim{\cal O}(10^{-9})$ at the CMB scale to $\sim{\cal
O}(10^{-2})$ at the small scales, which is ideally suitable for PBHs
formation. Also, the results of our model for the PBHs production are in good
agreement with the observational constraints like the observations of CMB
$\mu$-distortion, big bang nucleosynthesis (BBN), and pulsar timing array
(PTA) Inomata:2019-a ; Inomata:2016 ; Fixsen:1996 .

Figure 2: The scalar power spectrum ${\cal P}_{\cal R}$ in terms of wavenumber
$k$ for the quartic inflation model. The blue, red and green lines are
corresponding to the cases $A_{\rm H}$, $B_{\rm H}$ and $C_{\rm H}$,
respectively. The CMB observations exclude the light-green shaded area
akrami:2018 . The orange zone shows the $\mu$-distortion of CMB Fixsen:1996 .
The cyan area represents the effect on the ratio between neutron and proton
during the big bang nucleosynthesis (BBN) Inomata:2016 . The yellow region
demonstrates the constraint from the PTA observations Inomata:2019-a .

###  III.2 Natural inflation with field-dependent kinetic term

In this section, we study the PBHs formation in the framework of inflation
with field-dependent kinetic term for the natural potential given by

$V(\phi)=\lambda^{4}\left[1+\cos\left(\frac{\phi}{f}\right)\right],$ (19)

where $\lambda$ and $f$ are constant. Here we fix the parameter $\lambda$ by
using $P_{\cal R}(k_{*})=2.1\times 10^{-9}$ at the CMB scale
$k_{*}=0.05~{}{\rm Mpc}^{-1}$ akrami:2018 . With the help of slow-roll
approximation, we obtain the scalar field at the pivot scale as
$\phi_{*}\simeq 0.23$. Also, we set $N_{*}=0$ at the horizon exit. Note that
in the standard inflation setting, the prediction of the natural potential
(19) is not in very good agreement with the latest observations, because its
result takes place inside the 95$\%$ CL constraints of Planck 2018 TT+lowE
data akrami:2018 . This motivates us to investigate whether the results of the
natural potential in this scenario can be improved in light of the Planck 2018
observations.

For the case of natural potential (19), we have seven free parameters
($\alpha$, $\lambda$, $f$, $d$, $c$, $\phi_{c}$, $M$) for our model. By
choosing the proper values for $M$, $\alpha$, and $f$, one can adjust the
values of $n_{s}$ and $r$ to be compatible with the Planck observations. Also,
the parameter $\lambda$ is fixed by the amplitude of scalar power spectrum at
the pivot scale $k_{*}=0.05~{}{\rm Mpc}^{-1}$ akrami:2018 . As a result, the
free parameters $d$, $c$, and $\phi_{c}$ are only remaining parameters which
may affect the generation of PBHs.

For the natural inflation model we obtain the three parameter sets listed in
Table 3 which can lead to the PBHs generation. By setting $M=1.7$,
$\alpha=-1$, and $f=1$, our model is compatible with the Planck observations
for $n_{s}$ and $r$. In the cases $A_{\rm N}$ and $B_{\rm N}$, the results of
$n_{s}$ and $r$ are compatible with the Planck observations at the $68\%$ CL
akrami:2018 . Also in the case $C_{\rm N}$, we get $n_{s}=0.961$ which
satisfies the 95$\%$ CL constraint, but the value of $r=0.005$ falls in the
$68\%$ CL region of Planck 2018 data akrami:2018 .

Table 4 shows the values of power spectrum and the PBHs abundances
corresponding to the three cases of Table 3. In Fig. 3(a), evolution of the
scalar field as a function of the $e$-fold number is plotted for the case
$A_{\rm N}$ in Table 3. In this figure, the plateau-like region at
$\phi=\phi_{c}$ leads to a severe enhancement of the primordial curvature
perturbations (see Fig. 4). In this region, the slow-roll condition is
violated because $|\varepsilon_{2}|>1$, as shown in Fig. 3(c).

Figure 4 shows the scalar power spectrum for the all cases of Table 3. For the
PBHs formation, one needs the scalar power spectrum to increase seven orders
at the small scales, and Fig. 4 exhibits such growth. Furthermore, Fig. 4
shows that the scalar power spectrum of this model is in good agreement with
the observational bounds. These constraints include the observations of CMB
$\mu$-distortion, big bang nucleosynthesis (BBN), and pulsar timing array
(PTA) Inomata:2019-a ; Inomata:2016 ; Fixsen:1996 .

Table 3: The parameter sets for PBHs production in natural inflation model. The parameters $\phi_{c}$ and $c$ are in units of $M_{\rm{pl}}=1$. Here, $M=1.7$, $\alpha=-1$, and $f=1$. Sets |  $\phi_{c}$ |  $d$ | $c$
---|---|---|---
$\rm{Case~{}}A_{\rm N}$ |  $0.74$ |  $5.34\times 10^{9}$ |  $1.9\times 10^{-11}$
$\rm{Case~{}}B_{\rm N}$ |  $0.41$ |  $2.272\times 10^{9}$ |  $1.3\times 10^{-11}$
$\rm{Case~{}}C_{\rm N}$ |  $0.28$ |  $1.311\times 10^{9}$ |  $1.07\times 10^{-11}$
Table 4: The values of $n_{s}$, $r$, $k_{\text{peak}}$, ${\cal P}_{\cal R}^{\text{peak}}$, $M_{\text{PBH}}^{\text{peak}}$ and $f_{\text{PBH}}^{\text{peak}}$ for the cases listed in Table 3. Sets |  $n_{s}$ | $r$ |  $k_{\text{peak}}/\text{\rm Mpc}^{-1}$ |  ${\cal P}_{\cal R}^{\text{peak}}$ |  $M_{\text{PBH}}^{\text{peak}}/M_{\odot}$ |  $f_{\text{PBH}}^{\text{peak}}$
---|---|---|---|---|---|---
$\rm{Case~{}}A_{\rm N}$ |  $0.964$ |  $0.007$ |  $2.03\times 10^{12}$ |  $0.036$ |  $6.05\times 10^{-13}$ | $0.92$
$\rm{Case~{}}B_{\rm N}$ |  $0.965$ |  $0.006$ |  $4.8\times 10^{8}$ |  $0.044$ |  $1.60\times 10^{-5}$ | $0.035$
$\rm{Case~{}}C_{\rm N}$ |  $0.961$ |  $0.005$ |  $2.68\times 10^{5}$ |  $0.053$ |  $34.8$ | $0.0015$

(a) (b) (c)

Figure 3: Same as Fig. 3, but for the natural inflation model. The auxiliary
parameters are given by the parameter set of case $A_{\rm N}$ in Table 3 and
we take $N_{*}=0$ at the horizon exit.  Figure 4: Same as Fig. 2, but for the
natural inflation model. The blue, red and green lines are corresponding to
the cases $A_{\rm N}$, $B_{\rm N}$ and $C_{\rm N}$, respectively.

##  IV Reheating stage

After inflation, the scalar field starts to oscillate around the minimum value
of the potential. Consequently, inflaton is decayed to the particles of the
Standard Model. This process is known as the reheating, which joins the
supercooled universe at the end of inflation to the thermalized universe at
the radiation dominated (RD) epoch.

The decay rate of the inflaton can influence the reheating era duration. On
the other hand, the reheating temperature $T_{\rm reh}$ can affect the decay
rate of the inflaton. Hence, the large and small reheating temperatures lead
to the short and prolonged reheating eras, respectively. The horizon re-entry
can occur before the RD era if reheating phase has a prolonged period.
Depending on the re-entry horizon era, the mathematical formalisms applied to
determine the mass fraction and the energy density of GWs will be different.
Here, we discuss whether PBHs produce in the RD era or reheating epoch
mahbub:2020 .

In this paper, we use the method introduced in Dalianis:2019  to estimate the
time of horizon re-entry. For a scale $k^{-1}$ that exits the horizon $\Delta
N_{k}$ $e$-fold before the end of inflation, we can write

$\left(\frac{a_{k,\text{re}}}{a_{\text{end}}}\right)^{\frac{1}{2}(1+3w)}=e^{\Delta
N_{k}},$ (20)

where $a_{k,\text{re}}$ denotes the scale factor at the horizon re-entry
moment, and $w$ is the equation of state parameter. The number of $e$-folds
between the end of inflation and horizon re-entry can be defined as

$\tilde{N}_{k}\equiv\ln\left(\frac{a_{k,\text{re}}}{a_{\text{end}}}\right).$
(21)

Consequently, Eqs. (20) and (21) can be related to each other as follows

$\tilde{N}_{k}=\left(\frac{2}{1+3w}\right)\Delta N_{k},$ (22)

where $w>-1/3$. The number of $e$-folds throughout the reheating period is
$\tilde{N}_{\text{reh}}\equiv\ln\left(a_{\text{reh}}/a_{\text{end}}\right)=\left[3(1+w_{\text{reh}})\right]^{-1}\ln\left(\rho_{\text{end}}/\rho_{\text{reh}}\right)$,
in which $a_{reh}$ shows the value of the scale factor at the end of reheating
mahbub:2020 ; Dalianis:2019 . Also,
$\rho_{\text{end}}=3H_{\text{end}}^{2}M_{p}^{2}$ is the energy density at the
end of inflation epoch.

The values of $\tilde{N}_{k}$ and $\tilde{N}_{\text{reh}}$ can specify the
time of horizon re-entry. If $\tilde{N}_{k}>\tilde{N}_{\text{reh}}$, the scale
$k^{-1}$ re-enters during the RD era. On the other hand, the re-entering of
the scale $k^{-1}$ with $\tilde{N}_{k}<\tilde{N}_{\text{reh}}$ occurs during
the reheating stage. It is possible to assume $w_{\rm reh}\simeq 0$ because
the reheating phase can be recognized as an early matter-dominated era. The
amount of observable inflation $\Delta N$ is given by mahbub:2020

$\Delta N\simeq
57.3+\frac{1}{4}\ln\left(\frac{\varepsilon_{*}V_{*}}{\rho_{\text{end}}}\right)-\frac{1}{4}\tilde{N}_{\text{reh}}~{},$
(23)

where $\varepsilon_{*}$ and $V_{*}$ show the values of the first slow-roll
parameter and the potential at the pivot scale, respectively. For a scale such
as $k^{-1}$, which re-enters the horizon just at the end of the reheating era,
we have $\tilde{N}_{k}=\tilde{N}_{\text{reh}}$. Hence, $\Delta N_{\text{k}}$
and $\tilde{N}_{\text{reh}}$ can be associated with $\Delta
N_{\text{k}}=\tilde{N}_{\text{reh}}/2$ in which we used Eq. (22) and assumed
$w_{\text{reh}}\simeq 0$. Therefore, the scale $k^{-1}$ re-enters horizon
after reheating phase if $\Delta N_{\text{k}}>\tilde{N}_{\text{reh}}/2$.

Now for the both models studied in this paper, we can calculate the values of
$\Delta N_{\text{k}}^{\rm peak}$ and $\tilde{N}_{\text{reh}}$, where $\Delta
N_{\text{k}}^{\rm peak}$ describes the number of $e$-folds after the end of
inflation until the scale $k^{-1}_{\text{peak}}$ re-enters the horizon. The
co-moving wavenumber of the scale $k^{-1}$ that re-enters the horizon at the
end of reheating phase, $k^{-1}_{\rm reh}$, is given by $k_{\rm
reh}=e^{-\frac{\tilde{N}_{\rm reh}}{2}}k_{\rm end}$, where $k^{-1}_{\rm end}$
is the scale that exits the horizon at the end of inflation Dalianis:2019 .

Our results in Table 5 indicate that for the all cases in the both models we
have $\Delta N_{\text{k}}^{\rm peak}>\tilde{N}_{\text{reh}}/2$. This means
that in our models, the PBHs formation occurs in the RD era. For instance, in
Figs. 5(a) and 5(b) we plot the scalar power spectrum associated with the
cases $A_{\rm H}$ and $A_{\rm N}$, respectively. The shaded regions describe
the scales that re-enters the horizon during the reheating stage. As shown in
this figure, the duration of the reheating phase is short, and the peak scales
re-enter the horizon after reheating epoch. Hence, In the following section,
we use the mathematical formalisms, which are verified in the RD era to
calculate PBHs abundance and energy density of the induced GWs.

Table 5: The values of $N_{\text{peak}}$, $\Delta N_{\text{k}}^{\rm peak}$ and $\tilde{N}_{\text{reh}}$ for the all cases in Tables 1 and 3. Sets |  $\Delta N$ | $k_{\text{peak}}$ |  $N_{\text{peak}}$ |  $\Delta N_{\text{k}}^{\rm peak}$ |  $\tilde{N}_{\text{reh}}$
---|---|---|---|---|---
$\rm{Case~{}}A_{\rm H}$ |  54.5 |  $2.29\times 10^{12}$ |  31.6 |  16.42 |  6.98
$\rm{Case~{}}B_{\rm H}$ |  54.4 |  $4.52\times 10^{8}$ |  23 |  22.16 |  7.44
$\rm{Case~{}}C_{\rm H}$ |  54.3 |  $2.61\times 10^{5}$ |  15.5 |  26.56 |  7.85
$\rm{Case~{}}A_{\rm N}$ |  54 |  $2.03\times 10^{12}$ |  31.4 |  16.5 |  8.4
$\rm{Case~{}}B_{\rm N}$ |  53.5 |  $4.8\times 10^{8}$ |  23 |  22.04 |  9.2
$\rm{Case~{}}C_{\rm N}$ |  53.5 |  $2.68\times 10^{5}$ |  15.5 |  26.83 |  9.2

(a) (b)

Figure 5: The scalar power spectrum for the cases (a) $A_{\rm H}$ and (b)
$A_{\rm N}$, respectively, in quartic and natural inflation models. The shaded
areas denote the scales that re-enter the horizon during the reheating phase

##  V Abundance of primordial black holes

During the RD era, when primordial curvature perturbations re-enter the
horizon, the gravitational collapse may generate PBHs. The PBH mass at the
production time is given by $M=\gamma M_{\rm H}$, where $M_{\rm H}$ is the
horizon mass and $\gamma=0.2$ is the collapse efficiency parameter
Inomata:2017 ; Sasaki:2018 . The current ratio of the PBH mass to the total DM
reads Sasaki:2018

$f_{\rm PBH}(M_{\rm PBH})=1.68\times
10^{8}\left(\frac{\gamma}{0.2}\right)^{\frac{1}{2}}\left(\frac{g_{*}}{106.75}\right)^{-\frac{1}{4}}\left(\frac{M_{\rm
PBH}}{M_{\odot}}\right)^{-\frac{1}{2}}\beta(M_{\rm PBH}),$ (24)

where $g_{*}\simeq 106.75$ is the effective degrees of freedom at the
formation time of the PBH Motohashi:2017 . The mass fraction of PBH $\beta$ in
Eq. (24) can be estimated by Sasaki:2018 ; young:2014 ; harada:2013 ;
Musco:2013 ; Germani:2019 ; Shibata:1999 ; Polnarev:2007 ; Musco:2009

$\beta(M_{\rm PBH})=\gamma\frac{\sigma_{M_{\rm PBH}}}{\sqrt{2\pi}\delta_{\rm
th}}\exp{\left(-\frac{\delta_{\rm th}^{2}}{2\sigma_{M_{\rm
PBH}}^{2}}\right)},$ (25)

where $\delta_{\rm th}=0.4$ shows the threshold density contrast for PBH
formation. Also, $\sigma_{M_{\rm PBH}}$ is the variance of density contrast at
the comoving horizon scale and it is obtained as young:2014

$\sigma_{k}^{2}=\left(\frac{4}{9}\right)^{2}\int\frac{{\rm
d}q}{q}W^{2}(q/k)(q/k)^{4}P_{\cal R}(q),$ (26)

where $W(x)=\exp{\left(-x^{2}/2\right)}$ is the Gaussian window function. It
is possible to associate the mass of the PBHs and the corresponding wavenumber
as Motohashi:2017 ; mishra:2020 ; Sasaki:2018

$M_{\rm PBH}=1.13\times
10^{15}\left(\frac{\gamma}{0.2}\right)\left(\frac{g_{*}}{106.75}\right)^{-\frac{1}{6}}\left(\frac{k_{\rm
PBH}}{k_{*}}\right)^{-2}M_{\odot}.$ (27)

With the help of Eqs. (24) and (27), one can estimate the abundance and the
mass of PBHs. For the both quartic and natural inflationary models, the
results of $M_{\text{PBH}}^{\text{peak}}$ and $f_{\text{PBH}}^{\text{peak}}$
calculated at the peak scale of the PBH formation are listed in Tables 2 and
4. Also in Fig. 6, we plot variations of $f_{\text{PBH}}$ versus
$M_{\text{PBH}}$ for the both models. In this figure, the shaded areas reveal
the observational bounds on the abundance of PBH. Figure 6 shows that (i) the
PBHs with masses $\sim 10^{-13}M_{\odot}$, $\sim 10^{-5}M_{\odot}$, and $\sim
10M_{\odot}$ which are compatible with the observations can be produced in the
aforementioned models. (ii) Surprisingly enough is that for the cases $A_{\rm
H}$ and $A_{\rm N}$, we obtain $f_{\rm PBH}^{\rm peak}=0.91$ and $f_{\rm
PBH}^{\rm peak}=0.92$, respectively, which can explain most of DM in the
universe. (iii) For the cases ${\rm B_{H}}$ and ${\rm B_{N}}$, the abundance
peak of the PBHs can grow up to $f_{\text{PBH}}^{\text{peak}}=0.031$ and
$f_{\text{PBH}}^{\text{peak}}=0.035$, respectively, for the mass scales
$1.82\times 10^{-5}M_{\odot}$ and $1.60\times 10^{-5}M_{\odot}$. The results
obtained for the cases ${\rm B_{H}}$ and ${\rm B_{N}}$ are located at the
allowed region of microlensing events in the OGLE data. (iv) For the cases
${\rm C_{H}}$ and ${\rm C_{N}}$, the peaks of the PBH abundance are located at
$36.9M_{\odot}$ and $34.8M_{\odot}$ and enhance to
$f_{\text{PBH}}^{\text{peak}}=0.0017$ and
$f_{\text{PBH}}^{\text{peak}}=0.0015$, respectively. This result is compatible
with the upper limit bounds on the LIGO merger rate and consequently can
describe BHs observed by LIGO and Virgo collaboration.

(a) (b)

Figure 6: (a) The PBH abundance $f_{\rm PBH}$ versus the mass $M$ in the
quartic inflation model for the cases $A_{\rm H}$ (blue line), $B_{\rm H}$
(red line), and $C_{\rm H}$ (green line). (b) Same as Fig. 6(a) but in the
natural inflation model for the cases $A_{\rm N}$ (blue line), $B_{\rm N}$
(red line), and $C_{\rm N}$ (green line). The observational constraints on PBH
abundance are shown by the shaded areas. The brown shaded area illustrates the
allowed region of PBH abundance in the OGLE data OGLE . The red dots show the
upper bound on the PBH abundance due to the upper limit on the LIGO event
merger rate Ali:2017 . Other shaded regions present the recent observational
constraints containing accretion constraints from CMB CMB-a ; CMB-b , with
EROS/MACHO (EROS/MACHO) EROS , with the Kepler satellite (Kepler) kepler ,
microlensing events with Subaru HSC (Subaru HSC) HSC , white dwarf explosion
(WD) WD , femtolensing of gamma-ray burst (Femto) femto  and extragalactic
gamma rays from PBH evaporation (EG$\gamma$) EGG .

##  VI Induced Gravitational Waves

The induced GWs can produce together with PBHs at the horizon re-entry of the
primordial density perturbations Matarrese:1998 ; Mollerach:2004 ; Saito:2009
; Garcia:2017 ; Cai:2019-a ; Cai:2019-b ; Cai:2019-c ; Bartolo:2019-a ;
Bartolo:2019-b ; Wang:2019 ; Fumagalli:2020b ; Domenech:2020a ; Domenech:2020b
; Hajkarim:2019 ; Kohri:2018 ; Xu:2020 . The induced GWs can be tested by the
GWs observatory like LISA lisa . To investigate the induced GWs, we start with
the perturbed FRW metric in the conformal Newtonian gauge which is written as
Ananda:2007

$\displaystyle
ds^{2}=a(\eta)^{2}\left[-(1+2\Psi)d\eta^{2}+\left((1-2\Psi)\delta_{ij}+\frac{h_{ij}}{2}\right)dx^{i}dx^{j}\right],$
(28)

where $a$ is the scale factor and $\eta$ represents the conformal time. Also
$\Psi$ denotes the first-order scalar perturbations, and $h_{ij}$ is the
perturbation of the second-order transverse-traceless tensor. After inflation,
the universe should be thermalized in the reheating era. Hence, the inflaton
will be decayed into the light particles, and the universe will be dominated
by radiation, consequently. In the RD era, the effect of the inflation field
has almost ignorable on the cosmic evolution. Accordingly, to study the scalar
induced GWs during the RD era, we can use the standard Einstein equation. So,
the equation of motion for second-order tensor perturbations $h_{ij}$
satisfies Ananda:2007 ; Baumann:2007

$\displaystyle
h_{ij}^{\prime\prime}+2\mathcal{H}h_{ij}^{\prime}-\nabla^{2}h_{ij}=-4\mathcal{T}^{lm}_{ij}S_{lm}\;,$
(29)

where $\mathcal{H}\equiv a^{\prime}/a$ is the conformal Hubble parameter. The
quantity $\mathcal{T}^{lm}_{ij}$ denotes the transverse-traceless projection
operator, and $S_{ij}$, which is the GW source term, is obtained as

$\displaystyle
S_{ij}=4\Psi\partial_{i}\partial_{j}\Psi+2\partial_{i}\Psi\partial_{j}\Psi-\frac{1}{\mathcal{H}^{2}}\partial_{i}(\mathcal{H}\Psi+\Psi^{\prime})\partial_{j}(\mathcal{H}\Psi+\Psi^{\prime})\;.$
(30)

During the RD era, the scalar metric perturbation $\Psi$ takes the form
Baumann:2007

$\displaystyle\Psi_{k}(\eta)=\psi_{k}\frac{9}{(k\eta)^{2}}\left(\frac{\sin(k\eta/\sqrt{3})}{k\eta/\sqrt{3}}-\cos(k\eta/\sqrt{3})\right)\;,$
(31)

where $k$ indicates the comoving wavenumber, and the primordial perturbation
$\psi_{k}$ is given by

$\displaystyle\langle\psi_{\bf k}\psi_{\tilde{\bf
k}}\rangle=\frac{2\pi^{2}}{k^{3}}\left(\frac{4}{9}\mathcal{P}_{\mathcal{R}}(k)\right)\delta(\bf{k}+\tilde{\bf
k})\;.$ (32)

The energy density of the induced GWs in the RD era can be estimated by
Kohri:2018

$\displaystyle\Omega_{\rm{GW}}(\eta_{c},k)=\frac{1}{12}{\displaystyle\int^{\infty}_{0}dv\int^{|1+v|}_{|1-v|}du}\left(\frac{4v^{2}-(1+v^{2}-u^{2})^{2}}{4uv}\right)^{2}\mathcal{P}_{\mathcal{R}}(ku)\mathcal{P}_{\mathcal{R}}(kv)\left(\frac{3}{4u^{3}v^{3}}\right)^{2}(u^{2}+v^{2}-3)^{2}$
$\displaystyle\times\left\\{\left[-4uv+(u^{2}+v^{2}-3)\ln\left|\frac{3-(u+v)^{2}}{3-(u-v)^{2}}\right|\right]^{2}+\pi^{2}(u^{2}+v^{2}-3)^{2}\Theta(v+u-\sqrt{3})\right\\}\;,$
(33)

where $\Theta$ is the Heaviside theta function, and $\eta_{c}$ is the time in
which the growing of $\Omega_{\rm{GW}}$ is stopped. Also the scalar power
spectrum $\mathcal{P}_{\mathcal{R}}$ is estimated by the MS equation (6). The
current energy density of the induced GWs can be obtained as follows
Inomata:2019-a

$\displaystyle\Omega_{\rm{GW},0}h^{2}=0.83\left(\frac{g_{c}}{10.75}\right)^{-1/3}\Omega_{\rm{r},0}h^{2}\Omega_{\rm{GW}}(\eta_{c},k)\;,$
(34)

where $\Omega_{\rm{r},0}h^{2}\simeq 4.2\times 10^{-5}$ is the radiation
density parameter at the present time. Also $g_{c}\simeq 106.75$ indicates the
effective degrees of freedom in the energy density at $\eta_{c}$. The
frequency and comoving wavenumber can be related as follows

$\displaystyle f=1.546\times 10^{-15}\left(\frac{k}{{\rm Mpc}^{-1}}\right){\rm
Hz}.$ (35)

Now, with numerical solving of Eq. (VI) and using (34)-(35) one can estimate
the present energy density of induced GWs. The results of $\Omega_{\rm{GW},0}$
for the quartic and natural models are shown in Fig. 7. The blue, red, and
green curves correspond to the PBHs with the mass in the order of
$10^{-13}M_{\odot}$, $10^{-5}M_{\odot}$, and $10M_{\odot}$, respectively. The
shape and the amplitude of the current energy spectra of the induced GWs are
almost the same for all cases, but the peak frequencies are different as
depicted in Figs. 7(a) and 7(b).

Figure 7 shows that (i) for the cases $A_{\rm H}$ and $A_{\rm N}$ with
$M\sim{\cal O}(10^{-13})M_{\odot}$, the peak frequency of $\Omega_{GW,0}$ is
around the $\rm{mHz}$ band, which means both of them can place in the
sensitive region of space-based observatories like LISA, $\rm{TianQin}$, and
Taiji. (ii) In the cases $B_{\rm H}$ and $B_{\rm N}$ with PBH mass around
$10^{-5}M_{\odot}$, the $\Omega_{GW,0}$ has a peak in $f\sim 10^{-6}\rm{Hz}$
band. (iii) The energy density of induced GWs produced from the cases $C_{\rm
H}$ and $C_{\rm N}$ has a peak in frequencies in the order of
$10^{-10}\rm{Hz}$. Hence, the cases $B_{\rm H}$, $C_{\rm H}$, $B_{\rm N}$, and
$C_{\rm N}$ can be tested by SKA observation.

(a) (b)

Figure 7: (a) The induced GWs energy density parameter $\Omega_{\rm GW{,0}}$
versus the frequency in the quartic inflation model for the cases $A_{\rm H}$
(blue line), $B_{\rm H}$ (red line), and $C_{\rm H}$ (green line). (b) Same as
Fig. 7(a) but in the natural inflation model for the cases $A_{\rm N}$ (blue
line), $B_{\rm N}$ (red line), and $C_{\rm N}$ (green line). The black dashed
line denotes the broken power-law behaviour of $\Omega_{GW}$. The dashed
curves describe the sensitivity of GWs observatories, such as the European PTA
(EPTA) EPTA-a ; EPTA-b ; EPTA-c ; EPTA-d , the Square Kilometer Array (SKA)
ska , the Advanced Laser Interferometer Gravitational Wave Observatory (aLIGO)
ligo-a ; ligo-b , the Laser Interferometer Space Antenna (LISA) lisa ; lisa-a
, Taiji taiji  and TianQin tianqin .

Recent studies confirm that the energy density can be parameterized as
$\Omega_{\rm GW}(f)\sim f^{n}$ Xu:2020 ; Fu:2020 ; Kuro:2018 . In this regard,
for the case $A_{\rm H}$ we estimate $\Omega_{\rm GW}\sim f^{1.78}$ for
$f<f_{c}=3.54\times 10^{-3}{\rm Hz}$, and $\Omega_{\rm GW}\sim f^{-2.69}$ for
$f>f_{c}$. Also, for the case $A_{\rm N}$, our calculations indicate that
$\Omega_{\rm GW}\sim f^{1.68}$ for $f<f_{c}=3.73\times 10^{-3}{\rm Hz}$, and
$\Omega_{\rm GW}\sim f^{-2.69}$ for $f>f_{c}$. In addition, the results in
both cases confirm $\Omega_{\rm GW}\sim f^{3-2/\ln(f_{c}/f)}$ in the infrared
limit $f\ll f_{c}$, which is entirely compatible with the results of Yuan:2020
; shipi:2020 .

##  VII Conclusions

In this paper, we investigated the possibility of PBHs formation in the
context of inflation with field-dependent kinetic term for two types of the
inflationary potentials (i.e. quartic and natural potentials). Utilizing a
proper kinetic term may increase the primordial perturbations to ${\cal
O}(10^{-2})$ at small scales, which is required for PBHs formation. On the
other hand, the value of the $P_{\cal R}(k)$ should be consistent with the
Planck observation at the pivot scale. In our model, we use a kinetic function
that contains two parts as $G(\phi)=g_{I}(\phi)\big{(}1+g_{II}(\phi)\big{)}$.
The first part $g_{I}(\phi)$ guarantees that the model is compatible with the
Planck measurements at the CMB scale. Moreover, the fine-tuning of the
parameters of the second section $g_{II}(\phi)$ can severely increase the
$P_{\cal R}(k)$ at small scales.

We find three sets of parameters for each potential which are listed in Tables
1 and 3. Our calculations demonstrate that the scales re-enter the horizon
after the reheating stage. Therefore, we can apply the mathematical formalism
that is valid in the RD epoch. In the following, we estimated the PBHs
abundance for all cases in Tables 2 and 4. In our model, the noteworthy sets
are $A_{\rm H}$ and $A_{\rm N}$, which their $f_{\rm PBH}$ are in order ${\cal
O}(1)$ and can explain the total DM in the universe. For the cases $B_{\rm H}$
and $B_{\rm N}$, the peak values of $f_{\rm PBH}$ are as ${\cal O}(10^{-2})$,
which locate in the allowed region of ultrashort timescale in OGLE data. Also,
the cases of $C_{\rm H}$ and $C_{\rm N}$ are compatible with the upper limit
of LIGO, and $f_{\rm PBH}$ grows to ${\cal O}(10^{-3})$.

In the last section, we studied the propagation of secondary GWs in our model.
For the cases $A_{\rm H}$ and $A_{\rm N}$, the peaks of $\Omega_{\rm GW}$ are
formed in the mHz frequency range and can be examined by the observations of
LISA, Taiji, and TianQin. Moreover, as shown in the Fig. 7, in the cases
$B_{\rm H}$, $C_{\rm H}$, $B_{\rm N}$, and $C_{\rm N}$ the GWs energy density
parameter can be examined by the SKA observation. In addition, our numerical
calculations expose that the GWs energy density parameter can be parameterized
as a power-law function $\Omega_{\rm GW}(f)\sim f^{n}$. The results show that
$\Omega_{\rm GW}\sim f^{1.78}$ for $f<f_{c}=3.54\times 10^{-3}{\rm Hz}$, and
$\Omega_{\rm GW}\sim f^{-2.69}$ for $f>f_{c}$ in the case of $A_{\rm H}$, and
$\Omega_{\rm GW}\sim f^{1.68}$ for $f<f_{c}=3.73\times 10^{-3}{\rm Hz}$, and
$\Omega_{\rm GW}\sim f^{-2.69}$ for $f>f_{c}$ in the case of $A_{\rm N}$.
Also, in the infrared limit $f\ll f_{c}$, for the both cases, the power index
of the GWs energy density satisfies the relation $n=3-2/\ln(f_{c}/f)$, which
is well consistent with the analytical result obtained in Yuan:2020 ;
shipi:2020 .

### Acknowledgements

The authors thank the referee for his/her valuable comments.

## References

  * (1) S. Hawking, Mon. Not. R. Astron. Soc. 152, 75 (1971).
  * (2) B. J. Carr and S.W. Hawking, Mon. Not. R. Astron. Soc. 168, 399 (1974).
  * (3) B. P. Abbott et al. (LIGO Scientific Collaboration and Virgo Collaboration), Phys. Rev. Lett. 116, 061102 (2016).
  * (4) B. P. Abbott et al. (LIGO Scientific Collaboration and Virgo Collaboration), Phys. Rev. Lett. 116, 241103 (2016).
  * (5) B. P. Abbott et al. (LIGO Scientific Collaboration and Virgo Collaboration), Phys. Rev. Lett. 116, 221101 (2017).
  * (6) B. P. Abbott et al. (LIGO Scientific Collaboration and Virgo Collaboration), Astrophys. J. 851, L35 (2017).
  * (7) B. P. Abbott et al. (LIGO Scientific Collaboration and Virgo Collaboration), Phys. Rev. Lett. 119, 141101 (2017).
  * (8) P. Ivanov, P. Nasselsky and I.D. Novikov, Phys. Rev. D 50, 7173 (1994)
  * (9) M. Yu. Khlopov, S. G. Rubin, and A. S. Sakharov, Astropart. Phys. 23, 265 (2005).
  * (10) P. H. Frampton, M. Kawasaki, F. Takahashi, and T. T. Yanagida, JCAP 1004, 023 (2010).
  * (11) K. M. Belotsky et al. Mod. Phys. Lett. A 29, 1440005 (2014).
  * (12) S. Clesse and J. GarcBellido, Phys. Rev. D 92, 023524 (2015).
  * (13) B. Carr, F. Kuhnel, and M. Sandstad, Phys. Rev. D 94, 083504 (2016).
  * (14) K. Inomata, M. Kawasaki, K. Mukaida, Y. Tada, and T. T. Yanagida, Phys. Rev. D 96, 043504 (2017).
  * (15) M. Sasaki, T. Suyama, T. Tanaka, and S. Yokoyama, Phys. Rev. Lett. 117, 061101 (2016).
  * (16) H. Niikura, M. Takada, S. Yokoyama, T. Sumi, and S. Masaki, Phys. Rev. D 99, 083503 (2019).
  * (17) Y. Ali-Haïmoud, E. D. Kovetz, and M. Kamionkowski, Phys. Rev. D 96, 123523 (2017).
  * (18) Y. Ali-Haïmoud and M. Kamionkowski, Phys. Rev. D 95, 043534 (2017).
  * (19) V. Poulin, P. D. Serpico, F. Calore, S. Clesse, and K. Kohri, Phys. Rev. D 96, 083524 (2017).
  * (20) P. Tisserand et al. (EROS-2 Collaboration), Astron. Astrophys. 469, 387 (2007).
  * (21) K. Griest, A. M. Cieplak, and M. J. Lehner, Phys. Rev. Lett. 111, 181302 (2013).
  * (22) H. Niikura et al., Nat. Astron. 3, 524 (2019).
  * (23) P. W. Graham, S. Rajendran, and J. Varela, Phys. Rev. D 92, 063007 (2015).
  * (24) A. Barnacka, J. F. Glicenstein, and R. Moderski, Phys. Rev. D 86, 043001 (2012).
  * (25) B. J. Carr, K. Kohri, Y. Sendouda, and J. Yokoyama, Phys. Rev. D 81, 104019 (2010)
  * (26) A. Katz, J. Kopp, S. Sibiryakov, and W. Xue, JCAP 12, 005 (2018).
  * (27) P. Montero-Camacho, X. Fang, G. Vasquez, M. Silva, and C. M. Hirata, JCAP 08, 031 (2019).
  * (28) G. Sato-Polito, E. D. Kovetz, and M. Kamionkowski, Phys. Rev. D 100, 063521 (2019).
  * (29) Y. Akrami et al. (Planck Collaboration), A&A 641, A10 (2020).
  * (30) H. Motohashi and W. Hu, Phys. Rev. D 96, 063503 (2017).
  * (31) S. Passaglia, W. Hu, and H. Motohashi, Phys. Rev. D 99, 043536 (2019)
  * (32) M. Yu. Khlopov, Res. Astron. Astrophys. 10, 495 (2010).
  * (33) K. M. Belotsky et al, Mod. Phys. Lett. A, 29, 1440005 (2014).
  * (34) C. Germani and T. Prokopec, Phys. Dark Univ. 18, 6 (2017).
  * (35) H. Di and Y. Gong, JCAP 07, 007 (2018).
  * (36) J. M. Ezquiaga, J. Garc´ıa-Bellido, and E. R. Morales, Phys. Lett. B 776, 345 (2018).
  * (37) K. Inomata, M. Kawasaki, K. Mukaida, and T. T. Yanagida, Phys. Rev. D 97, 043514 (2018).
  * (38) G. Ballesteros and M. Taoso, Phys. Rev. D 97, 023501 (2018).
  * (39) O. Ozsoy, S. Parameswaran, G. Tasinato, and I. Zavala, JCAP 07, 005 (2018).
  * (40) S. Pi, Y. l. Zhang, Q. G. Huang and M. Sasaki, JCAP 05, 042 (2018).
  * (41) Y. F. Cai, X. Tong, D. G. Wang, and S. F. Yan, Phys. Rev. Lett. 121, 081306 (2018).
  * (42) C. Chen and Y. F. Cai, JCAP 10, 068 (2019).
  * (43) G. Ballesteros, J. B. Jim´enez, and M. Pieroni, JCAP 06, 016 (2019).
  * (44) A. Y. Kamenshchik, A. Tronconi, T. Vardanyan, and G. Venturi, Phys. Lett. B 791, 201 (2019).
  * (45) V. Atal, J. Garriga, and A. Marcos-Caballero, JCAP 09 073 (2019).
  * (46) K. M. Belotsky et al, Eur. Phys. J. C, 79, 246 (2019).
  * (47) I. Dalianis, A. Kehagias, and G. Tringas, JCAP 01, 037 (2019).
  * (48) C. Fu, P. Wu, and H. Yu, Phys. Rev. D 100, 063532 (2019).
  * (49) S. S. Mishra and V. Sahni, JCAP 04, 007 (2020).
  * (50) J. Lin, Q. Gao, Y. Gong, Y. Lu, C. Zhang, and F. Zhang, Phys. Rev. D 101, 103515 (2020).
  * (51) G. Ballesteros, J. Rey, M. Taoso, A. Urbano, JCAP 07, 025 (2020).
  * (52) G. Ballesteros, J. Rey, F. Rompineve, JCAP 06, 014 (2020).
  * (53) M. Braglia, D. K. Hazrad, F. Finelli, G. F. Smoot, L. Sriramkumari and A. A. Starobinsky, JCAP 08, 001 (2020).
  * (54) J. Fumagalli, S. Renaux-Petel, J. W. Ronayne, and L. T. Witkowski, arXiv:2004.08369.
  * (55) G. A. Palma, S. Sypsas, and C. Zenteno, Phys. Rev. Lett. 125, 121301 (2020).
  * (56) M. Braglia, X. Chen, and D. K. Hazra, JCAP 03, 005 (2021).
  * (57) I. Dalianis, S. Karydas, and E. Papantonopoulos, JCAP 06, 040 (2020).
  * (58) M. Sasaki, T. Suyama, T. Tanaka, and S. Yokoyama, Class. Quant. Grav. 35, 063001 (2018).
  * (59) S. Matarrese, S. Mollerach, and M. Bruni, Phys. Rev. D 58, 043504 (1998).
  * (60) S. Mollerach, D. Harari, and S. Matarrese, Phys. Rev. D 69, 063002 (2004).
  * (61) R. Saito and J. Yokoyama, Phys. Rev. Lett. 102, 161101 (2009.)
  * (62) J. Garcia-Bellido, M. Peloso, and C. Unal, JCAP 1709, 013 (2017).
  * (63) K. Kohri and T. Terada, Phys. Rev. D 97, 123532 (2018).
  * (64) R. G. Cai, S. Pi, and M. Sasaki, Phys. Rev. Lett. 122, 201101 (2019).
  * (65) R. G. Cai, S. Pi, S. J. Wang, and X. Y. Yang, JCAP 05, 013 (2019).
  * (66) Y. F. Cai, C. Chen, X. Tong, D. G. Wang, and S. F. Yan, Phys. Rev. D 100, 043518 (2019).
  * (67) N. Bartolo, V. De Luca, G. Franciolini, A. Lewis, M. Peloso, and A. Riotto, Phys. Rev. Lett. 122, 211301 (2019).
  * (68) N. Bartolo, V. De Luca, G. Franciolini, M. Peloso, D. Racco, and A. Riotto, Phys. Rev. D 99, 103521 (2019).
  * (69) S. Wang, T. Terada, and K. Kohri, Phys. Rev. D 99, 103531 (2019).
  * (70) J. Fumagalli, S. Renaux-Petel, and L. T. Witkowski, JCAP 08, 030 (2021).
  * (71) F. Hajkarim and J. Schaffner-Bielich, Phys. Rev. D 101, 043522 (2020).
  * (72) W. T. Xu, J. Liu, T. J. Gao, and Z. K. Guo, Phys. Rev. D 101, 023505 (2020).
  * (73) G. Domènech, and M. Sasaki, Int. J. Mod. Phys. D 29, 2050028 (2020).
  * (74) G. Domènech, and M. Sasaki, Phys. Rev. D 103, 063531 (2021).
  * (75) G. W. Horndeski, Int. J. Theor. Phys. 10, 363 (1974).
  * (76) T. Kobayashi, M. Yamaguchi, J. Yokoyama, Phys. Rev. Lett. 105, 231302 (2010).
  * (77) C. Burrage, C. d. Rham, D. Seery, and A. J. Tolley, JCAP 01, 014 (2011).
  * (78) Z. Teimoori and K. Karami, Astrophys. J. 864, 41 (2018).
  * (79) G. Tumurtushaa, Eur. Phys. J. C 79, 920 (2019).
  * (80) A. Nicolis, R. Rattazzi, E. Trincherini, Phys. Rev. D 79, 064036 (2009).
  * (81) C. Deffayet, G. Esposito-Farese, A. Vikman, Phys. Rev. D 79, 084003 (2009).
  * (82) C. Deffayet, S. Deser, G. Esposito-Farese, Phys. Rev. D 80, 064015 (2009)
  * (83) S. Hirano, T. Kobayashi and S. Yokoyama, Phys. Rev. D 94, 103515 (2016).
  * (84) J. Ohashi and S. Tsujikawa, JCAP 10, 035 (2012).
  * (85) D. J. Fixsen, E. S. Cheng, J. M. Gales, J. C. Mather, R. A. Shafer, and E. L. Wright, Astrophys. J. 473, 576 (1996).
  * (86) K. Inomata, M. Kawasaki, and Y. Tada, Phys. Rev. D 94, 043527 (2016).
  * (87) K. Inomata and T. Nakama, Phys. Rev. D 99, 043511 (2019).
  * (88) R. Mahbub, Phys. Rev. D 101, 023533 (2020).
  * (89) M. Shibata and M. Sasaki, Phys. Rev. D 60 084002 (1999).
  * (90) A. G. Polnarev and I. Musco, Class. Quant. Grav. 24, 1405 (2007).
  * (91) I. Musco, J. C. Miller and A. G. Polnarev, Class. Quant. Grav. 26, 235001 (2009).
  * (92) T. Harada, C.-M. Yoo, and K. Kohri, Phys. Rev. D 88, 084051 (2013).
  * (93) I. Musco and J. C. Miller, Class. Quant. Grav. 30, 145009 (2013).
  * (94) S. Young, C. T. Byrnes, and M. Sasaki, JCAP 07, 045 (2014).
  * (95) C. Germani and I. Musco, Phys. Rev. Lett. 122, 141302 (2019)
  * (96) P. Amaro-Seoane et al. (LISA Collaboration), arXiv:1702.00786.
  * (97) K. N. Ananda, C. Clarkson, and D. Wands, Phys. Rev. D 75, 123518 (2007).
  * (98) D. Baumann, P. J. Steinhardt, K. Takahashi, and K. Ichiki, Phys. Rev. D 76, 084019 (2007).
  * (99) R. D. Ferdman et al., Class. Quant. Grav. 27, 084014 (2010).
  * (100) G. Hobbs et al., Class. Quant. Grav. 27, 084013 (2010).
  * (101) M. A. McLaughlin, Class. Quant. Grav. 30, 224008 (2013).
  * (102) G. Hobbs, Class. Quant. Grav. 30, 224007 (2013).
  * (103) C. J. Moore, R. H. Cole, and C. P. L. Berry, Class. Quant. Grav. 32, 015014 (2015).
  * (104) G. M. Harry (LIGO Scientific Collaboration), Class. Quant. Grav. 27, 084006 (2010).
  * (105) J. Aasi et al. (LIGO Scientific Collaboration), Class. Quant. Grav. 32, 074001 (2015).
  * (106) K. Danzmann, Class. Quant. Grav. 14, 1399 (1997).
  * (107) W. R. Hu and Y. L. Wu, Natl. Sci. Rev. 4, 685 (2017).
  * (108) J. Luo et al. (TianQin Collaboration), Class. Quant. Grav. 33, 035010 (2016).
  * (109) C. Fu, P. Wu, and H. Yu, Phys. Rev. D 101, 023529 (2020).
  * (110) S. Kuroyanagi, T. Chiba, and T. Takahashi, JCAP 11, 038 (2018).
  * (111) C. Yuan, Z. C. Chen, and Q. G. Huang, Phys. Rev. D 101, 043019 (2020).
  * (112) R. G. Cai, S. Pi, and M. Sasaki, Phys. Rev. D 102, 083528 (2020).