Multi-Field versus Single-Field in the Supergravity Models of Inflation and Primordial Black Holes

We review the models unifying inflation and Primordial Black Hole (PBH) formation, which are based on the modified (Starobinsky-type) supergravity. We begin with the basic (Starobinsky) inflationary model of modified gravity and its alpha-attractor-type generalizations for PBH production, and recall how all those single-field models can be embedded into the minimal supergravity. Then, we focus on the effective two-field models arising from the modified (Starobinsky-type) supergravity and compare them to the single-field models under review. Those two-field models describe double inflation whose first stage is driven by Starobinsky’s scalaron and whose second stage is driven by another scalar belonging to the supergravity multiplet. The power spectra are numerically computed, and it is found that the ultra-slow-roll regime gives rise to the enhancement (peak) in the scalar power spectrum leading to an efficient PBH formation. The resulting PBH masses and their density fraction (as part of dark matter) are found to be in agreement with cosmological observations. The PBH-induced gravitational waves, if any, are shown to be detectable by the ground-based and space-based gravitational interferometers under construction.


Introduction
The oldest signal from the past to the present day is given by the Cosmic Microwave Background (CMB) that emerged when electromagnetic radiation decoupled from matter and began propagating approximately 380,000 years after the Universe was born. That time can be considered as the observational wall against any electromagnetic probe because photons did not propagate freely before that time (due to Thompson scattering). Therefore, extracting observational signals from the earlier Universe may only be possible with gravitational waves or neutrino sources.
Some detailed information about the CMB spectrum and its anisotropies and fluctuations is available due to several satellite missions in the past (COBE, WMAP, and Planck). The theoretical explanation of the CMB observations is provided by assuming an era of cosmological inflation in the very early Universe. The cosmological paradigm of inflation is an assumption of the accelerating expansion of the Universe (the primordial dark energy with a Graceful Exit) during a very short time, which produced the Large-Scale Structure (LSS) we observe today. The energy scale of inflation must be well beyond the electro-weak scale of the Standard Model (SM) of elementary particles. Therefore, inflation can also be viewed as the most powerful particle accelerator in nature, giving us the great (but small) window into high-energy physics well beyond the SM.
The standard (minimal) theoretical mechanisms of inflation employ a (canonical) scalar field called inflaton with its scalar potential defining a single-field inflationary model. The physical nature of inflaton and the origin of its scalar potential are unknown, so

Why Inflation?
A spatially homogeneous and isotropic (1 + 3)-dimensional universe at large scales (beyond 100 Mpc) is described by the Friedmann-Lemaître-Robertson-Walker (FLRW) metric where the function a(t) is called the cosmic scale factor in co-moving spacetime coordinates (t, r, θ, φ), k is the FLRW topology index, k = (−1, 0, 1), and dΩ 2 = dθ 2 + sin 2 θdφ 2 . The FLRW metric (1) admits the six-dimensional isometry group G that is SO(1, 3), E(3), or SO(4), respectively, acting on the orbits of the homogeneous space G/SO (3). The spatial three-dimensional sections are H 3 , E 3 , or S 3 , respectively. The Weyl tensor of a FLRW metric vanishes for any function a(t). It is safe to take k = 0 for most applications. Cosmological inflation is referred to the very early universe (after 10 −36 s or below 10 16 GeV, but well before 10 −10 s or well above the electro-weak scale 10 2 GeV), and it is usually defined as the primordial dark energy, i.e., as the accelerating universe with (we use dots to denote the derivatives with respect to the time t) •• a (t) > 0 (2) and a quick (Graceful) exit. Unlike the present dark energy (without an exit), it is more useful to define inflation by the inequality d dt in terms of Hubble function H = • a /a and Hubble radius H −1 /a. Equation (3) means the decreasing Hubble radius (or the causal distance) that is just needed to explain the origin of the FLRW universe after inflation due to homogenization of initial fluctuations.
There is a compelling theoretical and (indirect) observational evidence for cosmological inflation in the early Universe, so that we do not discuss its possible alternatives here. The idea of inflation originated as a solution to theoretical puzzles of the Standard Cosmology, because inflation predicts homogeneity of our Universe at large scales, its spatial flatness, its large size and entropy, and the almost scale-invariant spectrum of cosmological perturbations [10,11]. In more recent times, the cosmological paradigm of inflation received remarkable support from the observations and measurements of the Cosmic Microwave Background (CMB) radiation by the COBE, WMAP, and PLANCK satellite missions, see, e.g., the recent data from PLANCK [12].
Inflation can also be considered as the great amplifier of microscopic quantum fluctuations in the very early universe, leading to macroscopic seeds of the LSS formation. The standard (single-field) mechanism of inflation in gravitational theory uses a scalar field (called inflaton), whose potential energy drives inflation. The single-field (canonical inflaton ϕ) action (in Einstein frame) reads The inflaton scalar potential V(ϕ) should be flat enough to meet the slow-roll conditions during the inflationary era. However, the physical nature and fundamental origin of inflaton and its interactions to other particles are unknown.
The single inflaton action (4) is to be compared to a multi-field action for describing inflation, having the form (in Einstein frame) in terms of several (real) scalar fields φ a , a = 1, 2, . . . , n. The scalar kinetic terms in Equation (5) have the form of the Nonlinear Sigma-Model (NLSM) [13] with a metric G ab (φ) in the field space. The simplest models of inflation and PBH formation are based on the single inflaton action (4), see the next sections. Their motivation is based on simplicity and (expected) low non-Gaussianity of primordial fluctuations. However, there is no good reason that forbids interactions of the inflaton scalar to other fields. Moreover, it is highly welcome, because inflaton is expected to decay into other particles after inflation. Therefore, the action (5) is apparently more suitable for physical applications. It is always possible to embed a single-field inflationary model into a multi-field inflationary model after finetuning of the latter. However, the multi-field inflation has its own problems because of (i) its relative complexity and (ii) the need to fix the potential and the NLSM metric from a more fundamental theory.

Why Starobinsky's Model for Single-Field Inflation?
There are many viable single-field inflationary models with different motivation and inflaton potentials V(ϕ). Instead of reviewing all of them, only one basic model will be described below, and it is known in the literature as the Starobinsky model [2]. Instead of following the references in chronological order, we define the Starobinsky model from the modern perspective by the following modified gravity action: where we have introduced the reduced Planck mass M Pl = 1/ √ 8πG N ≈ 2.4 × 10 18 GeV, and the mass parameter m. We use the spacetime signature (−, +, +, +, ).
At first sight, the action (6) is very different from the action (4), and it looks rather ad hoc. However, it is not the case, while any viable (single field slow roll) inflationary model has to be close to the Starobinsky model, as we are going to demonstrate below.
The action (6) reduces to the standard (Einstein-Hilbert) gravitational action in the low curvature regime. In the high curvature regime relevant to inflation, the action (6) reduces to the no-scale and no-ghost R 2 gravity action with the positive dimensionless coupling constant M 2 Pl /(12 m 2 ). The key result in [2] can be interpreted as the existence of an inflationary (quasi de Sitter) solution to the R 2 gravity in the form of attractor with a Graceful Exit (for H m), also known as chaotic inflation. The (R + R 2 ) gravity model in Equation (6) is a representative of the modified F(R) gravity theories, which are known to be classically equivalent to the scalar-tensor theories of gravity defined by the action (4) after a (Legendre-Weyl) field redefinition, see, e.g., in [14,15] for details of the equivalence derivation in our notation. The specific correspondence is given by where the field redefinition χ(ϕ) is defined by (the primes denote differentiation with respect to χ) The inverse transformation reads (see, e.g., in [16]) defining the function F(R) in the parametric form for a given inflaton scalar potential V(ϕ). The required flatness of the scalar potential (needed for slow roll of inflaton during a sufficient duration of inflation) amounts to the smallness of the first term against the second one in the square brackets of Equations (10) and (11). For instance, simply ignoring the first term (or assuming a constant V) immediately leads to F ∝ R 2 .
The exact inflaton potential of the Starobinsky inflation follows from Equations (6), (8), and (9), and it is given by The potential (12) is a sum of the induced cosmological constant and the exponentially small corrections (for large values of φ). It has a plateau with the positive height that determines the scale of inflation.
The duration of slow roll inflation is usually measured by e-folds where ϕ * is the inflaton value at the horizon crossing, and ϕ end is the inflaton value at the end of inflation when one of the slow roll parameters, becomes no longer small or close to one. The only parameter m of the Starobinsky model is fixed by the observed CMB amplitude (known as the COBE normalization) as follows: The amplitude of scalar perturbations at the horizon crossing reads [17] The Starobinsky model (6) is known (since 1980) to be the excellent model of inflation and is in very good agreement with current measurements. As regards the main cosmological observables, the index of scalar perturbations is given by n s ≈ 1 + 2η V − 6ε V ≈ 0.9649 ± 0.0042 (with 68% CL), while the index of tensor perturbations or the tensor-toscalar ratio is given by r ≈ 16ε V < 0.064 (with 95% CL), according to the PLANCK data [12]. The Starobinsky inflation gives r ≈ 12/N 2 e ≈ 0.004 and n s ≈ 1 − 2/N e , where N e is supposed to be between 50 and 60 [18,19].
Besides the main properties mentioned above (universality for slow roll and agreement with observations), the Starobinsky inflationary model has the following attractive features: • the model is geometrical, being based on gravitational interactions; • its inflaton (called scalaron) is a physical excitation of higher-derivative gravity, so it is not a new fundamental scalar (minimal tool); • the model has no free parameters after fixing its only parameter M by observations; • the Starobinsky inflation obeys the Einstein criterium: "simple but not too simple"; • scalaron can be identified with the Nambu-Goldstone boson related to spontaneous breaking of the scale invariance [20] and the approximate flatness of the potential (12) during slow roll; • the UV-cutoff of quantized (R + R 2 ) gravity is given by M Pl H inf. , as is clear from expanding the non-renormalizable scalar potential (12) in powers of ϕ; this feature ensures reasonable protection of the Starobinsky inflation on the scale H inf. ∼ 10 13 GeV against quantum gravity corrections expected near the Planck scale M Pl ; and • rewriting the scalar potential (12) to the form of a mass term by the field redefinition leads to a non-canonical kinetic term of the field φ, with the pre-factor having a singularity at φ cr. = √ 3/2M Pl and specified by the critical exponent defining the universality class of the Starobinsky inflation.
Regarding the shortcomings of the Starobinsky model, note that • its fundamental origin in quantum gravity or in string theory is unknown; • it cannot be applied for large negative values of the scalar curvature, when F (R) < 0; and • a numerical analysis of (13) with the potential (12) yields (assuming the best fit N e ≈ 55) the trans-Planckian values of the inflaton field [17], The CMB data merely provide a small window into high-energy physics of inflation, while there are no observational constraints for the scales well beyond or much smaller the inflationary scale of about 10 13 GeV (in the case of the Starobinsky inflation).

Why Primordial Black Holes?
Primordial density fluctuations in the early Universe (during or after inflation) may also be responsible for seeding PBH when their amplitude is larger by the factor of 10 7 compared to the CMB amplitude. The very idea of PBH was proposed the long time ago by Zeldovich and Novikov [3] and also by Hawking [4]. After entering the horizon, growing up and then contracting within the gravitational radius, those large primordial fluctuations can form PBH that may survive in the current Universe and provide candidates for (non-particle) Dark Matter (DM) [6,21,22].
More recently, the renewed interest in PBH is mainly related to the possibility that part or the whole of DM may be composed of PBH, which gives a viable alternative to the particle physics explanations of DM as being composed of Weakly Interacting Massive Particles (WIMP), such as neutralino, or Gravitationally Interacting Massive Particles (GIMP), such as axion [23] or gravitino [24]. Given DM in the form of PBH, DM phenomenology should search for DM signals in cosmological data rather than in direct detection on colliders or indirect detection in astroparticle physics. Another reason for the considerable attention on PBH in the recent literature is the observational progress in lensing, cosmic rays, gravitational waves, and CMB radiation, see, e.g., in [7,8] for a review of observational constraints on PBH. In addition, PBH may offer a solution for many astrophysical puzzles, such as, e.g., the existence of SuperMassive Black Holes (SMBH). The SMBH already existed during the first Giga-year of our Universe [25], which gives the powerful argument for their primordial origin [26].
There are several mechanisms that may catalyze the formation of PBH in the early Universe: (i) gravitational instabilities induced by scalar fields [27] in single-field or multi-field inflation, (ii) bubble collisions from first order phase transitions [28][29][30], and (iii) formation of critical topological defects such as cosmic strings [31] and domain walls [32,33]. This paper is devoted to the case (i) only.
On the theoretical side, PBH are considered as a probe of very high energy physics and quantum gravity, "even if they never formed" [34]. Numerous phenomenological scenarios were proposed for PBH formation and, especially, for PBH generation after inflation in the early Universe, under the assumption that PBHs do contribute to DM, see, e.g., in [35][36][37][38][39][40][41] and the references therein. At present, the whole PBH DM appears to allow only two limited windows for the PBH masses, namely, either around 10 −15 or around 10 −12 of the solar mass. Those masses are far less than the masses of black holes whose mergers resulted in the gravitational waves observed by the LIGO detector [42]. Thus, it seems that DM is more likely to be composed of both DM particles and PBH. It is therefore important to study various mechanisms of PBH formation and estimate the PBH fraction of DM both in the single-field and multi-field inflationary models.
Due to the absence of observed non-Gaussianities and iso-curvature perturbations in the current CMB data [12], the single-field models were distinguished and discriminated in the literature both in relation to inflation and PBH formation. Though, as was argued above, the Starobinsky inflationary model is to be favored among all of them at present, its sharp prediction for the tensor-to scalar ratio (r) may be ruled out by future observations. Moreover, the Starobinsky inflation does not allow PBH production after inflation. It is therefore of importance to generalize the Starobinsky model in order to allow arbitrary values of the tensor-to-scalar ratio r and PBH production in the models with single inflaton field. This first problem (any value of r) can be solved by employing the so-called alphaattractor models reviewed in Section 3. The second problem (enhancement of the power spectrum of perturbations by the factor of 10 7 compared to the CMB amplitude, which is needed for PBH formation) can be achieved by modifying the inflaton scalar potential with a nearly inflection point [43][44][45]. However, details of the PBH production after single-field inflation are very much dependent upon a choice of the inflaton potential with an inflection point that requires significant fine-tuning of the model parameters.
As a matter of fact, there are no fundamental reasons for the absence of non-Gaussianities and iso-curvature perturbations, while they just have to be below the observational limits. Relaxing fine-tuning and making PBH production more natural are the basic reasons for considering multi-field models of inflation and PBH. For instance, PBH production may be a generic feature of two-field inflation with a sharp turn of inflationary trajectory [46]. The required growth of primordial fluctuations can be achieved by tachyonic instabilities, like in the waterfall phase of hybrid inflation [47,48]. In addition, inflaton interactions with other particles are of crucial importance for reheating and pre-heating which are beyond the scope of this paper.

Why Supergravity?
Despite clear merits of multi-field models of inflation and PBH formation, there is an obvious problem because multi-field models increase a number of the physical degrees of freedom and possible interactions, which reduce the predictive power of those models. It is, therefore, important to restrict them and propose some general principles. Supersymmetry (SUSY) is a very good choice of the fundamental principle that organizes particles into irreducible supermultiplets and restricts a number of independent coupling constants. In the context of gravity, one needs local supersymmetry, i.e., supergravity that is expected to be important in multi-field inflation by limiting the number of fields involved and severely restricting their interactions.
Supergravity is also a good framework to study the possible theoretical origin of PBH at a more fundamental level than General Relativity (GR). Supergravity is also the first step towards superstring theory and is considered as a viable candidate for a theory of quantum gravity. Because of the constraints imposed by local supersymmetry on possible couplings, a viable description of inflation and PBH in supergravity should lead to significant discrimination of phenomenological models.
In this paper, we adopt the most economical approach by minimizing the relevant number of the physical degrees of freedom needed for inflation and PBH production, by using the Starobinsky model as a guide, in the framework of modified supergravity that can be viewed as a locally supersymmetric extension of the (R + R 2 ) gravity. In practical terms, it implies the supergravitational origin of both inflation and PBH, by using only the (off-shell) N = 1 supergravity multiplet and its locally supersymmetric interactions in four spacetime dimensions, without adding other matter, i.e., with the truly minimal number of the physical degrees of freedom in SUSY. Our modified supergravity models can be reformulated in terms of the standard (Einstein) supergravity coupled to chiral and vector superfields, similarly to the well-known relation between modified F(R) gravity and scalar-tensor gravity (see Section 2.2). Our supergravity approach is thus quite similar in its spirit to the Starobinsky inflation based on gravitational interactions [20,49,50] but allows us to relate inflation with PBH and the DM genesis, in addition.
To the end of this section, brief comments on the general merits and demerits of SUSY and supergravity are in order, in light of no evidence for SUSY particles at the Large Hadron Collider (LHC) so far.
Despite the current absence of experimental confirmation, SUSY remains one of the leading candidates for new physics beyond the Standard Model (SM) of elementary particles. SUSY has to be (spontaneously or softly) broken at energies available for observation. However, the scale of SUSY breaking is unknown and can be much higher than the TeV scale accessible by the LHC.
Supergravity is the field theory with local SUSY that automatically implies the general coordinate invariance. The minimal N = 1 supergravity in four spacetime dimensions is chiral, that is, necessary for particle phenomenology and CP violation. SUSY and supergravity also have other attractive theoretical features, namely, • SUSY unifies bosons and fermions; • supergravity includes GR; • SUSY Grand Unified Theories (super-GUT) lead to the perfect unification of electroweak and strong interactions; • the spectrum of matter-coupled supergravities with spontaneously broken SUSY has the natural DM candidate given by the Lightest SUSY Particle (LSP) provided that R-parity is conserved; • low-energy SUSY helps to stabilize the fundamental scales (the hierarchy problem), such as the electro-weak scale and the GUT scale; • SUSY leads to cancellation of quadratic UV-divergences in quantum field theory; • supergravity is the only way to consistently describe spin-3/2 particles with gravity; and • supergravity arises as the low-energy effective action of superstrings.
In summary, SUSY and supergravity are alive and healthy, while they are not ruled out by observations despite the current absence of any signs of their presence at the LHC and in the sky.

Single-Field Models
In this section, we briefly review the single-field models of inflation and PBH formation. First, we introduce the alpha attractors as the generalizations of the Starobinsky inflation. Next, we give their minimal realizations in supergravity.

Power Spectrum and Generalized Alpha Attractors
To set up calculation tools, here we recall the textbook formulae describing the spectrum of perturbations in terms of the inflaton scalar potential in the single-field models of inflation (see also the original paper [51]).
The spectrum of scalar (density) perturbations is given by for the scales k = a(t k )H(t k ), where G N is Newton's constant, and the subscript k of the potential V and its derivatives (denoted by primes) refers to their corresponding values. In these terms, e-folds are given by N e (k) = ln(k end /k). The scalar spectral index (tilt or slope) is given by Similarly, the spectrum of tensor perturbations (gravitational waves) is given by while its spectral index is Accordingly, the tensor-to-scalar ratio reads As was already mentioned in Section 2.2, the Starobinsky model defines the universality class of single-field inflationary models with the critical parameter α Star. = √ 2/3 that corresponds to the (R + R 2 ) gravity. Replacing it with arbitrary α defines the different class of inflationary models known in the literature as the (E-model) alpha attractors [52] with the potential In modified gravity, this potential corresponds to a more complicated function F(R) that is different from the quadratic one, see Equations (10) and (11). The key feature of alpha attractors is the (leading) value of the tensor-to-scalar ratio (on CMB scale) The alpha attractors thus allow arbitrary suppression of r, whereas they have the same (leading) value of n s ≈ 1 − 2/N e (on CMB scale) as in the Starobinsky model.
The alpha attractors can be generalized even further [52] by considering the inflaton potentials of the form with a monotonically increasing (during slow roll) function f , the parameterα = 3 2 α, and =M −1 Pl . In those models, slow-roll inflation occurs for large positive values of the canonical inflation field ϕ with the approximate scalar potential (=1) where we have introduced the parameters f ∞ = f | ϕ→∞ and f ∞ = ∂ ϕ f ϕ→∞ . The constant in front of the second term in Equation (28) can be adjusted at will by a constant shift of the inflaton field, so that the potential (28) can be simplified to thus establishing the asymptotic equivalence to the alpha attractors on CMB scales. Having obtained significant freedom given by arbitrary function f in the potential (27), one can exploit it for generating new physics after inflation. For instance, engineering a nearly inflection point in the potential can be used for a significant enhancement of the power spectrum of perturbations at smaller (than Hubble) scales towards formation of primordial black holes and related dark matter. As was demonstrated in [53], it can be achieved via Taylor expansion of the function f , while keeping only the first three terms in it as follows: after fine-tuning of the parameters V 0 and c i for i = 1, 2, 3. Among plenty of other possibilities, the same purpose (formation of primordial black holes) can be achieved by using modulated potentials as, e.g., the one also studied in [53] with the potential with the tuned parameters V 0 , A, and f . The common properties of those models are (i) the generated PBH masses in the low mass window of 10 17 ÷ 10 20 g and (ii) slightly smaller values of n s (still compatible with Planck data within 3σ), see in [53] for details. It is, therefore, of interest to investigate whether those features are still present in multi-field models (Section 4).

Single-Field Models of Inflation and PBH in Supergravity
SUSY requires superpartners for each particle and equal numbers of bosonic and fermionic degrees of freedom, as well as supersymmeteric actions. To match observations and connect to the SM, SUSY has to be (spontaneously) broken. The formal technology of local supersymmetry is based on either (i) the superconformal tensor calculus [54] or (ii) the curved superspace [55]. Both approaches are equivalent but technically involved. The N = 1 superspace technology is geometrical and offers manifest SUSY of supersymmetric actions. Being applied to inflation, it means that the bosonic actions (4)-(6) have to be supersymmetrized. Note that the action (6) includes higher derivatives so that its supersymmetrization goes beyond the textbooks [54,55], and it is known as the modified (Starobinsky) supergravity. In addition, inflation has a positive energy (or a positive height of the potential) and thus breaks SUSY. Therefore, there must be a Goldstone spin-1/2 fermion (called goldstino) associated with spontaneous SUSY breaking during inflation.
Both inflaton and goldstino have to belong to (irreducible) supermultiplets or superfields. Usually, both superfields are chosen to be chiral with the maximal spin 1/2 [56,57], which requires complexification of the inflaton and the need of two chiral superfields having four real scalars and, thus, multi-field inflation in general. It is also possible to identify the inflaton chiral superfield with the goldstino chiral superfield, which leads to two-(scalar)-field inflation [58,59].) Slow-roll inflation is usually realized by engineering the scalar potential V in terms of a Kähler potential K and a superpotential W of the chiral superfields as (M Pl = 1) where the primes stand for the derivatives with respect to chiral superfields (we use the same notation for superfields and their leading field components, and suppress indices). Single-field inflation is possible by identifying inflaton with one of the scalars, while suppressing the other scalars during inflation. There is no need to learn supergravity theory for that. There are, however, several problems with the standard approach in supergravity. First, it is the so-called η-problem related to the e K factor in the scalar potential that may be too steep and thus does not have the required flatness during inflation or, equivalently, may lead to a large slow-roll parameter η. Second, there is no fundamental justification for a choice of the inflationary trajectory in the (scalar) multi-field space. Third, potential instabilities of the inflationary trajectory (because of inflaton mixing with other scalars) may easily spoil single-field inflation and fail to achieve the desired number of e-foldings. Though all those problems are solvable, they require careful engineering of the potentials K and W, which implies low predictive power of the standard approach to inflation in supergravity.
There exist, however, a different and truly minimal supergravity framework that is more suitable for embedding single-field inflationary models to supergravity [60][61][62][63], with the only restriction imposed by SUSY on the inflationary potential to be a real function squared. For instance, Ref. [64] has an example of the inflaton potential in the minimal supergravity framework with an inflection point needed for PBH formation. The problem of inflaton complexification in a chiral supermultiplet can be removed by assigning inflaton to a massive vector supermultiplet V that has single physical scalar. The scalar potential of a vector multiplet is given by the D-term instead of the F-term, while any desired values of the CMB observables (n s and r) are possible because the inflaton potential is given by the squared derivative of arbitrary real potential J(V), see below. The manifestly supersymmetric Lagrangian in curved superspace of supergravity is given by where the W α denotes the abelian superfield strength of the vector superfield V, and the abbreviation c.c. stands for the complex conjugated term. The bosonic part of the Lagrangian in Einstein frame (after Weyl rescalings and elimination of the auxiliary fields) reads [60,61] where C = V| is the real (scalar) inflaton given by the leading field component of the superfield V, F mn = ∂ m B n − ∂ n B m , J is a function of C, and g is the coupling constant. The input for model building is given by a real function J (instead of K and W). Ghost freedom requires J (C) > 0. For instance, the scalar potential of the Starobinsky inflationary model is obtained by choosing the J potential as in terms of the canonical inflaton ϕ. In the case of the generalized alpha attractors, one gets the nonlinear differential equations Differentiating the first equation with respect to C, substituting the second equation, and taking square yields for a given function f (tanh). It completes our description of the minimal embedding of any single-field model of inflation and PBH formation into supergravity, as long as its inflaton potential is a real function squared. We conclude this subsection with a comment about the relation between the minimal supergravity above and the standard (non-minimal) procedure in terms of chiral superfields [15].
The master function J(V) can be replaced by the functionJ(He 2VH ), where we have introduced the (charged) Higgs chiral superfield H and have chosen g = 1 for simplicity. The argument of the functionJ and, thus, the functionJ itself are both invariant under the gauge transformations whose gauge parameter Z is also a chiral superfield. The original theory of the massive vector multiplet governed by the master function J is recovered in the supersymmetric gauge H = 1. We can, however, choose another (Wess-Zumino-type) supersymmetric gauge V = V 1 , where V 1 describes the irreducible massless vector gauge multiplet minimally coupled to the dynamical Higgs-type chiral multiplet H. The standard Higgs mechanism appears with the canonical function J = 1 2 He 2VH that corresponds to a linear functionJ [62,63]. This phenomenon is known in the SUSY literature as the super-Higgs effect [55]. In the case under consideration, the supersymmetric U(1) gauge theory in terms of the superfields H and V 1 coupled to supergavity defines the analogue of Higgs inflation that is equivalent to the Starobinsky inflation by our construction, because both models arise in two different gauges of the same theory. The chiral (Higgs) superfield H is charged with respect to U(1), is coupled to the U(1) supersymmetric gauge theory, and is, in fact, the gauge (unphysical) degree of freedom that can be eaten up by the vector gauge multiplet to become massive.

Two-Field Models in Modified Supergravity
In this section, we reverse our strategy: instead of embedding the known viable inflationary models into supergravity (as was the case in the preceding sections for single-field inflation), we use an extension of the Starobinsky (R + R 2 ) gravity model to supergravity in order to generate new viable models of inflation and PBH formation. When using the modified supergravity alone (without SUSY matter described by chiral or vector superfields), more physical scalars (in addition to inflaton) are inevitable and necessarily lead to multi-field inflation. Here, we want to treat it as an advantage over single-field inflation in order to relax fine-tuning and allow more mechanisms for natural PBH generation.
Supergravity is expected to constrain content and interactions of all scalars present during inflation and PBH formation. In this section, we take the minimalistic approach by demanding all fields to belong to the single supergravity multiplet. In particular, all scalar particles (including inflaton) will be superpartners of graviton and gravitino, while their origin and interactions will be entirely supergravitational, just in the spirit of the Starobinsky model. In concrete terms, the (modified) Starobinsky-type supergravity is going to be used for deriving the NLSM metric G ab (φ) and the scalar potential V(φ) in the multi-field models (5), where a, b, . . . = 1, 2 . . . , 4.

Modified (Starobinsky-Type) Supergravity
First, we briefly describe our supergravity setup in curved superspace of the (oldminimal) N = 1 supergravity. We use the standard notation in [55]. Our starting actions are complete in the sense that they are manifestly supersymmetric by construction.
The standard superspace Lagrangian of chiral superfields, Φ i coupled to supergravity is given by (M Pl = 1) where E is the chiral density superfield, R is the chiral curvature superfield, D α ,Dα are the superspace covariant derivatives, D 2 ≡ D α D α , andD 2 ≡DαD˙α.
A (non-holomorphic) Kähler potential K(Φ i ,Φ i ) and a (holomorphic) superpotential W(Φ i ) define the model and uniquely determine its scalar sector in the form (5). A chiral superfield is expanded in terms of its field components as follows: After eliminating the auxiliary fields (F) and going to Einstein frame, the bosonic part of the Lagrangian (39) is given by where we have used the same notation for the superfields and their leading field components, together with The Θ-expansion of the chiral superfields E and R of supergravity in terms of their field components is given by where e ≡ det(e a m ), ψ mn ≡D m ψ n −D n ψ m andD m ψ n ≡ (∂ m + ω m )ψ n . The chiral superfield E is the SUSY extension of e = √ −g, and R is the SUSY extension of the (Ricci) scalar curvature R. The real vector b m and complex scalar X are known in the supergravity literature as the (old-minimal) set of auxiliary fields needed to complete the supergravity multiplet with a closed algebra of SUSY transformations. In the modified (Starobinskytype) supergravity, those "auxiliary" fields become dynamical (or propagating) because of the presence of higher derivatives in the initial Lagrangian (see below). The curved superspace techniques allow us to manifestly supersymmetrize the (R + R 2 ) gravity model as follows [65,66]: or, equivalently, as This (modified) supergravity action is governed by two potentials (N and F) of the single chiral superfield R having the scalar curvature R as its field component at Θ 2 . Therefore, no higher powers of R (beyond the linear and quadratic terms) can appear. The first (D-type) term and the second (F-type) term in Equation (45) are similar to the two terms in Equation (39), but they do not represent matter terms and depend upon the supergravity multiplet only.
The action (45) can be transformed into the standard (Einstein) matter-coupled supergravity action of the type (39) in terms of two chiral matter superfields, whose Kähler potential and the superpotential are expressed in terms of the input potentials N and F [1,49,65,66]. However, after that transformation, the supergravitational origin of the induced Kähler potential and the superpotential are lost.
In this section, we prefer to deal with the action (45) directly. Our modified supergravity setup goes beyond the supergravity textbooks because the standard (Einstein) supergravity actions are the extensions of Einstein-Hilbert term linear in R, whereas our Equations (45) and (46) are more general. They reduce to the the standard Einstein supergravity action in the special case of N = 0 and F = −3R. In other words, Equations (45) or (46) is the most general extension of the (R + R 2 ) gravity. Though an F-type term (except a constant) can be rewritten as the D-type term in Equations (45) and (46), we keep both of them because they represent different modifications of the standard (pure) supergravity.

The Effective Two-Field Models
Let us expand the functions N and F in Taylor series and keep a few leading terms as a probe of our modified supergravity as follows: where we have introduced three free parameters-ζ, γ, and δ-together with their normalization by the Starobinsky mass parameter M. Let us also ignore the vector field b m and the angular part of the R| = X field for simplicity by setting where σ is the real scalar field. It is straightforward (but tedious) to obtain the scalar part of the Lagrangian in Einstein frame. It takes the form of Equation (5) indeed [67], where the functions A, B, U are given by The model (50) has two scalars: scalaron ϕ and another scalar σ, while both have the supergravitational origin: the scalaron φ may be considered as spin-0 part of metric and σ as its superpartner. Moreover, their kinetic terms and the scalar potential are determined by Equations (50) and (51). One may worry about σ to become a ghost field due to the negative signs at some terms in the second Equation (51). However, the infinite wall in the scalar potential prevents σ from obtaining the values leading to the wrong sign of its kinetic term.
It is also straightforward to derive the equations of motion in our two-field model in the FLRW universe, when keeping only time dependence. We find [67] where we have added the Friedmann equation at the end. Our results for double inflation and PBH formation in the two special cases: the one with δ = 0 (called the γ model) and the one with γ = 0 (called the δ model) are given below, in accordance with the work in [67]. Every term in Equations (47) and (48) is well motivated, as follows: (i) the first term in F comes from the pure (standard) supergravity, (ii) the first term in N describes the R 2 supergravity, (iii) the second ζ-term in N is required for stabilization of the R 2 supergravity and of its scalar potential, and (iv) the remaining γ-and δ-terms are introduced for PBH formation.

The γ Models
As a representative of the γ models (δ = 0), let us choose the case with the parameters γ = 1 and ζ = −1.7774, see Figure 1. The scalar potential (for ϕ 1) has two valleys at σ = 0 and a single Minkowski minimum at σ = ϕ = 0. The first Slow-Roll (SR) inflation is possible along each of the valleys. The two valleys merge into the Minkowski minimum by passing through the near-inflection points leading to the intermediate Ultra-Slow-Roll (USR) stage, where Hubble friction dominates over the potential slope, followed by the second (SR) stage of inflation. The SR approximation must be violated for PBH production [45]. In the USR regime, the scalar field(s) actually roll down the potential faster than in the SR regime [68]. We plot our numerical solutions to the equations of motion in Figure 2. The total number of e-folds is set to ∆N = 60, and the end of the first stage of inflation is defined by the time when the SR parameter η first reaches 1 (see Figure 2e). The USR period USR SR in Figure 2e leads to the significant enhancement of the scalar power spectrum (see below). Inflation ends when = 1. The first stage lasts ∆N 1 ≈ 50 e-folds, whereas the second stage lasts for ∆N 2 ≈ 10. In our figures, the first stage of inflation is represented by the blue shaded region, whereas the second stage is marked by the green shaded region. The length of the second stage is controlled by the parameter ζ for a given γ. The CMB observables in this model (as the best fit) are given by n s ≈ 0.9600 and r max ≈ 0.004 .
The value of the spectral tilt n s is in tension with the Planck observations [12], whereas the value of the tensor-to-scalar ratio r is essentially the same as that in the Starobinsky model.
The power spectrum of curvature perturbations can be numerically computed at fixed ∆N 2 by using the standard transport method [69,70] with the Mathematica package [71]. We compute the spectrum around the pivot scale k * that leaves the horizon at the end of the first stage (we call this scale k ∆N 2 ). The inflaton mass is chosen to be 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 (we call it k 60 ). Our results for various values of γ are shown in Figure 3. The values of the parameters are collected in Table 1, where ζ is tuned to satisfy ∆N 2 ≈ 10. The changes in n s and r max are negligible for those parameters.  Figure 3. The power spectrum P ζ near the pivot scale k * = k ∆N 2 at ∆N 2 = 10 for some values of γ. Table 1. The parameters leading to the power spectrum in Figure 3 with ∆N 2 ≈ 10. The enhancement of primordial curvature perturbations needed for PBH formation is usually assumed to be about 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 the CMB scale. However, given a broad peak, the enhancement in the power spectrum may be one order of the magnitude less [36]. According to our numerical estimates (see Figure 3), 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.
The power spectrum also changes with ∆N 2 . To demonstrate that dependence, we consider the power spectrum at γ = 0.1 and γ = 1 for various values of duration of the USR stage, ∆N 2 = 10, 20, 23, 26 for each γ. The results are collected in Figure 4. In order to achieve the enhancement of the order of 10 6 , we must require ∆N 2 30. However, this pushes the spectral index outside of the 3σ (lower) limit of n s ≈ 0.946 again.
In Table 2, we collect the 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. The mass of PBH created as a result of the primordial power spectrum enhancement can be estimated from the peak data as follows [35]: where t peak is the time when the wavenumber corresponding to the power spectrum peak (k peak ) exits the horizon, and t 60 is the time when k 60 exits the horizon.
Our values of M PBH for some values of ∆N 2 from Equation (57) 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. On the one hand, the PBH with masses smaller than ∼ 10 16 g would have already evaporated until now via Hawking radiation. Therefore, 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. Thus, the PBH masses are restricted by O(10 16 g) < M PBH < O(10 19 g) even before considering the observational constraints on them.   Regarding the constraints on γ, the power spectrum in Figure 4 tells us for ∆N 2 > 20 that it is sufficient to have γ O(1) in order to produce the required enhancement in the spectrum. Therefore, the dependence of our results upon the value of γ is weak.
The PBH density fraction in DM can be estimated by using the standard (Press-Schechter) formalism [72]. 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., in [73,74] and the references therein): 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. It is usually assumed that δ c ≈ 1/3 [5], though it may be different. Then, the PBH-to-DM density fraction is estimated as [73,74] In order to numerically evaluate the functions (58) and (59), we have to normalize the values of k in terms of the observable scales today. We choose the scale k 60 to represent the largest currently observable scale around 10 −4 Mpc −1 . Then, our numerical evaluation shows that we need the δ c parameter smaller than 1/3, say, close to δ c = 0.275.

The δ Models
Having found the γ models to be in tension with the observed value of the spectral tilt n s (unless the total number of e-folds is significantly reduced to less than 60), we study the δ models in this subsection, when γ = 0 and δ = 0 in Equations (47) and (48). It breaks the R-symmetry and the reflection symmetry σ → −σ of the potential, see Figure 5. For any non-zero δ, there is a value of ζ leading to an inflection point. In contrast to the γ models, here we have a single valley for large positive ϕ and σ = 0, where the δ models reduce to the single-field Starobinsky model. When approaching ϕ = 0, the inflationary trajectory passes 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, where ζ is chosen to get ∆N 2 = 10. After numerically solving the field equations, we find the time dependence of ϕ, σ,H, N, , and η in Figure 6. The near-inflection point divides inflation into two stages with ∆N 1 = 50 (SR) and ∆N 2 = 10 (USR). We set the initial velocities to zero, with ϕ(0) = 6 and σ(0) = 0.05. Similarly to the γ models, the inflationary trajectory is (locally) stable against variations of the initial conditions.
Having fixed ∆N 2 = 10, we find that the spectrum has a non-trivial dependence on δ, see Figure 7. In the left plot, δ is changed 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.
The dependence upon ∆N 2 for the highest power spectrum peaks with δ = 0.1 and δ = 0.6 is displayed in Figure 8. As expected, the enhancement becomes larger with increasing ∆N 2 .
Regarding PBH masses and their density fraction, 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 (59) to be close to one and the corresponding density threshold to not deviate far away from the region 1/3 ≤ δ c ≤ 2/3, we find that the values δ = 0.094 and δ = 0.58 are suitable for efficient generation of PBH with their masses around 10 19 g, while avoiding their overproduction ( f 1). We estimate the generated PBH masses by using Equation (57) and summarize our results in Table 4. To get 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).  Let us now compare the specific models studied in the preceding subsections in order to tune their parameters to the best fit with the CMB observations, and then compare those models with the available observational constrains on PBH and DM [67,75]. We also want to estimate the PBH-to-DM density fractions in the γ-models and in the δ-models.
For those purposes, we pick up one γ-model (Case I) and two δ-models (Cases II and III, respectively) with different shapes of the power spectrum. The need of two δmodels is motivated by the existence of two suitable parameter regions, where δ 0.1 and δ 0.6 yield truly different shapes of the power spectrum (broad and narrow, respectively). The parameter sets of those three models are given in Table 5, and the corresponding (numerically computed) power spectra P ζ and PBH density fractions f (M) are shown in Figure 9. We use the normalization of the wavenumber k exit = 0.05 Mpc −1 , where k exit is the scale that leaves the horizon around 54 e-folds before the end of inflation. The parameter ζ is fixed by choosing the value of ∆N 2 at given γ and δ. In the cases I-III, we find ζ as −2.374, 0.032, and 0.102, respectively. Table 5. The parameters corresponding to the PBH density fraction in Figure 9. The n s and r are computed with the total ∆N = 54 e-folds before the end of inflation.   Table 5 with PBH as the DM, on the left plot (a). The k * represents the end of SR and the beginning of USR. The corresponding PBH density fractions are given on the right plot (b). The background observational constraints on PBH are taken from Refs. [9,76]. 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.
According to Table 5, the spectral tilt n s in case I is in tension by 3σ against the CMB data [12], whereas in cases II and III, the value of n s is within the current 3σ constraints. The PBH fraction in case II of Figure 9 implies that those δ-models are more flexible to accommodate slightly larger values of n s . It happens because the PBH fraction in 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 further the PBH masses. In turn, it will decrease ∆N 2 and, consequently, increase n s .

Gravitational Waves Induced by PBH Formation
One of the most visible opportunities to check the proposed models of inflation and PBH production by observations is via detection of the stochastic Gravitational Waves (GW) induced by the PBH formation. We follow the work in [75] here.
The present-day GW density function Ω GW is given by [77,78] where the constant c g ≈ 0.4 in the case of the Standard Model (SM), and c g ≈ 0.3 in the case of the Minimal Supersymmetric Standard Model (MSSM). The present-day value of the radiation density Ω r is equal to h 2 Ω r ≈ 2.47 × 10 −5 , according to measurements of CMB temperature [79]. Here, h is the reduced (present-day) Hubble parameter that we take as h = 0.67 (ignoring the Hubble tension). The variables x, y are related to the integration variables s, d as while the functions I c and I s of x(s, d) and y(s, d) are given by [77,78] I c = −4 ∞ 0 dη sin η 2T(xη)T(xη) + T(xη) + xηT (xη) T(yη) + yηT (yη) and where the T-function is defined by in terms of the conformal time η.
The integrations in I c and I s can be performed analytically, and the results are [77] where θ is the Heaviside step function. Using the definitions above, the GW density can be numerically computed from a given power spectrum. It suffices to consider the power spectra in the cases of Table 5 where PBH are part of DM because the cases with f tot = 1 have the quite similar power spectra, though with slightly larger peaks. By using the power spectra of Figure 9 (on the left side) we plot the density Ω GW (k) in terms of frequency k = 2π f in Figure 10 together with the expected sensitivity curves for several space-based GW experiments planned in the near future. We have used the power-law integrated curves [80] and have applied them to the LISA noise model [81,82]. As the alternative, the peak-integrated curves can be used [83]). To draw the sensitivity curves, we have used the parameters and the noise models for TianQin [84], Taiji [85], and DECIGO [86].
It follows from Figure 10 that the upcoming space-based GW experiments are expected to be sensitive enough to detect the stochastic GW background predicted by a large class of our two-field inflationary models where PBH may account for a significant fraction (or all) of DM. Figure 10 shows that our supergravity models produce GW peaking in the frequency range 10 −3 ÷ 10 −1 Hz expected to be accessible by LISA, TianQin, Taiji, and DECIGO gravitational interferometers.

Conclusions and Discussion
The proposed extensions of the Starobinsky inflationary model in the modified supergravity can produce PBH that could account for part (or whole) of Cold Dark Matter. The PBH generation can be efficiently catalyzed by primordial perturbations sourced by the Starobinsky scalaron coupled to another scalar of the (super)gravitational origin.
Our conclusions are mainly based on the theoretical considerations (by putting forward the two fundamental principles, namely, modified gravity and supergravity) before comparing them to cosmological observations, in the top-down approach. We started from the Starobinsky inflationary model serving as the theoretical tool, and argued about its universality (as the most basic model) for slow roll inflation. Then, the alpha-attractor models were reviewed as the natural extensions of the Starobinsky inflation in the different universality classes with different (arbitrary) tensor-to-scalar ratios.
All those single-field models of inflation and PBH were minimally embedded into the particular version of supergravity with inflaton in a massive vector multiplet after some fine-tuning of the parameters needed for efficient PBH formation. The formal theoretical advantages of supergravity are obvious: (a) it may be further embedded into a fundamental theory of quantum gravity such as superstrings; (b) the modified supergravity offers the supergravitational origin of inflaton and its potential, being (classically) equivalent to the standard (Einstein) matter-coupled supergravity; and (c) extra fields and extra couplings are severely restricted by SUSY, while the relation between them is not broken by renormalization due to SUSY preservation against quantum (gravity) corrections. Therefore, supergravity has higher predictive power for phenomenology of the early universe cosmology. In the reverse way, inflation and PBH DM can also be considered as the probes of supergravity for its role as a more fundamental approach to cosmology.
All that becomes even more important for multi-field models of inflation, PBH, and DM, where SUSY and supergravity can serve as the guiding fundamental principle for discrimination of phenomenological models. We find that it is the case for the effective two-field models derived from the modified (Starobinsky) supergravity by starting from the manifestly supersymmetric Lagrangian (46). After (Taylor) expanding its potentials N and F in powers of the scalar curvature superfield R, and keeping only the few leading terms, we arrived at our models defined by Equations (47) and (48), whose scalar part is described by Equations (50) and (51). Unlike the other two-field models in the literature, where Starobinsky's scalaron is paired either with dilaton [87,88] or a Higgs field [89][90][91], the modified supergravity predicts another scalar having the supergravitational origin and belonging to the single supergravity multipet together with graviton and gravitino.
The power spectra were numerically computed, and the PBH masses and their density fractions in DM were estimated in three special cases of those two-field models. We found that each model can simultaneously describe both the viable (Starobinsky-type) inflation and the efficient PBH production after inflation, with remarkably less fine-tuning of the parameters. Apparently, the δ models are preferable against the γ models from the viewpoint of the precision measurements on the CMB spectral index n s .
We also found that a significant part of DM is realizable in the form of PBH in our two-field models, with the PBH masses being in the range between 10 16 g and 10 20 g. The whole DM as the PBH is also possible, though after some more fine-tuning of the model parameters. Note that the PBH masses found are beyond the lower bound provided by Hawking radiation, but are much lower than the Solar mass and the black hole masses discovered by the LIGO experiment [92].
Unlike single-field inflationary models demanding an inflection point of the potential for enhancement of curvature perturbations, multi-field inflationary models offer other natural mechanisms of the power spectrum enhancement, which may be different from the double inflationary scenario adopted in this paper. For instance, it could be via processes of isocurvature mode amplification via sharp turns of the inflationary trajectory [38,88] or via first order phase transitions and related bubble collisions [28]. All that implies PBH production may be a generic feature of multi-field inflation, though more studies are needed to check it in the context of supergravity.
As is well known, PBH formation necessarily leads to GW because large scalar overdensities act as a source of stochastic GW background. Frequencies of those GW can be related to the expected PBH masses and the duration of the second stage of inflation [78]. These GW may be detected in the future ground-based experiments, such as the Einstein telescope [93] and the global network of GW interferometers, including advanced LIGO, Virgo, and KAGRA [92], as well as by the space-based GW interferometers such as LISA [81], TAIJI (old ALIA) [85], TianQin [84], and DECIGO [86].
Our results support the proposal that PBH may account for a large part or the whole DM, while those PBH and DM may have their origin in supergravity. Our supergravity models also predict the GW stochastic background radiation that is sensitive to the inflationary parameters and the PBH mass spectrum. Fine-tuning in our models amounts to fixing the parameter M ∼ 10 −5 M Pl as the (Starobinsky) scalaron mass and the dimensionless parameter ζ for the desired duration of the USR. The obtained PBH mass spectra are compatible with all astrophysical and cosmological constrains, while induced GW signals can be detected by the next space-based or land-based gravitational interferometers. The NANOGrav Collaboration data [94] hint to the PBH as DM [95][96][97], in agreement with our results in Figure 9.
Supergravity is often regarded as an extension of gravity at super-high energy scales. We find that the new scalars of modified supergravity can play the active 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, so that our models have the predictive power to be falsified in future experiments. It gives us the reason to believe that indirect footprints of SUSY may be detected from GW physics rather than from high-energy particle colliders.