# Supermassive Primordial Black Holes From Inflation

Dan Hooper  Aurora Ireland  Gordan Krnjaic  Albert Stebbins

###### Abstract

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 CMB spectral
distortion constraints. Such limits, however, can be avoided if the
distribution of primordial curvature perturbations is highly non-Gaussian. In
this study, we present a model of multi-field inflation — the curvaton model
supplemented with self-interactions — which can viably yield such dramatic
non-Gaussinities. Furthermore, we calculate the maximal abundance of black
holes that can be generated in this scenario and find this to be consistent
with the observed population of high-redshift SMBHs. This result is
particularly timely in light of recent evidence from the NANOGrav experiment
for a stochastic gravitational wave background consistent with SMBH mergers.

FERMILAB-PUB-23-390-T

##  1 Introduction

Supermassive black holes (SMBHs) are ubiquitous in our universe, being present
in the centers of nearly all massive galaxies. Quasars,111Quasars 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. powered by black holes with with masses $M\sim
10^{8-10}\,M_{\odot}$, 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\sim 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 $\Omega_{\rm SMBH}(z=0)\sim
10^{-6}$ [11]. This contrasts with the peak quasar mass density of
$\Omega_{\rm quasar}(z=2)\sim 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 mass accretion is Eddington limited to

$\dot{M}_{\rm Edd}\lesssim\frac{M_{\rm BH}}{\tau_{S}},$ (1.1)

where

$\tau_{\mathrm{S}}=\frac{\varepsilon\,\sigma_{\mathrm{T}}}{4\pi\,G\,m_{\mathrm{p}}}\approx
45\,\,{\rm Myr}\,\bigg{(}\frac{\varepsilon}{0.1}\bigg{)}\,$ (1.2)

is the Salpeter time [16], $\sigma_{\mathrm{T}}=8\pi\,\alpha^{2}/(3m_{e}^{2})$
is the Thomson cross section, $m_{p}$ is the proton mass, and $\varepsilon$ is
the radiative efficiency. At the maximum growth rate, a $10^{2}\,M_{\odot}$
black hole seed would require $\sim 0.8\,\mathrm{Gyr}$ to grow to $M\sim
10^{10}\,M_{\odot}$. Since $z\sim 6-7$ corresponds only to $\sim 0.7-0.9\,{\rm
Gyr}$ after the Big Bang, such a scenario would require that the largest high-
redshift black holes must have grown at high accretion rates almost
continuously during the first Gyr of our universe’s history. This contrasts
with the intermittent accretion that SMBHs undergo at lower redshifts. From
this perspective, it is surprising that so many highly-massive quasars have
been observed at such high redshifts [17, 18, 19, 20, 21].222Ref. [22]
proposes a mechanism relying on large initial clustering to generate massive
PBH in the absence of efficient accretion.

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.

There is a significant population of $M_{\rm BH}\gtrsim 10^{9}\,M_{\odot}$
quasars at $z\sim 6$ to $7$ [23, 4, 24, 25]. The $M_{\rm BH}\gtrsim
10^{10}\,M_{\odot}$ population has remained approximately constant since at
least $z\sim 5$ [26] (see the red curve of Fig. 1) which, curiously, requires
any growth in the number of SMBHs to be balanced by a decrease in the fraction
of SMBHs that are actively quasars. Also curious is the observation that the
maximum quasar mass has not increased much since $z\sim 5$. This leaves us
with two intriguing questions:

  1. 1.

If these quasars grew from small black hole seeds, how did they come to be so
massive on such a short timescale?

  2. 2.

Why did their growth rate dramatically slow down during the subsequent
$13\,{\rm Gyr}$?

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$-$\sigma$ relation333The
observed correlation between the velocity dispersion of a stellar bulge
$\sigma$ and the mass of the SMBH at its center. to evolve with redshift and
for the quasar luminosity function to steepen at the highest values [27, 28].
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 [29]. 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. Unlike
smaller primordial black holes (PBH), these objects would have formed at later
times, as governed by the size of the cosmological horizon. During the
radiation dominated era, the horizon contains the following amount of energy:

$M_{H}=\frac{4\pi}{3}\rho R_{H}^{3}\approx 3\times
10^{9}\,M_{\odot}\,\bigg{(}\frac{10\,{\rm
keV}}{T}\bigg{)}^{2}\bigg{(}\frac{3.36}{g_{\star}(T)}\bigg{)}^{1/2}\,,$ (1.3)

where $R_{H}=H^{-1}$ is the size of the horizon,
$\rho=\pi^{2}g_{\star}(T)T^{4}/30$ is the radiation density, and $g_{\star}$
is the number of relativistic species at temperature $T$. When a sufficiently
large density fluctuation collapses to form a PBH,444There exist other PBH
formation mechanisms beyond the collapse of overdensities seeded by inflation.
For example, PBH could form during first order cosmological phase transitions
[30] or due to the collapse of supercritical vacuum bubbles nucleated during
inflation [31]. The prospect of primordial SMBHs from these mechanisms is
discussed in Ref. [32] and Ref. [33], respectively. the mass of the resulting
black hole is typically an order one fraction of the horizon mass, $M_{\rm
BH}\simeq\gamma\,M_{H}$, where $\gamma\sim 0.2$ quantifies the efficiency of
collapse [34]. 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\sim{\rm
MeV}$ corresponds to $M_{H}\sim 10^{5}\,M_{\odot}$) but well before the onset
of matter domination or recombination ($T\sim{\rm eV}$ corresponds to
$M_{H}\sim 3\times 10^{17}\,M_{\odot}$).

Many inflationary models enable PBHs to form efficiently [35], including those
in which the inflaton undergoes a period of ultra-slow-roll [36, 37, 38, 39],
or whose potential features localized bumps, dips, or steps on small scales
[40, 41, 42]. More generically, a local enhancement of the power spectrum
$\mathcal{P}_{\zeta}$ requires a deviation from slow-roll evolution (see for
example, Refs. [43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57]).
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 $\mathcal{P}_{\zeta}(k_{\rm BH})\sim
10^{-2}$, seven orders of magnitude larger than observed on large scales,
$\mathcal{P}_{\zeta}(k_{\rm CMB})\simeq 2.1\times 10^{-9}$ [58]. Such an
amplification inevitably leads to large CMB spectral distortions. In
particular, for scales $k_{\rm BH}\mathrel{\raise
1.29167pt\hbox{$<$\kern-7.5pt\lower 4.30554pt\hbox{$\sim$}}}10^{5}\,{\rm
Mpc}^{-1}$, corresponding to black holes in the mass range $M_{\rm BH}\gtrsim
10^{3}M_{\odot}$, the predicted spectral distortions are in strong conflict
with COBE/FIRAS measurements [59, 60].

In principle, PBHs can form from smaller peaks in the power spectrum if the
tail of the $\zeta$ distribution is sufficiently non-Gaussian. Of course,
observations on large scales and measurements of the non-linearity parameters
$f_{\rm NL}$ and $g_{\rm NL}$ seem to indicate that $\zeta$ is very nearly
Gaussian [61, 62]. However, such perturbative measures of non-Gaussianity do
not capture the tail of the distribution, where fluctuations can be large and
rare [42]; it is precisely this tail that governs PBH formation.

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 primoridal SMBHs, so in Sec. 6 we introduce self-interactions to
augment the non-Gaussianity. We calculate the maximum PBH mass fraction that
is compatible with spectral distortion bounds in this model and find that it
is capable of accounting for the observed abundance of SMBHs. We conclude in
Sec. 7 with a summary of our results and offer some comments on directions for
future investigations.

##  2 Primordial Black Holes from Gaussian Perturbations

The initial abundance of PBHs is quantified in terms of their mass fraction at
the time of formation:

$\beta\equiv\frac{\rho_{\rm BH}}{\rho_{\rm tot}}\,.$ (2.1)

For PBHs that formed during radiation domination, this quantity is related to
the fractional energy density in black holes today, $\Omega_{\rm
BH}\simeq\beta\,(1+z_{f})\,\Omega_{R}$, 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:

$\Omega_{\rm BH}\simeq
10^{-6}\,\left(\frac{\beta}{10^{-10}}\right)\bigg{(}\frac{\gamma}{0.2}\bigg{)}^{1/2}\left(\frac{10^{8}\,M_{\odot}}{M_{\rm
BH}}\right)^{1/2}\left(\frac{g_{\star,s}(T_{f})}{3.91}\right)^{1/3}\left(\frac{3.36}{g_{\star}(T_{f})}\right)^{1/4},$
(2.2)

where $T_{f}$ is the temperature at the time of black hole formation.

In the standard555One can also calculate the black hole abundance using peak
theory [63, 64, 65, 66]. Unlike Press-Schechter, where the overdensity must
simply exceed the threshold, peak theory further demands that it be a local
maximum. This formalism has been demonstrated to be more appropriate when
perturbations exist on multiple scales. As we consider a sharply peaked
spectrum, the simpler Press-Schechter prescription suffices for our purposes.
treatment based on the Press-Schechter formalism [67], $\beta$ is computed by
integrating the probability distribution $P_{\delta}$ for the coarse-grained
density contrast $\delta=\delta\rho/\bar{\rho}$ over all values greater than
the critical threshold for collapse, $\delta_{c}$:

$\beta\simeq 2\int_{\delta_{c}}^{\infty}d\delta\,P_{\delta}[\delta]\,.$ (2.3)

The factor of $2$ is customarily introduced666It is unclear whether this
factor should still be included when considering asymmetric probability
distribution functions, as in the case of non-Gaussianities. We retain it
nevertheless since this is an ultimately inconsequential $\mathcal{O}(1)$
effect. to compensate for the undercounting that otherwise arises [67]. 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 PBH — and a simple implementation. One issue concerns the
uncertainty777It has been demonstrated that a more precise quantity which does
away with this uncertainty is the compaction function $\mathcal{C}$, defined
as twice the local mass excess over the comoving areal radius [68]. The value
of the compaction at its peak is interpreted as the threshold. See also Ref.
[69]. of the collapse threshold, which should in principle depend on the
overdensity profile.888Since the overdensity profile is modified in the
presence of non-Gaussianities, these can also change the threshold for
collapse [70, 71, 72, 73]. Given the uncertainties involved in these
calculations, we take a conservative position and determine $\delta_{c}$ using
Eq. (2.4).. 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 $\delta_{c}=c_{s}^{2}=1/3$ [34], where $c_{s}$ is the sound speed
of density perturbations. A more careful treatment that is applicable for an
arbitrary equation of state, $w$, finds [74]:

$\delta_{c}=\frac{3(1+w)}{5+3w}\sin^{2}\left(\frac{\pi\sqrt{w}}{1+3w}\right)\,,$
(2.4)

which has been shown to faithfully replicate the results of numerical
simulations (provided a sufficiently spherical profile) and yields
$\delta_{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 $\zeta$, 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, $\zeta$ is related to the
density contrast field by [75]

$\delta=-\frac{2(1+w)}{5+3w}\left(\frac{1}{aH}\right)^{2}e^{-2\zeta}\left(\nabla^{2}\zeta+\frac{1}{2}(\partial_{i}\zeta)^{2}\right)\,.$
(2.5)

Working to linear999A certain degree of non-Gaussianity inevitably arises in
the density contrast field due to this nonlinear relation. This implies that
even if the statistics of $\zeta$ were perfectly Gaussian, those of $\delta$
would not be. See, for example, Ref. [75] for discussion. order in radiation
domination, this simplifies to the Fourier space relation:

$\delta_{k}\simeq\frac{4}{9}\left(\frac{k}{aH}\right)^{2}\zeta_{k}\,,$ (2.6)

which implies that their power spectra, defined for arbitrary Fourier variable
$f_{k}$ via the 2-point function as
$\expectationvalue{f_{k}f_{k^{\prime}}}=\frac{2\pi^{2}}{k^{3}}\mathcal{P}_{f}(k)\delta^{(3)}(\vec{k}+\vec{k}^{\prime})$,
are related as

$\mathcal{P}_{\delta}(k)\simeq\frac{16}{81}\left(\frac{k}{aH}\right)^{4}\mathcal{P}_{\zeta}(k)\,.$
(2.7)

Note that at the time of horizon crossing $k\simeq aH$, when a perturbation
can collapse to form a black hole,101010There is additional uncertainty in the
calculation of PBH abundance coming from non-linear effects around the time of
horizon crossing. See Ref. [76]. the density contrast and curvature
perturbation are linearly related, $\delta\simeq\frac{4}{9}\zeta$. We can then
assume that peaks in $\delta$ also correspond to peaks in $\zeta$, and work
directly with the curvature perturbation. In terms of $\zeta$, the initial
abundance in the Press-Schechter formalism reads:

$\beta\simeq 2\int_{\zeta_{c}}^{\infty}d\zeta\,P_{\zeta}[\zeta]\,,$ (2.8)

where the collapse threshold $\zeta_{c}\simeq\frac{9}{4}\delta_{c}$ follows
from $\sigma_{\delta}^{2}\simeq\frac{16}{81}\sigma_{\zeta}^{2}$, since
$\mathcal{P}\sim\sigma^{2}$ for a sharply peaked spectrum. More precisely, the
variance of $\zeta$ smoothed on the scale $R\simeq(aH)^{-1}\simeq k^{-1}$ can
be computed from the power spectrum as [77]

$\sigma_{\zeta}^{2}(R)\equiv\expectationvalue{\zeta}^{2}_{R}=\int_{0}^{\infty}\frac{dk}{k}\tilde{W}^{2}(k,R)\mathcal{P}_{\zeta}(k)\,,$
(2.9)

where $\tilde{W}(k,R)$ is the Fourier transform of the (real space) window
function used to coarse-grain $\delta$. 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:
$\tilde{W}(k,R)=\exp\left(-k^{2}R^{2}/2\right)$, but see Ref. [78] for a
discussion of the resultant uncertainties.

Consider the case of a Gaussian distribution $P_{\zeta}=P_{\zeta}^{G}$, where

$P_{\zeta}^{G}=\frac{1}{\sqrt{2\pi}\sigma_{\zeta}}e^{-\zeta^{2}/2\sigma_{\zeta}^{2}}\,.$
(2.10)

In this case the mass fraction at formation evaluates to

$\beta=\text{erfc}\left(\frac{\zeta_{c}}{\sqrt{2}\sigma_{\zeta}}\right)\simeq\sqrt{\frac{2}{\pi}}\frac{\sigma_{\zeta}}{\zeta_{c}}\exp\left(-\frac{\zeta_{c}^{2}}{2\sigma_{\zeta}^{2}}\right)\,,$
(2.11)

where erfc is the complimentary error function and the second approximation
holds for $\zeta_{c}\gg\sigma_{\zeta}$, 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 $\beta=10^{-20}$, we see we need $\sigma_{\zeta}^{2}\simeq 0.01$,
corresponding to a peak in the power spectrum of $\mathcal{P}_{\zeta}\sim
0.01$, which is 7 orders of magnitude greater than the value measured on CMB
scales. Since the amplification of $\mathcal{P}_{\zeta}$ needed for PBH
formation depends on $\beta$ only logarithmically, this degree of enhancement
is a generic requirement for any non-vanishing initial abundance.

##  3 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 [79, 80, 81, 82, 83]. These deviations are known as
spectral distortions and in this context it is useful to distinguish between
three characteristic redshift intervals:

  * •

Thermalization era ($\bm{z>2\times 10^{6}})$: At high redshifts, Compton
scattering $\gamma e\to\gamma e$, double Compton scattering $\gamma
e\to\gamma\gamma e$, and Bremsstrahlung $ep\to ep\gamma$ maintain a blackbody
spectrum for the photons, and spectral distortions are exponentially
suppressed.

  * •

$\bm{\mu}$-era $(\bm{2\times 10^{5}<z<2\times 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, $\mu$. A $\mu$-distortion refers to a Bose-
Einstein distribution with $\mu\neq 0$.

  * •

$\bm{y}$-era ($\bm{z<2\times 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.111111Technically, the division between the two types of
spectral distortions is not entirely unambiguous, and inhomogeneities
dissipating around $z\sim 5\times 10^{4}$ can give rise to distortions of an
intermediate type, whose shape is not simply the sum of $\mu$\- and $y$-type
distortions [84].

A spectrum with a (positive) $y$-distortion can be expressed as an average of
blackbodies with slightly different temperatures [85, 86, 87, 88]. An average
of blackbodies with a mean temperature $\bar{T}$ and variance
$\bar{T}^{2}\,\Delta$ will (for $\Delta\ll 1$) be a $y$-distorted blackbody
characterized by the temperature $T=\bar{T}\,(1+\Delta^{2})$ and the parameter
$y=\Delta^{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 [89].

$\mu$\- and $y$-type spectral distortions are traditionally quantified in
terms of the parameters $\mu$ 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) [79, 83]:

$\mu\simeq
1.4\,\frac{\Delta\rho_{\gamma}}{\rho_{\gamma}}\qquad\textrm{and}\qquad y\simeq
0.25\,\frac{\Delta\rho_{\gamma}}{\rho_{\gamma}}\ .$ (3.1)

By introducing $k$-space window functions accounting for the effects of
thermalization and dissipation, $\mu$ and $y$ can be approximately calculated
from the spectrum of density perturbations, $P_{\zeta}$ [87, 83]:

$\displaystyle\mu\simeq 2.2\int_{k_{\rm
min}}^{\infty}\frac{dk}{k}\mathcal{P}_{\zeta}(k)\left[\exp\left(-\frac{k}{5400\,\text{Mpc}^{-1}}\right)-\exp\left(-\left[\frac{k}{31.6\,\text{Mpc}^{-1}}\right]^{2}\right)\right]\,,$
(3.2) $\displaystyle y\simeq 0.4\int_{k_{\rm
min}}^{\infty}\frac{dk}{k}\mathcal{P}_{\zeta}(k)\exp\left(-\left[\frac{k}{31.6\,\text{Mpc}^{-1}}\right]^{2}\right)\,,$
(3.3)

where $k_{\rm min}=1\,\text{Mpc}^{-1}$.

The strongest existing constraints on spectral distortions come from the
COBE/FIRAS instrument, which restricts $|\mu|\lesssim 9.0\times 10^{-5}$ and
$|y|\lesssim 1.5\times 10^{-5}$ at the 95% C.L. [59, 60]. The standard
cosmological model predicts spectral distortions of $\mu\sim 2\times 10^{-8}$
and $y\sim 10^{-6}$ [90], consistent with current limits. The models of
interest in this study have enhanced $P_{\zeta}$, and thus enhanced $\mu$ 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_{\rm BH}$:

$\mathcal{P}_{\zeta}(k)=\sigma_{\zeta}^{2}\,k\,\delta(k-k_{\rm BH}).$ (3.4)

Using the delta-function to perform the integral, $\mu$ and $y$ become
functions of $k_{\rm BH}$ and $\sigma_{\zeta}^{2}$ alone. From the horizon
crossing condition $k=aH$, and fact that
$H=1.66\sqrt{g_{\star}(T)}T^{2}/M_{\rm Pl}$ during radiation domination, we
can relate the wavenumber to the temperature,

$k_{\rm BH}=92\,{\rm Mpc}^{-1}\,\left(\frac{T}{10\,{\rm
keV}}\right)\,\left(\frac{g_{\star}(T)}{3.36}\right)^{1/2}\,\left(\frac{3.91}{g_{\star,S}(T)}\right)^{1/3}\,,$
(3.5)

which can then be related to the horizon mass using Eq. (1.3) to estimate the
mass of the resulting black hole:

$M_{\rm BH}\simeq 5\times
10^{8}\,M_{\odot}\,\left(\frac{92\,\text{Mpc}^{-1}}{k_{\rm
BH}}\right)^{2}\left(\frac{\gamma}{0.2}\right)\left(\frac{g_{\star}(T)}{3.36}\right)^{1/2}\left(\frac{3.91}{g_{\star,S}(T)}\right)^{2/3}\,.$
(3.6)

For the sharply peaked spectrum of Eq. (3.4), adopting $\gamma=0.2$, and for
black holes in the mass range of interest, the $\mu$\- and $y$-parameters can
be written as

$\displaystyle\mu\simeq
2.2\sigma_{\zeta}^{2}\,\Bigg{\\{}\exp\left[-\left(\frac{1.5\times
10^{5}\,M_{\odot}}{M_{\rm
BH}}\right)^{1/2}\right]-\exp\left[-\left(\frac{4.5\times
10^{9}\,M_{\odot}}{M_{\rm BH}}\right)\right]\Bigg{\\}}\,,$ (3.7)
$\displaystyle y\simeq 0.4\sigma_{\zeta}^{2}\,\exp\left[-\left(\frac{4.5\times
10^{9}\,M_{\odot}}{M_{\rm BH}}\right)\right]\,.$ (3.8)

Note that these results should also hold, for example, in the case of a log-
normal spectrum of sufficiently narrow width.

Figure 2: Constraints on the primordial power spectrum $\mathcal{P}_{\zeta}$
[91] coming from CMB temperature anisotropies (dark blue) [62], Lyman-$\alpha$
forest (light blue) [92], CMB spectral distortions (red) [59, 60], and pulsar
timing arrays (green) [93]. The cusp in the COBE/FIRAS excluded region
signifies the wavenumber where constraints from $\mu$\- and $y$-type
distortions are equally restrictive. Overlaid are illustrative sharply peaked
log-normal power spectra resulting in the formation of PBHs with $M_{\rm
BH}=10^{10},10^{6},10^{2},10^{-2},\,\text{and}\,10^{-6}M_{\odot}$ and an
initial abundance of $\beta=10^{-20}$, assuming Gaussian statistics for
$\zeta$.

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_{\rm BH}\ll 1.5\times
10^{5}\,M_{\odot}$ (corresponding to $k_{\rm BH}\gg 5400/\textrm{Mpc}$), the
black holes are formed during the thermalization era, and both $\mu$\- and
$y$-type spectral distortions are suppressed. For $1.5\times
10^{5}\,M_{\odot}\ll M_{\rm BH}\ll 4.5\times 10^{9}\,M_{\odot}$
($31.6/\textrm{Mpc}\lesssim k_{\rm BH}\lesssim 5400/\textrm{Mpc}$), the black
holes are forming during the $\mu$-era, leading primarily to $\mu$-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.
121212 The validity of the formulas used in this paper has been explored in
Ref. [94]. We believe they are valid across in our regime of interest, but may
differ at most by order one factors for large amplitude inhomogeneities.
Recall that, for the case of Gaussian statistics, a value of
$\sigma_{\zeta}^{2}\mathrel{\raise 1.29167pt\hbox{$>$\kern-7.5pt\lower
4.30554pt\hbox{$\sim$}}}10^{-2}$ is required in order to obtain a non-
negligible abundance of PBHs. For such large values of $\sigma_{\zeta}^{2}$,
spectral distortions exclude all black holes masses $M_{\rm
BH}\gtrsim\text{few}\times 10^{2}M_{\odot}$. 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_{\rm 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.

##  4 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 without necessarily violating spectral distortion constraints.131313Ref.
[95] proposes an alternative way to obtain the extreme statistics needed to
form primordial SMBHs compatible with spectral distortion bounds. Their multi-
field model results in effectively two different inflationary histories for
the casually disconnected Hubble patches, with a subdominant fraction
experiencing more expansion. By the $\delta N$ formalism this equates to large
curvature perturbations, such that these patches collapse to form PBH. Since
this occurs only in a tiny fraction of patches, there is no observable
generation of spectral distortions.

Fortunately, the distribution of primordial density perturbations is
generically predicted to be non-Gaussian. Firstly, there is intrinsic non-
Gaussianity that arises from the non-linear mapping between the curvature
perturbation $\zeta$ and the density contrast $\delta$, as can be seen in Eq.
(2.5). Thus, even if the probability distribution function for $\zeta$ were
exactly Gaussian, the distribution in $\delta$ 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
[75].

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 [95]

$P^{(n)}_{\delta}=\frac{1}{2\sqrt{2}\,\sigma_{0}\,\Gamma\left(1+\frac{1}{n}\right)}\exp\left[-\left(\frac{|\delta|}{\sqrt{2}\sigma_{0}}\right)^{\\!n}\,\right]\,,$
(4.1)

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:

$\sigma_{\delta}^{2}(\sigma_{0})\equiv\int_{-\infty}^{\infty}d\delta\,\delta^{2}\,P^{(n)}_{\delta}=\frac{\,\,\sigma_{0}^{2}\,\Gamma\left(1+\frac{3}{n}\right)}{3\Gamma\left(1+\frac{1}{n}\right)}~{}.$
(4.2)

Note that for $n=2$, this reduces to the Gaussian form of Eq. (2.10) with
$\sigma_{0}^{2}=\sigma_{\delta}^{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.”141414Formally,
the probability distribution $P_{x}$ of a random variable $x$ is said to be
“heavy” if its tail is not exponentially bounded,
$\lim_{x\rightarrow\infty}e^{\lambda x}\bar{F}(x)=\infty\,\forall\,\lambda>0$,
where $\bar{F}(x)=\int_{x}^{\infty}dx^{\prime}\,P_{x}(x^{\prime})$.

In Fig. 3, we show the shape of this class of probability distributions for
various choices of $n$.

Figure 3: Probability distribution functions Eq. 4.1 for various choices of
$n$ at fixed $\beta$ (left) and fixed variance $\sigma^{2}$ (right). We see
that smaller $n$ corresponds to heavier-tailed distributions.

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, $\mathcal{P}_{\zeta}\lesssim 10^{-4}$. In Fig. 4, we plot the
maximum PBH mass fraction at formation $\beta_{\rm max}$, for a variance that
saturates the spectral distortion constraints from COBE/FIRAS. Note that to
generate a present day abundance of $\Omega_{\rm BH}\gtrsim 10^{-20}$ with
$M_{\rm BH}\mathrel{\raise 1.29167pt\hbox{$<$\kern-7.5pt\lower
4.30554pt\hbox{$\sim$}}}10^{11}M_{\odot}$, we need a heavy tail with
$n\lesssim 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. [96, 97, 98, 99] 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. [47,
42] use the $\delta N$ formalism [100, 101, 102, 103, 104, 105] to compute the
non-perturbative distribution of curvature perturbations.151515Ref. [106]
presents another non-perturbative method for calculating the PBH abundance in
models where $\zeta$ is related to a Gaussian reference variable $\zeta_{G}$
by a “generalized local transformation”. For inflaton potentials with a small
step or bump-like feature that induces a period of off-attractor behavior,
these studies find that the tail of the distribution can become exponential,
without inducing any significant non-Gaussianities in the perturbative regime
[42]. 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 [101, 107] may be necessary for a proper
description of the dynamics. Combining this with the $\delta N$ formalism
[108], a number of studies [43, 45, 46, 109, 110, 111] have found that
prominent exponential tails arise generically from quantum diffusion.

While many single-field models have been found to yield exponential tails
[112], there is currently no known model which generates a heavier-tailed
distribution.161616Ref. [113] interprets the NANOGrav signal as evidence of
PBH mergers with $M_{\rm BH}\sim 10^{11}-10^{12}M_{\odot}$, and claims that
$\mu$-distortion constraints can be overcome for sufficiently non-Gaussian
single-field models. However they make no reference to $y$-type distortions,
which are more constraining for this mass range, and which we have verified
rule out this single-field scenario for any non-negligible abundance. However,
as shown in Fig. 4, a heavy tail
$P_{\delta}\sim\exp\left(-|\delta|^{n}\right)$ with $n\lesssim 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, the necessary
heavy tails can arise in certain multi-field scenarios. As we show in the next
section, the non-minimal curvaton model supplemented with self-interactions
can generate such a heavy-tailed distribution.

Figure 4: The maximum primordial black hole mass fraction at formation
$\beta_{\rm max}$ as a function of mass $M_{\rm BH}$ for a value of the
variance $\sigma_{\zeta}^{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_{\rm BH}\sim
3\times 10^{9}\,M_{\odot}$ correspond to the value of $M_{\rm BH}$ at which
$\mu$\- and $y$-type spectral distortions are equally restrictive. Contours of
constant $\Omega_{\rm BH}$, as computed using Eq. (2.2), are shown in dashed
gray.

##  5 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 [114, 115, 116]. The
perturbations of the curvaton are initially isocurvature, but become adiabatic
upon curvaton decay sometime after inflation ends [117]. Due to the non-
linearity inherent in this transfer, the full perturbation $\zeta$ 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 $\zeta$.

Non-Gaussianity in the curvaton model was first investigated using the $\delta
N$ formalism in Ref. [118], with implications for PBH formation analyzed in
Refs. [54, 119]. 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.

###  5.1 Curvaton Cosmology

We begin by reviewing the calculation of the curvature perturbation $\zeta$
and its statistics in the standard curvaton scenario [118], for which the
total potential is

$V(\phi,\chi)=V_{\phi}(\phi)+V_{\chi}(\chi)\,,$ (5.1)

where $V_{\phi}$ is the unspecified potential of the inflaton $\phi$, and
$V_{\chi}=m_{\chi}^{2}\chi^{2}/2$ is the quadratic potential of the curvaton
$\chi$. The curvaton mass is required to be light, satisfying $m_{\chi}\ll H$
throughout inflation, such that quantum perturbations are the dominant
influence on its evolution. This also implies that the background field
$\bar{\chi}$ will remain effectively fixed at some initial value
$\bar{\chi}_{*}$ 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
$\bar{\chi}_{*}\ll\sqrt{2V_{\phi}/m_{\chi}^{2}}$. Just like the inflaton, the
curvaton receives perturbations $\delta\chi_{*}\simeq H_{*}/2\pi$ set by the
Hubble rate at horizon exit $H_{*}$. Since the curvaton is a weakly coupled
field, we expect the perturbations $\delta\chi_{*}$ 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:

$\chi_{*}=\bar{\chi}_{*}+\delta\chi_{*}\,.$ (5.2)

The goal of this section will be to relate these initial Gaussian field
perturbations to the total curvature perturbation $\zeta$ via some mapping
$\zeta=f(\delta\chi_{*})$. This will be the key to constructing the
probability distribution function for $\zeta$, 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 $\rho_{R}\sim 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_{\chi}$, causing the curvaton to start oscillating about the minimum
of its potential. We denote the field value at which this occurs by
$\bar{\chi}_{0}$. During this oscillating phase, the curvaton redshifts like
matter with $\rho_{\chi}\propto a^{-3}$ and its energy density grows linearly
relative to radiation. Finally, when $H\sim\tau^{-1}$, where $\tau$ is the
$\chi$ lifetime, the curvaton decays to radiation and its isocurvature
perturbations become adiabatic perturbations, assuming the decay products
thermalize with the existing radiation.

###  5.2 The $\delta N$ Formalism

To calculate the distribution of the curvature perturbations in this model, we
employ the $\delta N$ formalism [100, 101, 102, 103, 104, 105, 120], which we
review in Appendix A. This technique identifies $\zeta$ on large scales ($k\ll
aH$) with the variation of inflationary $e$-folds across Hubble patches and
non-perturbatively captures its non-Gaussianities. The $\delta N$ formalism
was first used to study non-Gaussianity in curvaton models in Ref. [121, 118].
On a general hypersurface of uniform curvaton density, the conserved curvaton
curvature perturbation $\zeta_{\chi}$ is [121, 104]

$\zeta_{\chi}(t,\vec{x})=\delta
N(t,\vec{x})+\frac{1}{3}\ln\left(\frac{\rho_{\chi}(t,\vec{x})}{\bar{\rho}_{\chi}(t)}\right)\,,$
(5.3)

where $\delta N(t,\vec{x})$ is the perturbed number of $e$-folds,
$\rho_{\chi}(t,\vec{x})$ is the $\chi$ energy density, and
$\bar{\rho}_{\chi}(t)$ is its background value. In spatially flat slicing,
this becomes

$\zeta_{\chi}(t,\vec{x})=\frac{1}{3}\ln\left(\frac{\rho_{\chi}(t,\vec{x})}{\bar{\rho}_{\chi}(t)}\right)\,,$
(5.4)

and the curvaton energy density can be written as

$\rho_{\chi}(t,\vec{x})=e^{3\zeta_{\chi}(t,\vec{x})}\,\bar{\rho}_{\chi}(t)\,.$
(5.5)

In uniform total density slicing, Eq. (5.3) becomes

$\zeta_{\chi}(t,\vec{x})=\zeta+\frac{1}{3}\ln\left(\frac{\rho_{\chi}^{(u)}(t,\vec{x})}{\bar{\rho}_{\chi}^{(u)}(t)}\right)\,,$
(5.6)

where $\zeta$ is the total curvature perturbation.171717Note that this
generically has non-vanishing mean, $\expectationvalue{\zeta}\neq 0$, and so
when we later consider PBH formation, we will have to define a physical
$\zeta_{\rm phys}\equiv\zeta-\expectationvalue{\zeta}$ [54]. The expectation
value $\expectationvalue{\zeta}$ can be computed using the Gaussian
$P_{\delta\chi}$ as $\expectationvalue{\zeta}=\int
d\delta_{\chi}\,\zeta\,P_{\delta_{\chi}}$, with $\zeta$ expressed as a
function of $\delta_{\chi}$ given in Eq. (5.13). Since $\zeta$ and
$\zeta_{\chi}$ are gauge invariant quantities, Eq. (5.5) can be equated with
Eq. (5.6) to yield

$\rho_{\chi}^{(u)}(t,\vec{x})=e^{3(\zeta_{\chi}-\zeta)}\bar{\rho}_{\chi}^{(u)}(t)\,.$
(5.7)

A similar treatment for the radiation energy density in uniform total density
slicing gives

$\rho_{R}^{(u)}(t,\vec{x})=e^{4(\zeta_{R}-\zeta)}\bar{\rho}_{R}^{(u)}(t)\simeq
e^{-4\zeta}\bar{\rho}_{R}^{(u)}(t)\,,$ (5.8)

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=\tau^{-1}$, where $\tau$ is the curvaton lifetime. On
a uniform total density slice at $t=\tau$, the energy densities satisfy

$\rho_{R}^{(u)}(\tau,\vec{x})+\rho_{\chi}^{(u)}(\tau,\vec{x})=\bar{\rho}^{(u)}(\tau)\,,$
(5.9)

where $\bar{\rho}^{(u)}=\bar{\rho}^{(u)}_{R}+\bar{\rho}^{(u)}_{\chi}$ is the
total homogeneous energy density. Substituting Eqs. (5.7) and (5.8), this
becomes a 4th degree algebraic equation for $\zeta$ at $\tau$:

$e^{4\zeta}-\left(e^{3\zeta_{\chi}}\,\Omega_{\chi,{\tau}}\right)e^{\zeta}+\left(\Omega_{\chi,{\tau}}-1\right)=0\,.$
(5.10)

Alternatively, it is customary to introduce the parameter $r_{\tau}$, defined
as [114]:

$r_{\tau}=\frac{3\Omega_{\chi,{\tau}}}{4-\Omega_{\chi,{\tau}}}=\frac{3\bar{\rho}_{\chi}^{(u)}}{3\bar{\rho}_{\chi}^{(u)}+4\bar{\rho}_{R}^{(u)}}\bigg{|}_{\tau}\,,$
(5.11)

in terms of which the equation for $\zeta$ becomes:

$e^{4\zeta}-\frac{4r_{\tau}}{3+r_{\tau}}\left(e^{3\zeta_{\chi}}\right)e^{\zeta}+\frac{3r_{\tau}-3}{3+r_{\tau}}=0\,.$
(5.12)

The general solution is

$\zeta=\ln
X\,,\,\,\,X=\frac{B^{1/2}+\sqrt{Ar_{\tau}B^{-1/2}-B}}{(3+r_{\tau})^{1/3}}\,,$
(5.13)

where $A=e^{3\zeta_{\chi}}$ and we have defined

$\displaystyle B$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left[C^{1/3}+(r_{\tau}-1)(3+r_{\tau})^{1/3}C^{-1/3}\right]$
(5.14) $\displaystyle C$ $\displaystyle=$
$\displaystyle(Ar_{\tau})^{2}+\sqrt{(Ar_{\tau})^{4}+(3+r_{\tau})(1-r_{\tau})^{3}}\,.$
(5.15)

This gives the mapping $\zeta=\ln[X(\zeta_{\chi})]$ between $\zeta$ and
$\zeta_{\chi}$.

###  5.3 Calculating the Probability Distribution

Finally, to obtain $\zeta_{\chi}$ in terms of $\delta\chi_{*}$, 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
$\rho_{\chi}(t,\vec{x})=m_{\chi}^{2}\chi^{2}/2$. Expanding
$\chi(t,\vec{x})=\bar{\chi}(t)+\delta\chi(t,\vec{x})$, we can write this as

$\rho_{\chi}(t,\vec{x})=\frac{1}{2}m_{\chi}^{2}\bar{\chi}^{2}\left(1+\frac{\delta\chi}{\bar{\chi}}\right)^{2}=\bar{\rho}_{\chi}(t)(1+\delta_{\chi})^{2}\,,$
(5.16)

where $\delta_{\chi}=\delta\chi/\bar{\chi}$ is the curvaton contrast in
spatially flat slicing. Comparing with
$\rho_{\chi}=e^{3\zeta_{\chi}}\bar{\rho}_{\chi}$, we see we should identify
$e^{3\zeta_{\chi}}=(1+\delta_{\chi})^{2}$. Finally in order to relate the
contrast $\delta_{\chi}$ to its value at horizon exit, consider the equations
of motion for $\chi$ and the perturbation $\delta\chi$:

$\frac{d^{2}\bar{\chi}}{dt^{2}}+3H\frac{d\bar{\chi}}{dt}+m_{\chi}^{2}\bar{\chi}=0\,,$
(5.17a)
$\frac{d^{2}}{dt^{2}}(\delta\chi)+3H\frac{d}{dt}(\delta\chi)+\left(\frac{k^{2}}{a^{2}}+m_{\chi}^{2}\right)\delta\chi=0\,.$
(5.17b)

On superhorizon scales $k\ll aH$ is negligible and so these reduce to the same
equation. This implies $\delta\chi\sim\bar{\chi}$, such that
$\delta\chi/\bar{\chi}=\delta\chi_{*}/\bar{\chi}_{*}$ and

$A=e^{3\zeta_{\chi}}=(1+\delta_{\chi})^{2}=\left(1+\frac{\delta\chi_{*}}{\bar{\chi}_{*}}\right)^{2}\,.$
(5.18)

Combining with Eq. (5.13), we obtain the desired mapping between $\zeta$ and
the Gaussian initial field perturbations $\delta\chi_{*}$. We can now use
probability conservation to write

$P_{\zeta}[\zeta]=P_{\delta_{\chi}}\left[\delta_{\chi}^{+}(\zeta)\right]\bigg{|}\frac{d\delta_{\chi}^{+}}{d\zeta}\bigg{|}+P_{\delta_{\chi}}\left[\delta_{\chi}^{-}(\zeta)\right]\bigg{|}\frac{d\delta_{\chi}^{-}}{d\zeta}\bigg{|}\,,$
(5.19)

where $P_{\delta_{\chi}}$ is fully determined by the Gaussian variance
$\sigma_{0}^{2}$ and the roots $\delta^{\pm}_{\chi}(\zeta)$ satisfy

$\delta_{\chi}^{\pm}=-1\pm\sqrt{\left(\frac{3+r_{\tau}}{4r_{\tau}}\right)e^{3\zeta}+\left(\frac{3r_{\tau}-3}{4r_{\tau}}\right)e^{-\zeta}}\,,$
(5.20)

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_{\tau}$, which controls the heaviness of
the tail.

Figure 5: The probability distribution function for the curvature perturbation
$\zeta$ as given by Eq. (5.19) for various choices of $r_{\tau}$. The vertical
dashed line corresponds to the threshold value for collapse $\zeta_{c}$, and
the PBH mass fraction $\beta$ is obtained by integrating beyond this
threshold. We see that a smaller $r_{\tau}$ corresponds to a heavier-tailed
distribution, leading to a larger PBH abundance. Note that the reference value
$\sigma^{2}_{\zeta}=0.04$ is chosen to illustrate the heaviness of the tail as
$r_{\tau}$ 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.

##  6 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. These goals
can be achieved with 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 [54]. In Sec. 6.1 we review this scenario and calculate the power
spectrum that results 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 $\chi$
perturbations between horizon exit and the onset of curvaton oscillations.
While such interactions were studied perturbatively in Refs. [122, 123, 124,
125, 126, 127, 128, 129, 130], in Sec. 6.2 we go beyond these works to compute
the non-perturbative probability distribution function.

A curvaton model with these additional ingredients can generate an appreciable
SMBH population from the direct collapse of inflationary perturbations.

###  6.1 Non-Minimal Curvaton Scenario

The minimal curvaton model described in Sec. 5 provides non-Gaussian
statistics, but does not amplify $\mathcal{P}_{\zeta}$. This deficiency can be
remedied with a non-canonical kinetic term181818Such a term naturally arises
in many dilatonic and axionic models of inflation [131]. for the curvaton
[54]:

$\mathcal{L}\supset\frac{1}{2}f(\phi)^{2}(\partial\chi)^{2}\,.$ (6.1)

If $f(\phi)$ is chosen such that this kinetic term is suppressed at field
values $\phi=\phi_{*}$, the power spectrum will be enhanced on scales
corresponding to the horizon size at $\phi_{*}$. For concreteness, consider
the evolution of the inflaton and curvaton governed by the following system of
equations:

$\displaystyle\ddot{\phi}+3H\dot{\phi}+V^{\prime}(\phi)=ff^{\prime}\dot{\bar{\chi}}^{2}\simeq
0\,,~{}~{}\ddot{\bar{\chi}}+\left(3H+\frac{2f^{\prime}}{f}\dot{\phi}\right)\dot{\bar{\chi}}+\frac{m_{\chi}^{2}}{f^{2}}\bar{\chi}=0\,,$
(6.2)

where $f^{\prime}\equiv\partial_{\phi}f$ and the source term of the first
equation is negligible until the curvaton begins to oscillate. Similarly, the
curvaton perturbation evolves according to

$\displaystyle\frac{d^{2}}{dt^{2}}(\delta\chi)+\left(3H+\frac{2f^{\prime}}{f}\dot{\phi}\right)\frac{d}{dt}(\delta\chi)+\left(\frac{k^{2}}{a^{2}}+\frac{m_{\chi}^{2}}{f^{2}}\right)\delta\chi\simeq
0\,,$ (6.3)

which, to leading order, can be simplified inside the horizon, $k\gg aH$, to
yield

$\frac{d^{2}}{dt^{2}}\left(f\delta\chi\right)+\frac{k^{2}}{a^{2}}(f\delta\chi)\simeq
0\,,$ (6.4)

whose solution can be written as

$\delta\chi\simeq\frac{1}{\sqrt{2k}\,af}\exp\left(-ik\\!\int\frac{dt}{a}\right),$
(6.5)

which establishes the initial conditions in the adiabatic vacuum. As in Sec.
5, on superhorizon scales, $k\ll aH$, $\bar{\chi}$ and $\delta\chi$ evolve
according to the same equation, so their solutions have the same functional
form for all $t>t_{*}$:

$\frac{\delta\chi}{\bar{\chi}}=\frac{\delta\chi_{*}}{\bar{\chi}_{*}}\simeq\frac{H_{*}}{\sqrt{2k^{3}}f(\phi_{*})\bar{\chi}_{*}}\,,$
(6.6)

where $k=a(t_{*})H(t_{\star})$. Using Eq. (6.6), the power spectrum for
$\delta_{\chi}$ becomes

$\mathcal{P}_{\delta_{\chi}}(k)=\frac{k^{3}}{2\pi^{2}}\left|\frac{\delta\chi_{k}}{\bar{\chi}}\right|^{2}=\frac{1}{\bar{\chi}_{*}^{2}}\left(\frac{H_{*}}{2\pi
f(\phi_{*})}\right)^{2}~{},$ (6.7)

so, if $f(\phi)$ is chosen to have a dip at $\phi_{*}$, corresponding to
$k_{*}=k_{\rm BH}$, $\mathcal{P}_{\delta_{\chi}}$ will exhibit a peak at
$k_{\rm BH}$. However, for modes far away from $k_{\rm BH}$,
$f(\phi_{*})\approx 1$, recovering the nearly scale-invariant spectrum
required for consistency with CMB observations on larger scales. Assuming
$f(\phi)$ has such a localized feature, combining Eqs. (5.19), (6.7), and
(2.8), the PBH abundance becomes

$\displaystyle\beta=2\int_{\delta_{\chi,{c}}^{+}}^{\infty}d\delta_{\chi}\,P_{\delta_{\chi}}[\delta_{\chi}]+2\int_{-\infty}^{\delta_{\chi,{c}}^{-}}d\delta_{\chi}\,P_{\delta_{\chi}}[\delta_{\chi}]=\text{erfc}\left(\frac{\delta_{\chi,{c}}^{+}}{\sqrt{2}\sigma_{0}}\right)+\text{erfc}\left(\frac{|\delta_{\chi,{c}}^{-}|}{\sqrt{2}\sigma_{0}}\right)~{},$
(6.8)

where $\delta_{\chi,{c}}^{\pm}=\delta_{\chi}^{\pm}(\zeta_{c})$ are the roots
of Eq. (5.20) evaluated at the threshold.

Using Eq. (3.7), we obtain the maximum value of $\beta$ consistent with
spectral distortion constraints.191919Note that the variance
$\sigma_{\zeta}^{2}$ corresponds to the physical curvature perturbation,
$\sigma_{\zeta}^{2}=\langle\zeta_{\rm
phys}^{2}\rangle=\expectationvalue{\zeta^{2}}-\expectationvalue{\zeta}^{2}$,
which can be computed from $\zeta(\delta_{\chi})$ in Eq. (5.13) as
$\sigma_{\zeta}^{2}=\int
d\delta_{\chi}\,\zeta^{2}\,P_{\delta_{\chi}}-\left(\int
d\delta_{\chi}\,\zeta\,P_{\delta_{\chi}}\right)^{2}$.

Figure 6: The maximal PBH mass fraction at formation $\beta_{\rm max}$ in the
standard curvaton scenario for $\sigma_{\zeta}^{2}$ subject to spectral
distortion constraints for various choices of $r_{\tau}$ from Eq. (5.11)

The degree of non-Gaussianity in this scenario is governed by $r_{\tau}$, as
defined in Eq. (5.11). In the $r_{\tau}\rightarrow 1$ limit, the curvaton
dominates the energy density prior to its decay, so the relation between
$\zeta$ and $\zeta_{\chi}$ from Eq. (5.12) is approximately linear. Thus, in
this regime, the non-Gaussianity comes entirely from the non-linear
relationship between $\zeta_{\chi}$ and $\delta_{\chi}$; see Eq. (5.18). In
the opposite regime202020It may seem counterintuitive that a very subdominant
curvaton can still generate the curvature perturbation. The key is that
because the curvaton is subdominant during inflation, its perturbations from
quantum fluctuations are large relative to the background field value. Ref.
[132] finds that a subdominant curvaton can still viably produce the curvature
perturbation provided $r_{\tau}\gtrsim 10^{-3}$. that the curvaton is still
very subdominant when it decays ($r_{\tau}\ll 1$), the relation between
$\zeta$ and $\zeta_{\chi}$ 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.

###  6.2 Self-Interacting Curvaton Scenario

Introducing curvaton self-interactions allows for non-linear $\delta_{\chi}$
between horizon exit and the onset of oscillations, leading to even more
dramatic departures from Gaussianity.212121This 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 weak
self-interactions in curvaton models were first investigated in Refs. [122,
123, 124, 125, 126, 127, 128, 129, 130]. However, these studies restricted
themselves to computing the non-linearity parameters $f_{\rm NL}$ and $g_{\rm
NL}$, which do not capture the non-perturbative statistics in the tail of the
distribution. By employing the $\delta N$ formalism, we extend these studies
to compute the full probability distribution for $\zeta$.

We now allow the curvaton potential to be an arbitrary well-defined function
of $\chi$. In our setup, at $t=t_{\rm int}$, corresponding to
$V^{\prime}(\chi_{*})\sim H$, $\chi$ begins rolling towards the minimum of its
potential and its energy density is initially dominated by the interaction
terms. At a later time $t=t_{0}>t_{\rm int}$, the interaction terms become
subdominant, and the curvaton mass term drives field evolution, resulting in
matter-like scaling, $\rho_{\chi}\propto 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_{\tau}$ at the time of decay.

Recall that in the case of the quadratic potential, the curvaton density
contrast $\delta_{\chi}$ remained constant after horizon exit since the
background field $\bar{\chi}$ and perturbation $\delta\chi$ obeyed the same
equation of motion, shown in Eq. (5.17). This led to
$\delta\chi/\bar{\chi}=\delta\chi_{*}/\bar{\chi}_{*}$, which allowed
$\delta_{\chi}$ to be used as a Gaussian reference variable. Upon introducing
interactions this is no longer the case, as $\delta_{\chi}$ evolves non-
trivially between $t_{\rm osc}$ and curvaton decay at $t=\tau$. In the regime,
the equation of motion for the perturbation is:

$\displaystyle\frac{d^{2}}{dt^{2}}(\delta\chi)+3H\frac{d}{dt}(\delta\chi)+V^{\prime\prime}_{\chi}\delta\chi=0\,,$
(6.9)

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 $\delta\chi_{*}$, which can be related to the $\zeta_{\chi}$ by
first solving Eq. (6.9) along with the equation for the background field
$\bar{\chi}$, computing the total and background energy densities at decay,
and finally applying Eq. (5.4). The resulting
$\delta\chi_{*}^{i}=g_{i}(\zeta_{\chi})$ can then be mapped onto the total
curvature perturbation $\zeta$ via Eq. (5.12), since the relationship between
$\zeta_{\chi}$ and $\zeta$ is unchanged in the presence of $\chi$ 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 $\delta\chi_{*}$ and $\zeta_{\chi}$
at the time of curvaton decay. For $\chi$ values sufficiently close to the
minimum of its potential, the potential is approximately quadratic, and the
energy density is

$\rho_{\chi}\simeq\frac{1}{2}m_{\chi}^{2}\chi_{0}^{2}\,,$ (6.10)

where $\chi_{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 $\chi_{*})$. In terms
of background field values $\bar{\chi}_{0}\equiv\chi_{0}(\bar{\chi}_{*})$,
this can be expanded as:

$\chi_{0}=\bar{\chi}_{0}+\sum_{n=1}\frac{1}{n!}\bar{\chi}_{0}^{(n)}\delta\chi_{*}^{n}~{}~{},~{}~{}\bar{\chi}_{0}^{(n)}\equiv\frac{\partial^{n}\chi_{0}}{\partial\chi_{*}^{n}}\Big{|}_{\chi=\bar{\chi}_{*}},$
(6.11)

and the energy density can then be written as:

$\rho_{\chi}\simeq\bar{\rho}_{\chi}\left[1+\frac{1}{\bar{\chi}_{0}}\sum_{n=1}\frac{1}{n!}\bar{\chi}_{0}^{(n)}\delta\chi_{*}^{n}\right]^{2}\,,$
(6.12)

where $\bar{\rho}_{\chi}=\frac{1}{2}m_{\chi}^{2}\bar{\chi}_{0}^{2}$. Comparing
with Eq. (5.4), the bracketed quantity is identified with $e^{3\zeta_{\chi}}$.
Then also expanding

$\displaystyle\zeta_{\chi}=\zeta_{\chi,1}+\sum^{\infty}_{n=2}\frac{1}{n!}\zeta_{\chi,n}~{}~{}~{},$
(6.13)

we can write to leading order:

$\displaystyle\zeta_{\chi}(\delta\chi_{*})=\frac{2}{3}\left(\frac{\bar{\chi}^{\prime}_{0}}{\bar{\chi}_{0}}\right)\delta\chi_{*}+\frac{1}{3}\left(\frac{\bar{\chi}_{0}\bar{\chi}^{\prime\prime}_{0}}{{\bar{\chi}}^{\prime\,2}_{0}}\right)\left(\frac{\bar{\chi}^{\prime}_{0}}{\bar{\chi}_{0}}\right)^{2}\delta\chi_{*}^{2}+\frac{2}{9}\left(1-\frac{3}{2}\frac{\bar{\chi}_{0}\bar{\chi}^{\prime\prime}_{0}}{{\bar{\chi}}^{\prime\,2}_{0}}+\frac{1}{2}\frac{\bar{\chi}_{0}^{2}\bar{\chi}^{\prime\prime\prime}_{0}}{{\bar{\chi}}^{\prime\,3}_{0}}\right)\left(\frac{\bar{\chi}^{\prime}_{0}}{\bar{\chi}_{0}}\right)^{3}\delta\chi_{*}^{3}\,,~{}~{}~{}~{}~{}$
(6.14)

which can then be substituted into Eq. (5.13) to obtain $\zeta$ as a function
of the Gaussian reference variable $\delta\chi_{*}$. It is also possible to
invert this mapping to obtain the roots
$g_{i}(\zeta)=\delta\chi_{*}^{i}(\zeta)$. The probability distribution is
then:

$P_{\zeta}[\zeta]=\sum_{j}\left|\frac{g_{j}(\zeta)}{d\zeta}\right|P_{\delta\chi_{*}}[g_{j}(\zeta)]\,,$
(6.15)

where the sum runs over all real roots.

The resulting maximal PBH abundance consistent with spectral distortion
constraints is shown in Fig. 7 for a sample set of parameters. Comparison with
$\beta_{\rm max}$ in the standard curvaton scenario (see Fig. 6) 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 ($r_{\tau}\ll 1$) to yield an appreciable fraction of primordial
SMBHs while evading spectral distortion constraints.222222It is actually even
more natural to fulfill this condition here since energy density in the
curvaton field dilutes more quickly upon introducing interactions.

Figure 7: Maximal PBH mass fraction at formation $\beta_{\rm max}$ as a
function of PBH mass in the self-interacting curvaton model from Sec. 6.2,
where $\sigma_{\zeta}^{2}$ saturates spectral distortion bounds for different
values of $r_{\tau}$, as defined in Eq. (5.11). We take
$\zeta_{\chi}=c_{1}\,\delta\chi_{*}+c_{2}\,\delta\chi_{*}^{2}+c_{3}\,\delta\chi_{*}^{3}$
to take the perturbative form of Eq. (6.14), with sample parameters fixed as
$(c_{1},c_{2},c_{3})=(1,0.5,0.1)$. For a concrete potential which can realize
these values, see Appendix B. Contours of constant $\Omega_{\rm BH}$ today are
shown in dashed gray.

##  7 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_{\rm BH}\sim 10^{6}-10^{11}M_{\odot}$ 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\sim{\rm 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_{\zeta}\sim\exp\left(-|\zeta|^{n}\right)$ with $n\lesssim
0.6$; we are not aware of any single-field inflationary model that can realize
this behavior.

In this paper, we have shown that such heavy-tailed distributions can be
generated if a self-interacting curvaton field is also present as a spectator
during inflation. In the standard curvaton scenario, non-Gaussianity 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. By adding quartic self-interactions, we have found that the
non-Gaussian tail can be made significantly heavier, allowing for a viable
primordial SMBH population, as shown in Fig. 7.

Finally, this work is especially timely in light of recent pulsar timing array
observations [133, 134, 135]. In particular, the NANOGrav collaboration just
announced evidence of a signal in its 15-year data set consistent with a
stochastic gravitational wave background in the nHz frequency range [136]. 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 [137, 138].
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 [139], a highly clustered population could still
viably source the gravitational wave signal [140]. 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. [141,
142] for closely related work. There are also other signals, such as scalar-
induced secondary gravitational waves, which could offer complimentary
evidence of this scenario and deserve further study; see Refs. [143, 144] for
a potential connection between scalar-induced gravitational waves and the PTA
signal. Regardless, the recent detection of the gravitational wave background
provides us with motivation to better understand the cosmic origin of our
universe’s supermassive black holes.

## Acknowledgments

We would like to thank Wayne Hu and Keisuke Inomata for helpful conversations.
AI thanks the organizers and participants of the New Horizons in Primordial
Black Hole Physics (NEHOP) workshop, where part of this work was completed.
DH, GK and AS are supported by the Fermi Research Alliance, LLC under Contract
No. DE-AC02-07CH11359 with the U.S. Department of Energy, Office of Science,
Office of High Energy Physics.

## References

  * [1] Y. Shen, G.T. Richards, M.A. Strauss, P.B. Hall, D.P. Schneider, S. Snedden et al., _A catalog of quasar properties from sloan digital sky survey data release 7_ , _The Astrophysical Journal Supplement Series_ 194 (2011) 45.
  * [2] E. Bañados, B.P. Venemans, C. Mazzucchelli, E.P. Farina, F. Walter, F. Wang et al., _An 800-million-solar-mass black hole in a significantly neutral Universe at a redshift of 7.5_ , _Nature_ 553 (2018) 473 [1712.01860].
  * [3] N.P. Ross and N.J.G. Cross, _The near and mid-infrared photometric properties of known redshift $z\geq 5$ quasars_, _MNRAS_ 494 (2020) 789 [1906.06974].
  * [4] C.J. Willott et al., _The Canada-France High-z Quasar Survey: nine new quasars and the luminosity function at redshift 6_ , _Astron. J._ 139 (2010) 906 [0912.0281].
  * [5] E. Bañados et al., _The Pan-STARRS1 distant $z>5.6$ quasar survey: more than 100 quasars within the first Gyr of the universe_, _Astrophys. J. Suppl._ 227 (2016) 11 [1608.03279].
  * [6] L. Jiang, I.D. McGreer, X. Fan, M.A. Strauss, E. Bañados, R.H. Becker et al., _The Final SDSS High-redshift Quasar Sample of 52 Quasars at $z>5.7$_, _APJ_ 833 (2016) 222 [1610.05369].
  * [7] DES, WISE, VISTA collaboration, _Eight new luminous z $\geq$ 6 quasars discovered via SED model fitting of VISTA, WISE and Dark Energy Survey Year 1 observations_, _Mon. Not. Roy. Astron. Soc._ 468 (2017) 4702 [1701.04852].
  * [8] Y. Matsuoka, K. Iwasawa, M. Onoue, N. Kashikawa, M.A. Strauss, C.-H. Lee et al., _Subaru High-z Exploration of Low-luminosity Quasars (SHELLQs). X. Discovery of 35 Quasars and Luminous Galaxies at $5.7\leq z\leq 7.0$_, _APJ_ 883 (2019) 183 [1908.07910].
  * [9] J. Yang, F. Wang, X. Fan, X.-B. Wu, F. Bian, E. Bañados et al., _Filling in the Quasar Redshift Gap at $z\sim 5.5$ II: A Complete Survey of Luminous Quasars in the Post-Reionization Universe_, _arXiv e-prints_ (2018) arXiv:1810.11927 [1810.11927].
  * [10] F. Wang, J. Yang, X. Fan, X.-B. Wu, M. Yue, J.-T. Li et al., _Exploring Reionization-era Quasars. III. Discovery of 16 Quasars at $6.4\lesssim z\lesssim 6.9$ with DESI Legacy Imaging Surveys and the UKIRT Hemisphere Survey and Quasar Luminosity Function at $z\sim 6.7$_, _APJ_ 884 (2019) 30 [1810.11926].
  * [11] B. Mutlu-Pakdil, M.S. Seigar and B.L. Davis, _The local black holes mass function derived from the $m_{BH}-p$ and the $m_{BH}-n$ relations_, _The Astrophysical Journal_ 830 (2016) 117.
  * [12] F. Shankar, D.H. Weinberg and J. Miralda-Escudé, _Accretion-driven evolution of black holes: Eddington ratios, duty cycles and active galaxy fractions_ , _Monthly Notices of the Royal Astronomical Society_ 428 (2012) 421.
  * [13] X. Shen, P.F. Hopkins, C.-A. Faucher-Giguère, D.M. Alexander, G.T. Richards, N.P. Ross et al., _The bolometric quasar luminosity function at $z=0-7$_, _Monthly Notices of the Royal Astronomical Society_ 495 (2020) 3252.
  * [14] A. Sicilia, A. Lapi, L. Boco, F. Shankar, D.M. Alexander, V. Allevato et al., _The black hole mass function across cosmic time. II. heavy seeds and (super)massive black holes_ , _The Astrophysical Journal_ 934 (2022) 66.
  * [15] N. Banik, J.C. Tan and P. Monaco, _The formation of supermassive black holes from Population III.1 seeds. I. Cosmic formation histories and clustering properties_ , _Monthly Notices of the Royal Astronomical Society_ 483 (2018) 3592.
  * [16] E.E. Salpeter, _Accretion of Interstellar Matter by Massive Objects._ , _apj_ 140 (1964) 796.
  * [17] M. Volonteri, _Formation of Supermassive Black Holes_ , _Astron. Astrophys. Rev._ 18 (2010) 279 [1003.4404].
  * [18] Z. Haiman, _Constraints from gravitational recoil on the growth of supermassive black holes at high redshift_ , _Astrophys. J._ 613 (2004) 36 [astro-ph/0404196].
  * [19] S.L. Shapiro, _Spin, accretion and the cosmological growth of supermassive black holes_ , _Astrophys. J._ 620 (2005) 59 [astro-ph/0411156].
  * [20] M. Volonteri and M.J. Rees, _Quasars at z=6: The survival of the fittest_ , _Astrophys. J._ 650 (2006) 669 [astro-ph/0607093].
  * [21] T. Tanaka and Z. Haiman, _The Assembly of Supermassive Black Holes at High Redshifts_ , _Astrophys. J._ 696 (2009) 1798 [0807.4702].
  * [22] V. De Luca, G. Franciolini and A. Riotto, _Heavy Primordial Black Holes from Strongly Clustered Light Black Holes_ , _Phys. Rev. Lett._ 130 (2023) 171401 [2210.14171].
  * [23] SDSS collaboration, _A Survey of z > 5.8 quasars in the Sloan Digital Sky Survey I: Discovery of three new quasars and the spatial density of luminous quasars at $z\sim 6$_, _Astron. J._ 122 (2001) 2833 [astro-ph/0108063].
  * [24] D.J. Mortlock et al., _A luminous quasar at a redshift of z = 7.085_ , _Nature_ 474 (2011) 616 [1106.6088].
  * [25] X.-B. Wu, F. Wang, X. Fan, W. Yi, W. Zuo, F. Bian et al., _An ultraluminous quasar with a twelve-billion-solar-mass black hole at redshift 6.30_ , _Nature_ 518 (2015) 512 [1502.07418].
  * [26] B. Trakhtenbrot, H. Netzer, P. Lira and O. Shemmer, _Black-Hole Mass and Growth Rate at $z\simeq 4.8$: A Short Episode of Fast Growth Followed by Short Duty Cycle Activity_, _Astrophys. J._ 730 (2011) 7 [1012.1871].
  * [27] H. Netzer, _The Largest black holes and the most luminous galaxies_ , _Astrophys. J. Lett._ 583 (2003) L5 [astro-ph/0210548].
  * [28] P. Natarajan and E. Treister, _Is there an upper limit to black hole masses?_ , _Mon. Not. Roy. Astron. Soc._ 393 (2009) 838 [0808.2813].
  * [29] K. Inayoshi and Z. Haiman, _Is There a Maximum Mass for Black Holes in Galactic Nuclei?_ , _APJ_ 828 (2016) 110 [1601.02611].
  * [30] J. Liu, L. Bian, R.-G. Cai, Z.-K. Guo and S.-J. Wang, _Primordial black hole production during first-order phase transitions_ , _Phys. Rev. D_ 105 (2022) L021303 [2106.05637].
  * [31] J. Garriga, A. Vilenkin and J. Zhang, _Black holes and the multiverse_ , _JCAP_ 02 (2016) 064 [1512.01819].
  * [32] H. Davoudiasl, P.B. Denton and J. Gehrlein, _Supermassive Black Holes, Ultralight Dark Matter, and Gravitational Waves from a First Order Phase Transition_ , _Phys. Rev. Lett._ 128 (2022) 081101 [2109.01678].
  * [33] H.-L. Huang, Y. Cai, J.-Q. Jiang, J. Zhang and Y.-S. Piao, _Supermassive primordial black holes in multiverse: for nano-Hertz gravitational wave and high-redshift JWST galaxies_ , 2306.17577.
  * [34] B.J. Carr, _The Primordial black hole mass spectrum_ , _Astrophys. J._ 201 (1975) 1.
  * [35] S.M. Leach, M. Sasaki, D. Wands and A.R. Liddle, _Enhancement of superhorizon scale inflationary curvature perturbations_ , _Phys. Rev. D_ 64 (2001) 023512 [astro-ph/0101406].
  * [36] J. Garcia-Bellido and E. Ruiz Morales, _Primordial black holes from single field models of inflation_ , _Phys. Dark Univ._ 18 (2017) 47 [1702.03901].
  * [37] C. Germani and T. Prokopec, _On primordial black holes from an inflection point_ , _Phys. Dark Univ._ 18 (2017) 6 [1706.04226].
  * [38] G. Ballesteros and M. Taoso, _Primordial black hole dark matter from single field inflation_ , _Phys. Rev. D_ 97 (2018) 023501 [1709.05565].
  * [39] M.P. Hertzberg and M. Yamada, _Primordial Black Holes from Polynomial Potentials in Single Field Inflation_ , _Phys. Rev. D_ 97 (2018) 083509 [1712.09750].
  * [40] S.S. Mishra and V. Sahni, _Primordial black holes from a tiny bump/dip in the inflaton potential_ , _JCAP_ 2020 (2019) 007  [1911.00057].
  * [41] K. Inomata, E. McDonough and W. Hu, _Amplification of primordial perturbations from the rise or fall of the inflaton_ , _JCAP_ 02 (2022) 031 [2110.14641].
  * [42] Y.-F. Cai, X. Ma, M. Sasaki, D.-G. Wang and Z. Zhou, _Highly non-Gaussian tails and primordial black holes from single-field inflation_ , _JCAP_ 12 (2022) 034 [2207.11910].
  * [43] C. Pattison, V. Vennin, H. Assadullahi and D. Wands, _Quantum diffusion during inflation and primordial black holes_ , _JCAP_ 10 (2017) 046 [1707.00537].
  * [44] M. Biagetti, G. Franciolini, A. Kehagias and A. Riotto, _Primordial Black Holes from Inflation and Quantum Diffusion_ , _JCAP_ 07 (2018) 032 [1804.07124].
  * [45] J.M. Ezquiaga, J. Garcia-Bellido and V. Vennin, _The exponential tail of inflationary fluctuations: consequences for primordial black holes_ , _JCAP_ 03 (2020) 029 [1912.05399].
  * [46] C.A. Pattison, V. Vennin, D. Wands and H. Assadullahi, _Ultra-slow-roll inflation with quantum diffusion_ , _Journal of Cosmology and Astroparticle Physics_ 2021 (2021) .
  * [47] Y.-F. Cai, X.-H. Ma, M. Sasaki, D.-G. Wang and Z. Zhou, _One small step for an inflaton, one giant leap for inflation: A novel non-Gaussian tail and primordial black holes_ , _Phys. Lett. B_ 834 (2022) 137461 [2112.13836].
  * [48] J. Garcia-Bellido, A.D. Linde and D. Wands, _Density perturbations and black hole formation in hybrid inflation_ , _Phys. Rev. D_ 54 (1996) 6040 [astro-ph/9605094].
  * [49] S. Clesse and J. Garcia-Bellido, _Massive primordial black holes from hybrid inflation as dark matter and the seeds of galaxies_ , _Phys. Rev. D_ 92 (2015) 023524.
  * [50] A.R. Brown, _Hyperbolic Inflation_ , _Phys. Rev. Lett._ 121 (2018) 251601 [1705.03023].
  * [51] G.A. Palma, S. Sypsas and C. Zenteno, _Seeding primordial black holes in multifield inflation_ , _Phys. Rev. Lett._ 125 (2020) 121301 [2004.06106].
  * [52] J. Fumagalli, S. Renaux-Petel, J.W. Ronayne and L.T. Witkowski, _Turning in the landscape: A new mechanism for generating primordial black holes_ , _Phys. Lett. B_ 841 (2023) 137921 [2004.08369].
  * [53] M. Braglia, D.K. Hazra, F. Finelli, G.F. Smoot, L. Sriramkumar and A.A. Starobinsky, _Generating PBHs and small-scale GWs in two-field models of inflation_ , _JCAP_ 08 (2020) 001 [2005.02895].
  * [54] S. Pi and M. Sasaki, _Primordial Black Hole Formation in Non-Minimal Curvaton Scenario_ , 2112.12680.
  * [55] S.R. Geller, W. Qin, E. McDonough and D.I. Kaiser, _Primordial black holes from multifield inflation with nonminimal couplings_ , _Phys. Rev. D_ 106 (2022) 063535 [2205.04471].
  * [56] O. Özsoy and G. Tasinato, _Inflation and Primordial Black Holes_ , _Universe_ 9 (2023) 203 [2301.03600].
  * [57] Y. Cai, M. Zhu and Y.-S. Piao, _Primordial black holes from null energy condition violation during inflation_ , 2305.10933.
  * [58] Planck collaboration, _Planck 2018 results. VI. Cosmological parameters_ , _Astron. Astrophys._ 641 (2020) A6 [1807.06209].
  * [59] J.C. Mather et al., _Measurement of the Cosmic Microwave Background spectrum by the COBE FIRAS instrument_ , _Astrophys. J._ 420 (1994) 439.
  * [60] D.J. Fixsen, E.S. Cheng, J.M. Gales, J.C. Mather, R.A. Shafer and E.L. Wright, _The Cosmic Microwave Background spectrum from the full COBE FIRAS data set_ , _Astrophys. J._ 473 (1996) 576 [astro-ph/9605054].
  * [61] Planck collaboration, _Planck 2018 results. IX. Constraints on primordial non-Gaussianity_ , _Astron. Astrophys._ 641 (2020) A9 [1905.05697].
  * [62] Planck collaboration, _Planck 2018 results. X. Constraints on inflation_ , _Astron. Astrophys._ 641 (2020) A10 [1807.06211].
  * [63] J.M. Bardeen, J.R. Bond, N. Kaiser and A.S. Szalay, _The Statistics of Peaks of Gaussian Random Fields_ , _Astrophys. J._ 304 (1986) 15.
  * [64] A.M. Green, A.R. Liddle, K.A. Malik and M. Sasaki, _A New calculation of the mass fraction of primordial black holes_ , _Phys. Rev. D_ 70 (2004) 041502 [astro-ph/0403181].
  * [65] S. Young and M. Musso, _Application of peaks theory to the abundance of primordial black holes_ , _JCAP_ 11 (2020) 022 [2001.06469].
  * [66] C.-M. Yoo, T. Harada, S. Hirano and K. Kohri, _Abundance of Primordial Black Holes in Peak Theory for an Arbitrary Power Spectrum_ , _PTEP_ 2021 (2021) 013E02 [2008.02425].
  * [67] W.H. Press and P. Schechter, _Formation of galaxies and clusters of galaxies by selfsimilar gravitational condensation_ , _Astrophys. J._ 187 (1974) 425.
  * [68] M. Shibata and M. Sasaki, _Black hole formation in the Friedmann universe: Formulation and computation in numerical relativity_ , _Phys. Rev. D_ 60 (1999) 084002 [gr-qc/9905064].
  * [69] A. Escrivà, C. Germani and R.K. Sheth, _Universal threshold for primordial black hole formation_ , _Phys. Rev. D_ 101 (2020) 044022 [1907.13311].
  * [70] V. Atal, J. Cid, A. Escrivà and J. Garriga, _PBH in single field inflation: the effect of shape dispersion and non-Gaussianities_ , _JCAP_ 05 (2020) 022 [1908.11357].
  * [71] V. Atal, J. Garriga and A. Marcos-Caballero, _Primordial black hole formation with non-Gaussian curvature perturbations_ , _JCAP_ 09 (2019) 073 [1905.13202].
  * [72] C.-M. Yoo, J.-O. Gong and S. Yokoyama, _Abundance of primordial black holes with local non-Gaussianity in peak theory_ , _JCAP_ 09 (2019) 033 [1906.06790].
  * [73] A. Kehagias, I. Musco and A. Riotto, _Non-Gaussian Formation of Primordial Black Holes: Effects on the Threshold_ , _JCAP_ 12 (2019) 029 [1906.07135].
  * [74] T. Harada, C.-M. Yoo and K. Kohri, _Threshold of primordial black hole formation_ , _Phys. Rev. D_ 88 (2013) 084051 [1309.4201].
  * [75] G. Ferrante, G. Franciolini, A. Iovino, Junior. and A. Urbano, _Primordial non-Gaussianity up to all orders: Theoretical aspects and implications for primordial black hole models_ , _Phys. Rev. D_ 107 (2023) 043520 [2211.01728].
  * [76] V. De Luca, A. Kehagias and A. Riotto, _How Well Do We Know the Primordial Black Hole Abundance? The Crucial Role of Non-Linearities when Approaching the Horizon_ , 2307.13633.
  * [77] S. Young, C.T. Byrnes and M. Sasaki, _Calculating the mass fraction of primordial black holes_ , _JCAP_ 07 (2014) 045 [1405.7023].
  * [78] K. Ando, K. Inomata and M. Kawasaki, _Primordial black holes and uncertainties in the choice of the window function_ , _Phys. Rev. D_ 97 (2018) 103528 [1802.06393].
  * [79] R.A. Sunyaev and Y.B. Zeldovich, _Small scale fluctuations of relic radiation_ , _Astrophys. Space Sci._ 7 (1970) 3.
  * [80] Y.B. Zeldovich and R.A. Sunyaev, _The Interaction of Matter and Radiation in a Hot-Model Universe_ , _Astrophys. Space Sci._ 4 (1969) 301.
  * [81] A.F. Illarionov and R.A. Sunyaev, _Why the number of galactic X-ray stars is so small?_ , _Astron. Astrophys._ 39 (1975) 185.
  * [82] W. Hu and J. Silk, _Thermalization and spectral distortions of the cosmic background radiation_ , _Phys. Rev. D_ 48 (1993) 485.
  * [83] J. Chluba, A.L. Erickcek and I. Ben-Dayan, _Probing the inflaton: Small-scale power spectrum constraints from measurements of the CMB energy spectrum_ , _Astrophys. J._ 758 (2012) 76 [1203.2681].
  * [84] R. Khatri and R.A. Sunyaev, _Beyond y and $\mu$: the shape of the cmb spectral distortions in the intermediate epoch, $1.5\times 10^{4}<z<2\times 10^{5}$_, _JCAP_ 09 (2012) 016 [1207.6654].
  * [85] J. Chluba and R.A. Sunyaev, _Superposition of blackbodies and the dipole anisotropy: A possibility to calibrate CMB experiments_ , _Astronomy & Astrophysics_ 424 (2004) 389.
  * [86] A. Stebbins, _Cmb spectral distortions from the scattering of temperature anisotropies_ , 2007.  10.48550/ARXIV.ASTRO-PH/0703541.
  * [87] J. Chluba and R.A. Sunyaev, _The evolution of CMB spectral distortions in the early Universe_ , _Mon. Not. Roy. Astron. Soc._ 419 (2012) 1294 [1109.6552].
  * [88] C. Pitrou and A. Stebbins, _Parameterization of temperature and spectral distortions in future CMB experiments_ , _General Relativity and Gravitation_ 46 (2014) .
  * [89] B.J. Carr and J.E. Lidsey, _Primordial black holes and generalized constraints on chaotic inflation_ , _Physical Review D_ 48 (1993) 543.
  * [90] J. Chluba, _Which spectral distortions does $\Lambda$CDM actually predict?_, _Mon. Not. Roy. Astron. Soc._ 460 (2016) 227 [1603.02496].
  * [91] A.M. Green and B.J. Kavanagh, _Primordial Black Holes as a dark matter candidate_ , _J. Phys. G_ 48 (2021) 043001 [2007.10722].
  * [92] S. Bird, H.V. Peiris, M. Viel and L. Verde, _Minimally Parametric Power Spectrum Reconstruction from the Lyman-alpha Forest_ , _Mon. Not. Roy. Astron. Soc._ 413 (2011) 1717 [1010.1519].
  * [93] C.T. Byrnes, P.S. Cole and S.P. Patil, _Steepest growth of the power spectrum and primordial black holes_ , _JCAP_ 06 (2019) 028 [1811.11158].
  * [94] S. Ando, N. Hiroshima and K. Ishiwata, _Constraining the primordial curvature perturbation using dark matter substructure_ , _Phys. Rev. D_ 106 (2022) 103014 [2207.05747].
  * [95] T. Nakama, T. Suyama and J. Yokoyama, _Supermassive black holes formed by direct collapse of inflationary perturbations_ , _Phys. Rev. D_ 94 (2016) 103522 [1609.02245].
  * [96] M.H. Namjoo, H. Firouzjahi and M. Sasaki, _Violation of non-Gaussianity consistency relation in a single field inflationary model_ , _EPL_ 101 (2013) 39001 [1210.3692].
  * [97] X. Chen, H. Firouzjahi, M.H. Namjoo and M. Sasaki, _A Single Field Inflation Model with Large Local Non-Gaussianity_ , _EPL_ 102 (2013) 59001 [1301.5699].
  * [98] X. Chen, H. Firouzjahi, E. Komatsu, M.H. Namjoo and M. Sasaki, _In-in and $\delta N$ calculations of the bispectrum from non-attractor single-field inflation_, _JCAP_ 12 (2013) 039 [1308.5341].
  * [99] Y.-F. Cai, X. Chen, M.H. Namjoo, M. Sasaki, D.-G. Wang and Z. Wang, _Revisiting non-Gaussianity from non-attractor inflation models_ , _JCAP_ 05 (2018) 012 [1712.09998].
  * [100] D.S. Salopek and J.R. Bond, _Nonlinear evolution of long wavelength metric fluctuations in inflationary models_ , _Phys. Rev. D_ 42 (1990) 3936.
  * [101] A.A. Starobinsky, _Dynamics of Phase Transition in the New Inflationary Universe Scenario and Generation of Perturbations_ , _Phys. Lett. B_ 117 (1982) 175.
  * [102] M. Sasaki and E.D. Stewart, _A General analytic formula for the spectral index of the density perturbations produced during inflation_ , _Prog. Theor. Phys._ 95 (1996) 71 [astro-ph/9507001].
  * [103] M. Sasaki and T. Tanaka, _Superhorizon scale dynamics of multiscalar inflation_ , _Prog. Theor. Phys._ 99 (1998) 763 [gr-qc/9801017].
  * [104] D.H. Lyth, K.A. Malik and M. Sasaki, _A General proof of the conservation of the curvature perturbation_ , _JCAP_ 05 (2005) 004 [astro-ph/0411220].
  * [105] D.H. Lyth and Y. Rodriguez, _The Inflationary prediction for primordial non-Gaussianity_ , _Phys. Rev. Lett._ 95 (2005) 121302 [astro-ph/0504045].
  * [106] A.D. Gow, H. Assadullahi, J.H.P. Jackson, K. Koyama, V. Vennin and D. Wands, _Non-perturbative non-Gaussianity and primordial black holes_ , _EPL_ 142 (2023) 49001 [2211.08348].
  * [107] A.A. Starobinsky, _Stochastic de Sitter (Inflationary) Stage in the Early Universe_ , _Lect. Notes Phys._ 246 (1986) 107.
  * [108] V. Vennin and A.A. Starobinsky, _Correlation Functions in Stochastic Inflation_ , _Eur. Phys. J. C_ 75 (2015) 413 [1506.04732].
  * [109] Y. Tada and V. Vennin, _Statistics of coarse-grained cosmological fields in stochastic inflation_ , _JCAP_ 02 (2022) 021 [2111.15280].
  * [110] H. Firouzjahi, A. Nassiri-Rad and M. Noorbala, _Stochastic Ultra Slow Roll Inflation_ , _JCAP_ 01 (2019) 040 [1811.02175].
  * [111] A. Achucarro, S. Cespedes, A.-C. Davis and G.A. Palma, _The hand-made tail: non-perturbative tails from multifield inflation_ , _JHEP_ 05 (2022) 052 [2112.14712].
  * [112] S. Pi and M. Sasaki, _Logarithmic Duality of the Curvature Perturbation_ , _Phys. Rev. Lett._ 131 (2023) 011002 [2211.13932].
  * [113] V. Atal, A. Sanglas and N. Triantafyllou, _NANOGrav signal as mergers of Stupendously Large Primordial Black Holes_ , _JCAP_ 06 (2021) 022 [2012.14721].
  * [114] D.H. Lyth and D. Wands, _Generating the curvature perturbation without an inflaton_ , _Phys. Lett. B_ 524 (2002) 5 [hep-ph/0110002].
  * [115] T. Moroi and T. Takahashi, _Effects of cosmological moduli fields on cosmic microwave background_ , _Phys. Lett. B_ 522 (2001) 215 [hep-ph/0110096].
  * [116] N. Bartolo and A.R. Liddle, _The Simplest curvaton model_ , _Phys. Rev. D_ 65 (2002) 121301 [astro-ph/0203076].
  * [117] S. Mollerach, _Isocurvature Baryon Perturbations and Inflation_ , _Phys. Rev. D_ 42 (1990) 313.
  * [118] M. Sasaki, J. Valiviita and D. Wands, _Non-Gaussianity of the primordial perturbation in the curvaton model_ , _Phys. Rev. D_ 74 (2006) 103003 [astro-ph/0607627].
  * [119] A.D. Gow, T. Miranda and S. Nurmi, _Primordial black holes from a curvaton scenario with strongly non-Gaussian perturbations_ , 2307.03078.
  * [120] D. Langlois, F. Vernizzi and D. Wands, _Non-linear isocurvature perturbations and non-Gaussianities_ , _JCAP_ 12 (2008) 004 [0809.4646].
  * [121] N. Bartolo, S. Matarrese and A. Riotto, _On nonGaussianity in the curvaton scenario_ , _Phys. Rev. D_ 69 (2004) 043503 [hep-ph/0309033].
  * [122] K. Enqvist and T. Takahashi, _Signatures of Non-Gaussianity in the Curvaton Model_ , _JCAP_ 09 (2008) 012 [0807.3069].
  * [123] K. Enqvist, S. Nurmi, O. Taanila and T. Takahashi, _Non-Gaussian Fingerprints of Self-Interacting Curvaton_ , _JCAP_ 04 (2010) 009 [0912.4657].
  * [124] K. Enqvist, S. Nurmi, G. Rigopoulos, O. Taanila and T. Takahashi, _The Subdominant Curvaton_ , _JCAP_ 11 (2009) 003 [0906.3126].
  * [125] O. Taanila, _Primordial Perturbations from a Self-interacting Curvaton_ , Ph.D. thesis, Helsinki U., 2010.
  * [126] K. Enqvist, _The self-interacting curvaton_ , _Prog. Theor. Phys. Suppl._ 190 (2011) 62 [1012.1711].
  * [127] C.T. Byrnes, K. Enqvist and T. Takahashi, _Scale-dependence of Non-Gaussianity in the Curvaton Model_ , _JCAP_ 09 (2010) 026 [1007.5148].
  * [128] J. Fonseca and D. Wands, _Non-Gaussianity and Gravitational Waves from Quadratic and Self-interacting Curvaton_ , _Phys. Rev. D_ 83 (2011) 064025 [1101.1254].
  * [129] C.T. Byrnes, K. Enqvist, S. Nurmi and T. Takahashi, _Strongly scale-dependent polyspectra from curvaton self-interactions_ , _JCAP_ 11 (2011) 011 [1108.2708].
  * [130] T. Kobayashi and T. Takahashi, _Runnings in the Curvaton_ , _JCAP_ 06 (2012) 004 [1203.3011].
  * [131] G. Domènech and M. Sasaki, _Conformal Frame Dependence of Inflation_ , _JCAP_ 04 (2015) 022 [1501.07699].
  * [132] K. Enqvist, S. Nurmi, G. Rigopoulos, O. Taanila and T. Takahashi, _The Subdominant Curvaton_ , _JCAP_ 11 (2009) 003 [0906.3126].
  * [133] EPTA collaboration, _The second data release from the European Pulsar Timing Array III. Search for gravitational wave signals_ , 2306.16214.
  * [134] D.J. Reardon et al., _Search for an Isotropic Gravitational-wave Background with the Parkes Pulsar Timing Array_ , _Astrophys. J. Lett._ 951 (2023) L6 [2306.16215].
  * [135] H. Xu et al., _Searching for the Nano-Hertz Stochastic Gravitational Wave Background with the Chinese Pulsar Timing Array Data Release I_ , _Res. Astron. Astrophys._ 23 (2023) 075024 [2306.16216].
  * [136] NANOGrav collaboration, _The NANOGrav 15 yr Data Set: Evidence for a Gravitational-wave Background_ , _Astrophys. J. Lett._ 951 (2023) L8 [2306.16213].
  * [137] NANOGrav collaboration, _The NANOGrav 15-year Data Set: Bayesian Limits on Gravitational Waves from Individual Supermassive Black Hole Binaries_ , 2306.16222.
  * [138] NANOGrav collaboration, _The NANOGrav 15-year Data Set: Constraints on Supermassive Black Hole Binaries from the Gravitational Wave Background_ , 2306.16220.
  * [139] Y. Gouttenoire, S. Trifinopoulos, G. Valogiannis and M. Vanvlasselaer, _Scrutinizing the Primordial Black Holes Interpretation of PTA Gravitational Waves and JWST Early Galaxies_ , 2307.01457.
  * [140] P.F. Depta, K. Schmidt-Hoberg and C. Tasillo, _Do pulsar timing arrays observe merging primordial black holes?_ , 2306.17836.
  * [141] T. Shinohara, T. Suyama and T. Takahashi, _Angular correlation as a novel probe of supermassive primordial black holes_ , _Phys. Rev. D_ 104 (2021) 023526 [2103.13692].
  * [142] T. Shinohara, W. He, Y. Matsuoka, T. Nagao, T. Suyama and T. Takahashi, _Supermassive primordial black holes: a view from clustering of quasars at $z\sim 6$_, 2304.08153.
  * [143] G. Franciolini, A. Iovino, Junior., V. Vaskonen and H. Veermae, _The recent gravitational wave observation by pulsar timing arrays and primordial black holes: the importance of non-gaussianities_ , 2306.17149.
  * [144] S. Wang, Z.-C. Zhao, J.-P. Li and Q.-H. Zhu, _Exploring the Implications of 2023 Pulsar Timing Array Datasets for Scalar-Induced Gravitational Waves and Primordial Black Holes_ , 2307.00572.
  * [145] J.M. Bardeen, _Gauge Invariant Cosmological Perturbations_ , _Phys. Rev. D_ 22 (1980) 1882.

##  Appendix A The $\delta N$ Formalism

###  A.1 Review

The $\delta N$ formalism [100, 101, 102, 103, 104, 105, 120] is a technique
for computing the non-linear curvature perturbation $\zeta$ on large scales by
identifying it with the perturbed logarithmic expansion from some initial
state to a final state of fixed energy density. In a homogenous background,
the number of $e$-folds elapsed between two moments of times $t_{1}$ and
$t_{2}$ is

$\bar{N}(t_{2},t_{1})=\int_{t_{1}}^{t_{2}}dt\,H\,.$ (A.1)

Meanwhile, the amount of expansion in a perturbed universe is [104]

$N(t_{2},t_{1};\vec{x})=\int_{t_{1}}^{t_{2}}dt\,(H+\dot{\psi})\,,$ (A.2)

where $\psi$ is the curvature perturbation appearing in the decomposition of
the spatial 3-metric $g_{ij}=a^{2}(t)e^{2\psi(t,\vec{x})}\delta_{ij}$.232323We
ignore tensor perturbations. Note that this expression holds on superhorizon
scales, where spatial gradients can be neglected. We define $\delta N\equiv
N(t_{2},t_{1};x)-\bar{N}(t_{2},t_{1})$ to be the difference between the
perturbed and unperturbed expansion, and see that it equates to the change in
the curvature perturbation from an initial hypersurface at $t_{1}$ to a final
hypersurface at $t_{2}$

$\delta N\equiv
N(t_{2},t_{1};\vec{x})-\bar{N}(t_{2},t_{1})=\psi(t_{2},\vec{x})-\psi(t_{1},\vec{x})\,.$
(A.3)

Equivalently, the change in $\psi$ going from one choice of slicing to another
is the difference between the actual number of $e$-folds $N$ and the
homogenous background value $\bar{N}$.

Since the curvature perturbation $\psi$ is a gauge-dependent quantity, whose
value depends on the choice of slicing, it is convenient to combine this with
the gauge-dependent density perturbation $\delta\rho$ to form the gauge-
invariant curvature perturbation [145]:

$\zeta(t,\vec{x})=\psi(t,\vec{x})+\frac{1}{3}\int_{\bar{\rho}(t)}^{\rho(t,\vec{x})}\frac{d\rho}{(1+w)\rho}\,,$
(A.4)

where $w$ is the equation of state of the cosmological fluid,
$\rho(t,\vec{x})$ is the inhomogenous local energy density, and
$\bar{\rho}(t)$ is the homogenous background energy density. It was
demonstrated in Ref. [104] that this quantity is conserved on superhorizon
scales.

We would like to equate this gauge-invariant, conserved curvature perturbation
$\zeta$ with the perturbed logarithmic expansion $\delta N$. By choosing the
initial hypersurface at $t_{1}$ to be spatially flat, such that
$\psi(t_{1},\vec{x})=0$, from Eq. (A.3) we can equate $\delta
N(t_{2},t_{1};\vec{x}=\psi(t_{2},\vec{x})$. By choosing the finial
hypersurface at $t_{2}$ to be uniform density, such that
$\rho(t_{2},\vec{x})=\bar{\rho}(t_{2})$, from Eq. (A.4) we have
$\zeta(t_{2},\vec{x})=\psi(t_{2},\vec{x})$. Combining these expressions gives
the $\delta N$ formula

$\zeta(t_{2},\vec{x})=\delta N(t_{2},t_{1};\vec{x})\,,$ (A.5)

relating the curvature perturbation with the perturbed expansion between a
spatially flat and uniform energy density hypersurface.

###  A.2 Numerical Implementation

It is also possible to implement the $\delta N$ formalism numerically.
Following inflaton decay, the exact system of equations describing the
evolution of the curvaton and radiation bath is

$\ddot{\chi}+\left(3H+\frac{1}{\tau}\right)\dot{\chi}+V_{\chi}^{\prime}=0~{}~{}~{},~{}~{}~{}\dot{\rho}_{R}+4H\rho_{R}=\frac{\dot{\chi}^{2}}{\tau}~{}~{}~{},~{}~{}~{}H^{2}=\frac{8\pi}{3M_{\rm
Pl}^{2}}(\rho_{R}+\rho_{\chi})~{}.~{}~{}~{}$ (A.6)

We set the initial conditions at the end of inflation by specifying
$\chi_{i}=\chi_{*}$, $\dot{\chi}_{i}=0$, and $H_{i}=H_{*}$, which in turn
determines $\rho_{R,i}=\frac{3M_{\rm
Pl}^{2}}{8\pi}H_{i}^{2}-V_{\chi}(\chi_{i})$. The system of equations should be
evolved until a final time $t_{f}$ satisfying $H_{f}\ll\tau^{-1}$, such that
the curvaton has decayed completely. The number of $e$-folds elapsed is then
computed as $N=\ln(a_{f}/a_{i})$. This procedure should then be repeated for
the perturbed field value $\chi_{*}+\delta\chi_{*}$, with the fluctuation
determined by the size of Hubble at horizon exit $\delta\chi_{*}=H_{*}/2\pi$.
Evolving until the same final hypersurface of fixed energy density, the
curvature perturbation is computed as:

$\zeta=N(\chi_{*}+\delta\chi_{*})-N(\chi_{*})\,.$ (A.7)

Repeating for many different $\delta\chi_{*}$ gives a functional relation
between $\zeta$ and the Gaussian $\delta\chi_{*}$. We leave an in-depth
numerical study to future investigation.

##  Appendix B Weak Quartic Interactions

The sample parameters of Fig. 7 can be obtained rather generically for a
variety of weakly interacting models. For concreteness, consider the quartic
potential:

$V_{\chi}=\frac{1}{2}m_{\chi}^{2}\chi^{2}+\frac{\lambda}{4}\chi^{4}\,.$ (B.1)

In the limit of small interactions, we can find an explicit functional form
for the amplitude at the onset of quadratic oscillations $\chi_{0}(\chi_{*})$
by treating the quartic term as a perturbation. We expand the solution as:

$\chi_{0}=\chi_{0}^{(0)}+\lambda\chi_{0}^{(1)}\,.$ (B.2)

At zeroth order in $\lambda$, the equation for $\chi_{0}^{(0)}$ reads:

$\ddot{\chi}_{0}^{(0)}+\frac{3}{2t}\dot{\chi}_{0}^{(0)}+m_{\chi}^{2}\chi_{0}^{(0)}=0\,,$
(B.3)

where we have assumed radiation domination such that $H=1/2t$. The general
solution is

$\chi_{0}^{(0)}=A_{0}\,w_{1/4}(m_{\chi}t)+B_{0}\,w_{-1/4}(m_{\chi}t)\,,$ (B.4)

where $w_{1/4}(m_{\chi}t)\equiv J_{1/4}(m_{\chi}t)/(m_{\chi}t)^{1/4}$ and
$w_{-1/4}(m_{\chi}t)\equiv J_{-1/4}(m_{\chi}t)/(m_{\chi}t)^{1/4}$, with $J$
Bessel functions of the first kind. For convenience we take evolution to start
at $t_{\rm int}=0$; demanding regularity here instructs us to set $B_{0}=0$.
The initial condition $\chi_{0}|_{t=0}=\chi_{*}$ fixes
$A_{0}=2^{1/4}\,\Gamma(5/4)\,\chi_{*}\simeq 1.078\,\chi_{*}$, such that:

$\chi_{0}^{(0)}=1.078\,\chi_{*}\,\frac{J_{1/4}(m_{\chi}t)}{(m_{\chi}t)^{1/4}}\,.$
(B.5)

At first order in $\lambda$, the correction satisfies

$\ddot{\chi}_{0}^{(1)}+\frac{3}{2t}\dot{\chi}_{0}^{(1)}+m_{\chi}^{2}\chi_{0}^{(1)}=-\left(\chi_{0}^{(0)}\right)^{3}\,.$
(B.6)

This second order non-linear differential equation can be solved with the
method of variation of parameters, wherein we take the Ansatz:

$\chi_{0}^{(1)}=A_{1}(t)\,w_{1/4}(m_{\chi}t)+B_{1}(t)\,w_{-1/4}(m_{\chi}t)\,.$
(B.7)

The corresponding system of equations for $\dot{A}_{1}$ and $\dot{B}_{1}$ is

$\begin{split}&\dot{A}_{1}\,w_{1/4}(m_{\chi}t)+\dot{B}_{1}\,w_{-1/4}(m_{\chi}t)=0\\\
&\dot{A}_{1}\,\dot{w}_{1/4}(m_{\chi}t)+\dot{B}_{1}\,\dot{w}_{-1/4}(m_{\chi}t)=-\left(\chi_{0}^{(0)}\right)^{3}\,,\end{split}$
(B.8)

with integral solutions

$\begin{split}&A_{1}=A_{0}^{3}\int
dt\,\frac{w_{-1/4}(m_{\chi}t)\,w_{1/4}(m_{\chi}t)^{3}}{w_{1/4}(m_{\chi}t)\,\dot{w}_{-1/4}(m_{\chi}t)-\dot{w}_{1/4}(m_{\chi}t)\,w_{-1/4}(m_{\chi}t)}\,,\\\
&B_{1}=-A_{0}^{3}\int
dt\,\frac{w_{1/4}(m_{\chi}t)^{4}}{w_{1/4}(m_{\chi}t)\,\dot{w}_{-1/4}(m_{\chi}t)-\dot{w}_{1/4}(m_{\chi}t)\,w_{-1/4}(m_{\chi}t)}\,.\end{split}$
(B.9)

We should integrate these from the onset of interactions $t_{\rm int}=0$ to
the onset of oscillations, corresponding to $m_{\chi}t\sim 1$, with the
result:

$A_{1}\simeq-0.689\frac{\chi_{*}^{3}}{m_{\chi}^{2}}\,,\,\,\,B_{1}\simeq
0.535\frac{\chi_{*}^{3}}{m_{\chi}^{2}}\,.$ (B.10)

The full solution is then

$\chi_{0}\simeq
1.078\,\chi_{*}\frac{J_{1/4}(m_{\chi}t)}{(m_{\chi}t)^{1/4}}+\lambda\left(-0.689\,\frac{J_{1/4}(m_{\chi}t)}{(m_{\chi}t)^{1/4}}+0.535\,\frac{J_{-1/4}(m_{\chi}t)}{(m_{\chi}t)^{1/4}}\right)\frac{\chi_{*}^{3}}{m_{\chi}^{2}}\,.$
(B.11)

Finally evaluating this at $m_{\chi}t\sim 1$ gives

$\chi_{0}\simeq
0.81\,\chi_{*}-0.16\,\lambda\frac{\chi_{*}^{3}}{m_{\chi}^{2}}\,.$ (B.12)

This expression can then be substituted into Eq. (6.14) to obtain
$\zeta_{\chi}$. We see for example that the benchmark point of Fig. 7 can be
realized for $(\chi_{*},\,m_{\chi},\,\lambda)=(3.6,\,1,\,0.5)$.