WUCG-22-11

# Primordial black holes from Higgs inflation with a Gauss-Bonnet coupling

Ryodai<EMAIL_ADDRESS>and Shinji
<EMAIL_ADDRESS>Department of Physics, Waseda University,
3-4-1 Okubo, Shinjuku, Tokyo 169-8555, Japan

###### Abstract

Primordial black holes (PBHs) can be the source for all or a part of today’s
dark matter density. Inflation provides a mechanism for generating the seeds
of PBHs in the presence of a temporal period where the velocity of an inflaton
field $\phi$ rapidly decreases toward 0. We compute the primordial power
spectra of curvature perturbations generated during Gauss-Bonnet (GB)
corrected Higgs inflation in which the inflaton field has not only a
nonminimal coupling to gravity but also a GB coupling. For a scalar-GB
coupling exhibiting a rapid change during inflation, we show that curvature
perturbations are sufficiently enhanced by the appearance of an effective
potential $V_{\rm eff}(\phi)$ containing the structures of plateau-type, bump-
type, and their intermediate type. We find that there are parameter spaces in
which PBHs can constitute all dark matter for these three types of $V_{\rm
eff}(\phi)$. In particular, models with bump- and intermediate-types give rise
to the primordial scalar and tensor power spectra consistent with the recent
Planck data on scales relevant to the observations of cosmic microwave
background. This property is attributed to the fact that the number of
e-foldings $\Delta N_{c}$ acquired around the bump region of $V_{\rm
eff}(\phi)$ can be as small as a few, in contrast to the plateau-type where
$\Delta N_{c}$ typically exceeds the order of 10.

##  I Introduction

If there were over-density regions in the early universe, primordial black
holes (PBHs) may form as a result of the gravitational collapse during the
radiation-dominated epoch Zel’dovich and Novikov (1967); Hawking (1971); Carr
and Hawking (1974). Unlike astrophysical black holes, the PBHs can have a wide
range of masses and can be the source for all or a part of dark matter (DM)
Chapline (1975); Meszaros (1975) (see Refs. Khlopov (2010); Sasaki  _et al._
(2018); Carr and Kuhnel (2020); Green and Kavanagh (2021); Villanueva-Domingo
_et al._ (2021); Carr and Kuhnel (2022); Escrivà  _et al._ (2022); Karam  _et
al._ (2022) for recent reviews). Although PBHs have not been observationally
discovered yet, the detection of gravitational waves from binary black holes
Abbott  _et al._ (2016, 2019, 2021a, 2021b) suggested a possibility that they
may arise from non-stellar origins Bird  _et al._ (2016); Sasaki  _et al._
(2016); Clesse and García-Bellido (2017); Wang  _et al._ (2018). In addition,
PBHs have been also considered as possible seeds of supermassive black holes
in the center of galaxies Bean and Magueijo (2002). Various observations have
given the upper limit to their abundance $f$ as a function of mass $M$. In
particular, the mass window in which all DM can be explained by PBHs exists in
the range $10^{-16}M_{\odot}\lesssim M\lesssim 10^{-11}M_{\odot}$, where
$M_{\odot}$ is a solar mass (see Ref. Carr  _et al._ (2021) for a recent
study).

Inflation can provide a possible framework for generating the seed of PBHs on
scales smaller than the observed cosmic microwave background (CMB) temperature
anisotropies Ivanov  _et al._ (1994); Garcia-Bellido  _et al._ (1996); Bullock
and Primack (1997); Yokoyama (1997, 1998); Kawasaki  _et al._ (1998, 2006);
Kohri  _et al._ (2008); Saito  _et al._ (2008); Bugaev and Klimai (2008);
Alabidi and Kohri (2009); Drees and Erfani (2011, 2012); Martin  _et al._
(2013); Kohri  _et al._ (2013); Kawasaki  _et al._ (2013); Clesse and García-
Bellido (2015); Kawasaki and Tada (2016); Kawasaki  _et al._ (2016); Pi  _et
al._ (2018); Garcia-Bellido and Ruiz Morales (2017); Kannike  _et al._ (2017);
Germani and Prokopec (2017); Ando  _et al._ (2018); Ezquiaga  _et al._ (2018);
Motohashi and Hu (2017); Di and Gong (2018); Ballesteros and Taoso (2018);
Garcia-Bellido  _et al._ (2017); Hertzberg and Yamada (2018); Inomata  _et
al._ (2018); Cai  _et al._ (2018a); Drees and Xu (2021); Atal  _et al._ (2019,
2020); Mishra and Sahni (2020); Cheong  _et al._ (2021); Fu  _et al._ (2019);
Dalianis  _et al._ (2020); Ashoorioon  _et al._ (2021); Lin  _et al._ (2020);
Yi  _et al._ (2021); Palma  _et al._ (2020); Braglia  _et al._ (2020); Kefala
_et al._ (2021); Ballesteros  _et al._ (2020); Aldabergenov  _et al._ (2020,
2021); Inomata  _et al._ (2021, 2022a); Dalianis  _et al._ (2021); Cai  _et
al._ (2022a); Lin  _et al._ (2021); Kawai and Kim (2021a); Zhang (2022); Ahmed
_et al._ (2022); Cai  _et al._ (2022b); Pi and Wang (2022); Cheong  _et al._
(2022); Kawai and Kim (2022). If there is an intermediate stage in which the
velocity $\dot{\phi}$ of an inflaton field rapidly decreases toward 0 during
inflation, it is possible to enhance curvature perturbations at particular
scales. In the presence of an inflection point in the inflaton potential
$V(\phi)$ around which the derivative ${\rm d}V/{\rm d}\phi$ is close to 0,
the field velocity decreases as $\dot{\phi}\propto a^{-3}$ in the ultra-slow-
roll (USR) regime, where $a$ is a scale factor Garcia-Bellido and Ruiz Morales
(2017); Kannike  _et al._ (2017); Germani and Prokopec (2017); Ezquiaga  _et
al._ (2018); Motohashi and Hu (2017); Di and Gong (2018); Ballesteros and
Taoso (2018); Hertzberg and Yamada (2018); Drees and Xu (2021); Ballesteros
_et al._ (2020). In many of these models, we require a tuning of model
parameters to generate a plateau region of the potential. Moreover, if the
number of e-foldings acquired in the USR regime exceeds the order 10, the
scalar spectral index $n_{s}$ on CMB scales tends to be inconsistent with the
value constrained by the Planck data Akrami  _et al._ (2020).

There are also models containing one or more bumps/dips or steps in the
potential Atal  _et al._ (2019, 2020); Mishra and Sahni (2020); Kefala  _et
al._ (2021); Inomata  _et al._ (2021, 2022a); Dalianis  _et al._ (2021); Cai
_et al._ (2022a, b); Pi and Wang (2022). In this case, the scalar field
rapidly loses its kinetic energy around them, resulting in a strong
enhancement of curvature perturbations. It is also known that oscillating
features can appear in the scalar power spectrum especially for the step-type
potential. One advantage of these models is that the number of e-foldings
during transition can be as small as the order 1. This allows a possibility
for the compatibility of models with the observed values of $n_{s}$ and
tensor-to-scalar ratio $r$. There are also multi-field inflationary models
leading to the enhancement of curvature perturbations at particular scales
Yokoyama (1997); Garcia-Bellido  _et al._ (1996); Kawasaki  _et al._ (1998,
2013); Kohri  _et al._ (2013); Ando  _et al._ (2018); Aldabergenov  _et al._
(2020, 2021); Palma  _et al._ (2020); Braglia  _et al._ (2020); Cheong  _et
al._ (2022); Kawai and Kim (2022). In this case, we need to address whether
the presence of entropy perturbations does not contradict with CMB constraints
on the isocurvature mode. The possibility for producing the seed of PBHs
during preheating after inflation was also discussed in Refs. Green and Malik
(2001); Bassett and Tsujikawa (2001); Suyama  _et al._ (2005); Martin  _et
al._ (2020).

One of the advantages of PBHs as DM is that the origin of DM can be explained
within the framework of Standard Model (SM) of particle physics. A minimal
model of inflation without introducing additional scalar degrees of freedom to
those appearing in SM is known as Higgs inflation, in which the Higgs field
$\phi$ is nonminimally coupled to gravity Bezrukov and Shaposhnikov (2008);
Bezrukov  _et al._ (2009, 2011) (see Refs. Futamase and Maeda (1989); Fakir
and Unruh (1990) for early related works). Indeed, this model is perfectly
consistent with observational bounds on $n_{s}$ and $r$ constrained by the CMB
data Komatsu and Futamase (1999); Tsujikawa and Gumjudpai (2004); Linde  _et
al._ (2011); Ade  _et al._ (2014); Tsujikawa  _et al._ (2013). If we allow the
runnings of Higgs self-coupling $\lambda(\phi)$ and nonminimal coupling
$\xi(\phi)$, then it is possible to have an inflection point in the Higgs
potential. This gives rise to a plateau region in which the enhancement of
curvature perturbations occurs to generate the seed of PBHs Ezquiaga  _et al._
(2018); Drees and Xu (2021); Yi  _et al._ (2021); Lin  _et al._ (2021). In
this scenario the typical number of e-foldings acquired during the USR regime
is of order 10, which results in the values of CMB observables deviating from
those in standard Higgs inflation. The model can be consistent with the
current CMB observations, but it is typically outside the $1\sigma$ contour
constrained by the Planck data Ezquiaga  _et al._ (2018).

Recently, Kawai and Kim Kawai and Kim (2021a) proposed a single-field
inflationary scenario in which the inflaton field is coupled to a Gauss-Bonnet
(GB) curvature invariant $R_{\rm GB}^{2}$ of the form $\mu(\phi)R_{\rm
GB}^{2}$. A scalar-field dependent GB coupling $\mu(\phi)$ can give rise to an
inflection point $\phi=\phi_{c}$ in an effective potential of the inflaton. On
the other hand, we have to caution that a large contribution from the scalar-
GB coupling to the inflaton energy density modifies the primordial scalar
power spectrum on CMB scales Hwang and Noh (2005); Guo  _et al._ (2007); Satoh
and Soda (2008); Guo and Schwarz (2009, 2010); Kawai and Kim (2021b).
Moreover, the scalar-GB coupling leads to the propagation speeds of scalar and
tensor perturbations different from the speed of light Kobayashi  _et al._
(2011); Kase and Tsujikawa (2019). Since the Laplacian instabilities
associated with negative squared propagation speeds may arise, we need to make
sure whether the stability conditions are not violated during inflation.

In Ref. Kawai and Kim (2021a), the scalar-GB coupling
$\mu(\phi)=\mu_{0}\tanh[\mu_{1}(\phi-\phi_{c})]$ was proposed to generate the
seed of PBHs around the inflection point $\phi=\phi_{c}$. Since $\mu(\phi)$
approaches constants in the two asymptotic regimes $\phi\ll\phi_{c}$ and
$\phi\gg\phi_{c}$, the scalar-GB coupling is important only in the vicinity of
$\phi=\phi_{c}$. A temporal USR region can arise from the balance between the
GB term and the scalar potential. The enhancement of curvature perturbations
in such a transient epoch was studied for natural inflation Kawai and Kim
(2021a) and for $\alpha$-attractor Zhang (2022). In these papers, the authors
mostly focused on the USR regime realized by a plateau-type effective
potential. In this case, the CMB observables are subject to modifications by
the presence of a plateau-region with the number of e-foldings of order 10.
Hence it is nontrivial to produce the large amplitude of primordial scalar
perturbations responsible for the seed of PBHs, while satisfying CMB
constraints on $n_{s}$ and $r$.

In this paper, we will address this issue in Higgs inflation with a scalar-GB
coupling mentioned above. We will not incorporate the runnings of Higgs and
nonminimal couplings to focus on effects of the scalar-GB coupling on the
background and perturbations. We show that, besides the plateau-type effective
potential, it is possible to realize a bump- or step-type effective potential.
In this latter case, the field velocity $\dot{\phi}$ around $\phi=\phi_{c}$
decreases faster in comparison to the USR regime with a smaller number of
e-foldings of order 1. The primordial scalar power spectrum can also have a
sharp feature with a peak amplitude enhanced by a factor of $10^{7}$. In such
cases, PBHs can be the source for all DM in the mass range
$10^{-16}M_{\odot}\lesssim M\lesssim 10^{-13}M_{\odot}$. Moreover, the bump-
type effective potential can give rise to the values of $n_{s}$ and $r$ inside
the $1\sigma$ observational contour constrained by the Planck CMB data. There
are also intermediate-type potentials between plateau- and bump-types
consistent with the CMB constraints, while generating the seed of PBHs. Thus,
our inflationary scenario provides a versatile possibility for realizing
various shapes of the effective scalar potential. We note that each shape of
potentials was discussed separately in different contexts in the literature.

This paper is organized as follows. In Sec. II, we obtain the background
equations of motion in Higgs inflation with a scalar-GB coupling
$\mu(\phi)R_{\rm GB}^{2}$ and revisit the scalar and tensor power spectra
generated in slow-roll Higgs inflation with $\mu(\phi)=0$. In Sec. III, we
derive an effective potential $V_{\rm eff}(\phi)$ of the inflaton field and
classify it into three classes: (1) plateau-type, (2) bump-type, and (3)
intermediate-type. In Sec. IV, we compute the primordial scalar power spectra
for three sets of model parameters with which there are neither ghost nor
Laplacian instabilities. We show that the bump-type is favored over the
plateau-type for the consistency with CMB observables. In Sec. V, we calculate
the PBH abundance relative to the relic DM density and show that our model
produces a sufficient amount of PBHs that can be the source for all DM. Sec.
VI is devoted to conclusions. Throughout the paper, we use the natural units
($c=\hbar=1$).

##  II Inflationary model with a Gauss-Bonnet term

We begin with theories given by the action

${\cal S}=\int{\rm d}^{4}x\sqrt{-g}\left[\left(\frac{M_{\rm
Pl}^{2}}{2}+\frac{1}{2}\xi\phi^{2}\right)R-\frac{1}{2}g^{\mu\nu}\nabla_{\mu}\phi\nabla_{\nu}\phi-V(\phi)+\mu(\phi)R_{\rm
GB}^{2}\right],$ (1)

where $g$ is a determinant of the metric tensor $g_{\mu\nu}$, $M_{\rm Pl}$ is
the reduced Planck mass, $\xi$ is a nonminimal coupling constant, $\phi$ is a
scalar field with the covariant derivative operator $\nabla_{\mu}$, and $R$ is
the Ricci scalar. The scalar field has a potential of the form

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

where $\lambda$ is a positive coupling constant. The dynamics of nonminimally
coupled inflation with the potential (2) was originally addressed in Refs.
Futamase and Maeda (1989); Fakir and Unruh (1990) (see also Refs. Salopek  _et
al._ (1989); Makino and Sasaki (1991); Kaiser (1995); Komatsu and Futamase
(1999); Tsujikawa and Gumjudpai (2004)). It can also accommodate the Higgs
potential $V(\phi)=\lambda(\phi^{2}-v^{2})^{2}/4$ in the large field regime
$\phi^{2}\gg v^{2}$, where $v={\cal O}(10^{2})$ GeV Bezrukov and Shaposhnikov
(2008); Bezrukov  _et al._ (2009, 2011). Provided that the nonminimal coupling
is in the range

$\xi\gg 1\,,$ (3)

the self-coupling of order $\lambda=0.01\sim 0.1$ can be consistent with the
amplitude of observed CMB temperature anisotropies333If we consider quantum
corrections arising from the renormalization group running of the standard
model, the Higgs self-coupling $\lambda$ can be much smaller than 0.01 or even
negative at inflationary energy scales De Simone  _et al._ (2009); Hamada  _et
al._ (2014, 2015); Bezrukov  _et al._ (2015, 2018). In this paper, we do not
consider the runnings of coupling constants $\lambda$ or $\xi$..

The scalar field is coupled to a GB curvature invariant defined by

$R_{\rm
GB}^{2}\equiv{R^{2}-4R_{\mu\nu}R^{\mu\nu}+R_{\mu\nu\rho\sigma}R^{\mu\nu\rho\sigma}}\,,$
(4)

with a $\phi$-dependent coupling function $\mu(\phi)$, where $R_{\mu\nu}$ and
$R_{\mu\nu\rho\sigma}$ are the Ricci and Riemann tensors, respectively. The
action (1) belongs to a subclass of Horndeski theories with second-order field
equations of motion Horndeski (1974); Kobayashi  _et al._ (2011); Deffayet
_et al._ (2011); Charmousis  _et al._ (2012) (see Appendix. A). In this case,
there is only one propagating scalar degree of freedom besides two tensor
polarizations. As we will study in Sec. III, it is possible to enhance scalar
perturbations at particular scales for a specific choice of $\mu(\phi)$. The
action (1) corresponds to Higgs inflation corrected by the Higgs-GB coupling.
This allows a possibility for generating the seed of PBHs as the source for
all DM within the framework of SM of particle physics.

We note that the $\xi\to 0$ limit in the action (1) with the potential of
natural inflation corresponds to the model studied by Kawai and Kim Kawai and
Kim (2021a). In our model, the basic mechanism for the generation of seeds of
PBHs is similar to that advocated in Ref. Kawai and Kim (2021a). However,
natural inflation is in tension with the observation of CMB temperature
anisotropies Akrami  _et al._ (2020). Instead, we would like to construct an
explicit inflationary model consistent with CMB observations, while enhancing
curvature perturbations on scales relevant to PBHs. As we will show later,
this is indeed possible for Higgs inflation with $\xi\gg 1$ in the presence of
the scalar-GB coupling.

For the background, we consider a spatially flat Friedmann-Lemaître-Robertson-
Walker (FLRW) line element given by

${\rm d}s^{2}=-{\rm d}t^{2}+a^{2}(t)\delta_{ij}{\rm d}x^{i}{\rm d}x^{j}\,,$
(5)

where $a(t)$ is a time-dependent scale factor. On this background, the
Friedmann and scalar-field equations of motion are

$\displaystyle 3\left(M_{\rm
Pl}^{2}+\xi\phi^{2}\right)H^{2}=\frac{1}{2}\dot{\phi}^{2}+V(\phi)-6\left(\xi\phi+4H^{2}\mu_{,\phi}\right)H\dot{\phi}\,,$
(6)
$\displaystyle\ddot{\phi}+3H\dot{\phi}+V_{,\phi}-6\xi\phi\left(2H^{2}+\dot{H}\right)-24\mu_{,\phi}H^{2}\left(H^{2}+\dot{H}\right)=0\,,$
(7)

where a dot represents a derivative with respect to $t$, $H=\dot{a}/a$ is the
Hubble expansion rate, and we use the notations $\mu_{,\phi}={\rm d}\mu/{\rm
d}\phi$ and $V_{,\phi}={\rm d}V/{\rm d}\phi$. Taking the time derivative of
Eq. (6) and using Eq. (7), we obtain the closed-form differential equations

$\displaystyle\dot{H}=\frac{1}{2{\cal
D}}[8\mu_{,\phi}H^{2}(V_{,\phi}-18\xi\phi
H^{2}-24\mu_{,\phi}H^{4})+2\xi\phi(V_{,\phi}-12\xi\phi
H^{2})-\dot{\phi}^{2}(2\xi+1+8\mu_{,\phi\phi}H^{2})$
$\displaystyle\qquad\qquad+8H\dot{\phi}(\xi\phi+4\mu_{,\phi}H^{2})]\,,$ (8)
$\displaystyle\ddot{\phi}+3H\dot{\phi}+V_{{\rm eff},\phi}=0\,,$ (9)

where

$\displaystyle V_{{\rm eff},\phi}$ $\displaystyle\equiv$
$\displaystyle\frac{1}{{\cal D}}[(\xi\phi^{2}+M_{\rm
Pl}^{2})(V_{,\phi}-12\xi\phi
H^{2}-24\mu_{,\phi}H^{4})+3\dot{\phi}^{2}(2\xi+1+8\mu_{,\phi\phi}H^{2})(\xi\phi+4\mu_{,\phi}H^{2})$
(10)
$\displaystyle\quad+8H\dot{\phi}(\mu_{,\phi}V_{,\phi}-3\xi^{2}\phi^{2})-288H^{3}\dot{\phi}\mu_{,\phi}(\xi\phi+2H^{2}\mu_{,\phi})]\,,$
$\displaystyle{\cal D}$ $\displaystyle\equiv$
$\displaystyle(6\xi+1)\xi\phi^{2}+M_{\rm Pl}^{2}+8H\mu_{,\phi}(\dot{\phi}+6\xi
H\phi+12H^{3}\mu_{,\phi})\,.$ (11)

Notice that $V_{\rm eff}$ is an effective potential of the scalar field.
Numerically, we solve Eqs. (8) and (9) for $H$ and $\phi$ with the initial
conditions of $H$, $\phi$, and $\dot{\phi}$ consistent with the Hamiltonian
constraint (6).

###  II.1 Linear perturbations and stability conditions

To study the evolution of cosmological perturbations during inflation, we
consider a perturbed line element containing scalar perturbations $\alpha$,
$\psi$, $\zeta$ and tensor perturbations $h_{ij}$ as

${\rm d}s^{2}=-(1+2\alpha){\rm d}t^{2}+2\partial_{i}\psi{\rm d}t{\rm
d}x^{i}+a^{2}(t)\left[(1+2\zeta)\delta_{ij}+h_{ij}\right]{\rm d}x^{i}{\rm
d}x^{j}\,.$ (12)

In full Horndeski theories including the action (1) as a special case, the
linear perturbation equations of motion were already derived in the literature
Kobayashi  _et al._ (2011). On using the Hamiltonian and momentum constraints
to eliminate $\alpha$ and $\psi$ and integrating the action (1) by parts, the
second-order action of scalar perturbations is given by Kobayashi  _et al._
(2011); De Felice and Tsujikawa (2011a); Kawai and Kim (2021b)

${\cal S}_{s}^{(2)}=\int{\rm d}t{\rm
d}^{3}x\,a^{3}Q_{s}\left[\dot{\zeta}^{2}-\frac{c_{s}^{2}}{a^{2}}(\nabla\zeta)^{2}\right]\,,$
(13)

where

$Q_{s}=16\frac{\Sigma}{\Theta^{2}}Q_{t}^{2}+12Q_{t}\,,\qquad
c_{s}^{2}=\frac{1}{Q_{s}}\left[\frac{16}{a}\frac{{\rm d}}{{\rm
d}t}\left(\frac{a}{\Theta}Q_{t}^{2}\right)-4c_{t}^{2}Q_{t}\right]\,,$ (14)

with

$\displaystyle Q_{t}=\frac{1}{4}\left(8H\mu_{,\phi}\dot{\phi}+M_{\rm
Pl}^{2}+\xi\phi^{2}\right)\,,\qquad c_{t}^{2}=\frac{1}{4Q_{t}}\left(M_{\rm
Pl}^{2}+\xi\phi^{2}+8\mu_{,\phi\phi}\dot{\phi}^{2}+8\mu_{,\phi}\ddot{\phi}\right)\,,$
(15) $\displaystyle\Sigma=\frac{1}{2}\dot{\phi}^{2}-3H^{2}(M_{\rm
Pl}^{2}+\xi\phi^{2})-6\xi
H\phi\dot{\phi}-48H^{3}\mu_{,\phi}\dot{\phi}\,,\qquad\Theta=H(M_{\rm
Pl}^{2}+\xi\phi^{2})+\xi\phi\dot{\phi}+12H^{2}\mu_{\phi}\dot{\phi}\,.$ (16)

In the tensor sector, the reduced action is of the form

${\cal S}_{t}^{(2)}=\frac{1}{2}\int{\rm d}t{\rm d}^{3}x\hskip
2.84544pta^{3}Q_{t}\left[\dot{h}_{ij}^{2}-\frac{c_{t}^{2}}{a^{2}}(\nabla
h_{ij})^{2}\right],$ (17)

where $Q_{t}$ and $c_{t}^{2}$ are defined in Eq. (15). To avoid the ghost and
Laplacian instabilities of scalar and tensor perturbations, we require the
following conditions

$Q_{s}>0\,,\qquad c_{s}^{2}>0\,,\qquad Q_{t}>0\,,\qquad c_{t}^{2}>0\,.$ (18)

On using the background Eq. (6), we can express $Q_{t}$ and $Q_{s}$ in the
forms

$Q_{t}=\frac{1}{12H^{2}}\left[\frac{1}{2}\dot{\phi}^{2}+V(\phi)-6\xi
H\phi\dot{\phi}\right]\,,\qquad
Q_{s}=\frac{4\dot{\phi}^{2}}{\Theta^{2}}Q_{t}\left[2Q_{t}+3(\xi\phi+4H^{2}\mu_{,\phi})^{2}\right]\,.$
(19)

Provided that $\phi$ decreases during inflation in the region $\phi>0$, we
have $-6\xi H\phi\dot{\phi}>0$ for $\xi>0$. In such cases, both $Q_{t}$ and
$Q_{s}$ are positive and hence the ghost instabilities are absent. In the
absence of the scalar-GB coupling, both $c_{t}^{2}$ and $c_{s}^{2}$ are
equivalent to 1. However, the deviations of $c_{t}^{2}$ and $c_{s}^{2}$ from 1
arise in theories with $\mu(\phi)\neq 0$, so we need to numerically compute
$c_{t}^{2}$ and $c_{s}^{2}$ for a given coupling $\mu(\phi)$ to ensure the
absence of Laplacian instabilities.

###  II.2 Higgs slow-roll inflation with $\mu(\phi)=0$

We briefly revisit the background dynamics and perturbation spectra generated
during Higgs slow-roll inflation for $\mu(\phi)=0$. From Eqs. (6) and (9), we
have

$\displaystyle H^{2}=\frac{1}{3(M_{\rm
Pl}^{2}+\xi\phi^{2})}\left(\frac{1}{2}\dot{\phi}^{2}+\frac{1}{4}\lambda\phi^{4}-6\xi
H\phi\dot{\phi}\right)\,,$ (20)
$\displaystyle\ddot{\phi}+3H\dot{\phi}+\frac{(6\xi+1)\xi\phi\dot{\phi}^{2}}{(6\xi+1)\xi\phi^{2}+M_{\rm
Pl}^{2}}+\frac{\lambda\phi^{3}M_{\rm Pl}^{2}}{(6\xi+1)\xi\phi^{2}+M_{\rm
Pl}^{2}}=0\,.$ (21)

Let us consider the large coupling regime with $\xi\gg 1$ and $\xi\phi^{2}\gg
M_{\rm Pl}^{2}$. During slow-roll inflation, the dominant term in Eq. (20) is
the potential $V(\phi)=\lambda\phi^{4}/4$, while the dominant contributions to
Eq. (21) are second and fourth terms. Then, Eqs. (20) and (21) approximately
reduce to

$H^{2}\simeq\frac{\lambda\phi^{2}}{12\xi}\,,\qquad
H\dot{\phi}\simeq-\frac{\lambda\phi M_{\rm Pl}^{2}}{18\xi^{2}}\,.$ (22)

The field value $\phi_{f}$ at the end of inflation is determined by the
condition $\epsilon_{H}\equiv-\dot{H}/H^{2}=1$. Using the two equations in
(22), we have $\epsilon_{H}\simeq 2M_{\rm Pl}^{2}/(3\xi\phi^{2})$ and hence
$\phi_{f}=\sqrt{2/(3\xi)}M_{\rm Pl}$. The number of e-foldings counted
backward from the end of inflation can be estimated as

$N=\int_{\phi}^{\phi_{f}}\frac{H}{\dot{\phi}}{\rm
d}\phi\simeq\frac{3\xi\phi^{2}}{4M_{\rm Pl}^{2}}-\frac{1}{2}\,,$ (23)

where we exploited Eq. (22) in the second approximate equality. For $N\gg 1$,
we obtain the simple relation $\phi^{2}\simeq 4M_{\rm Pl}^{2}N/(3\xi)$.

For perturbations deep inside the Hubble radius (with the wavenumber $k\gg
aH$), they are initially in the Bunch-Davies vacuum state. On the inflationary
background, $\zeta$ and $h_{ij}$ approach constants after the sound horizon
crossing. The power spectra of scalar and tensor perturbations generated
during the quasi de Sitter period are given, respectively, by Kobayashi  _et
al._ (2011)

${\cal
P}_{\zeta}=\frac{H^{2}}{8\pi^{2}Q_{s}c_{s}^{3}}\biggr{|}_{c_{s}k=aH}\,,\qquad{\cal
P}_{h}=\frac{H^{2}}{2\pi^{2}Q_{t}c_{t}^{3}}\biggr{|}_{c_{t}k=aH}\,.$ (24)

In theories with $\mu(\phi)=0$, we have $c_{s}^{2}=1=c_{t}^{2}$ and hence both
${\cal P}_{\zeta}$ and ${\cal P}_{h}$ should be evaluated at $k=aH$.

Since $Q_{s}\simeq 3\xi\dot{\phi}^{2}/H^{2}$ and $Q_{t}\simeq\xi\phi^{2}/4$ in
the regimes $\xi\gg 1$ and $\xi\phi^{2}\gg M_{\rm Pl}^{2}$, the power spectra
(24) reduce to

$\displaystyle{\cal
P}_{\zeta}\simeq\frac{H^{4}}{24\pi^{2}\xi\dot{\phi}^{2}}\simeq\frac{\lambda\phi^{4}}{128\pi^{2}M_{\rm
Pl}^{4}}\simeq\frac{\lambda N^{2}}{72\pi^{2}\xi^{2}}\,,$ (25)
$\displaystyle{\cal
P}_{h}\simeq\frac{2H^{2}}{\pi^{2}\xi\phi^{2}}\simeq\frac{\lambda}{6\pi^{2}\xi^{2}}\,,$
(26)

where we used the approximate background Eq. (22). We note that the subscript
$k=aH$ is omitted in Eqs. (25) and (26). Then, we obtain the scalar spectral
index $n_{s}$ and the tensor-to-scalar ratio $r$, as

$\displaystyle n_{s}-1$ $\displaystyle=$ $\displaystyle\frac{{\rm d}\ln{\cal
P}_{\zeta}}{{\rm d}\ln
k}\biggr{|}_{k=aH}\simeq\frac{4\dot{\phi}}{H\phi}\simeq-\frac{8M_{\rm
Pl}^{2}}{3\xi\phi^{2}}\simeq-\frac{2}{N}\,,$ (27) $\displaystyle r$
$\displaystyle=$ $\displaystyle\frac{{\cal P}_{h}}{{\cal
P}_{\zeta}}\simeq\frac{12}{N^{2}}\,.$ (28)

Taking $N=60$ for scales relevant to the observed CMB temperature
anisotropies, we obtain $n_{s}=0.9667$ and $r=3.3\times 10^{-3}$. These values
are consistent with the bounds $n_{s}=0.9661\pm 0.0040$ (68 % CL) and
$r<0.066$ (95 % CL) constrained by the Planck 2018 data Akrami  _et al._
(2020). The Planck normalization ${\cal P}_{\zeta}=2.1\times 10^{-9}$ with
$N=60$ gives the constraint $\lambda/\xi^{2}=4.1\times 10^{-10}$. If
$\lambda=0.1$, then $\xi=1.6\times 10^{4}$.

The above results are valid for slow-roll inflation with $\mu(\phi)=0$. In the
presence of the scalar-GB coupling, the background dynamics and perturbation
spectra are subject to modifications. In subsequent sections, we will address
this issue along with the problem of generating the source for PBHs.

##  III Effective potentials with plateau and bump

Let us proceed to the case in which the scalar-GB coupling $\mu(\phi)R_{\rm
GB}^{2}$ is present. If $\mu(\phi)$ is a smooth function whose time variation
is small during inflation, the inflaton field $\phi$ can slowly evolve along
the potential. In such a case, the primordial power spectra of scalar and
tensor perturbations are given by Eq. (24). Since $Q_{s}$ is proportional to
$\dot{\phi}^{2}$, a smaller inflaton velocity generally leads to a larger
amplitude of ${\cal P}_{\zeta}$. In the context of slow-roll inflation,
however, this enhancement of ${\cal P}_{\zeta}$ is limited by a small time
variation of $\dot{\phi}$.

If the scalar-GB coupling generates a period in which the field velocity
$\dot{\phi}$ temporally approaches 0, it is possible to realize the large
enhancement of ${\cal P}_{\zeta}$ for scales smaller than those of the
observed CMB temperature anisotropies. One possible choice of $\mu(\phi)$ is a
dilatonic coupling of the form $\mu(\phi)=\mu_{0}e^{-\lambda\phi}$ Gross and
Sloan (1987); Metsaev and Tseytlin (1987); Gasperini  _et al._ (1997); Kawai
_et al._ (1998); Cartier  _et al._ (2001); Calcagni  _et al._ (2005); Guo  _et
al._ (2007). However, this type of continuously varying functions affects not
only the scalar perturbation on particular scales but also that on other
scales including CMB. Moreover, it can happen that the dominance of the
scalar-GB coupling over the potential and nonminimal couplings leads to the
violation of stability conditions (18).

Instead, we consider the step-like coupling given by Kawai and Kim (2021a);
Zhang (2022); Khan and Yogesh (2022)

$\mu(\phi)=\mu_{0}\tanh{[\mu_{1}(\phi-\phi_{c})]},$ (29)

where $\mu_{0}$, $\mu_{1}$, and $\phi_{c}$ are constants. Around the field
value $\phi=\phi_{c}$, this coupling rapidly changes from the asymptotic
constant $\mu(\phi)=-\mu_{0}$ (for $\phi\ll\phi_{c}$) to the other asymptotic
constant $\mu(\phi)=+\mu_{0}$ (for $\phi\gg\phi_{c}$). Since the GB curvature
invariant is a topological term, the scalar-GB coupling does not affect the
cosmological dynamics in the two asymptotic regimes with constant $\mu(\phi)$.
Provided that $\phi_{c}$ is in the range $\phi_{f}<\phi_{c}<\phi_{\rm CMB}$,
where $\phi_{\rm CMB}$ is the field value about 60 e-foldings before the end
of inflation (with the field value $\phi_{f}$), it should be possible to
enhance the scalar power spectrum for scales smaller than those of observed
CMB temperature anisotropies. In Sec. IV, we will study whether the sufficient
generation of seeds for PBHs is possible, while satisfying observational
constraints on $n_{s}$ and $r$ on CMB scales.

Provided that the field kinetic term is sufficiently small during inflation,
we can ignore the $\dot{\phi}$-dependent terms in Eq. (10). Then, the
$\phi$-derivative of the effective potential $V_{\rm eff}$ is approximately
given by

$V_{{\rm eff},\phi}\simeq\frac{(\xi\phi^{2}+M_{\rm
Pl}^{2})(V_{,\phi}-12\xi\phi
H^{2}-24\mu_{,\phi}H^{4})}{(6\xi+1)\xi\phi^{2}+M_{\rm
Pl}^{2}+48H^{2}\mu_{,\phi}(\xi\phi+2H^{2}\mu_{,\phi})}\,.$ (30)

From Eq. (6), the Hubble parameter is approximately given by

$H^{2}\simeq\frac{V(\phi)}{3(M_{\rm Pl}^{2}+\xi\phi^{2})}\,.$ (31)

Substituting Eq. (31) into Eq. (30) and integrating it with respect to $\phi$,
the effective potential $V_{\rm eff}$ can be numerically known as a function
of $\phi$.

In the following, we classify the effective potential into three classes: (1)
plateau-type, (2) bump-type, and (3) intermediate-type. The set 1, 2, 3 model
parameters shown in Table 1 are the typical examples of plateau-, bump-, and
intermediate-types, respectively. The nonminimal coupling constant is fixed to
be $\xi=5000$ in all cases.

Since PBHs are treated as the main component of DM in this paper, we consider
the PBH mass range $10^{-16}M_{\odot}\lesssim M\lesssim 10^{-11}M_{\odot}$.
This gives a constraint on the value of $\phi_{c}$. The two constants
$\mu_{0}$ and $\mu_{1}$ determine the types of $V_{\rm eff}(\phi)$ mentioned
above. Instead of the parameter $\mu_{0}$, we will use the combination

$\tilde{\mu}_{0}\equiv\frac{1}{4}\mu_{0}\mu_{1}\lambda\,,$ (32)

which appears later in Eq. (36). The values of $\tilde{\mu}_{0}$ are chosen to
be close or not far away from the right hand side (RHS) of Eq. (36), see the
last two columns in Table 1. The Higgs self-coupling $\lambda$ is determined
by the observed amplitude of primordial curvature perturbations on CMB scales.
As we see in Table 1, $\lambda$ is of order 0.01 for three sets of model
parameters.

Table 1: Three sets of model parameters. For dimensionfull parameters, the units are given inside the squared parenthesis. [.5]0 |  $\xi$ |  $\lambda$ |  $\phi_{c}\hskip 5.69046pt[M_{\rm Pl}]$ |  $\mu_{1}\hskip 5.69046pt[M_{\rm Pl}^{-1}]$ |  $\tilde{\mu}_{0}\hskip 5.69046pt[10^{8}M_{\rm Pl}^{-1}]$ |  RHS of Eq. (36)$\hskip 5.69046pt[10^{8}M_{\rm Pl}^{-1}]$
---|---|---|---|---|---|---
[.5]0 Set 1 |  5000 |  0.0244211 |  0.0380 |  1000 |  1.564709 |  1.556127
[.5]0 Set 2 |  5000 |  0.0110810 |  0.0760 |  5000 |  0.249909 |  0.176768
[.5]0 Set 3 |  5000 |  0.0152511 |  0.0600 |  1600 |  0.376125 |  0.366512

###  III.1 Plateau type

Thanks to the existence of the scalar-GB coupling, there is a stationary fixed
point $\phi=\phi_{*}$ at which $V_{{\rm eff},\phi}$ vanishes Kawai and Kim
(2021a). From Eq. (30) with Eq. (31), there is the following relation

$V_{,\phi}-\frac{4\xi\phi V}{M_{\rm
Pl}^{2}+\xi\phi^{2}}-\frac{8\mu_{,\phi}V^{2}}{3(M_{\rm
Pl}^{2}+\xi\phi^{2})^{2}}\biggr{|}_{\phi=\phi_{*}}=0\,.$ (33)

This is known as the USR regime in which the scalar-field equation (9) reduces
to Garcia-Bellido and Ruiz Morales (2017); Kannike  _et al._ (2017); Germani
and Prokopec (2017); Motohashi and Hu (2017)

$\ddot{\phi}+3H\dot{\phi}\simeq 0\qquad\text{(around $\phi=\phi_{*}$)}\,.$
(34)

The solution to this equation is given by

$\dot{\phi}\propto a^{-3}\propto e^{-3n}\qquad\text{(around
$\phi=\phi_{*}$)}\,,$ (35)

where $n=\ln a$ is the number of e-foldings counted forward. Hence
$\dot{\phi}$ rapidly decreases toward 0 in the USR regime.

Setting $\phi_{*}=\phi_{c}$ for the coupling (29), the moment at which
$V_{{\rm eff},\phi}$ vanishes coincides with the instant at transition of
$\mu(\phi)$. For the potential $V(\phi)=\lambda\phi^{4}/4$, the condition (33)
translates to

$\tilde{\mu}_{0}=\frac{3M_{\rm Pl}^{2}(M_{\rm
Pl}^{2}+\xi\phi_{c}^{2})}{2\phi_{c}^{5}}\,,$ (36)

which gives a constraint between $\mu_{0},\mu_{1},\phi_{c}$ and $\lambda,\xi$.
During the USR regime, the variation of the scalar field is of order
$1/\mu_{1}$. On using Eqs. (22) and (35), we can estimate the order of
$1/\mu_{1}$ as

$\frac{1}{\mu_{1}}\simeq\int\frac{|\dot{\phi}|}{H}{\rm d}n\simeq\int^{\Delta
N_{c}}_{0}\frac{2M_{\rm Pl}^{2}}{3\xi\phi_{c}}e^{-3n}{\rm d}n=\frac{2M_{\rm
Pl}^{2}}{9\xi\phi_{c}}\left(1-e^{-3\Delta N_{c}}\right)\,,$ (37)

where $\Delta N_{c}$ is the number of e-foldings acquired during the USR
phase.

Figure 1:  Shapes of the effective potential $V_{\rm eff}(\phi)$ around
$\phi=\phi_{c}$, where $V_{{\rm eff},c}=V_{\rm eff}(\phi_{c})$. Each plot
shows (1) plateau-type (left), (2) bump-type (middle), and (3) intermediate-
type (right).

Figure 2:  Evolution of $\phi$ during inflation (left) and $|\dot{\phi}|$
around $\phi=\phi_{c}$ (right) for the three sets of parameters given in Table
1. Here, $n$ is the number of e-foldings counted forward toward the end of
inflation, where $n_{\rm CMB}$ and $n_{c}$ are the e-foldings corresponding to
the field values $\phi_{\rm CMB}$ and $\phi_{c}$ respectively. In the right
panel, the dotted lines represent the evolution of $\dot{\phi}$ in the USR
regime, i.e., $\dot{\phi}\propto e^{-3n}$.

In the left panel of Fig. 1, we plot $V_{\rm eff}(\phi)$ in the vicinity of
$\phi=\phi_{c}$ for the Set 1 model parameters given in Table 1. In this case,
$\tilde{\mu}_{0}$ is chosen to be close to the RHS of Eq. (36). There is a
plateau of $V_{\rm eff}(\phi)$ in the region $|\phi-\phi_{c}|\lesssim 5\times
10^{-4}M_{\rm Pl}$, where $\phi_{c}=0.0380M_{\rm Pl}$. In the left panel of
Fig. 2, we show the evolution of $\phi$ versus the number of e-foldings for
the parameter Set 1 as a red line. The field value $\phi_{f}$ at the end of
inflation is numerically derived by the condition $\epsilon_{H}=1$, which
gives $\phi_{f}=0.0077M_{\rm Pl}$.

The initial field value for realizing the total number of e-foldings $n-n_{\rm
CMB}=60$ by the end of inflation corresponds to $\phi_{\rm CMB}\simeq
0.105M_{\rm Pl}$ (at which $n=n_{\rm CMB}$). In this case, the first slow-roll
stage of inflation is followed by the USR period starting at $n-n_{\rm
CMB}\simeq 36$. After the inflaton approaches the plateau of $V_{\rm
eff}(\phi)$ around $\phi=\phi_{c}$, the field derivative rapidly decreases as
$|\dot{\phi}|\propto e^{-3n}$, see the right panel of Fig. 2. For the model
parameters of Set 1, the number of e-foldings acquired during the USR epoch is
18.6. Finally, the scalar field exits from the USR regime, after which
$|\dot{\phi}|$ starts to increase.

As the plateau region of $V_{\rm eff}(\phi)$ gets wider, the number of
e-foldings $\Delta N_{c}$ acquired during the USR phase tends to be larger.
For $\Delta N_{c}$ exceeding the order 10, the CMB observables like $n_{s}$
and $r$ are subject to modifications in comparison to those derived for
$\mu(\phi)=0$. We will discuss this issue in Sec. IV.

###  III.2 Bump type

For the plateau-type effective potential, the scalar-GB coupling balances the
contributions arising from the potential and nonminimal couplings in the
scalar-field equation of motion. On the other hand, it should be possible that
the scalar-GB term temporarily becomes larger than the contributions from
other terms. This causes an instantaneous slowdown of the inflaton velocity in
a manner different from the USR case discussed in Sec. III.1. We call this
class as a bump-type model, in which the $\phi$ derivative of $V_{\rm
eff}(\phi)$ has the following feature

$V_{{\rm eff},\phi}(\phi)=\begin{cases}\hskip 5.69046pt>0\hskip 19.91684pt{\rm
for}~{}\phi\gg\phi_{c}\,,\\\ \hskip 5.69046pt<0\hskip 19.91684pt{\rm
around}~{}\phi=\phi_{c}\,,\\\ \hskip 5.69046pt>0\hskip 19.91684pt{\rm
for}~{}\phi\ll\phi_{c}.\end{cases}$ (38)

The parameter Set 2 in Table 1 gives rise to an effective potential of the
bump-type, which is illustrated in the middle panel of Fig. 1. In this case,
$\tilde{\mu}_{0}$ exhibits some deviation from the value on the RHS of Eq.
(36). The effective potential has a local maximum as well as a local minimum
in the vicinity of $\phi=\phi_{c}$. We note that similar toy models have been
studied in Refs. Atal  _et al._ (2019, 2020); Inomata  _et al._ (2022a); Cai
_et al._ (2022a, b) in different contexts.

If $\tilde{\mu}_{0}$ is larger than the RHS of Eq. (36), the scalar-GB
coupling dominates over the contributions from the potential and nonminimal
couplings. However, the period of the dominance must be sufficiently short to
end inflation properly. This requires that the parameter $\mu_{1}$ is quite
large. In the limit $\mu_{1}\rightarrow\infty$, $\mu_{,\phi}$ and
$\mu_{,\phi\phi}$ behave as

$\displaystyle\lim_{\mu_{1}\to\infty}\mu_{,\phi}=2\mu_{0}\delta(\phi-\phi_{c})\,,$
(39)
$\displaystyle\lim_{\mu_{1}\to\infty}\mu_{,\phi\phi}=\mu_{0}\mu_{1}\left[\delta(\phi-(\phi_{c}-\epsilon))-\delta(\phi-(\phi_{c}+\epsilon))\right]\,,$
(40)

where $\epsilon=\text{arctanh}(1/\sqrt{3})/\mu_{1}$. This type of step-like
behavior for large $\mu_{1}$ may induce Laplacian instabilities of
cosmological perturbations, so we will address this issue in Sec. IV by
computing the values of $c_{s}^{2}$ and $c_{t}^{2}$.

For the Set 2 model parameters corresponding to the bump-type effective
potential, we plot the evolution of $\phi$ and $|\dot{\phi}|$ as a green line
in Fig. 2. Unlike the plateau model, the field velocity decreases faster than
$|\dot{\phi}|\propto e^{-3n}$ due to the existence of the region $V_{{\rm
eff},\phi}<0$ around $\phi=\phi_{c}$. When the field reaches a local maximum
of $V_{\rm eff}(\phi)$, however, this period of the rapid decrease of
$\dot{\phi}$ soon comes to end. After this short epoch, the scalar field
quickly returns back to the slow-roll evolution. As we see in the left panel
of Fig. 2, the number of e-foldings $\Delta N_{c}$ acquired during the
transient phase around $\phi=\phi_{c}$ is only a few, which is much smaller
than $\Delta N_{c}$ in the USR case.

###  III.3 Intermediate type

Besides the two types of $V_{\rm eff}(\phi)$ discussed above, there is also an
intermediate case between the plateau- and bump-types. In this case the
effective potential is not exactly flat in the vicinity of $\phi=\phi_{c}$,
but it has a small peak and trough with a slight negative value of $V_{{\rm
eff},\phi}$ around $\phi=\phi_{c}$. The Set 3 parameters in Table 1 give rise
to such a shape of $V_{\rm eff}(\phi)$, see the right panel of Fig. 1.

As we plot as a blue line in Fig. 2, the field derivative initially decreases
in proportion to $|\dot{\phi}|\propto e^{-3n}$, which is followed by a
temporal period in which the decreasing rate of $|\dot{\phi}|$ becomes faster
than $|\dot{\phi}|\propto e^{-3n}$. In this latter regime, the scalar field
loses its velocity by climbing up the potential hill with $V_{{\rm
eff},\phi}<0$. After the field reaches the local maximum of $V_{\rm
eff}(\phi)$, $|\dot{\phi}|$ starts to grow toward the slow-roll regime. Since
the effective potential has neither an exactly flat region nor a sharp bump,
the number of e-foldings $\Delta N_{c}$ acquired around $\phi=\phi_{c}$ is
between those of plateau- and bump-types. In Set 3 model parameters, we have
$\Delta N_{c}\simeq 9.7$.

##  IV Generation of the seed for PBHs

In this section, we study how the power spectrum of curvature perturbations is
enhanced by the presence of the scalar-GB coupling (29). As we discussed in
Sec. III, the inflaton effective potential $V_{\rm eff}(\phi)$ can be
classified into three classes: (1) plateau-type, (2) bump-type, and (3)
intermediate-type. Examples of the model parameters corresponding to each type
of $V_{\rm eff}(\phi)$ are given in Table 1 as Sets 1, 2, 3, respectively. In
Sec. IV.1, we first discuss whether each model satisfies the stability
conditions (18). In Sec. IV.2, we compute the primordial scalar power spectra
by paying particular attention to the enhancement of curvature perturbations
around $\phi=\phi_{c}$. In Sec. IV.3, we confront our models with the observed
scalar spectral index and tensor-to-scalar ratio on CMB scales.

###  IV.1 Stability conditions

Let us first study whether the stability conditions (18) can be satisfied
during the whole stage of inflation. As we alluded in Sec. II.1, the no-ghost
parameter $Q_{t}$ in the tensor sector is positive for $\xi
H\phi\dot{\phi}<0$. For the three sets of model parameters in Table 1, the
field derivative $\dot{\phi}$ is negative without reaching 0, see the right
panel of Fig. 2. Then, we have $Q_{t}>0$ even during the transient epoch
around $\phi=\phi_{c}$. From Eq. (19), the other no-ghost parameter $Q_{s}$ is
also positive for $Q_{t}>0$. Thus, the ghost instabilities are absent for both
tensor and scalar perturbations, whose property is also confirmed numerically.

Figure 3:  Evolutions of $c_{t}^{2}$ (left) and $c_{s}^{2}$ (right) for the
three sets of model parameters given in Table 1.

The scalar-GB coupling generally gives rise to the tensor and scalar
propagation speeds different from 1. In the left panel of Fig. 3, we show the
evolution of $c_{t}^{2}$ versus the number of e-foldings $n-n_{\rm CMB}$ for
the three sets of model parameters in Table 1. Just after the entry to the
transient regime around $\phi=\phi_{c}$, $c_{t}^{2}$ exhibits the deviation
from 1. During the USR stage realized for the Set 1 model parameters, the
field velocity is significantly suppressed and hence $c_{t}^{2}$ approaches 1.
Just after the inflaton field exits from the USR regime, there is the temporal
deviation of $c_{t}^{2}$ from 1. In the subsequent slow-roll regime,
$c_{t}^{2}$ again goes back to 1. For the bump-type potential, which
corresponds to the Set 2 model parameters, the time variation of $c_{t}^{2}$
occurs instantly in the vicinity of $\phi=\phi_{c}$. In all the cases shown in
the left panel of Fig. 3, the deviation of $c_{t}^{2}$ from 1 is insignificant
($|c_{t}^{2}-1|\lesssim 0.04$) and hence the stability condition $c_{t}^{2}>0$
is always satisfied.

In the right panel of Fig. 3, we plot $c_{s}^{2}$ versus $n-n_{\rm CMB}$ for
the same sets of model parameters as those used in the left. For scalar
perturbations, the no-ghost quantity $Q_{s}$ is proportional to
$\dot{\phi}^{2}$. Since $Q_{s}$ appears in the denominator of $c_{s}^{2}$ in
Eq. (14), a smaller $\dot{\phi}^{2}$ does not imply a value of $c_{s}^{2}$
closer to 1. Indeed, during the transient regime around $\phi=\phi_{c}$, the
scalar propagation speeds deviate from 1 in all the three cases shown in Fig.
3. For Set 1, the minimum value of $c_{s}^{2}$ reached during the USR regime
is about 0.91. The bump-type potential corresponds to the Set 2 model
parameters, in which case $c_{s}^{2}$ temporally increases to the superluminal
region and then quickly evolves to the minimum around $c_{s}^{2}=0.90$. For
Set 3, $c_{s}^{2}$ decreases to the minimum around 0.35 and returns back to
the value close to 1 in the slow-roll regime. Since the positivity of
$c_{s}^{2}$ holds in all these cases, there are no Laplacian instabilities for
scalar perturbations.

Depending on the model parameters, there are cases in which the scalar sound
speed enters the region $c_{s}^{2}<0$. Since such models should be excluded by
the Laplacian instability, we will focus on the case $c_{s}^{2}>0$ in
subsequent sections.

###  IV.2 Primordial scalar power spectrum

To study the evolution of curvature perturbations during inflation, we
introduce the “sound-horizon” time defined by

$\tau_{s}=\int{\rm d}t\frac{c_{s}}{a}\,,$ (41)

where $c_{s}>0$. Then, the second-order action (13) of scalar perturbations
reduces to

${\cal S}_{s}^{(2)}=\int{\rm d}\tau_{s}{\rm
d}^{3}x\,a^{2}Q_{s}c_{s}\left[\zeta^{\prime 2}-(\nabla\zeta)^{2}\right]\,,$
(42)

where a prime represents the derivative with respect to $\tau_{s}$. We
decompose the curvature perturbation into the Fourier modes $\zeta_{k}$ as

$\zeta(\tau_{s},{\bm{x}})=\frac{1}{(2\pi)^{3}}\int{\rm
d}^{3}k\left[\zeta_{k}(\tau_{s},{\bm{k}})a({\bm{k}})+\zeta_{k}^{*}(\tau_{s},-{\bm{k}})a^{\dagger}(-{\bm{k}})\right]e^{i{\bm{k}}\cdot{\bm{x}}}\,,$
(43)

where ${\bm{k}}$ is a comoving wavenumber, and $a({\bm{k}})$ and
$a^{\dagger}({\bm{k}})$ are annihilation and creation operators, respectively.
Introducing the rescaled field

$u_{k}=Z_{s}\zeta_{k}\,,\qquad{\rm with}\qquad Z_{s}=a\sqrt{2Q_{s}c_{s}}\,,$
(44)

we obtain the following differential equation

$u_{k}^{\prime\prime}+\left(k^{2}-\frac{Z_{s}^{\prime\prime}}{Z_{s}}\right)u_{k}=0\,.$
(45)

For the perturbations deep inside the sound-horizon ($k^{2}\gg
a^{2}H^{2}/c_{s}^{2}$), $Z_{s}^{\prime\prime}/Z_{s}$ is suppressed relative to
$k^{2}$. For such modes, we choose a positive-frequency solution in a Bunch-
Davies vacuum state, i.e.,

$u_{k}=\frac{1}{\sqrt{2k}}e^{-ik\tau_{s}}\,,$ (46)

as an initial condition. In practice, we set the initial time for each
wavenumber that corresponds to 6 e-foldings before the sound-horizon crossing,
i.e., at $k=e^{6}aH/c_{s}$. Numerically we integrate Eq. (45) by the end of
inflation (characterized by the time $\tau_{s}=\tau_{sf}$) and compute the
scalar power spectrum given by

${\cal
P}_{\zeta}(k)=\frac{k^{3}}{2\pi^{2}}\left|\zeta_{k}(\tau_{sf},{\bm{k}})\right|^{2}\,.$
(47)

We evaluate ${\cal P}_{\zeta}(k)$ for three sets of model parameters given in
Table 1. At the 60 e-foldings before the end of inflation, the observed
amplitude ${\cal P}_{\zeta}=2.1\times 10^{-9}$ on CMB scales is used to give a
constraint among the model parameters. The two constants $\mu_{0}$ and
$\mu_{1}$ in Table 1 are chosen to realize the sufficient enhancement of
${\cal P}_{\zeta}$ around $\phi=\phi_{c}$.

In Fig. 4, we show the primordial scalar power spectrum ${\cal P}_{\zeta}(k)$
for three sets of model parameters. In comparison to the value of ${\cal
P}_{\zeta}$ on the CMB scale (characterized by the comoving wavenumber of
order $k_{\rm CMB}=10^{-3}~{}$Mpc), there is an enhancement of ${\cal
P}_{\zeta}$ by a factor of $10^{7}$. While the heights of peaks are similar
between the three sets of parameters, the widths of ${\cal P}_{\zeta}$ are
different from each other. The latter property is mostly attributed to the
fact that the number of e-foldings $\Delta N_{c}$ acquired around
$\phi=\phi_{c}$ are different between the plateau-, bump-, and intermediate-
types.

Figure 4:  The scalar power spectrum ${\cal P}_{\zeta}$ (normalized by its
value on the CMB scale ${\cal P}_{\zeta,{\rm CMB}}$) versus the comoving
wavenumber $k$ (normalized by the wavenumber $k_{\rm CMB}$ on the CMB scale).
Each plot corresponds to the three sets of model parameters presented in Table
1.

Let us try to understand how the enhancement of ${\cal P}_{\zeta}$ around
$\phi=\phi_{c}$ occurs in detail. For the perturbations outside the sound
horizon ($k^{2}\ll a^{2}H^{2}/c_{s}^{2}$), Eq. (45) is approximately given by

$u_{k}^{\prime\prime}-\frac{Z_{s}^{\prime\prime}}{Z_{s}}u_{k}\simeq 0\,.$ (48)

The solution to this equation can be generally expressed in the form

$\zeta_{k}(\tau_{s})=\frac{u_{k}(\tau_{s})}{Z_{s}(\tau_{s})}=A_{k}+B_{k}\int^{\tau_{s}}\frac{{\rm
d}\tilde{\tau}_{s}}{Z_{s}^{2}(\tilde{\tau}_{s})}\,,$ (49)

where $A_{k}$ and $B_{k}$ are constants. On using the approximations
$\Theta\simeq H\xi\phi^{2}$ and $Q_{t}\simeq V/(12H^{2})$ in the regimes
$\xi\phi^{2}\gg M_{\rm Pl}^{2}$ and $\xi\gg 1$, the no-ghost parameter $Q_{s}$
in Eq. (19) approximately reduces to

$Q_{s}\simeq\frac{4(3\xi^{2}+\lambda\mu_{,\phi}\phi)^{2}}{\lambda\xi^{2}\phi^{2}}\dot{\phi}^{2}\,,$
(50)

which is positive. Since we are considering the scalar-field evolution with
$\dot{\phi}<0$, the quantity $Z_{s}$ has the dependence

$Z_{s}\simeq-2\sqrt{\frac{2c_{s}}{\lambda}}\frac{3\xi^{2}+\lambda\mu_{,\phi}\phi}{\xi\phi}a\dot{\phi}\,.$
(51)

During the slow-roll regime in which the contribution from the scalar-GB
coupling is negligible, the field derivative $\dot{\phi}$ changes slowly with
$c_{s}\simeq 1$. In the limit that $\dot{\phi}$ is constant, we have
$Z_{s}\propto a$ and hence the last term of Eq. (49) decays in proportion to
$a^{-3}$ on the quasi de-Sitter background (where $a\simeq-(H\tau_{s})^{-1}$).
Then, after the sound horizon crossing, $\zeta_{k}$ approaches a constant
value $A_{k}$.

On the other hand, the field derivative exhibits a temporal rapid decrease
during the transition around $\phi=\phi_{c}$. Since the scalar-field solution
in the USR regime arising from the plateau-type effective potential is given
by $\dot{\phi}\propto a^{-3}$, neglecting the time dependence of $c_{s}$ leads
to the approximate relation $Z_{s}\propto a^{-2}$. Then, the last term in Eq.
(49) increases in proportion to $a^{3}$. This means that the decaying mode in
the slow-roll regime is replaced by the rapidly growing mode in the USR
regime. Since $\dot{\phi}$ also decreases rapidly for bump- and intermediate-
type potentials, the curvature perturbation can be strongly enhanced around
$\phi=\phi_{c}$ as well.

Figure 5:  Evolution of $Z_{s}$ (left) and
$(aH)^{-2}Z_{s}^{\prime\prime}/Z_{s}$ (right) around $n=n_{c}$ for three sets
of model parameters given in Table 1. Note that $Z_{s}$ is normalized by the
value $Z_{s,c}$ at $n=n_{c}$.

In the left panel of Fig. 5, we plot $Z_{s}$ versus $n-n_{c}$ for three sets
of model parameters in Table 1. During the initial slow-roll regime, the
evolution of $Z_{s}$ is approximately given by $Z_{s}\propto a$ in all these
three cases. For Set 1, which corresponds to the plateau-type potential, we
find that $Z_{s}$ decreases as $Z_{s}\propto a^{-2}$ in the vicinity of
$\phi=\phi_{c}$ as expected. As we see in the right panel of Fig. 3, the
deviation of $c_{s}^{2}$ from 1 is insignificant for Set 1 model parameters
and hence the evolution of $Z_{s}$ is hardly affected by the variation of
$c_{s}^{2}$.

The bump-type potential realized by Set 2 model parameters leads to a larger
decreasing rate of $Z_{s}$ around $\phi=\phi_{c}$ in comparison to the
plateau-type because of the stronger suppression of $\dot{\phi}$ (see the
right panel of Fig. 2). For Set 2, the scalar sound speed temporally enters
the region $c_{s}^{2}>1$ with a shorter transient period around
$\phi=\phi_{c}$ in comparison to Set 1. For Set 3, which corresponds to the
intermediate-type potential, the decreasing rate of $Z_{s}$ is between those
of the plateau- and bump-types. In all these cases, $Z_{s}$’s begin to
increase after reaching their minima, whose behavior is correlated with the
evolution of $\dot{\phi}$ plotted in Fig. 2.

In the superhorizon regime, the enhancement of $\zeta_{k}$ depends on the
integral of $1/Z_{s}^{2}(\tau_{s})$ with respect to $\tau_{s}=\int{\rm
d}t\,c_{s}/a$. On the other hand, we recall that we ignored the $k^{2}$ term
in Eq. (45) relative to $Z_{s}^{\prime\prime}/Z_{s}$ in the above argument. To
understand which modes of $\zeta_{k}$ are subject to the amplification, we
explicitly compute $Z_{s}^{\prime\prime}/Z_{s}$ as

$\frac{Z_{s}^{\prime\prime}}{Z_{s}}=\frac{(aH)^{2}}{c_{s}^{2}}\left[2-\epsilon_{H}+\frac{3\epsilon_{Q}}{2}+\frac{\epsilon_{c}}{2}-\frac{1}{4}(\epsilon_{c}+\epsilon_{Q})(2\epsilon_{H}-\epsilon_{Q}+\epsilon_{c})+\frac{1}{2}(\epsilon_{c}\eta_{c}+\epsilon_{Q}\eta_{Q})\right]\,,$
(52)

where

$\epsilon_{c}\equiv\frac{\dot{c_{s}}}{Hc_{s}}\,,\qquad\epsilon_{Q}\equiv\frac{\dot{Q_{s}}}{HQ_{s}}\,,\qquad\eta_{c}\equiv\frac{\dot{\epsilon_{c}}}{H\epsilon_{c}}\,,\qquad\eta_{Q}\equiv\frac{\dot{\epsilon_{Q}}}{H\epsilon_{Q}}\,.$
(53)

In the regime of slow-roll inflation, we have the approximate relation
$Z_{s}^{\prime\prime}/Z_{s}\simeq 2(aH)^{2}$ and hence $\zeta_{k}$ soon
approaches a constant after the Hubble radius crossing ($k\lesssim aH$).
During the transient regime around $\phi=\phi_{c}$, the quantities defined in
Eq. (53) can be larger than order 1. In the right panel of Fig. 5, we plot
${\cal Z}\equiv(aH)^{-2}Z_{s}^{\prime\prime}/Z_{s}$ versus $n-n_{c}$ for three
sets of model parameters. For Set 1, ${\cal Z}$ starts to evolve from the
value close to 2, temporally shows some decrease, and again increases to the
value around 4. This means that, for the wavenumbers in the range $k\lesssim
2aH$, there is the enhancement of $\zeta_{k}$ when the scalar field evolves
along the plateau region of $V_{\rm eff}(\phi)$. For the perturbations which
crossed the Hubble radius in the preceding slow-roll period, the last integral
in Eq. (49) has already decayed sufficiently around the time at which $\phi$
approaches $\phi_{c}$. Hence the enhanced modes of $\zeta_{k}$ are those
crossed the Hubble radius during the transient epoch around $\phi=\phi_{c}$.

For Set 2, we observe in Fig. 5 that ${\cal Z}$ quickly reaches the value
around 25 in the range $0\lesssim n-n_{c}\lesssim 2$. This means that
curvature perturbations up to the wavenumber $k\lesssim 5aH$, which include
subhorizon modes, can be amplified. In comparison to Set 1, the shorter
transient period around $\phi=\phi_{c}$ still limits the range of $k$ for
enhanced modes, see Fig. 4. Nevertheless, the height of peak of ${\cal
P}_{\zeta}(k)$ in Set 2 is similar to that in Set 1 thanks to the rapid
increase of ${\cal Z}$. For Set 3, the enhanced scalar power spectrum spans in
the ranges of $k$ between those of Sets 1 and 2. In this case, the rapid
decrease of $c_{s}^{2}$ seen in the right panel of Fig. 3 generates a sharp
peak in ${\cal Z}$. This gives rise to a shape of ${\cal P}_{\zeta}(k)$ whose
peak structure is different from those in other two cases. As we will see in
Sec. V, the PBH abundance generated by these three sets of primordial power
spectra can be sufficiently large to serve as almost all DM.

###  IV.3 CMB constraints

The existence of the transient regime with strongly suppressed values of
$\dot{\phi}$ modifies the spectra of scalar and tensor perturbations on scales
relevant to the observed CMB temperature anisotropies. Besides Eq. (45), we
also numerically solve the equation of tensor perturbations following from the
second-order action (17). We then compute the CMB observables like $n_{s}$,
$r$, and ${\cal P}_{\zeta}$ around the $N_{\rm CMB}=60$ e-foldings backward
from the end of inflation.

As we see in Fig. 2, the inflaton field stays nearly constant in the region
around $\phi=\phi_{c}$. In this transient regime, the acquired number of
e-foldings $\Delta N_{c}$ is different depending on the model parameters. For
larger $\Delta N_{c}$, the field value $\phi_{\rm CMB}$ around the CMB scale
tends to be smaller. Since $\phi$ stays nearly constant during the transient
regime, we can simply replace the relation $\phi^{2}\simeq 4M_{\rm
Pl}^{2}N/(3\xi)$ (which was derived for $\mu(\phi)=0$ in the limit $N\gg 1$)
with $\phi\to\phi_{\rm CMB}$ and $N\to N_{\rm CMB}-\Delta N_{c}$. Then, it
follows that

$\phi_{\rm CMB}\simeq 2\sqrt{\frac{N_{\rm CMB}-\Delta N_{c}}{3\xi}}M_{\rm
Pl}\,.$ (54)

In the USR case corresponding to the Set 1 model parameters in Fig. 2, we have
$\Delta N_{c}=18.6$ with $N_{\rm CMB}=60$ and hence $\phi_{\rm CMB}\simeq
0.105M_{\rm Pl}$ from Eq. (54). This shows good agreement with the numerical
value of $\phi_{\rm CMB}$.

From Eqs. (25) and (26), the scalar power spectrum on CMB scales is known by
the replacement $N\to N_{\rm CMB}-\Delta N_{c}$, while the tensor power
spectrum is hardly subject to modifications by the scalar-GB coupling.
Applying the change $N\to N_{\rm CMB}-\Delta N_{c}$ to Eqs. (25), (27), and
(28), we obtain

$\displaystyle{\cal P}_{\zeta}\bigr{|}_{\text{CMB}}$ $\displaystyle\simeq$
$\displaystyle\frac{\lambda(N_{\text{CMB}}-\Delta
N_{c})^{2}}{72\pi^{2}\xi^{2}}\,,$ (55) $\displaystyle n_{s}-1$
$\displaystyle\simeq$ $\displaystyle-\frac{2}{N_{\text{CMB}}-\Delta N_{c}}\,,$
(56) $\displaystyle r$ $\displaystyle\simeq$
$\displaystyle\frac{12}{(N_{\text{CMB}}-\Delta N_{c})^{2}}\,.$ (57)

For larger $\Delta N_{c}$, these observables exhibit more significant
deviations from those in standard Higgs inflation. The values of $\lambda$ in
Table 1 are chosen to match the amplitude ${\cal
P}_{\zeta}|_{\text{CMB}}=2.1\times 10^{-9}$ constrained by the Planck CMB
data.

Table 2: Numerical values of the scalar spectral index $n_{s}$ and the tensor-to-scalar ratio $r$ for three sets of model parameters given in Table 1 and Higgs inflation with $\mu(\phi)=0$. The number of e-foldings on the CMB scale is fixed to be $N_{\rm CMB}=60$.  [.5]0 |  $n_{s}$ |  $r$
---|---|---
[.5]0 Set 1 |  0.951415 |  0.00757526
[.5]0 Set 2 |  0.965142 |  0.00347235
[.5]0 Set 3 |  0.960055 |  0.00476115
[.5]0 Higgs inflation |  0.966527 |  0.00323724
Figure 6:  Dark grey and light grey areas represent the $1\sigma$ and
$2\sigma$ observational regions, respectively, constrained by the joint data
analysis of Planck 2018 + BK14 + BAO at $k=0.002$ Mpc-1. The red, green, and
blue points correspond to the theoretical predictions of $n_{s}$ and $r$ for
the model parameters of Sets 1, 2, 3, respectively, while the yellow star
represents those of Higgs inflation. The number of e-foldings on the CMB scale
is fixed to be $N_{\rm CMB}=60$.

For the Set 1 model parameters in Table 1, substituting $\Delta N_{c}=20$ and
$N_{\rm CMB}=60$ into the analytic estimations of Eqs. (56) and (57) gives
$n_{s}=0.950$ and $r=0.00750$. These are close to the numerically derived
values presented in the first column of Table 2. In Fig. 6, we show the
$1\sigma$ and $2\sigma$ observational bounds constrained from Planck 2018 data
combined with the data of B-mode polarizations available from the BICEP2/Keck
field (BK14) and baryon acoustic oscillations (BAO) Akrami  _et al._ (2020).
The theoretical values of $n_{s}$ and $r$ for Set 1 are outside the $2\sigma$
observational contour, so the plateau-type effective potential with $\Delta
N_{c}=18.6$ is disfavored from the data. Unless we choose an unusually large
value of $N_{\rm CMB}$ exceeding 68, the model with $\Delta N_{c}\simeq 20$
does not enter the inside of the $2\sigma$ contour.

For Set 2, the number of e-foldings acquired around the bump region of $V_{\rm
eff}(\phi)$ is as small as $\Delta N_{c}=2.3$, so the CMB observables are
similar to those in standard Higgs inflation, see the second and fourth
columns in Table 2. As we observe in Fig. 6, the model with Set 2 is inside
the $1\sigma$ observational contour. For Set 3, we numerically obtain the
values of $n_{s}$ and $r$ given in the third column of Table 2, with $\Delta
N_{c}=9.7$. In this case, the model is between $1\sigma$ and $2\sigma$
observational contours in the ($n_{s},r$) plane. For $N_{\rm CMB}=60$, the
number of e-foldings acquired around $\phi=\phi_{c}$ should be in the ranges

$\displaystyle 0\leq\Delta N_{c}<7.6\qquad\quad~{}(1\sigma)\,,$ (58)
$\displaystyle 7.6\leq\Delta N_{c}<12.5\qquad(2\sigma)\,,$ (59)

for the consistency with observations of $n_{s}$ and $r$ at $1\sigma$ and
$2\sigma$ confidence levels, respectively. Thus, the bump-type model with
$\Delta N_{c}=$ a few is favored over the plateau-type model with $\Delta
N_{c}>12.5$ from the viewpoint of CMB constraints.

Figure 7: $1\sigma$ (yellow) and $2\sigma$ (blue) confidence level parameter
spaces in the $(\lambda,\xi)$ plane constrained from the observed CMB
temperature anisotropies with ${\cal{P}}_{\zeta}=2.1\times 10^{-9}$ and
$N_{\rm CMB}=60$. The $1\sigma$ and $2\sigma$ regions correspond to $\Delta
N_{c}$ in the ranges (58) and (59), respectively. The red region is outside
the $2\sigma$ observational contour.

On using the Planck normalized value ${\cal P}_{\zeta|{\rm CMB}}=2.1\times
10^{-9}$ in Eq. (55) together with Eq. (56), we obtain the following relations

$\frac{\lambda}{\xi^{2}}\simeq\frac{1.49\times 10^{-6}}{(N_{\rm CMB}-\Delta
N_{c})^{2}}\simeq 3.73\times 10^{-7}(n_{s}-1)^{2}\,.$ (60)

Fixing $N_{\rm CMB}$ to be 60, the scalar spectral index is in the range
$n_{s}=1-2/(60-\Delta N_{c})\leq 0.9667$ for $\Delta N_{c}\geq 0$. This means
that, from Eq. (60), the ratio $\lambda/\xi^{2}$ is in the range
$\lambda/\xi^{2}\geq 4.1\times 10^{-10}$. As $\Delta N_{c}$ increases,
$\lambda/\xi^{2}$ gets larger. The $1\sigma$ and $2\sigma$ confidence regions
of $\Delta N_{c}$, which are given by Eqs. (58) and (59) respectively,
translate to

$\displaystyle 4.1\times 10^{-10}\leq\frac{\lambda}{\xi^{2}}<5.4\times
10^{-10}\qquad(1\sigma)\,,$ (61) $\displaystyle 5.4\times
10^{-10}\leq\frac{\lambda}{\xi^{2}}<6.6\times 10^{-10}\qquad(2\sigma)\,,$ (62)

In Fig. 7, these parameter spaces are plotted as yellow ($1\sigma$) and blue
($2\sigma$) regions in the $(\lambda,\xi)$ plane. The red region is outside
the $2\sigma$ observational contour. The allowed parameter spaces shown in
Fig. 7 will be useful to put further constraints on the values of $\lambda$
and $\xi$ from future collider experiments.

##  V PBH abundance

After inflation the perturbations reenter the Hubble radius, whose epoch
depends on the comoving wavenumber $k$. In Sec. IV, we showed that the
presence of the scalar-GB coupling can lead to the sufficient amplification of
curvature perturbations during inflation for particular wavelengths smaller
than the CMB scale ($k^{-1}\simeq 10^{3}$ Mpc). Such overdense regions can
collapse to form PBHs after the horizon reentry.

The horizon mass associated with the Hubble distance $H^{-1}$ is given by
$M_{H}=4\pi M_{\rm Pl}^{2}H^{-1}$. The mass of PBHs at its formation time
$t_{\text{form}}$ can be expressed as Green  _et al._ (2004)

$M=\gamma M_{H}=4\pi\gamma M_{\rm Pl}^{2}H(t_{\rm form})^{-1}\,,$ (63)

where $\gamma$ is the ratio of how much of the inner region of the Hubble
radius collapses into PBHs. We will consider the case in which the formation
of PBHs occurs during the radiation-dominated epoch. The Hubble parameter at
$t=t_{\rm form}$ is related to today’s Hubble constant $H_{0}$ as Kawai and
Kim (2021a)

$\frac{H(t_{\rm form})}{H_{0}}=\sqrt{\Omega_{r0}}\left[\frac{a_{0}}{a(t_{\rm
form})}\right]^{2}\left[\frac{g_{*0}}{g_{*}(t_{\text{form}})}\right]^{1/6}\,,$
(64)

where the subscript “0” represents today’s values and $g_{*}$ is the
relativistic degrees of freedom. The wavenumber at horizon reentry corresponds
to $k=a(t_{\text{form}})H(t_{\text{form}})$, whose relation can be used to
eliminate $a(t_{\text{form}})$ in Eq. (64). Solving Eq. (64) for $H(t_{\rm
form})$ and substituting it into Eq. (63), it follows that

$M(k)=10^{-13}M_{\odot}\left(\frac{\gamma}{0.2}\right)\left[\frac{g_{*}(t_{\text{form}})}{106.75}\right]^{-1/6}\left(\frac{k}{4.9\times
10^{12}~{}\text{Mpc}^{-1}}\right)^{-2}\,,$ (65)

where we used the values $g_{*0}=3.36$, $\Omega_{r0}=9\times 10^{-5}$, and
$H_{0}=10^{-42}$ GeV.

The abundance of PBHs can be estimated by using the Press-Schechter theory.
Assuming a Gaussian distribution for the coarse-grained density fluctuation
$\delta_{g}$, the probability that $\delta_{g}$ is higher than a certain
threshold value $\delta_{c}$ at $t=t_{\text{form}}$ is given by Green  _et
al._ (2004); Young  _et al._ (2014); Harada  _et al._ (2013); Germani and
Musco (2019)

$\beta(M(k))=\int^{\infty}_{\delta_{c}}\frac{1}{\sqrt{2\pi\sigma^{2}(k)}}\exp\left[-\frac{\delta_{g}^{2}}{2\sigma^{2}(k)}\right]{\rm
d}\delta_{g}\,,$ (66)

where

$\delta_{g}(\bm{x},R)\equiv\int W(|\bm{x}-\bm{y}|,R)\delta(\bm{y}){\rm
d}^{3}y\,.$ (67)

The window function $W(|\bm{x}-\bm{y}|,R)$ determines how the density contrast
$\delta=(\rho-\bar{\rho})/\bar{\rho}$ around $\bm{x}$ is coarse-grained, where
$\rho$ is the density and $\bar{\rho}$ is its background part. We choose the
Gaussian window function of the form

$W(|\bm{x}-\bm{y}|,R)=\frac{1}{(2\pi)^{3/2}R^{3}}\exp\left(-\frac{|\bm{x}-\bm{y}|^{2}}{2R^{2}}\right)\,,$
(68)

where the distance $R$ is taken to be $k^{-1}$. Then, the variance of
$\delta_{g}$ is given by

$\sigma^{2}(k)=\frac{16}{81}\int\exp\left[-\left(\frac{p}{k}\right)^{2}\right]\left(\frac{p}{k}\right)^{4}{\cal{P}}_{\zeta}(p)\,{\rm
d}\ln p\,.$ (69)

Since the PBH density decreases as $\rho_{\rm PBH}\propto a^{-3}$ after its
formation, today’s PBH density can be estimated as

$\rho_{{\rm PBH},0}=\rho_{\rm PBH}(t_{\rm form})\left[\frac{a_{0}}{a(t_{\rm
form})}\right]^{-3}=\gamma\beta\rho_{r}(t_{\rm
form})\left[\frac{a_{0}}{a(t_{\rm form})}\right]^{-3}\,.$ (70)

Since $\rho_{r}(t_{\rm form})$ is related to $H(t_{\rm form})$ as
$\rho_{r}(t_{\rm form})=3M_{\rm Pl}^{2}H(t_{\rm form})^{2}$, today’s density
parameter of PBHs corresponding to Eq. (70) is

$\Omega_{{\rm PBH},0}=\frac{\rho_{{\rm PBH},0}}{3M_{\rm
Pl}^{2}H_{0}^{2}}=\gamma\beta\frac{H(t_{\rm
form})^{2}}{H_{0}^{2}}\left[\frac{a(t_{\rm form})}{a_{0}}\right]^{3}\,,$ (71)

with total PBH density parameter $\int\Omega_{{\rm PBH},0}\,{\rm d}\ln M$.
Then, we obtain the ratio of the PBH abundance in a mass range $[M,M+{\rm
d}\ln M]$ to the entire density of cold DM (today’s density parameter
$\Omega_{{\rm CDM},0}$) as

$f(M)=\frac{\Omega_{{\rm PBH},0}}{\Omega_{{\rm
CDM},0}}\simeq\left[\frac{\beta(M)}{1.04\times
10^{-14}}\right]\left(\frac{\gamma}{0.2}\right)^{3/2}\left[\frac{g_{*}(t_{\text{form}})}{106.75}\right]^{-1/4}\left(\frac{\Omega_{\text{CDM},0}\,h^{2}}{0.12}\right)^{-1}\left(\frac{M}{10^{-13}M_{\odot}}\right)^{-1/2}\,,$
(72)

where $H_{0}=100\,h$ km sec-1 Mpc-1.

Figure 8:  The ratio of PBH abundances $f(M)$ relative to all the cold DM
density as a function of the PBH mass $M$ (in the unit of solar mass
$M_{\odot}$) for the three sets of model parameters presented in Table 1. The
colored areas correspond to excluded regions from the evaporation of black
holes (dark yellow) and observational constraints from the microlensing (dark
blue). These observational constraints are taken from Refs. Kavanagh (2019);
Green and Kavanagh (2021). The list of references used to derive these
constraints can be found in
https://github.com/bradkav/PBHbounds/blob/master/bounds/README.md.

In the following, we use the values $\gamma=0.2$ Carr (1975), $g_{*}(t_{\rm
form})=106.75$, $\Omega_{\text{CDM},0}h^{2}=0.12$, and $\delta_{c}=0.4$ for
the computations of $M$ and $f(M)$. With a given wavenumber $k$ at horizon
reentry, the PBH mass $M(k)$ is known from Eq. (65). The variance (69) is
affected by the primordial spectrum ${\cal P}_{\zeta}(k)$ enhanced at some
particular scales during inflation. This modifies the PBH mass function $f(M)$
through the change of $\beta(M(k))$ in Eq. (66).

In Fig. 8, we plot $f(M)$ versus $M$ for the three sets of model parameters
presented in Table 1. In these three cases, the mass functions $f(M)$ span in
the range $10^{-16}M_{\odot}\lesssim M\lesssim 10^{-13}M_{\odot}$. As we
observe in Fig. 4, the wavenumber $k$ corresponding to the peak positions of
${\cal P}_{\zeta}(k)$ is smallest for Set 1, while largest for Set 3. From Eq.
(65), the PBH mass $M$ decreases for larger $k$. Then, the mass $M$
corresponding to the peaks of $f(M)$ is smallest for Set 3, while largest for
Set 1. The maximum values of $f(M)$ are found to be $0.85$, $0.71$, and $1$
for Sets 1, 2, 3, respectively, so PBHs are the source for practically all
cold DM in these three cases.

In Fig. 8, we also show the regions excluded by the black hole evaporation and
by microlensing observations. The models with Sets 1, 2, 3 are also consistent
with such constraints. While we have considered the PBH mass range
$10^{-16}M_{\odot}\lesssim M\lesssim 10^{-13}M_{\odot}$, it is also possible
to produce PBHs with the mass $M>10^{-13}M_{\odot}$ by choosing different sets
of model parameters. The heights of peaks of $f(M)$ can be also lower than the
order 0.1 to be consistent with the microlensing data. Thus, our model allows
versatile possibilities for generating PBHs in broad mass ranges.

Finally, we explicitly derive a relation between $M$ and the critical field
value $\phi_{c}$. This is useful for estimating the mass of PBH and
determining the model parameters. After the transient regime the inflationary
dynamics rapidly approaches the slow-roll solution, so we can exploit Eq. (23)
as a good approximation. Then, we obtain

$\frac{3\xi\phi_{c}^{2}}{4M_{\rm Pl}^{2}}\simeq N_{\rm CMB}-\ln\frac{a(t_{\rm
PBH})}{a(t_{\rm CMB})}-\Delta N_{c}\,,$ (73)

where $t_{\rm CMB}$ and $t_{\rm PBH}$ represent the moments at which the CMB
and PBH scales leave the Hubble horizon, respectively. From Eq. (65), we can
estimate $a(t_{\rm PBH})/a(t_{\rm CMB})$ as

$\frac{a(t_{\rm PBH})}{a(t_{\rm CMB})}\simeq\frac{k}{k_{\rm
CMB}}\simeq\left[\frac{10^{-13}M_{\odot}}{M(k)}\right]^{1/2}\left(\frac{\gamma}{0.2}\right)^{1/2}\left[\frac{g_{*}(t_{\text{form}})}{106.75}\right]^{-1/12}\left(\frac{4.9\times
10^{12}~{}\text{Mpc}^{-1}}{k_{\rm CMB}}\right)\simeq\left[\frac{6\times
10^{17}M_{\odot}}{M(k)}\right]^{1/2}\,,$ (74)

where we used the values $\gamma=0.2$, $g_{*}(t_{\text{form}})=106.75$, and
$k_{\rm CMB}=0.002~{}\text{Mpc}^{-1}$ in the last equality. Substituting Eq.
(74) into Eq. (73), it follows that

$\phi_{c}\simeq M_{\rm Pl}\sqrt{\frac{4}{3\xi}\left\\{N_{\rm CMB}-\Delta
N_{c}+\frac{1}{2}\ln\left[\frac{M(k)}{6\times
10^{17}M_{\odot}}\right]\right\\}}\,.$ (75)

In the USR case corresponding to Set 1 model parameters, for example, we have
$\Delta N_{c}=18.6$ and $M\simeq 4\times 10^{-14}M_{\odot}$, and hence
$\phi_{c}\simeq 0.0383M_{\rm Pl}$ from Eq. (75). This is in good agreement
with the exact value of $\phi_{c}$, i.e., $\phi_{c}=0.0380M_{\rm Pl}$. For the
PBH mass range $10^{-16}M_{\odot}\lesssim M\lesssim 10^{-11}M_{\odot}$ in
which PBHs can be the source for all DM, the corresponding region of
$\phi_{c}$ is given by

$M_{\rm Pl}\sqrt{\frac{4}{3\xi}(21-\Delta N_{c})}\lesssim\phi_{c}\lesssim
M_{\rm Pl}\sqrt{\frac{4}{3\xi}(27-\Delta N_{c})}\,.$ (76)

We recall that, for the plateau-type potential, $\tilde{\mu}_{0}$ and
$\mu_{1}$ are related to $\phi_{c}$ according to Eqs. (36) and (37),
respectively. By using these relations with Eq. (75), the orders of parameters
$\tilde{\mu}_{0}$ and $\mu_{1}$ are known for given values of $M(k)$ and
$\Delta N_{c}$. We caution that the relations (36) and (37) lose their
accuracy for the bump-type potential, but they are still useful to estimate
the orders of $\tilde{\mu}_{0}$ and $\mu_{1}$. To know the precise values of
$\tilde{\mu}_{0}$ and $\mu_{1}$ with which curvature perturbations are
sufficiently enhanced for scales relevant to PBHs, the numerical analysis is
required as we performed in Sec. IV. Basically, $\tilde{\mu}_{0}$ and
$\mu_{1}$ determine the height of enhancement of $\zeta$ and the length of
transient period, respectively.

##  VI Conclusions

We studied a mechanism for producing the seed of PBHs in Higgs inflation in
the presence of a scalar-GB coupling. This provides a minimal scenario of
inflation within a framework of standard model of particle physics, while
allowing for the compatibility with observed CMB temperature anisotropies. In
comparison to original Higgs inflation, however, the scalar-GB coupling can
modify theoretical predictions of the scalar spectral index and tensor-to-
scalar ratio on CMB scales. We explored the possibility for enhancing
curvature perturbations at some particular scales to generate the seed for
PBHs, while the model is still compatible with the CMB observations.

The enhancement of curvature perturbations during inflation is possible when
the inflaton velocity $\dot{\phi}$ rapidly decreases toward 0 during some
transient epoch. If there is a period in which the scalar-GB coupling
$\mu(\phi)$ quickly changes, it is possible to generate particular shapes in
the inflaton effective potential $V_{\rm eff}(\phi)$. For this purpose, we
considered a coupling function of the form
$\mu(\phi)=\mu_{0}\tanh[\mu_{1}(\phi-\phi_{c})]$, where $\phi_{c}$ is the
field value at transition. Since $\mu(\phi)$ approaches constants in two
asymptotic regimes $\phi\ll\phi_{c}$ and $\phi\gg\phi_{c}$, the scalar-GB
coupling affects the dynamics of background and perturbations only in the
vicinity of $\phi=\phi_{c}$. Depending on the model parameters, the effective
potential can be classified into three classes: (1) plateau-type, (2) bump-
type, (3) intermediate-type.

A typical example of the plateau-type is the Set 1 model parameters given in
Table 1, in which case the scalar-GB coupling almost balances the potential
and nonminimal coupling terms around $\phi=\phi_{c}$. For this type the
effective potential has a nearly flat region, as seen in the left panel of
Fig. 1. During the USR regime, the field derivative decreases as
$\dot{\phi}\propto a^{-3}$. As we observe in the left panel of Fig. 2 for Set
1, the number of e-foldings $\Delta N_{c}$ acquired around $\phi=\phi_{c}$ is
as large as 20. The Set 2 model parameters in Table 1 give rise to a bump-type
potential illustrated in the middle panel of Fig. 1. In this case, the
existence of an explicit peak in $V_{\rm eff}(\phi)$ leads to the approach of
$\dot{\phi}$ toward 0 even faster than the USR case. For Set 2, the number of
e-foldings $\Delta N_{c}$ acquired around $\phi=\phi_{c}$ is as small as a
few. The Set 3 model parameters in Table 1 correspond to the intermediate-type
potential plotted in the right panel of Fig. 1. For this type, both the
decreasing rate of $\dot{\phi}$ and $\Delta N_{c}$ are between those of
plateau- and bump-types.

The evolution of scalar perturbations during inflation crucially depends on
the quantity $Z_{s}$ defined in Eq. (44). After the sound horizon crossing,
the solution to the Fourier-transformed curvature perturbation $\zeta_{k}$ is
expressed in the form (49). In all the three types of $V_{\rm eff}(\phi)$
mentioned above, $Z_{s}$ rapidly decreases in the vicinity of $\phi=\phi_{c}$.
Since the last integral in Eq. (49) becomes a rapidly growing mode during the
transient regime, there is the strong enhancement of $\zeta_{k}$ for the modes
exiting the sound horizon around $\phi=\phi_{c}$. For the effective potential
closer to the bump-type, the quantity $Z_{s}^{\prime\prime}/Z_{s}$, which
appears in the equation of motion for $u_{k}=Z_{s}\zeta_{k}$, can exceed the
order of $10(aH)^{2}$ during the rapid transition. Then, the bump-type
potential also leads to the amplification of some subhorizon modes. Still, the
shapes of primordial scalar power spectra ${\cal P}_{\zeta}(k)$ mostly depend
on the number of e-foldings $\Delta N_{c}$, in such a way that the peak tends
to be sharper for smaller $\Delta N_{c}$.

The primordial power spectra ${\cal P}_{\zeta}(k)$ plotted in Fig. 4
correspond to those of Set 1, 2, 3 model parameters. In all these cases, the
peak values of ${\cal P}_{\zeta}(k)$ are about $10^{7}$ times as large as the
amplitude of ${\cal P}_{\zeta}(k)$ on CMB scales. With these model parameters,
all the stability conditions given in Eq. (18) are consistently satisfied and
hence there are neither ghost nor Laplacian instabilities. We note that the
inflaton-GB coupling gives rise to the scalar propagation speed different from
1, whose property also affects the shapes of ${\cal P}_{\zeta}(k)$.

The existence of the transient epoch around $\phi=\phi_{c}$ modifies the
scalar power spectrum and tensor-to-scalar ratio on CMB scales. In standard
Higgs inflation, the expressions of ${\cal P}_{\zeta}$, $n_{s}$, and $r$ are
given, respectively, by Eqs. (25), (27), and (28), where $N$ is the number of
e-foldings counted backward from the end of inflation. In the presence of the
scalar-GB coupling, these observables are subject to the modification $N\to
N_{\rm CMB}-\Delta N_{c}$, where $N_{\rm CMB}$ corresponds to the number of
e-foldings on CMB scales. This means that, for smaller $\Delta N_{c}$, the
models can exhibit better compatibility with the observational bounds on
$n_{s}$ and $r$. As we see in Fig. 6, the bump-type with $\Delta N_{c}=2.3$
(Set 2) is well inside the $1\sigma$ observational contour constrained from
CMB and other data, while the plateau-type with $\Delta N_{c}=18.6$ (Set 1) is
outside the $2\sigma$ contour. The intermediate-type model with $\Delta
N_{c}=9.7$ (Set 3) is between $1\sigma$ and $2\sigma$ contours. These results
show that the bump-type is generally favored over the plateau-type from the
CMB constraints.

In Fig. 8, we plot the PBH fraction function $f(M)$ versus the mass $M$ for
three sets of model parameters in Table 1. In all these cases we have
$f(M)\simeq 1$ in the mass range $10^{-16}M_{\odot}\lesssim M\lesssim
10^{-13}M_{\odot}$. Thus, the GB corrected Higgs inflation allows the
possibility for generating the amount of PBHs serving as almost all DM, while
being compatible with the observed temperature anisotropies on CMB scales
especially for the bump-type effective potential. We note that, in our model,
it is also possible to generate PBHs in the larger mass range
$M>10^{-13}M_{\odot}$ with the fraction $f(M)\lesssim 0.1$.

The enhanced primordial curvature perturbations on particular scales may
induce gravitational waves at nonlinear order, which affects the spectrum of
gravitational wave background Kohri and Terada (2018) (see also Ref.
Papanikolaou  _et al._ (2021)). In addition, while a Gaussian distribution was
used for the coarse-grained density fluctuation in this paper, it was reported
that the distribution function of curvature fluctuations may have a peculiar
shape deviating from the Gaussian Namjoo  _et al._ (2013); Franciolini  _et
al._ (2018); Cai  _et al._ (2018b); Atal and Germani (2019); Atal  _et al._
(2019, 2020); Passaglia  _et al._ (2019); Taoso and Urbano (2021); Biagetti
_et al._ (2021); Davies  _et al._ (2022); Cai  _et al._ (2022a, b). Recently,
it was shown that the PBH scenario in single-field inflation with an USR
regime can induce large one-loop corrections to the scalar power spectrum on
CMB scales Kristiano and Yokoyama (2022); Inomata  _et al._ (2022b). Since the
analysis is limited to a canonical scalar field with a plateau-type potential,
the results of Ref. Kristiano and Yokoyama (2022) are not applied to our model
in which the enhancement of curvature perturbations occurs by the presence of
the inflaton-GB coupling. However, it is worth computing one-loop corrections
to ${\cal P}_{\zeta}$ in our model, especially for the bump-type effective
potential. We leave these interesting issues for future works.

## Acknowledgments

We thank Tomohiro Fujita for useful discussions and comments. We also thank
Bradley Kavanagh for giving us permission to use observational constraints
given in Fig. 8. ST is supported by the Grant-in-Aid for Scientific Research
Fund of the JSPS Nos. 19K03854 and 22K03642.

##  Appendix A Appendix: Correspondence with Horndeski theory

The inflationary model studied in this paper belongs to a subclass of
Horndeski theory given by the action Horndeski (1974); Deffayet  _et al._
(2011); Kobayashi  _et al._ (2011); Charmousis  _et al._ (2012)

${\cal S}=\int{\rm d}^{4}x\sqrt{-g}\,{\cal L}_{H}\,,$ (77)

where

$\displaystyle{\cal L}_{H}$ $\displaystyle=$ $\displaystyle
G_{2}(\phi,X)-G_{3}(\phi,X)\square\phi+G_{4}(\phi,X)\,R+G_{4,X}(\phi,X)\left[(\square\phi)^{2}-(\nabla_{\mu}\nabla_{\nu}\phi)(\nabla^{\mu}\nabla^{\nu}\phi)\right]+G_{5}(\phi,X)G_{\mu\nu}\nabla^{\mu}\nabla^{\nu}\phi$
(78)
$\displaystyle-\frac{1}{6}G_{5,X}(\phi,X)\left[(\square\phi)^{3}-3(\square\phi)\,(\nabla_{\mu}\nabla_{\nu}\phi)(\nabla^{\mu}\nabla^{\nu}\phi)+2(\nabla^{\mu}\nabla_{\alpha}\phi)(\nabla^{\alpha}\nabla_{\beta}\phi)(\nabla^{\beta}\nabla_{\mu}\phi)\right]\,,$

with $G_{\mu\nu}$ being the Einstein tensor. Four functions $G_{j}$’s
($j=2,3,4,5$) depend on the scalar field $\phi$ and its kinetic term
$X=-g^{\mu\nu}\nabla_{\mu}\phi\nabla_{\nu}\phi/2$, where we use the notation
$G_{j,X}\equiv{\rm d}G_{j}/{\rm d}X$. The action (1) can be accommodated by
the following coupling functions Kobayashi  _et al._ (2011)

$\displaystyle
G_{2}(\phi,X)=X-\frac{1}{4}\lambda\phi^{4}+8\mu^{(4)}(\phi)X^{2}(3-\ln|X|)\,,$
(79) $\displaystyle G_{3}(\phi,X)=4\mu^{(3)}(\phi)X(7-3\ln|X|)\,,$ (80)
$\displaystyle G_{4}(\phi,X)=\frac{M_{\rm
Pl}^{2}}{2}+\frac{1}{2}\xi\phi^{2}+4\mu^{(2)}(\phi)X(2-\ln|X|)\,,$ (81)
$\displaystyle G_{5}(\phi,X)=-4\mu^{(1)}(\phi)\ln|X|\,,$ (82)

where $\mu^{(i)}(\phi)\equiv{{\rm d}^{i}\mu/{\rm d}\phi^{i}}$. In full
Horndeski theories, the background and perturbation equations of motion on the
flat FLRW background were already derived in Refs. Kobayashi  _et al._ (2011);
De Felice and Tsujikawa (2011b); De Felice  _et al._ (2011). In this paper, we
applied those results to the theory given by the coupling functions (79)-(82).

## References

  * Zel’dovich and Novikov (1967) Y. B. Zel’dovich and I. D. Novikov,  _Soviet Astronomy_ 10, 602 (1967).
  * Hawking (1971) S. Hawking,  _Mon. Not. Roy. Astron. Soc._ 152, 75 (1971).
  * Carr and Hawking (1974) B. J. Carr and S. W. Hawking,  _Mon. Not. Roy. Astron. Soc._ 168, 399 (1974).
  * Chapline (1975) G. F. Chapline,  _Nature_ 253, 251 (1975).
  * Meszaros (1975) P. Meszaros,  _Astron. Astrophys._ 38, 5 (1975).
  * Khlopov (2010) M. Y. Khlopov,  _Res. Astron. Astrophys._ 10, 495 (2010), arXiv:0801.0116 [astro-ph] .
  * Sasaki  _et al._ (2018) M. Sasaki, T. Suyama, T. Tanaka, and S. Yokoyama,  _Class. Quant. Grav._ 35, 063001 (2018), arXiv:1801.05235 [astro-ph.CO] .
  * Carr and Kuhnel (2020) B. Carr and F. Kuhnel,  _Ann. Rev. Nucl. Part. Sci._ 70, 355 (2020), arXiv:2006.02838 [astro-ph.CO] .
  * Green and Kavanagh (2021) A. M. Green and B. J. Kavanagh,  _J. Phys. G_ 48, 043001 (2021), arXiv:2007.10722 [astro-ph.CO] .
  * Villanueva-Domingo  _et al._ (2021) P. Villanueva-Domingo, O. Mena, and S. Palomares-Ruiz,  _Front. Astron. Space Sci._ 8, 87 (2021), arXiv:2103.12087 [astro-ph.CO] .
  * Carr and Kuhnel (2022) B. Carr and F. Kuhnel,  _SciPost Phys. Lect. Notes_ 48, 1 (2022), arXiv:2110.02821 [astro-ph.CO] .
  * Escrivà  _et al._ (2022) A. Escrivà, F. Kuhnel, and Y. Tada, arXiv:2211.05767 [astro-ph.CO] .
  * Karam  _et al._ (2022) A. Karam, N. Koivunen, E. Tomberg, V. Vaskonen, and H. Veermäe, arXiv:2205.13540 [astro-ph.CO] .
  * Abbott  _et al._ (2016) B. P. Abbott _et al._ (LIGO Scientific, Virgo),  _Phys. Rev. Lett._ 116, 061102 (2016), arXiv:1602.03837 [gr-qc] .
  * Abbott  _et al._ (2019) B. P. Abbott _et al._ (LIGO Scientific, Virgo),  _Phys. Rev. X_ 9, 031040 (2019), arXiv:1811.12907 [astro-ph.HE] .
  * Abbott  _et al._ (2021a) R. Abbott _et al._ (LIGO Scientific, Virgo),  _Phys. Rev. X_ 11, 021053 (2021a), arXiv:2010.14527 [gr-qc] .
  * Abbott  _et al._ (2021b) R. Abbott _et al._ (LIGO Scientific, VIRGO, KAGRA), arXiv:2111.03606 [gr-qc] .
  * Bird  _et al._ (2016) S. Bird, I. Cholis, J. B. Muñoz, Y. Ali-Haïmoud, M. Kamionkowski, E. D. Kovetz, A. Raccanelli, and A. G. Riess,  _Phys. Rev. Lett._ 116, 201301 (2016), arXiv:1603.00464 [astro-ph.CO] .
  * Sasaki  _et al._ (2016) M. Sasaki, T. Suyama, T. Tanaka, and S. Yokoyama,  _Phys. Rev. Lett._ 117, 061101 (2016), [Erratum: Phys.Rev.Lett. 121, 059901 (2018)], arXiv:1603.08338 [astro-ph.CO] .
  * Clesse and García-Bellido (2017) S. Clesse and J. García-Bellido,  _Phys. Dark Univ._ 15, 142 (2017), arXiv:1603.05234 [astro-ph.CO] .
  * Wang  _et al._ (2018) S. Wang, Y.-F. Wang, Q.-G. Huang, and T. G. F. Li,  _Phys. Rev. Lett._ 120, 191102 (2018), arXiv:1610.08725 [astro-ph.CO] .
  * Bean and Magueijo (2002) R. Bean and J. Magueijo,  _Phys. Rev. D_ 66, 063505 (2002), arXiv:astro-ph/0204486 .
  * Carr  _et al._ (2021) B. Carr, K. Kohri, Y. Sendouda, and J. Yokoyama,  _Rept. Prog. Phys._ 84, 116902 (2021), arXiv:2002.12778 [astro-ph.CO] .
  * Ivanov  _et al._ (1994) P. Ivanov, P. Naselsky, and I. Novikov,  _Phys. Rev. D_ 50, 7173 (1994).
  * Garcia-Bellido  _et al._ (1996) J. Garcia-Bellido, A. D. Linde, and D. Wands,  _Phys. Rev. D_ 54, 6040 (1996), arXiv:astro-ph/9605094 .
  * Bullock and Primack (1997) J. S. Bullock and J. R. Primack,  _Phys. Rev. D_ 55, 7423 (1997), arXiv:astro-ph/9611106 .
  * Yokoyama (1997) J. Yokoyama,  _Astron. Astrophys._ 318, 673 (1997), arXiv:astro-ph/9509027 .
  * Yokoyama (1998) J. Yokoyama,  _Phys. Rev. D_ 58, 083510 (1998), arXiv:astro-ph/9802357 .
  * Kawasaki  _et al._ (1998) M. Kawasaki, N. Sugiyama, and T. Yanagida,  _Phys. Rev. D_ 57, 6050 (1998), arXiv:hep-ph/9710259 .
  * Kawasaki  _et al._ (2006) M. Kawasaki, T. Takayama, M. Yamaguchi, and J. Yokoyama,  _Phys. Rev. D_ 74, 043525 (2006), arXiv:hep-ph/0605271 .
  * Kohri  _et al._ (2008) K. Kohri, D. H. Lyth, and A. Melchiorri,  _JCAP_ 04, 038 (2008), arXiv:0711.5006 [hep-ph] .
  * Saito  _et al._ (2008) R. Saito, J. Yokoyama, and R. Nagata,  _JCAP_ 06, 024 (2008), arXiv:0804.3470 [astro-ph] .
  * Bugaev and Klimai (2008) E. Bugaev and P. Klimai,  _Phys. Rev. D_ 78, 063515 (2008), arXiv:0806.4541 [astro-ph] .
  * Alabidi and Kohri (2009) L. Alabidi and K. Kohri,  _Phys. Rev. D_ 80, 063511 (2009), arXiv:0906.1398 [astro-ph.CO] .
  * Drees and Erfani (2011) M. Drees and E. Erfani,  _JCAP_ 04, 005 (2011), arXiv:1102.2340 [hep-ph] .
  * Drees and Erfani (2012) M. Drees and E. Erfani,  _JCAP_ 01, 035 (2012), arXiv:1110.6052 [astro-ph.CO] .
  * Martin  _et al._ (2013) J. Martin, H. Motohashi, and T. Suyama,  _Phys. Rev. D_ 87, 023514 (2013), arXiv:1211.0083 [astro-ph.CO] .
  * Kohri  _et al._ (2013) K. Kohri, C.-M. Lin, and T. Matsuda,  _Phys. Rev. D_ 87, 103527 (2013), arXiv:1211.2371 [hep-ph] .
  * Kawasaki  _et al._ (2013) M. Kawasaki, N. Kitajima, and T. T. Yanagida,  _Phys. Rev. D_ 87, 063519 (2013), arXiv:1207.2550 [hep-ph] .
  * Clesse and García-Bellido (2015) S. Clesse and J. García-Bellido,  _Phys. Rev. D_ 92, 023524 (2015), arXiv:1501.07565 [astro-ph.CO] .
  * Kawasaki and Tada (2016) M. Kawasaki and Y. Tada,  _JCAP_ 08, 041 (2016), arXiv:1512.03515 [astro-ph.CO] .
  * Kawasaki  _et al._ (2016) M. Kawasaki, A. Kusenko, Y. Tada, and T. T. Yanagida,  _Phys. Rev. D_ 94, 083523 (2016), arXiv:1606.07631 [astro-ph.CO] .
  * Pi  _et al._ (2018) S. Pi, Y.-l. Zhang, Q.-G. Huang, and M. Sasaki,  _JCAP_ 05, 042 (2018), arXiv:1712.09896 [astro-ph.CO] .
  * Garcia-Bellido and Ruiz Morales (2017) J. Garcia-Bellido and E. Ruiz Morales,  _Phys. Dark Univ._ 18, 47 (2017), arXiv:1702.03901 [astro-ph.CO] .
  * Kannike  _et al._ (2017) K. Kannike, L. Marzola, M. Raidal, and H. Veermäe,  _JCAP_ 09, 020 (2017), arXiv:1705.06225 [astro-ph.CO] .
  * Germani and Prokopec (2017) C. Germani and T. Prokopec,  _Phys. Dark Univ._ 18, 6 (2017), arXiv:1706.04226 [astro-ph.CO] .
  * Ando  _et al._ (2018) K. Ando, K. Inomata, M. Kawasaki, K. Mukaida, and T. T. Yanagida,  _Phys. Rev. D_ 97, 123512 (2018), arXiv:1711.08956 [astro-ph.CO] .
  * Ezquiaga  _et al._ (2018) J. M. Ezquiaga, J. Garcia-Bellido, and E. Ruiz Morales,  _Phys. Lett. B_ 776, 345 (2018), arXiv:1705.04861 [astro-ph.CO] .
  * Motohashi and Hu (2017) H. Motohashi and W. Hu,  _Phys. Rev. D_ 96, 063503 (2017), arXiv:1706.06784 [astro-ph.CO] .
  * Di and Gong (2018) H. Di and Y. Gong,  _JCAP_ 07, 007 (2018), arXiv:1707.09578 [astro-ph.CO] .
  * Ballesteros and Taoso (2018) G. Ballesteros and M. Taoso,  _Phys. Rev. D_ 97, 023501 (2018), arXiv:1709.05565 [hep-ph] .
  * Garcia-Bellido  _et al._ (2017) J. Garcia-Bellido, M. Peloso, and C. Unal,  _JCAP_ 09, 013 (2017), arXiv:1707.02441 [astro-ph.CO] .
  * Hertzberg and Yamada (2018) M. P. Hertzberg and M. Yamada,  _Phys. Rev. D_ 97, 083509 (2018), arXiv:1712.09750 [astro-ph.CO] .
  * Inomata  _et al._ (2018) K. Inomata, M. Kawasaki, K. Mukaida, and T. T. Yanagida,  _Phys. Rev. D_ 97, 043514 (2018), arXiv:1711.06129 [astro-ph.CO] .
  * Cai  _et al._ (2018a) Y.-F. Cai, X. Tong, D.-G. Wang, and S.-F. Yan,  _Phys. Rev. Lett._ 121, 081306 (2018a), arXiv:1805.03639 [astro-ph.CO] .
  * Drees and Xu (2021) M. Drees and Y. Xu,  _Eur. Phys. J. C_ 81, 182 (2021), arXiv:1905.13581 [hep-ph] .
  * Atal  _et al._ (2019) V. Atal, J. Garriga, and A. Marcos-Caballero,  _JCAP_ 09, 073 (2019), arXiv:1905.13202 [astro-ph.CO] .
  * Atal  _et al._ (2020) V. Atal, J. Cid, A. Escrivà, and J. Garriga,  _JCAP_ 05, 022 (2020), arXiv:1908.11357 [astro-ph.CO] .
  * Mishra and Sahni (2020) S. S. Mishra and V. Sahni,  _JCAP_ 04, 007 (2020), arXiv:1911.00057 [gr-qc] .
  * Cheong  _et al._ (2021) D. Y. Cheong, S. M. Lee, and S. C. Park,  _JCAP_ 01, 032 (2021), arXiv:1912.12032 [hep-ph] .
  * Fu  _et al._ (2019) C. Fu, P. Wu, and H. Yu,  _Phys. Rev. D_ 100, 063532 (2019), arXiv:1907.05042 [astro-ph.CO] .
  * Dalianis  _et al._ (2020) I. Dalianis, S. Karydas, and E. Papantonopoulos,  _JCAP_ 06, 040 (2020), arXiv:1910.00622 [astro-ph.CO] .
  * Ashoorioon  _et al._ (2021) A. Ashoorioon, A. Rostami, and J. T. Firouzjaee,  _JHEP_ 07, 087 (2021), arXiv:1912.13326 [astro-ph.CO] .
  * Lin  _et al._ (2020) J. Lin, Q. Gao, Y. Gong, Y. Lu, C. Zhang, and F. Zhang,  _Phys. Rev. D_ 101, 103515 (2020), arXiv:2001.05909 [gr-qc] .
  * Yi  _et al._ (2021) Z. Yi, Y. Gong, B. Wang, and Z.-h. Zhu,  _Phys. Rev. D_ 103, 063535 (2021), arXiv:2007.09957 [gr-qc] .
  * Palma  _et al._ (2020) G. A. Palma, S. Sypsas, and C. Zenteno,  _Phys. Rev. Lett._ 125, 121301 (2020), arXiv:2004.06106 [astro-ph.CO] .
  * Braglia  _et al._ (2020) M. Braglia, D. K. Hazra, F. Finelli, G. F. Smoot, L. Sriramkumar, and A. A. Starobinsky,  _JCAP_ 08, 001 (2020), arXiv:2005.02895 [astro-ph.CO] .
  * Kefala  _et al._ (2021) K. Kefala, G. P. Kodaxis, I. D. Stamou, and N. Tetradis,  _Phys. Rev. D_ 104, 023506 (2021), arXiv:2010.12483 [astro-ph.CO] .
  * Ballesteros  _et al._ (2020) G. Ballesteros, J. Rey, M. Taoso, and A. Urbano,  _JCAP_ 07, 025 (2020), arXiv:2001.08220 [astro-ph.CO] .
  * Aldabergenov  _et al._ (2020) Y. Aldabergenov, A. Addazi, and S. V. Ketov,  _Eur. Phys. J. C_ 80, 917 (2020), arXiv:2006.16641 [hep-th] .
  * Aldabergenov  _et al._ (2021) Y. Aldabergenov, A. Addazi, and S. V. Ketov,  _Phys. Lett. B_ 814, 136069 (2021), arXiv:2008.10476 [hep-th] .
  * Inomata  _et al._ (2021) K. Inomata, E. McDonough, and W. Hu,  _Phys. Rev. D_ 104, 123553 (2021), arXiv:2104.03972 [astro-ph.CO] .
  * Inomata  _et al._ (2022a) K. Inomata, E. McDonough, and W. Hu,  _JCAP_ 02, 031 (2022a), arXiv:2110.14641 [astro-ph.CO] .
  * Dalianis  _et al._ (2021) I. Dalianis, G. P. Kodaxis, I. D. Stamou, N. Tetradis, and A. Tsigkas-Kouvelis,  _Phys. Rev. D_ 104, 103510 (2021), arXiv:2106.02467 [astro-ph.CO] .
  * Cai  _et al._ (2022a) Y.-F. Cai, X.-H. Ma, M. Sasaki, D.-G. Wang, and Z. Zhou,  _Phys. Lett. B_ 834, 137461 (2022a), arXiv:2112.13836 [astro-ph.CO] .
  * Lin  _et al._ (2021) J. Lin, S. Gao, Y. Gong, Y. Lu, Z. Wang, and F. Zhang, arXiv:2111.01362 [gr-qc] .
  * Kawai and Kim (2021a) S. Kawai and J. Kim,  _Phys. Rev. D_ 104, 083545 (2021a), arXiv:2108.01340 [astro-ph.CO] .
  * Zhang (2022) F. Zhang,  _Phys. Rev. D_ 105, 063539 (2022), arXiv:2112.10516 [gr-qc] .
  * Ahmed  _et al._ (2022) W. Ahmed, M. Junaid, and U. Zubair,  _Nucl. Phys. B_ 984, 115968 (2022), arXiv:2109.14838 [astro-ph.CO] .
  * Cai  _et al._ (2022b) Y.-F. Cai, X.-H. Ma, M. Sasaki, D.-G. Wang, and Z. Zhou, arXiv:2207.11910 [astro-ph.CO] .
  * Pi and Wang (2022) S. Pi and J. Wang, arXiv:2209.14183 [astro-ph.CO] .
  * Cheong  _et al._ (2022) D. Y. Cheong, K. Kohri, and S. C. Park,  _JCAP_ 10, 015 (2022), arXiv:2205.14813 [hep-ph] .
  * Kawai and Kim (2022) S. Kawai and J. Kim, arXiv:2209.15343 [astro-ph.CO] .
  * Akrami  _et al._ (2020) Y. Akrami _et al._ (Planck),  _Astron. Astrophys._ 641, A10 (2020), arXiv:1807.06211 [astro-ph.CO] .
  * Green and Malik (2001) A. M. Green and K. A. Malik,  _Phys. Rev. D_ 64, 021301 (2001), arXiv:hep-ph/0008113 .
  * Bassett and Tsujikawa (2001) B. A. Bassett and S. Tsujikawa,  _Phys. Rev. D_ 63, 123503 (2001), arXiv:hep-ph/0008328 .
  * Suyama  _et al._ (2005) T. Suyama, T. Tanaka, B. Bassett, and H. Kudoh,  _Phys. Rev. D_ 71, 063507 (2005), arXiv:hep-ph/0410247 .
  * Martin  _et al._ (2020) J. Martin, T. Papanikolaou, and V. Vennin,  _JCAP_ 01, 024 (2020), arXiv:1907.04236 [astro-ph.CO] .
  * Bezrukov and Shaposhnikov (2008) F. L. Bezrukov and M. Shaposhnikov,  _Phys. Lett. B_ 659, 703 (2008), arXiv:0710.3755 [hep-th] .
  * Bezrukov  _et al._ (2009) F. L. Bezrukov, A. Magnin, and M. Shaposhnikov,  _Phys. Lett. B_ 675, 88 (2009), arXiv:0812.4950 [hep-ph] .
  * Bezrukov  _et al._ (2011) F. Bezrukov, A. Magnin, M. Shaposhnikov, and S. Sibiryakov,  _JHEP_ 01, 016 (2011), arXiv:1008.5157 [hep-ph] .
  * Futamase and Maeda (1989) T. Futamase and K.-i. Maeda,  _Phys. Rev. D_ 39, 399 (1989).
  * Fakir and Unruh (1990) R. Fakir and W. G. Unruh,  _Phys. Rev. D_ 41, 1783 (1990).
  * Komatsu and Futamase (1999) E. Komatsu and T. Futamase,  _Phys. Rev. D_ 59, 064029 (1999), arXiv:astro-ph/9901127 .
  * Tsujikawa and Gumjudpai (2004) S. Tsujikawa and B. Gumjudpai,  _Phys. Rev. D_ 69, 123523 (2004), arXiv:astro-ph/0402185 .
  * Linde  _et al._ (2011) A. Linde, M. Noorbala, and A. Westphal,  _JCAP_ 03, 013 (2011), arXiv:1101.2652 [hep-th] .
  * Ade  _et al._ (2014) P. A. R. Ade _et al._ (Planck),  _Astron. Astrophys._ 571, A22 (2014), arXiv:1303.5082 [astro-ph.CO] .
  * Tsujikawa  _et al._ (2013) S. Tsujikawa, J. Ohashi, S. Kuroyanagi, and A. De Felice,  _Phys. Rev. D_ 88, 023529 (2013), arXiv:1305.3044 [astro-ph.CO] .
  * Hwang and Noh (2005) J.-c. Hwang and H. Noh,  _Phys. Rev. D_ 71, 063536 (2005), arXiv:gr-qc/0412126 .
  * Guo  _et al._ (2007) Z.-K. Guo, N. Ohta, and S. Tsujikawa,  _Phys. Rev. D_ 75, 023520 (2007), arXiv:hep-th/0610336 .
  * Satoh and Soda (2008) M. Satoh and J. Soda,  _JCAP_ 09, 019 (2008), arXiv:0806.4594 [astro-ph] .
  * Guo and Schwarz (2009) Z.-K. Guo and D. J. Schwarz,  _Phys. Rev. D_ 80, 063523 (2009), arXiv:0907.0427 [hep-th] .
  * Guo and Schwarz (2010) Z.-K. Guo and D. J. Schwarz,  _Phys. Rev. D_ 81, 123520 (2010), arXiv:1001.1897 [hep-th] .
  * Kawai and Kim (2021b) S. Kawai and J. Kim,  _Phys. Rev. D_ 104, 043525 (2021b), arXiv:2105.04386 [hep-ph] .
  * Kobayashi  _et al._ (2011) T. Kobayashi, M. Yamaguchi, and J. Yokoyama,  _Prog. Theor. Phys._ 126, 511 (2011), arXiv:1105.5723 [hep-th] .
  * Kase and Tsujikawa (2019) R. Kase and S. Tsujikawa,  _Int. J. Mod. Phys. D_ 28, 1942005 (2019), arXiv:1809.08735 [gr-qc] .
  * Salopek  _et al._ (1989) D. S. Salopek, J. R. Bond, and J. M. Bardeen,  _Phys. Rev. D_ 40, 1753 (1989).
  * Makino and Sasaki (1991) N. Makino and M. Sasaki,  _Prog. Theor. Phys._ 86, 103 (1991).
  * Kaiser (1995) D. I. Kaiser,  _Phys. Rev. D_ 52, 4295 (1995), arXiv:astro-ph/9408044 .
  * De Simone  _et al._ (2009) A. De Simone, M. P. Hertzberg, and F. Wilczek,  _Phys. Lett. B_ 678, 1 (2009), arXiv:0812.4946 [hep-ph] .
  * Hamada  _et al._ (2014) Y. Hamada, H. Kawai, K.-y. Oda, and S. C. Park,  _Phys. Rev. Lett._ 112, 241301 (2014), arXiv:1403.5043 [hep-ph] .
  * Hamada  _et al._ (2015) Y. Hamada, H. Kawai, K.-y. Oda, and S. C. Park,  _Phys. Rev. D_ 91, 053008 (2015), arXiv:1408.4864 [hep-ph] .
  * Bezrukov  _et al._ (2015) F. Bezrukov, J. Rubio, and M. Shaposhnikov,  _Phys. Rev. D_ 92, 083512 (2015), arXiv:1412.3811 [hep-ph] .
  * Bezrukov  _et al._ (2018) F. Bezrukov, M. Pauly, and J. Rubio,  _JCAP_ 02, 040 (2018), arXiv:1706.05007 [hep-ph] .
  * Horndeski (1974) G. W. Horndeski,  _Int. J. Theor. Phys._ 10, 363 (1974).
  * Deffayet  _et al._ (2011) C. Deffayet, X. Gao, D. A. Steer, and G. Zahariade,  _Phys. Rev. D_ 84, 064039 (2011), arXiv:1103.3260 [hep-th] .
  * Charmousis  _et al._ (2012) C. Charmousis, E. J. Copeland, A. Padilla, and P. M. Saffin,  _Phys. Rev. Lett._ 108, 051101 (2012), arXiv:1106.2000 [hep-th] .
  * De Felice and Tsujikawa (2011a) A. De Felice and S. Tsujikawa,  _JCAP_ 04, 029 (2011a), arXiv:1103.1172 [astro-ph.CO] .
  * Gross and Sloan (1987) D. J. Gross and J. H. Sloan,  _Nucl. Phys. B_ 291, 41 (1987).
  * Metsaev and Tseytlin (1987) R. R. Metsaev and A. A. Tseytlin,  _Nucl. Phys. B_ 293, 385 (1987).
  * Gasperini  _et al._ (1997) M. Gasperini, M. Maggiore, and G. Veneziano,  _Nucl. Phys. B_ 494, 315 (1997), arXiv:hep-th/9611039 .
  * Kawai  _et al._ (1998) S. Kawai, M.-a. Sakagami, and J. Soda,  _Phys. Lett. B_ 437, 284 (1998), arXiv:gr-qc/9802033 .
  * Cartier  _et al._ (2001) C. Cartier, J.-c. Hwang, and E. J. Copeland,  _Phys. Rev. D_ 64, 103504 (2001), arXiv:astro-ph/0106197 .
  * Calcagni  _et al._ (2005) G. Calcagni, S. Tsujikawa, and M. Sami,  _Class. Quant. Grav._ 22, 3977 (2005), arXiv:hep-th/0505193 .
  * Khan and Yogesh (2022) H. A. Khan and Yogesh,  _Phys. Rev. D_ 105, 063526 (2022), arXiv:2201.06439 [astro-ph.CO] .
  * Green  _et al._ (2004) A. M. Green, A. R. Liddle, K. A. Malik, and M. Sasaki,  _Phys. Rev. D_ 70, 041502 (2004), arXiv:astro-ph/0403181 .
  * Young  _et al._ (2014) S. Young, C. T. Byrnes, and M. Sasaki,  _JCAP_ 07, 045 (2014), arXiv:1405.7023 [gr-qc] .
  * Harada  _et al._ (2013) T. Harada, C.-M. Yoo, and K. Kohri,  _Phys. Rev. D_ 88, 084051 (2013), [Erratum: Phys.Rev.D 89, 029903 (2014)], arXiv:1309.4201 [astro-ph.CO] .
  * Germani and Musco (2019) C. Germani and I. Musco,  _Phys. Rev. Lett._ 122, 141302 (2019), arXiv:1805.04087 [astro-ph.CO] .
  * Kavanagh (2019) B. J. Kavanagh, “bradkav/pbhbounds: Release version,”  (2019).
  * Carr (1975) B. J. Carr,  _Astrophys. J._ 201, 1 (1975).
  * Kohri and Terada (2018) K. Kohri and T. Terada,  _Phys. Rev. D_ 97, 123532 (2018), arXiv:1804.08577 [gr-qc] .
  * Papanikolaou  _et al._ (2021) T. Papanikolaou, V. Vennin, and D. Langlois,  _JCAP_ 03, 053 (2021), arXiv:2010.11573 [astro-ph.CO] .
  * Namjoo  _et al._ (2013) M. H. Namjoo, H. Firouzjahi, and M. Sasaki,  _EPL_ 101, 39001 (2013), arXiv:1210.3692 [astro-ph.CO] .
  * Franciolini  _et al._ (2018) G. Franciolini, A. Kehagias, S. Matarrese, and A. Riotto,  _JCAP_ 03, 016 (2018), arXiv:1801.09415 [astro-ph.CO] .
  * Cai  _et al._ (2018b) Y.-F. Cai, X. Chen, M. H. Namjoo, M. Sasaki, D.-G. Wang, and Z. Wang,  _JCAP_ 05, 012 (2018b), arXiv:1712.09998 [astro-ph.CO] .
  * Atal and Germani (2019) V. Atal and C. Germani,  _Phys. Dark Univ._ 24, 100275 (2019), arXiv:1811.07857 [astro-ph.CO] .
  * Passaglia  _et al._ (2019) S. Passaglia, W. Hu, and H. Motohashi,  _Phys. Rev. D_ 99, 043536 (2019), arXiv:1812.08243 [astro-ph.CO] .
  * Taoso and Urbano (2021) M. Taoso and A. Urbano,  _JCAP_ 08, 016 (2021), arXiv:2102.03610 [astro-ph.CO] .
  * Biagetti  _et al._ (2021) M. Biagetti, V. De Luca, G. Franciolini, A. Kehagias, and A. Riotto,  _Phys. Lett. B_ 820, 136602 (2021), arXiv:2105.07810 [astro-ph.CO] .
  * Davies  _et al._ (2022) M. W. Davies, P. Carrilho, and D. J. Mulryne,  _JCAP_ 06, 019 (2022), arXiv:2110.08189 [astro-ph.CO] .
  * Kristiano and Yokoyama (2022) J. Kristiano and J. Yokoyama, arXiv:2211.03395 [hep-th] .
  * Inomata  _et al._ (2022b) K. Inomata, M. Braglia, and X. Chen, arXiv:2211.02586 [astro-ph.CO] .
  * De Felice and Tsujikawa (2011b) A. De Felice and S. Tsujikawa,  _Phys. Rev. D_ 84, 083504 (2011b), arXiv:1107.3917 [gr-qc] .
  * De Felice  _et al._ (2011) A. De Felice, T. Kobayashi, and S. Tsujikawa,  _Phys. Lett. B_ 706, 123 (2011), arXiv:1108.4242 [gr-qc] .