August 2021 IPMU21-0034

Exploring the parameter space of modified supergravity
for double inflation and primordial black hole formation

Ryotaro Ishikawa a and Sergei V. Ketov a,b,c

a Department of Physics, Tokyo Metropolitan University

1-1 Minami-ohsawa, Hachioji-shi, Tokyo 192-0397, Japan

b Research School of High-Energy Physics, Tomsk Polytechnic University

2a Lenin Avenue, Tomsk 634028, Russian Federation

c Kavli Institute for the Physics and Mathematics of the Universe (WPI)

The University of Tokyo Institutes for Advanced Study, Kashiwa 277-8583, Japan

<EMAIL_ADDRESS><EMAIL_ADDRESS>
Abstract

We study the parameter space of the effective (with two scalars) models of
cosmological inflation and primordial black hole (PBH) formation in the
modified $(R+R^{2})$ supergravity. Our models describe double inflation, whose
first stage is driven by Starobinsky’s scalaron coming from the $R^{2}$
gravity, and whose second stage is driven by another scalar belonging to the
supergravity multiplet. The ultra-slow-roll regime between the two stages
leads a large peak (enhancement) in the power spectrum of scalar
perturbations, which results in efficient PBH formation. Both inflation and
PBH formation are generic in our models, while those PBH can account for a
significant part or the whole of dark matter. Some of the earlier proposed
models in the same class are in tension (over $3\sigma$) with the observed
value of the scalar tilt $n_{s}$, so that we study more general models with
more parameters, and investigate the dependence of the cosmological tilts
$(n_{s},r)$ and the scalar power spectrum enhancement upon the parameters. The
PBH masses and their density fraction (as part of dark matter) are also
calculated. A good agreement (between $2\sigma$ and $3\sigma$) with the
observed value of $n_{s}$ requires fine tuning of the parameters, and it is
only realized in the so-called $\delta$-models. Our models offer the
(super)gravitational origin of inflation, PBH and dark matter together, and
may be confirmed or falsified by future precision measurements of the cosmic
microwave background radiation and PBH-induced gravitational waves.

##  1 Introduction

Symmetry principles play the very important role in fundamental physics,
especially well beyond the Standard Model (SM) of elementary particles, where
they provide advanced guidance to the phenomenological model building.
Supersymmetry (SUSY) is one of such fundamental symmetry principles that was
proposed many years ago as the nontrivial extension of the Poincaré symmetry.
No SUSY signature at the large hadron collider (LHC) rules out low-scale SUSY
but does not exclude high-scale SUSY well beyond 10 TeV. Local SUSY implies
the invariance under general coordinate transformations, so that supergravity
as the theory of local SUSY is more fundamental than general relativity. The
scale of cosmological inflation is expected at $10^{10}$ TeV that is well
beyond the electro-weak scale, where quantum gravity corrections are expected
to be significant also. Supergravity can serve as a bridge between classical
and quantum gravity, which makes it indispensable at the super-high energy
scales close to the GUT scale, see e.g., Refs. [1, 2] for the recent reviews.

The standard (concordance) model of cosmology (SCM) assumes cold dark matter
(CDM) and the positive cosmological constant as the dark energy (DE). However,
it does not tell us about physics behind CDM and DE. It is often assumed that
CDM is composed of unknown (electrically neutral and stable or meta-stable)
particles, though it is also possible that part or the whole CDM may also be
composed of primordial black holes (PBH) formed in the early Universe during
inflation [3, 4, 5, 6, 7]. We adopt the second possibility in this paper (PBH
as DM), see also Refs. [8, 9, 10, 11] for the recent reviews.

There is the well-known problem with the cosmological constant because its
interpretation as the vacuum energy contradicts quantum theory. A positive
cosmological constant implies a de Sitter vacuum that is incompatible with the
S-matrix approach and, therefore, is inconsistent with quantum gravity [12].
In particular, de Sitter vacua are believed to be forbidden in string theory
also [13].

Yet another well known indication on the necessity to go beyond the SCM is the
Hubble tension, namely, a disagreement between local measurements of Hubble
parameter from supernova and lensing, on the one side, and its value inferred
from the SCM fit to the cosmic microwave background (CMB) radiation [14].

To avoid (and mimic) a positive cosmological constant in the classical
description of gravity, one can employ the modified $f(R)$ gravity theory,
where DE can be effectively described by a non-trivial function $f(R)$ of the
(Ricci) scalar curvature in the gravitational action, in agreement with
gravitational tests in Solar system [15, 16, 17], see Ref. [18] for a
comprehensive review.

Cosmological inflation in the early Universe can also be successfully
described in the modified $f(R)$ gravity by the celebrated Starobinsky model
of $(R+R^{2})$ gravity [19], see also Ref. [20] for a recent review of
Starobinsky’s inflation. The Starobinsky inflation model is more than the best
fit to the current CMB observations. It is geometrical because it relies on
gravitational interactions, has no free parameters, and is driven by the
scale-invariant $R^{2}$ term in the action. Hence, Starobinsky’s inflaton
(called scalaron) has the clear gravitational origin as the physical degree of
freedom in the higher-derivative gravity, while it can be identified with the
Nambu-Goldstone boson related to spontaneous breaking of the scale invariance
[21].

However, the Starobinsky model has limited applicability because the CMB data
is just a small window into high-energy physics, while there are no
observational constraints for the scales both well beyond and well under the
(expected) inflationary scale of $10^{13}$ GeV, when being measured by the
Hubble parameter. There is also no good reason for the absence of other
scalars during inflation. Therefore, it is of importance to figure out
possible interactions of scalaron with other particles before and during
reheating. Specifying those scalaron interactions can be done either by a
trial-and-error procedure or by postulating some fundamental (symmetry)
principles. In this paper we assume Supersymmetry (SUSY) as the guiding
fundamental symmetry principle, which means local SUSY or supergravity in the
context of gravity. Despite the current absence of experimental confirmation,
SUSY is still one of the leading candidates for new physics beyond the SM. Of
course, SUSY has to be (spontaneously) broken at the energies available for
observations, see e.g., Refs. [22, 23] for specific mechanisms of spontaneous
SUSY breaking after inflation in supergravity.

Therefore, the natural questions arise: (i) whether the modified supergravity
can successfully describe inflation and PBH formation in the early Universe,
and (ii) how much fine tuning is needed for those purposes?

There exist many supergravity extensions of the Starobinsky inflation, see
e.g., Ref. [20] for a brief review of the existing approaches to embedding the
Starobinsky model into supergravity, and Ref. [24] for inclusion of DM and DE
too. Supergravity generically leads to multi-field inflation with several
scalar fields involved, whereas the Starobinsky inflation model is an example
of single-field inflation. Most of the known embeddings of inflation into
supergravity have a single inflaton and suppress other scalars. In this paper
we employ the generalized (or modified) Starobinsky-type supergravity [21, 24,
25, 26, 27] in order to describe multi-field inflation. We minimize the
particle spectrum by keeping all fields in the same (irreducible) supergravity
multiplet, while their interactions are severely restricted by local SUSY. We
give a supergravity realization of two-field double inflation with a sharp
turn of the inflationary trajectory, where one scalar is given by
Starobinsky’s inflaton and another scalar is its scalar superpartner. We show
that it leads to PBH generation and DM genesis after Starobinsky inflation.
Our models in this paper do not have a cosmological constant and have
Minkowski vacua, without addressing DE. However, they allow further extensions
that may lead to a viable DE description also, which is beyond the scope of
this investigation.

Our paper is organized as follows. In Sec. 2 we introduce our setup by
defining the modified Starobinsky-type supergravity and deriving the effective
two-field inflationary models, which coincides with the setup adopted in Refs.
[28, 29]. In Sec. 3 we review the two special cases already studied in Ref.
[28], having one parameter less, and called the $\gamma$ and $\delta$ models,
respectively. The main purpose of this investigation is to study the models
where both parameters $\gamma$ and $\delta$ are non-vanishing, as well as even
more general models beyond Ref. [28]. The significance of the new cases stems
from the observation that the $\gamma$-models with one less parameter are in
tension (over $3\sigma$) with the observed value of the CMB scalar spectral
index $n_{s}$. It was conjectured in Ref. [28] that one may get a much better
agreement in more general models with both non-vanishing $\gamma$ and $\delta$
parameters, or by adding more parameters. We present our new findings in Sec.
4 that contains our main results. Our Conclusion is given by Sec. 5.

The modified (Starobinsky) supergravity is used to generate the scalar
potential, instead of postulating it ad hoc. Once the scalar potential is
derived, the standard techniques are applied to compute the power spectra of
scalar and tensor perturbations and the inflationary observables. In the case
of two-field inflation, those techniques are described at length in Refs. [28,
30, 31], so that we skip their details here and focus on our new results.

##  2 Modified Starobinsky supergravity and the effective two-field models of
inflation and PBH formation

We use the standard notation of Ref. [32] in curved superspace of the (old-
minimal) supergravity. It guarantees that our starting actions are manifestly
$N=1$ locally supersymmetric in four spacetime dimensions. As regards the
supersymmetry transformation laws of the field components, see Ref. [32].

The chiral superspace Lagrangian of the chiral superfields $\mathbf{\Phi}^{i}$
coupled to supergravity reads (we take the reduced Planck mass $M_{\rm Pl}=1$
for simplicity)

${\cal L}=\int d^{2}\Theta 2{\cal
E}\left[{\textstyle{3\over\vphantom{2}\smash{\raise
0.60275pt\hbox{$\scriptstyle{8}$}}}}(\mkern 1.5mu\overline{\mkern-1.5mu\cal
D\mkern-1.5mu}\mkern 1.5mu^{2}-8{\cal R})e^{-K(\mathbf{\Phi}^{i},\mkern
1.5mu\overline{\mkern-1.5mu\mathbf{\Phi}\mkern-1.5mu}\mkern
1.5mu^{i})/3}+W(\mathbf{\Phi}^{i})\right]+{\rm h.c.}~{},$ (1)

where we have introduced the chiral density superfield $\cal E$, the chiral
curvature superfield $\cal R$, the superspace covariant derivatives ${\cal
D}_{\alpha},\mkern 1.5mu\overline{\mkern-1.5mu\cal D\mkern-1.5mu}\mkern
1.5mu_{\dot{\alpha}}$ , with ${\cal D}^{2}\equiv{\cal D}^{\alpha}{\cal
D}_{\alpha}$ and $\mkern 1.5mu\overline{\mkern-1.5mu\cal D\mkern-1.5mu}\mkern
1.5mu^{2}\equiv\mkern 1.5mu\overline{\mkern-1.5mu\cal D\mkern-1.5mu}\mkern
1.5mu_{\dot{\alpha}}\mkern 1.5mu\overline{\mkern-1.5mu\cal
D\mkern-1.5mu}\mkern 1.5mu^{\dot{\alpha}}$, the Kähler potential
$K(\Phi^{i},\mkern 1.5mu\overline{\mkern-1.5mu\Phi\mkern-1.5mu}\mkern
1.5mu^{i})$ and the superpotential $W(\Phi^{i})$. The $K$ and $W$ define the
model and uniquely determine its scalar sector.

Eliminating the auxiliary fields and going to Einstein frame via field
redefinitions yield the bosonic part of the Lagrangian (1) in the form

$e^{-1}{\cal L}={\textstyle{1\over\vphantom{2}\smash{\raise
0.60275pt\hbox{$\scriptstyle{2}$}}}}R-K_{i\bar{j}}\partial_{m}\Phi^{i}\partial^{m}\mkern
1.5mu\overline{\mkern-1.5mu\Phi\mkern-1.5mu}\mkern
1.5mu^{j}-e^{K}\left(K^{i\bar{j}}D_{i}WD_{\bar{j}}\mkern
1.5mu\overline{\mkern-1.5muW\mkern-1.5mu}\mkern 1.5mu-3|W|^{2}\right)~{},$ (2)

where we have used the same notation for the chiral superfields and their
leading (scalar) field components, together with

$\displaystyle
K_{i\bar{j}}\equiv{{\partial^{2}K}\over{\partial\Phi^{i}\partial\mkern
1.5mu\overline{\mkern-1.5mu\Phi\mkern-1.5mu}\mkern 1.5mu^{j}}}~{},\quad
K^{i\bar{j}}\equiv K^{-1}_{i\bar{j}}~{},\quad{\rm and}\quad
D_{i}W\equiv{{\partial W}\over{\partial\Phi^{i}}}+W{{\partial
K}\over{\partial\Phi^{i}}}~{}~{}.$ (3)

The expansions of the supergravity chiral superfields $\cal E$ and $\cal R$
with respect to the chiral anti-commuting coordinates $\Theta_{\alpha}$ define
their field components as follows:

$\displaystyle 2{\cal E}$ $\displaystyle=$ $\displaystyle
e\left[1+i\Theta\sigma^{m}\mkern
1.5mu\overline{\mkern-1.5mu\psi\mkern-1.5mu}\mkern
1.5mu_{m}+\Theta^{2}(6\mkern 1.5mu\overline{\mkern-1.5muX\mkern-1.5mu}\mkern
1.5mu-\mkern 1.5mu\overline{\mkern-1.5mu\psi\mkern-1.5mu}\mkern
1.5mu_{m}\mkern 1.5mu\overline{\mkern-1.5mu\sigma\mkern-1.5mu}\mkern
1.5mu^{mn}\mkern 1.5mu\overline{\mkern-1.5mu\psi\mkern-1.5mu}\mkern
1.5mu_{n})\right]~{},$ (4) $\displaystyle{\cal R}$ $\displaystyle=$
$\displaystyle X+\Theta\left(-{\textstyle{1\over\vphantom{2}\smash{\raise
0.60275pt\hbox{$\scriptstyle{6}$}}}}\sigma^{m}\mkern
1.5mu\overline{\mkern-1.5mu\sigma\mkern-1.5mu}\mkern
1.5mu^{n}\psi_{mn}-i\sigma^{m}\mkern
1.5mu\overline{\mkern-1.5mu\psi\mkern-1.5mu}\mkern
1.5mu_{m}X-{\textstyle{i\over\vphantom{2}\smash{\raise
0.60275pt\hbox{$\scriptstyle{6}$}}}}\psi_{m}b^{m}\right)+$ (5)
$\displaystyle+\Theta^{2}\left(-{\textstyle{1\over\vphantom{2}\smash{\raise
0.60275pt\hbox{$\scriptstyle{12}$}}}}R-{\textstyle{i\over\vphantom{2}\smash{\raise
0.60275pt\hbox{$\scriptstyle{6}$}}}}\mkern
1.5mu\overline{\mkern-1.5mu\psi\mkern-1.5mu}\mkern 1.5mu^{m}\mkern
1.5mu\overline{\mkern-1.5mu\sigma\mkern-1.5mu}\mkern
1.5mu^{n}\psi_{mn}-4X\mkern 1.5mu\overline{\mkern-1.5muX\mkern-1.5mu}\mkern
1.5mu-{\textstyle{1\over\vphantom{2}\smash{\raise
0.60275pt\hbox{$\scriptstyle{18}$}}}}b_{m}b^{m}+{\textstyle{i\over\vphantom{2}\smash{\raise
0.60275pt\hbox{$\scriptstyle{6}$}}}}\nabla_{m}b^{m}+\right.$
$\displaystyle+\left.{\textstyle{1\over\vphantom{2}\smash{\raise
0.60275pt\hbox{$\scriptstyle{2}$}}}}\mkern
1.5mu\overline{\mkern-1.5mu\psi\mkern-1.5mu}\mkern 1.5mu_{m}\mkern
1.5mu\overline{\mkern-1.5mu\psi\mkern-1.5mu}\mkern
1.5mu^{m}X+{\textstyle{1\over\vphantom{2}\smash{\raise
0.60275pt\hbox{$\scriptstyle{12}$}}}}\psi_{m}\sigma^{m}\mkern
1.5mu\overline{\mkern-1.5mu\psi\mkern-1.5mu}\mkern
1.5mu_{n}b^{n}-{\textstyle{1\over\vphantom{2}\smash{\raise
0.60275pt\hbox{$\scriptstyle{48}$}}}}\varepsilon^{abcd}(\mkern
1.5mu\overline{\mkern-1.5mu\psi\mkern-1.5mu}\mkern 1.5mu_{a}\mkern
1.5mu\overline{\mkern-1.5mu\sigma\mkern-1.5mu}\mkern
1.5mu_{b}\psi_{cd}+\psi_{a}\sigma_{b}\mkern
1.5mu\overline{\mkern-1.5mu\psi\mkern-1.5mu}\mkern 1.5mu_{cd})\right)~{},$

where $e\equiv{\rm det}(e^{a}_{m})$,
$\psi_{mn}\equiv\tilde{D}_{m}\psi_{n}-\tilde{D}_{n}\psi_{m}$ and
$\tilde{D}_{m}\psi_{n}\equiv(\partial_{m}+\omega_{m})\psi_{n}$.

The supergravity multiplet includes vierbein $e^{a}_{m}$, gravitino
$\psi_{m}$, the real vector field $b_{m}$ and the complex scalar $X$. The
chiral superfield ${\cal E}$ can be seen as the SUSY extension of the
spacetime density $e=\sqrt{-g}$, and the chiral superfield ${\cal R}$ can be
seen as the SUSY extension of the (Ricci) scalar curvature $R$. The real
vector $b_{m}$ and the complex scalar $X$ are known in the supergravity
literature as the (old-minimal) set of the ”auxiliary” fields needed to
complete the supergravity multiplet off-shell. In the higher-derivative
modified supergravity (see below), those ”auxiliary” fields become the
physical ones, i.e. they are dynamical or propagating.

The modified Starobinsky-type supergravity action in curved (full) superspace
is defined by [33, 34]

$S=\int d^{4}xd^{4}\theta E^{-1}N(\mathcal{R},\bar{\mathcal{R}})+\left[\int
d^{4}xd^{2}\Theta 2\mathcal{E}{\cal F}(\mathcal{R})+h.c\right]~{}~{}.$ (6)

The equivalent action in the curved chiral superspace reads

${\cal L}=\int d^{2}\Theta 2{\cal
E}\left[-{\textstyle{1\over\vphantom{2}\smash{\raise
0.60275pt\hbox{$\scriptstyle{8}$}}}}(\mkern 1.5mu\overline{\mkern-1.5mu\cal
D\mkern-1.5mu}\mkern 1.5mu^{2}-8{\cal R})N({\cal R},\mkern
1.5mu\overline{\mkern-1.5mu\cal R\mkern-1.5mu}\mkern 1.5mu)+{\cal F}({\cal
R})\right]+{\rm h.c.}$ (7)

This action is governed by two arbitrary potentials $N$ and ${\cal F}$ that
are similar to $K$ and $F$ in Eq. (1), respectively. Actually, the ${\cal
F}$-type term in Eq. (6) can be included into the $N$-type term, unless the
former is a constant. However, we prefer to distinguish the two structures
because they represent two different deformations of the simplest supergravity
action. It is worth mentioning that, unlike the chiral and anti-chiral
(matter) superfields in Eq. (1), the arguments of the functions $N$ and ${\cal
F}$ are the supergravity superfields, ${\cal R}$ and $\mkern
1.5mu\overline{\mkern-1.5mu\cal R\mkern-1.5mu}\mkern 1.5mu$, that are not
unconstrained chiral superfields but satisfy the supergravity constraints in
curved superspace. The chiral superfield $\mathcal{R}$ has the scalar
curvature $R$ as its field component at $\Theta^{2}$, so that no higher powers
of $R$ (beyond the linear and quadratic terms) appear in the action (6).
Therefore, supergravity apparently distinguishes the $(R+R^{2})$ model amongst
all the $f(R)$ gravity models! The higher powers of $R$ can, nevertheless,
appear when we allow the superspace derivatives of the supergravity
superfields ${\cal R}$ and $\mkern 1.5mu\overline{\mkern-1.5mu\cal
R\mkern-1.5mu}\mkern 1.5mu$ in the arguments of $N$ and ${\cal F}$. However,
it would lead to the appearance of the spacetime derivatives of the scalar
curvature in the Lagrangian and more physical scalars.

The action (6) can be transformed into the (dual) standard matter-coupled
supergravity action of the type (1) in terms of two chiral matter superfields
$T$ and $S$ in the manifestly supersymmetric way, whose Kähler potential and
the superpotential are related to the input potentials $N$ and $F$ as follows
[33, 34, 35]:

$K=-3\log\left[T+\mkern 1.5mu\overline{\mkern-1.5muT\mkern-1.5mu}\mkern
1.5mu-{\textstyle{1\over\vphantom{2}\smash{\raise
0.60275pt\hbox{$\scriptstyle{3}$}}}}N(S,\mkern
1.5mu\overline{\mkern-1.5muS\mkern-1.5mu}\mkern 1.5mu)\right]\quad{\rm
and}\quad W=3ST+{\cal F}(S)~{}.$ (8)

This Kähler potential is an example of the no-scale supergravity that
naturally arises in compactifications of heterotic strings [36]. Unlike the
non-SUSY duality between $f(R)$ gravity and scalar-tensor gravity, there is no
inverse transformation that could relate any choice of $K$ and $W$ to the
potentials $N$ and ${\cal F}$ of the generalized Starobinsky-type
supergravity. The particular form (8) distinguishes the induced Kähler
potentials $K$ and the superpotentials $W$ of the supergravitational origin
(from the modified supergravity).

Let us expand the functions $N$ and $\cal F$ in Taylor series and keep a few
leading terms as follows:

$\displaystyle N={{12}\over{M^{2}}}|{\cal
R}|^{2}-{{72}\over{M^{4}}}\zeta|{\cal R}|^{4}-{{768}\over{M^{6}}}\gamma|{\cal
R}|^{6}~{}~{},$ (9) $\displaystyle{\cal F}=-3{\cal
R}+{{3\sqrt{6}}\over{M}}\delta{\cal R}^{2}~{}~{},$ (10)

where we have introduced three dimensionless parameters $\zeta$, $\gamma$ and
$\delta$, and the Starobinsky mass $M$. This supergravity model reduces to the
standard (pure) supergravity in the very special case of $N=0$ and ${\cal
F}=-3{\cal R}$. In particular, the first term in Eq. (10) leads to the
Einstein-Hilbert term ${{{1}\over{2}}}eR$ in the Lagrangian, according to Eqs.
(4) and (5). Similarly, the first term in the $N$-potential (9) is needed to
generate the $R^{2}$ term in the Lagrangian with the Starobinsky mass
parameter $M$. We recall that the value of $M\sim 10^{-5}$ is fixed by CMB
observations. The second term in (9) is needed for stabilization of the
Starobinsky supergravity, see Ref. [27] for details. The other $\gamma$ and
$\delta$ terms represent further deformations of that Starobinsky-type
modified supergravity [28], whose impact on PBH production is reviewed in the
next Section.

Let us also ignore the vector field $b_{m}$ and the angular part of the scalar
field ${\cal R}|=X$ for simplicity (the scalar potenial is flat in the angular
direction), and rescale the real $X$ as

$X={{M\sigma}\over{\sqrt{24}}}~{}.$ (11)

The dual matter-coupled Einstein supergravity is described by Eq. (8) with the
functions

$\displaystyle N(S,\mkern 1.5mu\overline{\mkern-1.5muS\mkern-1.5mu}\mkern
1.5mu)$
$\displaystyle=3\left(|S|^{2}-{\textstyle{3\over\vphantom{2}\smash{\raise
0.60275pt\hbox{$\scriptstyle{2}$}}}}\zeta|S|^{4}-4\gamma|S|^{6}\right)~{},$
(12) $\displaystyle{\cal F}(S)$
$\displaystyle=3MS\left({\textstyle{\sqrt{6}\over\vphantom{2}\smash{\raise
0.60275pt\hbox{$\scriptstyle{4}$}}}}\delta
S-{\textstyle{1\over\vphantom{2}\smash{\raise
0.60275pt\hbox{$\scriptstyle{2}$}}}}\right)~{},$ (13)

after rescaling ${\cal R}=MS/2$. Equation (11) implies $S=\sigma/\sqrt{6}$.
The canonical scalaron $\varphi$ in the dual picture is given by

$e^{\sqrt{{\textstyle{2\over\vphantom{2}\smash{\raise
0.60275pt\hbox{$\scriptstyle{3}$}}}}}\varphi}=T+\mkern
1.5mu\overline{\mkern-1.5muT\mkern-1.5mu}\mkern
1.5mu-{\textstyle{1\over\vphantom{2}\smash{\raise
0.60275pt\hbox{$\scriptstyle{3}$}}}}N(S,\mkern
1.5mu\overline{\mkern-1.5muS\mkern-1.5mu}\mkern 1.5mu)~{}.$ (14)

It is straightforward (after a long calculation) to get the scalar part of the
Lagrangian in Einstein frame. It takes the form of two scalar fields coupled
to gravity as follows [28]:

$e^{-1}{\cal
L}={{1}\over{2}}R-{{1}\over{2}}(\partial\varphi)^{2}-{{3M^{2}}\over{2}}Be^{-\sqrt{{\textstyle{2\over\vphantom{2}\smash{\raise
0.60275pt\hbox{$\scriptstyle{3}$}}}}}\varphi}(\partial\sigma)^{2}-{{1}\over{4B}}\left(1-Ae^{-\sqrt{{\textstyle{2\over\vphantom{2}\smash{\raise
0.60275pt\hbox{$\scriptstyle{3}$}}}}}\varphi}\right)^{2}-e^{-2\sqrt{{\textstyle{2\over\vphantom{2}\smash{\raise
0.60275pt\hbox{$\scriptstyle{3}$}}}}}\varphi}U~{},$ (15)

where the functions $A,B,U$ are given by

$\displaystyle A$
$\displaystyle=1-\delta\sigma+{{1}\over{6}}\sigma^{2}-{{11}\over{24}}\zeta\sigma^{4}-{{29}\over{54}}\gamma\sigma^{6}~{},$
$\displaystyle B$
$\displaystyle={{1}\over{3M^{2}}}(1-\zeta\sigma^{2}-\gamma\sigma^{4})~{},$
(16) $\displaystyle U$
$\displaystyle={{M^{2}}\over{2}}\sigma^{2}\left(1+{{1}\over{2}}\delta\sigma-{{1}\over{6}}\sigma^{2}+{\textstyle{3\over\vphantom{2}\smash{\raise
0.60275pt\hbox{$\scriptstyle{8}$}}}}\zeta\sigma^{4}+{\textstyle{25\over\vphantom{2}\smash{\raise
0.60275pt\hbox{$\scriptstyle{54}$}}}}\gamma\sigma^{6}\right)~{}.$

Both scalars, scalaron $\varphi$ and extra scalar $\sigma$, have the
supergravitational origin and are related via local SUSY transformations by
construction. Their kinetic terms and the scalar potential are determined by
Eqs. (15) and (27). In particular, they have the same masses. One may worry
about $\sigma$ to become a ghost due to the negative signs at some terms in
the second equation (27). Actually, the infinite wall in the scalar potential
prevents $\sigma$ from obtaining the values leading to the wrong sign of its
kinetic term. Thus, there are no ghosts.

It is also straightforward (after a long calculation) to derive the equations
of motion of our two-field model in the FLRW universe, when keeping only time
dependence of the fields. One finds [28]

$\displaystyle\ddot{\varphi}+3H\dot{\varphi}+{{1}\over{\sqrt{6}}}(1-\zeta\sigma^{2}-\gamma\sigma^{4})e^{-\sqrt{{\textstyle{2\over\vphantom{2}\smash{\raise
0.60275pt\hbox{$\scriptstyle{3}$}}}}}\varphi}\dot{\sigma}^{2}+\partial_{\varphi}V=0~{},$
(17)
$\displaystyle\ddot{\sigma}+3H\dot{\sigma}-{{\zeta\sigma+2\gamma\sigma^{3}}\over{1-\zeta\sigma^{2}-\gamma\sigma^{4}}}\dot{\sigma}^{2}-\sqrt{{{2}\over{3}}}\dot{\varphi}\dot{\sigma}+{{e^{\sqrt{{\textstyle{2\over\vphantom{2}\smash{\raise
0.60275pt\hbox{$\scriptstyle{3}$}}}}}\varphi}}\over{1-\zeta\sigma^{2}-\gamma\sigma^{4}}}\partial_{\sigma}V=0~{},$
(18) $\displaystyle{\textstyle{1\over\vphantom{2}\smash{\raise
0.60275pt\hbox{$\scriptstyle{2}$}}}}\dot{\varphi}^{2}+{\textstyle{1\over\vphantom{2}\smash{\raise
0.60275pt\hbox{$\scriptstyle{2}$}}}}(1-\zeta\sigma^{2}-\gamma\sigma^{4})e^{-\sqrt{{\textstyle{2\over\vphantom{2}\smash{\raise
0.60275pt\hbox{$\scriptstyle{3}$}}}}}\varphi}\dot{\sigma}^{2}+\dot{H}=0~{},$
(19)

where the last equation describes the kinetic energy balance. These equations
are to be supplemented by the Friedman equation for completeness,

$V=3H^{2}+\dot{H}~{},$ (20)

where the full scalar potential $V$ is given by the last two terms in Eq.
(15).

The main cosmological parameters of inflation are given by the scalar tilt
$n_{s}$ and the tensor-to-scalar ratio $r$, whose observational values are
constrained by the Planck measurements of CMB as [37]

$n_{s}=0.9649\pm 0.0042\quad({\rm 68\%CL})\qquad{\rm and}\qquad
r<0.056\quad({\rm 95\%CL})~{}.$ (21)

##  3 Two special cases

In this Section we review the results of Ref. [28] for two special families of
our models with one less free parameter, and then compare them. It serves as
the pre-requisite for our new investigation in the next Section. Those special
cases are the one with $\delta=0$ (called the $\gamma$ models) and another one
with $\gamma=0$ (called the $\delta$ models). Simultaneously, this Section
introduces our methods of investigation of PBH production, which are used in
the next Section also. More technical details can be found in Refs. [30, 31].

###  3.1 The $\gamma$ models

Let us choose the case with the parameters $\gamma=1$ and $\zeta=-1.7774$, as
a representative of the $\gamma$ models ($\delta=0$) The scalar potential for
$\varphi\gg 1$ in Fig. 1 has two valleys at $\sigma\neq 0$, and a single
Minkowski minimum at $\sigma=\varphi=0$. The first slow-roll (SR) Starobinsky-
like inflation goes along either of the valleys. The inflationary trajectory
has a sharp turn by passing through one of the near-inflection points followed
by the second inflation stage toward the Minkowski minimum. In the very short
so-called ”ultra-slow-roll” (USR) regime between the two stages of inflation,
the scalar field(s) roll down the potential faster than in the SR regime [38].

Figure 1: The scalar potential $V/M^{2}$ for $\delta=0$, $\gamma=1$ and
$\zeta=-1.7774$.

The numerical solutions to the equations of motion are plotted in Fig. 2. The
total number of e-folds is set to $\Delta N=60$, and the end of the first
stage of inflation is defined by the time when the SR parameter $\eta$ first
reaches $1$, see Fig. 2(e). The USR period $\epsilon_{\rm USR}\ll\epsilon_{\rm
SR}$ in Fig. 2(e) between the two stages of inflation leads to a significant
enhancement of the scalar power spectrum, see Fig. 3. Inflation ends when
$\epsilon=1$. The first stage lasts $\Delta N_{1}\approx 50$ e-folds, whereas
the second stage lasts for $\Delta N_{2}\approx 10$. The length of the second
stage is controlled by the parameter $\zeta$ at a given $\gamma$.

(a)

(b)

(c)

(d)

(e)

Figure 2: (a) The solution to the field equations (17) and (18) with the
initial conditions $\varphi(0)=6,\sigma(0)=0.1$, the vanishing initial
velocities, and the parameters $\delta=0$, $\gamma=1$ and $\zeta=-1.7774$. The
blue shaded region represents the first stage of inflation, and the green
shaded region represents the second stage of inflation. (b) The trajectory of
the solution. (c) The corresponding Hubble function. (d) The e-folds. (e) The
SR parameters $\epsilon$ (red) and $\eta$ (blue).

The power spectrum of curvature perturbations is numerically computed at fixed
$\Delta N_{2}$ by using the standard transport method [39, 40] with the
Mathematica package [41], around the pivot scale $k_{*}$ that leaves the
horizon at the end of the first stage (we call this scale $k_{\Delta N_{2}}$).
The inflaton mass is chosen to be $0.5\times 10^{-5}M_{\rm Pl}$ by requiring
$P_{\zeta}\approx 2\times 10^{-9}$ for the mode $k$ that exits the horizon 60
e-folds before the end of inflation (we call it $k_{60}$). The results for
various values of $\gamma$ are shown in Fig. 3. The values of the parameters
are collected in Table 1, where $\zeta$ is tuned to satisfy $\Delta
N_{2}\approx 10$.

Figure 3: The power spectrum $P_{\zeta}$ near the pivot scale $k_{*}=k_{\Delta N_{2}}$ at $\Delta N_{2}=10$ for some values of $\gamma$. $\gamma$ | $0.1$ | $1$ | $10$ | $100$ | $1000$
---|---|---|---|---|---
$\zeta$ | $-0.31165$ | $-1.7774$ | $-8.91495$ | $-42.7976$ | $-201.722$
Table 1: The parameters leading to the power spectrum in Fig. 3 with $\Delta
N_{2}\approx 10$.

The enhancement of primordial curvature perturbations needed for PBH formation
is usually assumed to be ${\textstyle{P_{\zeta,{\rm
max.}}\over\vphantom{2}\smash{\raise
0.60275pt\hbox{$\scriptstyle{P_{\zeta,{\rm min.}}}$}}}}\equiv P_{\rm enh.}\sim
10^{7}$, where the subscripts refer to the values of $P_{\zeta}$ at the peak
and the base of the peak, respectively, in comparison to the CMB scale. Given
a broad peak, the enhancement may be less by one order of the magnitude [42].
According to Fig. 3, when choosing $\Delta N_{2}=10$, the enhancement $P_{\rm
enh.}\gtrsim 10^{6}$ is achieved for $\gamma\gtrsim 10$. Actually, it ranges
from $10^{6}$ at $\gamma=10~{}$ to $10^{7}$ at $\gamma=1000$.

In Table 2 the values of $n_{s}$ and $r_{\rm max.}$ are collected at the CMB
scales for the values of $\Delta N_{2}=10,17,20,23$, universally across the
considered values of $\gamma=0.1,1,10,100,1000$. The tensor-to-scalar ratio
$r$ is well within the observational limits in all those cases, but the scalar
tilt $n_{s}$ is outside the $1\sigma$ limit when $\Delta N_{2}=20$, and is
marginally outside the $3\sigma$ limit when $\Delta N_{2}=23$.

$\Delta N_{2}$ | $10$ | $17$ | $20$ | $23$
---|---|---|---|---
$n_{s}$ | $0.955$ | $0.946$ | $0.942$ | $0.936$
$r_{\rm max}$ | $0.006$ | $0.008$ | $0.009$ | $0.011$
Table 2: The approximate values of $n_{s}$ and $r$ for some choices of $\Delta
N_{2}$.

The mass of PBH created as a result of the primordial power spectrum
enhancement can be estimated from the peak data as follows [43]:

$\displaystyle M_{\rm PBH}\simeq{{M^{2}_{\rm Pl}}\over{H(t_{\rm
peak})}}\exp\left[2(N_{\rm end}-N_{\rm peak})+\int^{t_{60}}_{t_{\rm
peak}}\epsilon(t)H(t)dt\right]~{},$ (22)

where $t_{\rm peak}$ is the time when the wavenumber corresponding to the
power spectrum peak ($k_{\rm peak}$) exits the horizon, and $t_{60}$ is the
time when $k_{60}$ exits the horizon.

The values of $M_{\rm PBH}$ for some values of $\Delta N_{2}$ from Eq. (22)
are shown in Table 3 together with the corresponding values of the spectral
index. Those estimates are universal across the values of
$\gamma=0.1,1,10,100,1000$. On the one hand, PBH with masses smaller than
$\sim 10^{16}{\rm g}$ should have already evaporated until now via Hawking
radiation. Hence, we require $\Delta N_{2}\geq 20$. On the other hand, the
lower $3\sigma$ limit on the spectral index, $n_{s}\approx 0.946$, requires
$\Delta N_{2}<23$. Hence, PBH masses are restricted in our models by ${\cal
O}(10^{16}{\rm g})<M_{\rm PBH}<{\cal O}(10^{19}{\rm g})$ before imposing
current observational constraints on them.

$\Delta N_{2}$ | $10$ | $17$ | $20$ | $23$
---|---|---|---|---
$M_{\rm PBH}$, g | $10^{9}$ | $10^{15}$ | $10^{17}$ | $10^{20}$
$n_{s}$ | $0.955$ | $0.946$ | $0.942$ | $0.936$
Table 3: The PBH masses from Eq. (22) in the $\gamma$ model and the values of
the scalar spectral index $n_{s}$. In the Solar mass units, $1~{}{\rm
g}\approx 5.03\times 10^{-34}~{}M_{\odot}$.

The PBH density fraction in DM can be estimated by using the standard (Press-
Schechter) formalism [44]. The useful formulae include the PBH mass
$\tilde{M}_{\rm PBH}(k)$, the production rate $\beta_{f}(k)$, and the density
contrast $\sigma(k)$ coarse-grained over $k$ as follows (see e.g., Refs. [45,
46] and the references therein):

$\displaystyle\tilde{M}_{\rm PBH}\simeq 10^{20}\left({{7\times
10^{12}}\over{k~{}{\rm Mpc}}}\right)^{2}{\rm
g}~{},\quad\beta_{f}(k)\simeq{{\sigma(k)}\over{\sqrt{2\pi}\delta_{c}}}e^{-{{\delta^{2}_{c}}\over{2\sigma^{2}(k)}}}~{},$
(23)
$\displaystyle\sigma^{2}(k)={{16}\over{81}}\int{{dq}\over{q}}\left({{q}\over{k}}\right)^{4}e^{-q^{2}/k^{2}}P_{\zeta}(q)~{}.$

We have chosen the Gaussian window function for the density contrast, and have
introduced $\delta_{c}$ as a constant representing the density threshold for
PBH formation. It is often assumed that $\delta_{c}\approx 1/3$ [47], though
it may be different. The PBH-to-DM density fraction is estimated as [45, 46]

$\displaystyle{{\Omega_{\rm PBH}(k)}\over{\Omega_{\rm DM}}}\equiv
f(k)\simeq{{1.4\times 10^{24}\beta_{f}(k)}\over{\sqrt{\tilde{M}_{\rm
PBH}(k){\rm g}^{-1}}}}~{}~{}.$ (24)

In order to numerically evaluate the function (24), we have normalized the
values of $k$ in terms of the observable scales today, and the scale $k_{60}$
was chosen to represent the largest currently observable scale around
$10^{-4}~{}{\rm Mpc}^{-1}$. Then our numerical calculation reveals that the
$\delta_{c}$ parameter should be a bit smaller than $1/3$ (actually, close to
$\delta_{c}=0.275$ in the case under consideration).

The $\gamma$-models of inflation and PBH formation studied above have rather
low values of the CMB spectral index (with $n_{s}\approx 0.942$ as the best
fit) which are in tension ($3\sigma$) with the precision measurements of
Planck mission [37] in Eq. (21). The better values are found in Subsection
3.3.

###  3.2 The $\delta$ models

The $\delta$ models, defined by $\gamma=0$ and $\delta\neq 0$ in Eqs. (9) and
(10), break the R-symmetry and the reflection symmetry
$\sigma\rightarrow-\sigma$ of the potential, see Fig. 4. For any non-zero
$\delta$, there is the value of $\zeta$ leading to an inflection point.

Figure 4: The scalar potential for $\gamma=0$, $\delta=0.1$ and
$\zeta=0.033407$.

In contrast to the $\gamma$ models, there is a single valley for large
positive $\varphi$ and $\sigma=0$. When approaching $\varphi=0$, the
inflationary trajectory passes the (near-)inflection point and then falls to
the Minkowski minimum at $\varphi=\sigma=0$.

Let us take the parameter values $\delta=0.1$ and $\zeta=0.033407$, where
$\zeta$ is chosen to get $\Delta N_{2}=10$. A numerical solution to the field
equations yields the time dependence of $\varphi$, $\sigma$, $\tilde{H}$, $N$,
$\epsilon$ and $\eta$ shown in Fig. 5. The near-inflection point divides
inflation into two stages with $\Delta N_{1}=50$ and $\Delta N_{2}=10$,
respectively. The initial velocities are set to zero, with $\varphi(0)=6$ and
$\sigma(0)=0.05$. Similarly to the $\gamma$ models, the inflationary
trajectory is (locally) stable against variations of the initial conditions.

(a)

(b)

(c)

(d)

(e)

Figure 5: (a) The solution to the field equations (17) and (18) with the
initial conditions $\varphi(0)=6$ and $\sigma(0)=0.05$, the vanishing initial
velocities, and the parameter choice $\gamma=0$, $\delta=0.1$ and
$\zeta=0.033407$. (b) The trajectory of the solution ($\varphi$ – blue,
$\sigma$ – red). (c) The Hubble function. (d) The e-folds. (e) The slow-roll
parameters $\epsilon$ (red) and $\eta$ (blue).

The power spectrum enhancement also takes place, either with a smooth peak for
$\delta\sim 0.1$, or with a sharp peak for $\delta\sim 0.6$, in the scalar
power spectrum. When requiring the PBH density fraction (24) to be close to
one and the corresponding density threshold not to deviate far away from the
region $1/3\leq\delta_{c}\leq 2/3$, one finds that the values $\delta=0.094$
and $\delta=0.58$, respectively, are suitable for efficient generation of PBH
with their masses around $10^{19}{\rm g}$, while avoiding their overproduction
($f\lesssim 1$).

The results about the $\delta$-models are summarized in Table 4. To get
$M_{\rm PBH}\sim 10^{19}{\rm g}$ one should either set $\delta=0.094$ and
$\Delta N_{2}=20$ (in this case $P_{\rm enh}\approx 4.5\times 10^{7}$), or
$\delta=0.58$ and $\Delta N_{2}=23$ (in this case $P_{\rm enh}\approx
2.7\times 10^{8}$). In the former case, $n_{s}$ is within $2\sigma$ CL,
whereas in the latter case, $n_{s}$ is within $3\sigma$ but outside $2\sigma$
CL. Therefore, those $\delta$-models of inflation and PBH formation are better
than the $\gamma$ models, as regards matching the predicted values of the CMB
scalar perturbations index $n_{s}$ with Planck measurements [37].

| $\delta=0.09$ | $\delta=0.61$
---|---|---
$\Delta N_{2}$ | $10$ | $17$ | $20$ | $23$ | $10$ | $17$ | $20$ | $23$
$M_{\rm PBH}$, g | $10^{9}$ | $10^{15}$ | $10^{18}$ | $10^{20}$ | $10^{9}$ | $10^{15}$ | $10^{18}$ | $10^{20}$
$n_{s}$ | $0.9566$ | $0.9486$ | $0.9443$ | $0.9390$ | $0.9581$ | $0.9504$ | $0.9461$ | $0.9409$
$r_{\rm max}$ | $0.005$ | $0.007$ | $0.008$ | $0.010$ | $0.004$ | $0.006$ | $0.007$ | $0.008$
Table 4: The PBH masses from Eq. (22) for $\delta=0.09$ and $\delta=0.61$,
with the corresponding values of $n_{s}$ and $r_{\rm max}$.

###  3.3 Comparison with observational constraints on PBH and DM

The specific models studied in the preceding Subsections can be further
improved and compared to the current observational constraints on PBH and DM,
after fine tuning their parameters to the best fit within the modified
Starobinsky-type supergravity [28, 35].

For those purposes we pick up a $\gamma$-model (Case I) and two
$\delta$-models (Cases II and III) with the different shapes of the power
spectrum. The choice of two $\delta$-models is motivated by the existence of
two suitable parameter regions, where $\delta\simeq 0.1$ and $\delta\simeq
0.6$ yield truly different shapes of the power spectrum (broad and narrow
peaks), respectively. The parameters of those three models are given in Table
5, and the corresponding (numerically computed) power spectra $P_{\zeta}$ and
PBH density fractions $f(M)$ are shown in Fig. 6 with the normalization of the
wavenumber $k_{\rm exit}=0.05~{}{\rm Mpc}^{-1}$, where $k_{\rm exit}$ is the
scale that leaves the horizon around $54$ e-folds. The parameter $\zeta$ is
fixed by choosing the value of $\Delta N_{2}$ at given $\gamma$ and $\delta$.
In the cases I, II and III, the parameter $\zeta$ is given by $-2.374$,
$0.032$, and $0.102$, respectively.

| $\gamma$ | $\delta$ | $\Delta N_{2}$ | $\delta_{c}$ | $n_{s}$ | $r$
---|---|---|---|---|---|---
Case I | $1.5$ | $0$ | $20$ | $0.4$ | $0.942$ | $0.009$
Case II | $0$ | $0.09$ | $19$ | $0.47$ | $0.946$ | $0.008$
Case III | $0$ | $0.61$ | $20$ | $0.4$ | $0.946$ | $0.007$
Table 5: The parameters corresponding to the PBH density fraction in Fig. 6.
The $n_{s}$ and $r$ are computed with the total $\Delta N=54$ e-folds before
the end of inflation.

Figure 6: The power spectrum in the three examples of Table 5, see the left
side. The $k_{*}$ represents the end of SR and the beginning of USR. The
corresponding PBH density fractions are given on the right side. The
observational constraints on PBH in the background (on the right side) are
taken from Refs. [10, 11]. In both plots, the case I is denoted by the solid
line, the case II by the dashed line, and the case III by the dotted line. The
whole DM composed from PBH is possible in the cases I and III.

According to Table 5, the spectral tilt $n_{s}$ in the case I is in tension by
$3\sigma$ against the CMB data [37], whereas in the cases II and III the value
of $n_{s}$ is within the current $3\sigma$ constraints. Therefore, the
$\delta$-models are more suitable to accommodate the observed values of
$n_{s}$, though with the PBH fraction less than the whole DM. In particular,
the PBH fraction in the case II peaks at the center of the allowed window,
while it is still possible to move the peak further to the left, thus lowering
the PBH masses. It is also worth noticing that we have used the value of the
critical $\delta_{c}$ parameter beyond $1/3$, see Table 5.

##  4 Scanning the parameter space

In this Section we study the parameter space of the models defined by Eqs. (9)
and (10), with both $\gamma$ and $\delta$ to be non-vanishing. By searching
for the proper values of those parameters, we investigate whether it is
possible to move the value of the cosmological observable $n_{s}$ into the
most favorable region (within $1\sigma$), while keeping the efficient PBH
production. We also study some new models by changing the second (quadratic)
term in our Ansatz (10) to a higher power.

###  4.1 The $(\gamma,\delta)$ models

In the $(\gamma,\delta)$ models, the parameter values should be restricted in
order to get a sufficient length of the second stage of inflation, required
for efficient PBH production (as DM). In addition, with the parameter values
similar to those in the previous Section, the USR phase between the first and
the second plateaus of the inflationary trajectory requires extreme fine-
tuning that should be avoided. As a fully analytic approach is impossible, we
scan the parameter space randomly by doing numerical calculations in our
search for better agreement with the Planck measurements of $n_{s}$.

Figure 7: The power spectrum $P_{\zeta}$ around the pivot scale $k_{*}$ at
$\Delta N_{2}\approx 20$ for the selected values of $\gamma$ and $\delta$.

Figure 7 shows the power spectrum $P_{\zeta}$ around the pivot scale $k_{*}$
for some non-vanishing values of $\gamma$ and $\delta$ given in Table 6, which
were used (separately) in the preceding Section.

In Table 6 we also collect the values of the parameters $\gamma$, $\delta$ and
$\zeta$, where $\zeta$ has been tuned to get $\Delta N_{2}\approx 20$,
together with the values of the power spectrum enhancement $P_{\rm enh.}$
(relative to CMB) for them, as our first (naive) trials.

Like the special cases of the preceding Section, going well below $20$ for
$\Delta N_{2}$ leads to an increase of $n_{s}$ accompanied by a drastic
decrease in the PBH masses (significantly below the Hawking radiation limit)
and a much lower enhancement of the scalar power spectrum, see Fig. 8. We
assume $\Delta N_{2}\approx 20$ in what follows.

As can be seen from Table 6, the adjustment of the value of $\zeta$ to get
$\Delta N_{2}\approx 20$ is generically possible in a large region in the
parameter space $(\gamma,\delta)$, while the required value of $\zeta$ mainly
depends on the value of $\gamma$.

$\gamma$ | $0.1$ | $0.1$ | $0.5$ | $0.5$ | $1.0$
---|---|---|---|---|---
$\delta$ | $0.1$ | $0.5$ | $0.1$ | $0.5$ | $0.5$
$\zeta$ | $-0.3335$ | $-0.3934$ | $-1.1137$ | $-1.2280$ | $-1.9826$
$\Delta N_{2}$ | $20.23$ | $20.11$ | $20.35$ | $20.39$ | $20.06$
$P_{\rm enh.}$ | $4.4\times 10^{4}$ | $1.2\times 10^{3}$ | $1.7\times 10^{5}$ | $5.8\times 10^{3}$ | $1.2\times 10^{4}$
Table 6: The selected values of the parameters $\gamma$, $\delta$ and $\zeta$
leading to $\Delta N_{2}\approx 20$, and the power spectrum enhancement
$P_{\rm enh.}$ (relative to CMB), respectively.

As is clear from Table 6 and Fig. 7, we observe the sharp decrease in the
scalar power spectrum enhancement, well below $10^{6}$. Therefore, stricter
restrictions on the values of the parameters $\gamma$ and $\delta$ are
necessary, while their existence for efficient PBH production (as DM) seems to
be non-trivial. It signals a qualitative departure from the two special cases
considered in the preceding Section.

Figure 8: The power spectrum $P_{\zeta}$ near the pivot scale $k_{*}=k_{\Delta
N_{2}}$ for some values of $\Delta N_{2}$ less or equal $20$.

Afer studying the power spectra around $\gamma,\delta\approx 1.0$ with $\Delta
N_{2}\approx 20$, we find that their enhancement (peak) becomes larger when
$\gamma$ increases and $\delta$ decreases, though only slightly. The
sufficient enhancement of primordial curvature perturbations, needed for the
efficient PBH formation, is given by $P_{enh.}\sim 10^{6}$ at least [42]. It
follows from our numerical calculation and animation that the value of
$\gamma$ in the models with $\delta\geq 0.1$ should be as large as $10$ at
least. Accordingly, we are led to the representative values of $\gamma$ and
$\delta$ as $\gamma=10$, $\delta=0.2$ and adjust $\zeta=-9.15915$. The shape
of the (canonical) scalar potential in that $(\gamma,\delta)$ model is shown
in Fig. 9(a). As may have been expected, it combines the shapes of the
potentials observed in the $\gamma$\- and $\delta$-models, see Figs. 1 and 4,
respectively. In particular, the potential in Fig. 9(a) has the same pattern
of the reflection symmetry breaking seen in the $\delta$ models.

The evolution of the scalars is given by a solution to the equations of motion
in the FLRW spacetime, and it is displayed in Fig. 9. The evolution of
Starobinsky’s scalaron $\varphi$ and the SR parameter $\epsilon$ are shown by
the red lines, whereas the evolution of the field $\sigma$, the SR parameter
$\eta$ and the Hubble function $H$ are shown by the blue lines. We use the
initial conditions $\varphi(0)=6$ and $\sigma(0)=0.1$, with the vanishing
velocities. Should the initial condition for the $\sigma$-field be with the
opposite sign, say, $\delta(0)=-0.1$, there will be no second plateau, no
second stage of inflation and, hence, no PBH production.

(a)

(b)

(c)

(d)

(e)

Figure 9:  (a) The scalar potential of the model, and the two solutions
corresponding to the initial conditions $\sigma(0)=0.1$ (blue) and
$\sigma(0)=-0.1$ (red), respectively. The red solution does not have double
inflation and should be discarded. (b) The power spectrum around the pivot
scale. (c) The solution to the field equations (17) and (18) for $\varphi$
(red) and $\sigma$ (blue) with the initial conditions
$\varphi(0)=6,\sigma(0)=0.1$, the vanishing initial velocities, and the
parameters $\delta=0.2$, $\gamma=10$ and $\zeta=-9.15915$. (d) The Hubble
function. (e) The SR parameters $\epsilon$ (red) and $\eta$ (blue). The blue
shaded region represents the first stage of inflation, and the green shaded
region represents the second stage.

Our results for the CMB observables $n_{s}$ and $r$ at some values of
$\gamma\geq 10$ with $\delta=0.2$, all leading to the sufficient enhancement
of the power spectrum for the efficient PBH production, are collected in Table
7 together with the masses of the induced PBH.

$\gamma$ | $\zeta$ | $\Delta N_{2}$ | $n_{s}$ | $r$ | $M_{\mathrm{P}BH}$ | $P_{enh.}$
---|---|---|---|---|---|---
$10$ | $-9.15915$ | $20.00$ | $0.9407$ | $0.0097$ | $3.5\times 10^{16}$ | $3.6\times 10^{6}$
$15$ | $-12.06468$ | $19.94$ | $0.9409$ | $0.0096$ | $3.8\times 10^{16}$ | $5.0\times 10^{6}$
$20$ | $-14.66315$ | $20.02$ | $0.9407$ | $0.0097$ | $4.9\times 10^{16}$ | $5.1\times 10^{6}$
Table 7: The values of $n_{s}$ and $r$ together with the PBH masses
$M_{\mathrm{PBH}}$ and the corresponding enhancement $P_{\rm enh.}$ of the
scalar power spectrum for some values of $\gamma$ above $10$, with $\Delta
N_{2}\approx 20$.

It follows from Table 7 that the scalar tilt $n_{s}$ is still in tension with
the Planck measurements (by $3\sigma$ or more), despite the presence of an
extra (adjustable) parameter. Therefore, the naive generalization of the
special models by using two non-vanishing parameters $\gamma$ and $\delta$
within the Ansatz defined by Eqs. (9) and (10) does not lead to a perfect
match (within $1\sigma$) with Planck observations, contrary to the
expectations raised in Refs. [28, 35]. It implies that the ”pure”
$\delta$-models of inflation and PBH formation, defined in the preceding
Section, are special when compared to their hybrid extensions studied in this
Subsection, as regards PBH production (as DM) and observational constraints.
Those $\delta$-models are distinguished by a chiral deformation of the
modified supergravity potentials, so that we focus on generalizations of the
${\cal F}$ term in the next Subsection.

###  4.2 The new models with higher powers of ${\cal R}$ in the ${\cal F}$
term

Let us now revise our Ansatz (7) by changing the chiral potential (10),
because a more general non-chiral potential (9) apparently does not lead to a
significant improvement in the value of $n_{s}$, according to our studies in
the preceding Subsection.

The linear term in the function ${\cal F}$ is required to generate the
Einstein-Hilbert term for gravity, so that a generic chiral potential reads

$\displaystyle{\cal F}({\cal R})=-3{\cal R}+f({\cal R})~{},$ (25)

where the (complex) function $f$ of the chiral superfield ${\cal R}$ should be
quadratic at least (we assume it in the polynomial form).

The simplest new models are defined by using the function $f$ as a higher
power $(n>2)$ of ${\cal R}$,

$\displaystyle f({\cal R})\equiv{{C_{n}}\over{M^{n-1}}}{\cal
R}^{n}\quad\mathrm{and}\quad X={{M}\over{\sqrt{24}}}\sigma~{},$ (26)

where we have ignored the angular mode of $X=\left.{\cal R}\right|$ again.

The effective action for two scalars takes the same form (15), though with new
functions $A,B$ and $U$ of $\sigma$. After a straightforward calculation, we
find

$\displaystyle A$
$\displaystyle=1+{{1}\over{6}}\sigma^{2}-{{11}\over{24}}\zeta\sigma^{4}-{{29}\over{54}}\gamma\sigma^{6}-{{1}\over{6}}(f_{X}+\mkern
1.5mu\overline{\mkern-1.5muf\mkern-1.5mu}\mkern 1.5mu_{\mkern
1.5mu\overline{\mkern-1.5muX\mkern-1.5mu}\mkern 1.5mu})~{},$ $\displaystyle B$
$\displaystyle={{1}\over{3M^{2}}}(1-\zeta\sigma^{2}-\gamma\sigma^{4})~{},$
(27) $\displaystyle U$
$\displaystyle={{M^{2}}\over{2}}\sigma^{2}\left[1-{{1}\over{6}}\sigma^{2}+{{3}\over{8}}\zeta\sigma^{4}+{{25}\over{54}}\gamma\sigma^{6}+{{1}\over{3}}(f_{X}+\mkern
1.5mu\overline{\mkern-1.5muf\mkern-1.5mu}\mkern 1.5mu_{\mkern
1.5mu\overline{\mkern-1.5muX\mkern-1.5mu}\mkern
1.5mu})-{{1}\over{2}}\left({{f}\over{X}}+{{\mkern
1.5mu\overline{\mkern-1.5muf\mkern-1.5mu}\mkern 1.5mu}\over{\mkern
1.5mu\overline{\mkern-1.5muX\mkern-1.5mu}\mkern 1.5mu}}\right)\right]~{}.$

The $(\gamma,\delta)$ models are reproduced by choosing
$f(X)={\textstyle{3\sqrt{6}\over\vphantom{2}\smash{\raise
0.60275pt\hbox{$\scriptstyle{M}$}}}}\delta X$.

####  4.2.1 The model with a cubic chiral term

The model is defined by the cubic chiral function $f$ with the new parameter
$\lambda$ as follows:

$\displaystyle{\cal F}=-3{\cal R}+{{24}\over{M^{2}}}\lambda{\cal
R}^{3}~{}~{}.$ (28)

The corresponding functions $A,B$ and $U$ are given by

$\displaystyle
A=1+\left({{1}\over{6}}-\lambda\right)\sigma^{2}-{{11}\over{24}}\zeta\sigma^{4}-{{29}\over{54}}\gamma\sigma^{6},$
$\displaystyle
B={{1}\over{3M^{2}}}\left(1-\zeta\sigma^{2}-\gamma\sigma^{4}\right),$ (29)
$\displaystyle
U={{1}\over{2}}M^{2}\sigma^{2}\left[1+\left(\lambda-{{1}\over{6}}\right)\sigma^{2}+{{3}\over{8}}\zeta\sigma^{4}+{{25}\over{54}}\gamma\sigma^{6}\right]~{}.$

The proper values of the parameters $\gamma,\lambda$ and $\zeta$, which are
needed for efficient PBH production, are collected in Table 8, where the
parameter $\zeta$ has been tuned to get $\Delta N_{2}\approx 20$.

$\gamma$ | $1.0$ | $1.0$ | $10$ | $10$ | $20$
---|---|---|---|---|---
$\lambda$ | $0.1$ | $0.2$ | $0.1$ | $0.2$ | $0.2$
$\zeta$ | $-1.92564$ | $-2.06496$ | $-9.20329$ | $-9.48826$ | $-15.044427$
$\Delta N_{2}$ | $19.97$ | $19.99$ | $20.09$ | $19.95$ | $20.01$
$n_{s}$ | $0.9407$ | $0.9402$ | $0.9407$ | $0.9408$ | $0.9408$
$r$ | $0.0097$ | $0.0099$ | $0.0097$ | $0.0096$ | $0.0096$
${\cal P}_{\mathrm{e}nh.}$ | $3.0\times 10^{5}$ | $2.1\times 10^{6}$ | $3.1\times 10^{6}$ | $4.0\times 10^{6}$ | $3.8\times 10^{6}$
${M}_{\mathrm{PBH}}({\rm g})$ | $1.6\times 10^{16}$ | $1.5\times 10^{11}$ | $4.2\times 10^{16}$ | $3.7\times 10^{16}$ | $5.6\times 10^{16}$
Table 8: The parameters $\gamma$, $\lambda$ and $\zeta$ in the cubic model,
and the resulting PBH values.

As can be seen from Table 8, $\gamma\geq 10$ for $\lambda\sim 0.1$ is required
to get a sufficient enhancement of the power spectrum. The value of $n_{s}$ is
sensitive to large changes in the value of $\Delta N_{2}$ that is one of the
major factors contributing to the mass $M_{\rm{PBH}}$, as long as the value of
$\lambda$ is much smaller than $\gamma$. The obtained values of $n_{s}$ are
below $0.946$, i.e. outside the observational constrains by $3\sigma$ or more.

(a)

(b)

(c)

(d)

(e)

Figure 10:  (a) The shape of the potential and the inflation trajectory in the
cubic model with the parameters and the initial conditions given under (c).
(b) The power spectrum around the pivot scale. (c) The solution to the field
equations (17) and (18) in the cubic model (28) with the initial conditions
$\varphi(0)=6,\sigma(0)=0.1$, the vanishing initial velocities, and the
parameters $\lambda=0.2$, $\gamma=10$ and $\zeta=-9.48826$. (d) The Hubble
function. (e) The SR parameters $\epsilon$ (red) and $\eta$ (blue) in the
cubic model. The blue shaded region represents the first stage of inflation,
and the green shaded region represents the second stage.

As a representative of those models, let us take the one with $\gamma=10$,
$\lambda=0.2$ and $\zeta=-9.48826$. These values are determined by fixing the
value of the power spectrum enhancement as $10^{6}$, after setting the length
of the second stage of inflation by about 20 e-folds. Fig. 10 shows the shape
of the potential and the inflationary trajectory (a), the power spectrum (b),
the evolution of both scalar fields (c), the shape of the Hubble function (d),
and the slow roll parameters (e), in the cubic model with the parameters
chosen above.

Unlike the quadratic model, the potential becomes completely symmetric with
respect to a change of the sign of $\sigma$. As a result of that, the shape of
the trajectory is not affected by the sign of the initial condition on
$\sigma$. It implies that the inflationary trajectory is more stable against
changes in the initial value of the inflaton field.

Using the higher powers $n$ of ${\cal R}$ until $n=7$ in the chiral function
${\cal F}$ leads to the similar results. For instance, the PBH production in
the quartic model with $n=4$ appears to be similar to that in the quadratic
$\delta$-models, with the potential that is narrower along the $\sigma$-axis
and resembles the potentials in the $(\gamma,\delta)$ models with a higher
value of $\gamma$. It may have been expected because all those models merely
change some coefficients in the functions $A,B$ and $U$. In particular, the
value of $n_{s}$ is not significantly affected.

####  4.2.2 The model with $n=9$

One may wonder, what happens when the power of ${\cal R}$ in the holomorphic
potential ${\cal F}$ is increased further in Eq. (28) ? As is clear from Eq.
(4.2.1), until $n=9$ it would merely affect some of the coefficients in the
polynomials there, without producing new terms. Therefore, we should not
expect significant changes from the simplest $\delta$-models until $n=8$.
Going beyond $n=8$ leads to the appearance of the higher powers of $\sigma$
according to Eq. (27), when keeping the symmetry under the sign change of
$\sigma$.

The superpotential ${\cal F}$ for the $n=9$ model with the new parameter $\mu$
reads as follows:

$\displaystyle{\cal F}=-3{\cal R}+{{(24)^{4}}\over{M^{8}}}\mu{\cal
R}^{9}~{}~{}.$ (30)

(a)

(b)

(c)

(d)

(e)

Figure 11:  (a) The scalar potential of the model, and the inflationary
trajectory with the initial conditions given under (c). (b) The power spectrum
around the pivot scale. (c) The solution to the field equations (17) and (18)
in the $n=9$ model (30) for $\varphi$ (red) and $\sigma$ (blue) with the
initial conditions $\varphi(0)=6,\sigma(0)=0.1$, the vanishing initial
velocities, and the parameters $\gamma=1.0$, $\delta=0.1$ and
$\zeta=-2.41238$. (d) The Hubble function. (e) The SR parameters, $\epsilon$
(red) and $\eta$ (blue). The blue shaded region represents the first stage of
inflation, and the green shaded region represents the second stage.

As a representative of these models, let us take $\gamma=1.0$, $\lambda=0.1$
and $\zeta=-2.41238$, towards the desired value of $\Delta N_{2}\approx 20$
with at least the order $10^{6}$ enhancement on the scalar power spectrum. The
corresponding plots are given in Fig. 11.

After following the procedure described in the preceding Subsection, we find
the values collected in Table 9.

$\gamma$ | $1.0$ | $1.0$ | $1.0$ | $2.0$ | $5.0$
---|---|---|---|---|---
$\mu$ | $0.01$ | $0.1$ | $0.2$ | $0.1$ | $0.1$
$\zeta$ | $-1.858261$ | $-2.41238$ | $-2.87306$ | $-3.37249$ | $-5.8009472$
$\Delta N_{2}$ | $20.00$ | $20.07$ | $20.06$ | $19.99$ | $20.05$
$n_{s}$ | $0.9408$ | $0.9402$ | $0.9314$ | $0.9405$ | $0.9407$
$r$ | $0.0096$ | $0.0093$ | $0.0059$ | $0.0098$ | $0.0097$
${\cal P}_{\mathrm{e}nh.}$ | $4.3\times 10^{6}$ | $1.9\times 10^{6}$ | $3.7\times 10^{6}$ | $4.1\times 10^{6}$ | $2.7\times 10^{6}$
${M}_{\mathrm{PBH}}({\rm g})$ | $3.3\times 10^{16}$ | $6.6\times 10^{14}$ | $7.0\times 10^{11}$ | $3.2\times 10^{16}$ | $5.6\times 10^{16}$
Table 9: The parameters $\gamma$, $\mu$ and $\zeta$ in the $n=9$ model (30),
and the resulting PBH values.

As is clear from Table 9, the enhancement of the scalar power spectrum in the
$n=9$ model (30) can be as large as $10^{6}$ with the relatively smaller
values of $\gamma$ in comparison to the previous models, but the predicted
values of the tilt $n_{s}$ are once again outside the $3\sigma$ region of the
Planck measurements. The calculated value of $n_{s}$ increases when the value
of $\mu$ decreases, which leads us back to the $\gamma$-models. It appears to
be difficult to increase the masses of the generated PBH beyond $10^{16}$ g in
the $n=9$ model. The parameter $\zeta$ needs to be fine-tuned as before.

Therefore, it seems to be impossible to improve the values of cosmological
observables (mainly $n_{s}$) by changing the second term in Eq. (10) to
another term with a higher power.

##  5 Conclusion

Our extensions of the Starobinsky inflation model in the modified supergravity
are capable to produce PBH that may account for part or the whole Cold Dark
Matter, with the PBH masses in the range between $10^{16}$ g and $10^{20}$ g.
For example, the NANOGrav Collaboration data [48] hints to the PBH as DM [49],
in agreement with our results in Fig. 6.

The PBH generation can be efficiently catalyzed by primordial perturbations
sourced by Starobinsky’s scalaron coupled to another scalar of the
(super)gravitational origin, leading to the effective two-field inflation. Our
conclusions are based on theoretical considerations with the two fundamental
principles: modified gravity and supersymmetry, starting from the manifestly
(locally) supersymmetric Lagrangians, in the top-down approach.

Our models have a limited number of free parameters that have to be properly
tuned. Surprisingly enough, a naive increase in the number of the parameters
does not relax their tuning needed for viable inflation and efficient PBH
production (as DM). We found that it is highly non-trivial (and, perhaps,
impossible) to meet the $1\sigma$ Planck value of $n_{s}$ in our models for a
large part of the parameter space, though it is still possible to get an
agreement between $2\sigma$ and $3\sigma$ in the $\delta$-models with
$\gamma=\lambda=\mu=0$,

It is worth mentioning that we insisted on the duration of the second stage of
inflation to be close to 20 e-folds that allows the whole DM as PBH. Should
the lower $n_{s}$ values (close to $0.946$) be excluded by $5\sigma$, only a
small part (under 10%) of DM as PBH can be realized in our models.

Isocurvature perturbations play the key role for PBH formation in our approach
because their growth during USR is part of our amplification mechanism. Large
isocurvature perturbations appear when the effective isocurvature mass becomes
tachyonic and the inflationary trajectory makes a sharp turn. It is the
genuine multi-field inflation effect that cannot be realized in single-field
models of inflation, see e.g., Refs. [29, 31] for details. Our mechanism of
the amplification of curvature perturbations is different from the one
employed in Ref. [50] where the kinetic coupling between two scalars is
responsible for the tachyonic isocurvature mass.

As is well known, PBH formation necessarily leads to gravitational waves (GW)
because large scalar over-densities act as a source of the stochastic GW
background. Frequencies of those GW can be related to the expected PBH masses
and the duration of the second stage of inflation [51, 52]. Those GW may be
detected in the future ground-based experiments, such as the Einstein
telescope [53] and the global network of GW interferometers including the
advanced LIGO, Virgo and KAGRA [54], as well as the space-based GW
interferometers such as LISA [55], TAIJI (old ALIA) [56, 57], TianQin [58] and
DECIGO [59].

Our models of inflation and PBH production may be tested by observations via
detection of the stochastic GW induced by the PBH formation. Using the power
spectra on the left side of Fig. 6, the density $\Omega_{\rm GW}(k)$ in terms
of frequency $k=2\pi f$ was computed in Ref. [35], with the results collected
in Fig. 12 including the expected sensitivity curves for several space-based
GW experiments planned in the near future. To draw the sensitivity curves, the
parameters and the noise models for LISA [55], TianQin [58], Taiji [56], and
DECIGO [59] have been used. Our Fig. 12 here is different from the one given
in Ref. [35] because the updated sensitivity curve of the Taiji project [57]
has been used.

It follows from Fig. 12 that the upcoming space-based GW experiments may be
sensitive enough to detect the stochastic GW background predicted by our
models. According to Fig. 12, our supergravity models produce GW peaks in the
frequency range $10^{-3}\div 10^{-1}$ Hz expected to be accessible by the
LISA, TianQin, Taiji, and DECIGO gravitational interferometers.

It is yet to be seen to what extent future observations of CMB and GW would
allow us to reconstruct our models and distinguish them from others. One
possibility is to look at spectral distortions as a probe of primordial
density perturbations, as is proposed in Ref. [60]. Dissipation of those
perturbations through photon diffusion will distort the CMB spectrum at an
observable level [61, 62]. Another possibility may be looking at specific
oscillatory features of the stochastic GW background [63].

Figure 12:  The density of stochastic GW induced by the power spectrum
enhancement in the supergravity models studied in Subsection 3.3: the case I
(solid black curve), the case II (dashed black curve), and the case III
(dotted black curve). The expected sensitivity curves for the space-based GW
experiments are represented by the different colors.

Fine tuning in our best models amounts to fixing the parameter $M\sim
10^{-5}M_{\rm Pl}$ as the scalaron mass, the dimensionless parameter $\zeta$
for the desired duration of the USR phase between two stages of inflation, and
the parameter $\delta$ for agreement with the Planck measurements of $n_{s}$.
The obtained PBH masses and cosmological tilts agree with all astrophysical
and cosmological constrains, including $n_{s}$ (in the $\delta$-models). We
are not aware of any fundamental reason for the $\delta$-models to be
selected, beyond the phenomenological constraints.

Supergravity is usually regarded as a fundamental extension of gravity at
super-high energy scales. Our new findings demonstrate that the new scalars of
modified supergravity can play the important role during inflation, catalyze
PBH formation and produce GW radiation. The interactions of those scalars are
dictated by local SUSY and are not assumed ad hoc. Our models thus have the
predictive power that can be confirmed or falsified in future experiments,
also implying that indirect footprints of SUSY may be detectable from GW
physics rather than high-energy particle colliders!

## Acknowledgements

RI and SVK were supported by Tokyo Metropolitan University. SVK was also
supported by the World Premier International Research Center Initiative, MEXT,
Japan, and the Development Program of Tomsk Polytechnic University within the
assignment of the Ministry of Science and Higher Education of the Russian
Federation. The authors are grateful to M. Braglia, J. Chluba, H. Murayama, S.
Pi and A.A. Starobinsky for discussions and correspondence.

## References

  * [1] J. Ellis, M. A. G. Garcia, N. Nagata, D. V. Nanopoulos, K. A. Olive, and S. Verner, “Building models of inflation in no-scale supergravity,” Int. J. Mod. Phys. D 29 no. 16, (2020) 2030011, arXiv:2009.01709 [hep-ph].
  * [2] S. V. Ketov and M. Y. Khlopov, “Cosmological Probes of Supersymmetric Field Theory Models at Superhigh Energy Scales,” Symmetry 11 no. 4, (2019) 511.
  * [3] I. Novikov and Y. Zeldovic, “Cosmology,” Ann. Rev. Astron. Astrophys. 5 (1967) 627–649.
  * [4] S. Hawking, “Gravitationally collapsed objects of very low mass,” Mon. Not. Roy. Astron. Soc. 152 (1971) 75.
  * [5] B. J. Carr and S. W. Hawking, “Black holes in the early Universe,” Mon. Not. Roy. Astron. Soc. 168 (1974) 399–415.
  * [6] A. Dolgov and J. Silk, “Baryon isocurvature fluctuations at small scales and baryonic dark matter,” Phys. Rev. D 47 (1993) 4244–4255.
  * [7] J. D. Barrow, E. J. Copeland, and A. R. Liddle, “The Cosmology of black hole relics,” Phys. Rev. D 46 (1992) 645–657.
  * [8] B. J. Carr, “Primordial black holes as a probe of cosmology and high energy physics,” Lect. Notes Phys. 631 (2003) 301–321, arXiv:astro-ph/0310838.
  * [9] M. Sasaki, T. Suyama, T. Tanaka, and S. Yokoyama, “Primordial black holes—perspectives in gravitational wave astronomy,” Class. Quant. Grav. 35 no. 6, (2018) 063001, arXiv:1801.05235 [astro-ph.CO].
  * [10] B. Carr, K. Kohri, Y. Sendouda, and J. Yokoyama, “Constraints on Primordial Black Holes,” arXiv:2002.12778 [astro-ph.CO].
  * [11] B. Carr and F. Kuhnel, “Primordial Black Holes as Dark Matter: Recent Developments,” arXiv:2006.02838 [astro-ph.CO].
  * [12] G. Dvali, “$S$-Matrix and Anomaly of de Sitter,” Symmetry 13 no. 1, (2020) 3, arXiv:2012.02133 [hep-th].
  * [13] H. Ooguri, E. Palti, G. Shiu, and C. Vafa, “Distance and de Sitter Conjectures on the Swampland,” Phys. Lett. B 788 (2019) 180–184, arXiv:1810.05506 [hep-th].
  * [14] W. L. Freedman, “Cosmology at a Crossroads,” Nature Astron. 1 (2017) 0121, arXiv:1706.02739 [astro-ph.CO].
  * [15] W. Hu and I. Sawicki, “Models of f(R) Cosmic Acceleration that Evade Solar-System Tests,” Phys. Rev. D 76 (2007) 064004, arXiv:0705.1158 [astro-ph].
  * [16] S. A. Appleby and R. A. Battye, “Do consistent $F(R)$ models mimic General Relativity plus $\Lambda$?,” Phys. Lett. B 654 (2007) 7–12, arXiv:0705.3199 [astro-ph].
  * [17] A. A. Starobinsky, “Disappearing cosmological constant in f(R) gravity,” JETP Lett. 86 (2007) 157–163, arXiv:0706.2041 [astro-ph].
  * [18] A. De Felice and S. Tsujikawa, “f(R) theories,” Living Rev. Rel. 13 (2010) 3, arXiv:1002.4928 [gr-qc].
  * [19] A. A. Starobinsky, “A new type of isotropic cosmological models without singularity,” Phys. Lett. B 91 no. 1, (1980) 99 – 102.
  * [20] S. V. Ketov, “On the equivalence of Starobinsky and Higgs inflationary models in gravity and supergravity,” J. Phys. A 53 no. 8, (2020) 084001, arXiv:1911.01008 [hep-th].
  * [21] S. V. Ketov and A. A. Starobinsky, “Inflation and non-minimal scalar-curvature coupling in gravity and supergravity,” JCAP 08 (2012) 022, arXiv:1203.0805 [hep-th].
  * [22] Y. Aldabergenov and S. V. Ketov, “Removing instability of inflation in Polonyi–Starobinsky supergravity by adding FI term,” Mod. Phys. Lett. A 91 no. 05, (2018) 1850032, arXiv:1711.06789 [hep-th].
  * [23] I. Antoniadis, “Inflation from supersymmetry breaking,” Int. J. Mod. Phys. A 33 no. 31, (2018) 1844021.
  * [24] S. V. Ketov, “Inflation, Dark Energy and Dark Matter in Supergravity,” in Meeting of the Division of Particles and Fields of the American Physical Society.  9, 2019.  arXiv:1909.05599 [hep-th].
  * [25] S. V. Ketov and A. A. Starobinsky, “Embedding $(R+R^{2})$-Inflation into Supergravity,” Phys. Rev. D 83 (2011) 063512, arXiv:1011.0240 [hep-th].
  * [26] S. V. Ketov and T. Terada, “Old-minimal supergravity models of inflation,” JHEP 12 (2013) 040, arXiv:1309.7494 [hep-th].
  * [27] A. Addazi and S. V. Ketov, “Energy conditions in Starobinsky supergravity,” JCAP 03 (2017) 061, arXiv:1701.02450 [hep-th].
  * [28] Y. Aldabergenov, A. Addazi, and S. V. Ketov, “Primordial black holes from modified supergravity,” Eur. Phys. J. C 80 no. 10, (2020) 917, arXiv:2006.16641 [hep-th].
  * [29] S. V. Ketov, “Multi-Field versus Single-Field in the Supergravity Models of Inflation and Primordial Black Holes,” Universe 7 no. 5, (2021) 115.
  * [30] A. Gundhi and C. F. Steinwachs, “Scalaron-Higgs inflation,” Nucl. Phys. B954 (2020) 114989, arXiv:1810.10546 [hep-th].
  * [31] A. Gundhi, S. V. Ketov, and C. F. Steinwachs, “Primordial black hole dark matter in dilaton-extended two-field Starobinsky inflation,” Phys. Rev. D 103 no. 8, (2021) 083518, arXiv:2011.05999 [hep-th].
  * [32] J. Wess and J. Bagger, Supersymmetry and supergravity.  Princeton University Press, Princeton, NJ, USA, 1992\.
  * [33] S. Cecotti, “Higher derivative supergravity is equivalent to standard supergravity coupled to matter. 1.,” Phys. Lett. B190 (1987) 86–92.
  * [34] J. Gates, S.James and S. V. Ketov, “Superstring-inspired supergravity as the universal source of inflation and quintessence,” Phys. Lett. B 674 (2009) 59–63, arXiv:0901.2467 [hep-th].
  * [35] Y. Aldabergenov, A. Addazi, and S. V. Ketov, “Testing Primordial Black Holes as Dark Matter in Supergravity from Gravitational Waves,” Phys. Lett. B 814 (2021) 136069, arXiv:2008.10476 [hep-th].
  * [36] J. Ellis, D. V. Nanopoulos, and K. A. Olive, “A no-scale supergravity framework for sub-Planckian physics,” Phys. Rev. D 89 no. 4, (2014) 043502, arXiv:1310.4770 [hep-ph].
  * [37] Planck Collaboration, Y. Akrami et al., “Planck 2018 results. X. Constraints on inflation,” Astron. Astrophys. 641 (2020) A10, arXiv:1807.06211 [astro-ph.CO].
  * [38] H. Motohashi, A. A. Starobinsky, and J. Yokoyama, “Inflation with a constant rate of roll,” JCAP 09 (2015) 018, arXiv:1411.5021 [astro-ph.CO].
  * [39] D. J. Mulryne, D. Seery, and D. Wesley, “Moment transport equations for non-Gaussianity,” JCAP 01 (2010) 024, arXiv:0909.2256 [astro-ph.CO].
  * [40] D. J. Mulryne, D. Seery, and D. Wesley, “Moment transport equations for the primordial curvature perturbation,” JCAP 04 (2011) 030, arXiv:1008.3159 [astro-ph.CO].
  * [41] M. Dias, J. Frazer, and D. Seery, “Computing observables in curved multifield models of inflation—A guide (with code) to the transport method,” JCAP 12 (2015) 030, arXiv:1502.03125 [astro-ph.CO].
  * [42] C. Germani and I. Musco, “Abundance of Primordial Black Holes Depends on the Shape of the Inflationary Power Spectrum,” Phys. Rev. Lett. 122 no. 14, (2019) 141302, arXiv:1805.04087 [astro-ph.CO].
  * [43] S. Pi, Y.-l. Zhang, Q.-G. Huang, and M. Sasaki, “Scalaron from $R^{2}$-gravity as a heavy field,” JCAP 05 (2018) 042, arXiv:1712.09896 [astro-ph.CO].
  * [44] W. H. Press and P. Schechter, “Formation of galaxies and clusters of galaxies by selfsimilar gravitational condensation,” Astrophys. J. 187 (1974) 425–438.
  * [45] K. Inomata, M. Kawasaki, K. Mukaida, Y. Tada, and T. T. Yanagida, “Inflationary Primordial Black Holes as All Dark Matter,” Phys. Rev. D 96 no. 4, (2017) 043504, arXiv:1701.02544 [astro-ph.CO].
  * [46] K. Inomata, M. Kawasaki, K. Mukaida, and T. T. Yanagida, “Double inflation as a single origin of primordial black holes for all dark matter and LIGO observations,” Phys. Rev. D 97 no. 4, (2018) 043514, arXiv:1711.06129 [astro-ph.CO].
  * [47] B. J. Carr, “The Primordial black hole mass spectrum,” Astrophys. J. 201 (1975) 1–19.
  * [48] NANOGrav Collaboration, Z. Arzoumanian et al., “The NANOGrav 12.5 yr Data Set: Search for an Isotropic Stochastic Gravitational-wave Background,” Astrophys. J. Lett. 905 no. 2, (2020) L34, arXiv:2009.04496 [astro-ph.HE].
  * [49] V. De Luca, G. Franciolini, and A. Riotto, “NANOGrav Data Hints at Primordial Black Holes as Dark Matter,” Phys. Rev. Lett. 126 no. 4, (2021) 041303, arXiv:2009.08268 [astro-ph.CO].
  * [50] M. Braglia, D. K. Hazra, F. Finelli, G. F. Smoot, and A. A. Starobinsky, “Generating PBHs and small-scale GWs in two-field models of inflation,” arXiv:2005.02895 [astro-ph.CO].
  * [51] R.-G. Cai, S. Pi, and M. Sasaki, “Gravitational Waves Induced by non-Gaussian Scalar Perturbations,” Phys. Rev. Lett. 122 no. 20, (2019) 201101, arXiv:1810.11000 [astro-ph.CO].
  * [52] N. Bartolo, V. De Luca, G. Franciolini, A. Lewis, M. Peloso, and A. Riotto, “Primordial Black Hole Dark Matter: LISA Serendipity,” Phys. Rev. Lett. 122 no. 21, (2019) 211301, arXiv:1810.12218 [astro-ph.CO].
  * [53] ET Collaboration, M. Punturo et al., “The Einstein Telescope: a third-generation gravitational wave observatory,” Class. Quant. Grav. 27 no. 19, (9, 2010) 194002.
  * [54] KAGRA, LIGO Scientific, VIRGO Collaboration, B. Abbott et al., “Prospects for Observing and Localizing Gravitational-Wave Transients with Advanced LIGO, Advanced Virgo and KAGRA,” Living Rev. Rel. 21 no. 1, (2018) 3, arXiv:1304.0670 [gr-qc].
  * [55] LISA Collaboration, P. Amaro-Seoane et al., “Laser Interferometer Space Antenna,” arXiv:1702.00786 [astro-ph.IM].
  * [56] X. Gong et al., “Descope of the ALIA mission,” J. Phys. Conf. Ser. 610 no. 1, (2015) 012011, arXiv:1410.7296 [gr-qc].
  * [57] W.-H. Ruan, Z.-K. Guo, R.-G. Cai, and Y.-Z. Zhang, “Taiji program: Gravitational-wave sources,” International Journal of Modern Physics A 35 no. 17, (Jun, 2020) 2050075.
  * [58] TianQin Collaboration, J. Luo et al., “TianQin: a space-borne gravitational wave detector,” Class. Quant. Grav. 33 no. 3, (2016) 035010, arXiv:1512.02076 [astro-ph.IM].
  * [59] H. Kudoh, A. Taruya, T. Hiramatsu, and Y. Himemoto, “Detecting a gravitational-wave background with next-generation space interferometers,” Phys. Rev. D 73 (2006) 064006, arXiv:gr-qc/0511145.
  * [60] J. Chluba, M. H. Abitbol, N. Aghanim, Y. Ali-Haimoud, M. Alvarez, K. Basu, B. Bolliet, C. Burigana, P. de Bernardis, J. Delabrouille, E. Dimastrogiovanni, F. Finelli, D. Fixsen, L. Hart, C. Hernández-Monteagudo, J. C. Hill, A. Kogut, K. Kohri, J. Lesgourgues, B. Maffei, J. Mather, S. Mukherjee, S. P. Patil, A. Ravenni, M. Remazeilles, A. Rotti, J. A. Rubiño-Martin, J. Silk, R. A. Sunyaev, and E. R. Switzer, “New horizons in cosmology with spectral distortions of the cosmic microwave background,” Exp. Astronomy (May, 2021) , arXiv:1909.01593 [astro-ph.CO].
  * [61] J. Silk, “Cosmic black body radiation and galaxy formation,” Astrophys. J. 151 (1968) 459–471.
  * [62] R. A. Sunyaev and Y. B. Zeldovich, “Small scale entropy and adiabatic density perturbations - Antimatter in the Universe,” Astrophys. and Space Science 9 no. 3, (Dec., 1970) 368–382.
  * [63] M. Braglia, X. Chen, and D. K. Hazra, “Probing Primordial Features with the Stochastic Gravitational Wave Background,” JCAP 03 (2021) 005, arXiv:2012.05821 [astro-ph.CO].