We study the formation of primordial black holes PBHs with ultraslowroll inflation when stochastic effects are important. We use the inlineformulammlmath displayinlinemmlmi mathvariantnormalmmlmimmlmiNmmlmimmlmathinlineformula formalism and simplify the stochastic equations with an analytical constantroll approximation. Considering a viable inflation model we find the spatial profile of the PBH compaction function numerically for each stochastic patch without assumptions about Gaussianity or the radial profile. The stochastic effects that lead to an exponential tail for the density distribution also make the compaction function very spiky unlike as assumed in the literature. Naively using collapse thresholds found for smooth profiles the PBH abundance is enhanced by up to a factor of inlineformulammlmath displayinlinemmlmsupmmlmn10mmlmnmmlmn9mmlmnmmlmsupmmlmathinlineformula and the PBH mass distribution is spread over 3 orders of magnitude in mass. The results point to a need to redo numerical simulations of PBH formation with spiky profiles.
['Cosmology', 'Astrophysics', 'and Gravitation', 'Astrophysics - Cosmology and Nongalactic Astrophysics', 'General Relativity and Quantum Cosmology']
Primordial Black Hole Compaction Function from Stochastic Fluctuations in UltraslowRoll Inflation
{'Primordial black hole compaction function from stochastic fluctuations in ultra-slow-roll inflation': "Sami Raatikainen, 1 Syksy Rasanen, 2 and Eemeli Tomberg 3, 4  \n1 University of Helsinki, Helsinki Institute of Physics, P.O. Box 64, FIN-00014 University of Helsinki, Finland 2 University of Helsinki, Helsinki Institute of Physics and Department  \nof Physics, P.O. Box 64, FIN-00014 University of Helsinki, Finland  \n3 Laboratory of High Energy and Computational Physics, National Institute of  \nChemical Physics and Biophysics, Ravala puiestee 10, 10143 Tallinn, Estonia  \n4 Consortium for Fundamental Physics, Physics Department, Lancaster University, Lancaster LA1 4YB, United Kingdom (Dated: September 19, 2024)  \nWe study the formation of primordial black holes (PBH) with ultra-slow-roll inflation when stochastic effects are important. We use the ∆ N formalism and simplify the stochastic equations with an analytical constant-roll approximation. Considering a viable inflation model, we find the spatial profile of the PBH compaction function numerically for each stochastic patch, without assumptions about Gaussianity or the radial profile. The stochastic effects that lead to an exponential tail for the density distribution also make the compaction function very spiky, unlike assumed in the literature. Naively using collapse thresholds found for smooth profiles, the PBH abundance is enhanced by up to a factor of 10 9 , and the PBH mass distribution is spread over three orders of magnitude in mass. The results point to a need to redo numerical simulations of PBH formation with spiky profiles.  \nIntroduction.Primordial black holes (PBHs) [1-9] can constitute all of the dark matter or seed supermassive black holes [10-12]. Generating PBHs requires large density fluctuations. They may originate in the same process as the small fluctuations seen in the cosmic microwave background (CMB). The most successful proposal for this is cosmic inflation, in the simplest case driven by one scalar field, the inflaton [13-27]. The curvature perturbation ζ is inversely proportional to the inflaton velocity, so a feature in the inflaton potential that slows down the field amplifies ζ . A much studied case is ultra-slow-roll (USR) inflation, where the classical force due to the potential is small or opposes the inflaton velocity, and the velocity decreases exponentially [4, 28-43].  \nA simple way to determine the PBH abundance is to solve for the inflaton background evolution, calculate the linear perturbations around it to obtain the curvature power spectrum P ζ , and, assuming Gaussian statistics, use P ζ to find the fraction of values of ζ that exceed the curvature threshold ζ th corresponding to PBH formation [31-35, 43]. This approach has two shortcomings.  \nFirst, stochastic effects due to the stretching of shortwavelength modes into long wavelengths couple the evolution of the background to the perturbations [44-48]. USR is particularly sensitive to stochastic noise, because the classical field velocity is small and perturbations are large [41, 42, 49-62]. Stochastic kicks lead to an exponential tail for the fluctuations, which can dominate the  \nlinear theory Gaussian tail [49, 53, 56, 58-61]. 1  \nIt is therefore not enough to consider the power spectrum. In [56] we presented the first numerical calculation of the non-Gaussian tail of the probability distribution of ζ due to stochastic effects in USR, using the formation threshold ζ th = 1. This brings us to the second problem: ζ is not a good measure of PBH formation. The curvature perturbation is sensitive to long-wavelength fluctuations, whereas gravitational collapse and formation of a trapped surface is determined by local conditions [72-74].  \nThe maximum value of the compaction function C ≡ 2 G N ∆ M/R or its average has been found to be a good indicator of PBH formation [75-92]. Here G N is Newton's constant, ∆ M ( t, r ) is mass excess over the background, and R ( t, r ) is the areal radius at time t and coordinate radius r . For a Gaussian field, large overdensities are close to spherically symmetric, and we assume this also in our non-Gaussian case [93]. On scales larger than the Hubble radius, C is independent of time.  \nChecking whether C ( r ) exceeds the collapse threshold requires knowing the spatial profile of the perturbation. Previous studies have used the mean profile and its variance calculated for Gaussian random fields [64, 80, 83, 86, 89, 91, 94-96], or ad hoc profiles [72, 75, 77-79, 81, 82, 84, 88, 92, 96-98]. The former approach is unlikely to be accurate when the statistics are highly non-Gaussian, and the ad hoc profiles have been  \nsmoothly varying and have had only one or two maxima. When stochastic effects are significant, C ( r ) varies rapidly and has many maxima.  \nIn this Letter we for the first time calculate the compaction function C ( r ) directly from stochastic inflaton fluctuations, showing how the stochastic effects that generate the exponential tail of the perturbations also make their spatial profile spiky. We use C ( r ) to determine PBH abundance as a function of the collapse threshold, providing a complete route from the inflaton potential to PBH abundance including stochastic effects, without assumptions about the radial profile or Gaussianity.  \nThe compaction function.As in [56, 59], we use the ∆ N formalism, where ζ is equal to the difference between the number of e-folds N of a local patch of the universe and the mean number of e-folds ¯ N , ζ = N -¯ N ≡ ∆ N [99102]. By simulating one patch of comoving scale k , we obtain the coarse-grained curvature perturbation ζ <k ≡ ∫ d 3 p (2 π ) 3 / 2 ζ p e i p · x θ ( k -p ), where k ≡ | k | , p ≡ | p | . Given spherical symmetry, the Fourier modes satisfy  \nζ k = ζ k = √ 2 π 2 k 3 d ζ <k dln k . (1)  \nWe end the simulations on super-Hubble scales, where ζ k is time-independent. We probe a large number of k values in the same patch to reconstruct the function ζ ( r ) with the inverse Fourier transform:  \nζ ( r ) = 1 (2 π ) 3 / 2 ∫ d 3 k ζ k e i k · x = √ 2 π ∫ ∞ 0 d k k 2 ζ k sin( kr ) kr . (2)  \nA window function is sometimes introduced in (2) to cut off scales and define the PBH mass [80, 96, 103]. However, its physical motivation is unclear. Using a Gaussian window function would not have a large effect on our results. Our PBH masses are determined by where the stochastic process produces maxima of C ( r ), not fixed a priori. The patches just have to be small enough to be approximated as spherically symmetric, and large enough to capture all maxima that can lead to PBH formation. We also do not include the transfer function, which has a small impact on sub-Hubble evolution [86, 96, 98, 104].  \nOn super-Hubble scales in the comoving gauge, the compaction function in the gradient approximation is [77]  \nC ( r ) = 2 3 [1 -(1 + rζ ' ) 2 ] , (3)  \nwhere ' ≡ ∂ ∂r . The prefactor is gauge-dependent [75, 77, 89, 90]. It has been argued that the areal volume average of C is a better indicator of collapse [82, 83, 85, 87, 88, 91],  \n¯ C ( r ) ≡ 3 R ( r ) 3 ∫ R ( r ) 0 d ˜ R ˜ R 2 C (4) = -2 r 3 e 3 ζ ( r ) ∫ r 0 d˜ r ˜ r 2 e 3 ζ [2˜ rζ ' +3(˜ rζ ' ) 2 +(˜ rζ ' ) 3 ] ,  \nwhere R ( t, r ) = a ( t ) re ζ ( t,r ) , and a ( t ) is the background scale factor. In the literature, the integral has been taken to the maximum of C ( r ); we keep the integration endpoint free and study the profile of ¯ C ( r ).  \nThe collapse threshold C th or ¯ C th depends on the radial profile, varying from 0 . 4 to 2 / 3. It was observed in [77, 78] that a smoother transition between the overdensity and the surrounding background corresponds to a lower threshold, while a sharper peak leads to a higher threshold, as pressure gradients resist collapse. In [83] it was argued that if C ( r ) has several peaks, the threshold should be compared to the value of C ( r ) at the first peak from the origin, as other peaks dissipate. In contrast, [105] found that high-frequency modulation of a peak facilitates collapse, and [92] found that with two peaks, collapse can occur even if both are below the threshold.  \nIn our case, gradients are large due to rapid stochastic fluctuations; a typical realisation of C ( r ) is shown in figure 1. Because the compaction function is conserved on super-Hubble scales, the spikes persist until the fluctuation re-enters the Hubble radius, when the collapse starts. Then pressure gradients smooth out variations [98]. Determining the damping and finding the collapse criterion for such sharp profiles requires a dedicated numerical analysis. We simply find the probability distribution for the maxima of C ( r ) and ¯ C ( r ), called C max and ¯ C max respectively, and discuss the PBH abundance for conventional values of the collapse threshold.  \nFor comparison, we also find the abundance using the mean radial profile calculated assuming Gaussianity. The mean profile of a spherical fluctuation with maximum value ζ 0 is ⟨ ζ ( r ) ⟩ = ζ 0 σ 2 ζ ∫ d k k P ζ ( k ) sin( kr ) kr , where σ 2 ζ ≡ ∫ d k k P ζ ( k ) [106]. From ⟨ ζ ( r ) ⟩ we calculate C ( r ) and ¯ C ( r ). The Gaussian distribution of ζ 0 with variance σ 2 ζ then gives the distributions of C max and ¯ C max .  \nFrom the inflaton potential to the compaction function.To compute ζ stochastically, we split the inflaton field into long and short wavelength parts, ϕ = ¯ ϕ + δϕ , separated by the coarse-graining scale k σ ≡ σaH 0 , where H 0 is the initial value of the Hubble parameter, and σ ≪ 1 is a constant [44-48]. The long-wavelength part ϕ is treated as a background, and evolves according to the following stochastic equations to leading order in the gradient expansion (see e.g. [52, 59, 60]):  \n˙ ¯ ϕ = ¯ π + ξ ϕ , ˙ ¯ π = -( 3 -1 2 ¯ π 2 ) ( ¯ π + V , ¯ ϕ ( ¯ ϕ ) V ( ¯ ϕ ) ) + ξ π , (5)  \nwhere ˙ ≡ d d N , and ¯ π is the mean field momentum, and we use units where the Planck mass is unity. The dynamics are determined by the potential V . We consider the Higgs-inspired inflation model of [56, 59] that produces asteroid-mass PBHs while also fitting CMB observations [107]. The perturbations are amplified when the inflation rolls over a local maximum in V . The effect  \nof perturbations is described by the noise terms ξ ϕ and ξ π , caused by the short-wavelength modes crossing the coarse-graining scale. They are treated as linear perturbations, so the noise is Gaussian, with ⟨ ξ ϕ ( N ) ξ ϕ ( ˜ N ) ⟩ = P δϕ ( N,k σ ) δ ( N -˜ N ), where P δϕ ( N,k ) ≡ k 3 2 π 2 | δϕ k ( N ) | 2 . Correlators that involve ξ π behave similarly.  \nIn [56, 59], the coupled stochastic evolution of ¯ ϕ and δϕ was solved at every timestep. This is computationally expensive: it took over 10 6 CPU hours to resolve the probability distribution up to ζ = 0 . 95. To be able to reconstruct the spatial profile of C ( r ), we simplify the calculation by using results from [61].  \nOur PBH model, as is typical, first has a USR phase, when the inflaton climbs towards a local potential maximum, and then a dual constant-roll (CR) phase where it rolls down the other side [108]. The perturbations that exit the Hubble radius during USR get the largest amplitude. They exit the coarse-graining scale later, during the dual CR phase. We denote the wavenumber of the last mode that exits during USR by k PBH = 2 . 0 × 10 13 Mpc -1 . During the dual CR phase the perturbations are adiabatic, meaning that the stochastic noise only moves the inflaton back and forth on the CR trajectory, where the perturbation evolution turns out to be independent of the noise [60]. As long as this holds, the equations (5) have the following analytical solution in terms of discrete k modes [61]:  \nϕ ( N ) = ϕ 0 e ϵ 2 2 N ( 1 -ϵ 2 2 X <k σ ) , (6)  \nX <k ≡ -k ∑ ˜ k = k ini √ P ζ ( ˜ k ) d ln k ˆ ξ ˜ k , ⟨ ˆ ξ k ˆ ξ ˜ k ⟩ = δ k ˜ k . (7)  \nHere ϕ 0 is an integration constant, and ϵ 2 ≡ d d N ln ( -˙ H/H ) is the second slow-roll parameter during the dual CR phase (when it is constant). The power spectrum P ζ ( k ) is pre-solved from the Mukhanov-Sasaki equation. In our model it peaks at k peak ≈ 0 . 2 k PBH , and ϵ 2 = 0 . 807 [59], see figure 3 in Appendix A.  \nThe sum in (7) is taken over ˜ k values with fixed logarithmic step d ln k , which approaches zero in the continuum limit, starting from the initial coarse-graining scale k ini = 0 . 05 Mpc -1 = 2 . 5 × 10 -15 k PBH . The stochastic noises ˆ ξ k are independent Gaussian random variables.  \nOn the CR trajectory, the coarse-grained variable ∆ N <k = ζ <k is given in terms of X <k as [61]  \nζ <k = -2 ϵ 2 ln ( 1 -ϵ 2 2 X <k ) . (8)  \nThe resulting probability distribution for ζ <k matches the numerical result of [56, 59] excellently, demonstrating the power of this technique.  \nThe results (6)-(8) characterise a single stochastic realisation completely in terms of the variables ˆ ξ k , a tremendous simplification compared to solving (5) by force. The  \nmodel is fully described by P ζ ( k ) and the parameter ϵ 2 that controls non-Gaussianity. Combining (6)-(8) with (1) and (2) and using Itˆo's lemma [61, 109, 110] for differentiation of stochastic variables gives  \nrζ ' ( r ) = k end ∑ ˜ k = k ini [ -ˆ ξ ˜ k 1 -ϵ 2 2 X < ˜ k √ P ζ ( ˜ k ) d ln k + ϵ 2 4 ( 1 -ϵ 2 2 X < ˜ k ) 2 P ζ ( ˜ k ) d ln k ] × [ cos ( ˜ kr ) -sin ( ˜ kr ) ˜ kr ] . (9)  \nThe final mode with k end = 2 . 6 × 10 16 Mpc -1 = 1330 k PBH exits the Hubble radius about ten e-folds before the end of inflation. Plugging (9) into (3) gives the profile C ( r ) in one patch. We collect statistics for the patches by generating many stochastic realisations ˆ ξ k .  \nThe quantity C was first calculated from stochastic inflation in [111]. The authors used the first passage time formalism to find the coarse-grained ζ <k , which they used to approximate ζ ( r ) coarse-grained over a region of size ∼ k -1 . They computed the distribution of C max analytically in an example model, with the noise given in the slow-roll approximation. We obtain the radial profile in each patch, and compute the noise beyond the slow-roll approximation, which can make a large difference [56, 59].  \nResults.We generated 10 9 stochastic realisations of ζ ( r ), each with 10450 k -points from k ini to k end , separated by step length d ln k = 1 / 256. We recorded the values of C max and ¯ C max , the values of r for which they were reached, the value of ∆ N <k PBH = ζ <k PBH (denoted simply by ∆ N below), and the number of peaks above C th > 0 . 4 or ¯ C th > 0 . 4. We discuss resolution issues and the numerical algorithm in Appendices B and C.  \nThe left panel in figure 1 shows a typical realisation of ζ ( r ), chosen so that C max = 0 . 4, and the corresponding C ( r ) and ¯ C ( r ). The curvature perturbation is constant for small r , drops rapidly to zero around r ∼ 10 k -1 PBH , and oscillates stochastically before settling down. The oscillations are amplified in C ( r ), and its peak locations correspond to the big drop of ζ ( r ). The behaviour of ¯ C ( r ) is similar, but smoother and more robust to small fluctuations. For example, it has only one high peak instead of two. The plots also show the mean profile ⟨ ζ ( r ) ⟩ in the Gaussian case with the same value of ζ (0), and the corresponding C ( r ) and ¯ C ( r ).  \nIn [92] a C ( r ) profile with two peaks was studied, with the caveat that such cases might be rare. We see that several peaks is the norm: 39% of the simulations even have more than one peak above C th = 0 . 4. The number p of peaks above C th = 0 . 4 follows an exponential distribution, n ( p ) ∝ e -αp , with α = 0 . 79. See Appendix B for dependence on resolution. Perturbations in all patches are of type I [76, 112], i.e. R ' > 0.  \nFIG. 1. Left: Solid lines show one realisation of the curvature perturbation ζ ( r ), the compaction function C ( r ) and its average ¯ C ( r ). The realisation has been selected so that C max = 0 . 4. Dashed lines show the same quantities calculated from the mean profile assuming Gaussianity and the same maximum value of ζ ( r ). The horizontal line shows the collapse threshold 0.4. The vertical line marks the power spectrum peak scale.  \n<!-- image -->  \nMiddle: Probability distribution of C max and ¯ C max . The solid lines show the results of our numerical computation, while the dashed lines correspond to the Gaussian mean case. The vertical line marks the collapse threshold 0.4. Right: The mass distribution with the collapse threshold C th = 0 . 4. The vertical line marks the power spectrum peak scale.  \nThe middle panel in figure 1 shows the probability distribution of C max and ¯ C max . Averaging C suppresses the maximum values significantly. Also shown are the distributions calculated from the mean Gaussian profile. They fall off steeply, and do not reproduce the numerical distributions anywhere. Spherical symmetry, and thus our distributions, is not likely to be valid for small fluctuations, but they are not important for PBH formation.  \nThe initial PBH abundance β is given by an integral over P ( C max ) or P ( ¯ C max ) starting from the collapse threshold. For C th = 0 . 4, we get β = 0 . 01, while ¯ C th = 0 . 4 leads to β = 5 × 10 -9 . (The abundance for ¯ C max has an error of some tens of % because of the small number of realisations above ¯ C th ; in contrast, C is well sampled.) Using the Gaussian mean profile for the same thresholds instead gives β = 4 × 10 -15 and β = 9 × 10 -18 , respectively. In [56, 59], we calculated β for the same model from the ∆ N distribution using the threshold ζ th = 1, giving β = 3 × 10 -11 . (In the Gaussian approximation C max = 0 . 4 corresponds to ζ 0 ≈ 1, and ¯ C max = 0 . 4 to ζ 0 ≈ 1 . 1.) The corresponding Gaussian linear perturbation theory result was β = 3 × 10 -16 , tuned so that PBHs contribute all dark matter today.  \nIn figure 2 we show the correlation of ∆ N with C max and ¯ C max . It demonstrates that ∆ N is a poor predictor of C max and vice versa; many patches with C max > 0 . 4 even have ∆ N < 0. The average ¯ C max is better correlated with ∆ N , but even there the scatter is large, and we are limited by lack of statistics at large values of ¯ C max .  \nThe PBH mass is roughly the mass inside the Hubble radius given by the peak scale r . There are correc-  \ntions from evolution during the collapse and accretion [72, 81, 92, 113, 114], including critical behaviour [115118]. Neglecting these and assuming standard thermal history, the mass is (see e.g. [108])  \nM = 5 . 6 × 10 15 ( r × 0 . 05 Mpc -1 ) 2 M ⊙ . (10)  \nThe right panel of figure 1 gives the mass distribution normalised to unity using C th = 0 . 4. It peaks around M = 1 . 5 × 10 -11 M ⊙ , compared to the monochromatic approximation in [56, 59], where the scale k PBH gave the mass M = 7 × 10 -15 M ⊙ . It has significant support over three orders of magnitude in mass, much broader than the distribution found in [111], stretching to higher values than that determined from the maximum of the power spectrum, whereas critical collapse extends the spectrum to smaller masses. For the threshold ¯ C th = 0 . 4, there are not enough realisations to obtain a reliable mass distribution, but the masses are of the same order.  \nConclusions.We have for the first time constructed stochastic realisations of the spatial profiles of the curvature perturbation ζ and the compaction function C , starting from the inflaton potential. The stochastic effects that have been shown to generate an exponential tail for the distribution of the curvature perturbation also lead to highly oscillatory profiles for the individual patches. Our calculation produces the individual profiles of each patch, allowing us to find the radial profiles of both C and the average compaction function ¯ C .  \nPBHabundance has sometimes been calculated assuming Gaussian statistics and using the maximum of ζ as  \nFIG. 2. Correlation of ζ = ∆ N with C max (left) and ¯ C max (right).  \n<!-- image -->  \nthe collapse criterion. Naively applying the PBH collapse threshold for C or ¯ C that has been calibrated on smooth profiles to our spiky profiles enhances the PBH abundance by a factor of 10 9 for C and a factor of 10 2 for ¯ C , compared to using ζ . Compared to calculating the compaction function from the mean profile, assuming the field to be Gaussian, as often done in the literature, the enhancement factor is 10 13 for C and 10 9 for ¯ C .  \nHowever, the collapse criteria from simulations of smooth profiles cannot be reliably applied to the rapidly varying profiles we find, as shown by the fact that we have large differences between using C or ¯ C . On the one hand, having more than one spike can lead to collapse even if none of them exceed the threshold [92]. On the other hand, pressure gradients will dampen the spikes when the perturbation enters the Hubble radius and starts to collapse, which is expected to wash away some of the enhancement. Checking the collapse criterion anew with numerical simulations using the spiky profiles should establish which quantity correlates well with collapse, and the degree of profile-dependence of the collapse threshold. The enhancement factor can also be different for inflaton potentials corresponding to different PBH masses, as is the case for the exponential tail [59].  \nThe tail and spike enhancements could be important for determining whether the non-detection of PBHs rules out explaining the origin of gravitational waves possibly detected by pulsar timing arrays with second-order production from scalar perturbations [96, 113]. With this mechanism, detection of the amplitude of gravitational waves would fix the amplitude of the scalar power spectrum. Because of the enhancements we have discussed, the power spectrum will lead to a larger PBH abundance, without changing the amplitude of gravitational waves. This can make the constraints tighter. It has also been argued that generation of PBHs in single-field inflation  \nis excluded because the required large small-scale power spectrum would distort the power spectrum on CMB scales [119-135]. The effects we discuss would alleviate such tension, but it remains to be quantified how much smaller the power spectrum can be. We will report details of this and other extensions of our work elsewhere.  \nWe thank Daniel Figueroa who participated in the beginning of the project. The authors wish to acknowledge CSC - IT Center for Science, Finland, for computational resources. This work was supported by the Estonian Research Council grant PRG1055 and by the EU through the European Regional Development Fund CoE program TK133 'The Dark Side of the Universe.' E.T. was supported by the Lancaster-Manchester-Sheffield Consortium for Fundamental Physics under STFC grant: ST/T001038/1.  \n- [1] S. Hawking, Mon. Not. Roy. Astron. Soc. 152 , 75 (1971).\n- [2] G. F. Chapline, Nature 253 , 251 (1975).\n- [3] A. Dolgov and J. Silk, Phys. Rev. D47 , 4244 (1993).\n- [4] P. Ivanov, P. Naselsky, and I. Novikov, Phys. Rev. D50 , 7173 (1994).\n- [5] J. Yokoyama, Astron. Astrophys. 318 , 673 (1997), arXiv:astro-ph/9509027 [astro-ph].\n- [6] J. Garcia-Bellido, A. D. Linde, and D. Wands, Phys. Rev. D54 , 6040 (1996), arXiv:astro-ph/9605094 [astroph].\n- [7] K. Jedamzik, Phys. Rev. D55 , R5871(R) (1997), arXiv:astro-ph/9605152 [astro-ph].\n- [8] P. Ivanov, Phys. Rev. D57 , 7145 (1998), arXiv:astroph/9708224 [astro-ph].\n- [9] D. Blais, C. Kiefer, and D. Polarski, Phys. Lett. B535 , 11 (2002), arXiv:astro-ph/0203520 [astro-ph].\n- [10] B. Carr, K. Kohri, Y. Sendouda, and J. Yokoyama, Rept. Prog. Phys. 84 , 116902 (2021), arXiv:2002.12778 [astro-ph.CO].\n- [11] B. Carr and F. Kuhnel, Ann. Rev. Nucl. Part. Sci. 70 , 355 (2020), arXiv:2006.02838 [astro-ph.CO].\n- [12] A. M. Green and B. J. Kavanagh, J. Phys. G 48 , 043001 (2021), arXiv:2007.10722 [astro-ph.CO].\n- [13] A. A. Starobinsky, JETP Lett. 30 , 682 (1979).\n- [14] A. A. Starobinsky, Phys. Lett. 91B , 99 (1980).\n- [15] D. Kazanas, Astrophys. J. 241 , L59 (1980).\n- [16] A. H. Guth, Phys. Rev. D23 , 347 (1981).\n- [17] K. Sato, Mon. Not. Roy. Astron. Soc. 195 , 467 (1981).\n- [18] V. F. Mukhanov and G. V. Chibisov, JETP Lett. 33 , 532 (1981), [Pisma Zh. Eksp. Teor. Fiz. 33, 549 (1981)].\n- [19] A. D. Linde, Phys. Lett. 108B , 389 (1982).\n- [20] A. Albrecht and P. J. Steinhardt, Phys. Rev. Lett. 48 , 1220 (1982).\n- [21] S. W. Hawking and I. G. Moss, Phys. Lett. 110B , 35 (1982).\n- [22] G. V. Chibisov and V. F. Mukhanov, Mon. Not. Roy. Astron. Soc. 200 , 535 (1982).\n- [23] S. W. Hawking, Phys. Lett. 115B , 295 (1982).\n- [24] A. H. Guth and S. Y. Pi, Phys. Rev. Lett. 49 , 1110  \n(1982).  \n- [25] A. A. Starobinsky, Phys. Lett. 117B , 175 (1982).\n- [26] M. Sasaki, Prog. Theor. Phys. 76 , 1036 (1986).\n- [27] V. F. Mukhanov, Sov. Phys. JETP 67 , 1297 (1988), [Zh. Eksp. Teor. Fiz. 94N7, 1 (1988)].\n- [28] V. Faraoni, Am. J. Phys. 69 , 372 (2001), arXiv:physics/0006030 [physics].\n- [29] W. H. Kinney, Phys. Rev. D72 , 023515 (2005), arXiv:gr-qc/0503017 [gr-qc].\n- [30] J. Martin, H. Motohashi, and T. Suyama, Phys. Rev. D87 , 023514 (2013), arXiv:1211.0083 [astro-ph.CO].\n- [31] J. Garcia-Bellido and E. Ruiz Morales, Phys. Dark Univ. 18 , 47 (2017), arXiv:1702.03901 [astro-ph.CO].\n- [32] J. M. Ezquiaga, J. Garcia-Bellido, and E. Ruiz Morales, Phys. Lett. B776 , 345 (2018), arXiv:1705.04861 [astroph.CO].\n- [33] K. Kannike, L. Marzola, M. Raidal, and H. Veermae, JCAP 1709 , 020 (2017), arXiv:1705.06225 [astroph.CO].\n- [34] C. Germani and T. Prokopec, Phys. Dark Univ. 18 , 6 (2017), arXiv:1706.04226 [astro-ph.CO].\n- [35] H. Motohashi and W. Hu, Phys. Rev. D96 , 063503 (2017), arXiv:1706.06784 [astro-ph.CO].\n- [36] K. Dimopoulos, Phys. Lett. B775 , 262 (2017), arXiv:1707.05644 [hep-ph].\n- [37] Y. Gong and Y. Gong, JCAP 1807 , 007 (2018), arXiv:1707.09578 [astro-ph.CO].\n- [38] G. Ballesteros and M. Taoso, Phys. Rev. D97 , 023501 (2018), arXiv:1709.05565 [hep-ph].\n- [39] M. P. Hertzberg and M. Yamada, Phys. Rev. D 97 , 083509 (2018), arXiv:1712.09750 [astro-ph.CO].\n- [40] C. Pattison, V. Vennin, H. Assadullahi, and D. Wands, JCAP 1808 , 048 (2018), arXiv:1806.09553 [astroph.CO].\n- [41] M. Biagetti, G. Franciolini, A. Kehagias, and A. Riotto, JCAP 07 , 032 (2018), arXiv:1804.07124 [astro-ph.CO].\n- [42] J. M. Ezquiaga and J. Garcia-Bellido, JCAP 08 , 018 (2018), arXiv:1805.06731 [astro-ph.CO].\n- [43] S. Rasanen and E. Tomberg, JCAP 01 , 038 (2019), arXiv:1810.12608 [astro-ph.CO].\n- [44] A. A. Starobinsky, Lect. Notes Phys. 246 , 107 (1986).\n- [45] M. Morikawa, Phys. Rev. D 42 , 1027 (1990).\n- [46] S. Habib, Phys. Rev. D 46 , 2408 (1992), arXiv:grqc/9208006.\n- [47] M. Mijic, Phys. Rev. D 49 , 6434 (1994), arXiv:grqc/9401030.\n- [48] A. A. Starobinsky and J. Yokoyama, Phys. Rev. D 50 , 6357 (1994), arXiv:astro-ph/9407016.\n- [49] C. Pattison, V. Vennin, H. Assadullahi, and D. Wands, JCAP 1710 , 046 (2017), arXiv:1707.00537 [hep-th].\n- [50] D. Cruces, C. Germani, and T. Prokopec, JCAP 03 , 048 (2019), arXiv:1807.09057 [gr-qc].\n- [51] H. Firouzjahi, A. Nassiri-Rad, and M. Noorbala, JCAP 1901 , 040 (2019), arXiv:1811.02175 [hep-th].\n- [52] C. Pattison, V. Vennin, H. Assadullahi, and D. Wands, JCAP 07 , 031 (2019), arXiv:1905.06300 [astro-ph.CO].\n- [53] J. M. Ezquiaga, J. Garcia-Bellido, and V. Vennin, JCAP 03 , 029 (2020), arXiv:1912.05399 [astro-ph.CO].\n- [54] G. Ballesteros, J. Rey, M. Taoso, and A. Urbano, JCAP 2008 , 043 (2020), arXiv:2006.14597 [astro-ph.CO].\n- [55] A. De and R. Mahbub, Phys. Rev. D 102 , 123509 (2020), arXiv:2010.12685 [astro-ph.CO].\n- [56] D. G. Figueroa, S. Raatikainen, S. Rasanen, and E. Tomberg, Phys. Rev. Lett. 127 , 101302 (2021),  \n- arXiv:2012.06551 [astro-ph.CO].\n- [57] G. Rigopoulos and A. Wilkins, JCAP 12 , 027 (2021), arXiv:2107.05317 [astro-ph.CO].\n- [58] C. Pattison, V. Vennin, D. Wands, and H. Assadullahi, JCAP 04 , 080 (2021), arXiv:2101.05741 [astro-ph.CO].\n- [59] D. G. Figueroa, S. Raatikainen, S. Rasanen, and E. Tomberg, JCAP 05 , 027 (2022), arXiv:2111.07437 [astro-ph.CO].\n- [60] E. Tomberg, JCAP 04 , 042 (2023), arXiv:2210.17441 [astro-ph.CO].\n- [61] E. Tomberg, Phys. Rev. D 108 , 043502 (2023), arXiv:2304.10903 [astro-ph.CO].\n- [62] S. S. Mishra, E. J. Copeland, and A. M. Green, JCAP 09 , 005 (2023), arXiv:2303.17375 [astro-ph.CO].\n- [63] M. Biagetti, V. De Luca, G. Franciolini, A. Kehagias, and A. Riotto, Phys. Lett. B 820 , 136602 (2021), arXiv:2105.07810 [astro-ph.CO].\n- [64] N. Kitajima, Y. Tada, S. Yokoyama, and C.-M. Yoo, JCAP 10 , 053 (2021), arXiv:2109.00791 [astro-ph.CO].\n- [65] S. Hooshangi, M. H. Namjoo, and M. Noorbala, Phys. Lett. B 834 , 137400 (2022), arXiv:2112.04520 [astroph.CO].\n- [66] Y.-F. Cai, X.-H. Ma, M. Sasaki, D.-G. Wang, and Z. Zhou, Phys. Lett. B 834 , 137461 (2022), arXiv:2112.13836 [astro-ph.CO].\n- [67] V. De Luca and A. Riotto, Phys. Lett. B 828 , 137035 (2022), arXiv:2201.09008 [astro-ph.CO].\n- [68] G. Ferrante, G. Franciolini, A. Iovino, Junior., and A. Urbano, Phys. Rev. D 107 , 043520 (2023), arXiv:2211.01728 [astro-ph.CO].\n- [69] A. D. Gow, H. Assadullahi, J. H. P. Jackson, K. Koyama, V. Vennin, and D. Wands, EPL 142 , 49001 (2023), arXiv:2211.08348 [astro-ph.CO].\n- [70] S. Pi and M. Sasaki, Phys. Rev. Lett. 131 , 011002 (2023), arXiv:2211.13932 [astro-ph.CO].\n- [71] S. Hooshangi, M. H. Namjoo, and M. Noorbala, JCAP 09 , 023 (2023), arXiv:2305.19257 [astro-ph.CO].\n- [72] T. Nakama, T. Harada, A. G. Polnarev, and J. Yokoyama, JCAP 01 , 037 (2014), arXiv:1310.3007 [gr-qc].\n- [73] S. Young, C. T. Byrnes, and M. Sasaki, JCAP 1407 , 045 (2014), arXiv:1405.7023 [gr-qc].\n- [74] C.-M. Yoo, T. Harada, J. Garriga, and K. Kohri, PTEP 2018 , 123E01 (2018), arXiv:1805.03946 [astro-ph.CO].\n- [75] M. Shibata and M. Sasaki, Phys. Rev. D 60 , 084002 (1999), arXiv:gr-qc/9905064.\n- [76] T. Harada, C.-M. Yoo, and K. Kohri, Phys. Rev. D 88 , 084051 (2013), [Erratum: Phys.Rev.D 89, 029903 (2014)], arXiv:1309.4201 [astro-ph.CO].\n- [77] T. Harada, C.-M. Yoo, T. Nakama, and Y. Koga, Phys. Rev. D 91 , 084057 (2015), arXiv:1503.03934 [gr-qc].\n- [78] I. Musco, Phys. Rev. D 100 , 123524 (2019), arXiv:1809.02127 [gr-qc].\n- [79] S. Young, I. Musco, and C. T. Byrnes, JCAP 11 , 012 (2019), arXiv:1904.00984 [astro-ph.CO].\n- [80] S. Young, Int. J. Mod. Phys. D 29 , 2030002 (2020), arXiv:1905.01230 [astro-ph.CO].\n- [81] A. Escriv'a, Phys. Dark Univ. 27 , 100466 (2020), arXiv:1907.13065 [gr-qc].\n- [82] A. Escriv'a, C. Germani, and R. K. Sheth, Phys. Rev. D 101 , 044022 (2020), arXiv:1907.13311 [gr-qc].\n- [83] V. Atal, J. Cid, A. Escriv'a, and J. Garriga, JCAP 05 , 022 (2020), arXiv:1908.11357 [astro-ph.CO].\n- [84] C.-M. Yoo, T. Harada, and H. Okawa, Phys. Rev. D\n- 102 , 043526 (2020), arXiv:2004.01042 [gr-qc].\n- [85] A. Escriv'a, C. Germani, and R. K. Sheth, JCAP 01 , 030 (2021), arXiv:2007.05564 [gr-qc].\n- [86] I. Musco, V. De Luca, G. Franciolini, and A. Riotto, Phys. Rev. D 103 , 063538 (2021), arXiv:2011.03014 [astro-ph.CO].\n- [87] A. Escriv'a, Universe 8 , 66 (2022), arXiv:2111.12693 [grqc].\n- [88] A. Escriv'a, Y. Tada, S. Yokoyama, and C.-M. Yoo, JCAP 05 , 012 (2022), arXiv:2202.01028 [astro-ph.CO].\n- [89] C.-M. Yoo, Galaxies 10 , 112 (2022), arXiv:2211.13512 [astro-ph.CO].\n- [90] T. Harada, C.-M. Yoo, and Y. Koga, Phys. Rev. D 108 , 043515 (2023), arXiv:2304.13284 [gr-qc].\n- [91] A. Escriv'a, V. Atal, and J. Garriga, JCAP 10 , 035 (2023), arXiv:2306.09990 [astro-ph.CO].\n- [92] A. Escriv'a and C.-M. Yoo, JCAP 04 , 048 (2024), arXiv:2310.16482 [gr-qc].\n- [93] J. M. Bardeen, J. R. Bond, N. Kaiser, and A. S. Szalay, Astrophys. J. 304 , 15 (1986).\n- [94] C. Germani and I. Musco, Phys. Rev. Lett. 122 , 141302 (2019), arXiv:1805.04087 [astro-ph.CO].\n- [95] V. Atal, J. Garriga, and A. Marcos-Caballero, JCAP 09 , 073 (2019), arXiv:1905.13202 [astro-ph.CO].\n- [96] A. Kalaja, N. Bellomo, N. Bartolo, D. Bertacca, S. Matarrese, I. Musco, A. Raccanelli, and L. Verde, JCAP 10 , 031 (2019), arXiv:1908.03596 [astro-ph.CO].\n- [97] M. Kawasaki and H. Nakatsuka, Phys. Rev. D 99 , 123501 (2019), arXiv:1903.02994 [astro-ph.CO].\n- [98] C. Germani and R. K. Sheth, Universe 9 , 421 (2023), arXiv:2308.02971 [astro-ph.CO].\n- [99] M. Sasaki and E. D. Stewart, Prog. Theor. Phys. 95 , 71 (1996), arXiv:astro-ph/9507001.\n- [100] M. Sasaki and T. Tanaka, Prog. Theor. Phys. 99 , 763 (1998), arXiv:gr-qc/9801017.\n- [101] D. Wands, K. A. Malik, D. H. Lyth, and A. R. Liddle, Phys. Rev. D62 , 043527 (2000), arXiv:astroph/0003278 [astro-ph].\n- [102] D. H. Lyth, K. A. Malik, and M. Sasaki, JCAP 0505 , 004 (2005), arXiv:astro-ph/0411220 [astro-ph].\n- [103] K. Ando, K. Inomata, and M. Kawasaki, Phys. Rev. D97 , 103528 (2018), arXiv:1802.06393 [astro-ph.CO].\n- [104] V. De Luca, A. Kehagias, and A. Riotto, Phys. Rev. D 108 , 063531 (2023), arXiv:2307.13633 [astro-ph.CO].\n- [105] T. Nakama, JCAP 10 , 040 (2014), arXiv:1408.0955 [grqc].\n- [106] A. Dekel, Astron. Astrophys. 101 , 79 (1981).\n- [107] Y. Akrami et al. (Planck), Astron. Astrophys. 641 , A10 (2020), arXiv:1807.06211 [astro-ph.CO].\n- [108] A. Karam, N. Koivunen, E. Tomberg, V. Vaskonen, and H. Veermae, JCAP 03 , 013 (2023), arXiv:2205.13540 [astro-ph.CO].\n- [109] K. Itˆo, Proceedings of the Imperial Academy 20 , 519 (1944).\n- [110] V. Vennin and A. A. Starobinsky, Eur. Phys. J. C 75 , 413 (2015), arXiv:1506.04732 [hep-th].\n- [111] Y. Tada and V. Vennin, JCAP 02 , 021 (2022), arXiv:2111.15280 [astro-ph.CO].\n- [112] M. Kopp, S. Hofmann, and J. Weller, Phys. Rev. D 83 , 124025 (2011), arXiv:1012.4369 [astro-ph.CO].\n- [113] A. D. Gow, C. T. Byrnes, P. S. Cole, and S. Young, JCAP 02 , 002 (2021), arXiv:2008.03289 [astro-ph.CO].\n- [114] A. Escriv'a and A. E. Romano, JCAP 05 , 066 (2021), arXiv:2103.03867 [gr-qc].  \n- [115] J. C. Niemeyer and K. Jedamzik, Phys. Rev. Lett. 80 , 5481 (1998), arXiv:astro-ph/9709072.\n- [116] J. C. Niemeyer and K. Jedamzik, Phys. Rev. D59 , 124013 (1999), arXiv:astro-ph/9901292 [astro-ph].\n- [117] I. Musco, J. C. Miller, and A. G. Polnarev, Class. Quant. Grav. 26 , 235001 (2009), arXiv:0811.1452 [grqc].\n- [118] I. Musco and J. C. Miller, Class. Quant. Grav. 30 , 145009 (2013), arXiv:1201.2379 [gr-qc].\n- [119] K. Inomata, M. Braglia, X. Chen, and S. Renaux-Petel, JCAP 04 , 011 (2023), [Erratum: JCAP 09, E01 (2023)], arXiv:2211.02586 [astro-ph.CO].\n- [120] J. Kristiano and J. Yokoyama, Phys. Rev. Lett. 132 , 221003 (2024), arXiv:2211.03395 [hep-th].\n- [121] A. Riotto, (2023), arXiv:2301.00599 [astro-ph.CO].\n- [122] S. Choudhury, M. R. Gangopadhyay, and M. Sami, (2023), arXiv:2301.10000 [astro-ph.CO].\n- [123] J. Kristiano and J. Yokoyama, Phys. Rev. D 109 , 103541 (2024), arXiv:2303.00341 [hep-th].\n- [124] A. Riotto, (2023), arXiv:2303.01727 [astro-ph.CO].\n- [125] H. Firouzjahi, JCAP 10 , 006 (2023), arXiv:2303.12025 [astro-ph.CO].\n- [126] H. Motohashi and Y. Tada, JCAP 08 , 069 (2023), arXiv:2303.16035 [astro-ph.CO].\n- [127] H. Firouzjahi and A. Riotto, JCAP 02 , 021 (2024), arXiv:2304.07801 [astro-ph.CO].\n- [128] G. Franciolini, A. Iovino, Junior., M. Taoso, and A. Urbano, Phys. Rev. D 109 , 123550 (2024), arXiv:2305.03491 [astro-ph.CO].\n- [129] G. Tasinato, Phys. Rev. D 108 , 043526 (2023), arXiv:2305.11568 [hep-th].\n- [130] S.-L. Cheng, D.-S. Lee, and K.-W. Ng, JCAP 03 , 008 (2024), arXiv:2305.16810 [astro-ph.CO].\n- [131] J. Fumagalli, (2023), arXiv:2305.19263 [astro-ph.CO].\n- [132] S. Maity, H. V. Ragavendra, S. K. Sethi, and L. Sriramkumar, JCAP 05 , 046 (2024), arXiv:2307.13636 [astro-ph.CO].\n- [133] Y. Tada, T. Terada, and J. Tokuda, JHEP 01 , 105 (2024), arXiv:2308.04732 [hep-th].\n- [134] H. Firouzjahi, Phys. Rev. D 109 , 043514 (2024), arXiv:2311.04080 [astro-ph.CO].\n- [135] M. W. Davies, L. Iacconi, and D. J. Mulryne, JCAP 04 , 050 (2024), arXiv:2312.05694 [astro-ph.CO].  \nFIG. 3. The curvature power spectrum. The vertical lines mark the peak scale k peak and the CMB pivot scale k ∗ = 0 . 05 Mpc -1 .  \n<!-- image -->", 'APPENDIX A: CURVATURE POWER SPECTRUM': 'The inflationary model studied in this paper is the asteroid mass case previously studied in [56, 59-61], built from Higgs inflation by adding quantum corrections to the plateau potential and tuning the CMB region by hand. The curvature power spectrum is depicted in figure 3, and it has a peak at k = 0 . 2 k PBH = 3 . 3 × 10 12 Mpc -1 . After the peak, the power spectrum declines as k -ϵ 2 , where ϵ 2 = 0 . 807 is the second slow-roll parameter.', 'APPENDIX B: DETAILS ON RESOLUTION': 'We studied the effects of the k -resolution and the r -resolution used in our simulations on the PBH abundance, both separately and together.  \nTo test the effect of the k -resolution, we considered four step lengths d ln k : 1 / 32, 1 / 64, 1 / 128, and 1 / 256. Reducing the step length appears to increase the number of peaks in C ( r ) at large values of r , shifting the P ( C ) peak to larger values of C . The difference in the PBH abundance between the largest and the smallest step length is about a factor of 3. The relative difference decreases with each doubling of the k -resolution (70%, 45% and 30%, respectively). This suggests that the abundance may eventually saturate with a small enough step length.  \nThe scaling of the number of peaks n ( p ) ∝ e -αp above the threshold C th = 0 . 4 varies with the resolution as α = 0 . 79, 0 . 83, 1 . 0, and 1 . 4, each value corresponding to halving our maximum resolution. Again, the changes decrease with increasing resolution. For the quantity P ( ¯ C ), the effects of varying the k -resolution are negligible.  \nIn our simulations, we considered r -values ranging from 10 -1 k -1 PBH to 10 4 k -1 PBH . To test the effects of resolution in r , we considered three cases: 1024, 2048, and 4096 r -points logarithmically distributed in the given range. We also considered extending the upper limit of the r -range by an order of magnitude to 10 5 k -1 PBH . At the lowest k -resolution of d ln k = 1 / 32, the resolution or the range extension of r have no quantifiable effect. However, higher k -resolution works in synergy with the increased r -resolution to enhance the peaks at large r , further pushing the P ( C ) peak towards larger values of C . Doubling the r -resolution and extending the range of r by an order of magnitude at the highest k -resolution d ln k = 1 / 256 increases the PBH abundance by 28%. As with the k -resolution, increasing the r -resolution and extending the range have negligible effect on P ( ¯ C ).  \nThe computational cost of the simulations is directly proportional to both the k - and r -resolution, so increasing both ramps up the amount of computational resources needed to produce sufficient statistics. In this Letter we have settled on using the highest tested k -resolution (d ln k = 1 / 256) with 1024 r -points in the range spanning from 10 -1 k -1 PBH to 10 4 k -1 PBH .', 'APPENDIX C: NUMERICAL ALGORITHM': 'Before the simulations begin, we pre-compute the power spectrum P ζ ( k ) for our potential. Then, for each realisation, we execute the following numerical algorithm (simplified for readability):', 'Algorithm 1: Single realisation': 'Set X < 0 = 0.  \nfor all modes k do  \nSolve ζ <k from (8).  \nGenerate Gaussian random variable ˆ ξ  \nk . Solve X <k +1 from (7).  \nfor all points r do  \nSolve ζ ( r ) from (2).  \nSolve rζ ′ ( r ) from (9).  \nSolve C ( r ) from (3).  \nSolve ¯ C ( r ) from (4).  \nCount the number of C ( r ) > 0 . 4, ¯ C ( r ) > 0 .  \n4 peaks. Store C max , ¯ C max , the corresponding r values, ζ <k PBH , and the number of peaks.'}