#
Primordial black holes
and their gravitational-wave signatures

Eleni Bagui  , Sebastien Clesse111Project coordinator:
<EMAIL_ADDRESS>, Valerio De Luca  , Jose María Ezquiaga  ,
Gabriele Franciolini  , Juan García-Bellido222Project coordinator:
<EMAIL_ADDRESS>, Cristian Joana  , Rajeev Kumar Jain  ,
Sachiko Kuroyanagi  , Ilia Musco  , Theodoros Papanikolaou  , Alvise
Raccanelli  , Sébastien Renaux-Petel  , Antonio Riotto  , Ester Ruiz
Morales  , Marco Scalisi  , Olga Sergijenko  , Caner Unal  , Vincent
Vennin  , David Wands
(For the LISA Cosmology Working Group)

###### Abstract

In the recent years, primordial black holes (PBHs) have emerged as one of the
most interesting and hotly debated topics in cosmology. Among other
possibilities, PBHs could explain both some of the signals from binary black
hole mergers observed in gravitational wave detectors and an important
component of the dark matter in the Universe. Significant progress has been
achieved both on the theory side and from the point of view of observations,
including new models and more accurate calculations of PBH formation,
evolution, clustering, merger rates, as well as new astrophysical and
cosmological probes. In this work, we review, analyse and combine the latest
developments in order to perform end-to-end calculations of the various
gravitational wave signatures of PBHs. Different ways to distinguish PBHs from
stellar black holes are emphasized. Finally, we discuss their detectability
with LISA, the first planned gravitational-wave observatory in space.

##  1 Introduction

###  1.1 A brief history of primordial black holes (PBHs)

The idea that black holes could have been formed in the early Universe dates
back to the late 1960’s with the pioneering work of Zel’dovich and Novikov [1]
and to the early 1970’s with the work of Hawking [2]. This triggered Hawking’s
famous discovery of black hole evaporation [3]. Subsequently, Carr, Hawking’s
PhD student at the time, continued to investigate PBHs [4, 5]. Already in [2,
4, 6], it was suggested that such primordial black holes (PBHs) could
contribute to the suspected Dark Matter (DM) in the Universe or to the seeds
of supermassive black holes. The first formation scenarios in the context of
inflation were proposed in the 1990’s [7, 8, 9, 10, 11, 12, 13, 14, 15, 16,
17], but these usually led to (evaporating) PBHs of very small mass, except
for [11] who predicted solar mass PBHs that could account for the MACHO
microlensing events observed towards the Magellanic clouds [18, 19]. Other
mechanisms were also proposed, e.g., based on phase transitions [20, 21], an
early matter era [22, 23, 24], scalar field instabilities [25], collapse of
topological defects [26], modified gravity [27, 28, 29], string theory [30,
31, 32], bouncing cosmological scenarios [33, 34, 35, 36, 37] as well as
quantum gravity setups [38]. In the late 1990’s stellar-mass PBHs were
seriously considered as a dark matter candidate, following the possible
detection in the MACHO survey [18, 19] of several microlensing events towards
the Magellanic clouds. However, the EROS [39] and OGLE [40, 41, 42, 43]
surveys later set more stringent limits on the PBH abundance, and at the same
time, very stringent constraints from cosmic microwave background (CMB)
observations were claimed in [44].

Despite some pioneering numerical studies [45, 46, 47], only more recently it
has become possible to fully understand the mechanism of PBH formation with
detailed spherically-symmetric numerical simulations [48, 49, 50], showing
that superhorizon cosmological perturbations would collapse to PBHs after re-
entering the cosmological horizon, if their amplitude $\delta$ is larger than
a certain threshold value $\delta_{c}$. This quantity, measured at horizon
crossing, was initially estimated with a simplified Jeans length argument in
Newtonian gravity [5], giving $\delta_{c}\sim c_{s}^{2}$, where
$c_{s}^{2}=1/3$ is the sound speed of the cosmological radiation fluid
measured in units of the speed of light. More recently this argument has been
generalized within general relativity using a three-zone model [51, 52] which
gives $\delta_{c}\simeq 0.4$ for a radiation-dominated Universe. Numerical
simulations showed that the mechanism of critical collapse characterizing
perfect fluids [53] arises naturally in the context of PBH formation [16, 21,
54, 49]. In particular, the collapse is characterized by a relativistic wind
which progressively separates the collapsing perturbation from the expanding
background, approaching a self-similar critical solution characterizing the
threshold [55, 56]. Finally, a very recent detailed study has found a clear
relation between the value of the threshold, $\delta_{c}$, and the initial
configuration of the curvature (or energy density) profile used to set the
initial condition of the numerical simulations, with $0.4\leq\delta_{c}\leq
2/3$, where the shape is identified by a single parameter [57, 58]. This range
is reduced to $0.4\leq\delta_{c}\lesssim 0.6$ when the initial perturbations
are computed, using peak theory [59], from the primordial power spectrum of
cosmological perturbations [60].

Since 2016, the real game-changer that has rekindled the idea that PBHs may
exist and constitute a significant fraction to the total dark matter [61, 62,
63] has been the first gravitational-wave (GW) detection from a black hole
merger by LIGO/Virgo [64]. Initially, the merger rate of early PBH binaries
seemed to only allow a small fraction of the dark matter to be made of PBHs
[63], but N-body simulations later showed that these rates are in fact
significantly suppressed [65, 66]. Nowadays, the importance of the different
PBH binary formation channels, the possible abundance of PBHs, their viable
mass function, are a subject of intense activity and are hotly debated (for
recent reviews, see e.g. [67, 68]). Recent analysis based on the rate and mass
[69, 65, 70, 71, 72, 73, 74, 75, 76], or spin [77, 78, 79], or both [80, 81,
82, 83, 84, 85, 86, 87] distributions of compact binary coalescences observed
in the first, second [88] and third [89, 90] observing runs of the
LIGO/Virgo/Kagra collaboration (LVK), suggest that some black holes may be
primordial. Some are based on a Bayesian approach [69, 65, 71, 75, 81, 91, 92,
82, 83, 84, 87, 77, 79] that will ultimately allow us to compare statistically
the different possible origins and scenarios. However due to the broad variety
of PBHs and astrophysical models as well as the large theoretical and
observational uncertainties, it is still premature to firmly affirm (or deny)
a PBH origin of GW observations. In this context, the next generation of GW
detectors like Einstein Telescope [93], Cosmic Explorer [94] and the Laser
Interferometer Space Antenna (LISA) [95] will play a crucial role. Also, the
development of ultra-high frequency GW detectors may additionally complement
the search for GW signatures of asteroid-mass PBHs [96, 97] (see also e.g.
[98, 99, 100, 101, 102]).

Important progress has also been made on the limits on the allowed PBH
abundance set by various types of observations. Shortly after the first GW
detection, the CMB limits were re-analysed and found to be much less stringent
than initially thought [103], thereby re-opening the stellar mass region. At
the same time, the degree of validity of the microlensing limits has been
questioned [104, 105, 106, 107, 108, 70] for more realistic scenarios than
monochromatic and homogeneous distributions of PBHs. New probes, like ultra-
faint-dwarf galaxies [109, 110, 111] and X-ray sources towards the galactic
center [112], have revived the mass region between $1$ and $100M_{\odot}$,
while the constraining power of other probes like high-cadence microlensing
[113], neutron stars [114] and white dwarves [115], have been reduced, which
has reopened the asteroid-mass region [116]. Finally, a series of recent
microlensing events from OGLE [117, 118], HSC [113], as well as quasars [119,
120, 121, 122], may hint at the existence of an important population of
planetary-mass and solar-mass dark compact objects. The unexplained spatial
correlations between infrared and X-ray backgrounds at high redshift [123]
could also be due to an important population of stellar-mass PBHs [124, 125,
126].

Furthermore, since PBHs are formed by the collapse of large density
fluctuations, there is an associated stochastic GW background (SGWB) sourced
by these perturbations at second-order. It has been calculated [127, 128] that
if BHs detected by LVK have a primordial origin, there is an inevitable
accompanying SGWB peaking around pulsar-timing-array frequencies. In September
2020, NANOGrav claimed the possible detection of a SGWB at nanohertz frequency
using pulsar timing arrays [129], which may have been sourced by primordial
density perturbations which could be the origin of stellar-mass PBH formation
[130, 131, 132]. Remarkably, such claim was confirmed by more recent pulsar
timing array data released in 2023 by the NANOGrav [133, 134], EPTA (in
combination with InPTA) [135, 136, 137], PPTA [138, 139, 140] and CPTA [141]
collaborations, who found evidence for a Hellings and Downs angular
correlation, typical of an homogeneous spin-2 GW background and consistent
with the quadrupolar nature of GWs in GR [142]

On the theory side, a plethora of new models have been proposed in the last
four years. Most of them are still subject to fine-tuning issues related to
the amplitude of fluctuations required to form PBHs with a significant
abundance. Nevertheless, one should note a new class of models that do not
require a modification of the primordial power spectrum on the cosmological
scale, the PBH formation being related to a modification of the Gaussian tail
of the distribution of primordial perturbations [143, 144]. In these models, a
large population of PBHs could arise much more naturally. PBH masses of order
of the solar mass and around $30M_{\odot}$ could also arise naturally due to
the thermal history of the Universe and the change in the equation of state at
the QCD epoch [20, 21, 145, 70, 87, 76], which should strongly boost PBH
formation. Finally, besides explaining the dark matter or a significant
fraction of it, PBHs could also be related to baryogenesis and different
mechanisms have been proposed [146, 147, 148, 144, 149, 150].

Any firm detection would open a new window on the physics at play in the very
early Universe and a possible way to solve various long-standing astrophysical
and cosmological puzzles [69, 70]. Ultimately, the best way to distinguish
PBHs from stellar black holes would be to detect sub-solar-mass black holes.
Recently, four sub-solar black hole triggers have been found in the second
observing run of LVK [151, 152], and three in the third observing run [153].
These are not statistically significant enough to claim a firm detection, but
if confirmed, these would strongly support a PBH origin of some LVK black
holes and make a significant contribution to the dark matter in the Universe.

###  1.2 Why a(nother) review on PBHs?

There already exist a few recent and relatively complete review articles
dedicated to PBHs, see e.g. [154, 155, 67, 68, 156]. Earlier reviews that pre-
date the first GW detection also exist, see e.g. [157, 158]. All of them are
written by only a few authors. It has therefore been quite challenging for
them to include detailed discussions and analysis of the broad range of topics
covered by PBHs. Indeed, these include early Universe cosmology, numerical
relativity, dark matter, astrophysics, celestial mechanics, GW astrophysics,
etc. Those reviews therefore summarize some of the results obtained in these
topics but do not include model specific A-to-Z calculations of PBH related
observables. Very often, some of the most recent developments like the precise
determination of the critical threshold density leading to PBH formation,
which is highly model-dependent and varies with the evolution of the equation
of state of the Universe, or the rate suppression/boost linked to Poisson
clustering, are not always considered when models are compared to
observations.

The goal of this paper is to review these most recent developments and discuss
how they impact the GW signatures of PBHs. Compared to other recent reviews,
we aim at integrating together the most recent and accurate models, e.g., of
PBH formation, evolution, clustering, merger rates, in order to compute the GW
signatures of PBHs, and to explore their detectability. For this purpose, we
have developed in parallel a numerical code that implements all those models
and recent scenarios. The code will be released soon in a separate
publication, but it has been already used in this paper in order to produce
the key figures. The task of including these increasingly complex models has
been eased by the specific expertise of many authors, who have actively
contributed to these developments, covering a broad range of topics.

Finally, another original aspect of this work is the fact that we plan to have
a living review format, with bi-yearly updates of the present manuscript, in
order to track the rapid progress in this very active field of research and
provide at all times an up-to-date reference publication. In future, this
review could be expanded to include new sections on the detectability by PTAs
and ground-based GW detectors, or on the cosmological signatures of PBH
scenarios.

All of this encourages us to go beyond the present state, and provide a new
kind of review on PBHs.

###  1.3 Probing PBHs with LISA

This project is hosted by the cosmology working group of the LISA consortium.
Therefore our discussion has been focused on the GW signatures of PBHs in the
particular context of the LISA mission.

PBHs cover a very wide range of masses and therefore frequency ranges of GWs.
In particular, LISA will be sensitive to a broad band around the millihertz GW
frequencies [128, 159, 160, 161], which is complementary to the ground-based
GW detectors [93, 94] and electromagnetic probes [162]. The search for a SGWB
from PBH formation in the radiation era can be indirectly used for their
discovery or alternatively to constrain their existence in the
$10^{-12}~{}M_{\odot}$ mass range [163, 128, 159, 160, 161]. Moreover, the
coalescence of the heavy seeds of supermassive black holes (SMBH) in the late
Universe could leave clear signatures of their primordial origin, or even
detect individual events at high redshift which could not have arisen from
astrophysical mechanisms.

The great sensitivity of LISA at mHz frequencies [164, 165, 166, 167] opens
the possibility to detect the mergers of $10^{3}-10^{4}~{}M_{\odot}$ PBHs all
the way to $z\simeq 100$. Furthermore, the isotropic SGWB from the coalescence
of PBHs since recombination should have an amplitude and spectral shape that
will make it easily detectable by LISA. In this review, we will describe the
different features of known PBH scenarios that can be probed with LISA.

LISA will not only probe the GWs from individual PBH mergers and we will
review numerous other types of signatures of PBHs that can be probed with
LISA. Among them, the GW background from second order curvature fluctuations
will probe the very interesting asteroid-mass range. Other backgrounds include
the one from ultralight, evaporated PBHs, the ones from early and late PBH
binaries as well as from hyperbolic encounters. GW from individual sources
could also be observed by LISA, for instance from intermediate-mass PBH
mergers, including at high-redshifts, as above mentioned, or from inspirals
with extreme mass ratios which are expected for extended PBH mass
distributions. LISA will also have a key role if data are combined to LVK or
other detector observations at higher frequencies, for long-duration signals
including mergers with a subsolar component. Subsolar-mass PBHs could also
produce detectable quasi-monochromatic continuous waves in the LISA frequency
range. Finally, a more exotic possibility is to use LISA as a detector of
near-Earth asteroid-mass PBHs.

###  1.4 Outline of this review

This review is organized as follows. In Section 2, we present the different
classes of theoretical PBH formation scenarios, with a focus on recent
inflation models but including a discussion of other production mechanisms
such as (p)reheating, curvatons, phase transitions, topological defects and
primordial magnetic fields. Among inflationary scenarios, we especially focus
on recent proposals in which the usual Gaussian distribution of primordial
curvature fluctuations is modified in a more complex way than can be described
in terms of conventional non-Gaussian deviation parameters ($f_{\rm NL}$,
etc). Ultimately this could be one possible way to resolve the long-standing
fine-tuning issue for primordial fluctuations as the origin of PBHs.

Section 3 is dedicated to PBH formation and the computation of the PBH mass
distribution. We start with the standard formation formalism but then we
summarize and include in our calculations the results of the most recent and
accurate studies based on simulations in numerical relativity, such as non-
linear and non-Gaussian effects, the impact of the shape of the spectrum and
the curvature/overdensity profiles, the changes in the equation of state
linked to the thermal history, the mass and spin evolution after formation,
etc.

Once the PBH mass distribution is known, the next important step is the
calculation of PBH merger rates, covered in Section 4. PBH mergers are due to
two dominant channels: early binaries that are formed before matter-radiation
equality, and late binaries formed in PBH clusters. We use the latest
prescriptions, e.g., based on N-body simulations, in order to derive merger
rates for different mass scales and for some representative models, including
for very low mass ratios and intermediate-mass BH mergers. Finally, we
consider the case of hyperbolic encounters in PBH clusters that can also lead
to a large number of GW burst events.

In Section 5, we review the different sources of stochastic backgrounds
related to PBHs. One particularly relevant source of SGWB comes from the
curvature fluctuations that source GWs at second order in perturbation theory.
Other SGWBs are due to early or late PBH binaries as well as due to PBH
isocurvature fluctuations.

In Section 6, we review the possible ways to use the cross-correlations
between gravitational wave observations of individual mergers or the
stochastic background, and large-scale structure (LSS), in order to constrain
the existence and abundance of PBHs.

The current limits on the abundance of PBHs from various probes, including
GWs, are reviewed and discussed in Section 7. In particular, we aim to present
them in a model-dependent way. We emphasize and comment on the underlying
hypotheses and discuss how these limits change according to the different
classes of models.

Finally, the detectability of various GW signatures from PBHs with LISA is
discussed in Section 8, including individual merger events with extreme mass
ratios or intermediate-masses, hyperbolic encounters, stochastic backgrounds,
continuous waves, high-redshift signals, signal combination with other ground-
based detectors, etc.

At the end of each section, we provide an augmented discussion of recent
results, the current limitations and pave the way for future analyses by
listing some interesting perspectives. Our conclusions are presented in
Section 9.

##  2 Theoretical models

In this Section, we provide an overview of the principal mechanisms that can
lead to large curvature fluctuations and PBH formation. We start with
inflationary scenarios (see [168] for a recent review), either in single-field
(quasi-inflection point) or multi-field realisations (hybrid inflation, turns
in field space, gauge field interactions, etc.). Where possible, we provide
some general formulae for the calculation of the primordial power spectrum.
Then we consider recent models relying on quantum diffusion of the inflaton or
of a stochastic spectator field, generically leading to non-Gaussian
distributions of the curvature fluctuations. The third part of this Section is
dedicated to other models of PBH formation, due to curvaton fields,
preheating, phase transitions, early matter era, cosmic strings, domain walls
and primordial magnetic fields.

###  2.1 Single-field inflationary models

In the following, we provide a review of the main production mechanisms of
primordial black holes in the context of single-field inflationary models.

####  2.1.1 The basic idea and the slow-roll approximation

The inflationary epoch provides an ideal setting for producing PBHs in the
very early Universe, with single-field inflation being the minimal framework.
While a nearly scale-invariant spectrum of primordial curvature perturbations
is remarkably consistent with the observations of CMB temperature anisotropies
at large scales, the spectrum on smaller scales remains largely unconstrained.
When small-scale perturbations with sufficiently large overdensities re-enter
the expanding horizon during the radiation-dominated epoch, their gradients
induce a gravitational collapse that cannot be overcome even by the radiation
pressure of the expanding plasma, thereby producing PBHs with a mass of the
order of the horizon mass at the time of re-entry of a given wavemode, $k$.
The abundance of PBHs crucially depends on the collapse process and the nature
of primordial curvature perturbations, which are discussed later. If PBHs
start evaporating due to Hawking radiation right after their production, one
finds that PBHs with mass $m_{\rm PBH}\gtrsim 10^{-18}\,M_{\odot}$ will
survive until today and dominate the expansion of the Universe at matter-
radiation equality unless their initial abundance is strongly suppressed.

Since single-field slow-roll inflationary models typically predict a nearly
scale-invariant power spectrum of curvature perturbations on all scales, most
slow-roll single-field inflation models do not produce PBHs, because a
significant abundance of PBHs requires a substantial growth of the power
spectrum on small scales. Nevertheless, a large positive running of the scalar
primordial power spectrum can lead to PBH formation [9]. Another possibility
is to have a transition in the slope of the potential, with the scalar field
slowly rolling towards a region of almost constant potential, as one gets for
instance in the original hybrid inflation model along the valley. But the
difficulty is then to end inflation in a limited number of e-folds while, in
order to produce a large power, the potential is extremely flat. Another
possibility for producing PBHs in single-field inflation is through the
violation of the usual slow-roll conditions [10, 169, 170, 171, 172, 173],
although one should then pay a special attention to perturbative control
[173]. Moreover, not all single-field inflationary models can produce the
relevant abundance of PBHs as it crucially depends on the shape of the
potential and various model parameters. This typically leads to a high degree
of fine-tuning of model parameters.

The power spectrum of curvature perturbations in single-field models is given
at first order and in the slow-roll approximation, by

$\mathcal{P}_{\zeta}(k)=\frac{H_{*}^{2}}{\pi m_{\rm
Pl}^{2}\epsilon_{1*}}\left[1-2(C+2)\epsilon_{1*}+C\epsilon_{2*}-(2\epsilon_{1*}+\epsilon_{2*})\ln\left(\frac{k}{k_{*}}\right)\right]~{},$
(2.1)

where $C=0.578+2\ln 2-2$, $\epsilon_{1}\equiv-\dot{H}/H^{2}$ is the first
slow-roll Hubble-flow parameter, $\epsilon_{2}\equiv{\rm
d}\ln\epsilon_{1}/{\rm d}N$ is the second Hubble-flow parameter, $m_{\rm Pl}$
is the Planck mass (and $M_{\rm Pl}=\sqrt{8\pi}~{}m_{\rm Pl}$ denotes the
reduced Planck mass) and a star denotes a quantity evaluated at the time the
comoving wavenumber, $k_{*}$, crosses the Hubble radius during inflation when
$k_{*}=a_{*}H_{*}$. Therefore, one can naively expect that an enhancement of
$\mathcal{P}_{\zeta}$ can occur if there is a dynamical phase during the
inflationary expansion in which $\epsilon_{1}$ becomes much smaller than
unity. Such an epoch is often called an ultra slow-roll (USR) phase which can
be achieved by means of a plateau region in the inflationary potential, which
slows down the inflaton field more rapidly than in slow roll, before ending
inflation. An efficient and interesting way to produce this plateau is with
the introduction of an inflection point in the potential. However, one must
ensure that the inflaton does not stay too long at the inflection point
otherwise all the inflationary fluctuations that had successfully imprinted
the metric perturbations on large scales to explain the CMB anisotropies will
be expanded away. A schematic plot of the potential with an inflection point
is shown in Figure 1.

Figure 1: A schematic plot of the inflationary potential with an inflection
point. Adapted from Ref. [155].

####  2.1.2 Inflection point potentials

The simplest way to achieve an inflection or a near-inflection point in the
potential is to use a suitable polynomial potential or ratio of polynomials so
as to obtain the desired nearly scale invariant spectrum on CMB scales [169].
A very generic feature of such a potential is that although the inflaton
starts in a slow-roll regime, it dynamically enters into an USR stage close to
the inflection point for a short duration of e-folds that induces an
enhancement in the spectrum. Nevertheless, it is possible to evade the USR
regime in quasi-inflection point potentials [174]. However, one must always be
careful that this new feature in the potential does not affect the large-scale
dynamics which is already strongly constrained by CMB and large-scale
structure (LSS). These types of potentials arise naturally in models of Higgs
inflation [175, 174], MSSM inflation [176, 177], accidental inflation [178],
string inflation [30, 179], $\alpha$-attractor models [180, 181, 182],
punctuated inflation [183, 184], non-canonical inflation [185, 186, 187], as
an effective description of multifield inflation [188] and in more general
inflection point inflation [189]. However, to produce sufficient PBHs in these
models, one crucial difference is that, rather than imposing the initial
conditions on the inflaton evolution close to the inflection point, one must
assume that the field starts in the slow-roll regime at larger values of the
potential well above the inflection point, and then slow-rolls down towards
the minimum of the potential, crossing the inflection point, typically $30$ to
$40$ e-folds before the end of inflation. Such a dynamical evolution must
ensure that large-scale dynamics are consistent with the CMB constraints, PBHs
are produced at some intermediate smaller scales during the USR regime and
inflation naturally ends when the field finally rolls down towards the true
minimum of the potential. In this context, note that there has been an intense
discussion recently on loop corrections induced on CMB scales by enhanced
scalar perturbations on PBH scales [190, 191, 192, 193, 194, 195, 196, 197,
198]. Hints for significant one-loop contribution, to the point of undermining
the consistency of such set-ups, have appeared in computations including a
subset of vertices as derived from the interaction Hamiltonian. However, it
has been recently shown that these effects disappear due to cancellations once
all relevant interaction terms, at least at cubic order, are taken into
account [199], and this subject will likely deserve further studies. Note also
that enhanced fluctuations can truly lead, through loop effects, to an
increase of power on (intermediate) larger scales due to infrared rescattering
[200], in single-field models and beyond, with potentially important
consequences for the PBH population.

####  2.1.3 Example model: critical Higgs inflation

Critical Higgs Inflation (CHI) is probably among the best physically-motivated
single-field models of inflation which can give rise to a quasi-inflection
point [174] and a high plateau in the matter power spectrum. The non-minimally
coupled Higgs action is given by

$S\\!=\\!\\!\int\\!{\rm d}^{4}x\sqrt{g}\left[\left(\frac{M_{\rm
Pl}^{2}}{2}+\frac{\xi(\phi)}{2}\phi^{2}\right)\\!R-\frac{1}{2}(\partial\phi)^{2}-\frac{1}{4}\lambda(\phi)\phi^{4}\right]$
(2.2)

with the running of the couplings is parametrized by

$\displaystyle\lambda(\phi)$ $\displaystyle=$
$\displaystyle\lambda_{0}+b_{\lambda}\ln^{2}\left(\phi/\mu\right)\,,$ (2.3)
$\displaystyle\xi(\phi)$ $\displaystyle=$
$\displaystyle\xi_{0}+b_{\xi}\ln\left(\phi/\mu\right)\,,$ (2.4)

around the critical point $\phi=\mu$.

After standard metric and scalar field redefinitions,

$\displaystyle g_{\mu\nu}$ $\displaystyle\to$
$\displaystyle\left(1+M_{\rm{Pl}}^{-2}\xi(\phi)\phi^{2}\right)\,g_{\mu\nu}\,,$
(2.5) $\displaystyle\phi$ $\displaystyle\to$ $\displaystyle\varphi=\int{\rm
d}\phi\frac{\sqrt{1+M_{\rm{Pl}}^{-2}{\xi(\phi)\phi^{2}}\left(1+6\left(\xi(\phi)+\frac{1}{2}{\phi\xi^{\prime}(\phi)}\right)^{2}\right)}}{1+M_{\rm{Pl}}^{-2}{\xi(\phi)\phi^{2}}}\,,$
(2.6)

the effective inflationary potential becomes

$V(x)=\frac{V_{0}\,(1+a\,\ln^{2}x)\,x^{4}}{(1+c\,(1+b\,\ln x)\,x^{2})^{2}}\,,$
(2.7)

with $V_{0}=\lambda_{0}\mu^{4}/4$, $a=b_{\lambda}/\lambda_{0}$,
$b=b_{\xi}/\xi_{0}$ and $c=\xi_{0}\,M_{\rm{Pl}}^{-2}\mu^{2}$. The potential
has a flat plateau at large values of the field $x=\phi/\mu$, see Fig. 1,
where $V(x\gg x_{\rm c})\simeq
V_{0}\,\frac{a}{(b\,c)^{2}}=\frac{M_{\rm{Pl}}^{4}}{4}\,\frac{b_{\lambda}}{b_{\xi}^{2}}\ll
M_{\rm Pl}^{4}\,.$ The potential also has a short secondary plateau around the
critical point, $\phi_{\rm c}=\mu$, where the inflaton-Higgs slows down and
induces a large peak in the curvature power spectrum. This second plateau is
induced by a quasi-inflection point at $x=x_{\rm c}$, where $V^{\prime}(x_{\rm
c})\simeq 0,\,V^{\prime\prime}(x_{\rm c})\simeq 0$. As a consequence, the
number of $e$-folds has a sharp jump, $\Delta N$, at that point, plus a slow
rise towards larger field values, corresponding to CMB scales. This behavior
is very similar to the one discussed in Ref. [169].

In order to have large PBH production and good agreement with the CMB
constraints, the allowed range of CHI parameters can be found to be
$\lambda_{0}\approx(0.01–8)\times 10^{-7}$, $\xi_{0}\approx(0.5-15)$,
$x^{2}\mu^{2}\approx(0.05–1.2)$, $b_{\lambda}\approx(0.008–4)\times 10^{-6}$
and $b_{\xi}\approx(1–18)$, for the USR phase to occur $N\in(30,35)$ e-folds
before the end of inflation. The question arises whether these values,
corresponding to the model parameters at the critical scale $\mu$, are
consistent with the values of the SM parameters at the EW scale. Given the
latest values of the top quark mass $m_{\rm top}$ and strong coupling
$\alpha_{s}$ [201, 202], the values of $\lambda_{0}$ and $b_{\lambda}$ that we
consider for the Higgs quartic coupling, are consistent, within $2\sigma$,
with the existence of a critical point.

Figure 2: Power spectrum $\mathcal{P}_{\mathcal{R}=\zeta}(k)$ for the single-
field Critical Higgs Inflation model, with an inflection point at $N\approx
36$ satisfying the Planck 2018 constraints. Adapted from Ref. [169].

Taking into account the running of both the Higgs self-coupling $\lambda$ and
the non-minimal coupling to gravity $\xi$, there are regions of parameter
space allowed by the Standard Model for which the inflaton-Higgs potential
acquires a second plateau at smaller scales, around the critical point
$\lambda(\mu)\simeq\beta_{\lambda}(\mu)=0$. This plateau gives a super-slow-
roll evolution of the Higgs, inducing a high peak in the curvature power
spectrum which is very broad. When those fluctuations reenter the horizon
during the radiation era, they collapse to form PBHs with masses in the range
$0.1$ to $100\,M_{\odot}$.

In single-field models with a near-inflection point, the inflaton slows down
right before the end of inflation, creating a stronger backreaction and a
quick growth in curvature fluctuations, giving rise to a significant increase
in the power spectrum at scales much smaller than those probed by CMB and LSS
observations. Since the spectral amplitude is essentially inversely
proportional to the parameter $\epsilon_{1}$, whose exact evolution can
deviate significantly from the slow-roll approximation [169], the power
spectrum has to be calculated by integrating numerically the evolution of the
inflaton field $\phi$ following the exact (beyond slow-roll) equations. In
Fig. 2, we show a typical $\mathcal{P}(k)$ produced in quasi-inflection-point
models like CHI [174]. One can notice the difference between the exact power
spectrum and that obtained in the slow-roll approximation. The power spectrum
can be parametrized as a double step with different amplitudes and different
tilts:

${\mathcal{P}}_{\zeta}(k)=\left\\{\begin{array}[]{cc}A_{1}\left(\frac{k}{k_{1}}\right)^{n_{{\rm
s}1}-1}&{\rm for}\ \ k<k_{\rm c}\\\\[5.69054pt]
A_{2}\left(\frac{k}{k_{2}}\right)^{n_{{\rm s}2}-1}&{\rm for}\ \ k>k_{\rm
c}\end{array}\right.$ (2.8)

where $A_{i}$, $n_{{\rm s}i}$, $k_{i}$ are two different amplitudes, spectral
tilts and arbitrary pivot scales for the power spectrum, before and after the
cut at $k_{\rm c}$. It is also possible to introduce some small amount of
running of the tilt. This power spectrum neglects quantum diffusion, which
induces non-Gaussian contributions [203, 204] and exponential tails [143],
which can be responsible for a large probability of PBH formation even with a
relatively small power spectrum amplitude, as will be discussed in subsequent
sections.

Figure 3: Examples of primordial power spectra of curvature fluctuations,
leading to (stellar-mass) PBH formation: power-law, broken power-law, Gaussian
and log-normal models (solid lines), and particular examples of multifield and
axion-gauge models (dashed lines).

####  2.1.4 Reverse engineering approach

In single-field inflationary models, the commonly adopted strategy is to start
from some scalar potential that features an approximate stationary inflection
point, as shown in Fig. 1; the latter is controlled by various free
parameters, which are fine-tuned to guarantee the desired enhancement of the
power spectrum of curvature perturbations [205, 169, 206, 207, 208, 180, 209,
210, 211, 181, 28]. Following Ref. [212], a different, and much more powerful,
approach can be entertained. In a reverse engineering procedure, the model
building can be moved away from the scalar potential—that from this
perspective is the output of the analysis instead of the starting point—and
focused on the underlying inflationary dynamics (see refs. [213, 214, 215,
216] for a similar viewpoint).

The inflationary background dynamics can be captured by specifying a minimal
set of inputs. Consider the Hubble parameters $\epsilon_{1}$ and
$\eta\equiv\epsilon_{1}-\epsilon_{2}/2$, and the number of e-folds $N$,
defined by $\mathrm{d}N=H\mathrm{d}t$. Ref. [212] devised the following
analytical ansatz for the time-evolution of $\eta$

$\displaystyle 2\eta(N)=$ $\displaystyle\left[\eta_{\rm I}-\eta_{\rm
II}+(\eta_{\rm II}-\eta_{\rm I})\tanh\left(\frac{N-N_{\rm I}}{\delta N_{\rm
I}}\right)\right]+\left[\eta_{\rm II}+\eta_{\rm III}+(\eta_{\rm III}-\eta_{\rm
II})\tanh\left(\frac{N-N_{\rm II}}{\delta N_{\rm II}}\right)\right]+$
$\displaystyle\left[\eta_{\rm IV}-\eta_{\rm III}+(\eta_{\rm IV}-\eta_{\rm
III})\tanh\left(\frac{N-N_{\rm III}}{\delta N_{\rm III}}\right)\right]\,.$
(2.9)

The inflationary dynamics is therefore divided in different stages: i) a first
stage where $\eta_{\rm I}$ is constant and negative (necessarily taken to be
small to reproduce the conventional slow-roll dynamics) so that $\epsilon_{1}$
increases relatively gently; ii) a USR phase characterized by $\eta_{\rm
II}>(3+\epsilon_{1})/2\simeq 3/2$ , where negative friction makes the
parameter $\epsilon_{1}$ decrease abruptly iii) an optional intermediate stage
where $\epsilon_{1}$ remains constant, achieved if $\eta_{\rm III}\simeq 0$,
thus generating a second slow roll phase (i.e. an enhanced plateau in the
curvature power spectrum [131, 212]); iv) a final phase characterized by
$\eta_{\rm IV}<0$, that brings $\epsilon_{1}$ back to $O(1)$ values for
inflation to end. The parameters $\delta N_{\rm I,II,III}$ control the
sharpness of the transitions.

As shown in Refs. [212, 87], it is possible to choose the parameters
controlling the ansatz of Eq. (2.9) in such a way that the power spectrum of
curvature perturbations is compatible with CMB constraints at large scales as
well as featuring an enhanced peak or plateau at scales beyond the FIRAS
constraints (i.e. for $k\gtrsim 10^{4}/{\rm Mpc}$), where a sizeable PBH
abundance can be obtained, corresponding to the possible formation of masses
smaller than ${\cal O}(10^{4})M_{\odot}$.

Once the Hubble parameters are known, one can compute the inflationary
potential by means of

$\displaystyle\frac{V(N)}{V(N_{\rm ref})}$
$\displaystyle=\exp\left\\{-2\int_{N_{\rm
ref}}^{N}dN^{\prime}\left[\frac{\epsilon_{1}(3-\eta)}{3-\epsilon_{1}}\right]\right\\}\,,~{}~{}~{}~{}~{}~{}~{}\phi(N)=\phi(N_{\rm
ref})\pm\int_{N_{\rm ref}}^{N}dN^{\prime}\sqrt{2\epsilon_{1}}\,,$ (2.10)

where $N_{\rm ref}$ denotes an arbitrary reference e-fold time, and where in
the second equation we consider the minus sign having in mind a large-field
model in which the field value decreases as inflation proceeds. Combining
$V(N)$ and $\phi(N)$, we reconstruct the profile $V(\phi)$ of the inflationary
potential in field space [217]. Eq. (2.10) shows the convenience of modelling
the inflationary dynamics directly at the level of $\eta$ instead of
$V(\phi)$. This is because the Hubble parameters enter in the exponent in the
definition of $V(N)$, and thus allow for a much finer control on power
spectral features when performing the reverse engineering procedure. This
approach facilitates model building, even though it does not provide an
interpretation of the reconstructed potential in terms of high-energy theories
of inflation.

###  2.2 Multi-field inflationary models

As seen above, a common challenge of PBH formation within single-field
inflationary models is keeping control over the two phases that lead
respectively to the seeds for the CMB at large cosmic scales and for PBHs at
smaller scales.

Multi-field inflationary scenarios can help separate and better control these
two phases of perturbation generation, while also being very natural
candidates from the UV perspective. Furthermore, violation of the standard
slow-roll conditions can be efficiently realized in a context with more
degrees of freedom [105].

####  2.2.1 Hybrid inflation

In a setup with two scalar fields, the most studied scenario is perhaps hybrid
inflation [218, 219]. In this class of models, inflation takes place along an
almost flat valley in the direction of an initial slowly rolling scalar field,
$\phi$. Inflation ends with a waterfall instability when $\phi$ crosses a
critical point, $\phi_{\rm c}$, due to the presence of a second field, $\psi$.
The original hybrid potential can be written as

$V(\phi,\psi)=\Lambda^{4}\left[\left(1-\frac{\psi^{2}}{M^{2}}+\frac{\phi^{2}}{\mu^{2}}+2\frac{\phi^{2}\psi^{2}}{\phi_{\rm
c}^{2}M^{2}}\right)\right].$ (2.11)

The general shape of the hybrid potential is shown in Fig. 4.

Figure 4: Representation of a typical hybrid inflation potential, $V$, with a
possible trajectory in two-field space (dotted line). CMB perturbations are
created along the valley at $\psi=0$, during a first phase of inflation.
Curvature perturbations suitable for production of PBHs are generated in a
second flat part of the potential (red dotted line), when the mass-squared of
the field $\psi$ changes sign.

When $\mu M\ll M_{\rm Pl}^{2}$, the waterfall phase transition is nearly
instantaneous and lasts less than one e-fold. The global minima of the
potential are located at $\phi=0$ and $\psi=\pm M$. Inflation can take place
along the valley $\psi=0$ with the effective potential

$V(\phi,\psi)=\Lambda^{4}\left(1+\frac{\phi^{2}}{\mu^{2}}\right)$ (2.12)

that most often produces a blue-tilted power spectrum for the original
potential (a red-tilted spectrum can nevertheless be obtained for some
specific parameter values [105]). Other hybrid realisations, like F-term or
D-term inflation, can produce other shapes for the effective potential,
leading to a red-tilted spectrum, compatible with the Planck observations (see
e.g. [220], which realises this scenario via hybrid $\alpha$-attractors models
with scalar potentials with no tachyonic directions and, therefore, bounded
from below). PBHs of tiny mass can be produced during a fast waterfall phase
[11], which would have evaporated today, eventually leading to the reheating
of the Universe.

In the other regime $\mu M\gg M_{\rm Pl}^{2}$, the waterfall phase is slow,
lasting $N\gg 60$ e-folds [221], and it is then crucial to consider the exact
two-field dynamics to calculate the shape of the primordial power spectrum,
which is however not compatible with Planck CMB measurements [222].

Hybrid inflation becomes interesting for the production of massive PBHs in the
transient regime, when $\mu M\sim M_{\rm Pl}^{2}$, producing a significant but
not too large number of e-folds, $10\lesssim N\lesssim 40$. During the mild
waterfall phase, the multi-field perturbation dynamics leads to a strong,
exponential power spectrum enhancement (that can be sharper than the typical
maximal $\propto k^{4}$ increase in single-field models [217]) followed by a
relatively sharp decrease, when going from large to small scales. The position
and amplitude of the resulting broad peak is driven by a single combination of
the hybrid potential parameters, $\Pi\equiv M\sqrt{\phi_{\rm c}\mu_{1}}/M_{\rm
Pl}^{2}$ where $\mu_{1}$ comes from a linear expansion of the potential around
the critical point $V\propto(\phi-\phi_{\rm c})/\mu_{1}$. Numerical
calculations of the primordial power spectrum of curvature fluctuations using
either the exact multi-field perturbation dynamics or the $\delta N$
formalism, were first done in Ref. [105]. The resulting power spectrum is well
approximated by a log-normal in mass spectrum,

$\mathcal{P}_{\zeta}(k)\approx P_{\rm p}\exp\left[-\frac{(N_{k}-N_{\rm
p})^{2}}{2\sigma_{\rm p}^{2}}\right]\,,$ (2.13)

of maximal amplitude $P_{\rm p}$ and width $\sigma_{\rm p}$. Some examples of
power spectra are shown in Fig. 5. Finally, let us emphasize that the exact
amplitude of the power spectrum also depends on the quantum stochastic
dynamics of the auxiliary field $\psi$ close to the critical instability
point. In principle, this should lead to slightly different power spectrum
amplitudes in different regions of the Universe, and therefore different PBH
abundances. A natural outcome of the hybrid inflation scenario would thus be
that PBHs are formed in clusters whose mass and size depends on these quantum
stochastic dynamics. This can be also somehow used to reduce the fine-tuning
linked to the required power spectrum amplitude, since many realisations would
have been generically produced within the Universe.

Figure 5: Power spectrum of curvature perturbations for hybrid inflation
parameters values $M=0.1M_{\rm{Pl}}$, $\mu_{1}=3\times 10^{5}M_{\rm{Pl}}$ and
$\phi_{\rm c}=0.125M_{\rm{Pl}}$ (red)$,0.1M_{\rm{Pl}}$ (blue),
$0.075M_{\rm{Pl}}$ (green), and $0.05M_{\rm{Pl}}$ (cyan). Those parameters
correspond respectively to $\Pi^{2}=375/300/225/150$. The power spectrum is
degenerate for lower values of $M,\phi$ and larger values of $\mu_{1}$,
keeping the combination $\Pi^{2}$ constant. For larger values of $M,\phi_{\rm
c}$ the degeneracy is broken: power spectra in orange and brown are obtained
respectively for $M=\phi_{\rm c}=M_{\rm{Pl}}$ and
$\mu_{1}=300M_{\rm{Pl}}/225M_{\rm{Pl}}$. Dashed lines assume
$\psi_{c}=\psi_{0}$ whereas solid lines are obtained after averaging over 200
power spectra obtained from initial conditions on $\psi_{c}$ distributed
according to a Gaussian of width $\psi_{0}$. The power spectra corresponding
to these realizations are plotted in dashed light gray for illustration. The
$\Lambda$ parameter has been fixed so that the spectrum amplitude on CMB
anisotropy scales is in agreement with Planck data. The parameter
$\mu_{2}=10M_{\rm{Pl}}$ so that the scalar spectral index on those scales is
given by $n_{\rm s}=0.96$. Figure from Ref. [105].

####  2.2.2 Turns in multi-field inflation

Embeddings of inflation in high energy physics motivate scenarios in which
multiple scalar fields participate in the inflationary dynamics. An ubiquitous
phenomenon in this context is the one of turns in field space, corresponding
to a bending of the inflationary trajectory, or more accurately to its
deviation from a geodesic of the field space,333In the context of string
theory compactifications, it has recently been shown [223] that deviations
from geodesic trajectories must be negligible near the boundary of the moduli
space. Large turns can therefore be engineered just in the moduli space bulk.
quantified by a dimensionless parameter $\eta_{\perp}(t)$ measuring the
acceleration of the trajectory perpendicular to its direction [224, 225]. When
such a turn is sufficiently strong, with $\eta_{\perp}\gtrsim 1$ (see [226]
for a short overview), it leads to a transient tachyonic instability of the
scalar field fluctuations before Hubble crossing [227, 228, 229], similar in
essence to the one arising in axion gauge-field inflation (see section 2.2.4).
This gives rise to an exponential enhancement of the power spectrum on those
scales, compared to the ones that cross the Hubble radius before and after the
turn, and thus constitutes a natural mechanism to seed PBH formation [230,
231, 232, 181]. Interestingly, owing to its inherently multi-field origin, the
growth rate of the power spectrum can overcome the bound deduced in single-
field setups [217, 233, 234], i.e., the peak of the power spectrum can be
sharper. An important parameter is the duration of the turn in e-folds
$\delta\equiv\delta N$, which leads to qualitatively different behaviours for
broad ($\delta\gtrsim 1$) and sharp turns ($\delta\ll 1$). While broad turns
lead to a featureless bump of the power spectrum, sharp turns lead to a
localised peak of the envelope of the power spectrum, modulated by order-one
rapid oscillations, as described in simple cases by the analytical power
spectrum (assuming a constant value of $\eta_{\perp}$ during the turn) [230,
235]

$\frac{\mathcal{P}_{\zeta}(k)}{{\cal
P}_{0}}=\,\frac{e^{2\sqrt{(2-\kappa)\kappa}\,\eta_{\perp}\delta}}{2(2-\kappa)\kappa}\times\sin^{2}\left[e^{-\delta/2}\kappa\eta_{\perp}+\arctan\left(\frac{\kappa}{\sqrt{(2-\kappa)\kappa}}\right)\right]\,.$
(2.14)

Here, ${\cal P}_{0}$ denotes the amplitude of the power spectrum in the
absence of the transient instability, $\kappa=k/k_{\textrm{f}}$ is the
dimensionless wavenumber in units of the maximally enhanced one, with
$k_{\textrm{f}}$ the scale crossing the Hubble radius at the time of the sharp
turn, and the formula is valid for $\kappa\leq 2$ (see [235] for
generalisations). The oscillations there, which are to a very good
approximation periodic in $k$, are characteristic of a sharp feature, and go
hand in hand with a boosted power spectrum generated by an event localised in
time [235]. The mass function of PBHs is highly sensitive to the tail of the
probability density function (PDF) of the smoothed density contrast, which is
not yet known for the class of models described here. However, tentative
results for the mass function indicate that the oscillatory patterns are
washed out when assuming Gaussian statistics of the primordial curvature
fluctuations [231], and hence further investigation is required to assess the
observational consequences for PBHs of the oscillatory features. What is
known, however, is that the SGWB generated by scalar fluctuations at horizon
re-entry does keep an imprint of the oscillations of the scalar power
spectrum: they manifest as oscillations in the frequency profile of the SGWB,
providing a robust probe of small-scale features that would be inaccessible
otherwise [235, 236, 237, 238, 239].

If one uses the standard criterion that the power spectrum has to reach values
${\cal P}_{\zeta}\sim 10^{-2}$ to generate a substantial amount of PBHs, it is
interesting that the required enhancement from the baseline value ${\cal
P}_{0}\sim 10^{-9}$ necessitates a total angle
$\Delta\theta=\eta_{\perp}\delta$ swept by the trajectory during the turn
exceeding the value $\pi$ obtainable in flat field space, and hence requires
the inflationary field space to be curved [230]. Such large values of the
power spectrum are not without consequences however, and considerations of
backreaction and perturbative control are important to take into account to
delineate the space of theoretically viable models [235]. Indeed, the
criterion for perturbative control parametrically depends on the large
parameter $\eta_{\perp}$ in a way that imposes the bound
$\eta_{\perp}^{6}{\cal P}_{\zeta}\lesssim 1$ (or $\eta_{\perp}^{4}{\cal
P}_{\zeta}\lesssim 1$) for broad [240, 241] (or sharp) turns [235]. Thus, an
important open question is whether the inevitably present non-Gaussianities of
the curvature fluctuations in these models, enhanced in flattened
configurations [242, 240, 243], are such that a substantial amount of PBHs can
be generated with a value of the power spectrum less than $10^{-2}$ while
ensuring perturbative control.

####  2.2.3 Two-stage models

There are several concrete examples of PBH formation resulting from a turn in
the field space during inflation that separates two stages of inflation. The
first stage produces the primordial power spectrum on CMB scales while the
second stage is responsible for the enhancement of power on smaller scales,
leading to PBH formation. Oscillations are expected at the transition, as
explained in the previous Section. Pi et al. [244] have proposed a model
combining Starobinsky inflation with a scalaron field that becomes massive at
the end of the first stage of inflation and a curvature field that plays the
role of the inflaton in a second stage. Kawasaki et al [245] proposed a
scenario with a second phase of double inflation in a supergravity framework.
One should notice that non-Gaussianity is typically small in these models,
since most of the field trajectories are effectively single-field. Other
advantages of these models are that one does not especially need to tune the
position of a peak in the power spectrum, since it can be nearly scale
invariant, as generically predicted by slow-roll inflation, and that the
transition between the two regimes can be very sharp. Due to the features at
the transition in the power spectrum, such models can also produce multi-modal
distributions of black holes, as proposed in [246].

####  2.2.4 Axion-gauge scenario

PBHs can also arise in models in which rolling axions interact with gauge
fields, enhancing primordial density perturbations [247, 248, 249, 250, 128,
251, 252]. In natural inflation [253, 254], the couplings of the inflaton to
matter respect a shift symmetry, $\phi\to\phi+{\rm constant}$. Therefore, they
do not provide by themselves quantum corrections to the inflaton potential
(that arise from a breaking of the shift symmetry, as for instance from
instanton effects), which can therefore be technically naturally small and
flat. Also, UV-completed theories, such as string theory or supergravity
models of inflation, contain a large number of pseudo-scalar particles that
could be identified with the axion inflaton, or some other dynamically
relevant axion. The shift-symmetric coupling between the axion and the gauge
field is

${\cal L}\supset\frac{a}{f}F{\tilde{F}},$ (2.15)

where $a$ is the pseudoscalar axion, $F$ is the gauge field strength,
${\tilde{F}}$ its dual, and $f$ is the axion decay constant. Shift-symmetry is
satisfied due to the fact that $F{\tilde{F}}$ is a total derivative and the
constant term becomes a topological invariant. This interaction amplifies one
helicity of the vector field around horizon crossing. The amplified gauge
quanta might affect the evolution of the axion field [255] or, even if
produced in a smaller amount, source scalar density perturbations [256] and
parity violating tensor modes [257]. These effects are exponentially sensitive
to the parameter $\xi\equiv\frac{\dot{a}}{2fH}$, and visible effects are
obtained for $f\simeq 10^{-2}M_{\rm Pl}$, namely for an axion decay constant
parametrically equal to the string scale. On one hand, this enhances the
possibility that models of string cosmology have a rich phenomenology. On the
other hand, no effect is expected for models of a single axion inflaton, that
are required to have a super-Planckian excursion, unless some specific
mechanism allows excursions much greater than the scale of inflation, such is
the case for monodromy [258] or for aligned axion inflation [259].

Scalar perturbations are produced via inverse decay $A+A\to a$ (where $A$
denotes the gauge field), and therefore, they obey a non-Gaussian $\chi^{2}$
distribution. The amplified gauge fields also produce GWs via an analogous
$2\to 1$ process. We note that these primordial GWs, which are produced during
inflation, form a different stochastic GW background with respect to the
secondary GWs that are sourced by density perturbations in the process of
horizon re-entry in radiation or matter-dominated eras.

The power spectrum of density perturbations can be parametrized as [260, 261,
250]

$\mathcal{P}_{\zeta}(k)=\mathcal{P}_{\rm vacuum}(k)+\mathcal{P}_{\rm
sourced}(k)\;,$ (2.16)

where $\mathcal{P}_{\rm vacuum}\simeq 2\times 10^{-9}\,\left(k/k_{\rm
cmb}\right)^{n_{\rm s}-1}$ is the standard nearly scale invariant vacuum
signal (with $k_{\rm cmb}$ being some pivot wavenumber at CMB scales). As
already mentioned, the sourced signal is exponentially sensitive to the
parameter $\xi\equiv\dot{a}/{2fH}$, and it is therefore highly sensitive to
the speed of the axion, and, ultimately, to its potential. For axion inflation
with a single monomial potential the signal is typically blue, and one simply
needs to require that the production is sufficiently small so that no stable
PBH is produced via this mechanism [247]. This does not need to be the case
for a more complicated potential, or in the case in which the axion is not
identified with the inflaton [250, 128]. In this second case, one can imagine
a situation in which the axion rolls significantly only for a limited number
of e-folds [260, 261] and the sourced spectrum can be described as a nearly
log-normal bump, peaked at the scales that were exiting the horizon while $a$
was rolling. The infra-red part of the spectrum scales as $\mathcal{P}\propto
k^{3}$ due to causality, while, close to the peak, the spectrum is well
parametrized by [260]

$\mathcal{P}_{\rm sourced}(k)={\cal A}(\xi_{*},{\cal W})\,{\rm
exp}\left[-\frac{\ln^{2}(k/k_{*})}{2\sigma^{2}(\xi_{*},{\cal W})}\right]\;,$
(2.17)

where the amplitude ${\cal A}$, the central position $k_{*}$, and the width
$\sigma$ of the peak are functions of the maximum value of the particle
production parameter $\xi$ and of the amount of e-folds ${\cal W}$ during
which the axion is significantly rolling. An analogous parametrization also
holds for the amount of sourced GWs [260].

Figure 6: The stochastic GWs produced by axion inflation at four main scales
of interferometers for which we show the approximate sensitivities, nHz (PTAs
and SKA), $\mu$Hz (Gaia and Theia), mHz (LISA and BBO) and Hz (LVK, Einstein
Telescope, and Cosmic Explorer). The enhanced density perturbations may
produce PBHs which are a significant fraction of dark matter for
$1-100M_{\odot}$ and the totality in the $10^{-14}-10^{-11}M_{\odot}$ mass
range and, remarkably, the enhanced perturbations leave inevitable GW
backgrounds at the most sensitive regimes of GW detectors: for the first range
of BHs this corresponds to PTA frequencies (and possibly future SKA) and LISA
frequencies for the second range. The PBHs formed from fluctuations peaking at
Hz scales will be so light that they are expected to be part of thermal
history through Hawking radiation. The Figure is updated (with slight
modifications) from Refs. [128, 161, 262]

The sourced scalar perturbations are highly non-Gaussian [256]. For a nearly
scale-invariant signal (as is the case if the axion is the inflaton), their
bispectrum has a shape which is very close to the equilateral one. For the
case of short-duration production leading to Eq. (2.17), the bispectrum is a
maximum when the magnitude of the three momenta is parametrically equal to
$k_{*}$, and its scale-dependence can also be described accurately via a log-
normal shape [260]. The fact that the sourced scalar perturbations obey a
$\chi^{2}$ statistics makes PBH production much more significant compared to
Gaussian perturbations with the same power. Therefore, an equal amount of PBHs
is produced with a much smaller power spectrum than in the Gaussian case. This
results in a significantly smaller amount of stochastic GWs produced at second
order by the scalar perturbations, with a different prospect for detection at
LISA and PTA-SKA, see for example Fig. 6. This is discussed in more detail in
Subsection 5.1.

####  2.2.5 False-vacuum models

It is also possible to form PBHs in both single and multi-field scenarios if
some spatial regions end up in a false-vacuum. These regions can collapse into
black holes before the field can pass the potential barrier through quantum
tunneling and evolve towards the true vacuum. Several examples have been
proposed.

For instance, the Higgs field is suspected to be metastable. If it is a
spectator field during inflation, some regions may have ended in this region
of the potential, thereby leading to the formation of black holes [263].
Nevertheless, their abundance cannot be significant enough to contribute to
the dark matter. Another possibility is a tiny step in the inflaton potential
that can basically be overcome by most of field trajectories, but not all
[264]. This can radically change the tail of the Gaussian fluctuations, even
if the conventional non-linearity parameters ($f_{\rm NL},g_{\rm NL}$,…)
remain small. The regions with large curvature fluctuations later form PBHs.

It has also been proposed that a light spectator field during inflation,
having an asymmetric polynomial potential, experiences stochastic fluctuations
that displace the field from the global minimum of its potential and populate
a false vacuum state [265]. Lattice simulations have been used to show how
when the Hubble radius reaches the false-vacuum bubble size, the potential
energy is transferred to the kinetic and gradient energies of the bubble wall,
such that it begins to contract and possibly collapse into a black hole.
Finally, multi-field realisations have been proposed, for instance in [266,
267].

###  2.3 Non-Gaussian models

####  2.3.1 Motivations

In Subsection 2.1 we discussed how PBHs could be generated in single-field
models due to an enhancement of the power spectrum caused, for example, by a
second plateau in the potential. However, as we will discuss in more detail in
Subsection 3, what really matters to determine the PBH formation is the full
profile of the PDF of the primordial curvature perturbations, in particular
the probability to be above a given threshold. This implies that PBH abundance
is largely controlled by the tails of the PDF, which at the same time are
highly sensitive to any non-Gaussian corrections.

Non-Gaussian corrections can either enhance or suppress the production of PBHs
depending on their effect on the tail of the PDF. One way of accounting for
these effects is to compute the different statistical moments of the
distribution. Gaussian profiles are fully characterized by the second moment
(the power spectrum), while general PDFs can have contributions from any of
the higher moments. For illustrative purposes, we plot in Fig. 7 the effect of
the third and fourth moment as determined by the skewness and kurtosis,
respectively. As it can be clearly seen, a positive kurtosis has a large
impact in the tail of the PDF. Interestingly, one can analytically compute the
effect of non-zero higher order moments in the abundance of PBHs [268, 269].

Non-Gaussian corrections are known to be relevant in multi-field models of
inflation. For instance, models of axion inflation in which the gauge field
sources the curvature perturbations [250] display $\chi^{2}$ statistics for
the curvature perturbation. This has implications both for the PBH abundance
and for the GW signature, as studied in detail in [128]. In particular, a
strong increase of the PBH abundance associated with a greater tail of the PDF
does not correspond to a strong increase of the associated GW production.
Therefore, for equal PBH abundance, a smaller GW signal is produced with
respect to the Gaussian case.

Moreover, these non-Gaussian corrections can also enhance the abundance of
PBHs in single-field models due to the dynamics of the ultra-slow roll phase
[270, 271] and the effect of quantum diffusion beyond slow-roll [272, 203,
204, 273]444A violation of slow-roll is necessary for non-Gaussian effects to
be sufficiently large to change the PBH abundance in single-field models
[274].. In fact, in general terms, quantum diffusion introduces an exponential
tail in the PDF of curvature perturbations [275, 143]. For this reason, we
discuss next the effect of quantum diffusion in the production of PBHs.

Figure 7: Illustration of the effect of a positive skewness,
$\bar{\kappa}_{3}$, and a positive (excess) kurtosis, $\bar{\kappa}_{4}$, for
the probability density function $P(x)$ (left) and its logarithm (right) in
comparison with a Gaussian distribution ($\bar{\kappa}_{n>2}=0$). The dominant
effect in the tail is given by a positive kurtosis. Plots reproduced from
[204].

####  2.3.2 Quantum diffusion

Figure 8: Probability density functions of the curvature perturbations
generated by an inflection-point inflationary potential (see inset). Solid
lines correspond to the full distribution functions computed by mean of the
stochastic-$\delta N$ formalism, where different colours correspond to
different locations in the potential $v(\phi)$ where the scale under
consideration emerges from the Hubble radius and where $\Delta\phi_{\rm well}$
denotes the field range for which stochastic effects dominate over the
classical dynamics. The dotted lines correspond to the standard result, which
provides a good Gaussian approximation for the maximum of the distribution but
that however fails to describe the exponential tails, where PBHs are
nonetheless produced. This figure is adapted from Ref. [143].

Primordial curvature perturbations are expected to be well-described by a
quasi-Gaussian distribution when they are small and close to the maximum of
their probability distribution. This can be modelled using conventional
cosmological perturbation theory where the free field fluctuations during
inflation describe a Gaussian distribution and higher-order interactions lead
to non-Gaussian corrections. However, PBHs result from the collapse of large
density perturbations, which may be far from the peak of the distribution. In
this regime, we require a non-perturbative approach, such as the
stochastic-$\delta N$ formalism [276, 277, 278].

Stochastic inflation [279] splits the inflaton field into small-scale quantum
fluctuations and a large-scale effectively classical field above a coarse-
graining scale (typically around the Hubble length). The large-scale field
follows the non-linear evolution of a (locally) homogeneous and isotropic
cosmology [280, 281, 282], except that small-scale modes are constantly swept
up from the quantum vacuum state, crossing into the coarse-grained field,
leading to quantum diffusion
$\langle\Delta\varphi^{2}\rangle\approx(H/2\pi)^{2}$ per Hubble volume per
Hubble time. In regimes where PBHs are formed, the inflaton must undergo large
quantum fluctuations (for instance as the result of travelling through a very
flat potential), which implies that stochastic corrections are important and
must be taken into account [275, 204, 203, 143].

In order to derive the predicted abundance of PBHs in stochastic inflation, we
make use of the $\delta N$ formalism [279, 283, 284, 281, 285], in which the
primordial curvature perturbation is identified with the perturbation in the
number of e-folds $N=\int H\,\mathrm{d}t$ with respect to the average,
$\zeta={\cal N}-\langle{\cal N}\rangle$ (where $\mathcal{N}$ denotes the
duration of inflation starting from an initial field configuration, measured
in number of $e$-folds $N$). One thus has to extract the distribution
$P_{\phi}(\mathcal{N})$ of the number of e-folds of expansion $\mathcal{N}$
from a given inflaton field value $\phi$ to the end of inflation. This is the
program of the stochastic $\delta N$ formalism, in which it can be shown [278,
275, 143, 286] that the Fokker-Planck equation driving the stochastic
evolution of the coarse-grained field with respect to the integrated expansion
during inflation, $\partial P_{N}(\phi)/\partial
N=\mathcal{L}_{\mathrm{FP}}(\phi)\cdot P_{N}(\phi)$, can be recast in terms of
the adjoint Fokker-Planck operator acting on the duration distribution
$P_{\phi}(\mathcal{N})$, such that $\partial
P_{\phi}(\mathcal{N})/\partial\mathcal{N}=\mathcal{L}_{\mathrm{FP}}^{\dagger}(\phi)\cdot
P_{\phi}(\mathcal{N})$.

By solving this equation in different setups, one can show that quantum
diffusion leads unavoidably to non-Gaussian tails [275, 143]. In particular,
if the classical evolution gives a finite number of e-folds $N_{\rm c}$ from a
given field value to the end of inflation, then the large-${\cal N}$ limit
(where ${\cal N}\gg N_{\rm c}$) is inevitably dominated by quantum diffusion
since it requires many quantum kicks to evade the classical drift towards the
end of inflation. Quite generally, the probability distribution approaches an
exponential tail, $P_{\phi}({\cal
N})=\sum_{n}a_{n}e^{-\Lambda_{n}\mathcal{N}}\propto e^{-\Lambda_{0}{\cal N}}$
as ${\cal N}\to\infty$ [275, 143, 287], which cannot be described by a quasi-
Gaussian distribution (see Fig. 8). Thus predictions for the abundance of
PBHs, which will be described by the integrated probability distribution
function above a given threshold, may differ by several orders of magnitude
from results obtained assuming a Gaussian distribution [275, 143, 270, 288,
289, 290].

Let us also stress that these exponential tails cannot be properly described
by usual, perturbative parametrisations of non-Gaussian statistics (such as
those based on computing the few first moments of the distribution and the
non-linearity parameters $f_{\mathrm{NL}}$, $g_{\mathrm{NL}}$, etc.), which
can only account for polynomial modulations of Gaussian tails. A non-
perturbative approach, such as the one presented here, is therefore necessary.

It is also worth mentioning that the details of the tail (such as the value of
$a_{0}$ and $\Lambda_{0}$) are determined by the specifics of the inflationary
model under consideration. For instance, in the simplest setup where inflation
is realised with a single scalar field, if the potential energy $V$ is
constant across a field region of size $\Delta\phi$, then [275]
$\Lambda_{0}=V/(96M_{\rm{Pl}}^{2}\Delta\phi^{2})$ and
$\Lambda_{n}=(2n+1)^{2}\Lambda_{0}$. Wider “quantum wells” (i.e. larger values
of $\Delta\phi$) thus lead to smaller $\Lambda_{0}$ and to wider tails. This
result is valid if the initial velocity of the inflaton into the quantum well
is vanishing. Otherwise, a transient period of “ultra slow roll” takes place
until the inflaton relaxes to the slow-roll attractor, and on top of the
eigenvalues $\Lambda_{n}$ given above, a second set of decay rates
$\Lambda_{n}^{\prime}=\Lambda_{n}+3$ also arises [273]. When the initial
velocity is large compared to the velocity that would be required to cross the
well classically, the amplitudes $a_{n}^{\prime}$ associated with those
eigenvalues are larger than those for the rates $\Lambda_{n}$, hence the tails
are more suppressed. Another way to suppress the tails is to add a tilt in the
inflationary potential,
$V=24\pi^{2}M_{\rm{Pl}}^{4}v_{0}(1+\alpha\phi/M_{\rm{Pl}})$, for which the
decay rates are all shifted by a constant value [143], i.e.
$\Lambda_{n}^{\prime}=\Lambda_{n}+\alpha^{2}/(4v_{0})$. Unless $\alpha^{2}$ is
fine tuned to values less than the potential energy in Planckian units, the
exponential tails are heavily suppressed. As a last example, let us mention
the case where the potential features an inflection point, close to which one
can approximate
$V=24\pi^{2}M_{\rm{Pl}}^{4}v_{0}[1+\beta(\phi/{M_{\rm{Pl}}})^{3}]$. In such
models, if slow roll is not violated when approaching the inflection point
(which depends on the details of the inflationary potential away from that
point), one finds [143]
$\Lambda_{n}=(3/2)^{2/3}\pi^{2}(v_{0}\alpha)^{1/3}(n+1/2)^{2}$, such that
primordial black holes are always overproduced.

These examples show that the abundance of primordial black holes is ultimately
related to the details of the inflationary potential where quantum diffusion
plays an important role. PBHs would thus provide access to regions of the
inflationary potential that can otherwise not be probed with other
cosmological measurements such as the cosmic microwave background.

####  2.3.3 Stochastic spectator field

Another recent scenario [144] takes the advantage of the dynamics of a
stochastic spectator field to produce PBHs, while the primordial power
spectrum remains at the CMB level almost everywhere in our Universe, except in
subdominant regions where PBHs can form. If PBHs are linked to the generation
of the baryon asymmetry of the Universe, thereby connecting 1) the PBHs and
baryon abundance, 2) the baryon-to-photon ratio to the PBH density at the time
of their formation, the stochastic and quantum nature of the spectator field
fluctuations allows one to invoke an anthropic selection mechanism to resolve
the long-standing parameter fine-tuning issues related to PBH formation [148,
144, 149].

In this scenario, the stochastic spectator field experiences quantum
fluctuations during inflation that lead it to explore a wide range of the
potential, including its slow-roll part, in different Hubble patches. But the
field does not have any impact on the inflationary dynamics, thus being a
spectator field. Because the number of patches is huge, there necessarily
exist some in which the spectator field, at the time of Hubble exit, has the
required mean value for subsequent quantum fluctuations (but sub-Hubble at
that time) to later induce large curvature fluctuations. Our observable
Universe would correspond to one of these patches.

After inflation, the spectator field and its fluctuations within these patches
remain frozen during the radiation era until its potential energy starts to
dominate the density of the Universe, well after inflation but also possibly
well before the time of PBH formation. At this point, in the regions where the
field is in the slow-roll part of the potential, a very short second
inflationary phase (for at most a few e-folds) takes place, whereas in the
rest of the patch the field quickly rolls down its potential without
inflating. Such an extra expansion corresponds to the production of local non-
linear curvature fluctuations that are still super-horizon at that time. Only
later they re-enter the horizon and collapse into PBHs. But in the rest of the
Universe the curvature fluctuations are statistically Gaussian and behave as
expected in standard slow-roll inflation, unaffected by the spectator field.

The probability to form a PBH in a given region (that can be connected to the
function $\beta$ introduced in the next Section) has been calculated as [144]

$P_{\rm PBH}=\sqrt{\frac{2}{\pi}}\frac{\Delta\psi^{\rm
sr}}{\sqrt{H_{N}^{2}/4\pi^{2}+\langle\delta\psi^{2}\rangle_{N-1}}}\times\exp\left[-\frac{\langle\psi\rangle^{2}}{2(H_{N}^{2}/4\pi^{2}+\langle\delta\psi^{2}\rangle_{N-1})}\right],$
(2.18)

where $\Delta\psi^{\rm sr}$ is the size of the slow-roll region where the
field can induce $\mathcal{O}(1)$ curvature fluctuations. The subscript $N$ is
used to indicate that a quantity is evaluated at the e-fold time at which the
scale exited the horizon during inflation. It can be noticed that the
mathematical form of $P_{\rm PBH}$ is quite similar to the one of $\beta$ (see
next Section) in more standard models of Gaussian inflationary perturbations.
As a consequence, the shape of the PBH mass function is also very similar and
driven by the evolution of $H$ during inflation. There is also no need to
precisely fix the potential parameters in this model, the only condition being
that the potential is of plateau or small-field type and that it can dominate
the density of the Universe prior the QCD epoch (if one is interested to get
stellar-mass or heavier PBHs).

###  2.4 Curvaton

A different possibility to generate PBHs, while keeping the generation
mechanisms for curvature perturbations at different scales decoupled, is
represented by a simple modification of the original curvaton scenario (in the
original curvaton model [291, 292], one field, the inflaton, is responsible
only for sustaining the inflationary background evolution while the other
field, the curvaton, produces the primordial curvature perturbations). In this
setup with two fields, the primordial curvature perturbations on CMB scales
are produced by the inflaton, which acts very similarly to the standard
single-field scenario, while the curvaton field becomes responsible for
perturbations on smaller scales and thus formation of PBHs.

The power spectrum of curvature perturbations will be given as the sum of the
contributions of the inflaton and of the curvaton

$\mathcal{P}(k)=\mathcal{P}_{\text{infl}}(k)+\mathcal{P}_{\text{curv}}(k)\,.$
(2.19)

where $\mathcal{P}_{\text{curv}}(k)$ can be described by a lognormal function
such as in Eq. 2.17.

At large scales (reasonably $k\lesssim 1$ Mpc-1), the first term will be
dominating, with its magnitude being determined by CMB normalization, that is
$\mathcal{P}_{\rm infl}(k)\simeq 2\times 10^{-9}$. Moving to smaller scales
(larger $k$), the second contribution, given by the curvaton field, increases
and becomes dominant.

So, as long as the curvature perturbation on large cosmological scales is
dominated by inflaton fluctuations with an almost scale-invariant spectrum,
the curvaton field is not required to be light during inflation and its
fluctuations may have a steep blue spectrum. Curvaton fluctuations could then
give the dominant contribution to the primordial density perturbations on
small scales after the curvaton decays. The curvaton may have a steep blue
spectrum either due to interactions with the inflaton or other fields evolving
during inflation [293], a non-trivial kinetic term [294] or in an axion-like
model, where the curvaton is identified with the phase of a complex field
whose modulus decreases rapidly during inflation [295].

An important characteristic of any curvaton model is that the fluctuations
originate from non-adiabatic field fluctuations during inflation, a fact which
gives rise to a non-Gaussian distribution in the primordial curvature
perturbations, including a finite non-Gaussianity in the squeezed limit.
Although this local-type non-Gaussianity is typically small ($f_{\rm NL}\sim
1$), it can nonetheless have a significant effect on the abundance of PBHs, as
discussed in Subsection 2.3.

###  2.5 Preheating

Another way through which one can produce PBHs is the preheating instability.
During the preheating period after inflation, when the inflaton oscillates
coherently at its ground state and decays to other degrees of freedom, it is
widely argued that resonant amplification of quantum field fluctuations,
responsible for particle production, take place [296, 297]. These amplified
quantum fluctuations should be accompanied by a resonant amplification of the
scalar metric fluctuations (usually quoted as metric preheating [298, 299,
300, 301]), responsible for gravitational fluctuations in the curvature, since
the two are coupled through Einstein’s equations. One then expects that the
amplified metric perturbations give rise to large curvature/matter
perturbations, which in their turn collapse and form PBHs.

Historically, PBHs emanating from the preheating instability were speculated
in the context of multi-field inflation and in particular in the context of
two-field chaotic inflation [302, 303, 304, 305] and more recently in [306,
307], since in this case the parametric amplification of entropy/isocurvature
fluctuations can source the parametric amplification of the
adiabatic/curvature fluctuations in the regime of broad resonance. Thus, these
amplified curvature fluctuations could break the scale-invariance of the
primordial power spectrum and give birth to PBHs through non-linear
gravitational fluctuations inducing large density contrasts [308, 309].

It is predicted as well [310] (see also [311]) that in the context of single-
field inflation, there is a pronounced resonant instability structure in the
narrow regime where amplified metric perturbations can induce the production
of PBHs as recently studied in [312, 313]. At this point, one should point out
that the narrow resonant structure of metric preheating was shown in [314] to
be immune to the decay of the inflaton into a radiation fluid, ensuring in
this way the transition to the Hot Big Bang phase of the Universe.

##### The case of metric preheating

Most of the studies in which PBHs are produced during preheating treat them in
the context of multi-field and especially in the context of two-field
inflation. This makes it difficult to predict in an analytic way the matter
power spectrum responsible for the PBH formation, which is constructed
consequently numerically. On the contrary, in the context of single-field
inflation, one can in principle extract in a direct way the matter power
spectrum responsible for the PBH production. In what follows we emphasize on a
specific scenario studied in [312, 313], in which PBHs are produced during the
preheating phase due to the emergence of a resonant instability structure
concerning the equation of motion of the scalar metric perturbations that are
enhanced on small scales (metric preheating), which exit the Hubble radius
close to the end of inflation. Specifically, the modes concerned satisfy the
relation

$aH<k<a\sqrt{3Hm},$ (2.20)

where $k$ is the comoving scale concerned, $a$ is the scale factor, $H$ is the
Hubble parameter and $m$ is the mass parameter of the inflationary potential.
For these modes, as shown in [312], the density contrast scales linearly with
the scale factor, i.e. $\delta_{k}\sim a$, manifesting in this way an
effective “dust” behavior during the instability phase in which PBHs can be
easily produced.

Considering only the modes entering the instability band from above, namely
the left hand side of the inequality (2.20), the matter power spectrum
responsible for the PBH formation at the end of the preheating instability
phase can be written analytically as [312]

$\mathcal{P}_{\delta}\left(t_{\mathrm{\Gamma}},\frac{k}{k_{\mathrm{end}}}\right)=\left(\frac{k}{k_{\mathrm{end}}}\right)^{4}\left(\frac{\rho_{\mathrm{inf}}}{\rho_{\Gamma}}\right)^{2/3}\left(\frac{6}{5}\right)^{2}\mathcal{P}_{\zeta,\mathrm{end}}\left(\frac{k}{k_{\mathrm{end}}}\right),$
(2.21)

where $k_{\mathrm{end}}$ is the comoving scale exiting the Hubble radius at
the end of inflation, $\rho_{\mathrm{inf}}$ is the energy scale at the end of
inflation, $\rho_{\Gamma}$ is the energy scale at the end of the preheating
instability phase and
$\mathcal{P}_{\zeta,\mathrm{end}}\left(\frac{k}{k_{\mathrm{end}}}\right)$ is
the curvature power spectrum at the end of inflation which in the slow-roll
approximation reads as [315, 316]

$\mathcal{P}_{\zeta,\mathrm{end}}\left(\frac{k}{k_{\mathrm{end}}}\right)=\frac{H_{*}^{2}\left(k\right)}{8\pi^{2}M_{\rm{Pl}}^{2}\epsilon_{1*}\left(k\right)}\left[1+\left(\frac{k}{k_{\mathrm{end}}}\right)^{2}\right]\left[1-2\left(C+1\right)\epsilon_{1*}\left(k\right)-C\epsilon_{2*}\left(k\right)\right]\
\text{for}\ k<k_{\mathrm{end}}.$ (2.22)

The functions $H_{*}(k)$, $\epsilon_{1*}(k)$ and $\epsilon_{2*}(k)$ denote
respectively the values of $H$, $\epsilon_{1}$ and $\epsilon_{2}$ at the time
when the mode $k$ exits the Hubble radius during inflation. The parameter
$C\simeq-0.7296$ is a numerical constant.

The relevant range of the density contrast at the end of the preheating phase
is

$\displaystyle\left(\frac{\rho_{\mathrm{inf}}}{\rho_{\mathrm{\Gamma}}}\right)^{1/3}\left(\frac{k}{k_{\mathrm{end}}}\right)^{2}\left(\frac{3\pi}{2}\right)^{2/3}\left[\left(\frac{k}{k_{\mathrm{end}}}\right)^{3}\sqrt{\frac{\rho_{\mathrm{inf}}}{\rho_{\Gamma}}}-1\right]^{-2/3}<\delta_{k}(t_{\mathrm{\Gamma}})<1,$
(2.23)

while the relevant parameter space
$(\rho_{\mathrm{\Gamma}},\rho_{\mathrm{inf}})$ so that we do not have PBH
overproduction during preheating reads as

$\displaystyle\frac{\rho_{\Gamma}}{M_{\rm{Pl}}^{4}}<\frac{4}{125\sqrt{3}\pi^{5}}\left(\frac{\rho_{\mathrm{inf}}}{M_{\rm{Pl}}^{4}}\right)^{5/2}\
.$ (2.24)

Below, in Fig. 9 we see the behavior of
$\mathcal{P}_{\delta}\left(t_{\mathrm{\Gamma}},\frac{k}{k_{\mathrm{end}}}\right)$
[see Eq. 2.21] by fixing $\rho^{1/4}_{\mathrm{inf}}=10^{14}\mathrm{GeV}$ and
$\rho^{1/4}_{\mathrm{\Gamma}}=10^{7}\mathrm{GeV}$. At this point, it is
important to stress out that exactly at $k=k_{\mathrm{end}}$, Eq. (2.22) is
not trustful since there the slow-roll approximation breaks down.

Figure 9: The matter power spectrum at the end of the preheating instability
phase. The energy scales at the end of inflation and at the end of the
preheating instability phase are fixed as
$\rho^{1/4}_{\mathrm{inf}}=10^{14}\mathrm{GeV}$ and
$\rho^{1/4}_{\mathrm{\Gamma}}=10^{7}\mathrm{GeV}$.

The corresponding mass of the PBHs produced is related to the mass inside the
Hubble radius at the time when a scale $k$ re-enters the horizon, namely

$m_{\mathrm{PBH}}=\gamma
M_{\mathrm{H}}=\gamma\frac{\left(3M^{2}_{\mathrm{pl}}\right)^{3/2}}{\sqrt{\rho_{\mathrm{inf}}}}\left(\frac{k}{k_{\mathrm{end}}}\right)^{-3},$
(2.25)

where $\gamma$ is an efficiency parameter of order one depending on the
details of the gravitational collapse. Here, it is important to notice that in
(2.25), one clearly sees the dependence of the PBH mass on the energy scale at
the end of inflation and as a consequence on the underlying inflationary
model.

Regarding the characteristic PBH mass produced during preheating one can
choose $k/k_{\mathrm{end}}=0.1$, and using a relative high energy scale at the
end of inflation $\rho^{1/4}_{\mathrm{inf}}=10^{14}\mathrm{GeV}$, from Eq.
(2.25) one obtains $m_{\mathrm{PBH}}\simeq 10^{-27}M_{\odot}$. This mass,
around which the matter power spectrum peaks, is very small compared to the
mass range LISA can detect through the GW portal of PBH mergers. However, as
pointed out recently in [313], in which the PBH mass function is computed more
accurately with the use of the excursion-set formalism, the mass of PBHs
produced from the metric preheating instability, at which the PBH mass
function peaks, is higher than the PBH mass quoted above and can be as high as
the solar mass, lying in this way within the mass range probed by LISA.

At this point, one should also stress that ultralight PBHs, such as the ones
produced from the preheating instability, can be potentially detected by LISA
through the stochastic background of GWs induced at second order from PBH
isocurvature perturbations, as argued in [317]. For these GWs, the peak
frequency could lie within the LISA band despite the smallness of the PBH
mass.

###  2.6 Early matter era

An early prolonged period of matter domination – after the end of inflation –
could arise due to the domination of moduli fields which are a natural
prediction of any supergravity or superstring-inspired theory of the early
Universe. Since the curvature threshold for the PBH collapse depends on the
equation of state of the Universe, it might be easier to obtain the same
abundance of PBH during the moduli dominated epoch from a peak in the
inflationary power spectrum less strong than the typical $O(10^{-2})$ required
in the standard scenario of radiation domination. Several mechanisms have been
proposed to lead to such a production during an early matter dominated era,
see e.g. [22, 25, 24, 318, 319, 320, 321, 322, 323, 324, 325]. In such
scenarios, the mass and fraction of PBHs produced from a given inflationary
potential depend on whether the moduli dominated epoch is immediate or delayed
after the end of inflation. Recently, it has also been pointed out that even
if PBHs are formed during radiation domination, the presence of a reheating
epoch can still crucially affect the PBH abundance [326, 327].

In models of inflation based on moduli stabilisation in string theory
compactifications, for example, the inflaton may be identified with a modulus
field which is only gravitationally coupled to the visible sector, leading to
a slow-reheating phase, where the inflaton decay is via Planck-suppressed
operators, and hence a prolonged duration of matter-domination after inflation
[24, 323]. Ref. [323] finds that PBH formation is most effective for reheating
temperatures $\lesssim 10^{6}$ GeV, in which case PBHs can form a significant
fraction of the dark matter density today when primordial curvature
fluctuations from inflation are of order $10^{-4}$ in a matter-dominated era,
rather than $\sim 10^{-2}$ usually required in a radiation-dominated era.

###  2.7 Phase transitions

The formation and abundance of PBHs can also be affected by phase transitions
in the early Universe. During the radiation era, the formation of PBHs may be
facilitated if pre-existing super-horizon density fluctuations enter the
cosmological horizon during a first order phase transition which proceeds in
approximate equilibrium [25, 54]. During this phase transition the speed of
sound tends to zero and the pressure response of the fluid vanishes and does
not counter-balance the gravitational collapse of horizon-sized primordial
over-densities. Detailed numerical simulations using general-relativistic
hydrodynamics show that the overdensity threshold required for PBH formation
by fluctuations that enter the cosmological horizon around the phase
transition is smaller than the ordinary value for PBH formation during the
radiation epoch. This would produce an increase in the PBH mass spectrum
around the corresponding horizon mass in which the phase transition occurs.

It has been argued that the QCD phase transition could produce an enhancement
of the population of PBH’s around $1M_{\odot}$ [20], but lattice QCD
calculations have shown that the QCD phase transition in the Standard Model is
not first order but a continuous cross-over [328]. Nevertheless, it is
possible to collapse a significant amount of PBHs at the QCD transition if
there are large fluctuations and radiation pressure drops as the relativistic
quark and gluon degrees of freedom disappear from the plasma [145, 70, 73]

Alternative scenarios have been proposed, where the dynamics of the QCD phase
transition incorporates extra degrees of freedom, like lepton number violation
[329] and solitosynthesis [262]. In other models, the QCD transition could be
first order if the number of light quarks were larger than three at the
confinement. Models in which the dynamics of a light scalar field have been
tailored to suppress quark masses before the QCD transition, while recovering
the measured values after the transition, have been developed in [330].

Another possibility is that the collision of bubbles formed during a first
order phase transition in the early Universe could create a sufficiently large
local concentration of the energy stored in the walls to cause gravitational
collapse and the formation of a primordial black hole [331, 332, 333]. This
process could be affected by spatial clustering properties of bubble generated
in high-energy phase transitions [334, 335], Recent models [336] consider the
possibility that there may be particles whose mass increases significantly
during the phase transition, suppressing the transmission of the corresponding
particles through the advancing bubble walls. This effect can build a
sufficiently large overdensity in front of the walls that collapse into a
primordial black hole. The primordial black hole density and mass distribution
depend on the model parameters. Other models [337] consider fermion species
that gain big masses in the true vacuum, so that the corresponding particles
get trapped in the false vacuum as they do not have sufficient energy to
penetrate the bubble wall. After the first order phase transition, the
fermions are compressed into the false vacuum remnants to form Fermi-balls
that finally collapse to PBH due to the Yukawa attractive force.

Higher-order phase transitions periods in the early Universe may also lead to
the formation of PBHs. For example, domain walls can be generated during a
non-equilibrium second order phase transition in the vacuum state of a scalar
field with a flat direction during inflation [338]. The background de Sitter
fluctuations of such an effectively massless scalar field could provide
non–equilibrium redefinition of correlation length and give rise to the
islands of one vacuum in the sea of another one. After the phase transition
takes place in the Friedman– Robertson-Walker (FRW) epoch, the two vacua are
separated by a wall. Some of these closed and sufficiently large walls could
accumulate enough energy to collapse and form a massive black hole. The mass
spectrum of the PBHs which can be created in such a way depends on the details
of the potential of the scalar field that has the flat direction during
inflation and triggers the phase transition.

###  2.8 Domain walls and cosmic strings

Topological defects in the early Universe may also lead to the production of
PBHs. S. W. Hawking [339] was the first to show that PBHs could form from the
collapse of cosmic string loops. He argued that loops that collapse by a
factor $(G\mu)^{-1}$, where $\mu$ is the mass per unit length of the strings,
will inevitably collapse to form black holes. PBH production can be used to
place bounds on the cosmic string network [340, 26, 341, 342, 343, 344, 345,
346] or cosmic strings nucleated during inflation [347]. More recently it has
been argued that cosmic string cusps can collapse gravitationally into PBHs
with a mass function that could extend up to stellar-masses [348], although
this claim has been disputed in Ref. [349]

Moreover, domain walls, which arise in theories with a broken discrete
symmetry, can also form PBHs when such domain walls are smaller than a
critical radius fall within the cosmological horizon and collapse due to their
own tension [350] or are produced by tunnelling during inflation [351]. In
such scenarios, the mass function of PBHs in general has a spikelike
structure. Moreover, in QCD axion models, PBHs can form either by the collapse
of long-lived string-domain wall networks or from the collapse of closed
domain walls [350, 352]. The mass and abundance of PBHs formed in these
scenarios crucially depends on the QCD axion mass. However, interestingly,
various observational constraints on the PBH abundance in turn limit the QCD
axion parameter space.

###  2.9 Primordial magnetic fields

Primordial magnetic fields generated in the early Universe are considered one
of the feasible candidates to provide the required seeds for the observed
large-scale intergalactic magnetic fields. In addition, such primordial fields
also generally induce an anisotropic stress that can act as a source for the
evolution of the super-horizon curvature perturbations, leading to the
formation of PBHs with a broad possible range of masses. Since a large
abundance of PBHs requires a large amplitude of density perturbations, the
strength of these primordial magnetic fields on small scales should also be
large enough for this mechanism to work [353]. Based on this idea, strong
constraints on the amplitude of primordial magnetic fields have been
established, a few hundreds nano-Gauss on scales in the wave-number range
$10^{2}\,{\rm Mpc}^{-1}\leq k\leq 10^{18}\,{\rm Mpc}^{-1}$.

Alternatively, it has been proposed that PBHs can act as seeds of cosmic
magnetic fields due to their accretion disks, either due to a Biermann battery
[354, 355] mechanism or more exotic mechanisms such as monopoles [356, 357] or
Kerr-Newman PBHs [358]. In particular, in [355] it was shown that locally
isothermal disks around supermassive PBHs with masses $M>10^{10}M_{\odot}$,
can generate a seed primordial magnetic field of the order of
$10^{-30}\mathrm{G}$, which can be later amplified by various
dynamo/instability processes and provide the seed for the present day magnetic
field in intergalactic scales. Interestingly enough, this population of
magnetised supermassive PBHs can induce a stochastic GW background which can
be used to set tight constraints on the supermassive PBH abundances [See [359]
for more details].

###  2.10 Summary

In this Section, we have explored multiple mechanisms described in the
literature for generating primordial black holes, from gravitational collapse
of primordial fluctuations and inflation to preheating, phase transitions,
collisions of bubble walls and cosmic strings. Each of these mechanisms has
its own mass spectrum and individual characteristics, which makes it possible
to differentiate between them. The most predictive one is that of primordial
fluctuations, where specific models exist that can be put to test by present
and future observations. But we also emphasized the possible role of non-
Gaussian tails in the distribution of curvature fluctuations, which arise in
models with an ultra-slow-roll phase of inflation or in models with a
stochastic spectator field. This is an interesting, recent and active research
area since this could soften or even resolve the fine-tuning issues associated
to the formation of PBHs with significant abundances. In the next Sections, we
will explore the dynamics of gravitational collapse and how it depends on the
thermal history of the Universe, with its impact on the PBH mass function,
their spin distribution and merging rates.

##  3 PBH formation, mass function and clustering

###  3.1 From curvature perturbations to density contrast

PBHs form after cosmological perturbations re-enter the cosmological horizon.
Assuming spherical symmetry on superhorizon scales, the local region of the
Universe describing such perturbations has the following asymptotic form of
the metric

$\displaystyle{\rm d}s^{2}$ $\displaystyle=$ $\displaystyle-{\rm
d}t^{2}+a^{2}(t)\left[\frac{{\rm d}r^{2}}{1-K(r)r^{2}}+r^{2}{\rm
d}\Omega^{2}\right]=-{\rm d}t^{2}+a^{2}(t)e^{2\zeta(\hat{r})}\left[{\rm
d}\hat{r}^{2}+\hat{r}^{2}{\rm d}\Omega^{2}\right],$ (3.1)

where $a(t)$ is the scale factor, while $K(r)$ and $\zeta(\hat{r})$ are the
conserved comoving curvature perturbations, defined on a super-Hubble scale,
converging to zero at infinity where the universe is taken to be unperturbed
and spatially flat. The equivalence between the radial and the angular parts
of these two forms of the metric gives

$\left\\{\begin{aligned} &r=\hat{r}e^{\zeta(\hat{r})}\,,\\\
&\displaystyle{\frac{{\rm d}r}{\sqrt{1-K(r)r^{2}}}}=e^{\zeta(\hat{r})}{\rm
d}\hat{r}\,\end{aligned}\right.$ (3.2)

and the difference between the two Lagrangian coordinates $r$ and $\hat{r}$ is
related to the parametrisation of the comoving coordinate, fixed by the
curvature perturbation chosen in the metric (i.e. $K(r)$ or $\zeta(\hat{r})$).
From a geometrical point of view, the coordinate $\hat{r}$ considers the
perturbed region as a local FLRW separated universe with the curvature
perturbation $\zeta(\hat{r})$ modifying the local expansion, while the
curvature profile $K(r)$ is defined with respect to the background FLRW
solution ($K=0$) and measures more directly the spatial geometry of space-
time.

On superhorizon scales, where the curvature profile is time independent, these
two are related as

$K(r)r^{2}=-\hat{r}\zeta^{\prime}(\hat{r})\left[2+\hat{r}\zeta^{\prime}(\hat{r})\right]\,.$
(3.3)

In this regime, we can use the gradient expansion approach [48, 360, 280, 50],
based on expanding the time dependent variables, like energy density and
velocity profiles, as power series of a small parameter $\epsilon\ll 1$ up to
the first non zero order, where $\epsilon$ is conveniently identified with the
ratio between the Hubble radius and the length scale of the perturbation.
Although this approach reproduces the time evolution of the linear
perturbation theory, it also allows having non linear curvature perturbations
if the spacetime is sufficiently smooth on the scale of the perturbation (see
[285]). This is equivalent to saying that in this regime pressure gradients
are negligible, not playing an important role in the perturbation evolution,
which grows with the Universe expansion in a self similar way.

In this approximation, the energy density profile can be written as
555$K^{\prime}(r)$ denotes differentiation with respect to $r$ while
$\zeta^{\prime}(\hat{r})$ and $\nabla^{2}\zeta(\hat{r})$ denote
differentiation with respect to $\hat{r}$. [361, 57]

$\displaystyle\frac{\delta\rho}{\rho_{\rm b}}$ $\displaystyle\equiv$
$\displaystyle\frac{\rho(r,t)-\rho_{\rm b}(t)}{\rho_{\rm
b}(t)}=\frac{1}{a^{2}H^{2}}\frac{3(1+w)}{5+3w}\frac{\left[K(r)\,r^{3}\right]^{\prime}}{3r^{2}}$
(3.4) $\displaystyle=$
$\displaystyle-\frac{1}{a^{2}H^{2}}\frac{4(1+w)}{5+3w}e^{-5\zeta(\hat{r})/2}\nabla^{2}e^{\zeta(\hat{r})/2}\,,$

where $H(t)=\dot{a}(t)/a(t)$ is the Hubble parameter. The parameter $w$ is the
coefficient of the equation of state $p=w\rho_{\rm b}$, relating the total
(isotropic) pressure $p$ to the total background energy density $\rho_{\rm
b}$.

In the linear regime of curvature perturbations, when $\zeta\ll 1$, this
expression can be written as

$\frac{\delta\rho}{\rho_{\rm
b}}\simeq-\frac{1}{a^{2}H^{2}}\frac{2(1+w)}{5+3w}\nabla^{2}\zeta(\hat{r})=-\frac{k^{2}}{a^{2}H^{2}}\frac{2(1+w)}{5+3w}\zeta_{k}\,,$
(3.5)

where the second equality is obtained from the Fourier transformation
$\nabla^{2}\zeta=k^{2}\zeta_{k}$, showing that in the linear regime there is a
simple one to one mapping between the real space, where the perturbations
collapses, and the Fourier space where the power spectrum of cosmological
perturbations is defined. If the curvature perturbation $\zeta$ is a Gaussian
variable, the density contrast $\delta\rho/\rho_{\rm b}$ has also a Gaussian
distribution within the linear regime described by (3.5). However, because the
amplitude of the threshold $\delta_{\rm c}$ for PBHs is non linear, the linear
regime does not give an accurate description of the statistics of the density
contrast, and the mapping between the Fourier space and the real space
requires a more elaborate approach. As we will see later, this allows to
compute the threshold $\delta_{\rm c}$ from the shape of the power spectrum
$\mathcal{P}_{\zeta}$.

For this reason, the definition of the threshold in terms of $\zeta_{\rm c}$
computed from the linear approximation given by (3.5), does not give an
accurate description and should be dropped in favour of $\delta_{\rm c}$ which
is clearly defined in the following Section. We anticipate that the nonlinear
relation between the density contrast and the comoving curvature perturbations
causes the former to obey non-Gaussian statistics, even in absence of
primordial non-Gaussianity. Computations based on peak theory and on threshold
statistics show that this unavoidable non-Gaussianity makes the production of
PBHs more difficult. In case of a peaked power spectrum, this effect is
compensated by an increase of the amplitude by an order 2-3 factor [362, 363].

###  3.2 Perturbation amplitude

To define the amplitude of a cosmological perturbation, it is useful to
introduce the compaction function [48, 57], defined as

$\mathcal{C}\equiv 2\frac{\delta M(r,t)}{R(r,t)}\,,$ (3.6)

where $R(r,t)$ is the areal radius and $\delta M(r,t)$ is the difference
between the Misner-Sharp mass within a sphere of radius $R(r,t)$, and the
background mass $M_{\rm b}(r,t)=4\pi\rho_{\rm b}(r,t)R^{3}(r,t)/3$ within the
same areal radius but calculated with respect to a spatially flat FLRW metric.
In the superhorizon regime (i.e. $\epsilon\ll 1$) the compaction function is
time independent, and directly related to the curvature profile as

$\mathcal{C}(r)=\frac{3(1+w)}{5+3w}K(r)r^{2},$ (3.7)

where this expression can be written in terms of $\zeta(\hat{r})$ using (3.3).

The comoving length scale of the perturbation can be conveniently identified
with the distance where the compaction function reaches its peak (which is a
maximum for a positive perturbation and a minimum for a negative one) when
$r=r_{\rm m}$ [57](i.e. $\mathcal{C}^{\prime}(r_{\rm m})=0$), characterized by
a constraint relation in terms of the curvature profile

$\displaystyle K(r_{\rm m})+\frac{r_{\rm m}}{2}K^{\prime}(r_{\rm
m})=0,\quad\textrm{or}\quad\zeta^{\prime}(\hat{r}_{\rm m})+\hat{r}_{\rm
m}\zeta^{\prime\prime}(\hat{r}_{\rm m})=0\,.$ (3.8)

Given the curvature profile, the parameter $\epsilon$ of the gradient
expansion is therefore defined as

$\epsilon\equiv\frac{R_{\rm H}(t)}{R_{\rm b}(r_{\rm m},t)}=\frac{1}{aHr_{\rm
m}}=\frac{1}{aH\hat{r}_{\rm m}e^{\zeta(\hat{r}_{\rm m})}}\,,$ (3.9)

where $R_{\rm H}=1/H$ is the cosmological horizon and $R_{\rm b}(r,t)=a(t)r$
is the background component of the areal radius. With these definitions, the
expression written in Eq. (3.4) is valid with good accuracy for $\epsilon\ll
1$.

According to this, we can now consistently define the perturbation amplitude
as being the mass excess of the energy density within the scale $r_{\rm m}$
when $\epsilon=1$ ($aHr_{\rm m}=1$). This corresponds to defining the horizon
crossing time $t_{\rm H}$ with a linear extrapolation from the superhorizon
regime. Although this is not very accurate, it provides a well defined
criterion to measure and compare consistently the amplitude of different
perturbations, understanding how the threshold is varying with respect to the
initial curvature profiles.

The amplitude of the perturbation measured at $t_{\rm H}$, which we refer to
as $\delta_{\rm m}\equiv\delta(r_{\rm m},t_{\rm H})$, is given by the excess
of mass averaged over a spherical volume of radius $R_{\rm m}$, defined as

$\delta_{\rm m}=\frac{4\pi}{V_{R_{\rm m}}}\int_{0}^{R_{\rm
m}}\frac{\delta\rho}{\rho_{\rm b}}\,R^{2}{\rm d}R\,=\frac{3}{r_{\rm
m}^{3}}\int_{0}^{r_{\rm m}}\frac{\delta\rho}{\rho_{\rm b}}\,r^{2}{\rm d}r\,,$
(3.10)

where $V_{R_{\rm m}}={4\pi}R_{\rm m}^{3}/3$. The second equality is obtained
by neglecting the higher order terms in $\epsilon$, approximating $R_{\rm
m}\simeq a(t)r_{\rm m}$, which allows to simply integrate over the comoving
volume of radius $r_{\rm m}$. Inserting the expression for
$\delta\rho/\rho_{\rm b}$ given by (3.4) into (3.10), one obtains $\delta_{\rm
m}=\mathcal{C}(r_{\rm m})$, and with simple calculation seen in [57], one
obtains the fundamental relation

$\delta_{\rm m}=3\frac{\delta\rho}{\rho_{\rm b}}(r_{\rm m},t_{\rm H}),$ (3.11)

showing that at $r_{\rm m}$ the perturbation amplitude is not affected by the
location where this is measured.

###  3.3 Threshold for PBH formation

PBHs form when the perturbation amplitude $\delta_{\rm m}>\delta_{\rm c}$,
where the value of the threshold $\delta_{\rm c}$ depends on the shape of the
energy density profile and the equation of state. For a radiation-dominated
Universe ($p=\frac{1}{3}\rho$), it has been found that $2/5\leq\delta_{\rm
c}\leq 2/3$, with the shape of the cosmological perturbation simply
characterised by one dimensionless parameter corresponding to the width of the
peak of the compaction function [57, 58]

$\alpha=-\frac{\mathcal{C}^{\prime\prime}(r_{\rm m})r_{\rm
m}^{2}}{4\mathcal{C}(r_{\rm m})}\,.$ (3.12)

This allows to compute the threshold $\delta_{\rm c}$ as a function of the
shape parameter $\alpha$, with an analytic expression up to a few percent
precision [58]

$\delta_{\rm
c}\simeq\frac{4}{15}e^{-\frac{1}{\alpha}}\frac{\alpha^{1-\frac{5}{2\alpha}}}{\Gamma\left(\frac{5}{2\alpha}\right)-\Gamma\left(\frac{5}{2\alpha},\frac{1}{\alpha}\right)},$
(3.13)

where $\Gamma$ identifies the special Gamma-functions. This is consistent with
the numerical analysis made in [57], where it was shown that the effects on
the threshold of additional parameters are negligible, because the collapse in
spherical symmetry is basically characterized by the configuration of the
region within the forming apparent horizon ($r\leq r_{\rm m}$), whose initial
configuration is fully described by the shape parameter $\alpha$.

The mass of a PBH is determined by the amplitude $\delta$ of the perturbation
with respect to the corresponding threshold $\delta_{\rm c}$, according to the
scaling law of critical collapse

$m_{\rm PBH}=\mathcal{K}(\delta-\delta_{\rm c})^{\eta}M_{\rm H}$ (3.14)

where the critical exponent $\eta$ depends only on the equation of state
($\eta\simeq 0.36$ for a radiation-dominated Universe), while $\mathcal{K}$,
as for $\delta_{\rm c}$, depends on the initial configuration of the energy
density profile, roughly varying between $1$ and $10$. For the standard
scenario of a radiation-dominated Universe, expression (3.14) describes the
mass with good accuracy when $(\delta-\delta_{\rm c})\lesssim 10^{-2}$,
corresponding to $m_{\rm PBH}\lesssim M_{\rm H}$.

Although non linear cosmological density perturbations are described by a non
Gaussian random field, when the Universe is still radiation-dominated one can
compute the threshold $\delta_{\rm c}$ from the shape of the Gaussian
inflationary power spectrum following a simple prescription. The algorithm,
divided into a few simple steps, accounts for both the non linear effects
associated with the relation between the Gaussian curvature perturbation and
the density contrast as well as for those arising at horizon crossing. While a
more refined description of the various steps will be found in [60], we give
here a synthetic overview of the prescription one needs to follow.

  1. 1.

The power spectrum of the curvature perturbation: take the primordial power
spectrum $\mathcal{P}_{\zeta}$ of the Gaussian curvature perturbation and
compute, on superhorizon scales, its convolution with the transfer function
$T(k,\eta)$

$P_{\zeta}(k,\eta)=\frac{2\pi^{2}}{k^{3}}\mathcal{P}_{\zeta}(k)T^{2}(k,\eta).$
(3.15)

  2. 2.

The comoving length scale $\hat{r}_{\rm m}$ of the perturbation is related to
the characteristic scale $k_{*}$ of the power spectrum $P_{\zeta}$. Compute
the value of $k_{*}\hat{r}_{\rm m}$ by solving the following integral equation

$\int{\rm d}kk^{2}\\!\left[(k^{2}\hat{r}_{\rm
m}^{2}-1)\frac{\sin(k\hat{r}_{\rm m})}{k\hat{r}_{\rm m}}+\cos{(k\hat{r}_{\rm
m})}\right]\\!P_{\zeta}(k,\eta)=0\,.$ (3.16)

  3. 3.

The shape parameter: compute the corresponding shape parameter $\alpha$ of the
collapsing perturbation, including the correction from the non linear effects,
by solving the following equation

$\displaystyle F(\alpha)\left[1+F(\alpha)\right]\alpha$
$\displaystyle=2\alpha_{\text{\tiny G}}$ (3.17)

with

$\displaystyle\alpha_{\text{\tiny G}}$
$\displaystyle=-\frac{1}{4}\left[1+\hat{r}_{\rm m}\frac{\int{\rm
d}kk^{4}\cos{(k\hat{r}_{\rm m})}P_{\zeta}(k,\eta)}{\int{\rm
d}kk^{3}\sin{(k\hat{r}_{\rm m})}P_{\zeta}(k,\eta)}\right],$ $\displaystyle
F(\alpha)$
$\displaystyle=\sqrt{1-\frac{2}{5}e^{-\frac{1}{\alpha}}\frac{\alpha^{1-\frac{5}{2\alpha}}}{\Gamma\left(\frac{5}{2\alpha}\right)-\Gamma\left(\frac{5}{2\alpha},\frac{1}{\alpha}\right)}}\,.$
(3.18)

A numerical fit of the shape parameter $\alpha$ as a function of
$\alpha_{\text{\tiny G}}$, which represents the solution of Eq. (3.17) with a
percent accuracy, is given by [84]

$\alpha\simeq\left\\{\begin{aligned} &1.758\,\alpha_{\text{\tiny
G}}^{2.335}+1.912\,\alpha_{\text{\tiny
G}}\quad\quad&0.1\lesssim\,\alpha\lesssim 4.5,\\\ &4\,\alpha_{\text{\tiny
G}}^{2}+3.930\,\alpha_{\text{\tiny G}}\quad\quad\quad&\alpha\gtrsim
8.\,\,\,\,\end{aligned}\right.$ (3.19)

In the intermediate region where $4.5<\alpha<8$, no simple power law fit with
a percent accuracy was found, and one needs to solve Eq. (3.17) numerically.

  4. 4.

The threshold $\delta_{\rm c}$: compute the threshold as function of $\alpha$,
fitting the numerical simulations.

     * •

At _superhorizon scales_ making a linear extrapolation at horizon crossing
($aHr_{\rm m}=1$).

$\delta_{\rm c}\simeq\left\\{\begin{aligned}
&\alpha^{0.047}-0.50\quad&0.1\lesssim\,\alpha\lesssim\ 7\,\\\
&\alpha^{0.035}-0.475\quad\quad&7\lesssim\,\alpha\lesssim 13\\\
&\alpha^{0.026}-0.45\quad&13\lesssim\,\alpha\lesssim 30\end{aligned}\right.$
(3.20)

     * •

At _horizon crossing_ taking into account also the non linear effects.

$\delta_{\rm c}\simeq\left\\{\begin{aligned}
&\alpha^{0.125}-0.05\quad&0.1\lesssim\,\alpha\lesssim 3\\\
&\alpha^{0.06}+0.025\quad\quad&3\lesssim\,\alpha\lesssim 8\\\ &\quad\quad
1.15\quad&\alpha\gtrsim 8\end{aligned}\right.$ (3.21)

This prescription is using the fact that, if $\mathcal{P}_{\zeta}$ follows a
Gaussian distribution, $\Phi_{\rm L}\equiv-\hat{r}\zeta^{\prime}(\hat{r})$ is
also a Gaussian variable, and one can write $\delta_{\rm m}$ as

$\delta_{\rm c}=\frac{4}{3}\Phi_{\rm c,L}\left(1-\frac{1}{2}\Phi_{\rm
c,L}\right)$ (3.22)

where $\Phi_{\rm m,L}=\Phi_{\rm L}(\hat{r}_{\rm m})$, and $\Phi$ is such that
$0.37\lesssim\Phi_{\rm c,L}\leq 1$.

In Fig. 10, we show the threshold density contrast as a function of the shape
parameter (left panel). Also, in the right panel, we report the critical
collapse parameter $\mathcal{K}$ defined in Eq. (3.14). We recall that the
critical collapse exponent $\eta$ is independent from the shape of the
collapsing overdensity peak (i.e. independent from $\alpha$) and it is fixed
to be $\eta=0.358$ in a perfect radiation fluid. In the following, we will
discuss how these parameters are modified when thermal effects during the
evolution of the universe, causing a significant departure from the perfect
radiation fluid approximation.

Figure 10:  Threshold $\delta_{c}$ and critical collapse parameter ${\cal K}$
as a function of the shape parameter $\alpha$, assuming a universe dominate by
perfect radiation (i.e. $w=1/3$). Figure adapted from Ref. [60]

Before concluding this Section, we mention two relevant points. First, it has
been recently observed that corrections from the non-linear radiation transfer
function and the determination of the true physical horizon crossing tend to
decrease the PBH abundance [364]. Secondly, a recent study [365] which
investigated the formation of PBHs with an anisotropic relativistic perfect
fluid, using a covariant form of the EoS in terms of pressure and energy
density gradients. This indicates that the value of the threshold $\delta_{\rm
c}$ is increasing with the amplitude of the anisotropy, if this is not large,
while further studies are necessary for a nonlinear modification of the EoS.

###  3.4 Gravitational collapse of primordial curvature perturbations -
standard formalism

Although PBHs can be generated by many different mechanisms, we will focus
here on the collapse of cosmological perturbations. In this case, the fraction
$\beta_{\rm f}$ of PBHs formed in the early Universe is determined by the
probability that a given primordial curvature fluctuation $\zeta$,
characterized by an amplitude $\delta$ of the density contrast, is above a
certain threshold $\delta_{\rm c}$. The fraction of PBHs at formation is then
usually given by

$\beta_{\rm f}(m_{\rm PBH})=\int_{\delta_{\rm c}}^{\infty}P(m_{\rm
PBH},\delta){\rm d}\delta\,,$ (3.23)

and the abundance of PBHs is sensitive to the PDF $P(m_{\rm PBH},\delta)$ as
to the value of the threshold $\delta_{\rm c}$. In the previous section, we
have discussed in details how the amplitude $\delta$ is computed from the
curvature $\zeta$, understanding how this is related to the shape of the
cosmological power spectrum of curvature perturbations
$\mathcal{P}_{\zeta}(k)$. Later we will also discuss the effects of primordial
non-Gaussianity on the abundance.

In Eq. (3.23), we have expressed the PDF in terms of the mass $m_{\rm PBH}$ of
the PBHs formed instead of the number of $e$-folds because in first
approximation the mass of PBHs is related to the size of the causal horizon
collapsing, depending on the time of formation (see e.g. [11]). This implies
that

$m_{\rm PBH}\simeq\gamma M_{\rm H}$ (3.24)

where $M_{\rm H}$ is the mass of the cosmological horizon computed at horizon
crossing, while $\gamma$, varying roughly between $0.1$ and $1$, is the
efficiency parameter encapsulating the details of the gravitational collapse.

Given the fraction of PBHs $\beta_{\rm f}(M)$, one can compute the
contribution of PBHs to the energy density of the Universe and determine the
fraction of DM they represent today. In this sense it is convenient to compute
this quantity at the time of matter-radiation equality

$\Omega_{\rm PBH}^{\rm eq}=\int_{M_{\text{ev}}}^{M_{\rm eq}}\beta_{\rm
eq}(m_{\rm PBH}){\rm d}\ln m_{\rm PBH}$ (3.25)

where $M_{\text{ev}}$ is the lower bound due to Hawking evaporation, i.e. PBHs
with $M<M_{\text{ev}}$ would already been evaporated by now. We assume that
the fraction of PBH at the time of matter-radiation equality has grown because
of the cosmic expansion followed after the time of formation, which gives

$\beta_{\text{eq}}(m_{\rm PBH})=e^{(N_{\rm eq}-N_{\rm f})}\beta_{\rm f}(m_{\rm
PBH})\,.$ (3.26)

The PDF is characterized by the physics in the early Universe while the
threshold depends on the conditions at the time of formation. The analysis of
the gravitational collapse of curvature perturbations to form PBHs and the
appropriate threshold condition has been an active line of research in the
past years [49, 50, 366, 51, 367], with recent analysis on the dependence of
the threshold on the shape of the power spectrum [361, 368] and on the non-
linear relation between the curvature perturbations and the density contrast
[57, 60], which affects significantly the computation of the abundance with
respect to the standard calculation done with the linear Gaussian
approximation [363, 362].

Note that the threshold is also sensitive to the equation of state (EoS). For
example, the QCD phase transition is lowering the ratio between the pressure
and the energy density, reducing the threshold and therefore increasing the
production of PBHs formed during the QCD phase, characterized by a typical
mass of the order of $1M_{\odot}$ [20, 145].

Large curvature fluctuations can be produced by very different means in the
early Universe. If the curvature perturbations were Gaussian, then the PDF
would be described only by its variance $\sigma^{2}$ and the fraction of PBHs
measured at the time of formation is usually computed as

$\beta_{\rm f}(m_{\rm PBH})=2\int_{\delta_{\rm
c}}^{\infty}\frac{1}{\sqrt{2\pi}\sigma}e^{-\frac{\delta^{2}}{2\sigma^{2}}}{\rm
d}\delta=1-{\rm erf}\left(\frac{\delta_{\rm c}}{\sqrt{2}\sigma}\right)={\rm
erfc}\left(\frac{\delta_{\rm c}}{\sqrt{2}\sigma}\right)~{},$ (3.27)

where the factor $2$ outside the integral is sometimes introduced when the
abundance is computed using the excursion set method (see, for example, [369,
370]). This expression, however, does not take into account that when $\delta$
is larger then a certain value $\delta_{\rm max}$ (originally estimated
$\delta_{\rm max}\sim 1$ by Bernard Carr in 1975 [5]), the perturbation forms
a separate closed universe, topologically disconnected. A more accurate
version can nevertheless be obtained re-normalizing the previous expression as

$\beta_{\rm f}(m_{\rm PBH})=\frac{\displaystyle{\int_{\delta_{\rm
c}}^{\delta_{\rm
max}}\frac{1}{\sqrt{2\pi\sigma^{2}}}e^{-\frac{\delta^{2}}{2\sigma^{2}}}{\rm
d}\delta}}{\displaystyle{\int_{0}^{\delta_{\rm
max}}\frac{1}{\sqrt{2\pi\sigma^{2}}}e^{-\frac{\delta^{2}}{2\sigma^{2}}}{\rm
d}\delta}}=1-\frac{{\rm erf}\left(\displaystyle{\frac{\delta_{\rm
c}}{\sqrt{2}\,\sigma}}\right)}{{\rm erf}\left(\displaystyle{\frac{\delta_{\rm
max}}{\sqrt{2}\,\sigma}}\right)}=\frac{{\rm
erfc}\left(\displaystyle{\frac{\delta_{\rm c}}{\sqrt{2}\sigma}}\right)-{\rm
erfc}\left(\displaystyle{\frac{\delta_{\rm
max}}{\sqrt{2}\sigma}}\right)}{1-{\rm
erfc}\left(\displaystyle{\frac{\delta_{\rm max}}{\sqrt{2}\sigma}}\right)}\,.$
(3.28)

PBHs are produced through the tail of the PDF, ie. the area left under the PDF
curve where $\delta$ is larger than the critical threshold $\delta_{\rm c}$.
Therefore, PBH abundance is exponentially sensitive to the threshold and the
statistical properties of the primordial perturbations. Even tiny
modifications in the tail of the distribution can change the PBH abundance by
many orders of magnitude.

The variance of the field of density perturbations $\sigma$, according to the
Gaussian distribution of $\delta$, consistent with the linear approximation of
the relation between the energy density contrast and the curvature
perturbation (see next Subsection), is given by

$\sigma^{2}=\langle\delta^{2}\rangle=\int\limits_{0}^{\infty}\frac{\mathrm{d}k}{k}\mathcal{P}_{\delta}(k,r)=\frac{16}{81}\int\limits_{0}^{\infty}\frac{\mathrm{d}k}{k}(kr)^{4}\tilde{W}^{2}(k,r)T^{2}(k,r)\mathcal{P}_{\zeta}(k)\,,$
(3.29)

where ${\cal P}_{\delta}(k,r)$ and $\mathcal{P}_{\zeta}(k)$ are the density
and the curvature power spectra, while $\tilde{W}(k,r)$ is the Fourier
transform of the top-hat smoothing function and $T(k,r)$ is the linear
transfer function

$\displaystyle\tilde{W}(k,r)$
$\displaystyle=3\left[\frac{\mathrm{sin}(kr)-kr\mathrm{cos}(kr)}{(kr)^{3}}\right]$
(3.30) $\displaystyle T(k,r)$
$\displaystyle=3\left[\displaystyle{\frac{\mathrm{sin}\left(\frac{kr}{\sqrt{3}}\right)-\frac{kr}{\sqrt{3}}\mathrm{cos}\left(\frac{kr}{\sqrt{3}}\right)}{\left(\frac{kr}{\sqrt{3}}\right)^{3}}}\right]\,.$
(3.31)

The quantities $\delta_{\rm c}$ and $\sigma$, giving the critical threshold
$\nu_{\rm c}\equiv\delta_{\rm c}/\sigma$, are typically computed at super
horizon scale, when the curvature perturbation $\zeta$ for adiabatic
perturbations is time independent, making a linear extrapolation at horizon
crossing time, which is approximately the time when PBHs are formed (we will
discuss this later in more details). All this shows that a larger amplitude of
the power spectrum, combined with a shape decreasing the corresponding value
of $\nu_{\rm c}$, could increase significantly the fraction of PBHs because
the abundance of PBHs is exponentially sensitive to the value of $\delta_{\rm
c}$.

The linear transfer function666Ref. [364] has recently investigated the role
of non-linear corrections, coming from the non-linear radiation transfer
function and the determination of the true physical horizon crossing, on the
PBH abundance, showing that, while the critical threshold is larger than what
routinely assumed, the variance is unlikely to be changed with respect to the
linear one. ensures that the variance would always converge, because the
window function by itself, chosen consistently with the top-hat window
function in real space used in the definition of $\delta_{\rm c}$, does not
guarantee convergence if the power spectrum considered is sufficiently broad.
An alternative approach is to use a Gaussian window function, which provides
always convergence of the variance, but this is contaminating the value of
$\delta_{\rm c}$, introducing an error in the final computation of the
abundance [371].

###  3.5 The inevitable non-Gaussianity of the primordial black hole abundance

The non-linear relation between the curvature perturbation and the density
contrast in (3.4) induces an inescapable non-Gaussianity of the density
contrast, even assuming a Gaussian primordial initial condition [57, 363,
362]. This can be captured in the non-linear expression [363]

$\delta_{\rm m}=\delta_{\rm l}-\frac{3}{8}\delta_{\rm
l}^{2},\quad\,\,\,\,\delta_{\rm l}=-\frac{4}{3}r_{\rm m}\zeta^{\prime}(r_{\rm
m}).$ (3.32)

It highlights two key points. First of all, the probability of forming PBHs
does not depend on the comoving curvature perturbation itself $\zeta$, but on
its derivative $\zeta^{\prime}$. This is expected, given that on superhorizon
scales one can always add or subtract to the comoving curvature perturbation a
constant by a coordinate transformation and this may not influence any
physical result. Secondly, thanks to the conservation of the probability

$P(\delta_{\rm l})=P[\delta_{\rm m}(\delta_{\rm l})]\left|\frac{{\rm
d}\delta_{\rm m}}{{\rm d}\delta_{\rm l}}\right|,$ (3.33)

one can compute the abundance of PBHs simply using the Gaussian probability of
$\delta_{\rm l}$, integrating it from the critical amplitude for $\delta_{\rm
l}$

$\delta_{{\rm l},{\rm c}}=\frac{4}{3}\left(1-\sqrt{1-\frac{3}{2}\delta_{\rm
c}}\right),$ (3.34)

in terms of the critical amplitude $\delta_{\rm c}$.

###  3.6 Collapse of some initial non-Gaussian perturbations

It has been observed that nearly all models discussed in Sec. 2 generate a
considerable amount of primordial non-Gaussianity. We briefly summarize the
$\beta$ function in case of large local non-Gaussianity, which turns the PDF
into a $\chi^{2}$ distribution, explicitly $\delta_{\chi^{2}}=G^{2}-\langle
G^{2}\rangle$; and also in case of large cubic interactions, which turn the
PDF into a cubic-Gaussian distribution, explicitly $\delta_{G^{3}}=G^{3}$,
with $G$ a field obeying Gaussian statistics [372, 247, 373, 248].

$\displaystyle\beta_{{\rm f},\,\chi^{2}}(m_{\rm PBH})={\rm
erfc}\left(\sqrt{\frac{1}{2}+\frac{\delta_{\rm
c}}{\sqrt{2}\sigma}}\right)~{},$ $\displaystyle\beta_{{\rm f},\,G^{3}}(m_{\rm
PBH})={\rm erfc}\left(\left(\frac{\delta_{\rm
c}}{\sqrt{\frac{8}{15}}\sigma}\right)^{1/3}\right).$ (3.35)

We also note that the threshold density $\delta_{\rm c}$ depends slightly on
the non-Gaussianity [374, 288, 375].

###  3.7 Non-perturbative abundance from non-Gaussian perturbations

Refs. [376, 290] addressed the problem of computing the impact of a general
form of local primordial NG on PBH abundance beyond the perturbative approach.
In particular, let us focus on the functional form

$\displaystyle\zeta(\vec{x})=F(\zeta_{\text{\tiny G}}(\vec{x}))\,,$ (3.36)

where $F$ is a generic non-linear function of the Gaussian component
$\zeta_{\text{\tiny G}}$. There is a large number of relevant cases in the
literature of PBH formation where the NGs can be modelled as in Eq. (3.36). We
summarize few relevant examples in Tab. 1 (adapted from Ref. [376]).

Power-series expansion [377, 378, 379, 373, 361, 380, 381, 382, 383, 384, 375] | $\zeta=\zeta_{\text{\tiny G}}+\frac{3}{5}f_{\text{\tiny NL}}\zeta_{\text{\tiny G}}^{2}+\frac{9}{25}g_{\rm NL}\zeta_{\text{\tiny G}}^{3}+\dots$
---|---
Curvaton [385, 294] | $\zeta=\log\left[X(r_{\rm dec},\zeta_{\text{\tiny G}})\right]$
USR [386] | $\zeta=-\left(\frac{6}{5}f_{\text{\tiny NL}}\right)^{-1}\log\left(1-\frac{6}{5}f_{\text{\tiny NL}}\zeta_{\text{\tiny G}}\right)$
USR with an upward step [271] | $\zeta=-\frac{2}{|h|}\left[\sqrt{1-|h|\zeta_{\text{\tiny G}}}-1\right]$
Table 1:  Some examples of PBH formation scenarios characterised by local non-
Gaussianity of the form defined in Eq. (3.36). The function $X(r_{\rm
dec},\zeta_{\text{\tiny G}})$ given in the second row is defined in [385].

It is important to stress here that going beyond the perturbative approach
becomes a necessity at least when dealing with broad spectra of curvature
perturbations. This is because, as shown in Ref. [376], in these cases the
probability distribution function of $\zeta_{\text{\tiny G}}$ is evaluated
beyond the radius of convergence of the series expansion, leading to
inaccurate results if the series is truncated at finite order. Notice that in
non-attractor models the sign of non-Gaussianity is always positive [387],
thus enhancing the PBH abundance.

Using the relation between the curvature perturbation and the compaction
function in the gradient expansion (3.4), it is possible write

$\displaystyle\mathcal{C}(r)=-2\Phi\,r\,\zeta^{\prime}(r)\left[1+\frac{r}{2}\zeta^{\prime}(r)\right]=\mathcal{C}_{1}(r)-\frac{1}{4\Phi}\mathcal{C}_{1}(r)^{2}\,,~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\mathcal{C}_{1}(r)\equiv-2\Phi\,r\,\zeta^{\prime}(r)\,,$
(3.37)

where $\mathcal{C}_{1}(r)$ defines the so-called linear component of the
compaction function and we introduced the prefactor $\Phi$ that accounts for
potential thermal effects (see more in depth discussion below Eq. (3.46)).
Substituting in the previous expression the relation
$\zeta=F(\zeta_{\text{\tiny G}})$ (see eq. (3.36)), one obtains that the
linear component of the compaction function takes the form

$\displaystyle\mathcal{C}_{1}(r)=-2\Phi\,r\,\zeta_{\text{\tiny
G}}^{\prime}(r)\,\frac{dF}{d\zeta_{\text{\tiny G}}}=\mathcal{C}_{\text{\tiny
G}}(r)\,\frac{dF}{d\zeta_{\text{\tiny G}}}\,$ (3.38)

with $\mathcal{C}_{\text{\tiny G}}(r)\equiv-2\Phi\,r\,\zeta_{\text{\tiny
G}}^{\prime}(r)$. Consequently, the compaction function reads

$\displaystyle\mathcal{C}(r)=\mathcal{C}_{\text{\tiny
G}}(r)\,\frac{dF}{d\zeta_{\text{\tiny
G}}}-\frac{1}{4\Phi}\mathcal{C}^{2}_{\text{\tiny
G}}(r)\left(\frac{dF}{d\zeta_{\text{\tiny G}}}\right)^{2}\,.$ (3.39)

Crucially, $\mathcal{C}(r)$ depends on both the Gaussian linear component
$\mathcal{C}_{\text{\tiny G}}$ and the Gaussian curvature perturbation
$\zeta_{\text{\tiny G}}$. Both these random variables are Gaussian;
$\zeta_{\text{\tiny G}}$ is Gaussian by definition while
$\mathcal{C}_{\text{\tiny G}}$ is defined by means of the derivative of the
Gaussian variable $\zeta_{\text{\tiny G}}$.

Adopting threshold statistics on the compaction function, one can express the
non-Gaussian mass fraction as an integration of the joined PDF in the domain
which results in overthreshold perturbations, i.e.

$\displaystyle\beta_{\rm NG}$
$\displaystyle=\int_{\mathcal{D}}\mathcal{K}(\mathcal{C}-\mathcal{C}_{\rm
th})^{\gamma}\textrm{P}_{\text{\tiny G}}(\mathcal{C}_{\text{\tiny
G}},\zeta_{\text{\tiny G}}){\rm d}\mathcal{C}_{\text{\tiny G}}{\rm
d}\zeta_{\text{\tiny G}}\,,$ (3.40) $\displaystyle\mathcal{D}$
$\displaystyle=\left\\{\mathcal{C}_{\text{\tiny G}},\,\zeta_{\text{\tiny
G}}\in\mathbb{R}~{}:~{}~{}\mathcal{C}(\mathcal{C}_{\text{\tiny
G}},\zeta_{\text{\tiny G}})>\mathcal{C}_{\rm
th}~{}\land~{}\mathcal{C}_{1}(\mathcal{C}_{\text{\tiny G}},\zeta_{\text{\tiny
G}})<2\Phi\right\\}\,,$ (3.41)

where the multivariate normal distribution of $(\mathcal{C}_{\text{\tiny
G}},\zeta_{\text{\tiny G}})$ can be written as

$\textrm{P}_{\text{\tiny G}}(\mathcal{C}_{\text{\tiny G}},\zeta_{\text{\tiny
G}})=\frac{1}{(2\pi)\sigma_{c}\sigma_{r}\sqrt{1-\gamma_{cr}^{2}}}\exp\left(-\frac{\zeta_{\text{\tiny
G}}^{2}}{2\sigma_{r}^{2}}\right)\exp\left[-\frac{1}{2(1-\gamma_{cr}^{2})}\left(\frac{\mathcal{C}_{\text{\tiny
G}}}{\sigma_{c}}-\frac{\gamma_{cr}\zeta_{\text{\tiny
G}}}{\sigma_{r}}\right)^{2}\right]\,.$ (3.42)

To shorten the notation, we have followed Ref. [388] and defined the element
of the covariance matrix as

$\displaystyle\langle\mathcal{C}_{\rm G}\mathcal{C}_{\rm G}\rangle$