Primordial black holes from a curvaton scenario with strongly non-Gaussian perturbations

We investigate the production of primordial black holes (PBHs) in a mixed inflaton-curvaton scenario with a quadratic curvaton potential, assuming the curvaton is in de Sitter equilibrium during inflation with 〈χ〉 = 0. In this setup, the curvature perturbation sourced by the curvaton is strongly non-Gaussian, containing no leading Gaussian term. We show that for m 2/H 2 ≳ 0.3, the curvaton contribution to the spectrum of primordial perturbations on CMB scales can be kept negligible but on small scales the curvaton can source PBHs. In particular, PBHs in the asteroid mass range 10-16 M ⊙ ≲ M ≲ 10-10 M ⊙ with an abundance reaching F PBH = 1 can be produced when the inflationary Hubble scale H ≳ 1012 GeV and the curvaton decay occurs in the window from slightly before the electroweak transition to around the QCD transition.


Introduction
The study of primordial black holes (PBHs) formed via gravitational collapse of large primordial density fluctuations was initiated over 50 years ago [1,2].Already in [3] it was proposed that PBHs could form a cold dark matter component in the universe.The possibility that PBHs of mass M 10 −10 M would constitute all of the dark matter is already ruled out by constraints from lensing, dynamical effects, structure formation and gravitational waves [4].In the asteroid mass window, 10 −16 M M 10 −10 M the constraints are uncertain and it is possible that PBHs with masses in this range could constitute a dominant dark matter component [4].Even a subdominant PBH contribution to dark matter can have characteristic observational imprints, such as gravitational waves from PBH mergers testable with LIGO-Virgo-KAGRA data [5][6][7][8][9].Observational signals associated to PBHs are among the very few ways to probe small scale primordial perturbations and the process responsible for their generation.
In this work, we will study PBH formation in the curvaton setup with a quadratic potential, assuming χ = 0.This is the equilibrium configuration of a light spectator during de Sitter inflation [64], provided the potential is minimized at χ = 0.Even if one starts from an initial configuration with a non-zero χ , the distribution is rapidly driven towards the equilibrium during de Sitter inflation [64,65] 1 .For χ = 0, perturbations of the curvaton energy density obey a Gaussian squared distribution and are large, δρ 2  χ / ρ χ 2 = 2 [40].Consequently, the curvature perturbation component sourced by the curvaton has no leading Gaussian part and it can not be expanded in small perturbations.Therefore, one needs to consider the full non-linear and non-Gaussian solution for the curvature perturbation when investigating the PBH formation.The curvature perturbation ζ on CMB anisotropy scales is Gaussian to a high precision [67] and contributions sourced by the Gaussian squared δρ χ / ρ χ must be strongly suppressed on these scales.If the curvaton spectrum is sufficiently bluetilted, it can however give a dominant contribution to ζ on small scales k Mpc −1 where there are no constraints on non-Gaussianity.We focus on such a setup, investigating a mixed inflaton-curvaton scenario with a strongly blue tilted curvaton component which dominates the curvature perturbation on small scales and sources a fully non-Gaussian ζ with no leading Gaussian component.On the CMB scales, we require that the spectrum of ζ is dominated by a Gaussian inflaton component and contributions from the curvaton are suppressed to the 10 −11 level at the CMB pivot scale k * = 0.05 Mpc −1 .We note that PBHs in a partially related phenomenological setup with a Gaussian blue tilted spectator curvature perturbation component was investigated in [52].
We use the δN approach [68][69][70][71] to obtain a non-linear solution for the superhorizon scale curvature perturbation.We expand the curvaton field χ(x) in spherical harmonics and, following [61], truncate the expansion to the monopole when considering fluctuations relevant for PBHs.Within this approximation, we compute the probability distribution for C l (r) which determines the compaction function C(C l (r)).The compaction function C(r) equals the comoving gauge density contrast smoothed over a radius r.Using C c = 0.55 [72] as the PBH collapse threshold and modeling the PBH mass with the collapse parameters obtained in [73] for the power law spectrum, we explore the fraction of dark matter in PBHs f PBH and the mass distribution f (M ) as functions of the curvaton model parameters.We show that the curvaton scenario with the quadratic potential and the equilibrium configuration χ = 0 can lead to very efficient PBH production.In particular, we find that the scenario can produce asteroid mass PBHs, 10 −16 M M 10 −10 M , with f PBH = 1 when the inflationary Hubble scale H 10 12 GeV, the curvaton mass m 2 /H 2 0.3 and the curvaton decay occurs in the window from slightly before the electroweak transition to around the QCD transition.
The paper is organised as follows.In Sec. 2 we present the setup and in Sec. 3 we compute the probability distribution for the compaction function.In Sec. 4 we collect the expressions for the PBH abundance f PBH and the mass distribution f (M ).In Sec. 5 we present our main results and summarise the discussion in Sec. 6.

The mixed inflaton-curvaton setup
We investigate the PBH abundance in a mixed inflaton-curvaton scenario where the inflaton generates the Gaussian, nearly scale invariant curvature perturbations on CMB scales and curvaton-sourced perturbations dominate on small scales relevant for PBH formation.We assume the quadratic curvaton potential We further assume that the curvaton distribution during inflation follows the de Sitter equilibrium result with a vanishing mean value χ = 0 [64].For the quadratic potential, the curvaton χ is a Gaussian field and for m 2 /H 2 < 3/2 its spectrum at the end of inflation on superhorizon scales is given by the standard power-law expression [74] We denote the Hubble scale at the end of inflation by H ≡ H(t end ), and k end = a(t end )H(t end ) is the mode exiting the horizon at the end of inflation.Perturbations of the curvaton energy density δρ χ (x) obey Gaussian squared statistics, and the perturbations are large since δρ 2 χ (x) / ρ χ 2 = 2. Consequently, any contribution to the curvature perturbation sourced by δρ χ / ρ χ must be suppressed on the large scales k Mpc −1 probed by the CMB and LSS data [67].
Using the δN formalism [68][69][70][71], the superhorizon curvature perturbation ζ in the mixed scenario with inflaton sourced perturbations in the radiation component and the perturbed curvaton component obeys the non-linear equation [75] e 4ζ − Ω χ e 3ζχ e ζ + (Ω χ − 1) e 4ζr = 0 . (2.4) The individual curvature perturbations of the radiation ζ r and curvaton ζ χ fluids, and Ω χ are given in terms of spatially flat gauge quantities by the expressions, Here and in the rest of the text we use χ ≡ χ end to denote the curvaton field at the end of inflation.
The solution for the fourth order algebraic equation (2.4) can be written as [75] Here the only time dependent quantity is the curvaton density parameter Ω χ .Assuming the curvaton is subdominant at the onset of oscillations, which we define as H(t osc ) ≡ m, the density parameter for t t osc can be written as Here we used that ρ χ,osc ≈ (1/2)m 2 0.816χ 2 which follows from solving the curvaton equation of motion in a radiation dominated universe.We approximate the curvaton decay as an instant process at ) is obtained by evaluating Eq. (2.7) at the moment t dec .We assume the perturbation of the radiation component ζ r is entirely sourced by the inflaton and uncorrelated with the curvaton, ζ r ζ χ = 0. We further assume ζ r obeys Gaussian statistics with the power law spectrum where k * = 0.05 Mpc −1 .As explained above, we want to realise a scenario where the Gaussian inflaton sourced perturbations ζ r dominate the two point function of ζ on large scales and generate the observed CMB spectrum with P ζ (k * ) = A s = 2.10 × 10 −9 , n s = 0.965 [76].Correspondingly, contributions from ζ χ to the spectrum of ζ must be sufficiently suppressed.
To this end, we require which, as we show in Appendix B, allows us to obtain the observed CMB spectrum with A r ∼ 10 −9 and |n r − 1| ∼ 0.01, assuming Ω χ 0.9.Note that the underlying computation is somewhat non-trivial as ζ is a strongly non-linear function of the non-Gaussian ζ χ .The condition (2.10) constrains the curvaton mass from below as we show in detail in Appendix A, and implies a strongly blue tilted spectrum both for the curvaton field Eq. (2.2) and for ζ χ .On small scales, k Mpc −1 , the connected correlators of ζ are dominated by ζ χ which, as we show below, can lead to efficient formation of PBHs.For the maximal inflationary Hubble scale H ≈ 5.2 × 10 13 GeV consistent with the observational bound on the tensor-to-scalar ratio r T < 0.044 [77], Eq. (2.10) implies m 2 /H 2 0.29, assuming instant transition from inflation to radiation domination.The lower bound on m 2 /H 2 slowly grows as a function of decreasing H.For example, for H = 1.6 × 10 11 GeV, Eq. (2.10) implies m 2 /H 2 0.31, see Fig. 6 in Appendix B.

Probability distribution of density fluctuations
The central quantity in determining if an overdense region collapses into a primordial black hole is the compaction function C(r) [78][79][80] which for spherical overdensities equals the comoving density contrast coarse-grained with a top hat window function over a spherical volume of comoving radius r.Denoting the background equation of state by w = p / ρ , the expression for C(r) can be written as [81,82] where and C l (r) is determined by the curvature perturbation ζ as Here the prime denotes a derivative with respect to r.
Around the high density peaks relevant for the PBH formation, spherical symmetry can be expected to be a reasonable first approximation.We implement the approximation following [61] by expanding the Gaussian field χ in spherical harmonics and retaining only the leading monopole term of the expansion χ(r) = dk (2π) 3 j 0 (kr)χ k . (3.5) Here j 0 (z) = sin(z)/z.Since Eq. (3.5) is a linear map from χ(x), the field χ(r) and its derivative χ (r) are Gaussian fields with the joint probability distribution given by The components of the covariance matrix depend on r and can be written as where a prime denotes a derivative with respect to r and P(k) is the power spectrum of the full curvaton field χ(x) In our case P(k) is given by Eq. (2.2).We denote the variance of the full field χ(x) by σ 2 , We proceed to substitute χ(x) → χ(r) in Eqs.(2.6) and (2.7) in places where the field χ(x) appears with no brackets, but evaluate the background quantities that depend on χ 2 in Eq. (2.7) using the variance of the full field (3.11).Using χ(r) 2 < χ(x) 2 would give a higher probability for larger ζ χ values and hence enhance the PBH abundance but it is hard to quantify to what extent this is a spurious effect of the monopole truncation.We will therefore conservatively use χ(x) 2 in the background quantities.
In this setup, Eq. (3.3) takes the form where ∂ χ ζ ≡ ∂ζ ∂χ is obtained by differentiating Eq. (2.7).The probability distribution of C l is given by Carrying out the integral over χ , we obtain where the remaining integral needs to be computed numerically.

Expression for the PBH abundance
The mass of a PBH formed by a collapsing region with compaction C can be approximated by obtained by fitting to numerical simulations [72,[83][84][85][86][87][88].Here M H = 4π/3H −3 ρ is the mass within a Hubble volume at the collapse time, γ depends on the equation of state, and K and the collapse threshold C c depend on both the equation of state and the shape of the collapsing overdensity.In a radiation dominated universe γ ≈ 0.36.For monochromatic PBHs formed during radiation domination from a Gaussian ζ with the spectrum There are no existing numerical collapse simulations corresponding to our case with a fully non-Gaussian ζ on small scales.We adopt a phenomenological approach and set the collapse parameters equal to the Gaussian power-law results in radiation domination, γ = 0.36, C c = 0.55, and K = 4. Therefore, we will compute the PBH mass using We evaluate M H (r) at the horizon entry of the smoothing scale a(t r )H(t r )r = 1, In our setup, the curvaton contribution to the energy density will in general not be negligible at t r and the universe is therefore not fully radiation dominated.However, decreasing the pressure decreases C c and using the threshold C c for radiation domination we should be estimating the PBH abundance from below.In any case, our results should be regarded as order of magnitude estimates both due to the use of Eq. (4.2) and due to the monopole truncation Eq. (3.5).
The probability that a spherical overdensity coarse grained over the comoving radius r collapses into a PBH upon horizon entry equals the probability that C exceeds the threshold C c .The contribution of the PBHs to the total energy density at the collapse time can then be written as where is given by Eq. (3.14).The upper limit C = f (w)/2 is the largest fluctuation amplitude that forms a type I overdensity for which the areal radius is a monotonic function of the coordinate r [91].
The fraction of the present day dark matter energy density constituted by the PBHs reads where k eq = (aH) eq at the matter radiation equality, and we have omitted the O(1) factor (g * (t r )/g * (t eq )) −1/6 from the effective number of relativistic degrees of freedom.This can be recast as with the mass distribution function f (M ) given by where ), as obtained from Eq. (4.2).

Results
It is straightforward to compute the variances Eq. (3.7) using Eq.(2.2) and numerically perform the integral in Eq. (3.14) to find the probability distribution P C l (C l , r).Using Eqs.(4.4) and (4.5) we then get f PBH (r) as a function of the coarse-graining scale r.GeV, H = 8.9 × 10 12 GeV, and m 2 /H 2 = 0.31.
Figure 1 illustrates the typical shape of the function f PBH (r) and the curvaton density parameter Ω χ at the horizon crossing of r, which enters in the computation via Eq.(3.14).The abundance f PBH (r) has a clearly peaked structure although the curvaton spectrum Eq. (2.2) is of pure power law form.This is a generic feature in the setup and it arises from an interplay of two opposite effects.First, increasing the coarse-graining scale r makes the variance σ 2 χχ (r) smaller, see Eq. (3.7), and therefore suppresses the probability of large χ(r) values that can source PBHs.Second, the curvature perturbation ζ and Ω χ keep growing in time until the curvaton decay at H(t dec ) = Γ.For t r < (t dec ), increasing r corresponds to later horizon crossing times t r , making ζ(t r ) larger and enhancing the probability for PBH formation.This effect dominates to the left of the f PBH (r) peak in Fig. 1, and the peak corresponds to t r = t dec .To the right of the peak, ζ stays constant as t r > t dec .In this region, increasing r only acts to decrease σ 2 χχ (r) and therefore f PBH (r) starts to decrease.In the following, we will choose the coarse-graining scale r equal to the peak scale by setting r = r dec ≡ 1/(a(t dec )H(t dec )), and define f PBH = f PBH (r dec ).(5.1) Black hole evaporation sets very strong constraints for PBHs with M 10 −16 M .Constraints for M 10 −10 M PBHs from various different systems range downwards from f PBH = O(10 −2 ) depending on the mass [4].In the asteroid mass window 10 −16 M M 10 −10 M , the constraints are subject to significant uncertainties and it can be possible to have f PBH = 1 in this window [4].Interestingly, asteroid mass PBHs can be efficiently produced in the curvaton scenario.This is demonstrated in Fig. 4 which shows the PBH mass spectra f (M ) computed using Eq.(4.7) for Γ = 5.5 × 10 −18 GeV, Γ = 6.5 × 10 −16 GeV and Γ = 5.5 × 10 −14 GeV, with m 2 /H 2 = 0.31, and H chosen in each case such that f PBH ≈ 0.10.The corresponding curvaton decay temperatures are T dec ∼ 2 GeV for Γ = 5.5 × 10 −18 GeV, T dec ∼ 20 GeV for Γ = 6.5 × 10 −16 GeV, and T dec ∼ 200 GeV for Γ = 5.5 × 10 −14 GeV, approximating the curvaton decay into radiation and the thermalisation of the decay products as instant processes, and using the Standard Model g * (T ).The respective mass spectra in Fig. 4 are peaked at M ∼ 10 −11 M , M ∼ 10 −13 M and M ∼ 10 −15 M .In all three cases, the PBH abundance can be increased by slightly increasing H, see Fig.  Figures 2, 3 and 4 represent the main results of this work.In particular, they show that the curvaton scenario can generate a significant dark matter fraction consisting of asteroid mass scale PBHs when the inflationary Hubble scale is large enough H 10 12 GeV, the curvaton mass parameter m 2 /H 2 0.3 and the decay rate falls in the window 10 −18 GeV Γ 10 −13 GeV, corresponding to decay temperatures from slightly above the electroweak transition to around the QCD transition scale.If the curvaton decay occurs earlier, 10 −13 GeV Γ 10 −8 GeV, the scenario generates PBHs with M 10 −16 M and, according to the dependencies shown in Figs. 2 and 3, the observational constraints on f PBH constrains the viable parameter space for H from above and for m2 /H 2 from below.For Γ 10 −8 GeV there are no constraints as the the PBH abundance is exponentially suppressed for any H 5.2 × 10 13 GeV compatible with the observational bound on the tensor-to-scalar ratio r T < 0.044 [77], and m 2 /H 2 in the range compatible with Eq. (2.10) 2 .For Γ 10 −18 GeV, the scenario generates PBHs with M 10 −10 M but in this region our numerical integration of Eq. (3.14) starts to become inaccurate for configurations leading to f PBH 0.01.We are therefore not able to perform a detailed study of the Γ 10 −18 GeV region in this work.
Finally, Fig. 5 depicts the probability distribution P C l (C l , r dec ) given by Eq. (3.14) for the same choice of parameters as in Fig. 1, and in the third panel of Fig. 4. For comparison, we also show a Gaussian distribution with the variance equal to the variance computed from the full distribution for this set of parameters, dC l C 2 l P C l (C l , r dec ) ≈ 5.1×10 −4 .The full distribution deviates significantly from the Gaussian case and decays much slower as a function of C l .The slowly decaying tail of P C l (C l , r dec ) is essential for the PBH formation in our setup.We recall, that the fully non-Gaussian form of P C l (C l , r dec ) follows from the vanishing mean of the curvaton field χ = 0 which in turn is the equilibrium configuration during inflation for the χ 2 potential.The curvature perturbation component ζ χ ∝ lnρ χ / ρ χ = lnχ 2 / χ 2 has no leading Gaussian term and there is no suppression for fluctuations χ 2 / χ 2 ∼ 1. Together with the non-linear form of Eq. (2.7), this gives rise to the fully non-Gaussian distribution of C l seen in Fig. 5.

Summary
In this work we have investigated the PBH production in the mixed inflaton-curvaton scenario with a quadratic curvaton potential and a strongly blue-tilted curvaton spectrum, assuming the curvaton is in the de Sitter equilibrium χ = 0 during inflation.We require that the inflaton sourced Gaussian component dominates the spectrum of the curvature perturbation ζ on CMB scales and the curvaton sourced contribution is suppressed below the 2 × 10 −11 level at the pivot scale k * = 0.05 Mpc −1 .This constrains the curvaton mass from below, for example for the inflationary Hubble scale H = 10 13 GeV, the curvaton mass must satisfy m 2 /H 2 0.3.On small scales, however, the curvaton sourced component of ζ can take large values leading to PBH formation.
A key feature in this setup, is that the curvaton sourced part of ζ is fully non-Gaussian.It contains no leading Gaussian term, unlike in scenarios with a non-vanishing curvaton mean value | χ | H (a specific initial condition during inflation) which have been studied previously in the PBH context e.g. in [45][46][47][48][49].We use the δN formalism to obtain the full non-linear solution for ζ as a function of χ and, following [61], include only the monopole term of the spherical harmonics series of χ(x) when considering fluctuations relevant for PBHs.With this approximation, we compute the probability distribution for C l (r) which determines the compaction function C(C l (r)), and where r is the coarse-graining scale of perturbations.Comparing to a fiducial Gaussian distribution with the same variance C 2 l , we find that the full distribution P l (C l , r) decreases exponentially slower as a function of C l .We use C c = 0.55 as the threshold for the PBH collapse [73] and model the PBH mass with parameters obtained for the power-law spectrum in [72,73].
We find that the curvaton scenario with χ = 0 can lead to very efficient production of PBHs.In particular, the setup can generate asteroid mass PBHs, 10 −16 M M 10 −10 M , with an abundance equal to the observed dark matter abundance, f PBH = 1, when the Hubble scale at the end of inflation H 10 12 GeV, the curvaton mass m 2 /H 2 0.3 and the curvaton decay occurs in the window from slightly before the electroweak transition to around the QCD transition.If the curvaton decays before or after the aforementioned window, PBHs with masses below or above the asteroid mass range can be generated, respectively.The PBH abundance depends sensitively on H, m 2 /H 2 and the curvaton decay rate Γ, and the observational bounds on f PBH imply non-trivial constraints on these parameters when Γ 10 −8 GeV.For larger values of Γ the PBH production in the setup is exponentially suppressed.
The main uncertainties in our results arise from the monopole truncation Eq. (3.5) and the phenomenological choice of collapse parameters in Eq. (4.2).Changing the collapse parameter values would affect f PBH predicted for a fixed set of curvaton parameters.However, even order of magnitude changes of f PBH are compensated by just slight changes of H, m 2 /H 2 and Γ, as can be seen in Figs. 2 and 3.The error caused by the use of Eq. (3.5) is harder to quantify but we expect that the very efficient PBH production from the strongly non-Gaussian ζ is a robust conclusion.

A The spectrum of ζ χ
We use the stochastic formalism [64] and the spectral expansion method [92][93][94] to compute the infrared spectrum of ζ χ in our setup where χ = 0 and Eq.(2.6) cannot be expanded in small perturbations around a mean field.For the quadratic curvaton potential (2.1), the joint equal time two-point distribution of χ(t, r) in de Sitter equilibrium is given by the spectral sum [92][93][94] where and H n (x) = (−1) n e x 2 /2 d n dx n e −x 2 /2 are the Hermite polynomials.The eigenfunctions Using that χ 2 = 3H 4 /(8π 2 m 2 ), we have from Eq. (2.6) ζ χ = (1/3)ln(8π3 x 2 /3), and the connected part of its two-point function can be written as Truncating the series at the leading order we get The corresponding power spectrum is given by The existence of the Fourier transform requires 4m 2 < 9H 2 which we assume here.Assuming an instant transition from de Sitter inflation to radiation dominated epoch, and approximating the universe as radiation dominated until the curvaton decay  6 where the shaded orange area marks the region where Eq. (2.10) holds.GeV consistent with the observational bound on the tensor-to-scalar ratio r T < 0.044 [77].The shaded orange area between the solid and dashed lines depicts the phenomenologically viable region where P ζχ (k * ) < 2 × 10 −11 .
We will also need the two point functions of the powers ζ χ in the analysis below.These can be directly computed using the spectral expansion expression .
Truncating the series again at the leading order we obtain

. 5 )
Both ζ r and ζ χ are separately conserved, i.e. constant in time.For the quadratic potential Eq. (2.1), the curvaton component ζ χ can be written in terms of the field χ as

Figure 1 .
Figure 1.Left panel: The fraction of the dark matter abundance in PBHs f PBH as function of the coarse-graining scale r.Right panel: The curvaton density parameter Ω χ evaluated at a(t r )H(t r )r = 1 shown as a function of the coarse-graining scale r.Both panels are computed setting Γ = 5.5 × 10 −14GeV, H = 8.9 × 10 12 GeV, and m 2 /H 2 = 0.31.

Figure 2 .
Figure 2. PBH abundance f PBH as a function of the curvaton decay rate Γ and the inflationary Hubble scale H. Computed for m 2 /H 2 = 0.31 (left) and 0.35 (right).

Figures 2
Figures 2 and 3 illustrate the behaviour of f PBH as a function of the decay rate Γ, the inflationary Hubble scale H, and the curvaton mass parameter m 2 /H 2 .Varying the decay rate Γ alters the curvaton density parameter Ω χ at t r = t dec .As seen in Fig.1, f PBH is a strongly dependent function of Ω χ .Increasing Γ moves the decay to earlier times and decreases Ω χ for fixed H, m 2 /H 2 values.This explains the strong decrease of f PBH as a function of Γ in Figs.2 and 3.The figures further show that, for a fixed Γ, increasing m 2 /H 2 or decreasing H causes f PBH to decrease.The former is because larger m 2 /H 2 makes the curvaton spectrum Eq. (2.2) more blue-tilted, decreasing the variance σ 2 χχ (r), Eq. (3.7), and therefore suppressing the probability for large χ(r) values.The latter is because the mean curvaton energy at the end of inflation ρ χ ∝ m 2 χ 2 ∝ H 4[64], and decreasing H therefore decreases Ω χ for fixed Γ and m 2 /H 2 values.

Figure 5 .
Figure 5.The full probability distribution function P C l (C l , r dec ) and a Gaussian distribution with the same variance C 2 l ≈ 5.1 × 10 −4 .Parameters chosen as in Fig. 1.

Figure 7 .
Figure 7.The solid line shows A/A r defined in Eq. (B.7) as a function of the curvaton density parameter Ω χ .The dashed line shows the corresponding quantity computed in linear perturbation theory.