#  Computing the abundance of primordial black holes

Sam Young<EMAIL_ADDRESS>Instituut-Lorentz for Theoretical
Physics, Leiden University,
Niels Bohrweg 2, 2333 CA Leiden, The Netherlands

###### Abstract

An accurate calculation of their abundance is crucial for numerous aspects of
cosmology related to primordial black holes (PBHs). For example, placing
constraints on the primordial power spectrum from constraints on the abundance
of PBHs (or vice-versa), calculating the mass function observable today, or
predicting the merger rate of (primordial) black holes observable by
gravitational wave observatories such as LIGO, Virgo and KAGRA. In this
chapter, we will discuss the different methods used for the calculation of the
abundance of PBHs forming from large-amplitude cosmological perturbations,
assuming only a minimal understanding of modern cosmology. Different
parameters to describe cosmological perturbations will be considered
(including different choices for the window function), and it will be argued
that the compaction is typically the most appropriate choice. Different
methodologies for calculating the abundance and mass function are explained,
including _Press-Schechter_ -type and peaks theory approaches.

###### Contents

  1. I Introduction
  2. II Describing cosmological perturbations
    1. II.1 Smoothing functions
    2. II.2 The compaction function
  3. III Perturbation profiles
  4. IV Methods for the calculation
    1. IV.1 The _Press-Schechter_ -type approach
      1. IV.1.1 The critical scaling relationship
      2. IV.1.2 The total primordial black hole abundance
      3. IV.1.3 The mass function
    2. IV.2 The theory of peaks
    3. IV.3 Extensions to the theory of peaks
  5. V The effect of non-linearities
  6. VI Summary

##  I Introduction

Primordial black holes (PBHs) have a wealth of potential implications for
cosmology - from placing unique constraints on the primordial power spectrum,
solving numerous unexplained phenomena (such as the early formation of
supermassive black holes (SMBHs)), and the observed presence of dark matter in
the Universe. Drawing robust conclusions in such scenarios typically relies on
a calculation of quantities such as the abundance of PBHs from given initial
conditions, the masses at which they form, and their initial clustering.
Making this calculation as accurate and reliable as possible is therefore
crucial.

In this chapter, we will discuss the calculation of the abundance of PBHs from
the collapse of large density fluctuations during radiation domination, but
let us first turn our attention to the cosmological setting in which such PBHs
form. During (exponential) inflation, the physical Hubble horizon is constant,
and quantum fluctuations become classicial density perturbations as they exit
the horizon as the universe expands. Once inflation ends and the universe
becomes radiation dominated, the Hubble horizon starts to grow, and modes
which had previously exited the horizon re-enter. Figure 1 shows a schematic
plot of how perturbations exit and re-enter the horizon. Larger-scale
perturbations exit the horizon earlier, and then re-enter the horizon later.
Once a perturbation re-enters the horizon, if the amplitude of the
perturbation is large enough, gravity can overcome pressure forces, causing it
to quickly collapse and form a PBH Zel’dovich and Novikov (1967); Hawking
(1971); Carr (1975), with a mass of approximately the horizon mass, or else be
quickly damped out if the perturbation is not dense enough. This means that
larger-scale perturbations correspond to earlier times during inflation, and
more massive PBHs which form at later times.

Figure 1: _Horizon exit and re-entry._ The black line shows the Hubble horizon
as a function of time, remaining constant during (quasi-)exponential
inflation, before growing during the radiation dominated epoch following
reheating. The dashed red and dotted blue lines show how perturbations of a
small and large scale, respectively, first exit the horizon during inflation
before re-entering at a later time. PBHs can form upon horizon re-entry if the
amplitude of a perturbation is large enough.

Once an abundance of PBHs has formed, they behave on large scales like a
pressure-less fluid, and the PBH density grows inversely proportional to the
cube of the scale factor,

$\rho_{\mathrm{PBH}}\propto a^{-3},$ (1)

i.e. the PBH density on large scales evolves like matter, whilst the density
of the radiation fluid grows inversely to the fourth power of the scale
factor,

$\rho_{\mathrm{r}}\propto a^{-4}.$ (2)

Although not relevant for our discussion here, the final quantity which makes
up a significant fraction of the universe today is dark energy, which we will
assume is due to a cosmological constant $\Lambda$. The density of the
cosmological constant remains constant over time,

$\rho_{\Lambda}\propto\mathrm{const.}$ (3)

The density parameter $\Omega_{i}$ describes the density of a particular fluid
$\rho_{i}$ relative to the total density $\rho_{T}$,

$\Omega_{i}=\frac{\rho_{i}}{\rho_{T}}.$ (4)

Figure 2 shows how the density parameter changes over time for matter,
radiation and the cosmological constant. PBHs form very early during radiation
domination, at a time when the (relative) density of matter is very low.
Because the abundance of PBHs cannot be larger than the abundance of matter,
with the exception of very light PBHs which evaporate very quickly, the
abundance of PBHs in the universe at the time of formation is typically
assumed to be very small: PBHs are extremely rare at the time of formation.

Figure 2: _The evolution of cosmological densities._ The figure show the
relative energy density of different species within the universe. Initialy,
following reheating, the energy density of the universe was dominated by
radiation. As the universe expands, the radiation density decreases faster
than the matter density, and the universe enter a matter dominated epoch.
Finally, the universe becomes dominated by the cosmological constant, whose
energy density remains constant. Because PBHs act like matter, and form very
early during radiation domination, the fraction of the energy of the universe
contained in PBHs is typically expected to be very small (expect in models
where such PBHs evaporate extremely rapidly).

For PBHs to form, significantly larger perturbations are needed than are
observed on large scales in e.g. the cosmic microwave background radiation
(CMB). Perturbations are required to be approximately $50\%$ more dense than
the background Musco  _et al._ (2005, 2009); Nakama  _et al._ (2014a, b). A
consequence of this is that the primordial power spectrum must be orders of
magnitude larger on small scales than the value of $\sim 2\times 10^{-9}$
observed in the CMB Akrami  _et al._ (2020), with the typical values quoted in
the literature around $\mathcal{O}(10^{-2})$ (see e.g. Carr  _et al._ (2021)).
However, there can be significant deviation in this value depending on the
method used (e.g.Gow  _et al._ (2021)).

In this chapter, we will first discuss how cosmological perturbations can be
described and how they evolve over time. We will then discuss the typical
profile shapes of perturbations, and how this can affect the calculation of
the PBH abundance. We will then explain the different methods which can be
used to calculate the PBH abundance, beginning with a simple _Press-Schechter_
-type approach, before moving on to the theory of peaks and its extensions.
Finally, we will consider the effect of non-linearities in the density, which
become important as large perturbations are required for PBH formation.

##  II Describing cosmological perturbations

Many models of inflation can predict a large number of PBHs to form in the
early universe, including the running mass model Drees and Erfani (2011),
axion inflation Bugaev and Klimai (2014); Özsoy  _et al._ (2018), or a
waterfall transition during hybrid inflation Garcia-Bellido  _et al._ (1996);
Lyth (2012); Bugaev and Klimai (2012), and many others. The curvature
perturbation $\zeta$ is typically used to study cosmological perturbations
generated with the different models and to predict their observable
consequences (NB. Often $\mathcal{R}$ is used instead to denote the curvature
perturbation, whilst $\zeta$ is used to describe the curvature perturbation
only on a uniform density slicing), and is also what is typically constrained
by large-scale observations, such as the _Planck_ satellite observations of
the CMB Akrami  _et al._ (2020). The curvature perturbation $\zeta$ appears as
a perturbative quantity in the Friedmann-Lemaitre-Robertson-Walker (FLRW)
metric in the comoving uniform-density gauge as

$\mathrm{d}s^{2}=-\mathrm{d}t^{2}+a^{2}(t)\exp(2\zeta)\mathrm{d}\mathbf{X}^{2},$
(5)

where $t$ is time, $a(t)$ is the scale factor and $\mathbf{X}$ represents the
three comoving spatial coordinates. To use constraints on PBHs to constrain
models of inflation (or to make predictions of PBH abundances from a given
model), we therefore need to relate the primordial curvature perturbation
power spectrum $\mathcal{P}_{\zeta}(k)$ to the abundance of PBHs of different
masses.

However, the curvature perturbation itself is not a suitable parameter to
determine PBH formation. Whilst it is possible to state a threshold value for
$\zeta$ above which a fluctuation will collapse to form a PBH Shibata and
Sasaki (1999), this assumes that the curvature perturbation goes to zero
outside the scale of the perturbation. This is not generally the case: small-
scale perturbations will typically be superposed on top of large-scale
perturbations. Naively using the curvature perturbation as the formation
threshold would therefore suggest that larger-scale perturbations have a
significant impact on PBH formation. However, precisely because such
perturbations are well outside the horizon at the time of PBH formation, the
separate universe approach suggests that it should not affect the local
evolution of the universe Wands  _et al._ (2000). This can be seen by
separating the curvature perturbation into a small horizon-scale perturbation
$\zeta_{s}$ and a large super-horizon perturbation $\zeta_{l}$, and then
treating the large-scale super-horizon mode as a constant:

$\zeta=\zeta_{s}+\zeta_{l}.$ (6)

We then consider the term in the FLRW metric where the curvature perturbation
appears,

$a(t)\exp(\zeta)\rightarrow
a(t)\exp(\zeta_{s}+\zeta_{l})=\left[a(t)\exp(\zeta_{l})\right]\exp(\zeta_{s}).$
(7)

If we then consider that $\zeta_{l}$ is (very close to) constant on the
horizon scale, we can treat that term as simply a rescaling of the scale
factor, corresponding to a time shift $t^{\prime}$ in the background evolution
of the universe,

$a(t)\rightarrow a(t)\exp(\zeta_{l})=a(t+t^{\prime}).$ (8)

We conclude that, super-horizon modes in $\zeta$ should not affect whether a
PBH forms or not, but can have an effect at the time at which a PBH would
form.

Seemingly the most obvious parameter to use is then the density contrast
$\delta$, defined as

$\delta(t,\mathbf{x})\equiv\frac{\rho-\bar{\rho}}{\bar{\rho}}$ (9)

which, at first order111We will consider the effect of non-linearities arising
from higher order terms in section V., is related to $\zeta$ as

$\delta(t,\mathbf{x})=-\frac{2(1+\omega)}{5+3\omega}\left(\frac{1}{aH}\right)^{2}\nabla^{2}\zeta(\mathbf{x}),$
(10)

where $\mathbf{x}$ is the spatial coordinate, $\rho$ is the density,
$\bar{\rho}$ is its background value, $\omega$ is the equation of state (we
will consider that $\omega=1/3$ during radiation domination for the rest of
this chapter), $H(t)=\dot{a}(t)/a(t)$ is the Hubble parameter, and
$\nabla^{2}=\frac{\partial}{\partial x}+\frac{\partial}{\partial
y}+\frac{\partial}{\partial z}$ is the Laplacian. The comoving Hubble horizon
is given by $R_{\mathrm{H}}=(aH)^{-1}$. In Fourier space, this relation is
given by

$\delta(t,k)=\frac{2(1+\omega)}{5+3\omega}\left(\frac{k}{aH}\right)^{2}\zeta(t,k).$
(11)

On super-horizon scales, $\zeta$ is constant over time,
$\zeta(t,k)\rightarrow\zeta(k)$, and so the term
$\left(\frac{k}{aH}\right)^{2}$ contains the time-dependance of the density
contrast: perturbations becoming larger over time as the horizon grows, and
larger-scale modes are suppressed by a factor $k^{2}$.

On sub-horizon scales, it is necessary to account for the evolution of
perturbations, which, at linear order, is easily accounted for in Fourier
space by including a transfer function $T(R_{\mathrm{H}},k)$ which includes
the time-dependence,

$\delta(R_{\mathrm{H}},k)=\frac{2(1+\omega)}{5+3\omega}\left(\frac{k}{aH}\right)^{2}\zeta(k)T(R_{\mathrm{H}},k),$
(12)

with the transfer function given by Josan  _et al._ (2009),

$T(R_{\mathrm{H}},k)=3\frac{\sin\left(\frac{kR_{\mathrm{H}}}{\sqrt{3}}\right)-\frac{kR_{\mathrm{H}}}{\sqrt{3}}\cos\left(\frac{kR_{\mathrm{H}}}{\sqrt{3}}\right)}{\left(\frac{kR_{\mathrm{H}}}{\sqrt{3}}\right)^{3}},$
(13)

where the Hubble horizon scale $R_{\mathrm{H}}(t)$ is a function of time.
However, it will be more convenient for our purposes to simply use
$R_{\mathrm{H}}$ itself as the time parameter. On super-horizon scales, $k\ll
R_{\mathrm{H}}$, and the transfer function approach unity, recovering equation
(11), and is highly oscillatory on sub-horizon scales.

It is common to use the density contrast as the formation criterion: for any
perturbations where the density peaks above some threshold $\delta_{c}$, it is
considered that a PBH will form. However, simulations reveal that the
formation (or not) of a PBH is best determined by the average density within
the horizon, rather than the peak value (likely at the center of the
perturbation) Musco (2019), which brings us to the discussion of smoothing
functions (also referred to as window functions, or simply windows).

###  II.1 Smoothing functions

A smoothing function can be considered as an averaging function - it returns
the average value of a quantity within a specified volume, centred on a
certain coordinate. In order to calculate the smoothed density contrast
$\delta_{R}$, a smoothing function $W(R_{s},\mathbf{x})$ is convolved with the
density contrast:

$\delta_{R}(t,\mathbf{x},R_{s})=\int\mathrm{d}^{3}\mathbf{y}\delta(t,\mathbf{x}-\mathbf{y})W(R_{s},\mathbf{y}),$
(14)

where $R_{s}$ is the smoothing scale, and the integration is performed over
all space (or at least, the full extent of the smoothing function). In Fourier
space, the density is expressed simply as a multiplication with the smoothing
function:

$\delta_{R}(t,k)=\delta(t,k)\tilde{W}(k,R),$ (15)

where the tilde denotes a Fourier transform.

A smoothing function will typically define a volume of space to be integrated
over (potentially giving different weight to different regions), divided by
the volume. Conceptually, the simplest smoothing function is the top-hat
window function $T_{\mathrm{TH}}$,

$W_{\mathrm{TH}}(\mathbf{x},R_{s})=\frac{3}{4\pi R_{s}^{3}}\Theta_{H}(R-x),$
(16)

where $\Theta_{H}(x)$ is the Heaviside theta function, and
$x=\left|\mathbf{x}\right|$. We can insert this into equation (14) to obtain

$\delta_{R}(t,\mathbf{x},R_{s})=\frac{3}{4\pi
R_{s}^{3}}\int\mathrm{d}^{3}\mathbf{y}\delta(t,\mathbf{x}-\mathbf{y})\Theta_{H}(R-y).$
(17)

We further simplify by choosing to center our coordinates on $\mathbf{x}$ and
assuming spherical symmetry around this point,

$\delta_{R}(t,0,R_{s})=\frac{3}{4\pi
R_{s}^{3}}\int\limits_{0}^{R_{s}}\mathrm{d}\mathbf{r}4\pi r^{2}\delta(t,r),$
(18)

where we have used $\mathrm{d}^{3}y=4\pi r^{2}\mathrm{d}r$, which is valid
under spherical symmetry. It is clear, here, that the smoothed value
$\delta_{R}$ is simply $\delta$ integrated over a sphere of radius $R_{s}$,
divided by the volume of the sphere i.e. the volume-averaged value.

In addition to the top-hat smoothing function, another common choice is the
Gaussian smoothing function,

$W_{G}(\mathbf{x},R_{s})=\frac{1}{\left(\pi
R_{s}^{2}\right)^{3/2}}\exp\left(\frac{-x^{2}}{R_{s}^{2}}\right).$ (19)

We note that, often, there is an additional factor of $1/2$ included in the
exponent of this function which has been neglected here. This is a free
choice, but is made here such that the characteristic scales of the Gaussian
and top-hat window functions coincide (i.e. the integrand $4\pi
r^{2}W(r,R_{s})$ peaks at $R_{s}$ in both cases).

The Fourier transforms of these smoothing functions are given by

$\tilde{W}_{\mathrm{TH}}(k,R_{s})=3\frac{\sin\left(kR_{\mathrm{H}}\right)-\frac{kR_{\mathrm{H}}}{\sqrt{3}}\cos\left(\frac{kR_{\mathrm{H}}}{\sqrt{3}}\right)}{\left(\frac{kR_{\mathrm{H}}}{\sqrt{3}}\right)^{3}},$
(20)
$\tilde{W}_{\mathrm{G}}(k,R_{s})=\exp\left(\frac{-k^{2}R_{s}^{2}}{4}\right),$
(21)

Figure 3 shows the real- and Fourier-space transforms of the top-hat and
Gaussian smoothing functions. Which smoothing function to use is the subject
of ongoing debate, but _a priori_ there is currently no reason to choose one
smoothing over another - although both presented here have their own
advantages and disadvantages which will be discussed later.

Figure 3: _Smoothing functions._ The shape of top-hat and Gaussian smoothing
functions are shown in real-space (left) and Fourier-space (right).

The use of a smoothing function has several benefits compared to using the
unsmoothed density contrast. As previously mentioned, the threshold for
collapse is more accurately stated in terms of the smoothed density contrast.
Whilst the exact value depends on the shape of a the specific perturbation
being considered, the critical value varies in the range
$0.41\lesssim\delta_{R,c}\lesssim 0.67$ for extreme shapes of perturbations.
However, the critical value for the peak value of the (unsmoothed) density
contrast can, in principle vary between positive and negative infinity (e.g.
for perturbations with a very sharp peak or trough at the centre) Musco
(2019). In addition, smoothing the distribution on different scales means that
the abundance of PBHs forming on different scales (and, therefore, of
different masses) can be calculated Young (2019).

###  II.2 The compaction function

Very closely related to the smoothed density contrast, the compaction function
$C$ has, in recent years, become the standard parameter used to describe
perturbations, and is defined as

$C(r,t,\mathbf{x})=\frac{2\delta M(r,t,\mathbf{x})}{r},$ (22)

where $\delta M(r,\mathbf{x})=M(r,t,\mathbf{x})-\bar{M}(r,t)$ is the mass
excess within a sphere of areal radius $r$ centred on spatial coordinate
$\mathbf{x}$; $M(r,t,\mathrm{x})$ is the Misner-Sharp mass and the bar denotes
the background value in a region of unperturbed space. As we will see soon,
the compaction is constant on superhorizon scales, and the time dependance is
typically dropped. The mass excess $\delta M$ can be expressed as an integral
over the density

$\delta
M(r,t,\mathbf{x})=\int\mathrm{d}^{3}\mathbf{y}\left[\rho(t,\mathbf{x}-\mathbf{y})-\bar{\rho}(t)\right]\Theta_{H}(r-y).$
(23)

Using natural units, $c=G=1$, and assuming a flat universe, the Friedmann
equation gives an expression for the background density,

$\bar{\rho}(t)=\frac{3H^{2}(t)}{8\pi}.$ (24)

Using the definition for the density contrast, equation (9), allows us to
rewrite equation (23) as

$\delta
M(r,t,\mathbf{x})=\frac{3H^{2}}{4\pi}\int\mathrm{d}^{3}\mathbf{y}\delta(t,\mathbf{x}-\mathbf{y})\Theta_{H}(r-y).$
(25)

The compaction function can then be expressed as,

$C(r,t,\mathbf{x})=\frac{r^{2}}{r_{\mathrm{H}}^{2}}\int\mathrm{d}^{3}\mathbf{y}\delta(t,\mathbf{x}-\mathbf{y})\frac{3}{4\pi
r^{3}}\Theta_{H}(r-y),$ (26)

where $r_{\mathrm{H}}=H^{-1}$ is the physical Hubble radius. The term inside
the integral is the top-hat smoothing function, equation (16), so finally we
can write

$C_{\mathrm{TH}}(r,t,\mathbf{x})=\frac{r^{2}}{R_{\mathrm{H}}^{2}}\int\mathrm{d}^{3}\mathbf{y}\delta(t,\mathbf{x}-\mathbf{y})W_{\mathrm{TH}}(r,\mathbf{y}),$
(27)

which is exactly the same as the smoothed density contrast, equation (14),
where $r$ has taken the role of the smoothing scale, and multiplied by
$r^{2}/r_{\mathrm{H}}^{2}$. The compaction function can also be defined using
a Gaussian smoothing:

$C_{\mathrm{G}}(r,t,\mathbf{x})=\frac{r^{2}}{R_{\mathrm{H}}^{2}}\int\mathrm{d}^{3}\mathbf{y}\delta(t,\mathbf{x}-\mathbf{y})W_{\mathrm{G}}(r,\mathbf{y}).$
(28)

The compaction function is therefore extremely closely related to the smoothed
density contrast, and indeed, the two are often used interchangeably in the
literature (see e.g. reference Young (2019)). This is especially true as the
smoothed density contrast is often evaluated at horizon entry, such that
$r/r_{\mathrm{H}}=1$. However, the compaction is typically preferred as it is
conserved on super-horizon scales. This can be seen by inserting equation (9),
and replacing the physical radius $r=aR$ with the comoving radius $R$:

$C(R,\mathbf{x})=-R^{2}\frac{2(1+\omega)}{5+3\omega}\int\mathrm{d}^{3}\mathbf{y}\nabla^{2}\zeta(\mathbf{x}-\mathbf{y})W(R,\mathbf{y}),$
(29)

where the time-dependance has been dropped, because $\zeta$ is conserved on
super-horizon scales.

It is typically simpler to leave this integral as it appears, as the
expressions in Fourier space are quite simple (e.g. equation (15)). However,
it will be useful later if we perform the integration analytically for the
case of the top-hat smoothing, which can be completed under the assumption of
spherical symmetry. This assumption is typically justified because, as
discussed in section I, PBHs are extremely rare at the time of formation -
which therefore means that the perturbations from which they form are very
rare, due to the large amplitude required. The theory of peaks tells us that,
the larger a perturbation is, the more spherical symmetry it will have Bardeen
_et al._ (1986)222However, we note that it has been pointed out that, whilst
the perturbations to the density (or compaction) are large and rare, this does
not necessarily correspond to large, rare perturbations in $\zeta$ Young
(2022).. Under the assumption of spherical symmetry, the Laplacian
$\nabla^{2}\zeta$ becomes

$\nabla^{2}\zeta=\zeta^{\prime\prime}(r)+\frac{2}{r}\zeta^{\prime}(r),$ (30)

where $r$ is here the radial coordinate, and the prime denotes a derivative
with respect to $r$. The compaction can then be expressed as

$C(R)=-\frac{2(1+\omega)}{5+3\omega}\frac{3}{R}\int\limits_{0}^{R}\mathrm{d}rr^{2}\left[\zeta^{\prime\prime}(r)+\frac{2}{r}\zeta^{\prime}(r)\right]=-\frac{6(1+\omega)}{5+3\omega}R\zeta^{\prime}(R).$
(31)

To avoid confusion, we stress here that, the $\zeta^{\prime}$ term has arisen
as the surface term of an integral over a spherical volume - therefore
corresponding to $\zeta^{\prime}$ smoothed over a spherical shell.

##  III Perturbation profiles

The single most important parameter when considering the abundance of PBHs
which will form from the collapse of cosmological perturbations is the
threshold amplitude of the perturbations required for collapse, $C_{c}$.
Carr’s initial order-of-magnitude estimate, calculated using a Jeans’ Length
argument, was that the threshold for collapse is equal to the equation of
state, $C_{c}\sim\omega$ Carr and Hawking (1974); Carr (1975), which is equal
to $1/3$ during radiation domination. Since then, many more complex methods
have been developed, utilising either analytic methods (e.g. Harada and
Jhingan (2016); Escrivà  _et al._ (2021)), or numerical simulations (e.g.
Musco  _et al._ (2005); Polnarev and Musco (2007); Musco  _et al._ (2009);
Nakama  _et al._ (2014c); Escrivà (2020); Escrivà  _et al._ (2020)). The
precise value depends on the initial profile shape of the perturbation which
forms a PBH, which is a subject of immense scrutiny. In this section, we will
briefly turn our attention towards the profile shape one might therefore
expect for cosmological perturbations.

We will first consider how a single perturbation may be described using the
various parameters discussed in section II, and we will consider only the
initial perturbation on super-horizon scales, before any complicated evolution
close to horizon entry. We therefore set the transfer function equal to unity,
$T(t,k)=1$. We will take a spherically symmetric perturbation in $\zeta$ in an
otherwise empty universe, given by a Gaussian profile:

$\zeta(r)=A\exp\left(-\frac{r^{2}}{r_{p}^{2}}\right),$ (32)

where $A$ describes the amplitude of the perturbation, $r_{p}$ determines the
scale of the perturbation, and $r$ is the comoving radial coordinate. The same
perturbation can be described in terms of the density contrast $\delta$ and
results in the ubiquitous Mexican-hat profile,

$\delta(r)=-\frac{4(1+\omega)}{5+3\omega}\left(\frac{1}{aH}\right)^{2}A\frac{\exp\left(-\frac{r^{2}}{r_{p}^{2}}\right)\left(2r^{2}-3r_{p}^{2}\right)}{r_{p}^{4}},$
(33)

as well as the compaction $C$ (for a top-hat smoothing),

$C(r)=-\frac{12(1+\omega)}{5+3\omega}A\frac{r^{2}}{r_{p}^{2}}\exp\left(-\frac{r^{2}}{r_{p}^{2}}\right).$
(34)

Figure 4 shows the profile shapes in each of these parameters, where we have
made the following choices,
$\left[A,r_{p},(aH)^{-1}\right]=\left[0.5,1,1\right]$. $A$ is chosen to
correspond approximately to the amplitude for PBH formation, whilst the choice
of $r_{p}$ is arbitrary, and $(aH)^{-1}$ is chosen so that the amplitude of
the density is similar to the others (technically, in the super-horizon limit,
we should have chosen $(aH)^{-1}\ll r_{p}$). For $\zeta$ and $\delta$, the
profile simply represents the value of that parameter at the coordinate $r$,
and this is also often considered the case for $C$. However, it is more
accurate to consider the shown profile corresponds to the value of the
compaction at the centre, when a smoothing scale $R_{s}=r$ has been used.

Figure 4: _Parameterising a perturbation._ The same perturbation is shown,
described with several different parameters: the curvature perturbation
$\zeta$, the density contrast $\delta$ and the compaction $C$. For $\zeta$ and
$\delta$ the plot shows the value at the radial coordinate $r$. For the
compaction, the plot shows the amplitude of the compaction at the centre of
the perturbation, when smoothed with a smoothing scale $r$.

For any profile shape, the compaction is generically zero at $r\rightarrow 0$
and $r\rightarrow\infty$, and peaks at some value of $r$. The value at which
this occurs is used to define the characteristic scale $r_{m}$ of the
perturbation: $C^{\prime}(r_{m})\equiv 0$, and can be found by solving

$\zeta^{\prime}(r_{m})+r_{m}\zeta(r_{m})=0.$ (35)

The characteristic scale $r_{m}$ will be important later as it corresponds to
the scale (and mass) at which a PBH forms. In this case, with the profile
given in equation (32), we find that $r_{m}$ coincides with our choice of
parameter for the scale, $r_{m}=r_{p}$.

We will now turn our attention to the consideration of what profile shapes we
might expect to see in the density in different scenarios - as this will
affect the amplitude a perturbation needs to have in order to form a PBH. The
expected profile shape can be calculated by considering the amplitude of
different Fourier modes in the universe Bardeen  _et al._ (1986), which is
described by the power spectrum. The power spectrum $P_{\zeta}(k)$ is the
Fourier transform of the 2-point correlation function,

$\langle\zeta(\mathbf{k}_{1})\zeta(\mathbf{k}_{2})\rangle=(2\pi)^{3}\delta_{D}^{(3)}(\mathbf{k}_{1}+\mathbf{k}_{2})P_{\zeta}(k),$
(36)

where $\delta_{D}^{(3)}$ is the (3D) Dirac-delta function. We will consider
the dimensionless power spectrum $\mathcal{P}_{\zeta}(k)$, given by

$\mathcal{P}_{\zeta}(k)=\frac{k^{3}}{2\pi^{2}}P_{\zeta}(k).$ (37)

The effect of the power spectrum on perturbation profiles is considered in
detail in the context of PBH formation by references Germani and Musco (2019);
Young (2019). In this case, we will consider density profiles, and so consider
the density power spectrum, $\mathcal{P}_{\delta}(k)$, which is related to the
curvature perturbation power spectrum, $\mathcal{P}_{\zeta}(k)$, as

$\mathcal{P}_{\delta}(R_{\mathrm{H}},k)=\frac{4(1+\omega)^{2}}{(5+3\omega)^{2}}\left(kR_{\mathrm{H}}\right)^{4}\mathcal{P}_{\zeta}(k)T^{2}\left(R_{\mathrm{H}},k\right),$
(38)

where, again, we use the comoving Hubble horizon $R_{\mathrm{H}}$ as the
parameter for time, and will here consider the initial perturbations long
before horizon entry, such that $k\ll R_{\mathrm{H}}$, and
$T(R_{\mathrm{H}},k)\approx 1$.

Assuming spherical symmetry, the profile shape close to the centre of a
perturbation can be stated in terms of the radial coordinate $r$,

$\delta(R_{\mathrm{H}},r)=\delta(R_{\mathrm{H}},0)\frac{\xi(R_{\mathrm{H}},r)}{\xi(R_{\mathrm{H}},0)},$
(39)

where $\delta(R_{\mathrm{H}},0)$ is the amplitude at the centre of the
perturbation.The function $\xi(R_{\mathrm{H}},r)$, which describes the average
profile shape, is given by

$\xi(R_{\mathrm{H}},t)=\frac{1}{(2\pi)^{3}}\int\frac{\mathrm{d}k}{k}\frac{\sin(kr)}{kr}\mathcal{P}_{\delta}(R_{\mathrm{H}},k),$
(40)

where equation (38) relates the density and curvature power spectra. We will
consider several forms for the curvature perturbation power spectrum, which
cover a wide range of models in the literature:

  * •

the (Dirac-)delta power spectrum

$\mathcal{P}_{\zeta}=\mathcal{A}\delta_{\mathrm{D}}(k-k_{*}),$ (41)

where $\delta_{\mathrm{D}}$ is the Dirac-delta function. It is notable that
this doesn’t correspond to a realistic or physical power spectrum, it is often
used as a reasonable approximation to sufficiently narrow spectra.

  * •

the narrow-peak power spectrum

$\mathcal{P}_{\zeta}=\mathcal{A}\exp\left(-\frac{(k-k_{*})^{2}}{2\Delta^{2}}\right),$
(42)

which contains a Gaussian-shaped peak. It is also possible to utilise a
lognormal peak instead.

  * •

The power-law power spectrum

$\mathcal{P}_{\zeta}=\mathcal{A}\left(\frac{k}{k_{*}}\right)^{n_{s}-1}\Theta_{H}(k-k_{min})\Theta_{H}(k_{*}-k),$
(43)

where the Heaviside functions provide both a UV and IR cut-off. We will here
consider the scale-invariant power spectrum, $n_{s}=1$.

Whilst these describe very different power spectra, one thing that they all
share in common is a cut-off at some scale in both the UV and IR. Without the
UV cut-offs, the integral in equation (40) diverges - and, as we will see
later, the variance of the distribution also diverges. The IR cut-off is not
so significant, as the $k^{4}$ term when converting to the density power
spectrum, as in equation (38), naturally provides a sharp IR cut-off.

Figure 5 shows the resulting density power spectra, using equation (38), and
average density profiles, calculated using equation (39). We can see that,
whilst the curvature perturbation power spectra are very different, the
density power spectrum are well characterised by a sharp peak. As a result,
the average density profile in each case is also very similar (although note
that the characteristic scale $r_{m}$ of the average profiles does vary
slightly, and it has been normalised to 1 in the figure), and matches closely
to the Mexican-hat profile already considered, equation (33). Because of this,
the threshold amplitude for collapse also shows little variation, and we can
safely use the threshold for the Mexican-hat profile shape, $C_{c}\approx
0.5$. There is, however, some dependance on the shape of the power spectrum,
which is discussed in reference Germani and Musco (2019).

Figure 5: _Power spectra and average profile shape._ The left plot shows
several different possible choices for the density power spectrum, normalised
to have an amplitude of unity at the peak value. The right plot shows that the
average profile associated with such power spectrum, with the Mexican hat
profile shape also included with the green dotted line.

We note that this critical value corresponds to a top-hat smoothing, and, as
we will see later, the choice of smoothing function can have a significant
impact on the calculated PBH abundance. If the choice is made to use a
Gaussian smoothing function for the calculation, then the critical value for
the Gaussian-smoothed compaction should be used: $C_{G,c}\approx 0.18$. When
properly accounted for, this introduces an uncertainty in the amplitude of the
power spectrum required to form a given abundance of PBHs of order $10\%$
Young (2019).

##  IV Methods for the calculation

We will first discuss some generalities with the calculation, before
discussing the simplest method, a _Press-Schechter_ -type approach, and
finally the theory of peaks and its extensions. Throughout this chapter, we
will assume a Gaussian distribution for the curvature perturbation $\zeta$. At
first-order, the density $\delta$, compaction $C$ and curvature perturbation
$\zeta$ are linearly related, implying both $\delta$ and $C$ also follow a
Gaussian distribution. Since we will use $C$ for the calculation, we will
consider the probability density function (PDF) of $C$ given by,

$P(C)=\frac{1}{\sqrt{2\pi\sigma_{0}^{2}}}\exp\left(-\frac{C^{2}}{2\sigma_{0}^{2}}\right),$
(44)

where $\sigma_{0}^{2}$ is its variance. The variance can be calculated by
taking the zeroth-order (n=0) moment of the power spectrum,

$\sigma_{n}^{2}(R_{s},R_{\mathrm{H}})=\int\limits_{0}^{\infty}\frac{\mathrm{d}k}{k}k^{2n}\mathcal{P}_{C}(R_{s},R_{\mathrm{H}},k)\\\
=\int\limits_{0}^{\infty}\frac{\mathrm{d}k}{k}k^{2n}\frac{4(1+\omega)^{2}}{(5+3\omega)^{2}}\left(kR_{s}\right)^{4}\mathcal{P}_{\zeta}(k)\tilde{W}^{2}(R_{s},k)T^{2}(R_{\mathrm{H}},k),$
(45)

where $R_{S}$ is the smoothing scale, and $R_{\mathrm{H}}$ is the Hubble
horizon. For the _Press-Schechter_ -type approach, we will only need
$\sigma_{0}^{2}$, whilst higher order terms will be needed for the theory of
peaks.

Several things are note-worthy at this point. The first is that, whilst PBH
formation is an extremely non-linear process, $T(R_{\mathrm{H}},k)$ is the
linear transfer function. It is, therefore, preferable to complete
calculations in the super-horizon limit $R_{s}\ll R_{\mathrm{H}}$, such that
$\delta\ll 1$, and $T(R_{\mathrm{H}},k)\approx 1$. The second is that the
integrand contains the factor $k^{4}$, which can cause the integral to diverge
if there is no UV cut-off sharper than $k^{-4}$. This is the case for a
Gaussian smoothing, equation (21), which has an exponential cut-off, but not
necessarily for the top-hat smoothing, equation (20), which (though highly
oscillatory) drops as $k^{-4}$. Therefore, when the top-hat function is used,
the calculation is often performed at horizon crossing, $R_{S}=R_{\mathrm{H}}$
on a somewhat _ad hoc_ basis, and the effect of the transfer function then
guarantees that the integral converges. For further reading on this topic, it
is discussed in more detail in Young (2019), and the effect of non-linearities
at horizon crossing has been studied in De Luca  _et al._ (2023).

The initial abundance of PBHs is typically stated in terms of $\beta$, the
energy fraction of the universe collapsing into black holes at a given time
$t$,

$\beta(t)=\left.\Omega_{\mathrm{PBH}}(t)\right|_{\mathrm{formation}}=\left.\frac{\rho_{\mathrm{PBH}}}{\rho_{r}}\right|_{\mathrm{formation}},$
(46)

where $\rho_{\mathrm{PBH}}$ is the mean energy density of PBHs, $\rho_{r}$ is
the mean radiation density, and we are assuming complete radiation domination.
It will be helpful for consider $\beta$ as the fraction of the universe
collapsing to form PBHs at a time $t$ when the Hubble horizon is $R_{s}$ and
the horizon mass is $M_{\mathrm{H}}$, where these 3 parameters can be used
interchangeably. Note that the smoothing scale $R_{s}$ is used here, rather
than the (actual) Hubble horizon $R_{\mathrm{H}}$ at the time of calculation.
This is because the calculation typically works by comparing the _initial_
amplitude of perturbations to the critical amplitude, long before such
perturbations enter the horizon. The scale of the perturbations forming PBHs,
corresponding to the smoothing scale, is therefore much larger than the Hubble
scale at this time. However, as discussed in the previous paragraph, it is
often necessary to set $R_{s}=R_{\mathrm{H}}$.

Throughout this chapter, we will assume that PBHs form spontaneously upon
horizon re-entry. Accounting for the fact that PBHs take some time to form
after horizon re-entry would mean a larger value for $\beta$ (because
$\rho_{r}$ drops faster than $\rho_{\mathrm{PBH}}$ as the universe expands),
which would then be cancelled out by the relative changes in the PBH and
radiation densities as the universe expands when the final PBH abundance is
calculated at a later time. For simplicity, this is therefore neglected here,
but may be important if one is considering the formation of very light PBHs
which evaporate extremely rapidly.

###  IV.1 The _Press-Schechter_ -type approach

The idea behind the _Press-Schechter_ -type approach is that we can calculate
$\beta$, by calculating the fraction of the volume of the universe where the
compaction is above the threshold value. This gives us the collapse constraint

$n(\bar{\nu})=\delta_{\mathrm{D}}\left(\nu-\bar{\nu}\right)\Theta_{H}(\nu-\nu_{c}),$
(47)

where $\delta_{\mathrm{D}}$ is the Dirac-delta function,
$\nu_{c}=C_{c}/\sigma_{0}$, and

$\nu(\mathbf{x})=\frac{C(\mathbf{x})}{\sigma_{0}},$ (48)

has unit variance. This returns the number of regions (either 0 or 1) where
$\nu$ takes a value between $\bar{\nu}$ and $\bar{\nu}+\mathrm{d}\nu$, above
the critical value at point $\mathbf{x}$. To calculate $\beta$, we will want
the expectation value of $n$:

$\mathcal{N}(\bar{\nu})=\langle
n\rangle=\int\limits_{-\infty}^{\infty}\mathrm{d}\nu
n(\bar{\nu})P(\nu)=\begin{cases}P(\bar{\nu})\quad\quad\bar{\nu}\geq\nu_{c}\\\
0\quad\quad\quad\bar{\nu}<\nu_{c}\\\ \end{cases}.$ (49)

We can then calculate $\beta$ as

$\beta(t)=2\int\limits_{\nu_{c}}^{\infty}\mathrm{d}\nu\frac{M_{\mathrm{PBH}}}{M_{\mathrm{H}}}\frac{1}{\sqrt{2\pi}}\exp\left(-\frac{\nu^{2}}{2}\right),$
(50)

where the fraction $M_{\mathrm{PBH}}/M_{\mathrm{H}}$ gives the fraction of
each Hubble patch which collapses to form a PBH. The time dependance is
somewhat hidden here,

The simplest scenario is to assume that PBHs form with a fixed fraction of the
horizon mass, $M_{\mathrm{PBH}}/M_{\mathrm{H}}=\alpha$, whilst simulations of
PBH formation show that the PBH mass depends on both the scale and amplitude
of the perturbation from which it formed, which will be discussed in section
IV.1.1. The factor 2 arises as a “fudge-factor” from requiring that, for
$\nu_{c}\ll 1$, all of the energy universe should collapse to form PBHs,
$\beta=1$. Without the addition of the factor of 2, this integral would give
$\beta=0.5$ (for $\alpha=1$).

For $M_{\mathrm{PBH}}/M_{\mathrm{H}}=\alpha$, the expression for $\beta$ is
simply the complementary error function $\mathrm{erfc}$,

$\beta(t)=\alpha\mathrm{erfc}\left(\frac{\nu_{c}}{\sqrt{2}}\right).$ (51)

Using the large $x$ limit of $\mathrm{erfc}(x)$ gives,

$\beta(t)\approx\frac{\alpha}{\sqrt{2\pi}\nu_{c}}\exp\left(-\frac{\nu_{c}^{2}}{2}\right).$
(52)

Whilst not very accurate, this analytic expression does give an order-of-
magnitude estimate for the amplitude of the power spectrum required to form a
given abundance of PBHs, and illustrate the exponential dependence of $\beta$
on $\sigma_{0}$ (and thus, on the power spectrum). Because PBHs form from rare
perturbations in the extreme tail of the distribution, $\nu_{c}\gg 1$, small
changes in the power spectrum result in large changes to $\beta$.

####  IV.1.1 The critical scaling relationship

When a perturbation collapses to form a PBH, the final mass depends on how
much of the surrounding matter has been pulled in. For this to happen,
gravitational attraction towards the centre needs to overcome pressure forces
and the initial Hubble expansion of the surrounding material. The larger the
density at the centre, the greater the force of gravity, and the more material
is pulled into the final PBH. The result is that PBH mass depends on the
amplitude of the perturbation from which it formed, as well as the scale.
Simulations of PBH formation reveal that the PBH follows a power law known as
the critical scaling relation Musco  _et al._ (2009); Musco and Miller (2013);
Young  _et al._ (2019),

$M_{\mathrm{PBH}}=\mathcal{K}M_{\mathrm{H}}\left(C-C_{c}\right)^{\gamma},$
(53)

where the exact values depend on the specific profile considered, but we will
take $\mathcal{K}=4$, and $\gamma=0.36$ as representative values during
radiation domination. The horizon mass when the comoving Hubble horizon is
$R_{H}$, assuming radiation domination until the time of matter-radiation
equality, is given by

$M_{\mathrm{H}}(R_{\mathrm{H}})=M_{\mathrm{eq}}\left(\frac{g_{*}}{g_{\mathrm{\star,eq}}}\right)^{-1/6}\left(\frac{R_{\mathrm{H}}}{R_{\mathrm{eq}}}\right)^{2}M_{\mathrm{eq}}\approx
0.5\left(\frac{g_{*}}{10.75}\right)^{-1/6}\left(\frac{R_{\mathrm{H}}}{10^{-6}\mathrm{Mpc}}\right)^{2}M_{\odot},$
(54)

where $g_{*}$ is the number of relativistic degrees of freedom (although the
dependance on this is extremely weak and is typically neglected), the
subscript $\mathrm{eq}$ denotes the value at matter-radiation equality, and
$M_{\odot}$ is a solar mass. For the parameter choices in Nakama  _et al._
(2017), $M_{\mathrm{eq}}=2.8\times 10^{17}M_{\odot}$. Here, we have given the
horizon mass $M_{\mathrm{H}}$ as a function of the Hubble horizon $R_{H}$ at a
given time, but for our purposes it will be more convenient to consider what
the horizon mass will be at the time when the smoothing scale enters the
horizon, which is achieved simply by replacing $R_{\mathrm{H}}$ with $R_{s}$.
Therefore, the expression we will use for the horizon mass is therefore

$M_{\mathrm{H}}(M_{\mathrm{H}})\approx
0.5\left(\frac{g_{*}}{10.75}\right)^{-1/6}\left(\frac{R_{s}}{10^{-6}\mathrm{Mpc}}\right)^{2}M_{\odot}.$
(55)

Utilising equation (53), the expression for the initial PBH abundance becomes

$\beta(t)=2\int\limits_{\nu_{c}}^{\infty}\mathrm{d}\nu\mathcal{K}\left(\nu\sigma_{0}-\nu_{c}\sigma_{0}\right)^{\gamma}\frac{1}{\sqrt{2\pi}}\exp\left(-\frac{\nu^{2}}{2}\right),$
(56)

where the integral is now performed numerically.

Whilst not needed at this point, it will be helpful later to express $\beta$
as an integral over the PBH mass. Equations (53) and (48) can be used to
express $\nu$ as a function of the PBH mass

$\nu(M_{\mathrm{PBH}})=\frac{C(M_{\mathrm{PBH}})}{\sigma_{0}}=\frac{1}{\sigma_{0}}\left(\frac{M_{\mathrm{PBH}}}{\mathcal{K}M_{\mathrm{H}}}\right)^{1/\gamma}+\nu_{c},$
(57)

and differentiated to give

$\frac{\mathrm{d}\nu}{\mathrm{d}M_{\mathrm{PBH}}}=\frac{1}{\gamma\mathcal{K}M_{\mathrm{H}}\sigma_{0}}\left(\frac{M_{\mathrm{PBH}}}{\mathcal{K}M_{\mathrm{H}}}\right)^{\frac{1}{\gamma}-1}.$
(58)

Performing a change of variable, we can then express $\beta$ as

$\beta(t)=2\int\limits_{0}^{\infty}\mathrm{d}M_{\mathrm{PBH}}\left[\frac{1}{\gamma\mathcal{K}M_{\mathrm{H}}\sigma_{0}}\left(\frac{M_{\mathrm{PBH}}}{\mathcal{K}M_{\mathrm{H}}}\right)^{\frac{1}{\gamma}-1}\right]\frac{M_{\mathrm{PBH}}}{M_{\mathrm{H}}}\frac{1}{\sqrt{2\pi}}\exp\left(-\frac{\nu^{2}(M_{\mathrm{PBH}})}{2}\right).$
(59)

Including the effect of the critical scaling relationship means that, even
when only considering perturbations of a single scale (i.e. PBHs forming at a
single time), PBHs form with a range of masses. It also means that, in
addition to affecting the total abundance of PBHs which form, the amplitude of
the power spectrum also affects the masses at which they form. Figure 6 shows
the abundance of PBH forming at different scales, normalised by their overall
abundance, with values for the variance
$\sigma_{0}^{2}=\\{0.003,0.0065,0.01$}. These were chosen to correspond
approximately with typical values for the abundance, and give
$\beta\approx\\{4\times 10^{-20},4\times 10^{-10},5\times 10^{-7}\\}$. For
larger $\sigma_{0}^{2}$, the distribution peaks at larger masses - although
still very similar to the horizon mass. NB. This plot is independent of the
shape of the power spectrum. In the _Press-Schechter_ formalism, the power
spectrum only enters into the calculation via the $\sigma_{0}^{2}$ term -
which we have here simply chosen values for.

Figure 6: _Critical scaling and extended mass functions._ The figure shows the
masses of PBHs forming at a specific time (or, equivalently, at a specific
scale), normalised by the total abundance. The critical scaling relationship
between the PBH mass and the amplitude of perturbations means that PBHs form
at a range of masses, even when they only form at one specific time. The
larger the amplitude of perturbations, the larger the mass of the PBHs which
form.

####  IV.1.2 The total primordial black hole abundance

We have so far considered the abundance of PBHs at a single time of formation,
$\beta$. In order to calculate the total abundance of PBHs at a later time
(e.g. today), we need to sum the amount of PBHs forming at all times, and
account for the evolution of their density. Assuming that the universe behaves
purely as radiation domination until the time of matter-radiation equality,
then we can express the abundance of PBHs which formed when the horizon mass
was $M_{\mathrm{H}}$ as a function of $\beta$,

$\beta_{\mathrm{eq}}(M_{\mathrm{H}})=\frac{a_{\mathrm{eq}}}{a(M_{\mathrm{H}})}\beta(M_{\mathrm{H}})=\left(\frac{M_{\mathrm{eq}}}{M_{\mathrm{H}}}\right)^{1/2}\beta(M_{\mathrm{H}}).$
(60)

The total PBH abundance at matter-radiation equality is therefore found by
integrating over all times (scales) at which PBHs form,

$\Omega_{\mathrm{PBH}}=\int\limits_{M_{\mathrm{min}}}^{M_{\mathrm{max}}}\frac{\mathrm{d}M_{\mathrm{H}}}{M_{\mathrm{H}}}\left(\frac{M_{\mathrm{eq}}}{M_{\mathrm{H}}}\right)^{1/2}\beta(M_{\mathrm{H}}).$
(61)

The integral should technically be performed over all scales, but in practice
it is only necessary between the minimum and maximum scales being considered.
Typically, this range is set by the power spectrum under consideration, where
PBH formation occurs over a limited range of scales near the peak of the power
spectrum (or, between UV and IR cut-offs).

The total PBH abundance is often also expressed in terms of the fraction of
(cold) dark matter composed of PBHs

$f_{\mathrm{PBH}}=\frac{\Omega_{\mathrm{PBH}}}{\Omega_{\mathrm{CDM}}},$ (62)

where $\Omega_{\mathrm{CDM}}$ is the dark matter density parameter. This is
also useful because it is a time-independent quantity (after PBHs have stopped
forming).

####  IV.1.3 The mass function

In addition to the total abundance, a key potential observable of PBHs is the
mass function $f_{\mathrm{PBH}}$, which describes how many PBHs there are at
different masses. There are several different ways in which this can be
parameterised, but we will use the definition,

$\psi\left(M_{\mathrm{PBH}}\right)=\frac{1}{\Omega_{\mathrm{CDM}}}\frac{\mathrm{d}\Omega_{\mathrm{PBH}}}{\mathrm{d}(\ln
M_{\mathrm{PBH}})}=M_{\mathrm{PBH}}\frac{\mathrm{d}f_{\mathrm{PBH}}}{\mathrm{d}M_{\mathrm{PBH}}},$
(63)

which satisfies the normalisation condition,

$\int\mathrm{d}(\ln M_{\mathrm{PBH}})\psi(M_{\mathrm{PBH}})=f_{\mathrm{PBH}}.$
(64)

Using equations (59), (61) and (62) gives us the final expression for the mass
function

$\psi(M_{\mathrm{PBH}})=\frac{1}{\Omega_{\mathrm{CDM}}}\int\limits_{M_{\mathrm{min}}}^{M_{\mathrm{max}}}\frac{\mathrm{d}M_{\mathrm{H}}}{M_{\mathrm{H}}}\left(\frac{M_{\mathrm{eq}}}{M_{\mathrm{H}}}\right)^{1/2}\frac{2}{\gamma\mathcal{K}M_{\mathrm{H}}\sigma_{0}}\left(\frac{M_{\mathrm{PBH}}}{\mathcal{K}M_{\mathrm{H}}}\right)^{\frac{1}{\gamma}-1}\\\
\times\frac{M_{\mathrm{PBH}}}{M_{\mathrm{H}}}\frac{1}{\sqrt{2\pi}}\exp\left(-\frac{\nu^{2}(M_{\mathrm{PBH}})}{2}\right),$
(65)

where the expression for $\nu(M_{\mathrm{PBH}})$ is given in equation (57).
Recall that the variance $\sigma_{0}^{2}$ is a function of the smoothing
scale, which gives the integrand a dependence on the horizon mass
$M_{\mathrm{H}}$.

###  IV.2 The theory of peaks

The theory of peaks, developed by Bardeen, Bond, Kaiser, and Szalay (BBKS)
Bardeen  _et al._ (1986), extends upon the _Press-Schechter_ formalism by
stating that compact objects (PBHs in this case) will form at peaks in the
density above some threshold value. Our collapse constraint becomes

$n(\mathbf{x})=\delta_{D}^{(3)}(\mathbf{\eta})\Theta_{H}(\lambda_{3})\delta_{\mathrm{D}}^{(3)}\left(\nu-\bar{\nu}\right)\Theta_{H}(\nu-\nu_{c}),$
(66)

where the index $(3)$ indicates a 3-dimensional function, $\sigma_{n}$ are
moments of the power spectrum given by equation (45), $\mathbf{\eta}$ is given
by

$\mathbf{\eta}(\mathbf{x})=\frac{1}{\sigma_{1}}\mathbf{\nabla}C(\mathbf{x})$
(67)

$\zeta_{ij}$ is given by

$\zeta_{ij}(\mathbf{x})=\frac{1}{\sigma_{2}}\frac{\partial^{2}C(\mathbf{(}x))}{\partial
x_{i}\partial x_{j}},$ (68)

and finally, $\lambda_{3}$ is the smallest eigenvalue of $-\zeta_{ij}$. The
collapse constraint $n(\mathbf{x})$ returns the number of peaks (either 0 or
1) of height between $\bar{\nu}$ and $\bar{\nu}+\mathrm{d}\nu$, above the
critical value $\nu_{c}$ in the infinitesimal volume $\mathrm{d}^{3}x$,
divided by $\mathrm{d}^{3}x$. This is analogous to the collapse constraint
used in the _Press-Schechter_ -type approach, where the additional
$\delta_{D}^{(3)}(\mathbf{\eta})$ term constrains the spatial derivative of
$C$ to be zero, whilst the $\Theta_{H}(\lambda_{3})$ term specifies peaks
rather than troughs.

The collapse constraint depends on the spatial coordinate $\mathbf{x}$ only
through a set of continuous Gaussian fields evaluated at the same position,
which we will denote $\mathbf{y}$,

$\mathbf{y}=\\{\nu,\eta_{i},\zeta_{ij}\\}.$ (69)

In order to calculate the abundance of PBHs, we need to calculate the average
value for $n$, calculated by integrating over the PDF $P(\mathbf{y})$

$\mathcal{N}=\int\mathrm{d}\mathbf{y}\frac{\sigma_{2}^{3}}{\sigma_{1}^{3}}\left|\mathrm{det}\zeta_{ij}\right|nP(\mathbf{y}).$
(70)

The term
$\frac{\sigma_{2}^{3}}{\sigma_{1}^{3}}\left|\mathrm{det}\zeta_{ij}\right|$ is
included to account for the transformation from the Gaussian variables of the
constraint to the physical coordinates $\mathbf{x}$. Whilst not
straightforward to calculate, the integral gives the result

$\mathcal{N}(\nu)=\frac{\sigma_{2}^{3}}{3^{3/2}(2\pi)^{2}\sigma_{1}^{3}}\exp\left(-\frac{\nu^{2}}{2}\right)\Theta_{H}(\nu-\nu_{c})\\\
\times\int\limits_{0}^{\infty}\mathrm{d}uf(u)\frac{1}{\sqrt{2\pi}(1-\gamma_{02}^{2})}\exp\left(-\frac{(u-\gamma_{02}\nu)^{2}}{2(1-\gamma_{02}^{2})}\right),$
(71)

where $\gamma_{02}=\sigma_{1}^{2}/\sigma_{0}\sigma_{2}$ is the correlation
function, and the function $f(u)$ is given by

$f(u)=\sqrt{\frac{2}{5\pi}}\left[\left(\frac{u^{2}}{2}-\frac{8}{5}\right)\exp\left(-\frac{5u}{2}\right)+\left(\frac{31u^{2}}{4}+\frac{8}{5}\right)\exp\left(-\frac{5u}{8}\right)\right]\\\
+\frac{1}{2}\left(u^{3}-3u\right)\left[\mathrm{erf}\left(\sqrt{\frac{5}{2}}u\right)+\mathrm{erf}\left(\sqrt{\frac{5}{2}}\frac{u}{2}\right)\right],$
(72)

where $\mathrm{erf}(x)$ is the error function

$\mathrm{erf}(x)=\frac{2}{\pi}\int\limits_{0}^{x}\mathrm{d}t\exp(-t^{2}).$
(73)

For the full derivation, the interested reader is directed towards the
original BBKS paper Bardeen  _et al._ (1986), as well as section II and III of
Lazeyras  _et al._ (2016).

This expression is somewhat complicated, but fortunately can be greatly
simplified when one considers that PBHs form from high peaks, and one can then
make the assumption $v\gg 1$, which is known as the high-peak limit. The
integrand appearing in equation (71) is the function $f(u)$ multiplied by a
Gaussian PDF with mean $\gamma_{02}\nu$ and variance $(1-\gamma_{02}^{2})$,
where we generically have $\gamma_{02}\nu\gg(1-\gamma_{02}^{2})$. This means
that the PDF can be well approximated by a Dirac-delta function,
$\delta_{D}(u-\gamma_{02}\nu)$. Under this approximation, the integral is then
equal to $f(\gamma_{02}\nu)$, with $\gamma_{02}\nu\gg 1$. The exponential
terms in equation (72) are then very small, whilst the $\mathrm{erf}$ are
close to unity. We then have

$f(u)\approx u^{3}-3u\approx u^{3}.$ (74)

Inserting this into equation (71) gives the us the well-known (and much
simpler) number density of peaks in the high-peak limit:

$\mathcal{N}_{hi-
pk}(\nu)=\frac{1}{3^{3/2}(2\pi)^{2}}\frac{\sigma_{1}^{3}}{\sigma_{0}^{3}}\exp\left(-\frac{\nu^{2}}{2}\right).$
(75)

From this expression, we can easily that the factor
$\sigma_{1}^{3}/\sigma_{0}^{3}$ gives the expression units of $k^{3}\propto
1/\mathrm{volume}$, since $\mathcal{N}$ is a number density (rather than
simply a number, as in section IV.1).

We can now proceed with the calculation of the PBH abundance, which proceeds
in much the same way as before, except that, since $\mathcal{N}$ now gives a
number density, we now need to multiply by a volume to find the fraction of
the universe collapsing into PBHs at a given time (or equivelently, scale or
horizon mass):

$\beta(R_{s})=V(R_{s})\int\limits_{\nu_{c}}^{\infty}\mathrm{d}\nu\frac{M_{\mathrm{PBH}}}{M_{\mathrm{H}}}\frac{1}{3^{3/2}(2\pi)^{2}}\frac{\sigma_{1}^{3}}{\sigma_{0}^{3}}\exp\left(-\frac{\nu^{2}}{2}\right),$
(76)

where, again, the critical scaling relationship, equation (53), can be used to
relate the PBH mass to the amplitude of the perturbation which formed it. The
volume term $V(R_{s})$ is the volume of the window function being used,

$V(R_{s})=\begin{cases}(2\pi)^{3/2}R_{s}^{3}\quad\quad\mathrm{Gaussian~{}smoothing},\\\
\frac{4\pi}{3}R_{s}^{3}\quad\quad\quad\quad\mathrm{Top-hat~{}smoothing}.\\\
\end{cases}$ (77)

The total abundance and mass function can then be calculated using the same
method as in section IV.1.

It is noteworthy that the dominant term in the expression for $\beta$ remains
the unchanged exponential term. However, using the theory of peaks, rather
than _Press-Shechter_ generally predicts a higher abundance of PBHs, due to
the additional $v^{3}$ term appearing in the expression. However, when one
calculates the amplitude of the power spectrum required to produce a specified
PBH abundance, the difference is typically only of order $10\%$, although
significant changes can be seen in the mass function Gow  _et al._ (2021);
Young  _et al._ (2014).

###  IV.3 Extensions to the theory of peaks

When studying the large-scale structure of the universe, a similar approach is
taken to that considered here. Initial fluctuations in the density of the
universe grow over time and provide the seeds for the formation of compact
objects (such as galaxies and galaxy clusters, in this case), and one can
calculate the later abundance of e.g. galaxies by considering the initial
distribution of perturbations and the amplitude required for collapse. A
common problem encountered in this case is the cloud-in-cloud phenomenon,
which arises when a small-scale compact object is calculated to form inside a
large-scale compact object, resulting in double counting. A common method used
to solve this problem for large-scale structure is the excursion set approach,
which works by not adding the stipulation that haloes are only considered to
collapse to form a compact object if there are no haloes on larger scales
which will form compact objects at that scale.

The excursion set has been considered in the case of PBH formation, and it is
found that the cloud-in-cloud phenomenon, whereby a smaller mass PBH would be
absorbed by bigger mass PBHs forming later, is basically absent for a Gaussian
distribution Moradinezhad Dizgah  _et al._ (2019). This is due to the rarity
of PBHs at formation, meaning that only a very small fraction of small mass
PBHs form at locations where larger mass PBHs later form. 333However, this may
be different for non-Gaussian distributions, which could mean that PBHs form
at small masses at a location where there is a large-scale perturbation which
will collapse to a large mass PBH (although this does not occur for non-
Gaussianities arising as a result of non-linearities, which are considered in
section V).

However, a related problem to the cloud-in-cloud problem can be relevant for
PBHs, which is most easily seen if we consider a single perturbation, such as
that shown in figure 4. For a smoothing scale $R_{s}=r=0$, the compaction $C$
is also zero, and as the smoothing scale increases, $C$ peaks at the
characteristic scale $r_{m}$ before dropping back to zero. If the peak value
of $C$ is above the threshold for collapse, then the perturbation will
collapse to form a PBH, with a mass determined by the amplitude of $C$ at the
peak and the scale $r_{m}$ (as in equation (53)). However, this does mean that
the compaction is above the threshold for a range of smoothing scales (unless
the peak occurs at exactly $C=C_{c}$, in which case a vanishingly small PBH
forms). The calculations that we have considered thus far would, therefore,
consider that PBHs form at each of these scales, and sum over the relevant
range of scales. However, in reality, only a single PBH would form,
corresponding to the scale at which $C$ peaks.

This can be accounted for by extending the theory of peaks to stipulate that
PBHs form not only at locations where the compaction peaks, but also at the
(smoothing) scale at which the compaction peaks:

$n(\mathbf{x})=\delta_{D}\left(\frac{\partial C}{\partial
R}\right)\Theta_{H}\left(\frac{\partial^{2}C}{\partial
C^{2}}\right)\delta_{D}^{(3)}(\mathbf{\eta})\Theta_{H}^{(3)}(\lambda_{3})\delta_{\mathrm{D}}^{(3)}\left(\nu-\bar{\nu}\right)\Theta_{H}(\nu-\nu_{c}).$
(78)

The calculations necessary were developed simultaneously by Germani & Sheth
Germani and Sheth (2020) and Young & Musso Young and Musso (2020). The effect
of accounting for the improved the theory of peaks calculated was considered
by Gow et al Gow  _et al._ (2021), and it was found that the change to power
spectrum amplitude to produce a given PBH abundance only varies by a few
percent, and that the calculated mass function is also largely unchanged. For
this reason, the extensions to the theory of peaks are typically not included
in calculations.

##  V The effect of non-linearities

Until this point, we have been considering a linear relation between the
curvature perturbation $\zeta$ and the density contrast $\delta$, resulting
also in a linear relation to the compaction $C$. However, for PBHs to form,
large amplitude perturbations are needed, and it is therefore necessary to go
beyond linear order Young  _et al._ (2019).

The first effect we need to account for is a coordinate transform between the
comoving uniform density gauge in which we have defined $\zeta$ and the
comoving synchronous gauge in which we calculate $C$ and determine PBH
formation. We will still consider that PBHs form at high peaks which are
spherically symmetric, and therefore only consider the radial coordinate. The
radial coordinate $\hat{r}$ in the comoving synchronous gauge is then related
to the radial coordinate $r$ appearing in the comoving uniform density gauge
by

$\hat{r}=a(t)\exp\left(\zeta(r)\right)r.$ (79)

For small amplitude perturbations, $\hat{r}$ increases monotonically with $r$,
and these are known as type I perturbations. However, as larger $\zeta$ is
considered, there comes a point where this is longer true, meaning that
increasing coordinate radius $r$ does not coincide with increasing physical
distance. Perturbations where this occurs are known as type II perturbations,
and their evolution (and their subsequent collapse to PBHs) is not currently
well understood. We can find the point at which the switch between type I and
II perturbations occurs by searching for a stationary point
$\mathrm{d}\hat{r}/\mathrm{d}r=0$. Stationary points occur when we have

$-r\zeta^{\prime}(r)=1.$ (80)

For any given perturbation, this function peaks at the scale $r_{m}$ (see the
discussion in section III), and we therefore have the condition for type II
perturbations,

$-r_{m}\zeta^{\prime}(r_{m})>1.$ (81)

The second effect we need to consider is the non-linear relationship between
$\zeta$ and $\delta$, and we will again assume that $\zeta$ follows a Gaussian
distribution. In the super-horizon regime, $\delta$ can be expressed as a
function of $\zeta$ as

$\delta(\hat{r})=-\frac{4(1+\omega)}{5+3\omega}\left(\frac{1}{aH}\right)^{2}\exp\left(-\frac{5\zeta(r)}{2}\right)\nabla^{2}\exp\left(\frac{\zeta(r)}{2}\right),$
(82)

where $\hat{r}$ is the radial coordinate in the comoving synchronous gauge,
and in spherical symmetry

$\nabla^{2}=\frac{\partial^{2}}{\partial
r^{2}}+\frac{2}{r}\frac{\partial}{\partial r}.$ (83)

We note that equation (82) is exact to all orders in $\zeta$, but is correct
to second order in $\epsilon=1/HL$, where $L$ is the length scale
characterizing the perturbation, and this expression is therefore valid in the
super-horizon regime.

Combining these equations and inserting them into the equation for the
compaction with a top-hat window function, equation (27) gives

$C(r)=-\frac{2}{ae^{\zeta(r)}r}\frac{3H^{2}}{8\pi}\int\limits_{0}^{r}\mathrm{d}\left(ae^{\zeta(r)}r\right)\left[4\pi\left(ae^{\zeta(r)}r\right)^{2}\right]\\\
\times\frac{2(1+\omega)}{5+3\omega}\left(\frac{1}{aH}\right)^{2}e^{-2\zeta(r)}\left[\zeta^{\prime\prime}(r)+\frac{2}{r}\zeta^{\prime}(r)+\frac{1}{2}\left(\zeta^{\prime}(r)\right)^{2}\right],$
(84)

where the prime denotes a derivative with respect to $r$. Fortunately, this
complicated looking integral has a simple analytic solution:

$C(r)=-\frac{3(1+\omega)}{5+3\omega}\left[2r\zeta^{\prime}(r)-\left(r\zeta^{\prime}(r)\right)^{2}\right].$
(85)

The linear component of this expression,

$C_{1}(r)=\frac{6(1+\omega)}{5+3\omega}r\zeta^{\prime}(r),$ (86)

is exactly the same as the expression for the compaction seen in equation
(31), and the variance $\langle C_{1}^{2}\rangle=\sigma_{0}^{2}$ is calculated
exactly as before, as in equation (45). Setting $\omega=1/3$ during radiation
domination, we have a simple quadratic relation between the compaction $C$ and
its linear component $C_{1}$:

$C=C_{1}-\frac{3}{8}C_{1}^{2},$ (87)

which takes a maximum value $C=2/3$ when $C_{1}=(4/3)r\zeta^{\prime}(r)=4/3$.
If we consider the characteristic scale $r_{m}$, where the compaction peaks,
this corresponds to the switch between type I and II perturbations, seen in
equation (81).

For the large, rare perturbations which form PBHs, peaks in the non-linear
compaction field will correspond to peaks in the linear compaction field Young
_et al._ (2019), and we can therefore simply apply the same calculation
methods we have used before, with some small modifications. We will consider
that PBHs form for type I perturbations with $C>C_{c}$. Solving equation (87)
for $C_{1}$ gives

$C_{1\pm}=\frac{4}{3}\left(1\pm\sqrt{\frac{2-3C}{2}}\right),$ (88)

which gives two solutions for any value $C<2/3$, where $C_{1-}$ and $C_{1+}$
correspond to type I and II perturbations respectively. Going forwards, we
will consider only type I solutions, and will drop the minus sign in the
subscript. Taking $C_{c}\approx 0.5$ (for Mexican hat-type perturbations)
means that we will consider PBHs to form in the range,

$C_{1,c}\approx 0.67<C<\frac{4}{3}.$ (89)

Due to the exponential suppression in the abundance of large peaks, the upper
limit will have hardly any effect on the total abundance of PBHs (although we
will need to account for it in the calculation nonetheless, to avoid the
appearance of imaginary solutions).

Neglecting the upper bound for the moment, therefore, we can obtain a simple
estimate for the effect of non-linearities using equation (52), using
$C_{1,c}\approx 0.67$ instead of $C_{c}=0.5$. The result is a severe reduction
in the number of PBHs forming, typically by numerous orders of magnitude. In
order to achieve the same PBH abundance as before, then $\sigma_{0}^{2}$ and
the amplitude of the power spectrum needs to be larger by a factor,

$\left(\frac{C_{1,c}}{C_{c}}\right)^{2}=\left(\frac{4}{3C_{c}}\left(1\pm\sqrt{\frac{2-3C_{c}}{2}}\right)\right)^{2}\approx
1.78,$ (90)

where we have used $C_{c}=0.5$ in the final approximation on the right-hand
side. This means that, when constraints on the PBH abundance are translated
into constraints on the primordial power spectrum, they are almost a factor of
2 weaker than if non-linearities were ignored.

The main differences we will need to make in the calculation of the PBH
abundance concerns the critical scaling relationship, and the limits of
several integrals. Inserting equation (87) into the equation for the critical
scaling relation, equation (53), gives the PBH mass as a function of $C_{1}$

$M_{\mathrm{PBH}}=\mathcal{K}M_{\mathrm{H}}\left(C_{1}-\frac{3}{8}C_{1}^{2}-C_{c}\right)^{\gamma},$
(91)

where the values are unchanged (we previously chose $\mathcal{K}=4$ and
$\gamma=0.36$ as representative values during radiation domination). Because
the compaction takes a maximum value $C_{max}$ at $C_{1,max}=4/3$, we now have
a maximum mass PBH which can form for a given horizon mass,

$M_{\mathrm{PBH,max}}=\mathcal{K}M_{\mathrm{H}}\left(C_{max}-C_{c}\right)^{\gamma}\approx
2.1M_{\mathrm{H}},$ (92)

or, equivalently, the minimum horizon mass which can form a PBH of a certain
mass,

$M_{\mathrm{H,min}}=\frac{M_{\mathrm{PBH}}}{\mathcal{K}\left(C_{max}-C_{c}\right)^{\gamma}}\approx
0.48M_{\mathrm{PBH}}.$ (93)

Inverting equation (91) gives us $\nu$ as a function of the PBH mass,

$\nu(M_{\mathrm{PBH}})=\frac{C_{1}(M_{\mathrm{PBH}})}{\sigma_{0}}=\frac{2}{3\sigma_{0}}\left(2-\sqrt{4-6C_{c}-6\left(\frac{M_{\mathrm{PBH}}}{\mathcal{K}M_{\mathrm{H}}}\right)^{1/\gamma}}\right),$
(94)

where we have kept only the solution for type I perturbations, and
$\sigma_{n}$ is given by equation (45). We will also require the derivative,

$\frac{\mathrm{d}\nu}{\mathrm{d}M_{\mathrm{PBH}}}=\frac{\sqrt{2}\left(\frac{M_{\mathrm{PBH}}}{\mathcal{K}M_{\mathrm{H}}}\right)^{\frac{1}{\gamma}-1}}{\gamma\mathcal{K}M_{\mathrm{H}}\sqrt{2-3C_{c}-3\left(\frac{M_{\mathrm{PBH}}}{\mathcal{K}M_{\mathrm{M_{H}}}}\right)^{\frac{1}{\gamma}}}}.$
(95)

Using the theory of peaks, the initial abundance of PBHs collapsing to form
PBHs is given by,

$\beta(r)=\frac{4\pi
r^{3}}{3}\int\limits_{\nu_{c}}^{\frac{4}{3\sigma_{0}}}\mathrm{d}\nu\mathcal{K}\left(\nu\sigma_{0}-\frac{3}{8}\left(\nu\sigma_{0}\right)^{2}-C_{c}\right)^{\gamma}\frac{1}{3^{3/2}(2\pi)^{2}}\frac{\sigma_{1}^{3}}{\sigma_{0}^{3}}\exp\left(-\frac{\nu^{2}}{2}\right),$
(96)

We have here included the upper limit for PBH formation, corresponding to
$C_{\mathrm{max}}=4/3$, although this typically has a minimal effect on the
PBH abundance.

Again, the total PBH abundance $f_{\mathrm{PBH}}$ and the mass function
$\psi_{\mathrm{PBH}}$ can then be obtained following the same method as
described in section IV.1;

$f_{\mathrm{PBH}}=\frac{\Omega_{\mathrm{PBH}}}{\Omega_{\mathrm{CDM}}}=\frac{1}{\Omega_{\mathrm{CDM}}}\int\limits_{M_{\mathrm{min}}}^{M_{\mathrm{max}}}\frac{\mathrm{d}M_{\mathrm{H}}}{M_{\mathrm{H}}}\left(\frac{M_{\mathrm{eq}}}{M_{\mathrm{H}}}\right)^{1/2}\beta(M_{\mathrm{H}}),$
(97)

where equation (54) relates the horizon mass $M_{\mathrm{H}}$ to the
(smoothing) scale $r$. The integral runs over the entire range of horizon
masses for which PBH formation occurs.

The PBH mass function is, likewise, calculated in the same fashion as
previously,

$\psi(M_{\mathrm{PBH}})=M_{\mathrm{PBH}}\frac{\mathrm{d}f_{\mathrm{PBH}}}{\mathrm{d}M_{\mathrm{PBH}}}=\frac{M_{\mathrm{PBH}}}{\Omega_{\mathrm{CDM}}}\int\limits_{M_{\mathrm{H,min}}}^{M_{\mathrm{max}}}\frac{\mathrm{d}M_{\mathrm{H}}}{M_{\mathrm{H}}}\left(\frac{M_{\mathrm{eq}}}{M_{\mathrm{H}}}\right)^{1/2}\\\
\times\frac{4\pi
r^{3}(M_{\mathrm{H}})}{3}\frac{\mathrm{d}\nu}{\mathrm{d}M_{\mathrm{PBH}}}\frac{M_{\mathrm{PBH}}}{M_{\mathrm{H}}}\frac{1}{3^{3/2}(2\pi)^{2}}\frac{\sigma_{1}^{3}}{\sigma_{0}^{3}}\exp\left(-\frac{\nu^{2}(M_{\mathrm{PBH}})}{2}\right).$
(98)

It is here important to note that we now integrate from the minimum horizon
mass which forms a PBH of mass $M_{\mathrm{PBH}}$, given by equation (93). It
is now important not to neglect this limit, as the peak height $\nu$ becomes
imaginary for smaller values.

##  VI Summary

We have discussed how cosmological perturbations collapse to form PBHs at
horizon entry, with larger mass PBHs corresponding to larger scale
perturbations formed at earlier times during cosmological inflation. Because
PBHs form deep in the radiation dominated epoch, and are constrained to be
less abundant than dark matter, they are very rare upon formation, making up
only a small fraction of the energy content of the universe. The abundance of
PBHs depends upon the amplitude of perturbations in the early universe, with
larger amplitude perturbations resulting in a higher PBH abundance. This means
that the PBH abundance can be calculated from the primordial power spectrum,
providing a unique method to constrain the early universe via constraints on
the PBH abundance. An accurate calculation of the PBH abundance also means
that predictions can be made for other interesting factors, such as the mass
function and initial clustering.

There are numerous ways in which perturbations can be parameterised, including
the curvature perturbation, the density contrast, but in the context of PBH
formation, it is typically most useful to consider perturbations described by
the compaction. A perturbation is considered to collapse to form a PBH if the
compaction peaks above a certain threshold value, $C_{c}$. The compaction
function is conserved on super-horizon scales, and the statistics are
relatively simple, even when non-linearities are accounted for. However, as
perturbations re-enter the horizon, their evolution becomes significantly more
complicated - although a linear transfer function is often used to describe
the evolution on horizon and sub-horizon scales, the validity of which can be
questioned De Luca  _et al._ (2023). It is therefore preferable to calculate
the PBH abundance in the super-horizon regime by comparing the initial
amplitude of perturbations (long before horizon entry) to the critical value.
However, depending on the power spectrum being considered, this is not always
possible.

The most important parameter for determining the PBH abundance is the the
threshold value for collapse. The critical value is sufficiently large that no
PBHs are expected to form in the universe, if the power spectrum is not
significantly larger on small scales than is visible in the CMB. The exact
value of the threshold depends upon the profile shape of a given perturbation,
but since perturbations are generally expected to be close to the Mexican-hat
profile, it is safe to use the critical value corresponding to that profile,
$C_{c}=0.5$.

Several methods for calculating the PBH abundance and mass function have been
described, with the theory of peaks being the most common method. Extensions
to the theory of peaks do exist, but these typically do not offer a
significant improvement in the accuracy of such calculations, and can be far
more complex.

Throughout this chapter, we have assumed an initial Gaussian distribution for
the curvature perturbation $\zeta$, but in section V the effect of non-
linearity between $\zeta$ and the density contrast is considered. This non-
linearity results in the compaction following a non-Gaussian distribution
function, which significantly increases the amplitude of the power spectrum
required to form a given abundance of PBHs.

## References

  * Zel’dovich and Novikov (1967) Y. B. Zel’dovich and I. Novikov, Soviet Astronomy 10, 602 (1967).
  * Hawking (1971) S. Hawking, Mon. Not. Roy. Astron. Soc. 152, 75 (1971).
  * Carr (1975) B. J. Carr, Astrophys. J. 201, 1 (1975).
  * Musco  _et al._ (2005) I. Musco, J. C. Miller, and L. Rezzolla, Class. Quant. Grav. 22, 1405 (2005), arXiv:gr-qc/0412063 .
  * Musco  _et al._ (2009) I. Musco, J. C. Miller, and A. G. Polnarev, Class. Quant. Grav. 26, 235001 (2009), arXiv:0811.1452 [gr-qc] .
  * Nakama  _et al._ (2014a) T. Nakama, T. Harada, A. G. Polnarev, and J. Yokoyama, JPS Conf. Proc. 1, 013115 (2014a).
  * Nakama  _et al._ (2014b) T. Nakama, T. Harada, A. G. Polnarev, and J. Yokoyama, in  _10th International Symposium on Cosmology and Particle Astrophysics_ (2014) arXiv:1401.7740 [gr-qc] .
  * Akrami  _et al._ (2020) Y. Akrami _et al._ (Planck), Astron. Astrophys. 641, A10 (2020), arXiv:1807.06211 [astro-ph.CO] .
  * Carr  _et al._ (2021) B. Carr, K. Kohri, Y. Sendouda, and J. Yokoyama, Rept. Prog. Phys. 84, 116902 (2021), arXiv:2002.12778 [astro-ph.CO] .
  * Gow  _et al._ (2021) A. D. Gow, C. T. Byrnes, P. S. Cole, and S. Young, JCAP 02, 002 (2021), arXiv:2008.03289 [astro-ph.CO] .
  * Drees and Erfani (2011) M. Drees and E. Erfani, JCAP 1104, 005 (2011), arXiv:1102.2340 [hep-ph] .
  * Bugaev and Klimai (2014) E. Bugaev and P. Klimai, Phys. Rev. D90, 103501 (2014), arXiv:1312.7435 [astro-ph.CO] .
  * Özsoy  _et al._ (2018) O. Özsoy, S. Parameswaran, G. Tasinato, and I. Zavala, JCAP 1807, 005 (2018), arXiv:1803.07626 [hep-th] .
  * Garcia-Bellido  _et al._ (1996) J. Garcia-Bellido, A. D. Linde, and D. Wands, Phys. Rev. D54, 6040 (1996), arXiv:astro-ph/9605094 [astro-ph] .
  * Lyth (2012) D. H. Lyth, JCAP 1205, 022 (2012), arXiv:1201.4312 [astro-ph.CO] .
  * Bugaev and Klimai (2012) E. Bugaev and P. Klimai, Phys. Rev. D85, 103504 (2012), arXiv:1112.5601 [astro-ph.CO] .
  * Shibata and Sasaki (1999) M. Shibata and M. Sasaki, Phys. Rev. D 60, 084002 (1999), arXiv:gr-qc/9905064 .
  * Wands  _et al._ (2000) D. Wands, K. A. Malik, D. H. Lyth, and A. R. Liddle, Phys. Rev. D 62, 043527 (2000), arXiv:astro-ph/0003278 .
  * Josan  _et al._ (2009) A. S. Josan, A. M. Green, and K. A. Malik, Phys. Rev. D 79, 103520 (2009), arXiv:0903.3184 [astro-ph.CO] .
  * Musco (2019) I. Musco, Phys. Rev. D 100, 123524 (2019), arXiv:1809.02127 [gr-qc] .
  * Young (2019) S. Young, Int. J. Mod. Phys. D 29, 2030002 (2019), arXiv:1905.01230 [astro-ph.CO] .
  * Bardeen  _et al._ (1986) J. M. Bardeen, J. R. Bond, N. Kaiser, and A. S. Szalay, Astrophys. J. 304, 15 (1986).
  * Young (2022) S. Young, JCAP 05, 037 (2022), arXiv:2201.13345 [astro-ph.CO] .
  * Carr and Hawking (1974) B. J. Carr and S. W. Hawking, Mon. Not. Roy. Astron. Soc. 168, 399 (1974).
  * Harada and Jhingan (2016) T. Harada and S. Jhingan, PTEP 2016, 093E04 (2016), arXiv:1512.08639 [gr-qc] .
  * Escrivà  _et al._ (2021) A. Escrivà, C. Germani, and R. K. Sheth, JCAP 01, 030 (2021), arXiv:2007.05564 [gr-qc] .
  * Polnarev and Musco (2007) A. G. Polnarev and I. Musco, Class. Quant. Grav. 24, 1405 (2007), arXiv:gr-qc/0605122 .
  * Nakama  _et al._ (2014c) T. Nakama, T. Harada, A. G. Polnarev, and J. Yokoyama, JCAP 01, 037 (2014c), arXiv:1310.3007 [gr-qc] .
  * Escrivà (2020) A. Escrivà, Phys. Dark Univ. 27, 100466 (2020), arXiv:1907.13065 [gr-qc] .
  * Escrivà  _et al._ (2020) A. Escrivà, C. Germani, and R. K. Sheth, Phys. Rev. D 101, 044022 (2020), arXiv:1907.13311 [gr-qc] .
  * Germani and Musco (2019) C. Germani and I. Musco, Phys. Rev. Lett. 122, 141302 (2019), arXiv:1805.04087 [astro-ph.CO] .
  * De Luca  _et al._ (2023) V. De Luca, A. Kehagias, and A. Riotto, Phys. Rev. D 108, 063531 (2023), arXiv:2307.13633 [astro-ph.CO] .
  * Musco and Miller (2013) I. Musco and J. C. Miller, Class. Quant. Grav. 30, 145009 (2013), arXiv:1201.2379 [gr-qc] .
  * Young  _et al._ (2019) S. Young, I. Musco, and C. T. Byrnes, JCAP 11, 012 (2019), arXiv:1904.00984 [astro-ph.CO] .
  * Nakama  _et al._ (2017) T. Nakama, J. Silk, and M. Kamionkowski, Phys. Rev. D 95, 043511 (2017), arXiv:1612.06264 [astro-ph.CO] .
  * Lazeyras  _et al._ (2016) T. Lazeyras, M. Musso, and V. Desjacques, Phys. Rev. D 93, 063007 (2016), arXiv:1512.05283 [astro-ph.CO] .
  * Young  _et al._ (2014) S. Young, C. T. Byrnes, and M. Sasaki, JCAP 07, 045 (2014), arXiv:1405.7023 [gr-qc] .
  * Moradinezhad Dizgah  _et al._ (2019) A. Moradinezhad Dizgah, G. Franciolini, and A. Riotto, JCAP 11, 001 (2019), arXiv:1906.08978 [astro-ph.CO] .
  * Germani and Sheth (2020) C. Germani and R. K. Sheth, Phys. Rev. D 101, 063520 (2020), arXiv:1912.07072 [astro-ph.CO] .
  * Young and Musso (2020) S. Young and M. Musso, JCAP 11, 022 (2020), arXiv:2001.06469 [astro-ph.CO] .