Scalaron–Higgs inflation reloaded: Higgs-dependent scalaron mass and primordial black hole dark matter

We propose an extension of the scalaron-Higgs model by a non-minimal coupling of the Standard Model Higgs boson to the quadratic Ricci scalar resulting in a Higgs-dependent scalaron mass. The model predicts a successful stage of effective single-field Starobinsky inflation. It features a multi-field amplification mechanism leading to a peak in the inflationary power spectrum at small wavelengths which enhances the production of primordial black holes. The extended scalaron-Higgs model unifies inflationary cosmology with elementary particle physics and explains the origin of cold dark matter in terms of primordial black holes without assuming any new particles.


Introduction
In view of the countless number of inflationary models, predictability as well as theoretical motivation become even more important. In this respect, Starobinsky's quadratic f (R) model and the model of Higgs inflation stand out. Starobinsky's model is the natural geometric extension of Einstein's theory in which the higher derivatives in the quadratic curvature scalar lead to the emergence of an additional massive scalar propagating degree of freedom, the scalaron, driving inflation [1]. In contrast, in Higgs inflation, the inflaton is identified with the Standard Model (SM) Higgs boson nonminimally coupled to the Ricci scalar [2]. Quantum corrections dominated by the heavy SM particles [3] establish the connection between particle physics at the electroweak scale and cosmology at the energy scale of inflation via the renormalization group (RG) running [4][5][6][7][8], see [9] for a recent review on the Higgs field in cosmology. Both models lead to a e-mail: anirudh.gundhi@phd.units.it b e-mail: christian.steinwachs@physik.uni-freiburg.de (corresponding author) the same predictions for the spectral inflationary observables [3,10,11]. This is a particular manifestation of a more general equivalence between f (R) theories and scalar-tensor theories for different field parametrizations at the classical and quantum level [12][13][14].
In contrast to single-field models of inflation [43][44][45][46][47][48], multi-field models offer additional mechanisms for the amplification of the adiabatic power spectrum due to the multidimensional potential landscape and the curved field space geometry. In particular, such an amplification can result in the formation of peaks in the adiabatic power spectrum at small wavelengths. For peaks with sufficiently large amplitudes, the production of Primordial Black Holes (PBHs) is greatly enhanced and might explain the observed Cold Dark Matter (CDM) content of the Universe, see [49][50][51] for a review. The production of PBHs in multi-field inflation has recently been realized in a variety of models, see e.g. [42,[52][53][54][55][56][57][58]. The formation of PBHs by an effective single-field ultraslow roll mechanism has also been investigated in a finetuned RG-analysis of the scalaron-Higgs model [31]. Aside from offering an explanation for the origin of the observed CDM, PBHs offer a unique opportunity to constrain the adiabatic power spectrum at small wavelengths not accessible to Cosmic Microwave Background (CMB) measurements.
In this article, we extend the scalaron-Higgs model by a Higgs-dependent scalaron mass crucial for the successful realization of the multifield "isocurvature pumping" amplification mechanism described in [18,42]. This results in the formation of a single peak in the adiabatic power spectrum triggering the formation of PBHs at small wavelengths. We find simple scaling relations among the model parameters, which allow to adjust the amplitude and the location of the peak in the power spectrum, such that a maximum total PBH-CDM mass fraction can be realized in all observationally viable mass windows.
The model proposed in this article provides a unified description of inflationary cosmology, CDM, and elementary particle physics without assuming any new physics except for a non-minimal coupling of the SM Higgs boson to the modified gravitational R + R 2 sector.

Extended scalaron-Higgs model
We assume the SM embedded into curved spacetime with a modified graviton-Higgs sector given by The curvature-Higgs-dependent function in (1) reads The structure of the action (1) with the function (2) has been investigated in [42] in the context of an abstract dilaton field ϕ. In contrast, in this article, we study the implications of identifying ϕ with the SM Higgs field. This model provides a natural extension of the scalaron-Higgs model considered in [15][16][17][18]. We formulate the model as a two-field scalartensor theory. 1 The auxiliary field χ emerges when formulating the f (R, ϕ) theory as a two-field scalar-tensor theory [18]. Transforming in addition to the the Einstein frame (EF) by performing the non-linear field redefinitions the action (1) is written as two-field action in the EF Here,χ is the scalaron, effectively emerging from the higher derivatives present in R 2 term in (2). A more detailed presentation of the transition from (1) to (4) can be found in [18]. The local scalar field coordinates I (x) and the metric G I J on the scalar field-space manifold are defined by The scalar two-field potentialŴ ( ) in the EF readŝ We have introduced the parametrization The extended scalaron-Higgs model is defined by Here, U (ϕ) corresponds to an effective Higgs-dependent gravitational constant with the reduced Planck mass M P = 1/ √ 8π G N ≈ 2.4 × 10 18 GeV (in natural units c =h = 1) and the non-minimal coupling ξ , while M(ϕ) leads to an effective Higgs-dependent scalaron mass m(ϕ) with the constant scalaron mass m 0 and non-minimal coupling ζ . The potential V (ϕ) is defined by the SM Higgs potential with the quartic self-coupling λ and the symmetry breaking scale ν ≈ 246 GeV. The explicit form of the functions in (9)-(11) may also be motivated by the their limits for small and large values of the Higgs field. In the limit ϕ/M P → 0, the original Starobinsky model [1] is recovered, while for ϕ/M P → ∞, the model features an asymptotic scale invariance. The main generalization compared to the scalaron-Higgs model investigated in [18] is the Higgs-dependent function (10). In view of ν/M P ≈ 10 −16 , we neglect the constant ν for the inflationary analysis. Under these assumptions, the inflationary EF two-field potential (6) reduces tô

Covariant multi-field formalism
Following the general treatment in [18], we formulate the inflationary dynamics of the background and the perturbations in terms of the covariant multi-field formalism. The line element of the perturbed flat Friedmann-Lemaître-Robertson-Walker (FLRW) universe reads Here, t is the cosmic time, a(t) is the scale factor, i, j, . . . = 1, 2, 3 are spatial indices, δ i j = diag(1, 1, 1) is the flat spatial metric and E i j := ψδ i j + E ,i j . The scalar metric perturbations A(t, x), B(t, x), ψ(t, x), and E(t, x) combine with the scalar perturbations δ I (t, x). The Friedmann equations and the Klein-Gordon equations for the homogeneous scalar field multiplet I (t) read The dot is shorthanded for ∂ t . The Hubble parameter H (t) and the covariant time derivative D t are defined by The connection I J K is defined with respect to (5) and the unit vector along the background trajectory readŝ The unit vectorŝ I orthogonal toσ I satisfies The unit vectorŝ I is proportional to the acceleration vector ω I which defines the turn rate ω, Instead of the perturbations δ I (t, x), we work with the gauge-invariant Mukhanov-Sasaki variable [59][60][61], The equation for the Fourier modes of the perturbation δ I g (t, k) is found to be [61][62][63], Following the conventions introduced in [18], I J and the effective mass tensor M I J are defined by Here R I J K L is the Riemannian curvature tensor associated with the curved scalar field space manifold. Projecting (21) alongσ I andŝ I defines the adiabatic and isocurvature perturbations Inserting δ I g = Q σσ I + Q sŝ I into (22), the dynamical equations for the Fourier modes Q σ (t, k) and Q s (t, k) in the large wavelength limit k a H read The effective masses σ σ and m 2 s are defined by projecting (23) and (24) respectively and include additional contributions of the turn rate The operator f (d/dt) in (26) is defined by Only if the combination of ω and Q s is sufficiently large, Q σ is sourced by the "isocurvature pumping" mechanism discussed in [18]. This amplification may lead to a peak in the adiabatic power spectrum that is crucial for the formation of PBHs [42]. The power spectra of the scalar perturbations with the first two Hubble slow-roll parameters defined by Fig. 1 Bird's eye view of the EF two-field potential (12). The red lines sketch the inflationary trajectories running along the valleys in the direction of the arrows. Left: In Scenario I, ϕ 0 is a global attractor as the two valleys ϕ ± v merge with ϕ 0 . Right: In Scenario II, the valleys ϕ 0 and ϕ ± v never merge

Inflation and peak formation
The landscape of the two-field potential (12) is characterized by three valleys at ϕ 0 = 0 and ϕ ± v , which are solutions ϕ(χ) of the valley equation The model features two different scenarios which are characterized by the parameter combination The two scenarios shown in Fig. 1 are distinguished by the conditions [42], x ≥ 1 Scenario II.
We restrict our analysis to Scenario I. In contrast to Scenario II, the global attractor nature of the ϕ 0 = 0 solution in Scenario I ensures that the inflationary background trajectory is independent of the initial conditions. In Scenario I, the first stage of inflation proceeds along the ϕ 0 = 0 attractor. It starts atχ i and lasts until shortly before the critical valueχ c is reached, at which the local ϕ 0 minimum turns into an unstable maximum [42], During that stage the inflationary dynamics reduces to an effective single-field model with Starobinsky potential Consequently, ifχ c is sufficiently small, the predictions for modes probed by the CMB radiation are that of Starobinsky's model. For the CMB modes (38), the scalar and tensor power spectra only feature a weak logarithmic k dependence parametrized by the power-law ansatz Here, k * is a pivot scale which first crosses the Hubble horizon N * efolds before the end of inflation. The CMB predictions for A R , n R and the tensor-to-scalar ratio r = A h /A R in Starobinsky's model read [1], Planck data [64] constraints A * R , n * R at k * = 0.05 Mpc −1 , and the tensor-to-scalar ratio r * at k * = 0.002 Mpc −1 , For N * = 50 ÷ 60, the predictions (40) for n * R and r * are in perfect agreement with (42) and (43). For N * = 60, the normalization (41) fixes the scalaron mass to be According to (36), by tuning the ratio ξ/ζ for fixed m 2 0 , the value ofχ c can be made sufficiently small such that all CMB modes (38) cross the horizon before the inflationary trajectory passes the critical pointχ c . In this way consistency with the observational constraints (41)-(43) on the spectral CMB observables is ensured.
It is important to note that the relations (40) only approximately hold in our analysis. This is primarily because the slow-roll dynamics along ϕ ± v , after crossing the critical point χ c , is faster compared to that in Starobinsky inflation. This implies thatχ(N * ) >χ Star (N * ), whereχ Star (N * ) represents the numerical value ofχ in pure Starobinsky inflation at N = N * . Therefore, in our model, k * would 'feel' a flatter part of the Starobinsky potential at the time of horizon crossing. The spectral index would be closer to one (closer to scale invariance), and therefore slightly higher than that predicted in Starobinsky inflation. We found this deviation to be nonnegligible only for the largest values of λ ≈ 10 −3 − 10 −2 considered in our model. Nevertheless, the spectral index would still be compatible with CMB constrains by choosing N * to be closer to 50 efolds for these large values of λ. In addition, the scalaron mass must be slightly re-adjusted compared to the value (44) for different values of λ. However, as for the spectral index, these numerical changes are insignificant for all parameter values considered in this work. We therefore fix N * = 60 and m 0 = 1.18 × 10 −5 M P for all our numerical analyses.
The details of the dynamics in the vicinity of the critical pointχ c are crucial for the formation of a peak in the adiabatic power spectrum. Before reachingχ c , the unit vectorσ I points inχ-direction and δϕ is directly related to the isocurvature perturbation Q s which is suppressed due to a large and positive m 2 s , defined in (28). Atχ c , the effective isocurvature mass m 2 s becomes tachyonic, and, according to (27), leads to a strong growth of the isocurvature modes Q s . However, already shortly before the critical point is reached, the classical restoring force in ϕ-direction (proportional to m 2 s ) decreases and becomes comparable to the unavoidable zeropoint fluctuations δϕ driving the trajectory away from the ϕ 0 attractor. During this transition region aroundχ c the diffusive quantum effects become important (even dominant) and the inflationary dynamics must be described in terms of a probability density function (PDF) P(N , ϕ) within the stochastic formalism [65]. 2 The PDF gives the probability of the field having the value ϕ at time N and is determined by the Fokker-Planck equation [42], Here, dN := −H dt counts the number of efolds such that N = 0 at the end of inflation. A solution to the Fokker-Planck equation (45) is provided by a Gaussian ansatz with time dependent variance S(N ) = ϕ 2 (N ), The resulting first-order equation for the variance, effectively determines the time evolution of the background field ϕ(N ) via the identification [70], The stochastic phase during which the quantum diffusive term H 2 /4π 2 dominates the classical term 2m 2 ϕ S/(3H 2 ) in (47) lasts for a period N estimated by [42], As discussed in [42], for generating a sufficiently sharp peak in the adiabatic power spectrum, the duration N must be smaller than one efold N 1.
Once the inflationary trajectory has been driven away from the ϕ 0 solution, it turns and falls into one of the ϕ ± v valley. According to (26), the combination of the non-zero turn rate ω and the amplified isocurvature modes Q s leads to a sourcing of the adiabatic modes Q σ responsible for the formation of a peak in P R (k). The peak formation at the modes k p corresponding to PBH masses in the two different mass windows M I I PBH and M I I I PBH (defined in Sect. 5) is shown in Fig. 2 for λ = 10 −2 . After the fall, the background trajectory ultimately settles in one of the outer ϕ ± v valleys in which the inflationary dynamics again reduces to an effective single-field model with a second phase of slow-roll inflation. Inflation ends at χ f determined by the condition ε H (ϕ ± v ,χ f ) = 1, close to the global minimum atχ = 0. The exact inflationary background dynamics is obtained numerically by patching the following three stages as described in detail in [42]: (Stage 1) Effective single-field Starobinsky inflation along ϕ 0 forχ >χ c . (Stage 2) Stochastic phase in the vicinity ofχ c . (Stage 3) Fall into the ϕ ± v valley and subsequent slow-roll inflation along ϕ ± v untilχ f .

PBH production and CDM
The generation of a peak in P R (k) centered at k p , triggers the production of PBHs with a mass distribution f (M PBH ) centered at M PBH corresponding to k p .
We assume that a PBH directly forms once an overdensity greater than some critical value δ c enters the horizon The PBH mass is given by the critical scaling [71,72], Here, M H (t f ) is the horizon mass at the time of formation t f and the parameters K , δ c and γ in (51) are determined numerically [72][73][74][75]. As in [42], following the discussion below Eq. (2.16) in [76], we fix these parameters to be K = 10, δ c = 0.25 and γ = 0.36 consistent with the choice of the window function, as specified in the following discussion. In the Press-Schechter formalism, the PBH mass fraction β at t f is calculated by [77], The Gaussian probability of generating an overdensity with amplitude δ at t f is given by The variance smoothed over a scale The Gaussian window function in (54) reads [78], Trading We take the matter density parameter m = 0.315 and the CDM density parameter c = 0.264 [79]. The total integrated PBH-CDM mass fraction today is defined as Following the analysis in [42], the mass distribution as a function of the PBH mass M PBH is defined by F PBH := f (M PBH )d ln M PBH and obtained by using (51) A simple estimate for the approximate PBH mass as a function of the peak scale k p is obtained as [80], The PBH mass windows for which F PBH = 1 is compatible with observational constraints are [51,81], According to [51,[82][83][84][85], a maximal contribution F PBH 10 −2 ÷ 10 −3 in M I I I PBH seems to be favored. The relation (59) implies that the mass windows (60)-(62) are related to the following peak scales k p , We present the PBH mass distribution f (M PBH ) for the LIGO mass window M I I I PBH and the mass window M I I PBH in Fig. 3. Both mass distributions are obtained for the inflationary power spectra shown in Fig. 2 and are consistent with all observational constraints. The numerically generated f (M PBH ) (red dots) are fitted well by a log-normal distribution (blue line) defined by For completeness, we mention that there are also parameter combinations leading to a mass distribution f (M PBH ), which is compatible with observational constraints and leads to F PBH = 1 in the mass window M I PBH . Since recent data from the NANOGrav Collaboration [86] suggests that PBHs may constitute a large part (if not all) of CDM with a mass distribution centered in the mass window 10 −15 M ÷10 −11 M [87], the distinction between the mass windows M I PBH and M I I PBH might become obsolete.

Constraints on parameters
We have demonstrated that there are parameter combinations for which the extended scalaron-Higgs model describes a successful phase of inflation yielding predictions for the spectral observables (40) that are in perfect agreement with Planck data (41)- (43). At the same time, it explains the observed CDM in terms of PBHs. The compatibility with CMB measurements requires fixing m 0 according to (44). The quartic Higgs coupling λ is determined by SM physics. The remaining free parameters ζ and ξ are fixed by the properties of the peak in P R (k) which ultimately determines the PBH fraction of the observed CDM. In this section, we derive general scaling relations among the parameters λ, ζ , and ξ , for a given F PBH with underlying mass distribution f (M PBH ) centered at a given M PBH . First, we derive a direct relation between the model parameters ζ and ξ and the PBH mass around which f (M PBH ) is centered. Using the definition of the number of efolds N * − N = ln a/a * , for modes which cross the horizon at k = a H and k * = a * H with constant H ≈ H * respectively, we obtain N * − N = ln k/k * . During the phase of effective Starobinsky inflation, there is a simple relation between N and the field valueχ given by N (χ) ≈ F(χ). The peak in P R (k) is centered around the modes k ≈ k p ± k which cross the horizon in the vicinity ofχ c . Hence, we can express k p in terms ofχ c defined in (36) by the relation Since M PBH is related to the peak scale k p via (59) we finally obtain M PBH in terms of the model parameters, the pivot scale k * and the total number of efolds N * , For N * = 60, k * = 0.002 Mpc −1 , and m 0 as in (44), we obtain the linear scaling relation between ζ and ξ with a M PBH -dependent proportionality coefficient Next, we obtain a relation involving λ and ξ from the requirement that F PBH acquires a specific value for a mass distribution f (M PBH ) centered around a given M PBH . In view of the complex inflationary dynamics around peak formation, described in Sect. 4, and the details involved in calculating F PBH described in Sect. 5, going beyond an order of magnitude estimate based on various simplifying assumptions seems to be illusive.
We start by noting that P R (k) is related to the perturbation Q σ (N , x) in position space by According to (25), Q σ is related to δχ and δϕ via Q σ = G I Jσ I δ J , with I and G I J defined in (5). During most of the inflationary dynamics along ϕ 0 and the later part of the dynamics in ϕ ± v , the inflaton vectorσ I points in theχ direction and Q σ exclusively receives contribution from δχ . Only during the short peak formation stage in the vicinity of χ c , where the trajectory turns andσ I has a non-zero component in ϕ-direction, Q σ also receives contribution from δϕ. Here, we assume that during this periodσ I points in the ϕdirection such that (σ ϕ ) 2 = F. 3 For modes k ≈ k p ± k, which cross the horizon during this period, (70) reduces to For a simplified treatment we take the sharp peak limit P R (k) ≈ A p δ(ln k − ln k p ) such that (71) becomes Although the slow-roll dynamics along ϕ ± v slightly differs from that of the effective Starobinsky inflation along ϕ 0 , for an order of magnitude estimate we use the background relations of Starobinsky inflation F(χ) ≈ N and ε H (N ) ≈ 1/N 2 evaluated at N c := N (χ c ), so that (72) reduces to 4 As discussed in Sect. 4, close toχ c quantum diffusive effects dominate and a stochastic treatment is required during which ϕ(N ) is identified with δϕ 2 (N , x) 1/2 . But even after the stochastic phase, during the fall from ϕ 0 to ϕ ± v , both δϕ and ϕ continue to grow together -δϕ becauseŴ ,ϕϕ is still negative, and ϕ because it moves away from ϕ = 0 to larger field values until it reaches ϕ ± v . However, just before the background trajectory settles in the ϕ ± v valley,Ŵ ,ϕϕ turns positive at the inflection pointŴ ,ϕϕ = 0 during the fall. This leads to a sudden stop of the growth of δϕ, while ϕ still continues to grow until ϕ = ϕ ± v . Hence δϕ 2 is bounded from above by the maximum distance between the two valleys ϕ ± v (χ max ) − ϕ 0 = ϕ ± v (χ max ) attained atχ max <χ c , A strong amplification of P R (k) requires a large δϕ and hence a large |ϕ ± v (χ max )|. The inequality in (74) can be parametrized by δϕ 2 ≈ α 2 |ϕ ± v (χ max )| 2 with α ∈ [0.1, 1]. 5 The analytic expression for |ϕ ± v (χ max )| 2 is found from (32) and the second equation in (74) as 6 4 To produce PBHs in the mass windows M I I I PBH and M I I PBH , we find N c ≈ 40 and N c ≈ 25, respectively. 5 Geometrically, the inflection point which lies between ϕ 0 and ϕ ± v cannot be too close to ϕ 0 = 0. In addition, the inertia of the background dynamics carries the trajectory along ϕ 0 even after reaching the bifurcation point shown in the left plot of Fig. 1, such that the fall into ϕ ± v happens only after the valleys reach a sufficient separation, justifying the lower bound on α. 6 The criterion to determineχ max in (74) only applies to Scenario I. Only in this scenario, the valleys re-emerge at the bifurcation pointχ c turn and again move towards ϕ = 0. with x defined in (33) and the function L(x) defined by Since (34) implies x < 1, the function (76) takes arguments from x ∈ [0, 1) which means that (76) takes values in the interval L(x) ∈ [0, 1). For an order of magnitude estimate we approximate L(x) = x/4 + O(x 2 ) and obtain Combining (73) with (77), we obtain the analytic estimate for the peak amplitude We argued in Sect. 4 that N ≈ m 2 0 /(M 2 P ζ ) corresponds to the duration of the stochastic phase and that this phase must be sufficiently short N 1 in order to produce a narrow peak in P R (k). Since N c = O(10), the magnitude of the total prefactor in (78) is estimated to be of order N c N α 2 /8 ≈ 10 −2 , leading to the condition Since a significant F PBH ≈ 1 requires a peak amplitude A p ≈ 10 −2 ÷10 −3 [42,50], it is clear that x cannot be much smaller than one and we finally obtain the estimate Using (33), this yields the approximate scaling relation Inserting the scalaron mass m 0 from (44), we obtain The precise value of F PBH for a given power spectrum P R (k) ≈ A p δ(ln k − ln k p ) is exponentially sensitive to the peak amplitude A p , as can be seen from the relations (52)- (57). This is the main reason why any attempt to obtain a precise analytical relation for F PBH in terms of the model parameters is hard to realize. Nevertheless, in view of (79), the amplification only depends on x such that the same amplification is achieved for different values of ξ and λ as long as they are related by the scaling relation (82). In general the exact numerical factor in the quadratic scaling law (82), depends on the values of F PBH and the PBH mass M PBH at which the mass distribution f (M PBH ) peaks, but the scaling law λ ∝ ξ 2 will be the same for all mass windows and total mass fractions. As a side note, the relation (82) coincides with the CMB normalization condition found in pure Higgs inflation [2][3][4][5][6]. This coincidence is surprising, as in our model the CMB normalization condition (41) is satisfied by m 0 alone and the parameters λ and ξ are not directly related to CMB physics at large wavelengths but rather determine the PBH formation resulting from a peak in P R at small wavelengths.
Finally, we check the analytical estimates (69) and (81) by an exact numerical analysis. We systematically perform a parameter scan for different values of λ, ξ and ζ such that a mass distribution f (M PBH ) in the window M I I PBH centered around M PBH = 10 −11 M with F PBH ≈ 1, as shown in the right plot of Fig. 3, is realized. The parameters λ, ξ and ζ that permit such realizations, are related to each other by the scaling relations shown in Fig. 4, which are remarkably close to the analytical estimates (69) and (81). The linear and quadratic fits to the numerically found scaling relations in Fig. 4 are given by ζ = 4.23 × 10 −12 ξ, λ = 2.47 × 10 −9 ξ 2 .
In addition to the correct functional form of the scaling relations, also the numerical coefficients in (83) agree well with those predicted by the analytic estimates (69) and (82), thereby numerically confirming them. All parameters of the extended scalaron-Higgs model are fixed. The parameter m 0 is fixed by the CMB constraint (44) on the scalar inflationary power spectra at large wavelengths, independently of the value for the quartic Higgs coupling λ. In contrast, the non-minimal couplings ζ and ξ are ultimately determined in terms of λ by the scaling relations (69) and (81). The relations are determined by the requirement that the peak in P R (k) leads to a significant F PBH with a PBH mass distribution f (M PBH ) centered around M PBH .
In the SM, the tree-level value of the quartic Higgs coupling λ ≈ 10 −1 is determined by the symmetry breaking scale ν and the Higgs mass M h . In view of the huge energy gap separating the electroweak energy scale and the inflationary energy scale, the RG improvement becomes crucial to determine the value of the running Higgs coupling λ during inflation [4][5][6]. The contributions to the beta function of λ are dominated by quantum loops of the heavy SM particles. The system of the RG equations is highly sensitive to the precise conditions of the RG flow at the electroweak scale, in particular to the Higgs mass M h and the Yukawa top-quark mass, which are constraint by recent experimental bounds [88], M t = 172.9 ± 0.4 GeV.
While an analysis of the full RG system of the extended scalaron-Higgs model would be required for a precise determination of the running λ at the inflationary energy scale, already the pure SM running may be sufficient to derive a lower bound on λ. The RG running based on the SM beta functions drives λ to small values at high energies during inflation. Depending on (84) and (85), λ might even become negative and trigger an instability of the RG improved effective Higgs potential [89,90]. Assuming a stable λ > 0, by continuously varying the masses (84) and (85), the value of λ can in principle be made arbitrarily small at some energy scale. However, the smallest value of λ which can be sustained over the course of the inflationary dynamics is found to be λ ≈ 10 −6 [18,[91][92][93]. Hence, once the total mass fraction F PBH ≤ 1 and the PBH mass M PBH around which the mass distribution f (M PBH ) is centered are specified, the values of ζ and ξ are fixed in terms of λ via the scaling relations (69) and (82). The value of λ = 10 −2 ÷ 10 −6 during inflation, in turn, is fixed by SM physics at the electroweak scale. In this way, our unified model incorporates the physics of the SM at the electroweak scale, explains the presently observed CDM content of the Universe by PBHs and leads to inflationary predictions in agreement with measurements of the CMB radiation.

Conclusions
The extended scalaron-Higgs model proposed in this article is a viable model of inflation, which, at the same time, explains the origin of the presently observed CDM by PBHs.
One of the main features of the model is a Higgsdependent scalaron mass which arises from a non-minimal coupling of the SM Higgs field to the quadratic scalar curvature invariant. Compared to the scalaron-Higgs model [15][16][17][18], the additional non-minimal coupling ζ introduces one more parameter. With this additional parameter, the physics of the early Universe and the physics of the SM at the electroweak scale are described in one unified model which explains the observed CDM content without assuming any new particle, except for the scalaron which effectively emerges from the modified gravitational sector.
In addition, due to the global attractor nature of the ϕ 0 solution, the scenario considered in this article has the appealing feature that its predictions do not depend on the initial conditions of the inflationary background trajectory. A correct description of the background dynamics requires a stochastic treatment in the vicinity of the critical pointχ c . The inclusion of these diffusive quantum effects are crucial for an accurate quantitative treatment of the multi-field "isocurvature pumping" mechanism and leads to the formation of a peak in P R at small wavelengths responsible for a significant production of PBHs [42] .
We find that the extended scalaron-Higgs model can produce an observationally viable mass distribution f (M PBH ) with F PBH ≈ 10 −2 ÷ 10 −3 in the LIGO mass window (62) and F PBH ≈ 1 in the mass windows (60) and (61). We find simple scaling relations (69) and (81) between the nonminimal coupling parameters ζ , ξ , and λ, for CDM to be described by PBHs of a given mass. Together with the CMB normalization condition (41), which fixes the scalaron mass m 0 via (40), these scaling relations uniquely determine all parameters of the model in terms of λ. The quartic Higgs coupling λ, in turn, is determined by SM physics. The RG analysis of the SM suggests that a positive λ can take values λ ≈ 10 −2 ÷ 10 −6 at the energy scale of inflation. For these values of λ, the model permits a viable phase of inflation and an explanation of the observed CDM through PBH production in one single unified model without assuming any new physics.
The predictions of the model on wavelengths probed by the CMB are identical to that of Starobinsky's model. Thus, a measurement of the tensor-to-scalar ratio higher than that predicted by Starobinsky's model would rule out this model.
Another characteristic feature of the model is that it can only produce a single peak in the adiabatic power spectrum, such that F PBH can only receive contributions from PBHs within a narrow mass interval. Therefore, the model can also be tested against the possibility of F PBH collecting significant contributions from PBHs in different mass intervals.
Finally, since all parameters of the model are fixed by the SM (λ), the CMB (m 0 ), and a significant F PBH (ξ ) for f (M PBH ) centered at M PBH (ζ ), any additional prediction of the model, such as the production of gravitational waves accompanying the formation of PBHs [50,80,94,95], which may be detected by the space-based gravitational interferometer LISA [57,96,97], could potentially rule out this model.
Acknowledgements AG and CFS thank Sergey V. Ketov for discussions. AG thanks the University of Trieste and INFN for financial support and Angelo Bassi for supporting this collaboration.

Data Availability Statement
This manuscript has no associated data or the data will not be deposited. [Authors' comment: This is a theoretical work and no experimental data has been listed.] Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/. Funded by SCOAP 3 .