Primordial black hole dark matter in dilaton-extended two-field Starobinsky inflation

We investigate the production of primordial black holes and their contribution to the presently observed dark matter in a dilaton two-field extension of Starobinsky's quadratic $f(R)$ model of inflation. The model features a multi-field amplification mechanism which leads to the generation of a sharp peak in the inflationary power spectrum at small wavelengths responsible for the production of primordial black holes. This mechanism is significantly different from single-field models and requires a stochastic treatment during an intermediate phase of the inflationary dynamics. We find that the model leads to a successful phase of effective single-field Starobinsky inflation for wavelengths probed by the cosmic microwave background radiation and explains the observed cold dark matter content in the Universe by the formation of primordial black holes.

We investigate the production of primordial black holes and their contribution to the presently observed dark matter in a dilaton two-field extension of Starobinsky's quadratic f (R) model of inflation. The model features a multi-field amplification mechanism which leads to the generation of a sharp peak in the inflationary power spectrum at small wavelengths responsible for the production of primordial black holes. This mechanism is significantly different from single-field models and requires a stochastic treatment during an intermediate phase of the inflationary dynamics. We find that the model leads to a successful phase of effective single-field Starobinsky inflation for wavelengths probed by the cosmic microwave background radiation and explains the observed cold dark matter content in the Universe by the formation of primordial black holes.

I. INTRODUCTION
Primordial Black Holes (PBHs) could provide an explanation of the origin of Cold Dark Matter (CDM) without assuming new particles, generate the seeds of Large Scale Structure (LSS), and probe very high energy physics including quantum gravity, see e.g. [1][2][3][4][5] for a review, cosmological and astrophysical constraints on PBHs, and references therein to original papers.
The phenomenologically most relevant PBH formation scenario is due to large density fluctuations generated during inflation. Aside from the relevance of PBHs in explaining (part of) the currently observed CDM, they offer a unique opportunity to probe the power spectrum of perturbations at smaller wavelengths. Therefore, they provide constraints complementary to those from the observations of the Cosmic Microwave Background (CMB) at large wavelengths and have the potential to limit the number of phenomenologically viable inflationary models. Among the theoretically best motivated models, whose predictions for the inflationary spectral observables are in perfect agreement with recent PLANCK data [6], are Starobinsky's quadratic model of f (R) gravity [7] and the model of Higgs inflation [8], see [9] for a recent review. A combination of the two individual models leads * anirudh.gundhi@phd.units.it; ketov@tmu.ac.jp; christian.steinwachs@physik.uni-freiburg.de to the unified scalaron-Higgs two-field model of inflation considered in [10][11][12][13][14][15][16]. Another multi-field extension of Higgs inflation and Starobinsky's R + R 2 model is based on a coupling to a dilaton field considered in [17,18].
In this paper we investigate the formation of PBHs in a two-field dilaton extension of Starobinsky's model. Although the structure is similar to the scalaron-Higgs model [13][14][15][16], regarding the coupling of the dilaton field, we consider the operators R and R 2 on equal footing. The non-minimal coupling to the linear Einstein-Hilbert term leads to an effective dilaton-dependent Planck mass, while the coupling to the quadratic R 2 term leads to an effective dilaton-dependent scalaron mass. The latter is a key element in the successful realization of the multi-field amplification mechanism that leads to an enhancement of the power spectrum of scalar perturbations crucial for the formation of PBHs.
Our paper is organized as follows. We formulate our model in the Jordan frame and perform the transition to the two-field scalar-tensor formulation in the Einstein frame in Sect. II. We summarize the formalism required for the covariant description of the background dynamics and the evolution of the perturbations in Sect. III. The arXiv:2011.05999v2 [hep-th] 20 Nov 2020 properties of the Einstein frame two-field potential landscape are discussed in Sect. IV. In Sect. V, we explain how inflation proceeds in three subsequent stages involving an intermediate stochastic phase connecting two effective single-field phases of slow-roll inflation. Sect. VI is devoted to the details of the peak formation mechanism, while we present numerical results of the inflationary dynamics and observables in Sect. VII. In Sect. VIII, we summarize the PBH formation mechanism and the formalism required for the calculation of the PBH mass distribution and the total PBH fraction of the currently observed CDM. Finally, in Sect. IX, we present our numerical results for the PBH mass distribution in the three mass windows which permit a significant PBH fraction of CDM. We summarize our main results and discuss various future applications of our model in Sect. X. A derivation of the PBH mass distribution function is presented in Appendix A. An approximate analytical estimate for the power spectrum peak amplitude required for a significant PBH production is provided in Appendix B.

II. THE MODEL
Starobinsky's model is defined by the action [7] S Star [g] = M 2 The higher derivatives entering (1) via the marginal R 2 operator lead to an additional scalar propagating degree of freedom, the scalaron [7,32], which becomes manifest in the scalar-tensor representation of (1). 1 The action (1) is characterized by two dimensional parameters: the reduced Planck mass M P = 1/ √ 8πG N ≈ 2.4 × 10 18 GeV, with Newton's constant G N and the scalaron mass m 0 . In the regime of large curvatures R/m 2 0 1, the marginal R 2 operator in (1) dominates and an inflationary quasi de Sitter phase and an almost scale invariant spectrum of perturbations is naturally realized due to the asymptotic scale invariance. In contrast, in the regime of small curvatures R/m 2 0 1, the linear Einstein-Hilbert term in (1) dominates and naturally realizes a graceful exit from inflation. The particular degeneracy structure of f (R) theories ensures that the scalaron is neither a tachyon nor a ghost, provided M 2 P and m 2 0 are both positive [33,34]. Starobinsky's geometric single-field model of inflation (1) is in excellent agreement with recent Planck data [6]. The main inflationary predictions are the power spectra of tensor perturbations h ij and scalar curvature perturbation R. Due to their weak logarithmic dependence on 1 Curved spacetime is described by a four-dimensional pseudo-Riemannian manifold with local coordinates µ, ν, . . . = 0, 1, 2, 3, metric field gµν (x), and metric-compatible connection ∇µ. Our sign conventions are sig(g) = (−, +, +, +, ), R ρ µσν = ∂σΓ ρ µν − ∂ν Γ ρ µσ + . . ., and Rµν = R ρ µρν .
the wave number k, they are parametrized (neglecting the higher order k-dependent terms) by the power laws The power spectra (2) are fully characterized by the two amplitudes A h and A R , the two spectral indices n h and n R , and the pivot scale k * chosen within the CMB window of scales accessible to PLANCK [35], The moment when k * first crosses the Hubble horizon is chosen to correspond to 50 ≤ N ≤ 60 efolds. Since no primordial gravitational waves have been measured yet, it is more convenient to express the amplitude A h in terms of the tensor-to-scalar ratio r = A h /A R . Due to the single-field consistency condition r = −8 n h , the parametrization of the two spectra (2) is effectively determined only by the three parameters A R , n R and r. PLANCK data [6] constrains the values for A R and n R at k * = 0.05 Mpc −1 , and for r at k * = 0.002 Mpc −1 , Starobinsky's model (1) predicts a scalar amplitude The normalization condition (4) fixes the only free parameter, the scalaron mass as Since the spectral observables are independent of the model parameter, they provide predictions of (1), The numerical values in (8) and (9) have been obtained for N = 60. Even though Starobinsky's inflationary model (1) has maximal predictive power, it does not exclude interactions with other fields. In particular, in string theory inspired cosmological models, as well as in modified supergravity, the coupling of gravity to a dilaton field ϕ arises naturally, see e.g. [36][37][38][39]. A generic coupling of the action (1) to a canonically normalized dilaton field ϕ with generic potential V (ϕ) can be written in the form with the function Compared to (1), the dimensional parameters in (11) have been promoted to generic functions of ϕ. The nonminimal coupling function U (ϕ) replaces the constant Planck mass M P , while the the mass function M (ϕ) replaces the scalaron mass m 0 . Following the general derivation in [16], the Jordan frame (JF) action (10) can be written as a two-field model in the classically equivalent 2 Einstein frame (EF) by performing the non-linear field redefinitions In terms of the field covariant formulation, the EF twofield action resembles a non-linear sigma-model [16], The scalars Φ I (x) are the local coordinates of the scalar field space manifold with the target metric G IJ , In terms of the parametrizations the scalar two-field potentialŴ (Φ) readŝ In order to specify a concrete dilaton extension of Starobinsky's geometric model, we assume the following explicit form of the functions U (ϕ), m 2 (ϕ), and V (ϕ), The choices (18)- (20) are motivated by several considerations: First, we demand that Starobinsky's model (1) is recovered in the limit ϕ → 0, justifying the presence of the field independent constants M P and m 2 0 in (18) and (20) and the absence of a field independent contribution (a cosmological constant Λ) to (20). Second, we assume an additional internal Z 2 symmetry ϕ → −ϕ leaving only ϕ-even operators. Third, we assume that the invariance under the global scale transformations g µν → α −2 g µν and ϕ → αϕ with constant parameter α is asymptotically realized for large field values ϕ/M P 1, justifying the absence of higher order ϕ monomials in (18)- (20). For (18)- (20),Ŵ defined in (17) reduces tô In addition to the two mass parameters M P and m 0 , present in the Starobinsky model (1), the EF two-field potential (21) is characterized by the three dimensionless parameters ξ, ζ, and λ. We do not include a dilaton mass term m 2 D ϕ 2 in (20), which in view of the EF potential (21), can be safely neglected as long as m 2 D ζM 2 P or m 2 D ξm 2 0 . Moreover, since in Sect. IX we find that a significant PBH production requires ξ 1, a dilaton mass term can be neglected as long as m 2 D m 2 0 . It would be interesting to study the impact of a large dilaton mass m 2 D m 2 0 , which, however, goes beyond the scope of our present work.

III. COVARIANT MULTI-FIELD FORMALISM AND INFLATIONARY OBSERVABLES
Following the general treatment in [16], we use the covariant multi-field formalism 3 to formulate the inflationary dynamics of the background and perturbations.

A. Background dynamics
The homogeneous and isotropic background dynamics of the metric is determined by the flat Friedmann-Lemaître-Robertson-Walker (FLRW) line element Here, t is the cosmic Friedmann time, a(t) is the scale factor, i, j, . . . = 1, 2, 3 are spatial indices and δ ij = diag(1, 1, 1) is the spatial metric. The Friedmann equations and the Klein-Gordon equations for the homogeneous scalar fields Φ I (t) read Here the dots denote the derivatives with respect to the cosmic time t. The Hubble parameter H(t) and the covariant time derivative D t are defined as The Christoffel connection Γ I JK is defined with respect to the field space metric (14). The unit vector tangential to the inflationary trajectory readŝ The two-field background dynamics is decomposed into a direction alongσ I and a direction along the unit vector s I orthogonal toσ I , The unit vectorŝ I is proportional to the acceleration vector ω I which defines the turn-rate ω, Projecting (23)-(25) alongσ I andŝ I leads to the set of background equations The derivatives ofŴ alongσ I andŝ I are defined bŷ

B. Perturbations
The perturbed FLRW line element reads Here E ij := ψδ ij + E ,ij and the scalar metric perturbations A(t, x), B(t, x), ψ(t, x), and E(t, x) combine with the perturbation of the scalar fields δΦ I (t, x). Instead of δΦ I (t, x), we work with the gauge-invariant multi-field Mukhanov-Sasaki variables [45][46][47], The equation for the Fourier modes of the perturbation δΦ I g (t, k) is found to be [47][48][49], Following the conventions introduced in [16], Ω I J and the effective mass tensor M I J are defined by The effective mass tensor M IJ includes the Riemannian curvature tensor R IJKL associated with the curved scalar field space manifold, as well as the curvature of the multifield potentialŴ . The tensor modes h(t, k) (suppressing tensor indices) satisfy the simple mode equation, Projecting (34) alongσ I andŝ I defines the adiabatic and isocurvature perturbations Inserting δΦ I g = Q σσ I + Q sŝ I into (35), the dynamical equations for the Fourier modes Q σ (t, k), Q s (t, k) and h(t, k) in the large wavelength limit k aH read h + 3Hḣ = 0.
The projections of (36) and (37) include the additional effective contributions of the turn rate, The operator f (d/dt) in (40), acting on the product ωQ s , is defined as [16] f Equation (40) shows that the adiabatic mode Q σ is sourced by the product ωQ s . The turn rate ω is determined by the background dynamics (31), while the isocurvature mode Q s is obtained by solving the homogeneous equation (41). Only if the combination of ω and Q s is sufficiently large, Q σ is sourced by the isocurvature mode, leading to the amplification of the adiabatic power spectrum. This "isocurvature pumping" amplification mechanism of Q σ was already described in [16] and, as we describe in detail in Sect. VI, is crucial for the production of PBHs and the numerical results obtained in Sect. IX. As in single-field models of inflation, the deviation from de Sitter spaceḢ = 0 is quantified by the Hubble slowroll parameters The power spectra of the scalar perturbations (46) and the tensor perturbation h(t, k) read Note that a simple power law ansatz for the power spectra (47)- (49), as in (2), is no longer adequate when the power spectra feature peaks at small wavelengths.

IV. TWO-FIELD POTENTIAL LANDSCAPE
As is shown in Fig. 1, the landscape of the two-field potential (21) is dominated by three valleys separated by two hills symmetrically located around ϕ = 0 . The action (10) with the potential (21) can also be viewed as an extension of the scalaron-Higgs potential having one additional parameter ζ.
The location of the three valleys and the two hills of the potential (21) are determined by the five roots of the valley equationŴ ,ϕ = 0.
The condition (50) follows from the background dynamics of ϕ andχ, The equations (51) and (52) in turn follow from (25) withχ and ϕ being functions of the number of efolds N defined by dN = −Hdt, i.e. we count N backwards with N = 0 at the end of inflation. Primes denote derivatives with respect to N . For fixedχ, a necessary condition to reach a stationary point in phase space (ϕ, ϕ ) isŴ , ϕ |χ = 0. Asχ takes different values during the background evolution, the classical trajectory is obtained by solving (50) for ϕ(χ). Equation (50) has five solutions The solutions (53) correspond to the central valley at ϕ 0 , the two outer valleys at ϕ ± v , and the two hills at ϕ ± h . At the onset of the inflationary dynamics the initial valueχ i must be sufficiently largeχ i /M P 1, in order to guarantee that inflation lasts at least N ≈ 60 efolds. Inflation ends in one of the outer valleys ϕ ± v close to the global minimum at (χ, ϕ) = (0, 0) when ε H (ϕ ± v ,χ f ) = 1. For inflationary background trajectories, which run along the ϕ 0 valley, there is a field valueχ c at which the local ϕ 0 minimum turns into an unstable maximum. The critical pointχ c is determined by the condition The solution of (54) only depends on m 0 /M P and ξ/ζ, For a fixed m 2 0 and a suitable ratio ξ/ζ, the value ofχ c can be made sufficiently small, such that all CMB modes (3) cross the horizon before the inflationary trajectory crosses the critical point (55). Depending on the parameters of the model, two qualitatively different scenarios are possible. They are shown in the top and bottom rows of Fig. 1, respectively. In Scenario I, the two outer ϕ ± v valleys merge with the central ϕ 0 valley before the critical point is reachedχ >χ c and the resulting landscape is that of a single global attractor at ϕ = 0 (top left plot in Fig. 1). Atχ =χ c , when the second derivative W, ϕϕ turns negative, the local minimum along the ϕ direction turns into an unstable local maximum and the two valleys symmetrically located around ϕ 0 emerge again (top right plot in Fig. 1).
In Scenario II, the two ϕ ± v valleys always run parallel to the central ϕ 0 valley and at no stage merge with it into a single global ϕ 0 attractor. Nevertheless, for initial conditions ϕ i /M P 1, atχ i χ c , the local minimum at ϕ = 0 keeps the inflationary trajectory trapped along the ϕ 0 solution untilχ reachesχ c . As in Scenario I, after the stable local minimum in the ϕ direction turns into an unstable local maximum, the trajectory falls into one of the two lower-lying adjacent ϕ ± v valleys. The parameter combination, which distinguishes between the two scenarios, is obtained by observing that in Scenario I the ϕ ± v valleys merge with ϕ ± h forχ >χ c and re-emerge forχ <χ c . In Scenario II, the valleys ϕ 0 and ϕ ± v run parallel to each other for the entire inflationary dynamics.
As is shown in Fig. 2, mathematically this observation is reflected by the fact that in Scenario I, during a certain period before reachingχ c , the valley solutions ϕ ± v and ϕ ± h become complex. In contrast, in Scenario II, the solutions ϕ ± v are real for all values ofχ. Only the two maxima ϕ ± h , which earlier separated the ϕ 0 valley from the ϕ ± v valleys, disappear (or become complex) for χ <χ c . Demanding that the ϕ ± v valley solutions are real for all values ofχ leads to the constraint λ ≤ 6ξ 2 m 2 0 /M 2 P . Hence, Scenario I and Scenario II are obtained for the parameter combinations As shown in the top right plot in Fig. 1, the merger of the ϕ ± v valleys with the ϕ ± h hills in Scenario I implies that the potential landscape reduces to a broad global attractor around ϕ 0 for someχ >χ c . In general this would make the predictions of this scenario independent of the initial field values ϕ i 4 . In contrast, in Scenario II the background trajectory always remains in one of the two ϕ ± v valleys over the course of the entire inflationary dynamics for initial values ϕ i = ϕ ± v (χ i ). In both scenarios, the dynamics in the vicinity ofχ c is important due to a rapid growth of isocurvature modes which, in combination with a non-zero turn rate, leads to a strong amplification of the adiabatic modes. This amplification mechanism is described in more detail in Sect. VI. As explained in Sect. V B, it requires a careful stochastic treatment including diffusive quantum effects.
Finally, let us compare the potential landscape (21) to that in the model of scalaron-Higgs inflation [16], which is recovered from (21) in the limit ζ → 0. In this limit, instead of the central valley, the potential features a central hill at ϕ = 0, where quantum diffusive effects almost immediately push any background trajectory running along the ϕ = 0 line into one of the outer valleys at the very onset of inflation. This prevents multi-field effects to have any observable consequences, see [16] for a detailed discussion.

V. THREE STAGES OF THE INFLATIONARY BACKGROUND DYNAMICS
The dynamics is naturally divided into three stages: Stage 1: Effective single-field slow-roll Starobinsky inflation and CMB observables at large wavelengths.
Stage 2: Stochastic transition regime with dominant quantum diffusion and tachyonic isocuravture mass.
Stage 3: PBH peak formation at small wavelengths and effective singe-field slow-roll dynamics in a ϕ ± v valley.

A. Stage 1: Starobinsky inflation and CMB
During the first stageχ i ≥χ >χ c , the inflationary dynamics effectively reduces to a single-field slow-roll phase along the ϕ 0 valley. Forχ >χ c , due to the large and positive value of the isocurvature mass (directly related to the curvature of the potential in the ϕ direction), any deviation of the background trajectory from ϕ = 0 and any growth of the δϕ perturbation will be immediately suppressed.
Along the ϕ 0 valley, the two-field potential (21) reduces to the EF potential of the Starobinsky model (1), Thus, the inflationary predictions are identical to those of the Starobinsky model (1) for wavelengths that can be probed by CMB measurements, provided the scalaron mass is fixed to the value in (8), and the ratio ξ/ζ in (55) is chosen such thatχ c is sufficiently small to ensure that all CMB modes (3) cross the horizon during Stage 1. In particular, m 0 is no longer a free parameter but is fixed at this stage, in order to ensure consistency with CMB measurements.

B. Stage 2: Stochastic dynamics and tachyonic isocuravture mass
The formalism of linear perturbation theory, in which the field ϕ(N, x) is decomposed into a homogeneous background fieldφ(N ) and a perturbation δϕ(N, x), can be safely applied to situations in which δϕ φ. Any linear perturbation δϕ around the classical solution ϕ 0 is strongly suppressed by a large and positive curvature of the potentialŴ ,ϕϕ . Along ϕ 0 the unit vector tangential to the inflationary trajectoryσ I points in χ direction. Consequently, in view of (39), along ϕ 0 the perturbation δϕ is directly related to the isocurvature perturbation Q s andŴ ,ϕϕ to the effective isocurvature mass m 2 s defined in (43). While the inflationary trajectory alongφ = ϕ 0 = 0 is classically stable forχ χ c , in the vicinity of the critical pointχ c , the restoring classical force, which keeps the trajectory focused to the ϕ 0 attractor, is no longer sufficiently strong to counteract the diffusive force that originates from the unavoidable quantum zero-point fluctuations the trajectory experiences in the ϕ direction. Forχ <χ c , the isocurvature mass turns negative and the solution ϕ 0 becomes unstable. At the same time, the perturbation δϕ starts to grow and even dominates over the classical solution ϕ 0 .
In situations where the quantum diffusive force dominates the classical background drift, the formalism, in which the background dynamics of the scalar fields is considered to be independent of the time evolution of the quantum fluctuations, breaks down. Instead, the application of the stochastic formalism introduced in [50], which properly takes into account the back reaction of quantum fluctuations on the coarse grained classical background dynamics, is required. In the stochastic formalism, the dynamics of ϕ during the transition stage aroundχ c is determined by a probability density function (PDF) P (ϕ, N ) that specifies the probability of the field having the value ϕ at a time N .
The time evolution of P (ϕ, N ) is described by the Fokker-Planck equation, see e.g. [51], The right-hand-side of the Fokker-Planck equation (59) is characterized by two terms: a classical drift term with the coefficient D(ϕ, N ) and a quantum diffusion, or "fluctuation", term with the coefficient F(ϕ, N ).
In the context of the inflationary dynamics in the vicinity ofχ c , the drift coefficient D in (59) corresponds to the rate of change of the classical (averaged) field, while the fluctuation coefficient F corresponds to the rate of change of the variance. For the decomposition ϕ(N, x) =φ(N ) + δϕ(N, x) into a homogeneous backgroundφ(N ) and a Gaussian random fluctuation δϕ(N, x) with δϕ = 0, the coefficients are obtained as In the following we omit the bar over a background quantity. Assuming slow-roll in theχ and ϕ directions 5 , the equation of motion (52) reduces to the single-field dynamics of ϕ depending only parametrically onχ(N ), Since we are interested in the dynamics around ϕ = 0, Taylor expansion ofŴ , ϕ yields the linearized equation which determines the drift coefficient with the effective N -dependent mass The variance δϕ 2 after coarse graining over the horizon scale is given in terms of the power spectrum P ϕ (k) ≡ k 3 |δϕ k | 2 /(2π 2 ) of the scalar perturbation δϕ, The modes k crossing the Hubble horizon in the time interval −∆N contribute to the increment of the integral (64) by increasing its upper bound and satisfy k = aH = a * e (N −N * ) H at their respective moments of horizon crossing. Here a * is the value of the scale factor at some earlier time, N > N * . Treating H as a constant, 5 In fact, slow-roll is automatically realized by the phase of Starobinsky inflation in Stage 1 which sets the initial conditions for the subsequent evolution in Stage 2. Within the slow-roll approximation ϕ ≈ 0, ε H 3 and the term proportional to F ,χχ ϕ /F can be neglected in (52) we obtain d ln k ≈ −dN , which results in the diffusive random "quantum kicks" the coarse-grained field experiences from the Fourier modes δϕ k continuously crossing the Hubble horizon, Inserting (62) and (65) into (59), P (ϕ, N ) satisfies the Fokker-Planck equation The Fokker Planck equation (66) is solved by a Gaussian ansatz with a time dependent variance S(N ) = ϕ 2 (N ), Inserting (67) into (66) yields an equation for S, Hence, we may interpret the equation (68) for the variance S as determining the time evolution of effective amplitude of the scalar field by identifying [52], We apply the stochastic formalism for a period in which the quantum diffusive term H 2 /4π 2 dominates the classical term 2m 2 ϕ S/(3H 2 ). The stochastic stage starts at some pointχ >χ c at which m 2 ϕ falls from large positive values towards zero, and ends when it takes large negative values for someχ <χ c . The onset of the stochastic stage can be determined by the condition 2m 2 ϕ /3H 2 < 1. This ensures that for any inevitable initial increment of S(N ) ∼ H 2 /4π 2 the diffusive term dominates the classical one. Similarly, the end of stochastic phase can be roughly estimated by the condition 2m 2 ϕ /3H 2 ≈ −1. 6 If the isocurvature mass becomes more negative, the classical dynamics will dominate again and we resort back to the original equation (52). Thus, the duration ∆N of Stage 2, can be estimated by calculating the time it takes 2m 2 ϕ /3H 2 to change from 1 to −1. For the potential (21), the ratio 2m 2 ϕ /3H 2 is given by that, in turn, provides the estimate 6 A more precise discussion is presented in Sect. VII Along ϕ 0 , the difference ∆F can be estimated from Starobinsky inflation ∆F ≈ ∆N , leading to C. Stage 3: Peak formation and slow-roll After the isocurvature mass first crosses zero and becomes tachyonic atχ c , it grows again (but negatively) while the trajectory still follows the ϕ 0 solution (which is now a hill). Eventually its magnitude becomes again comparable with the magnitude of the quantum diffusive kicks (65) and the classical dynamics takes over again. At this point Stage 3 starts, the trajectory turns/falls into one of the ϕ ± v valleys and the growing isocurvature modes source the adiabatic modes, ultimately leading to a peak in the adiabatic power spectrum.
The wavenumber k p , at which the power spectrum features a peak, crosses the horizon shortly after the trajectory passesχ c . Therefore, a peak at a fixed k p requires thatχ c = 0. In addition, a narrow peak in the power spectrum centered at k p requires that stage two lasts for less than an efold ∆N ≤ 1, which according to (72) Hence, according to (55), a non-zerô χ c implies ξ 1. As we will see in Sect. VII, the numerical results confirm that the parameter combinations, which lead to a suitable peak in the power spectrum, respect these expectations.
Although, Stage 2 lasts for less than one efold, the stochastic treatment is essential as it provides the initial conditions for the dynamics in Stage 3 during which the power spectrum is amplified. The remaining inflationary dynamics until the end of inflation takes place inside the ϕ ± v valley and is again that of an effective single-field model with a slow-roll potential.

VI. PEAK FORMATION MECHANISM
The peak formation mechanism during Stage 3 is based on the "isocurvature pumping" effect described in [53] and explained in detail in the context of scalaron-Higgs inflation in [16]. Qualitatively, this mechanism can be easily understood from the dynamical equations for the scalar perturbations (40) and (41). Although the exact dynamics of the isocuravture modes (41) is more complicated in general, let's assume for purely illustrative purposes thaẗ Q s ≈ 0. In this case the solution to (41) would read For a positive isocuravture mass m 2 s > 0, the amplitude of the isocurvature modes Q s is exponentially damped. Conversely, for a tachyonic isocurvature mass m 2 s < 0, the amplitude of the isocurvature modes Q s is exponentially amplified. The total amplification depends on the amplification factor m 2 s /(3H 2 ) and the duration ∆N = N 2 − N 1 .
Next, let us look at the equation for the adiabatic modes (40). They are sourced by the product of the turn rate ω and the isocurvature modes Q s . Hence, only if ω = 0 and the amplitude of Q s is sufficiently large, the source term will have a sizable impact on the adiabatic modes. 7 Thus, in order for this multi-field amplification mechanism to be realized in a concrete model, several factors have to work together in a synchronized way. First, a sufficiently long inflationary phase with a tachyonic effective isocurvature mass is required for the isocurvature modes to grow. Next, this amplification must be transmitted to the adiabatic modes, which requires that the inflationary trajectory follows a curved path in the field space to yield a non-zero turn rate. To avoid an overamplification, the sourcing must terminate at some point, either via the vanishing turn rate or via the exponential suppression of the isocurvature modes which requires the effective isocurvature mass to turn positive again.
The two-field potential (21) satisfies all these requirements. During the transition from Stage 1 to Stage 2, the effective isocurvature mass turns tachyonic and the isocurvature modes start growing. At the same time, the stable ϕ 0 valley turns into an unstable hill, such that the trajectory eventually turns/falls into one of the ϕ ± v valleys, and thereby permits a sourcing of adiabatic modes and ultimately leads to the formation of a peak in the adiabatic power spectrum. During the fall the effective isocurvature mass turns positive again and leads to an exponential suppression of isocurvature modes during the subsequent slow-roll phase of inflation in ϕ ± v . This growth and subsequent damping of the isocuravture modes for Scenario I is shown in Fig. 3. It is obtained for the same parameter values as for those in Fig. 4 and Fig. 5 obtained in Sect. VII, so that its relation to the temporary tachyonic nature of the effective isocurvature mass and the peak in the power spectrum can be directly compared. 7 The case in which the effective isocurvature mass in (43) is dominated by the (positive) contribution of the turn rate (assuming a Gaussian profile function for η ⊥ (N ) = ω(N )/H(N )) arising from the curvature of the scalar field space manifold was studied in [22]. Even if under such assumption the analysis can be carried out in a model-independent way, it does not seem to be applicable to most realistic models in which the model-dependent potential dominates the effective isocurvature mass and is responsible for its tachyonic instability.  Fig. 4 shows that the isocuravture perturbation Qs starts growing when the isocurvature mass turns negative.
The efficiency and magnitude of the amplification is highly sensitive to the time the trajectory spends on the hill during its unstable dynamics in the vicinity of the critical pointχ c . Aroundχ c , the inflationary trajectory is no longer protected by the strong restoring classical force in ϕ direction and, therefore, is exposed to the diffusive quantum kicks driving it away from ϕ 0 = 0. Although this phase of unstable motion along the ϕ 0 direction is very short, the magnitude of the quantum kicks it gives to the background trajectory is crucial for determining the time for which the iscocurvature pumping lasts. Hence, the stochastic treatment is crucial for the quantitative treatment of the amplification mechanism, and is described in detail in Sect. V B. Note that this peak formation mechanism is very different from that in single-field models of inflation, see e.g. [25,27,54] for models in which the amplification is based on an intermediate ultra slow-roll phase resulting from a potential featuring an inflection point.

VII. NUMERICAL TREATMENT
As described in Sect. V and Sect. VI, the individual stages of the background dynamics can be easily understood qualitatively. However, a precise calculation including the dynamics of the perturbations has to be carried out numerically. The entire background dynamics is obtained by patching the numerical solutions of the equations in the individual stages in such a way that the preceding stages provide the initial conditions for the subsequent stages.
During Stage 1 and Stage 3, we numerically solve the exact equations of motion (52) and (51) for both scalar fields. In Stage 2, in which the stochastic formalism is used to describe the ϕ dynamics, we numerically solve the equations of motion Since the actual fall into the valley occurs during Stage 3, during Stage 2 both scalar fieldsχ and ϕ can be safely considered as slowly rolling with the inflaton unit vectorσ still pointing in theχ direction. Consequently, we neglect the ϕ 2 term in (51) during Stage 2. In order to patch the numerical solutions obtained in the different stages, we have to find the transition moments between them. During the first phase along ϕ 0 , the steep positive curvature of the potential along the ϕ direction provides a strong restoring force which immediately erases the effect of the continuous quantum kicks trying to drive ϕ away from ϕ = 0. The moment N 1 of the transition between Stage 1 and Stage 2 can, therefore, be inferred from the moment at which for the first time the effect of a quantum kick H/(2π) on ϕ = √ S will not be erased, i.e. when the drift term in (74) becomes comparable to the diffusive term for S(N 1 ) = H 2 (N 1 )/(4π 2 ). The resulting condition is solved numerically for N 1 as Since S = ϕ 2 is effectively zero before N 1 , the complete set of initial conditions which result from patching Stage 1 and Stage 2 read Stage 2 lasts until the curvature of the potential becomes dominant again (but this time with a negative sign). The time N 2 , at which Stage 2 ends, is determined numerically from the condition The initial conditions for Stage 3 are The numerical analysis confirms the analytic estimate (72) that the second stage typically lasts for less than one efold and that the values acquired by ϕ at the beginning of Stage 3 are very small. This a posteriori justifies the assumptions of slow-roll along ϕ in (61), the Taylor expansion of the potential in (62), and the Gaussian solution (67) to the Fokker-Planck equation (66). The numerical solutions for ϕ(N ),χ(N ) and a(N ) are then used in the equations for the perturbations (35), which are solved numerically with Bunch-Davis initial conditions imposed in the deep subhorizon regime. Finally, the power spectra (47)- (49) are computed numerically for the pivot scale k * = 0.002 Mpc −1 crossing the horizon at N = 60. We discuss the numerical results separately for Scenario I and Scenario II.

A. Scenario I
For the parameter values (56), for which Scenario I is realized, the effective single-field stage with the Starobinsky potential (58) along ϕ 0 will be inevitably attained, independently of the initial value ϕ i .
A successful phase of Starobinsky inflation during Stage 1 requires m 2 0 /M 2 P = 1.18 × 10 −5 to satisfy the COBE normalization (4). The parameters are further constrained by demanding that PBHs are produced within a given PBH mass window (constraining the peak location in the power spectrum) and that these PBHs lead to a significant amount of CDM as observed today (constraining the amplitude of the peak). For instance, as explained in the context of the peak formation mechanism in Sec. VI, demanding that the power spectrum peaks on specific scales constraints the ratio ξ/ζ that, in turn, determinesχ c and, hence, also the efolding number around which the sourcing takes place. Altogether, these constraints fix three out of the four parameters m 0 , ζ, ξ, and λ.
The numerical results in Fig. 4 and Fig. 5 are based on the parameter combination m 0 /M P = 1.18 × 10 −5 , ζ = 2.6 × 10 −10 , ξ = 50, and λ = 10 −5 . The inlay in the top left plot of Fig. 4 shows ϕ as a function of N for both trajectories and illustrates that the inflationary trajectory (red) does not immediately turn/fall into the ϕ + v valley (blue) which re-emerges atχ c . Due to its inertia, it stays on the unstable hill along ϕ 0 for a short period. The exact moment, at which the trajectory (red) turns/falls before it catches up with the valley (blue), depends on the diffusive quantum kicks it experiences and requires the stochastic formalism discussed in Sect. V B.
The top right plot in Fig. 4 shows ϕ as a function ofχ. It illustrates how the inflationary trajectory runs along ϕ 0 , performs a sharp turn/fall and, after several mild oscillations, tracks the valley solution ϕ + v . The bottom left plot in Fig. 4 shows the ratio of the effective isocurvature mass (43) and the squared Hubble parameter as a function of N . It becomes tachyonic when it first crosses zero atχ c , corresponding to the moment at which the curvature of the potential in ϕ direction m 2 s ∝Ŵ ,ϕϕ changes sign. It turns positive again as the trajectory is pushed away from ϕ 0 and approaches the ϕ + v valley where m 2 s turns positive again. The oscillations in m 2 s disappear after the trajectory settles down completely in the ϕ + v valley. The bottom right plot in Fig. 4 shows the slow-roll parameter ε H as a function of N . It remains well within the slow-roll regime ε H 1 during all three stages of the inflationary background evolution and at no point enters a regime of ultra-slow roll (in particular not during the turn/fall). This illustrates again that the multi-field amplification mechanism described in Sect. VI is essentially different from that based on an intermediate phase of ultra slow-roll in single-field models of inflation. After mild oscillations, ε H settles down to a higher value in ϕ + v compared to its value along ϕ 0 . This implies that the inflationary dynamics in Stage 3 is slightly faster than that in Stage 1 due to the steeper slope ofŴ (χ, ϕ + v (χ)) inχ direction compared toŴ (χ, ϕ 0 ).
The formation of a peak in P R and the absence of any feature in P h is shown in Fig. 5. The left plot of Fig. 5 shows the weak logarithmic k dependence of the power spectrum P R for large wavelengths (small k) during the first slow-roll phase along ϕ 0 in Stage 1 with the amplitude P R ≈ 10 −9 required for the consistency with CMB measurements. At smaller wavelengths (larger k), P R experiences a strong amplification leading to a peak centered around k p /k * ≈ 10 15 with amplitude P R ≈ 10 −2 . This peak corresponds to the modes which cross the horizon during the turn/fall of the inflationary trajectory. For modes that cross the horizon during the slow-roll phase in Stage 3, the amplitude of P R ≈ 10 −10 is slightly smaller than that for the modes that cross the horizon during Stage 1.
The amplification in P R is entirely due to the sourcing of scalar perturbations. This can be seen from the right plot of Fig. 5, which shows that the tensor power spectrum P h remains constant for lower k in Stage 1, and drops as H decreases when the inflationary trajec-tory settles in the ϕ + v valley. No amplification or any other feature is visible, which is again explained by the fact that the isocurvature amplification mechanism only affects the scalar power spectrum P R , cf. (40).

B. Scenario II
While the Starobinsky phase along ϕ 0 during Stage 1 in Scenario I is similar to that in Scenario II, the transition to Stage 3 is more violent in Scenario II, as the inflationary background trajectory falls much deeper from the central ϕ 0 valley into one of the ϕ ± v valleys. All plots in this Section are obtained for m 0 /M P = 1.18 × 10 −5 , ζ = 6 × 10 −10 , λ = 10 −5 and ξ = 200. The top left plot of Fig. 6 shows the two-field potential superimposed by the numerically calculated exact inflationary trajectory (red) and the analytic valley equation ϕ ± v (blue). The inlay plot shows that the trajectory sharply bends and falls into the ϕ + v valley. The top right plot of Fig. 6 shows the parametric dependence of ϕ onχ and the oscillations in ϕ direction before settling in the ϕ + v valley. The effective isocurvature mass m 2 s /H 2 shown in the bottom left plot of Fig. 6 remains positive during the slow-roll phase in Stage 1, crosses zero at the critical pointχ c (shown by the lower inlay plot), grows negative while the trajectory is still on the unstable hill ϕ 0 , crosses zero again, peaks during the steep fall into the ϕ + v valley, undergoes heavy subsequent oscillations until the trajectory settles in the ϕ + v valley, and acquires a large positive value in the second slow-roll phase along ϕ + v well inside Stage 3. The upper inlay plot shows the peak amplitude of the oscillations in the effective isocurvature mass.
The bottom right plot of Fig. 6 shows the slow-roll parameter ε H (N ), that remains almost constant during Stage 1, grows rapidly during the fall and performs subsequent oscillations until the trajectory has settled in the ϕ + v valley. During Stage 3, it remains constant but is larger than in Stage 1 due to the steeper slope ofŴ (χ, ϕ + v (χ)). While overall ε H (N ) is well within the slow-roll regime, in contrast to Scenario I, the slow-roll condition is temporarily violated during Stage 3 when the trajectory falls into the valley, as shown by the inlay plot. This again emphasizes the difference between the isocurvature pumping mechanism and the ultra slow-roll mechanism for peak formation. The log-log plot of P R in Fig. 7 shows the weak logarithmic k dependence for large wavelengths in Stage 1 with amplitude P R ≈ 10 −9 consistent with CMB measurements. The strong amplification leads to a peak centered around k p /k * ≈ 10 15 with amplitude P R ≈ 10 −2 . Well inside stage 3, the trajectory settles into the ϕ + v valley in which the value of the potential is considerably smaller than during Stage 1 and the overall power P R ≈ 10 −11 along ϕ + v is about two orders of magnitude lower compared to that in Stage 1.
The log-log plot of P h in Fig. 7 shows that P h remains constant for large wavelengths during Stage 1 and sharply drops when the inflationary trajectory falls into the ϕ + v valley. It shortly oscillates while the trajectory settles in the ϕ + v valley and continuously decreases as the trajectory approaches the end of inflation along ϕ + v .

VIII. PBH DARK MATTER FROM INFLATION
The formation of PBHs in the early Universe was proposed more than 50 years ago [55,56]. Since PBHs do not form by the gravitational collapse of a star or the merger of two neutron stars, PBH masses have a much wider mass spectrum than permitted by the Chandrasekar mass bound.
Therefore, lighter black holes decay earlier. Universe, they are referred to as PBHs. 8 Despite the wide range of theoretically permitted PBH masses, there are strong observational constraints on the PBH mass spectrum which is probed by many different sources on a huge range of scales. The absence of any strong gamma ray burst that would have been emitted by the black hole evaporation [58], provides a direct lower bound on the allowed PBH mass. Tighter constraints come from Hawking radiation of decaying PBHs that would have heated the Universe and thereby delayed the formation of chemical elements during the Big Bang Nucleosynthesis (BBN). In addition, decaying PBHs might leave relics of Planck mass. 9 Since the total energy density of such hypothetical relics cannot exceed the present CDM density, this also leads to constraints on the PBH mass spectrum. Further constraints unrelated to PBH evaporation come from gravitational lensing experiments, pulsar timing experiments, as well as CMB and LSS measurements. A recent overview of observational constraints can be found in [4].

A. PBH formation
There are various formation mechanisms for PBHs, such as a sudden pressure reduction during phase transitions [60], cosmic strings [61] and bubble collisions [62]. The most relevant formation mechanism, however, is provided by the collapse of large density perturbations. The probability of such perturbations with a large amplitude can be enhanced by peaks in the inflationary power spectrum of curvature perturbation [63][64][65][66]. Since CMB measurements only cover large wavelengths (3), PBHs offer a unique opportunity to probe inflationary models on smaller wavelengths [35]. A PBH forms whenever the overdensity exceeds a critical value δ c in a given (spherical) Hubble volume V H (t) := 4πr 3 H (t)/3. The Hubble radius is defined as r H (t) := 1/H(t) and the background density in a flat FLRW universe is given byρ(t) = 3H 2 (t)/(8πG N ). In the simplest picture, a PBH with r S = r H immediately forms once a perturbation δ > δ c with wavelength λ δ enters the Hubble horizon λ δ = r H . 10 In addition, the 8 Even though PBHs could have also formed during inflation, their number density would have been strongly diluted due to the rapid expansion of the Universe. 9 A consistent description of the final stage of black hole evaporation may require a more fundamental quantum theory of gravity, see e.g. [59] for a review about various covariant approaches to quantum gravity. 10 A more detailed study of the formation time and its impact on the PBH abundance can be found in [67][68][69][70][71].
PBH mass M PBH is assumed to be directly proportional to the horizon mass M H at the time of formation t f , The critical density δ c is determined by the Jeans-length criterion applied in an expanding FLRW universe. The critical density δ c , as well as the proportionality factor K, depend only on the FRLW background dynamics. For an equation of state p = ωρ, with constant ω determined by the dominant energy density fraction at the time of formation, δ c ≈ ω and K ≈ ω 3/2 [56]. In the radiation dominated phase ω = 1/3 this leads to δ c ≈ 0.33 and K ≈ 0.19.
In this simplified picture, the PBH mass only depends on the horizon mass and, hence, on the time of formation but not on the amplitude δ. However, numerical collapse simulations show that the PBH mass satisfies a more complicated critical scaling relation and depends both on the time of formation t f as well as on the amplitude of the overdensity δ [72,73], The parameters K, δ c and γ in (85) are determined numerically [74][75][76][77]. An analytic approach to determine the threshold δ c was proposed recently in [78].

B. PBH abundance
In order to calculate the PBH abundance, it is useful to define the fraction of the mass in the universe, which collapsed into PBHs at the time of formation, In the Press-Schechter formalism 11 β is calculated as [83] β Here P (δ, t f ) is the PDF of generating an overdensity with amplitude δ at the moment of formation t f . Assuming that the perturbations δ are independent random variables, they follow Gaussian statistics. 12 The lower integration bound in (87) is determined by the critical collapse density δ c . 13 The probability density of having an overdensity with amplitude δ is given by Hence those perturbations forming PBHs are very rare and lie in the tail of the Gaussian PDF (88). Calculating β from (87) requires calculating the variance σ 2 (t f ) = δ 2 (t f ) in (88). The Fourier transform of the density contrast δ(t, x) is given by In a homogeneous and isotropic FLRW universe the variance σ 2 (t) is completely determined by the power spectrum P δ (t, k) via the two-point correlation function Assuming that the overdensities δ k (t) arise from the comoving curvature perturbations R k (t) amplified during inflation, we need to relate P δ (t, k) with the inflationary power spectrum P R (t, k) that is defined by the two-point correlation function The linear relation between δ k in the radiation dominated era and R k is given by 14 The transfer function T (t, k) describes the sub-horizon dynamics k > aH of δ k (t) after horizon re-entry, while T (t, k) = 1 for superhorizon scales k < aH. Thus, the variance (90) is obtained as The integral (94) diverges at the upper integration bound for small wavelengths λ = 1/k. This is avoided by smoothing δ(t, x) with a unit normalized window function W (x − y, R) at a smoothing scale R, Physically, the coarse graining induced by the smoothing means that at every point x, the smoothed overdensity δ R (t, x) represents the average of δ(t, x) over a spherical region of radius R centered at x, i.e. the substructures in the overdensity δ(t, x) below the resolution scale R are smoothed out in δ R (t, x) by the averaging procedure. We choose a modified Gaussian window function W G in (95). Following the conventions in [71], the window function in Fourier space reads 15 The window function (96) strongly damps out contributions from modes much larger than the "smoothing mode" k R = 1/R. Since we assume that a PBH forms when the modes δ k (t) re-enter the horizon at t = t f , the smoothing mode should be identified with the comoving Hubble radius at formation With the choice of the smoothing scale (97), the subhorizon modes k a(t f )H(t f ) are strongly suppressed by the window function such that the transfer function in (94) effectively becomes T (t, k) = 1 for the superhorizon modes k a(t f )H(t f ). The variance (94) at t f , smoothed at the horizon scale, acquires the form In order to have a sizable mass fraction (87), the smoothed variance (98) must be sufficiently large. Since σ R (t f ) is strongly damped by the (k/k R ) 4 factor for modes k k R and by the window function W 2 (k/k R ) for modes k k R , a sufficiently large σ R (t f ) can only be realized for a power spectrum which features a strong amplification at k ≈ k R . 16 Thus, in case the inflationary power spectrum P R (t f,k ) features a sharp peak at k p , this peak scale should be arranged to be close to k R by tuning the parameters of the inflationary model. Moreover, the horizon mass M H (t f ) at the time of formation is related to the peak scale k p by [95], 15 Note the additional factor of 1/2 in the argument of the exponential in (96). For a comparison of the impact of different window function see e.g. [94]. 16 Note that smoothing has a negligible effect on sharp peaks in P R , see Appendix B and [88] for more details.
Here g is the effective number of relativistic degrees of freedom, and t eq and k eq are the time of matter-radiation equality and the mode which crosses the horizon at that time, respectively. Neglecting for a moment the improved critical scaling relation (85) and assuming that the horizon mass is directly proportional to the PBH mass, (99) gives a good first estimate for the approximate PHB mass as a function of the peak scale k p , In the transition from (99) to (100), we have used M PBH = 0.19M H , g(t f ) = 106.75, g(t eq ) = 3.36, h = 0.674 and Ω m = 0.315 [35]. Conversely, if we are interested in a particular PBH mass range, the relation (100) gives an estimate for the peak location in P R (t, k).
For example, if we demand that all CDM we observe today is made of PBHs, only very small windows for the PBH mass are observationally allowed, and determine the scale k p at which P R (t f , k) must be strongly amplified.

C. PBHs as CDM
In the case when sufficiently large number of PBHs formed in the radiation dominated era, they could make up a large fraction of the presently observed CDM content in the Universe [5,[96][97][98][99][100]. Particularly interesting is the PBH mass window accessible to gravitational wave experiments performed by LIGO [54,[101][102][103]. The PBH mass distribution f (M PBH ) derived in (A23) measures the fraction of the presently observed CDM contributed by PBHs with mass M PBH .
In case CDM is entirely made of PBHs, the total fraction F PBH defined in (A22) equals one. In case only a fraction of CDM is made of PBHs F PBH < 1. Thus, the trivial constraint on the PBH abundance is given by F PBH ≤ 1 because the total PBH density cannot exceed the present CDM density. Since the mass spectrum of PBHs is already considerably constrained on a broad range of scales, there are only a few PBH mass windows in which F PBH = 1 can be realized [4,5]: In addition to M I PBH and M II PBH , the detection of binary black hole merges at LIGO/Virgo has renewed interest in the possibility of a primordial origin of CDM for PBHs in the mass range Even if the possibility of explaining all the observed CDM by PBHs in the mass window M III PBH seems to be ruled out observationally, as F PBH 10 −2 ÷ 10 −3 [4,5,103], a peak leading to the production of PBHs in the mass range M III PBH would provide an inflationary explanation for the observed merger events. While there are interesting scenarios connected to a mixed contribution of PBHs and additional CDM particles such as WIMPs (see [5] for an overview), we primarily focus on the possibility to explain all of the observed CDM content by PBHs with F PBH = 1 in the mass windows M I PBH and M II PBH . Nevertheless, in Sect. IX, we demonstrate that there are the parameter combinations that permit a generation of the observationally allowed distribution f (M PBH ) in all mass windows (101)-(103), i.e. F PBH ≈ 10 −1 ÷ 1 in the mass windows M I PBH and M II PBH , and F PBH ≈ 10 −3 ÷ 10 −2 in the mass window M III PBH . The mass intervals (101)-(103) directly translate into k intervals in which the peak featured in P R , centered at k p , must lie: k III p ≈ 10 6 Mpc −1 .
Conversely, once a particular PBH mass M PBH has been chosen, the constraint F PBH = 1 directly translates into a constraint on the amplitude A p of the peak in P R at k p . For a given PBH mass M PBH , an analytical estimate of the peak amplitude A p leading to F PBH = 1 is provided in Appendix B.

IX. NUMERICAL RESULTS
We present our numerical results for the PBH mass distribution (A23) in the three mass windows (101)- (103). We show that for all three mass windows, there are parameter combinations which lead to F PBH = 1. Moreover, in all cases it is also possible to find parameters such that F PBH < 1, which is particularly relevant for the observationally most interesting LIGO mass window (103). Therefore, we first present our results for the LIGO mass range, show that the log-normal function (B2) provides an excellent fit to the peak of the numerically generated power spectrum, and finally compare the peak amplitude A p to the results obtained in [94]. For the mass windows (101)-(103), we illustrate that the numerical value for F PBH is highly sensitive to the parameter ξ.
In both plots f (MPBH) is evaluated for g = 10.75. Fig. 8 shows that the chosen model parameters generate a distribution f (M PBH ) consistent with observational constraints. Hence, the two-field extended Starobinsky model (10) with the potential (21) can explain the origin of the black holes involved in the LIGO observed merger events as PBHs formed due to the enhancement of the inflationary power spectrum on wavelengths smaller than the ones probed by CMB. The fit-parameters, obtained by fitting the mass distribution f (M PBH ) in the left plot of Fig. 8 to a log-normal distribution, are close to the values reported in [113] for mass distributions that would give rise to the events in the recently released GWTC-2 event catalog of the LIGO-Virgo Collaboration [114] and can easily be made compatible by a slight modification of the model parameters ζ, ξ and λ.
The total fraction F PBH is highly sensitive to the model parameters, in particular to ξ, as can be seen by comparing the two plots in Fig 8 where F PBH = 4.1 × 10 −3 for the left plot and F PBH = 0.9 for the right plot. Therefore, by fine tuning the model parameters any numerical value F PBH ≤ 1 can be obtained.
As illustrated in Fig. 9, the log-normal function (B2) closely fits the peak of the exact numerically obtained power spectrum, implying that the peak is well characterized by three parameters: the peak scale k p , the peak amplitude A p and the peak width ∆ p . In addition, by using the log-normal fit, the subsequent numerical integrations can be performed much more efficiently. Finally, the log-normal fit to the peak allows a direct comparison with [94]. In Fig. 1 of [94], f (M PBH ) for the LIGO mass window was also computed using the Press-Schechter formalism with a Gaussian window function and by assuming a log-normal shape for the peak in the power spectrum. In [94] the value of the total fraction was fixed to F PBH = 2 × 10 −3 , which is of the same order as F PBH = 4.1 × 10 −3 obtained in left plot in Fig. 8.
Moreover, a comparison with the results tabulated in Table I of [94] shows (first column) that A p must lie between 4.14 × 10 −3 ≤ A p ≤ 8.92 × 10 −3 in order to have F PBH = 2 × 10 −3 for 0.3 ≤ ∆ p ≤ 1. The values A p = 7.2 × 10 −3 and ∆ p = 0.709 obtained in Fig. 8 leading to F PBH = 4.1 × 10 −3 therefore provide an additional consistency check for the numerical evaluation of but only within the smaller mass windows M I PBH and M II PBH resulting from splitting (107), see [5]. We explicitly show that for appropriate parameter values, mass distributions f (M PBH ) with F PBH = 1 can be realized in both mass windows M I PBH and M II PBH . Note that recent data from the NANOGrav Collaboration [115] lends further support to the proposal that PBHs may constitute a large part (or the whole) of CDM with the dominant contribution to the mass function in the range 10 −15 ÷ 10 −11 M [116]. Thus, in light of the observational ambiguity regarding the strict upper bound on F PBH in the respective mass windows M I PBH and M II PBH , our main objective is not to derive stringent constraints on the model parameters for all the mass windows but to demonstrate that there are parameter values that lead to observationally viable mass distributions f (M PBH ).
In Fig. 10, we show an exemplary parameter combination for which P R peaks at k p ≈ 10 15 Mpc −1 and generates a significant amount of CDM in mass window M I PBH . The left plot in Fig. 10 shows a mass distribution f (M PBH ) which leads to F PBH = 0.69 for ξ = 38. The right plot in Fig. 10 illustrates the sensitivity of f (M PBH ) on ξ. For a ξ larger by only 0.5%, the amplification mechanism is already too strong and leads to the observationally unacceptable large value of F PBH = 1.4.
Similarly, Fig. 11 shows f (M PBH ) for two parameter combinations in the mass window M II PBH . The left plot in Fig. 11 shows an observationally viable scenario with F PBH = 0.5, while the right plot leads to a an unacceptable value F PBH = 2.1. Summarizing, the discussion  presented in this section illustrates that the model parameters can be adjusted such that a significant fraction of CDM (including all CDM) is made of PBHs in all the three mass windows (101)-(103).

X. CONCLUSIONS
The two-field dilaton extension of Starobinsky's inflationary model predicts a successful phase of inflation in perfect agreement with recent Planck measurements. At the same time, it is capable of predicting the presently observed dark matter content in our Universe in the form of PBHs. The generation of gravitational waves and CDM from PBHs is also possible in the Starobinsky supergravity theory with two-field double inflation [38,39].
It is interesting to compare our model to the model of scalaron-Higgs inflation [16]. First, we identify the non-minimally coupled scalar field extending Starobinsky's geometric model with a dilaton field, whereas in the scalaron-Higgs model this scalar field is identified with the Standard Model Higgs field. Second, we assume a dilaton-dependent scalaron mass M 2 (ϕ) in (11), which effectively introduces an additional parameter ζ.
For positive values of ζ the stable inflationary dynamics along ϕ 0 reduces to that of an effective single-field model with the same predictions as Starobinsky's model at the scales probed by the CMB. However, the inflation-ary trajectory remains stable only up to the critical point χ c , at which ϕ 0 turns into an unstable hill. The trajectory subsequently falls into one of the outer ϕ ± v valleys. As explained in Sect. VI, this feature of the multi-field dynamics leads to an amplification of the adiabatic power spectrum at small wavelengths. In contrast, the two-field potential in scalaron-Higgs inflation features an unstable hill-top along ϕ 0 for all values ofχ and, therefore, leads to an immediate fall of the inflationary background trajectory already at the very onset of inflation.
Even though the scalaron-Higgs model predicts the same values for the inflationary observables as Starobinsky's model at large wavelengths, it does not support the isocurvature sourcing mechanism for smaller wavelengths required for significant PBH production. The detailed conditions required for a successful realization of the multi-field amplification mechanism include the growth of isocuravture perturbations resulting from an intermediate phase in which the effective isocurvature mass becomes tachyonic and the sourcing of adiabatic modes resulting from a curved trajectory in the scalar field space geometry. This "isocurvature pumping" mechanism, already discussed in [16], is a genuine multi-field effect and is essentially different from other amplification mechanisms in the single-field models of inflation.
We emphasized the necessity of the stochastic treatment within the short (less that one efold) transition phase, where quantum diffusive effects dominate the inflationary dynamics, and set the initial conditions for the subsequent classical dynamics in which the peak in the adiabatic power spectrum is generated. Thus, the stochastic treatment is crucial for the production of PBHs and for a precise prediction of their contribution to the presently observed CDM. The importance of the stochastic formalism in the context of PBHs produced in multi-field models of inflation is also discussed in [51,117,118].
All numerical results presented in Sect. IX are obtained for the parameter values satisfying the condition (56), and, therefore, are realized in Scenario I. In Scenario I, the valley along ϕ 0 is a global attractor, so that all its predictions are independent of the initial conditions for the inflationary background trajectory.
Even if the exact inflationary power spectrum is known numerically, the calculation of the PBH mass distribution depends on the details of the formalism used. For generic power spectra with multiple peaks of different amplitudes and shapes, a choice of the window functions and the formalism (Press-Schechter [83] vs. peaks theory [79,80,82]) do, in general, also affect the result for the PBH mass distribution [71]. The PBH collapse process does not only depend on the time at which the PBHs form, but also on the amplitude of the density perturbations. This is taken into account by the critical scaling relation for the PBH mass [74]. In the present work, we assumed that PBHs form immediately once an overcritical density fluctuation enters the horizon, took into account the critical scaling relation (85) for the PBH mass, worked with the Press-Schechter formalism (87), and used a Gaussian window function (96). However, since the dilaton-extended Starobinsky model predicts an almost scale invariant power spectrum with a single sharp peak, we neither expect our results to strongly depend on the specifics of the Press-Schechter formalism nor on the choice of the Gaussian window function. This expectation is further supported by the comparison with the results of [94].
We investigated three different PBH mass windows (101)- (103) where the observational constraints permit a sizable contribution of PBHs to the presently observed CDM content of our universe. This also includes the possibility to explain all CDM by PBHs in the mass windows (101) and (102), which is supported by the recent data from the NANOGrav Collaboration [115].
Of particular interest is the mass window (103), as the sources of the binary merger events observed by the LIGO Collaboration [101] may be identified with PBHs and thereby provide an additional observational window to the inflationary dynamics. While the maximum value of the fraction F PBH in the LIGO mass window (103) is controversially discussed [105][106][107][108][109][110][111][112], the extended dilaton two-field Starobinsky model can account for any observationally viable F PBH by a suitable combination of parameters.
Our model allows for several future applications. It would be interesting to study the effect of non-Gaussianities, both primordial ones and those arising from the non-linear relation between the comoving curvature perturbation and the density perturbation. Moreover, the formation of primordial black holes inevitably leads to the production of gravitational waves [95], which may be detected by the space-based gravitational interferometer LISA [23,119,120] and could provide additional constraints on the model. Finally, it may be possible to realize the same mechanism of PBH production with the Standard Model Higgs field instead of the dilaton in a similar extension of the scalaron-Higgs model [121].
The total fraction of CDM made from PBHs is obtained by integrating (A16) over the logarithmic mass intervals (or, equivalently, by integrating (A15) over all epochs when PBHs could have formed) as In order to obtain the distribution (A16) as a function of the PBH mass, we proceed as follows. The mass fraction β(M H ) as a function of M H is given by [72], with the Gaussian PDF (88) and variance σ 2 R (M H ), . (A19) The explicit dependence of M PBH on the horizon mass M H and the amplitude of the energy density contrast δ is given by the critical scaling (85), For fixed M H , we solve (A20) for δ as a function of M PBH , Hence, for fixed M H , the dδ integral in (A18) can be traded for the d ln M PBH integral with dδ = (µ 1/γ /γ)d ln M PBH . After doing that and using (A19) in (A18) and (A18) in (A16), (A17) takes the form Under several simplifying assumptions, a quick analytic estimate of β(M H ) can be obtained as follows. First, we calculate the variance σ 2 R by assuming that the inflationary power spectrum can be written as a sum of the constant CMB part A CMB and the part P PBH (k) responsible for the PBH production 18 , P R (k) = A CMB + P PBH (k). (B1) Assuming further that the PBH part of the power spectrum has a single symmetric peak centered at the peak scale k p , we parametrize the shape of the peak by a lognormal distribution, i.e. by a Gaussian in ln(k/k p ) with the standard deviation ∆ p centered around k p , In the case where an exact numerical treatment features a single peak, the parameters A p , ∆ p , and k p are extracted by fitting the log-normal distribution (B2) to the numerically obtained power spectrum. Here we focus on a derivation of the simple analytic estimate of the peak value A p required for f (M PBH ) = 1. We assume the simple scaling M PBH = KM H (i.e. not the critical scaling (85)) and use (A15) to related β with f . We further 18 The weak logarithmic scale dependence of P R at large scales is neglected as it has a negligible effect on the peak analysis.
assume that the peak is sufficiently sharp and can be approximated such that we can consider the limit ∆ p → 0 where the peak is described by a Dirac delta function, 19 P PBH (k) = A p δ (ln k − ln k p ) .
Since the generation of a significant number of PBHs requires A p A CMB , we can safely neglect the constant amplitude A CMB ≈ 10 −9 for the derivation of the PBH abundance and obtain P R (k) ≈ A p δ (ln k − ln k p ) . (B4) Using (B4) in (98) for a peak scale k p = k R 20 , we obtain a simple relation for the variance .
(B6) 19 See also [69,94,123,124] for related discussions as regards the impact and comparison of different peak widths. 20 As explained in Sect. VIII B, a strong amplification with a sufficiently large σ R can only be realized if the peak scales kp is not significantly different from the smoothing scale k R . Since we also assume that the time of formation coincides with the time the density perturbation enters the horizon, k R is identified with the comoving Hubble radius at the time of formation as in (97). (B7) A quick estimate of A p , required to obtain f (M PBH ) = 1 for a given PBH mass M PBH , can now be obtained easily. For example, let us take the values for the proportionality factor and the critical density found in [56] for a radiation dominated universe as K = 0.19 and δ c = 0.33, and choose a PBH mass M PBH = 5 × 10 −12 M . Then, according to (A15), the mass fraction should be β = 1.69 × 10 −14 in order to have f (M PNH ) = 1. Using (B6) and (B5), the corresponding peak value for the power spectrum can be deduced to be A p ≈ 1.7 × 10 −2 . Compared to the value of A p ≈ 5 × 10 −3 required to obtain F PBH ≈ 1 in Fig. 11 (the peak values in f (M PBH ) almost coincide with the total integrated fraction F PBH ), the analytic estimate roughly coincides up to an order of magnitude, with the main difference resulting from the omission of the critical scaling relation (85) in the analytic estimate.