Supermassive Primordial Black Holes From Inflation

There is controversy surrounding the origin and evolution of our universe's largest supermassive black holes (SMBHs). In this study, we consider the possibility that some of these black holes formed from the direct collapse of primordial density perturbations. Since the mass of a primordial black hole is limited by the size of the cosmological horizon at the time of collapse, these SMBHs must form rather late, and are naively in conflict with constraints from CMB spectral distortions. These limits can be avoided, however, if the distribution of primordial curvature perturbations is highly non-Gaussian. After quantifying the departure from Gaussianity needed to evade these bounds, we explore a model of multi-field inflation -- a non-minimal, self-interacting curvaton model -- which has all the necessary ingredients to yield such dramatic non-Gaussianities. We leave the detailed model building and numerics to a future study, however, as our goal is to highlight the challenges associated with forming SMBHs from direct collapse and to identify features that a successful model would need to have. This study is particularly timely in light of recent observations of high-redshift massive galaxy candidates by the James Webb Space Telescope as well as evidence from the NANOGrav experiment for a stochastic gravitational wave background consistent with SMBH mergers.


Introduction
Supermassive black holes (SMBHs) are ubiquitous in our universe, being present in the centers of nearly all massive galaxies.Quasars, 1 powered by black holes with with masses M ∼ 10 8−10 M ⊙ , are also found in large numbers in the high-redshift universe.At present, over 170 quasars have been observed at z > 6, with the most distant at z = 7.54, and several hundred others at z = 5 − 6 [2][3][4][5][6][7][8][9][10].Fig. 1 shows the quasar abundance as a function of redshift over several ranges of black hole mass.While this magnitude limited quasar catalog is believed to be nearly complete out to z ∼ 5, only a small fraction of SMBHs are, in fact, quasars.A more complete census of SMBHs is possible in the local universe where one finds Ω SMBH (z = 0) ∼ 10 −6 [11].This contrasts with the peak quasar mass density of Ω quasar (z = 2) ∼ 10 −8 .While specific models for SMBH population evolution have been proposed [12][13][14], the limited available data leaves a great deal of uncertainty.
It is usually assumed that SMBHs grow over time from relatively low-mass seeds (possibly the remnants of Population III stars [15]) through the process of accretion.The rate of 1 Quasars are the most luminous active galactic nuclei (AGN).In this paper, quasar refers to an AGN that is sufficiently luminous to appear in a quasar catalog such as SDSS DR7 [1].While this definition is redshift-dependent, the most luminous quasars should be consistently present in the catalog up to the DR7 redshift limit of 5.In the text, the phrase "quasar mass" refers to the mass of the black hole that powers the quasar.
Figure 1.Estimates of the quasar black hole comoving mass density in units of the z = 0 critical density as a function of redshift, plotted for the four different mass classes indicated.Black hole masses and redshifts are taken from the DR7 SDSS quasar survey [1], which covers 20% of the sky and the redshift range 0 < z < 5.The results shown represent the moving average of 100 quasars sorted by redshift.

mass accretion is Eddington limited to
where is the Salpeter time [16], σ T = 8π α 2 /(3m 2 e ) is the Thomson cross section, m p is the proton mass, and ε is the radiative efficiency.At the Eddington-limited rate, a 10 2 M ⊙ black hole seed would require ∼ 0.8 Gyr to grow to M ∼ 10 10 M ⊙ .Since z ∼ 6 − 7 corresponds to only ∼ 0.7 − 0.9 Gyr after the Big Bang, such a scenario would require these large high-redshift black holes to have grown at high accretion rates almost continuously throughout the first Gyr of our universe's history.Interestingly, there seems to be quite a significant population of quasars at z ∼ 6 to 7 in this mass range [4,[17][18][19].
The continuous accretion required to explain this rapid growth contrasts with the intermittent accretion observed of SMBHs at lower redshifts.Further, from Fig. 1 we see that the most massive M BH ≳ 10 10 M ⊙ population has remained approximately constant since at least z ∼ 5 [20].This would require any growth in the number of SMBHs to be balanced by a decrease in the fraction of SMBHs that are actively quasars.From this perspective, it is surprising that so many highly-massive quasars have been observed at such high redshifts [21][22][23][24][25].These observations prompt two intriguing questions: 1.If these quasars grew from small black hole seeds, how did they come to be so massive on such a short timescale?
2. Why did their growth rate dramatically slow down during the subsequent 13 Gyr?
Various solutions to the first question have been put forth, including an enhanced role of mergers as well as the possibility of super-Eddington accretion [26].However frequent mergers would require more heavily clustered populations in the early universe, and further can have both positive and negative impacts on SMBH growth, since they can also kick SMBHs out of the material-rich centers of galaxies.A certain degree of super-Eddington growth is expected to occur in high-redshift galaxies containing large reservoirs of gas and efficient angular momentum transport due to turbulence.However as increased accretion produces powerful jets and outflows which drive material away, it is uncertain how long it can be sustained.Feedback effects from transient periods of super-Eddington growth are actually expected to impede SMBH growth within a few Myr [27].
Regarding the second question, it has been suggested that the suppressed growth of the most massive black holes after the first Gyr could perhaps be attributed to galaxy-scale feedback.This, however, would require the M -σ relation2 to evolve with redshift and for the quasar luminosity function to steepen at the highest values [28,29].Alternatively, it has been proposed that there might be a maximum mass that black holes can reach through accretion, resulting from the fragmentation of the accretion disks that could have otherwise facilitated rapid black hole growth [30].Despite these suggestions, there remain many open questions concerning the origin and evolution of our universe's most massive black holes.
In this study, we take these questions to motivate another possibility: that our universe's most massive black holes did not acquire most of their mass through accretion, but are instead predominantly primordial in origin 3 .Unlike the small primordial black holes (PBH) typically considered as dark matter candidates, these massive objects would have formed at "late" times, as governed by the size of the cosmological horizon.During the radiation dominated era, the horizon contains the following amount of energy: where R H = H −1 is the size of the horizon, ρ = π 2 g ⋆ (T )T4 /30 is the radiation density, and g ⋆ is the number of relativistic species at temperature T .When a sufficiently large density fluctuation collapses to form a PBH, 4 the mass of the resulting black hole is typically an order one fraction of the horizon mass, M BH ≃ γ M H , where γ ∼ 0.2 quantifies the efficiency of collapse [38].From Eq. (1.3), we conclude that the PBHs in the mass range of interest here form after the start of Big Bang nucleosynthesis (T ∼ MeV corresponds to M H ∼ 10 5 M ⊙ ) but well before the onset of matter domination or recombination (T ∼ eV corresponds to Many inflationary models enable PBHs to form efficiently [39], including those in which the inflaton undergoes a period of ultra-slow-roll [40][41][42][43], or whose potential features localized bumps, dips, or steps on small scales [44][45][46].More generically, a local enhancement of the power spectrum P ζ requires a deviation from slow-roll evolution (see for example, Refs.[47][48][49][50][51][52][53][54][55][56][57][58][59][60][61]).Regardless of the mechanism, PBH formation requires a significant enhancement in the amplitude of the primordial power spectrum.Assuming Gaussian statistics, this amplitude must be P ζ (k BH ) ∼ 10 −2 , seven orders of magnitude larger than observed on large scales, [62].Such an amplification inevitably leads to large CMB spectral distortions.In particular, for scales k BH < ∼ 10 5 Mpc −1 , corresponding to black holes in the mass range M BH ≳ 10 3 M ⊙ , the predicted spectral distortions are in strong conflict with COBE/FIRAS measurements [63,64].
In principle, PBHs can form from smaller peaks in the power spectrum if the tail of the ζ distribution is sufficiently non-Gaussian.Of course, observations on large scales and measurements of the non-linearity parameters f NL and g NL seem to indicate that ζ is very nearly Gaussian on CMB scales [65,66].However the scales on which PBHs are formed can be disconnected from CMB scales, and so these bounds need not apply.
In this paper, we quantify the degree of non-Gaussianity that would be required to viably produce primordial SMBHs and present a model of inflation that contains the necessary features.We begin in Sec. 2 by reviewing the calculation of the PBH abundance in the Press-Schechter formalism.In Sec. 3, we examine how measurements of spectral distortions constrain the primordial power spectrum, and demonstrate that the appreciable formation of SMBHs from Gaussian density fluctuations is excluded on the basis of these constraints.Sec. 4 considers departures from Gaussianity, and quantifies the heaviness of the distribution needed to evade these bounds.In Sec. 5, we review the calculation of the curvature perturbation and its statistics in the standard curvaton scenario, which has been shown to be capable of producing appreciable non-Gaussianity.This model, however, cannot viably produce a significant abundance of primordial SMBHs, so in Sec.6 we introduce self-interactions to augment the non-Gaussianity and speculate on the maximum PBH abundance in this model.We conclude in Sec.7 with a discussion of our results and some comments on directions for future investigations.

Primordial Black Holes from Gaussian Perturbations
The initial abundance of PBHs is quantified in terms of their mass fraction at the time of formation: For PBHs that formed during radiation domination, this quantity is related to the fractional energy density in black holes today, , where z f is the redshift at the time of black hole formation.Using Eq. (1.3) and entropy conservation, this can be written more conveniently in terms of the initial PBH mass as: where T f is the temperature at the time of black hole formation.
In the standard5 treatment based on the Press-Schechter formalism [72], β is computed by integrating the probability distribution P δ for the coarse-grained density contrast δ = δρ/ρ over all values greater than the critical threshold for collapse, δ c : (2. 3) The factor of 2 is customarily introduced6 to compensate for the undercounting that otherwise arises [72].This prescription has the benefit of having an intuitive interpretation -we are summing over the fraction of regions with a sufficient overdensity to collapse to form a PBHand a simple implementation.One issue concerns the uncertainty7 of the collapse threshold, which should in principle depend on the overdensity profile 8 .Assuming a spherical profile, we expect a perturbation to be able to collapse during radiation domination if the size of the overdensity at maximum expansion exceeds the Jeans length.This led to Carr's original estimate of δ c = c 2 s = 1/3 [38], where c s is the sound speed of density perturbations.A more careful treatment that is applicable for an arbitrary equation of state9 , w, finds [81]: which has been shown to faithfully replicate the results of numerical simulations and yields δ c = 0.414 during radiation domination.We adopt this value throughout our analysis.
While the coarse-grained density contrast is the proper object to consider when computing the probability of PBH formation, it is often convenient to work directly with the curvature perturbation ζ, since its statistics are more easily computed from underlying inflationary models.When large perturbations exist only on one scale, as is the case for the sharply peaked power spectra we consider, only a minimal amount of error is incurred by making this approximation.On superhorizon scales, ζ is related to the density contrast field by [82] (2.5) Working to linear10 order in radiation domination, this simplifies to the Fourier space relation: which implies that their power spectra, defined for arbitrary Fourier variable f k via the 2-point function as , are related as Note that at the time of horizon crossing k ≃ aH, when a perturbation can collapse to form a black hole, 11 the density contrast and curvature perturbation are linearly related, δ ≃ 4 9 ζ.We can then assume that peaks in δ also correspond to peaks in ζ, and work directly with the curvature perturbation.In terms of ζ, the initial abundance in the Press-Schechter formalism reads: where the collapse threshold ζ c ≃ 9 4 δ c follows from σ 2 δ ≃ 16 81 σ 2 ζ , since P ∼ σ 2 for a sharply peaked spectrum.More precisely, the variance of ζ smoothed on the scale R ≃ (aH) −1 ≃ k −1 can be computed from the power spectrum as [84] where W (k, R) is the Fourier transform of the (real space) window function used to coarsegrain δ.It is unclear what functional form for the window function most accurately reproduces the actual relation between the PBH abundance and power spectrum, but popular choices in the literature include the (volume normalized) Gaussian, as well as real and k-space top hats.
For simplicity, we use the former, W (k, R) = exp −k 2 R 2 /2 , but see [85,86] for a discussion of the resultant uncertainties.
Consider the case of a Gaussian distribution P ζ = P G ζ , where (2.10) In this case the mass fraction at formation evaluates to where erfc is the complimentary error function and the second approximation holds for ζ c ≫ σ ζ , which is generically true for all cases of physical interest.Since the threshold is fixed, the initial abundance is determined solely by the variance of the power spectrum.In order to have β = 10 −20 , we see we need σ 2 ζ ≃ 0.01, corresponding to a peak in the power spectrum of P ζ ∼ 0.01, which is 7 orders of magnitude greater than the value measured on CMB scales.Since the amplification of P ζ needed for PBH formation depends on β only logarithmically, this degree of enhancement is a generic requirement for any non-vanishing initial abundance.

Spectral Distortions
The small scales relevant for primordial SMBH formation are well below the angular resolution probed by current CMB measurements.Nevertheless, inhomogeneities on these scales will generate isotropic deviations from the usual blackbody spectrum [87][88][89][90][91].These deviations are known as spectral distortions and in this context it is useful to distinguish between three characteristic redshift intervals: • Thermalization era (z > 2 × 10 6 ): At high redshifts, Compton scattering γe → γe, double Compton scattering γe → γγe, and Bremsstrahlung ep → epγ maintain a blackbody spectrum for the photons, and spectral distortions are exponentially suppressed.
• µ-era (2 × 10 5 < z < 2 × 10 6 ): During this era, photon number changing processes, double Compton scattering and Bremsstrahlung, become ineffective at maintaining a blackbody spectrum.Compton scattering, however, continues to redistribute photon energies to maintain a Bose-Einstein distribution, parameterized by both a temperature, T , and a chemical potential, µ.A µ-distortion refers to a Bose-Einstein distribution with µ ̸ = 0.
• y-era (z < 2 × 10 5 ): Compton scattering becomes ineffective at redistributing photon energies during this era, so there are no processes to maintain a Bose-Einstein distribution.Spectral distortions that are generated at these redshifts are characterized by a departure from an equilibrium distribution, and often yield so-called y-distortions. 12 spectrum with a (positive) y-distortion can be expressed as an average of blackbodies with slightly different temperatures [93][94][95][96].An average of blackbodies with a mean temperature T and variance T 2 ∆ will (for ∆ ≪ 1) be a y-distorted blackbody characterized by the temperature T = T (1+∆ 2 ) and the parameter y = ∆ 2 /2.y-type distortions can be generated through a variety of mechanisms, including the Compton scattering of CMB photons with a population of electrons with a different temperature.In this paper, we will be interested in y-distortions that are generated through photon diffusion.Limits on such spectral distortions allow us to constrain the amplitude of inhomogeneities on very small scales [97].
µ-and y-type spectral distortions are traditionally quantified in terms of the parameters µ and y, which are related to the fractional increase in energy per photon (relative to a blackbody spectrum with the same number density of photons) [87,91]: By introducing k-space window functions accounting for the effects of thermalization and dissipation, µ and y can be approximately calculated from the spectrum of density perturbations, P ζ [91,95]: where k min = 1 Mpc −1 .The strongest existing constraints on spectral distortions come from the COBE/FIRAS instrument, which restricts |µ| ≲ 9.0×10 −5 and |y| ≲ 1.5×10 −5 at the 95% C.L. [63,64].The standard cosmological model predicts spectral distortions of µ ∼ 2 × 10 −8 and y ∼ 10 −6 [98], consistent with current limits.The models of interest in this study have enhanced P ζ , and thus enhanced µ and y, which can be constrained by CMB measurements.
For concreteness, consider the case of a power spectrum that is sharply peaked at a single scale, k BH : Using the delta-function to perform the integral, µ and y become functions of k BH and σ 2 ζ alone.From the horizon crossing condition k = aH, and fact that H = 1.66 g ⋆ (T )T 2 /M Pl during radiation domination, we can relate the wavenumber to the temperature, which can then be related to the horizon mass using Eq. ( 1.3) to estimate the mass of the resulting black hole: For the sharply peaked spectrum of Eq. (3.4), adopting γ = 0.2, and for black holes in the mass range of interest, the µ-and y-parameters can be written as ) Note that these results should also hold, for example, in the case of a log-normal spectrum of sufficiently narrow width.In Eqs.(3.7) and (3.8), we can identify the impact of the various eras described earlier in this section.In particular, for M BH ≪ 1.5 × 10 5 M ⊙ (corresponding to k BH ≫ 5400/Mpc), the black holes are formed during the thermalization era, and both µ-and y-type spectral distortions are suppressed.For 1.5 × 10 5 M ⊙ ≪ M BH ≪ 4.5 × 10 9 M ⊙ (31.6/Mpc ≲ k BH ≲ 5400/Mpc), the black holes are forming during the µ-era, leading primarily to µ-type spectral distortions.Larger black holes form later, yielding primarily y-type spectral distortions.
We can use Eqs.(3.7) and (3.8) to quickly estimate whether a scenario featuring primordial SMBHs is consistent with spectral distortion constraints.Recall that, for the case of Gaussian statistics, a value of σ 2 ζ > ∼ 10 −2 is required in order to obtain a non-negligible abundance of PBHs.For such large values of σ 2 ζ , spectral distortions exclude all black holes masses M BH ≳ few×10 3 M ⊙ .This conclusion is consistent with Fig. 2, where we compare the bounds from COBE/FIRAS with the power spectra predicted for several values of M BH .Therefore, the existence of a non-negligible abundance of primordial SMBHs requires the presence of significant non-Gaussianities in the distribution of the primordial curvature perturbations.[99] coming from CMB temperature anisotropies (dark blue) [66], Lyman-α forest (light blue) [100], CMB spectral distortions (red) [63,64], and pulsar timing arrays (green) [101].The cusp in the COBE/FIRAS excluded region signifies the wavenumber where constraints from µ-and y-type distortions are equally restrictive.Overlaid are illustrative sharply peaked log-normal power spectra resulting in the formation of PBHs with M BH = 10 10 , 10 6 , 10 2 , 10 −2 , and 10 −6 M ⊙ and an initial abundance of β = 10 −20 , assuming Gaussian statistics for ζ.

Departures from Gaussianity
For Gaussian density perturbations, constraints from spectral distortions severely limit the abundance of primordial SMBHs that could have formed in the early universe.In this case, the variance must be small in order to limit spectral distortions, while a large variance is required to generate a non-negligible abundance of PBHs.This tension could be resolved, however, if the distribution of curvature perturbations features a heavier tail than that of a Gaussian.Such non-Gaussianities thus potentially allow primordial SMBHs to form 13 without necessarily violating spectral distortion constraints. 14ortunately, the distribution of primordial density perturbations is generically predicted to be non-Gaussian.Firstly, there is intrinsic non-Gaussianity that arises from the nonlinear mapping between the curvature perturbation ζ and the density contrast δ, as can be seen in Eq. (2.5).Thus, even if the probability distribution function for ζ were exactly Gaussian, the distribution in δ would not be.Secondly, and more crucially, large departures from Gaussianity are generically found in models which produce a local enhancement in the primordial power spectrum [82].
To quantify the degree of non-Gaussianity that would be required to generate primordial SMBHs without violating spectral distortion constraints, it is instructive to consider a class of probability distribution functions of the form [104] where n parameterizes the heaviness of the distribution's tail.The variance of the density contrast is set by the second moment of the distribution: Note that for n = 2, this reduces to the Gaussian form of Eq. (2.10) with σ 2 0 = σ 2 δ .For n = 1 the tail falls off exponentially, while for 0 < n < 1 it falls off even more slowly; we will refer to this class of distributions with n < 1 as "heavy-tailed." 15n Fig. 3, we show the shape of this class of probability distributions for various choices of n.As expected, we see that smaller n gives rise to heavier tails.This raises the question of how small n must be in order to efficiently form primordial SMBHs while keeping the peak of the power spectrum within bounds of spectral distortions, P ζ ≲ 10 −4 .In Fig. 4, we plot the maximum PBH mass fraction at formation β max , for a variance that saturates the spectral distortion constraints from COBE/FIRAS.Note that to generate a present day abundance of Ω BH ≳ 10 −20 with M BH < ∼ 10 11 M ⊙ , we need a heavy tail with n ≲ 0.6.
We now turn to the question of what inflationary models could yield such dramatic departures from Gaussianity.In the context of single-field inflation, Refs.[105][106][107][108] perturbatively studied the local non-Gaussianity that arises in models which deviate from the slow-roll attractor, as in ultra-slow-roll inflation.Going beyond perturbation theory, Refs.[46,51] used the δN formalism [109][110][111][112][113][114] to compute the non-perturbative distribution of curvature perturbations. 16For inflaton potentials with a small step or bump-like feature that induces a period of off-attractor behavior, these studies found that the tail of the distribution can become exponential, without inducing any significant non-Gaussianities in the perturbative regime [46].This result highlights the fact that perturbative measures of non-Gaussianity are not generally adequate to describe the large, rare fluctuations that lead to PBH formation.Finally, many of the mechanisms for enhancing local curvature perturbations rely on a temporary reduction in the inflaton's velocity.When the slow-roll classical drift vanishes, the field dynamics can receive large corrections from quantum diffusion, and the stochastic inflationary formalism [110,116] may be necessary for a proper description of the dynamics.Combining this with the δN formalism [117], a number of studies [47,49,50,[118][119][120] have found that prominent exponential tails arise generically from quantum diffusion.
While many single-field models have been found to yield exponential tails [121], there is currently no known model which generates a heavier-tailed distribution. 17However, as shown in Fig. 4, a heavy-tailed P δ ∼ exp (−|δ| n ) with n ≲ 0.6 is needed to yield a non-negligible population of primordial SMBHs while satisfying bounds from CMB spectral distortions.While it is unclear whether primordial SMBHs can appreciably form in any viable single-field models, it is plausible that the necessary heavy tails can arise in certain multi-field scenarios.In particular, curvaton models have long been known to generate density perturbations with sizeable non-Gaussianities in the case that the curvaton remains subdominant at the time of its decay [123].While the standard curvaton with quadratic potential is incapable of producing the dramatic non-Gaussianities needed for our scenario, there is good reason to believe that these can be augmented for a self-interacting curvaton model.We explore these possibilities in the following sections.

Standard Curvaton Scenario
Curvaton models introduce a second light, unstable spectator field that is present during inflation and that is responsible for generating the dominant contribution to the primordial curvature perturbations [124][125][126][127].The perturbations of the curvaton are initially isocurvature, but become adiabatic upon curvaton decay sometime after inflation ends [128].Due to the non-linearity inherent in this transfer, the full perturbation ζ can become quite non-Gaussian.In particular, when the curvaton is still very subdominant at decay, the inefficient conversion can yield a very heavy-tailed distribution for ζ.
Non-Gaussianity in the curvaton model was first investigated using the δN formalism in Ref. [123], with implications for PBH formation analyzed in Refs.[58,102,129].As we will see, the standard curvaton model with only a quadratic potential cannot produce sufficient non-Gaussianity to generate a non-negligible abundance of SMBHs without violating spectral distortions constraints.ζ that saturates the spectral distortion constraints from COBE/FIRAS, as estimated according to Eq. (3.7).We assume the distribution function given in Eq. (4.1) and consider Gaussian (n = 2, black), exponential (n = 1, red), and power law (n = 0.6 blue, n = 0.25 green) behavior in the tail.Note that the cusps which appear near M BH ∼ 3 × 10 9 M ⊙ correspond to the value of M BH at which µ-and y-type spectral distortions are equally restrictive.Contours of constant Ω BH , as computed using Eq.(2.2), are shown in dashed gray.

Curvaton Cosmology
We begin by reviewing the calculation of the curvature perturbation ζ and its statistics in the standard curvaton scenario [123], for which the total potential is where V ϕ is the unspecified potential of the inflaton ϕ, and V χ = m 2 χ χ 2 /2 is the quadratic potential of the curvaton χ.The curvaton mass is required to be light, satisfying m χ ≪ H throughout inflation, such that quantum perturbations are the dominant influence on its evolution.This also implies that the background field χ will remain effectively fixed at some initial value χ * during inflation, where a star denotes the value of quantities at horizon exit.For the curvaton energy density to remain subdominant throughout inflation, we demand χ * ≪ 2V ϕ /m 2 χ .Just like the inflaton, the curvaton receives perturbations δχ * ≃ H * /2π set by the Hubble rate at horizon exit H * .Since the curvaton is a weakly coupled field, we expect the perturbations δχ * to be described by a Gaussian random field.Thus we can write the curvaton at horizon exit as the sum of a background field and a linear perturbation, with no higher order terms: (5. 2) The goal of this section will be to relate these initial Gaussian field perturbations to the total curvature perturbation ζ via some mapping ζ = f (δχ * ).This will be the key to constructing the probability distribution function for ζ, since the statistics of a non-Gaussian variable are completely determined by the statistics of a Gaussian reference variable when the mapping between them is specified.When inflation ends and the inflaton decays, the universe enters into an era of radiation domination, with ρ R ∼ a −4 .At this point, the curvaton energy density is subdominant and its fluctuations are still isocurvature in nature.As the Hubble rate decreases, it eventually drops below m χ , causing the curvaton to start oscillating about the minimum of its potential.We denote the field value at which this occurs by χ0 .During this oscillating phase, the curvaton redshifts like matter with ρ χ ∝ a −3 and its energy density grows linearly relative to radiation.Finally, when H ∼ τ −1 , where τ is the χ lifetime, the curvaton decays to radiation and its isocurvature perturbations become adiabatic perturbations, assuming the decay products thermalize with the existing radiation.

The δN Formalism
To calculate the distribution of the curvature perturbations in this model, we employ the δN formalism [109][110][111][112][113][114]130], which we review in Appendix A. This technique identifies ζ on large scales (k ≪ aH) with the variation of inflationary e-folds across Hubble patches and non-perturbatively captures its non-Gaussianities.The δN formalism was first used to study non-Gaussianity in curvaton models in Refs.[123,131].On a general hypersurface of uniform curvaton density, the conserved curvaton curvature perturbation ζ χ is [113,131] where δN (t, ⃗ x) is the perturbed number of e-folds, ρ χ (t, ⃗ x) is the χ energy density, and ρχ (t) is its background value.In spatially flat slicing, this becomes and the curvaton energy density can be written as ρ χ (t, ⃗ x) = e 3ζχ(t,⃗ x) ρχ (t) . (5.5) In uniform total density slicing, Eq. ( 5.3) becomes where ζ is the total curvature perturbation. 18Since ζ and ζ χ are gauge invariant quantities, Eq. (5.5) can be equated with Eq. (5.6) to yield (5.7) A similar treatment for the radiation energy density in uniform total density slicing gives where we have assumed for simplicity that the main contribution to the curvature perturbation comes from the curvaton.In order to derive analytic results, we work in the instantaneous decay approximation such that the curvaton decays when H = τ −1 , where τ is the curvaton lifetime.On a uniform total density slice at t = τ , the energy densities satisfy where ρ(u) = ρ(u) R + ρ(u) χ is the total homogeneous energy density.Substituting Eqs.(5.7) and (5.8), this becomes a 4 th degree algebraic equation for ζ at τ : (5.10) Alternatively, it is customary to introduce the parameter r τ , defined as [124]: .11) in terms of which the equation for ζ becomes: (5.12) The general solution is where A = e 3ζχ and we have defined

Calculating the Probability Distribution
Finally, to obtain ζ χ in terms of δχ * , we return to Eq. (5.4), which gave the curvature perturbation in spatially flat slicing as a function of perturbed and background energy densities.
For the simple quadratic potential of this model, we have ρ χ (t, ⃗ x) = m 2 χ χ 2 /2.Expanding χ(t, ⃗ x) = χ(t) + δχ(t, ⃗ x), we can write this as where δ χ = δχ/ χ is the curvaton contrast in spatially flat slicing.Comparing with ρ χ = e 3ζχ ρχ , we see we should identify e 3ζχ = (1 + δ χ ) 2 .Finally in order to relate the contrast δ χ to its value at horizon exit, consider the equations of motion for χ and the perturbation δχ: On superhorizon scales k ≪ aH is negligible and so these reduce to the same equation.This implies δχ ∼ χ, such that δχ/ χ = δχ * / χ * and (5.18) Combining with Eq. ( 5.13), we obtain the desired mapping between ζ and the Gaussian initial field perturbations δχ * .We can now use probability conservation to write where P δχ is fully determined by the Gaussian variance σ 2 0 and the roots δ ± χ (ζ) satisfy which arise from solving Eq. (5.12) and substituting Eq. (5.18).In Fig. (5), we plot the probability distribution for the curvature perturbation as given by Eq. (5.19) for a few choices of r τ , which controls the heaviness of the tail.

Heavy Tails and Primordial Black Holes in Curvaton Models
Viably forming an appreciable number of primordial SMBHs requires both amplified power on small scales and a departure from Gaussianity.Thus we require the following two additional ingredients: • Enhanced Power: One mechanism for enhancing the power spectrum is to introduce a non-canonical kinetic term for the curvaton, which depends on the inflaton's field value [58].In Sec.6.1, we review this scenario and calculate the power spectrum resulting from such a kinetic term.
• Large Non-Gaussianity: The non-Gaussianity in a curvaton model can be amplified through self-interactions, which lead to non-linear growth of χ perturbations between horizon exit and the onset of curvaton oscillations [132][133][134][135][136][137][138][139][140].In Sec.6.2, we review the evolution of the non-quadratic curvaton and calculate the resulting probability distribution function.for various choices of r τ .The vertical dashed line corresponds to the threshold value for collapse ζ c , and the PBH mass fraction β is obtained by integrating beyond this threshold.We see that a smaller r τ corresponds to a heavier-tailed distribution, leading to a larger PBH abundance.Note that the reference value σ 2 ζ = 0.04 is chosen to illustrate the heaviness of the tail as r τ is varied, but the minimal scenario presented in Sec. 5 cannot realize such a large value without violating the observational limits on the power spectrum shown in Fig. 2.

Non-Minimal Curvaton Scenario
The minimal curvaton model described in Sec. 5 provides non-Gaussian statistics, but does not amplify P ζ .This deficiency can be remedied with a non-canonical kinetic term 19 for the curvaton [58]: If f (ϕ) is chosen such that this kinetic term is suppressed at field values ϕ = ϕ * , the power spectrum will be enhanced on scales corresponding to the horizon size at ϕ * .For concreteness, consider the evolution of the inflaton and curvaton governed by the following system of equations: where f ′ ≡ ∂ ϕ f and the source term of the first equation is negligible until the curvaton begins to oscillate.Similarly, the curvaton perturbation evolves according to which, to leading order, can be simplified inside the horizon, k ≫ aH, to yield whose solution can be written as which establishes the initial conditions in the adiabatic vacuum.As in Sec. 5, on superhorizon scales, k ≪ aH, χ and δχ evolve according to the same equation, so their solutions have the same functional form for all t > t * : where k = a(t * )H(t ⋆ ).Using Eq. (6.6), the power spectrum for δ χ becomes so, if f (ϕ) is chosen to have a dip at ϕ * , corresponding to k * = k BH , P δχ will exhibit a peak at k BH .However, for modes far away from k BH , f (ϕ * ) ≈ 1, recovering the nearly scale-invariant spectrum required for consistency with CMB observations on larger scales.Assuming f (ϕ) has such a localized feature, combining Eqs.(5.19), (6.7), and (2.8), the PBH abundance becomes where ) are the roots of Eq. (5.20) evaluated at the threshold.Using Eq. (3.7), we obtain the maximum value of β consistent with spectral distortion constraints. 20The degree of non-Gaussianity in this scenario is governed by r τ , as defined in Eq. (5.11).In the r τ → 1 limit, the curvaton dominates the energy density prior to its decay, so the relation between ζ and ζ χ from Eq. (5.12) is approximately linear.Thus, in this regime, the non-Gaussianity comes entirely from the non-linear relationship between ζ χ and δ χ ; see Eq. (5.18).In the opposite regime21 that the curvaton is still very subdominant when it decays (r τ ≪ 1), the relation between ζ and ζ χ is highly non-linear, and so the degree of non-Gaussianity is large.This is reflected in Figs. 5 and 6.Clearly this scenario is incapable of producing primordial SMBHs while satisfying spectral distortion constraints.

Self-Interacting Curvaton Scenario
Introducing curvaton self-interactions allows for non-linear evolution of δ χ between horizon exit and the onset of oscillations, which can lead to even more dramatic departures from Gaussianity.This scenario is also physically well-motivated since the curvaton needs to decay for the isocurvature perturbations to be converted to adiabatic perturbations.The effects of self-interactions in curvaton models were first investigated in Refs.[132][133][134][135][136][137][138][139][140].However, these studies restricted themselves to computing the non-linearity parameters f NL and g NL , which do not capture the non-perturbative statistics in the tail of the distribution.Unfortunately following the curvaton's evolution in the non-quadratic regime and computing the exact resulting curvature distribution require extensive numerics beyond the scope of this study.We offer here a schematic picture, but leave the detailed model building to future studies.We now allow the curvaton potential V (χ) to be an arbitrary well-defined function of χ.The cosmological evolution of the curvaton proceeds similarly to the case of the quadratic potential, but with a few crucial differences.At a time t = t int corresponding to V ′ (χ * ) ∼ H, χ begins rolling towards the minimum of its potential.We choose parameters such that this occurs while the curvaton's energy density is dominated by the interaction terms.At a later time t = t 0 > t int , these interaction terms become subdominant and the curvaton mass term drives field evolution, resulting in matter-like scaling ρ χ ∝ a −3 .Note that with self-interactions, the curvaton energy density generically falls off faster than in the quadratic case, resulting in a smaller r τ at the time of decay.
Recall that in the case of the quadratic potential, the curvaton density contrast δ χ remained constant after horizon exit since the background field χ and perturbation δχ obeyed the same equation of motion, shown in Eq. (5.17).This led to δχ/ χ = δχ * / χ * , which allowed δ χ to be used as a Gaussian reference variable.Upon introducing interactions this is no longer the case, as δ χ evolves non-trivially between t int and the onset of quadratic oscillations at t 0 .In this regime, the equation of motion for the perturbation is: where it is understood that the second derivative of the potential should be evaluated at the background field value.
A natural choice of Gaussian reference variable is the initial curvaton perturbation δχ * , which can be related to the ζ χ by first solving Eq. (6.9) along with the equation for the background field χ, computing the total and background energy densities at decay, and finally -18 -applying Eq. (5.4).The resulting δχ i * = g i (ζ χ ) can then be mapped onto the total curvature perturbation ζ via Eq.(5.12), since the relationship between ζ χ and ζ is unchanged in the presence of χ self-interactions.Although an exact solution requires the use of numerical techniques, an approximate relation can be derived in the limit of weak interactions.
We are interested in the relation between δχ * and ζ χ at the time of curvaton decay.For χ values sufficiently close to the minimum of its potential, the potential is approximately quadratic, and the energy density is where χ 0 is the amplitude at the onset of oscillations.Since non-linear evolution takes place between horizon exit and oscillation, this initial amplitude is a function of initial field value χ * ).In terms of background field values χ0 ≡ χ 0 ( χ * ), this can be expanded as: and the energy density can then be written as: where ρχ = 1 2 m 2 χ χ2 0 .Comparing with Eq. (5.4), the bracketed quantity is identified with e 3ζχ .Then also expanding we can write to leading order: which can then be substituted into Eq.(5.13) to obtain ζ as a function of the Gaussian reference variable δχ * .Note that though we have written this expression to fixed order, higher order terms can become significant in the presence of large non-Gaussianities.It is also possible to invert this mapping to obtain the roots g i (ζ) = δχ i * (ζ).The probability distribution is then: where the sum runs over all real roots.In Fig. 7, we plot the maximal PBH abundance consistent with spectral distortion constraints for a sample set of parameters.Comparison with β max in the standard curvaton scenario (see Fig. 6) reveals the dramatic amplification of non-Gaussianity that self-interactions can afford.We leave the determination of a potential that could actually yield these parameter values to future investigation.⇣ saturates spectral distortion bounds for differnt values of r ⌧ .We assume ⇣ the take the form of Eq. ( 5.41) with the potential chosen such that ⇣ ,2 = 0.5⇣ ,1 and ⇣ ,3 = 0.1⇣ ,1 .Contours of constant ⌦ BH today are shown in dashed gray.
Comparison with max in the standard curvaton scenario (see Fig. 5) reveals the dramatic amplification of non-Gaussianity that even weak interactions can afford.Note that while the non-linear growth of modes between horizon exit and the onset of oscillations significantly boosts the potential amount of non-Gaussianity, the curvaton still needs to be very subdominant at the time of its decay to produce a sufficiently heavy-tailed distributioni.e.we must have r ⌧ ⌧ 1.It is actually even more natural to fulfill this condition here since energy density in the curvaton field dilutes more quickly upon introducing interactions.

Conclusions
Much remains to be understood about the origin and evolution of our universe's most massive supermassive black holes.The inferred number and population distribution of SMBH in the high-redshift universe is perhaps surprising, and challenges the standard assumption that they formed from initially low mass seeds which gradually grew through accretion and possibly mergers.In this work we have instead motivated the possibility that some fraction of this population may be primordial in origin, having formed from the collapse of overdensities seeded by inflation.The main obstruction to this scenario comes from measurements of spectral distortions of the CMB, which constrain the amplitude of the primordial power spectrum those from CMB spectral distortions in our mass range of interest.In particular, Ref. [110] finds uses dark matter substructure -specifically the observed number of dwarf spheroidal galaxies and stellar streamsto constrain P ⇣ .While the limits derived in this paper likely depend on the specific form of the bump in the potential (and consequent form of the variance), as well as modeling of the evolution of the host halos and subhalos, this could potentially rule out some fraction of the parameter space for SMBH.

Discussion and Conclusions
Much remains to be understood about the origin and evolution of our universe's most massive black holes.The inferred population of supermassive black holes with M BH ∼ 10 6 − 10 11 M ⊙ in the high-redshift universe is perhaps surprising, and challenges the standard assumption that these objects formed from low mass seeds which grew through the processes of accretion and mergers.In this study, we have taken this as motivation to consider the possibility that some of our universe's supermassive black holes may be primordial in origin, having formed from the direct collapse of overdensities seeded by inflation.
Forming primordial SMBH from direct collapse requires an enhanced power spectrum on small scales (k ∼ 10 − 10 4 Mpc −1 ), which results in dangerous CMB spectral distortions.Since the CMB exhibits a nearly perfect blackbody spectrum, such distortions exclude the possibility that a population of primordial SMBH could have formed from Gaussian density perturbations.However, if the distribution of primordial curvature perturbations were highly non-Gaussian, it is possible that primordial black holes may have formed from smaller peaks in the power spectrum.To evade limits from spectral distortions, the tail of the probability distribution must be very heavy, falling off as P ζ ∼ exp (−|ζ| n ) with n ≲ 0.6; we are not aware of any single-field inflationary model that can realize this behavior.
In this paper, we have explored the possibility of generating curvature perturbations with heavy-tailed probability distributions in the event that a self-interacting curvaton field is also present as a spectator during inflation.In the standard curvaton scenario, non-Gaussianity -20 -arises from the inefficient conversion of isocurvature perturbations into adiabatic perturbations when the curvaton decays.However, the degree of non-Gaussianity in this minimal realization is insufficient to yield an appreciable primordial SMBH population.Introducing curvaton self-interactions results in non-linear evolution of the curvaton contrast δ χ between horizon exit and the onset of quadratic oscillations, potentially resulting in a heavier-tailed distribution, as seen in Fig. 7.We leave the detailed model building and determination of a potential that can realize this benchmark point to future investigation.Provided such a potential can be found, the non-minimal self-interacting curvaton could viably yield SMBHs from the direct collapse of primordial perturbations, without violating spectral distortion constraints.
We note that beyond the bounds from CMB spectral distortions, there are two further more speculative constraints that potentially need to be addressed in a full model.First, as pointed out by Ref. [143], the small scale curvature perturbation impacts structure evolution as traced by dark matter halos.Presuming a bump in the small scale power spectrum, they investigated host halo and subhalo evolution.By comparing the results of their simulations with the observed number of dwarf spheroidal galaxies and stellar streams, they were able to place bounds on P ζ to greater sensitivity than COBE/FIRAS for a region of wavenumbers.The result would exclude some, but not all, of the SMBH mass range we consider.These bounds are less robust than those coming from CMB spectral distortion and are subject to modeling uncertainties, but nevertheless should be noted.Second, the large non-Gaussianities described in this model would likely lead to heavy clustering of quasars.Ref. [144] calculated the observed angular correlation function using the 92 quasars at z ∼ 6 reported by the Subaru High-z Exploration of Low-Luminosity Quasars project and compared this with that predicted in a PBH model which used a light spectator field to source non-Gaussian primordial perturbations.They found the amplitude of the angular correlation function predicted by the model to be much larger than that observed, a result which would seem to preclude PBHs as the sole progenitors of SMBHs in this redshift range.More concretely, they were able to restrict the fraction of PBH originating SMBHs to ≲ 10 −4 .
These potential constraints notwithstanding, this work is especially timely in light of two distinct recent observations.The James Webb Space Telescope (JWST) has reported a number of surprisingly luminous high-redshift galaxy candidates [145][146][147][148] whose existence poses a challenge to the standard ΛCDM paradigm.These massive early galaxies could conceivably be explained if primordial black holes accelerated galaxy formation in the early universe [149].Meanwhile, the NANOGrav collaboration and several other pulsar timing array experiments have just announced evidence of a signal consistent with the stochastic gravitational wave background in the nHz frequency range [150][151][152][153].The leading astrophysical interpretation of this signal is that it consists of gravitational waves from supermassive black hole binary mergers.However, some aspects of this data, such as the frequency scaling of the spectral density parameter, are not particularly well-fit by this interpretation [154,155].Given that the distribution of SMBH binaries would be different if these objects were of a primordial origin, one avenue for future investigation would be to compute the gravitational wave signal predicted in this scenario.While for homogenously distributed primordial SMBHs this possibility is in conflict with constraints on total abundance from large scale structure and the UV galaxy luminosity function [156], a highly clustered population could still viably source the gravitational wave signal [157].Given that a clustered population can arise in the presence of non-Gaussianities, it would be interesting to see what our model predicts in this context; see Refs.[144,158] for closely related work.There are also other signals, such as scalar-induced -21 -secondary gravitational waves, which could offer complimentary evidence of this scenario and deserve further study; see Refs.[159,160] for a potential connection between scalar-induced gravitational waves and the PTA signal.Regardless, these recent observations provide us with strong motivation to better understand the cosmic origin of our universe's supermassive black holes.

Figure 3 .
Figure 3. Probability distribution functions Eq. 4.1 for various choices of n at fixed β (left) and fixed variance σ 2 (right).We see that smaller n corresponds to heavier-tailed distributions.

Figure 4 .
Figure 4.The maximum primordial black hole mass fraction at formation β max as a function of mass M BH for a value of the variance σ 2ζ that saturates the spectral distortion constraints from COBE/FIRAS, as estimated according to Eq. (3.7).We assume the distribution function given in Eq. (4.1) and consider Gaussian (n = 2, black), exponential (n = 1, red), and power law (n = 0.6 blue, n = 0.25 green) behavior in the tail.Note that the cusps which appear near M BH ∼ 3 × 10 9 M ⊙ correspond to the value of M BH at which µ-and y-type spectral distortions are equally restrictive.Contours of constant Ω BH , as computed using Eq.(2.2), are shown in dashed gray.

Figure 5 .
Figure 5.The probability distribution function for the curvature perturbation ζ as given by Eq.(5.19)   for various choices of r τ .The vertical dashed line corresponds to the threshold value for collapse ζ c , and the PBH mass fraction β is obtained by integrating beyond this threshold.We see that a smaller r τ corresponds to a heavier-tailed distribution, leading to a larger PBH abundance.Note that the reference value σ 2 ζ = 0.04 is chosen to illustrate the heaviness of the tail as r τ is varied, but the minimal scenario presented in Sec. 5 cannot realize such a large value without violating the observational limits on the power spectrum shown in Fig.2.

Figure 6 .
Figure 6.The maximal PBH mass fraction at formation β max in the standard curvaton scenario for σ 2ζ subject to spectral distortion constraints for various choices of r τ from Eq.(5.11)

Figure 6 .
Figure 6.Maximal PBH mass fraction at formation max as a function of PBH mass in the selfinteracting curvaton model, where 2⇣ saturates spectral distortion bounds for differnt values of r ⌧ .We assume ⇣ the take the form of Eq. (5.41) with the potential chosen such that ⇣ ,2 = 0.5⇣ ,1 and ⇣ ,3 = 0.1⇣ ,1 .Contours of constant ⌦ BH today are shown in dashed gray.