Primordial black holes from modified supergravity

The modified supergravity approach is applied to describe a formation of Primordial Black Holes (PBHs) after Starobinsky inflation. Our approach naturally leads to the two-(scalar)-field attractor-type double inflation, whose first stage is driven by scalaron and whose second stage is driven by another scalar field which belongs to a supergravity multiplet. The scalar potential and the kinetic terms are derived, the vacua are studied, and the inflationary dynamics of those two scalars is investigated. We numerically compute the power spectra and we find the ultra-slow-roll regime leading to an enhancement (peak) in the scalar power spectrum. This leads to an efficient formation of PBHs. We estimate the masses of PBHs and we find their density fraction (as part of Dark Matter). We show that our modified supergravity models are in agreement with inflationary observables, while they predict the PBH masses in a range between $10^{16}$ g and $10^{23}$ g. In this sense, modified supergravity provides a natural top-down approach for explaining and unifying the origin of inflation and the PBHs Dark Matter.


Introduction
The prospect that Dark Matter (DM) is composed of Primordial Black Holes (PBHs) is an intriguing and highly motivated alternative to any particle physics explanations such as Weak Interacting Massive Particles (WIMPs), gravitino or axion dark matter. Indeed, such a possibility reverses the strategy for DM phenomenology: DM signals may appear in cosmological data rather than colliders, direct detection searches or indirect detection in astroparticle physics. The idea of PBHs was proposed by Zeldovich and Novikov [1], and then by Hawking [2] who realized that some primordial density fluctuations may lead to PBH seeds in the early Universe. There are several mechanisms that may catalyze the formation of PBHs: (i) gravitational instabilities induced from scalar fields [3] such as axion-like particles or multi-field inflation, (ii) bubble-bubble collisions from first order phase transitions (see Refs. [4,5,6] for recent discussions), and (iii) formation of critical topological defects such as cosmic strings [7] and domain walls [8,9] in the early Universe.
After accretion, some PBHs may survive in the Universe today and provide candidates for (non-particle) Dark Matter (DM) [10]. More recently, PBHs attracted considerable attention in the literature, related to observational progress in lensing, cosmic rays, and Cosmic Microwave Background (CMB) radiation, see Refs. [11,12] for a review of observational constraints on PBHs and their prospects for being a fraction of or a whole DM.
On the theoretical side, PBHs are considered as a probe of very high energy physics and quantum gravity "even if they never formed" [13]. Numerous phenomenological scenarios were proposed for PBH formation and, especially, for PBH generation after inflation in the early Universe, under the assumption that PBHs contribute to DM (see e.g., Refs. [14,15,16,17,18,19,20] and the references therein). Indeed, the whole PBH DM case leaves only two limited windows for allowed PBH masses around either 10 −15 or 10 −12 of the Solar mass.
Therefore, It is of interest to study a possible theoretical origin of PBHs at a more fundamental level than General Relativity (GR) by using string theory, as a candidate for quantum gravity, and supergravity as the first step in that direction. Moreover, because of the constraints imposed by local supersymmetry on possible couplings, a viable description of PBHs in supergravity may lead to significant discrimination of phenomenological models of inflation and PBHs.
Due to the absence of large non-Gaussianities and isocurvature perturbations in the current observational CMB data [21], single-field inflationary models were distinguished and discriminated within the large landscape of inflation mechanisms. The Starobinsky R 2 inflation [22] seems to be favored as the best phenomenological fit. Then a required growth (by a factor of 10 7 compared to the CMB amplitude) of the amplitude of fluctuations to be responsible for PBH seeds can be achieved by modifying the inflaton scalar potential with a nearly inflection point [23,24,25]. Details of the PBH production after single-field inflation are very much dependent upon a choice of inflaton potential. This requires a significant fine-tuning for PBHs as a candidate of DM. Standard (Einstein) supergravity can accommodate single-field inflationary models in the (new) minimal setup, with the only restriction to the inflaton potential as a real function squared [26,27,28,29]. 1 Since there are no fundamental reasons for the absence of non-Gaussianities and isocurvature perturbations (they just have to be below the observational limits), multi-field inflationary models were also extensively studied. The required growth of primordial fluctuations can be achieved by tachyonic instabilities, say, in the waterfall phase of hybrid inflation [31,32]. Moreover, the PBH production may be a generic feature of two-field inflation [33]. On the other side, multi-field inflation considerably extends a number of physical degrees of freedom and possible interactions, which reduce predictive power.
Thus, supersymmetry is expected to be even more important in multi-field inflation by limiting the number of fields involved (in the minimal setup) and severely restricting their interactions.
As a guiding principle, in this paper we elaborate on a possible "supergravitational" origin of both inflation and PBHs, by using only supergravity fields and their locally supersymmetric interactions, without adding extra matter fields. Only the minimal number of the physical degrees of freedom associated with an N = 1 full supergravity multiplet is used. In its spirit, our approach is similar to Starobinsky inflation based on gravitational interactions only (see Ref. [34] for a recent review of Starobinsky inflation in gravity and supergravity). The Starobinsky inflation is based on the modified (R + ζR 2 ) gravity, which can be further extended to modified supergravity in the minimal setup [35,36,37] leading to the effective two-field double inflation. We will show that the emerging double-field inflationary model is suitable for a formation of PBH seeds after the first inflation. In this sense, Starobinsky supergravity naturally relates inflation with the dark matter genesis.
Our paper is organized as follows. In Sec. 2 we introduce the general modified supergravity setup and give a specific example of the bosonic terms arising in the simplest non-trivial model. In Sec. 3 we introduce the duality transformations between the modified supergravity and the standard supergravity (in Jordan and Einstein frames) in terms of the field components (of the bosonic part) and in terms of the superfields. In Sec. 4 we study the vacuum structure of our basic model and the effective inflationary dynamics of its two scalars. Section 5 is devoted to an investigation of two-field inflation in our basic model defined by keeping only the leading terms in a generic modified supergravity action. We demonstrate consistency of the basic model with CMB observations but also find the necessity of extreme fine-tuning of initial conditions for PBH generation. In Sec. 6 we extend our basic model by two subleading terms within the same modified supergravity framework, and study in detail the two modifications of our basic model, corresponding to activation of only one of the two subleading terms. We numerically compute the power spectra, and estimate PBH masses and their density fraction, in both cases. We find that our extended models are capable to simultaneously describe viable (Starobinsky-type) inflation and PBH production after inflation, with limited fine-tuning of the parameters, exhibiting an attractor-type behavior. In Sec. 7, we give our conclusions and comments. Some technical details are summarized in Appendices A and B.

Modified supergravity setup
Let us consider a modified supergravity theory with the general Lagrangian (in curved superspace of the old-minimal supergravity in four spacetime dimensions, with M Pl = 1) [38,36] which is parametrized by two arbitrary functions, a non-holomorphic real potential N and a holomorphic potential F, of the covariantly chiral scalar curvature superfield R of the old-minimal supergravity. 2 Some relevant details about supergravity in superspace are collected in Appendix A. It should be mentioned that the master Eq. (1) goes beyond the supergravity textbooks and describes a modified supergravity because the standard (Einstein) supergravity actions are the extensions of Einstein-Hilbert term, whereas Eq. (1) is more general and reduces to the pure Einstein supergravity action only in the very special

Superfield dual version
As was demonstrated in Ref. [36], the dual superfield Lagrangian of Eq. (1) is obtained by introducing the Lagrange multiplier (chiral) superfield T as 4 Varying it with respect to T gives back the original Lagrangian (1) by identifying the chiral superfield S with R. When using instead the superspace identity the Lagrangian (16) can be rewritten to Given the functions N (R, R) and F(R) according to Eq. (2), the Lagrangian (18) can be rewritten to the standard form, where the Kähler potential K and the superpotential W are given by (cf. Ref. [40]) after rescaling S → M S/2 and using the parameter ζ ≡ M 4 ξ/144. It is straightforward to derive the corresponding bosonic terms in field components. We find where Φ i = (T, S), i = 1, 2, and the Kähler metric reads with P ≡ T + T −Ñ . The inverse Kähler metric is given by Because of the non-vanishing non-diagonal elements of the Kähler metric, the kinetic part of the Lagrangian mixes the derivatives of S and T , In order to bring the Lagrangian to the form (14), where the contributions of b m 5 and the angular part of X are ignored, we set ImT = 0, S = |S|, and denote |S| = σ/ √ 6. The kinetic mixing between ∂S and ∂T can be eliminated by using P = T + T −Ñ as the independent (real) scalar instead of ReT . We get the canonical normalization of its kinetic term in the parametrization P = exp 2 3 ϕ as follows: that exactly matches the kinetic part of Eq. (14). It is also straightforward (albeit tedious) to check that the scalar potential of Eq. (22) also coincides with that of Eq. (15) after using the field redefinitions diagonalizing the scalar kinetic matrix above.

Critical points
To study vacuum equations in our basic model, we denote e − 2 3 ϕ ≡ x and rewrite the scalar potential (15) as The equations for critical points read where the primes denote the derivatives with respect to σ. A simple solution to these equations is It gives rise to the vanishing potential (27) for σ 0 = ϕ 0 = 0. There is another solution by taking U = 0 and obtaining On the other hand, the condition U = 0 is solved by Equating Eqs. (31) and (32) leads to an equation on the parameter ζ with a solution ζ = 1/54 ≈ 0.019 provided that the "plus" branch is chosen in Eq. (31). It means, when ζ = 1/54, we have three Minkowski minima: σ 0 and ±|σ 1 | where σ 2 1 is given by Eqs. (31) or (32). When ζ = 1/54, the two minima at ±σ 1 are not given by Eqs. (31) and (32), being more general solutions to the vacuum equations (28) and (29). In particular, when 0 < ζ < 1/54, the minima at ±σ 1 are AdS, while for 1/54 < ζ < 0.027 the minima are uplifted to metastable de-Sitter (dS). When ζ ≈ 0.027, there are two inflection points, whereas for ζ > 0.027 all the critical points, except of σ = 0, disappear. The scalar potential V /M 2 is shown in Fig. 1a (at ζ = 1/54) and Fig. 1b (at ζ = 0.027).
Let us comment on the scalar masses for the model (20) (21). Expanding around the Minkowski vacuum at ϕ = σ = 0, we find M ϕ = M σ = M , where M ϕ and M σ are the masses of ϕ and σ, respectively. As for the axion ImT , after its proper normalization we find that it also has the mass M . However, the last scalar θ has vanishing mass around the minimum, so it must be generated by additional means. This, together with the more detailed analysis of the dynamics of ImT and θ deserves a separate investigation that we leave to future works. Here we will focus on the two scalars ϕ and σ.

Two-field inflationary dynamics
Having derived the Lagrangian with the two-field scalar potential from the modified supergravity, in this Section we investigate its suitability for describing cosmological inflation in agreement with CMB observations.

Field equations
The Lagrangian in Eqs. (14) and (15) takes the form of a Non-Linear Sigma-Model (NLSM) minimally coupled to gravity, where φ A = {ϕ, σ}, A = 1, 2, and the NLSM metric is given by Varying the Lagrangian (33) with respect to the scalar fields yields equations of motion in the form where ≡ ∇ m ∇ m is the spacetime Laplace-Beltrami operator, and Γ C AB are the Christoffel symbols of the NLSM target space. The non-vanishing Christoffel symbols are After using these results, and the Friedmann-Lemaitre-Robertson-Walker (FLRW) spacetime metric g mn = diag(−1, a 2 , a 2 , a 2 ) with the time-dependent scale factor a(t), the equations of motion take the form where the dots stand for the time derivatives. The Friedmann equations for the system (33) read where the Hubble function has been introduced, H ≡ȧ/a. For numerical computations it is useful to rescale time ast ≡ M t (when usingt, the dots will denote the derivatives with respect tot) with the rescaled Hubble functionH = H/M .

Inflationary parameters
In this Subsection we employ the covariant formalism that is well known in the literature, see e.g., Refs. [42,43,44] and the references therein, with the slow-roll parameter In a two-field analysis, it is useful to define the field-space velocity and acceleration (turn rate) unit vectors as , respectively, where the absolute value of a field-space vector a A is defined by |a| ≡ √ G AB a A a B , and the acceleration vector ω A is defined by Another useful quantity is the effective mass matrix, where R A CDB is the Riemann tensor of the NLSM scalar manifold, with the non-vanishing components With the above definitions we can introduce the adiabatic and isocurvature parameters respectively, where η ΣΣ plays the role of the second slow-roll parameter, while η ΩΩ is proportional to the effective isocurvature mass.
The transfer functions are defined as follows: where The transfer functions describe the evolution of perturbations on superhorizon scales, i.e. from the moment of horizon exit t 1 (of the k-mode of interest) until some later time t 2 .
The inflationary observables (CMB tilts) can be computed as (by assuming that isocurvature modes are suppressed) n s = 1 − 6 + 2η ΣΣ and r = 16 As T RS is real, the maximum value of the tensor-to-scalar ratio is r max = 16 , while it can be computed without the transfer functions. In Appendix B we estimate both T RS and T SS and find them negligible. Therefore, isocurvature effects can be ignored at CMB scales indeed.
According to the latest PLANCK data [21], the observed values of n s and r are n s = 0.9649 ± 0.0042 (1σ CL) and r < 0.064 (2σ CL) .

Inflationary solutions
Let us first consider the case of ζ = 1/54 ≈ 0.019 with three Minkowski minima in Fig. 1a. We numerically solve the field equations (37), (38) and (39) with the initial conditions ϕ(0) = 6, σ(0) = 3, and the vanishing initial velocities, so let us call it the solution (I). The scalar field solutions are plotted in Figure 2a, and their trajectories in the scalar potential are plotted in Figure 2b. It can be seen that σ quickly drops to its minimum σ = 0, so that the trajectory becomes similar to that in the single-field Starobinsky inflation. In fact, this is a generic feature when the initial velocities are zero (or almost zero), ϕ(0) 6 and |σ(0)| σ max , where σ max = 1/ √ ζ is the upper bound on σ where the potential is infinite. When ζ = 1/54 we find σ max ≈ 7.35.
The solution (I) leads to the spectral tilt and the tensor-to-scalar ratio as n s ≈ 0.9662 and r max ≈ 0.003, respectively, which are consistent with the observed values and the theoretical (Starobinsky) predictions of chaotic single-field inflation.
As the initial value σ(0) approaches σ max and/or as the initial velocities become nonnegligible, the trajectory starts to curve. Also a smaller value of ϕ(0) makes it easier to curve the trajectory.
As regards PBH production after inflation, let us consider the field-space trajectory going through the saddle point of the potential, that is a maximum in the σ-direction and a   (37) and (38) with the initial conditions ϕ(0) = 6, σ(0) = 3 and the vanishing initial velocities. The blue shaded region represents the time period of the last 60 e-foldings.
(local) minimum in the ϕ-direction. Then the saddle point divides inflation into two stages.
We found a set of initial conditions that leads to such trajectory with Let us call the corresponding solution as the solution (II). We include its plots in Figures  3a and 3b. The time-dependence of the Hubble function and the e-foldings number defined byṄ = H are shown in Figures 3c and 3d. The total number of e-foldings is around 40, though it can be larger for larger values of ϕ(0) with more fine-tuning of the initial velocities. Thus, in order to achieve the two-stage inflation, where the field-space trajectory passes through the saddle point, we have to fine-tune the initial conditions as in Eq. (52), though the last choice is not unique. The reason is, when ϕ is large, the potential takes the shape of a valley with the minima at σ = 0, so a generic behavior of σ is to quickly relax at σ = 0, and let ϕ drive the entire inflationary period. The same remains true if we change the shape of the potential as in Figure 1b by changing ζ (the only difference is the saddle point to be replaced by an inflection point).

Generalized attractor-type models
Having learned the lessons in the previous Sections, we conclude that our basic ansatz in Eq. (2) for functions N (R, R) and F(R) is too restrictive because it requires extreme fine tuning of the initial conditions for PBHs production. Therefore, we generalize our ansatz by adding the next-order corrections as where we have introduced two new parameters γ and δ with their normalization chosen for later convenience. We keep α = 0, β = −3 and ζ ≡ M 4 ξ/144, ignore b m and the angular mode of R| = X, and set X = M σ/ √ 24, as in the previous Sections. In the framework of the dual matter-coupled supergravity (19), the γ-term resides in the Kähler potential that can be affected by quantum corrections, whereas the δ-term resides in the superpotential that does not receive (perturbative) quantum corrections.     (37) and (38) with the initial conditions (52). The red spot in (3b) represents one of the two saddle points of the potential.
After repeating the procedure outlined in Subsection 3.1, we obtain the Einstein frame Lagrangian as follows: where the functions A, B, U are given by When γ = δ = 0, all that reduces to Eqs. (14) and (15), as it should. Similary to the basic model, there is the infinite wall in the scalar potential, which prevents σ from obtaining values leading to the wrong sign of its kinetic term. The relevant field equations of the generalized model are The generalized model defined by Eqs. (53) and (54) appears to be rather complicated for a detailed numerical analysis, so we study only two special cases, the one with δ = 0 (dubbed the γ-extension) and the one with γ = 0 (dubbed the δ-extension), in what follows.

The γ-extension
As a representative of the γ-extension (δ = 0), we choose the parameters γ = 1 and ζ = −1.7774, see Fig. 4. This choice is interesting because the scalar potential (for ϕ 1) has two valleys where σ = 0, and a single Minkowski minimum at σ = ϕ = 0. The first Slow-Roll (SR) inflation is possible along either of the valleys. The valleys merge into the Minkowski minimum by passing through inflection points (or near-inflection points) followed by the second, Ultra-Slow-Roll (USR), inflationary stage 6 .
After numerically solving the equations of motion (57)-(60) we plot the solutions in Fig. 5. The total number of e-foldings is set to ∆N = 60, and the end of the first stage of inflation is defined by the time when η ΣΣ first crosses unity (see Fig. 5e). We could also define the end of the first stage from the local maximum of , which nearly coincides with the former definition. As may be expected from the USR period USR SR and can be seen in Fig. 5e, it leads to an enhancement in the scalar power spectrum indeed. Inflation ends when = 1, as usual. With the chosen parameters, the first stage lasts ∆N 1 ≈ 50 e-foldings, whereas the second stage lasts for ∆N 2 ≈ 10: in the subsequent Figures the first stage of inflation is represented by the blue shaded region, whereas the second stage is marked by the green shaded region, whenever is relevant. The length of the second stage is controlled by the parameter ζ for a given γ. 6 Actually, despite the name, during an ultra-slow-roll regime, the scalar field(s) roll down the potential faster than during a slow-roll regime (see e.g. Ref. [45]).
The parameter space. The parameter choice leading to a scalar potential with the suitable properties (as described above) is not unique, and for any γ greater than ∼ 0.004 there is a value of ζ that leads to a similar shape of the potential (with two inflection points, unique Minkowski minimum, etc.). For a given γ, one can solve the system of equations where H is the Hessian determinant of the potential, in order to obtain the value of ζ leading to the desired inflection points. Then, by fine-tuning ζ around that value, one can change a duration of the USR stage ∆N 2 .
In order to see how γ changes the shape of the scalar potential, let us evaluate the ratio V inflec. /V ∞ as a function of γ, where V inflec. is the value of the potential at an inflection point, and V ∞ is the asymptotic value of the potential when ϕ → ∞ and σ is at its local minimum, which corresponds to the SR stage. This ratio represents the depth of the inflection points relative to the SR valleys, and it does not significantly change the curvature of the inflationary path in the ϕ − σ plane. The plot of V inflec. /V ∞ versus γ, as well as the trajectory in the γ − ζ plane, which solves Eqs. (62), are shown in Fig. 6. After taking all that into account, we conclude the control over the overall shape of the potential is limited due to the attractor-type behavior of V inflec. /V ∞ at large γ.
When γ = 1, a solution to Eq. (62) gives ζ ≈ −1.774 that, in turn, leads to ∆N 2 ≈ 6.3. However, in our example we slightly departed from that value of ζ and set ζ = −1.7774 in order to obtain ∆N 2 ≈ 10. Our strategy is to compute power spectra for different choices of γ while keeping ∆N 2 ≈ 10. The latter condition fixes the value of ζ. Next, we examine the impact of a variation of ∆N 2 .
The power spectrum at fixed ∆N 2 . We numerically compute the power spectrum of curvature perturbations by using the transport method introduced in Refs. [46,47] with the Mathematica package described in Ref. [48]. We compute the spectrum around the pivot scale k * that leaves the horizon at the end of the first stage, i.e. ∆N 2 e-folds before the end of inflation (let us call this scale k ∆N 2 ). The inflaton mass is adjusted in each case around 0.5 × 10 −5 M Pl by requiring P ζ ≈ 2 × 10 −9 for the mode k that exits the horizon 60 e-folds before the end of inflation -call it k 60 .
The power spectrum for various values of γ is shown in Fig. 7. The parameters considered are collected in Table 1, where ζ is tuned to satisfy ∆N 2 ≈ 10. A change of n s and r max (still given by Eq. (61)) is negligible for those parameters.  Table 1: The parameters used in a computation of the power spectrum in Fig. 7 with ∆N 2 ≈ 10.
As is often adopted in the literature, the desired enhancement of primordial curvature perturbations should be at least P ζ,max. P ζ,min. ≡ P enh. ∼ 10 7 (where the subscripts max. and min. refer to the values at the peak and at the base of the peak, respectively) compared to CMB scales, in order to efficiently produce PBHs. However, the authors of Ref. [15] argued by using peak theory that, given a broad peak, the required enhancement in the power spectrum drops by one order of the magnitude to ∼ 10 6 . According to our numerical estimates (see Fig. 7), when choosing ∆N 2 = 10, the enhancement of the order P enh. 10 6 is achieved for γ 10, namely, it ranges from the order of 10 6 at γ = 10 to the order of 10 7 at γ = 1000. However, the power spectrum enhancement also grows with ∆N 2 (see below), so that smaller values of γ may also become viable for sufficiently large ∆N 2 .
Changing ∆N 2 . Let us examine how the power spectrum changes with ∆N 2 . To demonstrate that dependence, we consider the power spectrum at γ = 0.1 and γ = 1 with various values of duration of the USR stage, ∆N 2 = 10, 20, 23, 26 for each γ. The results are collected in Fig. 8. For example, when γ = 0.1 and ∆N 2 = 26, the enhancement can be estimated as P enh. ≈ 3.2 × 10 5 . Hence, in order to achieve the enhancement of the order of 10 6 , we must require ∆N 2 30. But this would push the spectral index outside of the 3σ (lower) limit of n s ≈ 0.946 (cf. Refs. [49,14]). On the other hand, at γ = 1 we have P enh. ≈ 1.4 × 10 8 already at ∆N 2 = 20 and, therefore, the values of γ 1 are favored for efficient production of PBHs.
In Table 2 we collect the approximate values of n s and r max. (at the CMB scales) for the values of ∆N 2 = 10, 20, 23, 26, universally across the considered values of γ = 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σ limit when ∆N 2 = 20, and is marginally outside the 3σ limit when ∆N 2 = 23.
PBH masses and their density fraction. The mass of a PBH created by lateinflationary overdensities was estimated in Ref. [14] as follows:   where t peak is the time when the wavenumber corresponding to the power spectrum peak (k peak ) exits the horizon, whereas t 60 is the time when k 60 exits the horizon (the beginning of observable inflation). The formula is independent of the period between t peak and the time of PBHs formation during the radiation-dominated era. We estimate the values of M PBH for various values of ∆N 2 by using Eq. (63) The results are shown in Table 3 together with the corresponding values of the spectral index. Our estimates are universal across the values of γ = 0.1, 1, 10, 100, 1000. PBHs with masses smaller than ∼ 10 16 g would have already evaporated by now via Hawking radiation. Thus, we require ∆N 2 > 20. On the other hand, the lower 3σ limit on the spectral index, n s ≈ 0.946, requires ∆N 2 < 23, so that viable PBH masses are restricted by O(10 16 g) < M PBH < O(10 19 g) even before considering observational constraints on PBHs masses.  As regards the constraints on γ, the power spectrum in Fig. 8 tells us for ∆N 2 > 20 that it is sufficient to have γ O(1) in order to produce the required enhancement in the spectrum.
We also estimate the PBHs density fraction by using Press-Schechter formalism [50]. The useful formulae include the PBH massM PBH (k), the production rate β f (k), and the density contrast σ(k) coarse-grained over k as follows (see e.g., Refs. [51,52] and references therein):M PBH 10 20 7 × 10 12 k Mpc respectively, where we have chosen the Gaussian window function for the density contrast and have introduced δ c as a constant representing the density threshold for PBH formation, which is usually estimated as δ c ≈ 1/3 [53] for simplicity (its more precise value depends upon details of the power spectrum). In terms of the above functions, the PBHs-to-DM density fraction can be estimated as follows [51,52]: In order to numerically evaluate the functions (6.1) and (65), we need to normalize the values of k in terms of the observable scales today. We choose the scale k 60 to represent the largest observable scale today, which is around 10 −4 Mpc −1 . 7 Then a numerical evaluation shows, in order to obtain a substantial density fraction, we need a relatively small δ c (of course, given a fixed δ c , we can instead adjust the parameters of the model to obtain the desired density fraction). For example, with the parameters γ = 1 and ∆N 2 = 22, setting δ c = 0.275 results in the PBHs fraction shown in Fig. 9, that we have superimposed on the observational constraints of Ref. [54] (see also Ref. [55]). According to Fig. 9, our peak is located at the center of the low-mass window of 10 16 − 10 23 grams. The constraints of Fig. 9 are imposed by assuming a monochromatic PBH mass spectrum. In our case the mass distribution is narrow, albeit not strictly monochromatic.
It should be noticed that the exact location of the peak in Fig. 9 depends on the normalization of the wavenumber k that, in turn, depends on the amount of post-inflationary expansion of the universe.

The δ-extension
As the next possibility, we study the δ-extension in this Subsection, when γ = 0 and δ = 0 in Eqs. (53) and (54). 8 It breaks the R-symmetry and the reflection symmetry σ → −σ of the potential, see Fig. 10. For any non-zero δ, there is a value of ζ that leads to an inflection point: for a positive δ the inflection point is at σ = −|σ inflec | (as in our example of Fig. 10), and for a negative δ the inflection point is at σ = +|σ inflec |.
In contrast to the γ-extension, here we have a single valley for large positive ϕ and σ = 0, so that in this limit the model reduces to a single-field Starobinsky model. As one approaches ϕ = 0, the potential inclines towards the (near-)inflection point which could, in principle, guide the inflationary trajectory towards passing through the (near-)inflection point before falling to the Minkowski minimum at ϕ = σ = 0.   Let us consider, for example, the parameter values δ = 0.1 and ζ = 0.033407 (ζ is chosen to get ∆N 2 = 10). After solving the corresponding field equations, we show the time dependence of ϕ, σ,H, N , and η ΣΣ in Fig. 11. The near-inflection point divides inflation into two stages with ∆N 1 = 50 (slow-roll) and ∆N 2 = 10 (ultra-slow-roll). We set initial velocities to zero, with ϕ(0) = 6 and σ(0) = 0.05. Similarly to the γ-extension, the inflationary trajectory is stable against variations of the initial conditions, as long as they are not very large.
The parameter space. When demanding the presence of a (near-)inflection point, the parameters must satisfy Eq. (62). The plot of V inflec. /V ∞ (V ∞ is taken for ϕ 1 and σ = 0) versus δ, and the solution ζ(δ) to Eq. (62), are displayed on the left side of Fig. 12. On the right side of Fig. 12 we show the profile of the potential with ϕ at its local minimum satisfying ∂ ϕ V = 0, for several choices of δ. In particular, our plot shows, when δ → 0, the inclination of the potential towards the inflection point becomes smaller until it vanishes when δ = 0 (in such case the potential coincides with the one in Fig. 1b). Therefore, when δ is very small, the inclination of the potential becomes insufficient for guiding the inflationary trajectory through the inflection point. Instead, the trajectory tends to the σ = 0 path (when δ = 0, the trajectory exactly follows the σ = 0 path).
We plot the inflationary trajectories in the ϕ − σ plane for various choices of δ (with ζ being fixed by requiring the existence of an inflection point) in Fig. 13. The colored spots represent the inflection points. As can be seen in Fig. 13 for δ = 0.01 and δ = 0.05, the trajectory misses the corresponding inflection point by a large margin and, therefore, avoids the USR regime. On the other hand, when δ = 0.1, the trajectory stops near the inflection point and then oscillates a few times before going to the minimum at ϕ = σ = 0. This indicates the possibility of an USR stage, and it happens in Fig. 11 for this parameter  The scalar power spectrum at fixed ∆N 2 and δ 0.1. Let us fix ∆N 2 = 10, and consider the power spectrum for several values of δ. We find that the spectrum has a non-trivial dependence on δ, see Fig. 14. In the left plot, δ is varied from 0.1 to 0.2, and we observe the spectrum enhancement to become smaller as δ grows. In the right plot, once δ reaches 0.2, the enhancement starts growing with δ and develops a sharper peak.
As regards larger values of δ, our numerical results show, when δ 0.6, it becomes increasingly more difficult to maintain the USR stage and to achieve ∆N 2 > 10, in particular. It may be due to the need of an extreme fine-tuning of the parameter ζ when δ is large.
Changing ∆N 2 . Amongst the values of δ studied above, let us pick up those with the highest power spectrum peaks, namely, δ = 0.1 and δ = 0.6, and plot the spectrum for ∆N 2 = 10, 20, 23, 26. The results are displayed in Fig. 15 with the plots on the left side and the right side corresponding to δ = 0.1 and 0.6, respectively. As expected, the enhancement becomes larger with increasing ∆N 2 . PBH masses and their density fraction. To be specific, let us consider two different examples: a smooth peak for δ ∼ 0.1, and a sharp peak for δ ∼ 0.6, in the power spectrum. Requiring the PBH density fraction (65) to be around unity and the corresponding density threshold not to deviate far away from the region 1/3 ≤ δ c ≤ 2/3, we find the choices δ = 0.094 and δ = 0.58 are suitable for efficient generation of PBHs with their masses around 10 19 g, while avoiding their overproduction (f 1).
We estimate the PBH masses by using Eq. (63) and summarize our results in Table  4. Hence, for M PBH ∼ 10 19 g we can either set δ = 0.094 and ∆N 2 = 20 (in this case P enh ≈ 4.5 × 10 7 ), or δ = 0.58 and ∆N 2 = 23 (in this case P enh ≈ 2.7 × 10 8 ). In the former case, n s is within 2σ CL, whereas in the latter case, n s is within 3σ (but outside 2σ) CL. In those two examples, we compute the PBH density fraction given by Eq. (65) and show the results in Fig. 16. The observational constraints in Fig. 16 are included for the reference purposes only, as they assume the monochromatic PBH mass function.

Conclusions and comments
In this paper, we analyze several supergravity extensions of the Starobinsky inflationary model. We explore possibilities of PBHs genesis that could account for part of Cold Dark Matter. We find that PBHs generation can be efficiently catalyzed by primordial perturbations sourced by the Starobinsky scalaron coupled to a new supersymmetric "modulus" (scalar) field. Let us summarize our strategy. We rely on theoretical considerations before comparing them with cosmological observations, as a top-down approach. As our starting point, we adopt the Starobinsky inflationary model serving as the theoretical tool and pointing out the need of modified gravity in a more fundamental approach, i.e. beyond considering the Starobinsky model as merely the best phenomenological fit to CMB observations. We extend the modified (R + ζR 2 ) gravity to the modified supergravity, where the latter is considered as the candidate (or as the approximation) of a more fundamental theory of quantum gravity. Amongst the theoretical advantages of modified supergravity are (i) the use of supergravitational couplings only, (ii) predicted new physical degrees of freedom, and (iii) its formal equivalence to the standard (matter-coupled) supergravities. However, unlike the standard supergravities coupled to matter, modified supergravity can be limited to the supergravity fields alone, where the new physical scalar naturally accompanies the inflaton (scalaron), together with metric and gravitino. It happens because modified supergravity is a higher-derivative theory, so that the "auxiliary" scalar of the standard (off-shell) supergravity multiplet becomes dynamical. We find that modified supergravity naturally leads to the two-field inflationary models with restricted couplings and a small number of free parameters. Therefore, local supersymmetry has predictive power for phenomenology of the early universe cosmology via the double-inflation scenario. Indeed, the second field coupled to the Starobinsky scalaron is not introduced ad hoc but is predicted by the supergravity extension of the Starobinsky model. Our strategy is to use those models for a viable description of Starobinsky inflation together with the PBH production after inflation. Cosmological inflation and the PBH production can be considered as probes of supergravity for its use as a more fundamental approach, and vice versa: modified supergravity provides a theoretical input for the discrimination of phenomenological models of inflation and PBHs.
We summarize our main results as follows.
A generic modified supergravity Lagrangian in the manifestly supersymmetric form (with all couplings included) is given by Eq. (1). After (Taylor) expanding its potentials N and F in powers of the scalar curvature superfield R and keeping only the leading terms (needed for minimal embedding of R + ζR 2 gravity), we arrive at our basic model defined by Eq. (2), whose relevant bosonic terms (in Jordan frame) are given by Eqs. (3) and (4). As the next step, we perform the duality transformation of the derived bosonic terms to Einstein frame, and arrive at the two-scalar NLSM minimally coupled to gravity with the derived NLSM metric and the scalar potential, given by Eqs. (14) and (15). We also provide the manifestly supersymmetric (complete) duality transformation in terms of the superfields, and compute the corresponding Kähler potential and the superpotential, given by Eqs. (20) and (21) in the case of the basic model as an example. Then, we study the critical points (vacua) of the derived scalar potentials and the inflationary dynamics of two scalars in the context of two-field inflation, and we find consistency of the basic model with CMB observations. However, we also observe that such a scenario can work only with an extreme fine-tuning of initial conditions for efficient formation of PBHs.
To overcome the problem, we add the next (subleading) terms to our basic model within the same modified supergravity master Lagrangian (1). There are two such terms, see Eqs. (53) and (54), so we study them separately. We numerically compute the power spectra, estimate PBH masses and their density fraction, in both cases. We find that any of the extended models can simultaneously describe viable (Starobinsky-type) inflation and the PBH production after inflation, with limited fine-tuning of the parameters, exhibiting an attractor-type behavior. Actually, the PBH production is less sensitive to changes of the parameter γ in the γ-extension of the N -potential. Thus, the γ-term seems to be preferable over the δ-term. We confront our theoretical predictions for PBHs (as part of DM) with current observations in Figs. 9 and 16; in the cases of the γ-and δ-extensions, respectively. Having both the γ-and δ-terms in the Lagrangian is expected to render our supergravity model even more flexible in accommodating the PBH DM, being compatible with the constrains on the inflationary parameters.
Of course, modified supergravity does not pretend on the status of an ultimate fundamental theory. However, there are indications that it may be embedded into superstrings considered as an ultra-violet complete theory of quantum gravity. Here it is worthwhile to mention that (i) modified supergravity always leads to the no-scale Kähler potential (20) that often arises in superstring compactifications (see e.g., Ref. [57]), and (ii) there is a possibility of interpreting (some) modified supergravity theories as the D3-brane worldvolume theories in type II superstrings [58,59]. Thus, the exploration of cosmological predictions from modified supergravity provides a remarkable bridge between quantum gravity on one side and phenomenology of inflation and PBHs on the other side.
PBH formation necessarily leads to Gravitational Waves (GWs) because large scalar overdensities act as a source for stochastic GWs background. Frequencies of those GWs can be directly related to expected PBHs masses and duration of the second stage of inflation [60]. Those GWs may be detected in the future ground-based experiments, such as the Einstein telescope [61] and the global network of GWs interferometers including advanced LIGO, Virgo and KAGRA [62], as well as in the space-based GWs interferometers such as LISA [63], TAIJI (old ALIA) [64], TianQin [65] and DECIGO [66].