Primordial black hole constraints with Hawking radiation - a review

Jérémy<EMAIL_ADDRESS>
aUniv. Lyon, Univ. Claude Bernard Lyon 1, CNRS/IN2P3, IP2I Lyon, UMR 5822,
F-69622, Villeurbanne, France

Primordial black holes are under intense scrutiny since the detection of
gravitational waves from mergers of solar-mass black holes in 2015. More
recently, the development of numerical tools and the precision observational
data have rekindled the effort to constrain the black hole abundance in the
lower mass range, that is $M<10^{23}$g. In particular, primordial black holes
of asteroid mass $M\sim 10^{17}-10^{23}\,$g may represent 100% of dark matter.
While the microlensing and stellar disruption constraints on their abundance
have been relieved, Hawking radiation of these black holes seems to be the
only detection (and constraining) mean. Hawking radiation constraints on
primordial black holes date back to the first papers by Hawking. Black holes
evaporating in the early universe may have generated the baryon asymmetry,
modified big bang nucleosynthesis, distorted the cosmic microwave background,
or produced cosmological backgrounds of stable particles such as photons and
neutrinos. At the end of their lifetime, exploding primordial black holes
would produce high energy cosmic rays that would provide invaluable access to
the physics at energies up to the Planck scale. In this review, we describe
the main principles of Hawking radiation, which lie at the border of general
relativity, quantum mechanics and statistical physics. We then present an up-
to-date status of the different constraints on primordial black holes that
rely on the evaporation phenomenon, and give, where relevant, prospects for
future work. In particular, non-standard black holes and emission of beyond
the Standard Model degrees of freedom is currently a hot subject.

###### Contents

  1. 1 Introduction
  2. 2 Primordial black holes
    1. 2.1 Historical overview
    2. 2.2 Formation mechanisms
      1. 2.2.1 Extended mass distribution
      2. 2.2.2 Spin distribution
  3. 3 Hawking radiation
    1. 3.1 Historical overview
      1. 3.1.1 Thermodynamics aspects
      2. 3.1.2 Quantum mechanics aspects
      3. 3.1.3 Hawking papers
    2. 3.2 Theoretical aspects
    3. 3.3 Basic formulas
    4. 3.4 Numerical recipes
  4. 4 Non-evaporated primordial black holes
    1. 4.1 Lensing constraints
    2. 4.2 Dynamical constraints
    3. 4.3 Accretion constraints
    4. 4.4 Gravitational wave constraints
  5. 5 Baryogenesis
    1. 5.1 Hawking model
    2. 5.2 Weinberg model
    3. 5.3 Dolgov model
    4. 5.4 Nagatani model
    5. 5.5 Fujita model
    6. 5.6 Other ideas and prospects
  6. 6 Big bang nucleosynthesis
    1. 6.1 Nucleon inter-conversion
    2. 6.2 Hadro-dissociation
    3. 6.3 Photo-dissociation
    4. 6.4 Lithium problem
    5. 6.5 Modern studies
  7. 7 Cosmic microwave background
    1. 7.1 CMB distortions
    2. 7.2 Reionization
    3. 7.3 CMB anisotropies
    4. 7.4 21 cm signal
  8. 8 Cosmological backgrounds
    1. 8.1 Photons
    2. 8.2 Neutrinos
    3. 8.3 Gravitons
    4. 8.4 Present backgrounds
  9. 9 Cosmic rays
    1. 9.1 Electrons & positrons
      1. 9.1.1 Direct detection
      2. 9.1.2 Indirect detection
      3. 9.1.3 511 keV line
    2. 9.2 Antimatter
  10. 10 Final bursts
    1. 10.1 Photons
    2. 10.2 Other particles
    3. 10.3 Direct observation
  11. 11 Beyond monochromatic and Schwarzschild constraints
    1. 11.1 Extended mass functions
    2. 11.2 Spinning black holes
    3. 11.3 Non-standard black holes
      1. 11.3.1 Higher-dimensional black holes
      2. 11.3.2 Other metrics
      3. 11.3.3 Black hole remnants
    4. 11.4 Non-standard particles
  12. 12 Conclusion

##  1 Introduction

Dark matter (DM) is currently one of the deepest problems of the cosmological
standard model. Indeed, whereas the precision era of cosmology has probed the
bases of this model, many unanswered questions remain. Even if the concept of
DM—an invisible matter which nevertheless has visible gravitational effects—is
well established, there are numerous solutions to explain its nature (see e.g.
[1] for a recent review and [2] for an historical account). Primordial black
holes (PBHs) were one of the candidates proposed to explain DM, in the early
1970’s. PBHs are black holes (BHs) that formed in the early universe from the
collapse of overdensities resulting from quantum fluctuations, and are
predicted also by other mechanisms such as phase transitions.

PBHs are not the result of stellar collapse, thus their mass can span a very
wide range of masses, in fact from the Planck mass to the mass of the universe
enclosed by the Hubble horizon today. On the other hand, their abundance is
tightly constrained by the numerous effects they have on the cosmological
history. Heavy PBHs, with mass $M>10^{23}\,$g, may act as gravitational
lenses, disrupt binaries, accrete matter or coalesce. This gives corresponding
microlensing, stellar surveys, X-ray observations and gravitational wave (GW)
constraints. Overall, PBHs with $M>10^{23}\,$g cannot represent 100% of DM and
the constraints will certainly tighten with the upcoming GW instruments (LISA,
Einstein telescope) or microlensing surveys (Roman, Rubin). For a recent
review of the PBH constraints as a DM candidate, see [3, 4], for a review of
the general prospects on PBH DM, we refer the reader to the Snowmass white
paper [5], and for a dedicated review of GW perspectives see [6].

The lightest PBHs are more interesting as their abundance is yet unconstrained
in the window $10^{17}-10^{23}\,$g, so that they could constitute 100% of DM.
These PBHs are light enough so that the phenomenon of Hawking radiation (HR)
becomes important. HR is a semi-classical process predicted by Hawking, that
makes BHs lose their mass by emitting all the quantum fields with a quasi-
thermal spectrum. Theorized in the 1970’s, this effect remains unobserved,
mainly because it is so faint for heavy BHs: the temperature of the Hawking
spectrum is inversely proportional to the BH mass. PBHs in the “HR window”
$10^{17}-10^{23}\,$g however would radiate huge amounts of particles of all
types, which could be the only way of assessing their abundance. The “HR
window” coincides with MeV astronomy, a quite unexplored domain that has
become the stake of many instrument proposals (AMEGO, ASTROGRAM) [7].
Prediction of the PBH signals in this mass window would then allow to close
the last possibility for PBH DM [5].

Finally, PBHs with mass $M<5\times 10^{15}\,$g would have already evaporated
by now, hence they would not participate in the DM density. While their final
explosion is searched for in $\gamma$-ray and even neutrino experiments [8],
more promising prospects come from the signature they may have left in the
past cosmological eras. The smallest PBHs, with mass $M<10^{9}\,$g, would have
evaporated before big bang nucleosynthesis (BBN). They may have produced the
baryon asymmetry or generated the particle DM abundance: BHs evaporate in all
(beyond) Standard Model ((B)SM) particles. Slightly heavier PBHs with mass
$M\sim 10^{9}-10^{12}\,$g evaporate during BBN, and their energetic
evaporation products enter the chain of nuclear reactions that created the
light elements; their abundance is then constrained by the successful
predictions of standard BBN. PBHs with $M\sim 10^{12}-10^{17}\,$g modify the
recombination history of the universe. As the cosmic microwave background
(CMB) is one of the better known cosmological probe, with precision data from
the Planck satellite and prospective ones (CMB-S4, PIXIE), any exotic energy
injection is limited by observation. In particular, the EDGES collaboration
has cleared the way for surveys of the Dark Ages after recombination,
observing the $21\,$cm sky for the first time. This may result in the
strongest constraints on PBHs up to $M\sim 10^{18}\,$g, with future
instruments such as SKA. PBHs emit particles during their whole lifetime, and
the spectrum of those backgrounds today encodes some of the thermal history of
the universe, on top of the genuine characteristics of HR. The photon
background, the easiest to observe, goes back to the CMB ($t\sim
380\,000\,$yr); the neutrino background would have been emitted since neutrino
decoupling ($t\sim 1\,$s) and would then constitute a complementary probe of
the early universe with BBN; the putative graviton spectrum would have
propagated unperturbed since inflation, which would represent the ultimate
probe of the first instants of the universe. Most of these topics are reviewed
in [9], but the constraints evolve constantly and no review focusing on HR
limits has yet been published.

With the development of precision observations, it further becomes more and
more important to predict the precise shape of the PBH signal, which is the
aim of several recent numerical tools such as BlackHawk [10, 142], CosmoLED
[11], nuHawkHunter [12] and ULYSSES [13, 14]. These tools rely on particle
physics codes that compute the interactions of the Hawking radiated particles
from the MeV to the Planck scale. Moreover, intense work is performed to
obtain the signals (and constraints) on non-standard (P)BHs, such as regular
solutions to the Einstein equations, or BSM particles, such as interacting DM.
Realistic PBH formation mechanisms are now taken into account when
constraining their abundance, with e.g. extended mass distributions or
spinning BHs.

In this review, we present first the historical development of the PBH (as DM)
(Section 2) and HR (Section 3) ideas, with a recall of the fundamental
equations required to compute the HR spectra. We briefly review the
constraints on (heavy) PBHs not linked to evaporation (Section 4). Then, and
this is the main focus of the present paper, we give a complete account of the
different constraints, in cosmological chronological order: baryogenesis
(Section 5), BBN (Section 6), CMB (Section 7), diffuse backgrounds of photons,
neutrinos and gravitons (Section 8), cosmic rays such as electrons/positrons
and antiprotons (CRs, Section 9), and the final burst of PBHs (Section 10). In
the last part (Section 11), we review the recent developments regarding
extended mass functions, spinning PBHs and non-standard BHs/BSM degrees of
freedom. We finally conclude with general prospects (Section 12). A master
figure summarizing all the constraints is given in Fig. 14. Throughout the
review, we use the natural system of units defined as $G=\hbar=k_{\rm
B}=c=4\pi\epsilon_{0}=1$.

##  2 Primordial black holes

PBHs are BHs formed in the early universe. In this Section, we first go
through a brief historical overview of the PBH ideas, and then sketch the main
formation mechanisms and mass ranges expected for PBHs, with particular focus
on formation during a dust phase, and extended mass and angular momentum
distributions.

###  2.1 Historical overview

Using only Newtonian mechanics, Michell and Laplace independently predicted at
the end of the XVIIth century that if there were an object sufficiently dense
that the escape velocity at its surface was faster than the speed of light, it
would be totally black and absorb anything that would come gravitationally
bound to it. That is precisely what a “black hole” would do. Schwarzschild was
the first one to derive the solution of the spherically symmetric and static
Einstein equations around a pointlike mass $M$ in 1916 [15] (hereafter SBH,
for Schwarzschild BH), right after Einstein published the general relativity
framework

${\rm d}s^{2}=-G(r){\rm d}t^{2}+\frac{1}{F(r)}{\rm d}r^{2}+H(r){\rm
d}\Omega\,,$ (1)

with

$\left\\{\begin{array}[]{ll}F(r)=G(r)=1-\dfrac{r_{\rm S}}{r}\,,\\\
H(r)=r^{2}\,,\end{array}\right.$ (2)

and $r_{\rm S}\equiv 2M$ is the Schwarzschild radius, illustrated in Fig. 1
(left panel). At space infinity, that is $r\gg r_{\rm S}$, the metric becomes
asymptotically flat, and reduces to the Minkowsky metric

${\rm d}s^{2}=-{\rm d}t^{2}+{\rm d}r^{2}+r^{2}{\rm d}\Omega^{2}\,,$ (3)

which is the empty and flat solution to the Einstein equations. There is a
singularity of the metric at $r=0$, corresponding to a real curvature
singularity as the Ricci scalar $R$ goes to infinity, and a fake singularity
($R$ is finite) at $r=r_{\rm S}$, which can be resolved by a change of
coordinates. The singularity is an intricate problem, but the general
relativity calculations show that the worldlines of free-falling objects
toward the singularity exchange time and space characteristics precisely at
$r=r_{\rm S}$, denoted as the BH horizon of events. Hence, any signal emitted
by a probe that would have crossed the horizon would never reach an observer
at infinity but remain enclosed in the BH. The interior of the BH horizon is
causally separated from the exterior by the event horizon, which effectively
hides away the metric singularity and satisfies the “cosmic censorship
conjecture” [16] (CCC). Interestingly, “invisible collapsed objects”, (i.e.
BHs), were proposed very early as a potential candidate for DM [17, 18].

The Maxwell equations were included in that framework to obtain the metric
around a pointlike mass $M$ with electric charge $Q$ by Reissner & Nordström
[19, 20] (hereafter RNBH, for Reissner-Nordström BH) with the metric
coefficients

$\left\\{\begin{array}[]{ll}F(r)=G(r)=1-\dfrac{r_{\rm S}}{r}+\dfrac{r_{\rm
Q}^{2}}{r^{2}}\,,\\\ H(r)=r^{2}\,,\end{array}\right.$ (4)

where $r_{Q}\equiv Q^{2}$. This metric exhibits one real singularity at $r=0$
and two fake coordinate singularities at

$r_{\pm}\equiv r_{\rm S}\dfrac{1\pm\sqrt{1-Q^{*2}}}{2}\,,$ (5)

where $r_{-}$ is a Cauchy horizon and $r_{+}$ is the BH event horizon.
$Q^{*}\equiv Q/M$ is the reduced (dimensionless) BH charge.

The case of a rotating BH with angular momentum $a\equiv J/M$, which is only
axisymmetric, is more mathematically involved and was solved by Kerr in 1963
[21] (hereafter KBH, for Kerr BH). The metric is here

${\rm d}s^{2}=\big{(}{\rm d}t-a\sin^{2}\theta{\rm
d}\phi\big{)}^{2}\,\frac{\Delta}{\Sigma}-\left(\frac{{\rm
d}r^{2}}{\Delta}+{\rm d}\theta^{2}\right)\Sigma-\big{(}(r^{2}+a^{2}){\rm
d}\phi-a{\rm d}t\big{)}^{2}\,\frac{\sin^{2}\theta}{\Sigma}\,,$ (6)

where $\Sigma(r)\equiv r^{2}+a^{2}\cos^{2}\theta$ and $\Delta(r)\equiv
r^{2}-2Mr+a^{2}$. The Cauchy and event horizons are located at

$r_{\pm}\equiv r_{\rm S}\dfrac{1\pm\sqrt{1-a^{*2}}}{2}\,,$ (7)

where $a^{*}\equiv a/M$ is the reduced (dimensionless) angular momentum or “BH
spin”.

The combination of these two solutions gives the Kerr–Newman metric around a
charged, rotating mass [22, 23], with horizons located at $r_{\pm}=r_{\rm
S}(1\pm\sqrt{1-Q^{*2}-a^{*2}})/2$. A comparison of the radii of a SBH with
those of a KBH or RNBH for different values of $a^{*}$ and $Q^{*}$ are given
in Fig. 1 (right panel).

One immediately remarks that when $a^{*},Q^{*}\rightarrow 0$, the SBH is
recovered, while in the opposite limit $a^{*},Q^{*}\rightarrow 1$, the Cauchy
and event horizons radii collapse to the value $r_{\rm S}/2$. In this case,
denoted as “extremal”, the BH is ill-defined as the coordinate singularity at
the BH center $r=0$ is not hidden behind a horizon—the CCC is violated; this
limit is shown as a black vertical line on the right panel of Fig. 1.

Figure 1: Left: Radius of a SBH as a function of its mass. The vertical solid
line at $M=M_{\rm Pl}$ represents the classical limit, the vertical dashed
line at $M\sim 0.1\,$g represents the inflation limit, the vertical dashed
line at $M\sim 10^{38}\,$g represents the BBN limit and the vertical dotted
line at $M\sim 10^{50}\,$g represents the “incredulity limit” (see text).
Right: Event (solid) and Cauchy (dashed) horizons radii $r_{\pm}$ normalized
to the Schwarzschild radius (dotted) $r_{\rm S}$ for a KBH ($x=a^{*}$) and
RNBH ($x=Q^{*}$). The vertical line at $x=0$ represents the Schwarzschild
limit and the vertical solid line at $x=1$ represents the extremal limit (see
text).

BHs result in general from the collapse of some energy density inside its own
Schwarzschild radius. At the end of its evolution, a star shrinks because most
hydrogen has burned into helium and the gravitational forces win against the
failing nuclear reactions radiation pressure. For a certain parameter space of
stellar conditions, the star can collapse into a “stellar” BH (sometimes
denoted as “astrophysical” BH). Tolman, Oppenheimer & Volkoff [24, 25] (TOV)
identified an inferior mass for these BHs, with modern studies giving $M_{\rm
TOV}\sim 2\,M_{\odot}$ [26]. Upper limits on their mass are more dubious but
it is generally not expected that genuine stellar mass BHs of more than $\sim
100\,M_{\odot}$ can form.

Nowhere else in the present universe is there matter denser than the interior
of stars, where BHs could form. However, the universe was order of magnitudes
denser in the past. Then, it is not unimaginable that in the early universe,
patches of the universe which were overdense relatively to the average density
due to some statistical fluctuation collapsed into early non-stellar BHs.
These are generally denoted as “primordial black holes”. The very first
mention of PBHs dates back to 1966: Zel’dovich & Novikov, in a paper entitled
“The hypothesis of cores retarded during expansion and the hot cosmological
model” [27], speak about “retarded cores” [28] whose evolution, because of
gravitational collapse, would escape the overall expansion of the universe,
and describe these cores to have radius $R<2M$ which is precisely the radius
of a SBH.

This original idea seems to have gone quite unnoticed until Hawking published
the paper “Gravitationally collapsed objects of very low mass” [29]. The idea
is fundamentally the same as Zel’dovich & Novikov, i.e. large perturbations
collapsing after inflation inside their Schwarzschild radius to give BHs
(explicitly named here). The collapse is described as a classical process,
hence a lower limit is imposed on the mass of the BH formed $M>M_{\rm Pl}\sim
2\times 10^{-5}\,$g for general relativity to be reliable. Hawking concludes
that paper by stating that:

> […] it is tempting to suppose that the major part of the mass of the
> Universe is in the form of collapsed objects. This extra density could
> stabilize clusters of galaxies which, otherwise, appear mostly not to be
> gravitationally bound. [29]

Some years later, Carr & Hawking [30] show that the concerns raised by
Zel’dovich & Novikov [27] about catastrophic early PBH accretion are not
established, however:

> The obvious place to look for such giant black holes [that would have
> accreted rapidly] would be in clusters of galaxies where they might provide
> the missing mass necessary to bind clusters gravitationally. [30]

Meszaros [31, 32] finally argues that galaxy formation could benefit from the
existence of PBHs, the latter seeding early structures and giving the “hidden
mass in clusters required by virial theorem arguments.” These claims, together
with that of [17, 18] concerning astrophysical BHs, constitute the prelude of
the “DM = PBHs” scenario.

###  2.2 Formation mechanisms

Carr & Hawking [30] first quantitatively described the process of PBH
formation, which was numerically computed in [33] (for a recent review, see
[34]). A region of locally enhanced density collapses if its radius is smaller
than the associated Jeans length. The collapse must supplant the pressure
forces (which depend on the equation of state of the cosmological fluid) and
the universe expansion. BH formation then depends on the amplitude of the
perturbation, which must be superior than some threshold
$\Delta\rho/\rho>\delta$.

At time $t$, the collapsing region must be of order of the particle horizon
size $M_{\rm hor}\sim H(t)^{-1}$ for causal reasons. This is the standard
scenario, which relies on the assumption that there was a spectrum of
primordial inhomogeneities. Modern calculations show that a PBH formed during
radiation domination (RD) has an initial mass $M_{\rm i}$ linked to its
formation time $t_{\rm f}$ by [9]

$M_{\rm i}(t_{\rm f})=\dfrac{\gamma}{2H(t_{\rm f})}\sim
10^{15}\,\text{g}\left(\dfrac{t_{\rm f}}{10^{-23}\,{\rm s}}\right),$ (8)

where $\gamma$ is a parameter linked to the collapse mechanism. For spherical
collapse during RD, $\gamma\sim 0.2$ [9]. In the usual model of PBH formation,
PBHs form after the inflation period, hence there is an inferior limit on the
size of the particle horizon derived from CMB observations $H(t_{\rm f})^{-1}$
[35] translating into an inferior limit on the mass of PBHs

$M_{\rm min}\sim 0.1\,\text{g}\,.$ (9)

However, this constraint applies only on conventional inflationary scenarios,
_i.e._ the standard slow-roll models of inflation with Einstein gravity (see
_e.g._ the review [36] and references therein). In more sophisticated
scenarios, the scale of inflation can not be determined by CMB observations.
An absolute lower limit on the PBH mass is the Planck mass $M_{\rm Pl}$ at
which quantum gravity effects should be sizeable and our understanding of
general relativity breaks down. PBHs with size smaller than the Planck length
$\ell_{\rm Pl}\sim 10^{-35}\,$m would have formed at times before the Planck
time $t_{\rm Pl}\sim 10^{-43}\,$s which is forbidden by the trans-Planckian
conjecture [37].

An upper limit on the mass of PBHs at formation is derived from the fact that
they should form before the onset of BBN not to spoil the baryon-to-photon
ratio

$M_{\rm BBN}\sim M(t_{\rm f}\approx 1\,\text{s})\sim 10^{5}\,M_{\odot}\,.$
(10)

There are also trivial “incredulity” reasons to believe that there is no giant
PBH of $M\gtrsim M_{\rm incr}\sim 10^{50}\,$g occupying most of the Hubble
sphere today, first because there are other complex structures in the sky,
second because if there were a giant PBH, we would be falling onto it due to
its gravitational attraction and would see a strong Doppler dipole feature in
the CMB [38, 39]. The limits $M_{\rm Pl}$, $M_{\rm min}$, $M_{\rm BBN}$ and
$M_{\rm incr}$ are represented by vertical lines on the left panel of Fig. 1.

The list of PBH formation scenarios is extensive, and we refer the reader to
the recent reviews [9, 4] for a complete description. They include:

  * 1.

collapse of primordial overdensities (the standard scenario);222This scenario
seems deprecated due to the tremendous increase of the fluctuation spectrum
needed at small scales compared to the CMB, see the very explicit Fig. 1 of
[4]. This has motivated more exotic formation scenarios, and refutes the usual
claim that PBH DM is compelling because it does not require new physics.

  * 2.

collapse of inflaton fields;

  * 3.

collapse of topological defects;

  * 4.

collapse from bubble collision due to an early first order phase transition.

In particular, some formation mechanisms (e.g. grand unified theory—GUT—phase
transition) explain both PBH formation and later cosmological events (e.g. PBH
baryogenesis mediated by GUT bosons [40]). In general, an alternative
cosmological model results in very different PBH constraints [41]. Let us
detail some formation mechanisms related to EMDE, or which result in extended
PBH mass or spin distributions.

For example, formation of PBHs during an early matter-dominated era (EMDE) was
reviewed in detail in [42]. It is shown that PBH formation is easier in this
case due to reduced pressure forces. The relationship between a PBH mass and
its formation time, and thus the PBH constraints then depend on the EMDE
parameters and are in general weaker than the corresponding RD constraints
[43]. Modern studies of PBH formation during an EMDE have been performed by
[44, 45, 46, 47, 48, 49, 50].

####  2.2.1 Extended mass distribution

One prominent feature of the first review of PBH constraints by Carr [51] is
that they use the Press-Schechter formalism designed for scale-invariant
perturbations in the early universe and the formation of structures. Carr
predicts that the mass spectrum of PBHs is a power-law with exponent ${\rm
d}n\propto M^{-f(w)}{\rm d}M$ where $f(w)$ is some function of the equation of
state parameter $P=w\rho$ and where the collapse into a PBH happens once the
density perturbation exceeds some threshold $\Delta\rho/\rho>\delta$. In the
RD era, one obtains $f(w)=5/2$. This exponent should be modified for PBH
formation during an EMDE.

Hence, the mass spectrum of PBHs was originally believed to be extended and to
span a rather large range of masses (from $M_{\rm Pl}$ upward). This link
between primordial scale-invariant fluctuations of density and the PBH
abundance shows that:

> PBHs are unique since they alone could be expected to survive the
> dissipative effects which erase all other imprints of conditions in the
> first second of the universe. [51]

Hence, PBH constraints are sometimes converted into primordial inhomogeneities
(see e.g. Fig. 19 of [9]). The monochromatic distribution of PBHs often used
to derive abundance constraints is then only a convenient approximation to a
more realistic extended mass distribution. This poses the mathematical
question of how to convert the constraints obtained for monochromatic mass
distributions to extended ones?

Yokoyama [52] is the first, to our knowledge, to try and use an analytical
procedure of conversion. This follows the discovery of the critical behaviour
of the PBH collapse process by Niemeyer & Jedamzik [53, 54]. The critical
collapse mechanism leads to the formation of PBHs with mass $M<M_{\rm hor}$;
hence at each time, a distribution of PBHs masses would arise. The total
distribution would then be the sum of these instantaneous formations. The
rough computation by Yokoyama shows that using an extended mass function may
broaden the excluded PBH parameter space.

At about the same time, calculations showed that a monochromatic peak in the
primordial fluctuation power spectrum results in fact in a log-normal PBH mass
distribution, centered around the characteristic horizon scale at formation
[55]. This distribution was further refined using peak theory [56, 57]. In
fact, a Gaussian log-normal function can mimic any peak in the PBH
distribution resulting from a particular mechanism of formation [58]. The most
common PBH mass distributions are described e.g. in [59].

Concerning conversion methods, pioneering analytical and numerical work has
been done in [60, 61, 59, 62]. Refs. [60, 61] propose to use a sufficiently
fine binning procedure to separate the extended distribution into independent
Dirac functions, and Ref. [62] implicitly computes the equivalent
monochromatic amplitude resulting in the same constraint as the extended
distribution at some mass $M$. The most commonly used method is that of Carr
et al. [59] (hereafter the “Carr method”) which decomposes any measure as an
infinite sum over correlated PBH contributions from different masses; the
difficulty is to find the $n$-th order correlation functions. These works
confirm that the constraints on an extended distribution should be more
stringent than the expected constraint resulting from the addition of
monochromatic ones. The precision of the conversion has been examined in some
papers, showing mitigated results [63]. For HR constraints, we advocate the
use of a numerical procedure such as BlackHawk which can directly compute the
full HR spectrum for any PBH distribution before applying the constraints.

In the following, PBH mass distributions at formation are denoted as ${\rm
d}n/{\rm d}M$ with the total number of PBHs $N_{\rm PBH}$ and their density
$\rho_{\rm PBH}$ given by

$N_{\rm PBH}=\int_{M_{\rm min}}^{M_{\rm max}}\dfrac{{\rm d}n}{{\rm d}M}\,{\rm
d}M\,,\quad\text{and}\quad\rho_{\rm PBH}=\int_{M_{\rm min}}^{M_{\rm
max}}M\,\dfrac{{\rm d}n}{{\rm d}M}\,{\rm d}M\,.$ (11)

####  2.2.2 Spin distribution

PBHs formed during RD are believed to have negligible spin [64], which has
also been proven for the critical collapse mechanism [65]. On the other hand,
PBHs formed during an EMDE could have sizeable to near-extremal spin [45, 66,
67, 68]. In any case, some of the PBHs should form with a significant spin
[69], and PBHs can also be spun up either through early accretion processes
[70] or through hierarchical mergers [71].

The problem is then the same as for extended mass distributions. The
constraints on SBHs could be converted to KBHs using some analytical
procedure; or realistic extended spin distributions could be reduced to their
peak/average value (e.g. [72]). We instead advocate for the use of the
BlackHawk capacity to simulate extended distributions of both mass and spin to
compute accurate constraints, as was done for the first time in [73].

In the following, the distribution of PBH secondary parameter $x$ ($x$ can be
the PBH spin $a^{*}$, electric charge $Q^{*}$, etc.) is denoted as ${\rm
d}\tilde{n}/{\rm d}x$ and is normalized as

$\int_{x_{\rm min}}^{x_{\rm max}}\dfrac{{\rm d}\tilde{n}}{{\rm d}x}\,{\rm
d}x=1\,.$ (12)

A joint distribution may be denoted as ${\rm d}^{2}n/{\rm d}M{\rm d}x$.

##  3 Hawking radiation

In this Section, we review the physics of Hawking radiation, with historical
ideas and papers. We describe the thermodynamics and quantum mechanics bases,
then detail the HR calculation and finally give the basic formulas, with
analytical and numerical results.

###  3.1 Historical overview

####  3.1.1 Thermodynamics aspects

In the 1970’s, some work was done in order to conciliate BH mechanics with
known theories of thermodynamics. Indeed, BHs seemed to violate the 2nd law of
thermodynamics as they can swallow a great amount of information with only
modification of mass $M$, charge $Q$ and angular momentum $J$. As BHs are
described only by these three quantities, as stated by the “no-hair” theorem
set out by Israel [74, 75, 76], there is a loss of entropy known as the
“information paradox” [77]. Thus, Bekenstein [78] searched for a generalized
second law for BHs, and associated their area to an effective entropy that
should not decrease during BH evolution. This claim was further supported by
an analog version of the first law of thermodynamics

${\rm d}M=\Theta{\rm d}A+\vec{\Omega}\cdot\vec{{\rm d}L}+\Phi{\rm d}Q\,,$ (13)

where $\vec{{\rm d}J}$ (resp. ${\rm d}Q$) is the change in angular momentum
(resp. charge) of the BH while $\vec{\Omega}$ (resp. $\Phi$) plays a role
analog to angular frequency (resp. electric potential). Thus, the term
$\Theta$ in front of the entropy change would be identified as an effective
temperature whose expression depends on the surface gravity $\kappa$ of the
BH.

Bardeen, Carter & Hawking [79] go one step further by giving explicitly the 4
laws of BH mechanics:

  * 0th law: the temperature of a BH $T=\kappa/2\pi$ is constant over its event horizon;

  * 1st law: the infinitesimal evolution of the BH mass is given by ${\rm d}M=T{\rm d}A/4+\vec{\Omega}\cdot\vec{{\rm d}J}+\Phi{\rm d}Q$;

  * 2nd law: the total entropy composed of that of a black hole $S=A/4$ plus that of the rest of the universe can only increase;

  * 3rd law: the temperature of a BH cannot be decreased to absolute 0 by a finite number of operations.

In this paper, the entropy is found to be $S=A/4$ in the first law, so that
the effective temperature is identified as

$T=\dfrac{\kappa}{2\pi}\,.$ (14)

The third law flows from the CCC. Then follows some noticeable statement in
adequacy with the views at that time:

> In fact the temperature of a black hole is absolute zero [because it] cannot
> be in equilibrium with blackbody radiation at any non-zero temperature. [79]

precisely because a BH can only accrete matter and radiation. For a complete
modern review of the thermodynamics of BHs, we refer the reader to [80, 81].

####  3.1.2 Quantum mechanics aspects

At about the same epoch, theoretical work was pursued in order to check
whether BHs were stable objects regarding quantum mechanics, that is, if BHs
could develop diverging perturbations that would challenge their survival.
This was particularly interesting for KBHs as these can experience stimulated
emission of bosonic fields through “superradiance” effects [82], which would
efficiently extract mass and angular momentum from the black hole and could
lead to “black hole bombs” [83].333For a complete review of the superradiance
effect we refer the interested reader to [84]. Teukolsky & Press published a
series of 3 papers [85, 86, 87] in which they study the stability of the Kerr
metric against a spin $s$ equal to 1, $\nicefrac{{1}}{{2}}$ and 2 wave
scattering. The separated equations for the scalar case were derived earlier
by Carter [88].444This problem was also explored by Starobinskii & Churilov
for bosonic fields in [89, 90] and Unruh for fermionic field in [91]. They
obtained the fundamental separated radial and angular equations for the
propagation of a wave of spin $s$ and energy $E$ outside a BH of mass $M$, in
semi-classical general relativity (i.e. classical metric and quantum fields).
The radial equation is therefore named the “Teukolsky equation” and the
angular part is solved by spin-weighted spheroidal harmonics [92]. Resolution
of the Teukolsky equation with the correct boundary conditions yields the
transmission (absorption) and reflection coefficients of a wave over a BH
horizon and thus the superradiant amplification of incident bosonic waves,
while fermions are not superradiant.

####  3.1.3 Hawking papers

Confronted with the possible existence of PBHs with size and mass way below
the classical stellar processes—in fact in the quantum mechanics
regime—Hawking wondered how quantum mechanics effects would come into play.
This led to the famous discovery that BHs are not black but radiate a steady
flux of particles like thermal bodies. This was first proposed in a paper
published in Nature in 1974 entitled “Black hole explosions?” [93]. In this
paper, Hawking claims that:

> It seems that any black hole will create and emit particles such as
> neutrinos or photons at just the rate as one would expect if the black hole
> was a body with a temperature of [$T=\kappa/2\pi$]. [93]

which is precisely the temperature (14) encountered when reconciling BH
mechanics with thermodynamics. Hawking immediately suspected a deep
fundamental link between this temperature (blackbody emission) and the
effective temperature in BH thermodynamics (entropy). The link was identified
in the context of the microcanonical ensemble in [94, 95]. Follows:

> As a black hole emits this thermal radiation one would expect it to lose
> mass. This in turn would increase the surface gravity [$\kappa\propto 1/M$]
> and so increase the rate of emission. [93]

Hence, any BH would have a finite lifetime that can be estimated if its
luminosity follows the blackbody law $L\sim AT^{4}$. Straightforward
integration shows that a PBH with the critical mass $M_{*}\sim 10^{15}\,$g
would have evaporated by today if formed just after inflation.555Given the
roughness of this estimation and the uncertainty concerning the age of the
universe at that time, it is remarkable that the estimation by Hawking falls
within 50% of the modern value $M_{*}=5\times 10^{14}\,$g. Due to the inverse
power-law dependency of the emissivity with the mass, a very powerful
“explosion” is to be expected at the end of the BH evolution.

The details of HR are given in a subsequent longer paper from 1975 entitled
“Particle creation by black holes” [96]. The calculations are based on quantum
mechanics in curved spacetime. We don’t give them here but refer the reader to
[97] which follows step by step the original Hawking derivation in a
comprehensive way. Basically, the flux of particle originates in the fact that
an observer close to the BH is freely falling with an acceleration given by
the surface gravity $\kappa$, while an observer far away from the BH horizon
is at rest in the BH frame. Thus, they define different local bases and vacuum
quantum states and the conversion between these bases results in a net thermal
flux of particles at infinity with rate

$\dfrac{{\rm d}^{2}N}{{\rm d}t{\rm
d}E}=\dfrac{1}{2\pi}\dfrac{\Gamma_{Eslm}}{e^{E/T}-(-1)^{2s}}\,,$ (15)

where $\Gamma_{Eslm}$ is precisely the absorption coefficient encountered by
Teukolsky & Press, here interpreted as a spontaneous emission coefficient. It
is often denoted as “greybody” (or “graybody”) factor (GF) as it encodes the
departure of the BH from a pure blackbody: not all radiation is absorbed but
some part is reflected (or, equivalently, spontaneously emitted).

The rate of emission is given for massless scalar fields from a SBH but
extension to higher spin $s$ fields for a Kerr–Newmann BH is straightforwardly
obtained by adding the corresponding angular velocity and electric potential
into the Boltzmann factor

$\dfrac{{\rm d}^{2}N}{{\rm d}t{\rm
d}E}=\dfrac{1}{2\pi}\dfrac{\Gamma_{Eslm}(\Omega,\Phi)}{e^{(E-m\Omega-q\Phi)/T(M,a^{*},Q)}-(-1)^{2s}}\,,$
(16)

where $\Omega=4\pi a/A$ and $\Phi=4\pi Qr_{+}/A$ [98]. The tendency of BHs to
emit aligned momentum and same-charge particles is apparent from the form of
the effective chemical potential in Eq. (16).

For the emission of a particle of non-zero rest mass $\mu$ to be kinematically
allowed, the energy of emission must be $E>\mu$. However, the rate of emission
is exponentially suppressed at high energies, hence:

> As the temperature rose, it would exceed the rest mass of particles such as
> the electron and the muon and the black hole would begin to emit them also.
> [96]

Therefore, if the number of kinematically allowed degrees of freedom (dofs) of
emission increases monotonically with energy, such as in the Hagedorn model,
then the latest emission would indeed be explosive.

The evolution of BHs by HR is considered as quasi-static, so that backreaction
can be ignored, and the successive particles are emitted independently.
Furthermore, to avoid problems with Planck scale naked singularities, Hawking
assumed that when reaching this scale the system shall “disappear altogether”
in a final flash of energy. Many interpretations were subsequently given as
for this spontaneous emission. Hawking proposed that it could be grasped
considering the spontaneous creation of pairs of particles and antiparticles
just outside the horizon with a Boltzmann distribution given by the
temperature of the BH, with one half escaping to infinity with positive energy
and the second half falling inside the BH with negative energy, therefore
decreasing its mass. Another analogy proposed is the spontaneous decay of the
BH because of particles tunneling out of the horizon [99] (see in particular
the discussion by Parikh & Wilczek [100]). For other HR derivation and
interpretations, we refer the interested reader to [80].

###  3.2 Theoretical aspects

The Hawking process was thus firmly established on statistical mechanics,
quantum mechanics and general relativity grounds. That new process and the
corollary finite lifetime of BHs has immediate consequences on PBH
constraints. PBHs with mass $M<M_{*}$ cannot participate in the DM density
today, and their evaporation may leave observational signatures, as discussed
in the first ever PBH review by Carr [51].666A review was also published on
the Soviet side [101] but it went almost unnoticed. In fact, the evaporation
of small PBHs may be the only reasonable way to detect them.

Very early arguments by Carter and Gibbons [102, 103] showed that even if
formed with non-zero charge and angular momentum, a BH would lose them more
rapidly than its mass and end up as a simple SBH for most of its lifetime
$\tau$. For rotating PBHs, even if the temperature $T\rightarrow 0$ as the
angular momentum approaches the extremal limit $a^{*}\rightarrow 1$, enhanced
emission rates (wrongly associated to superradiance) should compensate so that
the lifetime $\tau(a^{*}\lesssim 1)\sim\tau(a^{*}=0)$. Hence, if formed
initially rotating, a PBH should keep a sizeable angular momentum for most of
its lifetime; this in fact contradicts the above statement. The situation is
very different for charged BHs since a rapid discharge process called the
Schwinger effect [104] neutralizes the BH independently of its HR rate in a
timescale much shorter than the mass loss, due to spontaneous pair-production
of charged particles (see also [105]). In fact, it is estimated that only
supermassive BHs formed with appreciable charge would still be charged
today—they are protected from the Schwinger effect by the non-zero mass of the
electron. Hence, for quite a long time, HR by PBHs was considered only for
SBHs. This was further supported by the simple formation mechanism proposed
for PBHs; they originate from the collapse of neutral radiation dominated
material.

The GFs represented the complicated part of the Hawking process analysis.
Apart from some high- and low-energy limits, the calculation of GFs requires
numerical tools. Page was the first to compute them numerically in 1976-1977
[98, 106, 107], using the programs of Teukolsky & Press [85, 86, 87]. The
first paper deals with SBHs and defines the “Page coefficients” as the
fundamental functions for computing the BH evolution, integrating over the
mass (and spin) loss through HR. The second paper deals with KBHs and Page
remarks that depending on the set of accessible dofs, the spin evolution may
be drastically different. A coupling between the fundamental dofs spin and the
BH angular momentum is identified, which causes enhanced emission of higher
spin particles for rapidly spinning BHs. The third paper deals with charged
BHs and massive particles, with the conclusion that a non-zero rest mass $\mu$
acts more or less as a kinematic cut-off on HR at $E>\mu$.

Once the “primary” emission rates of fundamental particles are known, remains
the difficulty of their subsequent interactions and decays. This is deeply
related to the model of particle physics considered, and has dramatic
consequences on the PBH constraints, as noted in the second review by Carr
[108]. The theory of HR for coloured particles and gauge bosons [109] was
rapidly implemented. However, the relative independence of emitted particles
remained a subject of debate until MacGibbon & Webber [110, 111] (hereafter
MG&W) settled the fundamental theory of BH “secondary” HR, with the hypothesis
that BHs radiate mostly independent particles that subsequently behave like
the products of e+-e- collisions in accelerators.777This work is based on
earlier discussion by Carter et al. [112] and in particular by Oliensis & Hill
[113]. The first paper deals with the instantaneous emission of particles and
the hadronization of the coloured ones; use is made of the particle physics
code HERWIG [114] to obtain the secondary spectra. The second paper deals with
the total emission integrated over the BH lifetime, with the conclusion that
most of the energy radiated by BHs is in the form of QCD jets originating in
the hadronization and decay of quarks and gluons. The composition of the
primary and secondary PBH HR products gives therefore privileged access to
particle physics at any energy up to the Planck limit $E\sim 10^{19}\,$GeV.

The MG&W model was challenged by others, mainly in the 1990’s, with very
different PBH constraint predictions:

  * 1.

Cline & Hong [115] argue that the correct model of HR must lie between the two
extreme Hagedorn and Gell-Mann models for hadrons (or equivalently between the
“elementary” and “composite” particle models);

  * 2.

Heckler [116, 117] claims that taking into account number non-conserving
interactions of the $2\rightarrow 3$ kind, the density of an expanding plasma
around an exploding PBH can grow very fast (see also the computations by
Daghigh & Kapusta [119]);

  * 3.

Belyanin et al. [120] predict that a full magneto-hydrodynamic plasma should
develop around evaporating PBHs, based or earlier ideas from Rees & Blandford
[121, 122];

  * 4.

Nagatani [123] argues that symmetry restoration is possible in a sphere around
an evaporating BH when its temperature overcomes the associated
threshold.888Hawking had suspected just the contrary in [124].

All of these competing models rely on the formation of some kind of
“photosphere” or “fireball” around PBHs, due to non-negligible interactions
between emitted particles [125]. They predict very different HR signatures in
the astronomical observations, in particular concerning the PBH final bursts.
The party has been definitely settled by MacGibbon et al. [126] in 2008 in a
detailed study of all the relevant interactions in the SM framework. The MG&W
model is now commonly used by all HR papers, and we stick to it hereafter,
with comments on other models at relevant places.

###  3.3 Basic formulas

In this Section we present the basic formulas used to compute the HR rates in
the MG&W model. Recall that the fundamental equation for the emission rate of
a single spin $s$ dof with angular momentum numbers $(l,m)$, rest mass $\mu$,
electric charge $q$ and energy $E$ and for a BH metric described by the mass
$M$ and the set of secondary parameters $\\{x_{j}\\}$ ($\\{x_{j}\\}$ can
contain the BH spin, its electric charge, etc.) is

$\dfrac{{\rm d}^{2}N}{{\rm d}t{\rm
d}E}=\dfrac{1}{2\pi}\dfrac{\Gamma_{Esqlm}(M,\\{x_{j}\\})}{e^{E^{\prime}/T(\\{x_{j}\\})}-(-1)^{2s}}\Theta(E-\mu)\,.$
(17)

The GFs $\Gamma$ and the BH temperature $T$ depend in general of the BH
metric. In the high-energy limit (“geometrical optics” limit or GO), they are
the same for all particle spins and in the Schwarzschild case one obtains

$\Gamma_{\rm GO}(E)\underset{\lim
E\rightarrow+\infty}{=}\sum_{l,m}\Gamma_{Eslm}=27E^{2}M^{2}\,.$ (18)

The GO limit for general spherically symmetric BHs is given by [127]. The
temperature is obtained in general from $T=\kappa/2\pi$, which becomes in the
KBH case

$T_{\rm
K}=\dfrac{1}{2\pi}\left(\dfrac{r_{+}-M}{r_{+}^{2}+a^{2}}\right)\quad\underset{a^{*}\rightarrow
0}{\longrightarrow}\quad T_{\rm S}=\dfrac{1}{8\pi M}\,.$ (19)

The temperature of a SBH, compared to some typical energy scales, is given in
Fig. 2, while the temperature of a general spherically symmetric BH is
obtained in [128].

Figure 2: Temperature of a SBH, compared to the energy scales of the CMB
($T\sim 2.7\,$K), the QCD ($\Lambda_{\rm QCD}\sim 150\,$MeV) and QED
($\Lambda_{\rm QED}\sim 100\,$GeV) phase transitions and the putative GUT
scale ($E_{\rm GUT}\sim 10^{14}-10^{16}\,$GeV). The mass scales are the same
as in Fig. 1.

The energy of the particle may be corrected by effective chemical potentials
due to e.g. spin or charge coupling: $E^{\prime}=E-m\Omega-q\Phi$. The rate of
emission of a primary particle $i$ with internal multiplicity $g_{i}$ (color,
helicity) is

$\dfrac{{\rm d}^{2}N_{i}}{{\rm d}t{\rm d}E}=g_{i}\sum_{l,m}\dfrac{{\rm
d}^{2}N}{{\rm d}t{\rm d}E}\,.$ (20)

The emission of a secondary particle $j$ must be computed by convolving the
primary spectra with particle physics branching ratios
$\text{Br}_{i\rightarrow j}(E,E^{\prime})$

$\dfrac{{\rm d}^{2}N_{j}}{{\rm d}t{\rm
d}E}=\int_{0}^{+\infty}\sum_{i}\text{Br}_{i\rightarrow
j}(E,E^{\prime})\,\dfrac{{\rm d}^{2}N_{i}}{{\rm d}t{\rm d}E^{\prime}}\,{\rm
d}E^{\prime}\,.$ (21)

Hence, a distribution of BHs produces the following instantaneous flux of
secondary particles

$\dfrac{{\rm d}^{2}n_{j}}{{\rm d}t{\rm d}E}=\int{\rm d}\\{x_{j}\\}\int_{M_{\rm
min}}^{M_{\rm max}}{\rm d}M\,\dfrac{{\rm d}^{2}n}{{\rm d}M{\rm
d}\\{x_{j}\\}}\,\dfrac{{\rm d}^{2}N_{j}}{{\rm d}t{\rm d}E}\,.$ (22)

This flux may be integrated over some period of cosmological time
$(t_{1},t_{2})$, taking the universe redshift dilution $\mathfrak{r}\equiv
a(t_{2})/a(t)$ into account, through

$\dfrac{{\rm d}n_{j}}{{\rm
d}E}(t_{2})=\int_{t_{1}}^{t_{2}}\mathfrak{r}(t)\left.\dfrac{{\rm
d}^{2}n_{j}}{{\rm d}t{\rm
d}E^{\prime}}\right|_{E^{\prime}=E/\mathfrak{r}(t)}\,{\rm d}t\,,$ (23)

where the time evolution of the PBH parameters should also be taken into
account.

The evolution of KBHs is encoded in the “Page coefficients”

$\begin{pmatrix}-f(M,a^{*})/M^{2}\\\
-g(M,a^{*})a^{*}/M\end{pmatrix}\equiv\dfrac{{\rm d}}{{\rm
d}t}\begin{pmatrix}M\\\ J\end{pmatrix}=\int_{0}^{+\infty}\sum_{i}\dfrac{{\rm
d}^{2}N_{i}}{{\rm d}t{\rm d}E}\begin{pmatrix}E\\\ m\end{pmatrix}{\rm d}E\,,$
(24)

where one should remember that the emission rates depend non-trivially on the
particle angular momentum (an analog coefficient could be defined for e.g. the
BH electric charge). The following differential equations are obtained
straightforwardly

$\left\\{\begin{array}[]{l}\dfrac{{\rm d}M}{{\rm
d}t}=-\dfrac{f(M,a^{*})}{M^{2}}\,,\\\ \dfrac{{\rm d}a^{*}}{{\rm
d}t}=\dfrac{a^{*}(2f(M,a^{*})-g(M,a^{*}))}{M^{3}}.\end{array}\right.$ (25)

and the lifetime of a BH is computed as

$\tau(M,a^{*})=\int_{0}^{M}\frac{M^{2}}{f(M,a^{*})}{\rm d}M\,.$ (26)

Assuming constant $f(M,a^{*})$ raises $\tau\sim M^{3}/(3f(M,a^{*}))$. With the
SM emission only, one obtains $\tau(M_{*})=t_{0}$ the age of the universe for
the critical mass $M_{*}\sim 5\times 10^{14}\,$g.

###  3.4 Numerical recipes

There are several methods to compute the quantities given by Eqs. (17)–(23).
This can be done analytically under some assumptions, e.g. in the geometrical
optics (GO) high-energy limit for HR rates [129]. One can also compute the
numerical rates and fit them with heuristic analytical formulas, see e.g.
[13]. As the GFs are not analytically calculable for general BH metrics at all
energies, numerical resolution is required. In general, the Teukolsky equation
is complicated to solve numerically as such. It can be reduced to a
Schödinger-like wave equation with short-ranged potentials, as shown in the
general spherically symmetric case by [128] or in the Kerr case by [130, 131,
132, 133]. This equation may be solved numerically in the WKB approximation
[100], through a path-integral formalism [134] or by brute force computing.

Several numerical public codes are available to compute the HR rates.
Historically, the first programs were released in the context of HR by small
higher-dimensional BHs evaporating in the LHC detectors, and took the form of
event generators. These include Charybdis [135, 136], Catfish [137], BlackMax
[138, 139, 140] and QBH [141]. The present author (together with A. Arbey)
released the C code BlackHawk [10, 142] which is more dedicated to PBH studies
with the possibility to simulate general mass and spin distributions and to
study the cosmological evolution of the HR emission. More recently, the PBH
evaporation in the context of DM production by PBH evaporation and higher-
dimensional PBHs were implemented in respectively the Python code ULYSSES [13,
14] and CosmoLED [11]. Comparison of the PBH codes, when possible, shows that
they are in reasonable quantitative agreement with each other.

The secondary spectra are computed by particle physics codes, which differ
following the energy scale considered. For example, the public code BlackHawk,
available at https://blackhawk.hepforge.org/, relies on the Python package
Hazma [143] at low energies ($E\lesssim$ GeV), PYTHIA [144] or HERWIG [145] at
LHC energies ($E\sim$ GeV$-$TeV) and on HDMSpectra [146] at high energies (TeV
upward). An example of photon spectra is given in Fig. 3 for each of these
cases. On that figure, the bumpy feature in the secondary spectra at $\sim
100\,$MeV is due to the primary pion decays. The codes are not reliable for
$E_{\rm sec}/E_{\rm prim}\ll 10^{-6}$, hence the HDMSpectra line stops at
$\sim 1\,$GeV. As predicted by MG&W, most of the emission is realized at
energies $E\ll T$ due to the jet hadronization and decay. In the context of DM
studies, Cirelli & collaborators constructed tabulated branching ratios for DM
decay that can be directly transported into PBH secondary spectrum
calculations [147, 148]. The material is available at
http://www.marcocirelli.net/PPPC4DMID.html and is generally denoted as
PPPC4DMID.

Figure 3: Photon primary (dashed) and secondary (solid) spectra from a SBH
with temperature $T=100\,$MeV (left, Hazma), $100\,$GeV (center, PYTHIA) and
$100\,$TeV (right, HDMSpectra) as computed with BlackHawk. The vertical lines
represent the BH temperature.

##  4 Non-evaporated primordial black holes

In this Section, we provide a very brief account of the constraints on the PBH
density which are not based on HR. This is only a brief summary as complete
recent accounts can be found e.g. in [3, 4]. Theses reviews focus on PBHs as a
DM candidate. Indeed, PBHs are non-baryonic and non-relativistic, while not
being a new object of particle physics like WIMPs or axions. They appear
therefore as a convenient DM candidate and the constraints on their abundance
today $\Omega_{\rm PBH}\equiv\rho_{\rm PBH}/\rho_{\rm c}$ (where $\rho_{\rm
c}$ is the critical density) are often set as the fraction of DM $f_{\rm
PBH}\equiv\Omega_{\rm PBH}/\Omega_{\rm DM}$ they can represent, with the upper
limit $f_{\rm PBH}<1$ not to overclose the universe. A more universal
parameter is the fraction of the universe density into PBHs at the time of
their formation $\beta\equiv\Omega_{\rm PBH}(t_{\rm f})/\Omega_{\rm
tot}(t_{\rm f})$, with the particular combination $\beta^{\prime}$ that could
encode non-standard cosmologies [9]; we stick to $\beta^{\prime}$ hereafter to
express the PBH constraints. PBHs with mass over the evaporation limit
$M>M_{*}$ are constrained mostly by their gravitational effects, such as
gravitational lensing, dynamics in structures, accretion and gravitational
waves. The constraints that rely on PBH evaporation are presented in detail in
the next Sections.

###  4.1 Lensing constraints

PBHs of any mass behave as gravitational lenses that deviate light rays due to
their strong gravitational field (for a review of this effect and its
implications in cosmology, see [149]). PBHs should have a sensible proper
motion with respect to the sources they lense the light of. Their radius is so
small that they only behave as very transient lenses, that is called
millilensing, microlensing, femtolensing or picolensing depending on the
duration of the effect.

The visible effect on Earth is that the light from remote sources is suddenly
brightened by the passing of a PBH due to the focusing effect [149]. This
applies to the light of quasars, supernovae, stars, and $\gamma$-ray bursts
(GRBs). The most famous studies are that of MACHO [150], EROS [151],999The
names of these two collaborations are obviously irrespective. OGLE
[152]101010OGLE might have observed microlensing events compatible with
planetary mass PBH origin [153]. and Subaru-HSC [154], which exclude $f_{\rm
PBH}=1$ in the range $\sim 10^{23}-10^{35}\,$g. The original Subaru-HSC
constraint was initially much more stringent and excluded PBHs as DM down to
$\sim 10^{21}\,$g [155] but it was shown recently that this was due to a
simplistic treatment of the source spatial extension [154].

There is a disputed constraint on GRB femtolensing [156, 157] in the
$10^{17}-10^{19}\,$g mass range which would close a good deal of the HR
window, and an ongoing search for fast radio burst lensing [158]. Microlensing
constraints are the ones with the most expected improvement in the Snowmass
PBH white paper [5].

###  4.2 Dynamical constraints

PBHs are massive objects, and thus they exert gravitational attraction on the
nearby material. This results: $i)$ in gravitational friction that make them
settle down in compact objects that they can subsequently absorb or destruct,
$ii)$ in gravitational disturbance of bound systems like stellar binaries, and
$iii)$ in gravitational heating of the interstellar medium (ISM) gas. The most
severe constraints in the $10^{36}-10^{47}\,$g come from gas heating and
galaxies, globular clusters and stellar binaries stability [159, 160] (see
[38] for a review). “Stupendously” large PBHs would also accrete surrounding
material with too much efficiency if formed in the early universe; depleting
structure formation [38]. There could be constraints on the disruption of
neutron stars and white dwarfs (the most compact non-BH objects) in the range
$10^{18}-10^{23}\,$g [161, 162] (effectively closing all the HR window), but
these were recently debunked [163, 164, 165] and await more rigorous
theoretical treatment.

###  4.3 Accretion constraints

PBHs may accrete gas and dust material into disks like any other celestial
object. These disks would radiate in the X-ray band and heat the surrounding
ISM or the CMB [166, 167], or the PBHs could be so massive that their
gravitational potential distorts the CMB [168] ($\mu$-distortions). These
constraints exclude all DM into PBHs in the super-Solar mass range $M\gtrsim
M_{\odot}$. Accretion indirectly takes a part in all other constraints as it
modifies the PBH mass spectrum.

###  4.4 Gravitational wave constraints

As GWs propagate in space with no measurable alteration, they are the most
precise probes of gravitational events back to the early universe. The
stochastic GW background encompasses the GWs generated by PBH formation, by
the statistical distribution of PBHs (scalar-induced GWs), by their sudden
evaporation and modification of the equation of state, and by their continuous
mergers since the early era. This may be the set of constraints totaling the
greatest number of publications since 2015 and the first observation of the
merger of BHs by the LIGO instrument [169, 170]. It was immediately suggested
that the BH components of the binary merger were PBHs [171, 172], and the
background of GWs from PBH mergers was predicted in [173]. The associated
constraints are model-dependent but exclude PBHs in the $10^{30}-10^{34}\,$g
mass range [174, 175] (observations by the LIGO/VIRGO instruments and the
NANOGrav experiment). For a review about the GW perspectives of PBH studies,
please refer to [6].111111For a more general review about GWs and the early
universe, see [176]. GWs from mergers are also one of the only positive
evidence for PBHs, as even if they represent just a percent of the DM density,
their mergers are still in the range of detectability at LIGO/VIRGO and future
detectors [177].

Peculiar GW events such as very heavy $M\gtrsim 100\,M_{\odot}$ or light
$M\lesssim M_{\rm TOV}$ BH binary mergers (e.g. [178, 179]) would be smoking-
guns towards PBH origin (see however [180]). A stochastic GW background should
also be generated by PBH formation [181, 182, 183], and most interestingly it
depends on the spin of the PBHs, giving access to a double-check of the PBH
origin of BH mergers [184]. It could be the only way of constraining the PBH
abundance in the remote $M\lesssim 10^{9}\,$g mass range for which PBHs
evaporate before BBN. The forthcoming LISA instrument is expected to shed some
light on these yet unmeasured features [185].

##  5 Baryogenesis

The universe is not symmetric in (anti)baryon content. The baryon-to-photon
ratio $\eta\sim 6\times 10^{-10}$ [186], as measured at BBN or CMB epochs
shows that if initially produced in the same quantities, most of the baryons
and antibaryons should have annihilated into photons. However, a slight
overabundance of baryons of the order of $\eta$ must have been produced, so
that after annihilation some baryon content was still present. This baryon-
antibaryon asymmetry resulted in baryogenesis, the genesis of the baryon
content of the universe. For baryogenesis to take place, the three conditions
of Sakharov [187] must be fulfilled:

  1. Cond. 1: C and CP symmetry violation;

  2. Cond. 2: baryonic number violation;

  3. Cond. 3: thermodynamic equilibrium violation.

In the standard cosmological model with only SM fields, these conditions are
not met with sufficient amplitude to explain baryogenesis. Thus, extensions to
the standard model must assumed. One possibility resides in the capability of
BHs to fulfill conditions 1 and 3 immediately. Condition 1 is satisfied
because as a PBH explosion manifestly violates the T symmetry, it should
violate C/CP as well in order to conserve the CPT symmetry resulting from the
global Lorentz invariance [108]. CP symmetry violation is also present in the
SM as seen e.g. in neutral kaon decay [188]. Furthermore, when BHs evaporate,
the temperature of the quasi-thermal flux they produce is uncorrelated to the
background radiation temperature. Therefore, they are sources of out-of-
equilibrium processes, which fulfills condition 3.

Condition 2 is trickier: baryonic number violation is present in the SM but
not with a magnitude sufficient to explain the value of $\eta$. Furthermore,
any asymmetry generated in the early universe must survive until BBN and CMB
epochs, thus it must not be diluted away by other mechanisms, such as wash-out
or entropy injection [189, 190]. In general, it is assumed that baryogenesis
should occur after the background temperature has cooled to $T_{\rm EW}\sim
100\,$GeV, for the electroweak (EW) interactions to freeze out. In the
standard RD model, this translates to a time $t_{\rm EW}\sim 10^{-12}\,$s
which is the lifetime of a $M\sim 10^{5}\,$g PBH. PBH participation in
baryogenesis may proceed by five main scenarios:

  * 1.

the “Hawking model” where baryogenesis is the direct result of an asymmetry in
the GFs (e.g. an effective chemical potential);

  * 2.

the “Weinberg model” where heavy (GUT) bosons emitted by HR of PBHs decay
asymmetrically into (anti)baryons;

  * 3.

the “Dolgov model” where a mass asymmetry between (anti)baryons cause
asymmetric re-absorption probability by PBHs;

  * 4.

the “Nagatani model” where an EW symmetry-restored region exists around high-
temperature PBHs and baryon asymmetry develops by sphaleron processes;

  * 5.

the “Fujita model” where the heavy decaying particles are right-handed
neutrinos (RHNs) instead of GUT bosons.

These five models fulfill condition 2, as showed below. We have listed them in
chronological order of appearance.

The natural mass range for baryogenesis is $M\ll 10^{9}\,$g so that PBHs
evaporate before BBN and their abundance is not constrained. Indeed, the
initial abundance of PBHs would have to be very large in order to produce the
correct asymmetry [191], PBHs may even have dominated the energy density of
the universe before their evanescence [196], an era denoted as BHD (for black
hole domination). However, a very extended PBH distribution like the original
scale-invariant one poses immediate difficulty as it would conflict with the
BBN constraints at $M>10^{9}\,$g [192].

In general, the asymmetry proceeds from a baryon non-conservation process
either in the PBH evaporation directly or in the PBH evaporation product
decay. Thus, the generated asymmetry is related to the PBH abundance by a
proportionality factor $\epsilon_{{\rm B}\overline{\rm B}}$, which depends on
the PBH mass and on the asymmetry mechanism

$\eta\propto\left\\{\begin{array}[]{ll}\epsilon_{{\rm B}\overline{\rm
B}}(M)\,n_{\rm PBH}\,,&\text{RD}\\\ \epsilon_{{\rm B}\overline{\rm
B}}(M)\,,&\text{BHD}.\end{array}\right.$ (27)

Reviews about PBH baryogenesis can be found e.g. in [193, 194, 195].

###  5.1 Hawking model

In the seminal paper [96], Hawking expressed the idea that PBHs may be formed
from symmetric (anti)baryonic content, and evaporate preferentially into
baryons rather that antibaryons, due to e.g. a chemical potential in the GFs
leading to asymmetric branching ratios. On the other hand, PBHs could be
formed by collapsing baryonic material, but as they have no memory of their
progenitor composition due to the no-hair theorem, they would carry no
baryonic charge. Hence, PBHs would radiate (anti)baryons in equal numbers and
contribute only in the radiation density, making the baryon-to-photon ratio
effectively diminish. These two mechanisms fulfill condition 2.

Barrow [196] shows first that the correct PBH baryogenesis model depends on
the capability of PBHs to dominate the energy density of the universe. If they
do, then baryogenesis can result from the “Hawking model”. In this case, the
asymmetry does not depend on the initial PBH abundance because they are at the
origin of a second reheating (the first reheating is the decay of the putative
inflaton field) and of all the radiation content of the universe. If they do
not, then it should proceed through e.g. the “Weinberg model”. There is a
region of overlap that depends greatly on the assumptions on $\epsilon_{{\rm
B}\overline{\rm B}}$.

Hook [197] proposes a realization of the “Hawking model” of baryogenesis via a
chemical potential in the HR rate that is not the same for (anti)particles.
This chemical potential is assumed to be present everywhere in the universe in
[197], but it may also be located in the BH vicinity due to modified EW
interactions in curved spacetime [198, 199] and it could depend on the BH mass
[200].

The chemical potential could also be generated by baryon/lepton number
interactions with the space-time curvature [201, 202]. Ref. [201] is the first
to consider baryogenesis with an extended mass function originating in the
critical collapse mechanism.

###  5.2 Weinberg model

The first complete scenario of baryogenesis based on PBH evaporation in a GUT
may also be the first appearance of the idea that PBHs can emit particle dofs
that have not yet been discovered in astronomical and laboratory experiments.
It was presented by Turner & Schramm in [203] and assumes that PBHs can
evaporate into a GUT-scale ($\sim 10^{16}\,$GeV) boson $X$ (that could be the
Higgs boson), whose decay would strongly violate CP. The subsequent asymmetric
X/$\overline{\text{X}}$ decays into B/$\overline{\text{B}}$ would generate a
net baryon density $n_{\rm B}$ (fulfilling condition 2) which would then be
diluted by the emission of the rest of the PBH mass under the form of
radiation (entropy) down to the observed value of $\eta$. This is based on the
idea by Weinberg [204, 205] but provides a natural origin for the GUT bosons
in the evaporation of PBHs. The first numerical calculation of the relevant
parameter space in the PBH mass/GUT boson characteristics was done in [206].

The GUT “Weinberg model” was greatly improved by Barrow & Ross [196, 207,
194]. They remark that if the PBH spatial distribution is not homogeneous,
then inhomogeneous baryon-to-photon ratio can develop. This could be at the
origin of structure formation. They also discuss the monopole problem: PBHs
could produce dramatic magnetic monopole abundance by HR, which is strictly
constrained. This would result in a strong reduction of the available
parameter space for PBH baryogenesis.

10 years later, Barrow et al. [208, 209] review the PBH GUT baryogenesis
mechanism. They present a version of it where the heavy GUT boson is the Higgs
boson. Baumann et al. [210] spotted a typo in the numerical calculations of
Barrow et al. and shrink accordingly the available parameter space in [210].
The inclusion of this scenario in the full MSSM was done in [211].

###  5.3 Dolgov model

This model is the reverse of the “Hawking model” and it does not really rely
on HR mechanics, even if absorption and emission coefficients are not
independent; it is given here for completeness. Toussaint et al. [191]
proposed that the absorption of surrounding particles by PBHs could be
asymmetric: if B/$\overline{\text{B}}$ have different interactions with the
gravitational field (e.g. different rest masses), then their absorption
coefficients are different and an initially symmetric medium would become
asymmetric because of PBH absorption.

Dolgov’s model [212] is based on the “Weinberg model”. Suppose that the PBH
radiated heavy bosons decay into mass-asymmetric mesons. These are
subsequently recaptured by PBHs with a different rate due to different
coupling with the PBH metric. The baryonic asymmetry is computed analytically
thanks to a resolution of the Teukolsky equation with an approximate potential
for the massive GUT fields in [213].

Ambrosone et al. [214] note that the “Dolgov model” can in fact be
complementary to the other baryogenesis models because the very high PBH
fraction necessitated by those may result in effective re-absorption of the
asymmetry produced. If this re-absorption is asymmetric, the baryon asymmetry
is preserved.

###  5.4 Nagatani model

A new model of baryogenesis by PBHs is introduced by Nagatani in 1999 [123].
They argue that the if PBHs dominated the early universe, their energetic
radiation might heat up the PBH environment so much that the EW (Higgs)
symmetry might be restored in a spherical region around PBHs. Then, a domain
wall should form between the restored symmetry patch and the outer broken
symmetry universe. Baryogenesis would occur inside that domain wall (or in the
restored region) via sphaleron processes [215] which fulfill condition 2. With
very fine tuning, this model could also explain reheating and DM [216]. The
“Nagatani model” seems deprecated since PBHs may not form symmetry-restoring
fireballs [124, 126].

###  5.5 Fujita model

Baryogenesis through PBH evaporation into asymmetric RHNs is considered by
Fujita et al. [217] as a mean to explain the coincidence between the (dark)
matter cosmological densities. This baryogenesis model is based on an old
tentative to avoid recourse to ad hoc GUT bosons [218]. The parameter space
for baryogenesis depends on the “light” ($T>\mu_{\rm\nu}$) or “heavy”
($T<\mu_{\rm\nu}$) mass of RHNs and is incompatible with DM generation except
if entropy non-conservation is allowed (even so, very fine-tuning is needed
[219]).

Recent numerical calculations done with the public code ULYSSES [13, 14] have
shown that interplay between PBH generated RHNs and thermal ones from Higgs
decay complicate the scenario [220]. The parameter space depends further on
the neutrino model [221, 222]. Dilution of the asymmetry by entropy injection
may be avoided with very light PBHs $M\gtrsim 10\,$g [223, 224] but in this
case the wash-out problem is still present.

###  5.6 Other ideas and prospects

A small amount of asymmetry could also originate in the randomness of the HR
process, with a fluctuation of the baryon number of the emitted flux of the
same form as that of the emitted charge [108].

Lindley [225, 226] argued that PBH baryogenesis models are un-natural since
they would require “entirely unjustified” fine-tuning in the PBH abundance and
mass in order to obtain the correct asymmetry. They argue rather that the PBHs
could be used as a second-hand source of asymmetry, while the first-hand
source should remain e.g. the thermal GUT mechanism.

Majumdar et al. [227, 228] were the first to study the impact of PBH early
accretion on the baryogenesis mechanisms. The overall conclusion is that the
initial increase in the PBH mass delays their evanescence, which would help
them decay after the EW freeze-out and thus avoid the wash-out problem.
However, this effect can be sizeable only for $M>10^{3}\,$g PBHs [229].

If PBHs generate the baryon asymmetry of the universe, this constitutes a
major positive evidence for their existence. The distinctive signature of such
models is that they predict a stochastic GW background of multiple origins:
PBH evaporation into gravitons, or PBH formation/Poisson distribution (scalar-
induced GWs) [181].

Overall, the PBH baryogenesis scenario is still perfectly viable, with the
downside that fine-tuning is required in every model. As all these processes
rely on the existence of BSM physics or non-standard BH behavior, they provide
no direct constraint on the PBH fraction $\beta^{\prime}$ in the early
universe. Constraints could come from the GW background (see Section 8.3),
from the density of PBH remnants (see Section 11.3.3) or from the concomitant
production of DM (see Section 11.4); these could allow observational
discrimination of the different models. Note that the full GFs for HR of heavy
particles (Weinberg or Fujita models) are never used, with the remarkable
exception of Ref. [220]. Other PBH baryogenesis scenarios would benefit from
precision calculations of HR.

##  6 Big bang nucleosynthesis

BBN corresponds to the epoch when the light elements where built up by fusion
and fission of nucleons into nuclei (by order of importance H, 3He, 4He, 7Li
and 6Li). It began when the universe was $t\sim 2\,$min old, and lasted
effectively until it was $t\sim 20\,$min old. The complete BBN mechanisms
extend out of this time window with the freezing of the neutron-to-proton
ratio at $t\sim 1\,$s when p$\,\longleftrightarrow\,$n inter-conversion by
neutrino capture became ineffective and the dissociation or decay of some
light elements (e.g. 3H, 7Be) into stable elements at $t\gg 20\,$min. These
light elements enter in the composition of the first generation of stars where
they are further processed into heavier elements up to Fe.

The standard BBN scenario works remarkably well, and is the more robust probe
of the very early universe so far. The light element abundances were observed
in distant (old) gas clouds, and correspond perfectly to the prediction of
numerical codes such as e.g. the public one AlterBBN [230, 231] that compute
the nuclear reaction chains taking the universe expansion into account. The
only free parameter of the standard scenario is the baryon-to-photon ratio,
and if set to the CMB value of $\eta\sim 6\times 10^{-10}$, all the relevant
abundances are obtained. This success remains however mitigated by the 7Li
problem, namely the discrepancy of a factor of a few between the predicted 7Li
abundance and its value as deduced from the Spite plateau in metal-poor
stellar abundances. Recent element abundances and discussions on BBN can be
found in the review [186]. PBHs modify the BBN in several ways, described in
the following:

  * 1.

modification of the neutron-to-proton ratio prior to BBN;

  * 2.

dissociation of elements by injection of hadrons;

  * 3.

dissociation of elements by injection of photons;

  * 4.

modification of the expansion history by injection of entropy.

As the BBN yields are finely tuned by the different reactions at stake, they
tolerate only minimal exotic interference. This results in very stringent
limits on PBHs with mass $M\sim 10^{9}-10^{13}\,$g.

The BBN constraints from PBHs long remained subject to large uncertainties
because of the precision of the nuclear reaction rates, the need for full
numerical Boltzmann equation resolution as the PBH products are non-thermal,
and the uncertainties of the high-energy PBH spectrum. It is only with the
MG&W model in the 1990’s that studies have correctly taken into account the
hadronic injection in the form of QCD jets. A summary of the modern
constraints is provided in Fig. 4.

Figure 4: BBN constraints on PBHs. The coloured solid lines are that of Ref.
[9], compared to the older limits of Ref. [232] (coloured dashed lines). For
completeness, we also plot in grey the photodissociation limits obtained by
Refs. [233] (dotted) and [234] (dashed), and the limit of [235] (dot-dashed)
obtained by conversion of DM constraints. The shaded area is robustly
excluded.

###  6.1 Nucleon inter-conversion

PBHs of $M\gtrsim 10^{9}\,$g can produce energetic neutrinos and hadrons after
the neutron-to-proton ratio has frozen out. Hence, they modify the initial
conditions of BBN. This has a sizeable effect on the final yields as all the
supplementary neutrons that are produced will ultimately end up in helium,
with a corollary reduction of the deuterium abundance.

Vainer & Nasel’skii were the first to examine the effect of neutrino injection
in 1977 [236, 237] (followed by Zel’dovich et al. [238] for the nucleon
injection). They rely on the Wagoner code [239, 240, 241] to compute the light
element abundances. All these effects modify the abundance of D at the end of
BBN, which is tightly constrained by observations.

Kohri & Yokoyama recomputed the nucleon interconversion constraints in 1999
[242]. Their paper is inspired by the work performed for particle DM decay
during BBN by Dimopoulos and collaborators [243, 244, 245]. They used new
observational data on the D, 4He and 7Li abundances, and on the neutron
lifetime $\tau_{\rm n}$. They studied the alteration of (early) BBN by hadron
injection by $10^{8}-3\times 10^{10}\,$g PBHs at $t\sim 10^{-3}-10^{4}\,$s. As
the timescale of particle interaction in the BBN plasma is
$\mathcal{O}(10^{-8})\,$s, they considered all hadrons with that lifetime in
their cascade reactions: charged pions $\pi^{\pm}$, charged kaons K±, the long
neutral kaon KL and the (anti)nucleons. They estimated the rate of emission of
these particles from jets in the MG&W model, for the first time with unstable
particles. The other electromagnetic (EM) particles are assumed to decay or
thermalize with the plasma. They showed in particular that the constraints are
highly dependent of the baryon-to-photon ratio $\eta$, which is left free in
their study (as well as $M$ and $n_{\rm PBH}$), and on the choice of the D
data, which suffers from high uncertainties.

###  6.2 Hadro-dissociation

Hadrodissociation means dissociation of light elements by hadronic injection,
namely light mesons (pions, kaons) and nucleons (protons, neutrons) radiated
by PBHs.

Zel’dovich et al. [238] analytically studied the effect of PBHs on the
abundance of D. They find that energetic hadronic injection can destroy helium
and overproduce deuterium. Numerical calculations are performed with the
Wagoner code by Miyama & Sato [246]. Vainer et al. [247] further discussed
hadrodissociation of light nuclei by energetic neutrons from PBHs and conclude
that the effect is globally an increase in the production of D and a decrease
of 4He. The energetic protons, due to their non-zero electric charge, would
thermalize in the vicinity of the PBHs.

Sedel’nikov returned to the problem of standard BBN alterations by evaporating
PBHs into energetic hadrons in 1996 [248]. They took into account the fact
that the PBH p/$\overline{\rm p}$ yield comes uniquely from jet decays,
following the MG&W model. They obtained constraints on the abundance of PBHs
slightly more stringent that what was obtained before in the
$10^{9}-10^{13}\,$g mass range. They also discussed the effect of adding BSM
dofs on the evaporation of the PBHs: it mainly reduces their lifetime and
their branching ratio into p/$\overline{\rm p}$.

###  6.3 Photo-dissociation

Vainer & Nasel’skii [236] argued that the very high energy $\gamma$-rays
emitted at the end of PBH evaporation could break 4He into D by photo-
dissociation. Lindley [249] further showed that high energy photons can
trigger EM cascades and then produce more high energy photons, leading to the
late photodissociation of D after the BBN yields have frozen out.

The photo-dissociation constraints have been revisited by Acharya & Khatri
[233]. They considered only EM cascades from $10^{11}-10^{16}\,$g PBHs, and
followed the procedure of [250] for DM decay. Full GFs are taken into account,
secondary particles are computed with PYTHIA and results are presented for the
first time for a log-normal mass distribution of PBHs.

The photo-dissociation constraints were studied by Luo et al. [234] with a
full spectrum Boltzmann treatment (based on the DM study [251]). They computed
the constraints for critical collapse mass function, for the first time.
BlackHawk is used to compute the precise spectrum of secondary particles and
the evolution of PBHs. The Kawano code is used with the REACLIB database for
the BBN computations. Only a very narrow mass range around
$10^{12}-10^{13}\,$g is constrained.

###  6.4 Lithium problem

Jedamzik [252] studied the general impact of “something” decaying during BBN
(DM, PBHs), with a modified version of the more Kawano code. Under some
circumstances, namely fine-tuned neutron injection during BBN, the 7Li problem
is resolved. As can be expected, this is more complicated with full Boltzmann
injection of un-thermalized particles generated by PBH evaporation. In fact,
it has been proven that a monochromatic distribution of Schwarzschild PBHs
cannot solve the lithium problem because their temperature is unequivocally
related to their lifetime, forbidding fine-tuning [9].121212Unpublished work
performed by the author with J. Pradler reached the same conclusions. An
extended mass distribution may alleviate that conclusion by destroying 7Li
without overproducing D.

###  6.5 Modern studies

Carr et al. [232] were the first to compute the BBN constraints from PBHs in
detail in the MG&W model, taking into account all above phenomena. This work
is a direct continuation of Ref. [242]. The new constraints rely on a modified
version of the Kawano code v4.1 to include PBH injection and the relevant
interactions are borrowed from DM studies (e.g. [253] and references therein).
The full energy spectrum of the PBH products are taken into account and
convoluted with the nuclear chains with an iterative method.

The BBN constraints are revisited by Keith et al. [235], based on the latest
Planck measurements [186] and of direct conversion of constraints from the
more recent particle decay theory [254, 255]; however only pure blackbody HR
is used. They discuss the impact of a dark sector: the evaporation rate is
faster, and the branching ratio into SM particles is weaker. The expansion
history is modified by PBHs and their (dark) decay products, which is taken
into account using AlterBBN.

In their most recent review, Carr et al. [9] update their BBN constraints from
Ref. [232] and discuss the results of [233, 235] which obtain stronger limits.
They argue that the hadron and photon effects must be taken into account
together as they may cancel each other, as explained e.g. in [254, 255] on
which they base their calculation. Their results is thus slightly less
stringent while they use the same recent data from Planck. Their constraints
are given as solid colored lines in Fig. 4, as we consider they are the most
robust to date. Older constraints from [232] are shown for comparison,
together with that of Refs. [233, 234, 235].

Overall, BBN is a very well understood process that strongly constrains the
light PBH abundance. It does not seem that PBHs can have any positive impact
on that cosmological era, thus we expect that constraints will become more and
more stringent with additional light element data, in particular deuterium
observations.

##  7 Cosmic microwave background

The measure of the CMB by Penzias & Wilson in 1965 [256] was one of the chief
evidence in favor of the hot Big Bang model, together with the BBN yields. The
CMB has been observed by the satellites COBE [257], WMAP [258, 259] and more
recently Planck [260, 261], to reach an extraordinary sensitivity. Future
measures with the CMB-S4 experiments should increase further the understanding
on the CMB epoch [262, 263, 264].

The isotropic CMB light corresponds to the first light emitted by atoms when
nuclei and electrons combined and photons decoupled form the plasma. The CMB
had a pure blackbody distribution at $E=13.6\,$eV originally, redshifted today
to a $T\sim 2.7\,$K background. This light was emitted at $t\sim
380\,000\,$yr. The decoupling process encompasses all the $\Lambda$CDM
physical ingredients: the spectrum of the CMB is extremely sensitive to the
composition of the universe, the distribution of the energy density between
(dark) matter and radiation. The spatial fluctuations of the CMB measured by
Planck give stringent constraints on the amount of energy injection, as well
as the departure from pure Planckian spectrum measured by COBE/FIRAS. The CMB
epoch is followed by the Dark Ages, namely the period before the burning of
the first stars and galaxies. During this period, some reionization process
should have occurred, with disputed origin. Another signature is the $21\,$cm
background which contains information about an hyperfine transition of the
hydrogen atom. The reionization fraction of the universe and the $21\,$cm
“temperature” $T_{21}$ give access to the Dark Ages period. CMB physics has
benefited from numerous public computational tools, such as CLASS [265, 266,
267, 268], and many modules adapted for non-standard cosmology, including PBHs
(HyRec [269, 270], DarkHistory [271]).

The simplest idea linking PBHs and the CMB would be that PBHs of temperature
$T\sim 2.7\,$K distributed isotropically generate the observed background with
steady state radiation. Nartikar & Rana [272] showed that this is unrealistic.
However, evaporating PBHs modify the genuine CMB in several ways, described
chronologically in the following (see [273] for a review):

  * 1.

energy injection may distort the CMB spectrum and make it depart from pure
blackbody;

  * 2.

energetic EM particles may reionize the universe, leading to:

  * 3.

damping of the CMB anisotropies;

  * 4.

heating of the $T_{21}$ temperature (modified optical depth of Lyman-$\alpha$
signal).

All these processes rely on energy injection in the form of EM particles, that
is photons, charged leptons and charged mesons (see the seminal work by
Nasel’skii [274, 275]). The other particles are either unstable on the CMB
timescales or not produced by PBHs that survive until then. Hence, PBHs that
can modify the CMB must have initial mass $M\gtrsim 10^{11}\,$g. The fact that
CMB disturbance relies most on the electron injection makes these constraints
competitive only for $M<10^{18}\,$g.

###  7.1 CMB distortions

The CMB is characterized by a power spectrum very close to a pure blackbody of
initial temperature $13.6\,$eV corresponding to the ionization energy of the
hydrogen atom. After the redshift due to the expansion of the universe from
$t\sim 380\,000\,$yr to $t_{0}\sim 13.9\,$Gyr, the CMB is measured today as a
microwave background of temperature $T_{\rm CMB}\sim 2.7\,$K [276]. This light
is due to the recombination of electrons onto protons to form neutral hydrogen
atoms. Exotic EM injection, e.g. by evaporating PBHs, contributes to or
scatters on the CMB photons. Hence, the pure blackbody can be distorted in
several ways [273]:

  * 1.

$g$ distortions: a shift of the blackbody temperature away from $13.6\,$eV;

  * 2.

$\mu$ distortions: an effective chemical potential in the CMB spectrum;

  * 3.

$y$ distortions: a departure from equilibrium of the Compton scattering
process.

These are thermal distortions, but the CMB can be distorted by non-thermal
mechanisms, such as PBH injection of energy. In this case a dedicated
treatment is required. Ref. [277] hence distinguishes “yim” distortions from
“ntr” ones. All of these are severely constrained by CMB observations, the
most precise to date being COBE/FIRAS. As the CMB distortions would be caused
at the time of recombination, it concerns PBHs of mass $M\sim
10^{11}-10^{13}\,$g.

Zel’dovich et al. [189] were the first to estimate PBH constraints from CMB
distortions. PBHs evaporating after the temperature of the universe cools down
to $T\sim m_{\rm e}$, emit energetic radiation that cannot reach thermal
equilibrium with the plasma, resulting in a distortion of the CMB. Nasel’skii
& Shevelev [278] then considered the Boltzmann equations describing the
inverse Compton scattering of e± on the CMB and the induced distortion, but
they were limited by numerical capabilities at that time. Tashiro & Sugiyama
[56] studied the $\mu$ and $y$ distortions of the CMB by PBHs in the
$10^{11}-10^{13}\,$g range, using the WMAP data, considering only the primary
photon spectrum. The constraints were revisited by Acharya & Khatri [277],
based on the DM study [279] and the theory on EM cascades in [280] (see also
[281]). They used COBE/FIRAS data that result in constraints less stringent
than the BBN ones but predict major improvement with future PIXIE measures.
The PBH mass evolution during evaporation is taken into account and secondary
particles are also considered thanks to the PPPC4DMID package. The CMB
distortions have been corrected by Chluba et al. [282] in the case of very
large energy injection that depart from simple linear analysis. This results
in mild modifications of the PBH constraints from [277]. The constraints from
Refs. [273, 277] and the corrections from Ref. [282] are shown in Fig. 5.

Figure 5: CMB distortion constraints on PBHs. The solid lines represent the
more robust constraints by Acharya & Khatri [277] (AK20, solid green is the
ntr limit, while dashed green is the yim limit) and Chluba et al. [282] (C20,
solid red is the high-injection corrected version of dashed red). The older
constraints from Tashiro & Sugiyama [56] (TS08) is given in dashed blue for
comparison, while the prospective limits from PIXIE are shown as discontinuous
grey lines (the PIXIE design is not definitive, which explains the differences
between Refs. [277, 282], resp. AK20 and C20).

###  7.2 Reionization

In a series of 3 papers, Gibilisco [283, 284, 285] computed the reionization
of the universe due to EM energy injection by PBHs. The FIRAS/COBE data impose
strong constraints on the fraction of the universe reionized after the CMB. In
the first paper, they used the old Page spectra for the injected energy and
obtained very strong limitations $\Omega_{\rm PBH}\lesssim 10^{-8}-10^{-12}$
in the $M\sim 10^{14}\,$g range. They considered PBH remnants in the second
paper. In the third one, they finally used the MG&W model to compute the total
amount of energy injected into the CMB: both direct EM and QCD decays are
taken into account. Interestingly, the constraints are presented with a
dependency on the Hubble parameter. Indeed, the mass of PBHs evaporating at
the reionization epoch depends on the age of the universe. Reionization by
PBHs was reviewed by Belotsky et al. [286] (see also [287, 288]). They
computed the reionization due to $10^{16}-10^{17}\,$g PBHs, saturating the PBH
bounds from $\gamma$-rays. The distinctive signature would be that such a
reionization should be homogeneous, whereas that from the first stars would be
inhomogeneous. They proposed a way of alleviating the constraints on PBHs by
using their idea of PBH clusters [289]. Clustering would allow to “disregard”
the Galactic $\gamma$-ray limit. The monochromatic and power-law mass
functions considered result in a negligible reionization, while the two-peaked
function of [290] provides enhanced effect. Reionization by axions produced by
PBHs is studied by Schiavone et al. [291]. They obtain the first constraints
on axions by PBHs with the RecFast code by adapting the study of [292].

Reionization has two major observational consequences. The first one is to
dampen the CMB anisotropies, the second one is to modify the Lyman-$\alpha$
optical depth, and thus the $21\,$cm signal. Both are used below to obtain
constraints on PBHs.

###  7.3 CMB anisotropies

CMB anisotropy damping proceeds from EM energy injection at the late CMB
epoch, it then concerns PBHs with $M>10^{13}\,$g. In fact, it is a direct
corollary of the reionization mechanism as the free electrons interact with
the CMB photons and dampen their power spectrum [293]. Any EM energy injection
that alters the ionization fraction sufficiently late hence modifies the CMB
anisotropy spectrum.

Dorosheva & Nasel’skii [294, 295], and Nasel’skii & Polnarëv [296] assessed
how the CMB anisotropies could help constrain the DM candidates. They
presented a complete computation of the $\Delta(T)/T$ anisotropy damping of
the CMB, related to non-standard reionization by either decaying massive
particles or PBHs, both envisaged as decaying “hidden mass” candidates. Carr
et al. [232] were the first to apply this constraint to PBHs, for which they
obtained a limit in the $M\sim 10^{13}-10^{14}\,$g mass range adapted from
that of decaying particles [297], based on WMAP data.

This first detailed numerical study was performed by Poulin et al. [293]. They
used the time-dependent evaporation products of PBHs and the code CLASS to
study the effect of exotic energy injection. The bounds are competitive with
the EGXB ones in the $10^{14}-10^{17}\,$g mass range with the Planck 2015 data
(see also [298]). However, the HR spectrum was simplified to the extreme by
considering the GO limit at $E>3T$ and no emission below. The constraints are
alleviated if the $\Lambda$CDM parameters are left free during the
constraining procedure [299, 63]. The CMB constraint in the
$10^{15}-10^{17}\,$g mass range was recomputed by Poulter et al. [63] using
Planck 2018 data.131313See also [273] which relied on the “on-the-spot”
approximation, meaning immediate thermalization and resulting in stronger
constraints; but used complete secondary estimation with PYTHIA. Explicit
extended mass functions (uniform and log-normal) are used for the first time
but only primary particles were considered.141414Interestingly, they compared
the results obtained with these genuine mass functions to the conversion
method of Ref. [59] and a good agreement for log-normal mass function is
found, whereas the uniform mass function shows $\mathcal{O}(1)$ discrepancy.
Stöcker et al. [300] used the same simple primary HR parametrization as [293],
however, they computed the secondary particles with PYTHIA, with particular
attention given to the QCD phase transition at $E\sim\Lambda_{\rm QCD}$. The
CMB constraints were revisited by Acharya & Khatri [233], based on the
particle DM study [250]. The constraints of Refs. [293, 300, 233] are shown on
Fig. 6.

Prospective CMB anisotropy constraints on PBHs were computed by Cang et al.
[301] using e.g. the CMB-S4 designs (see Fig. 6). Several kinds of extended
mass functions were considered and full HR rates were extracted from
literature. HyRec was used to produce mock data for these instruments and
extended limits are obtained from the Carr method. Finally, the case of
spinning PBHs was treated in the subsequent paper [302] thanks to BlackHawk,
with extended mass functions as well (Carr method). They are the first to use
numerical GFs, thus we have chosen to represent their constraint in Fig. 6 as
the robust limit. CMB-S4 is predicted to improve the constraints a little in
the higher mass range.

Figure 6: CMB anisotropy constraints on PBHs. The solid lines represent the
results of Poulin et al. [293] (P17, blue), Stöcker et al. [300] (S18, green),
Acharya & Khatri [233] (AK20, red) and Cang et al. [302] (C22, purple). The
prospective CMB-S4 limit is from Ref. [301] (C21, grey dashed).

###  7.4 21 cm signal

PBH late EM energy injection provokes heating of the intergalactic medium
(IGM) during the Dark Ages. This in turn affects the optical depth of hydrogen
as the equilibrium of the hyperfine states HI and HII depends on the IGM
temperature. Hence, the redshifted $21\,$cm data contain the whole thermal
history of the hydrogen gas after recombination, during the Dark Ages. This is
a unique probe of an otherwise dark period. Other processes impact the
$21\,$cm temperature, such as particle DM [308], baryon cooling [309] and
X-ray heating by stars [310].

At the time of Gibilisco’s papers [283, 284, 285], no Lyman-$\alpha$ data were
available, only upper limits. The $21\,$cm signal as a constraint on PBHs was
evaluated by Mack & Wesley [303] with first comparison to prospective
instrument sensitivity. They noted that:

> […] future 21 cm observations can provide better constraints on PBHs than
> are currently available [due] to a coincidence between the Hawking
> temperature of PBHs that evaporate during the Dark Ages and a window of low
> optical thickness of the IGM to photon absorption. [303]

This explains why the $21\,$cm signal has received so much attention lately.
Ref. [303] modified RecFast to include the full PBH spectrum interaction with
the IGM temperature and the reionization process, and obtained constraints on
PBHs from prospective $21\,$cm surveys: LOFAR, SKA, EDGES, 21CMA, MWA.

The first EDGES data were finally released in 2018 [304], with a trough in the
$21\,$cm temperature deeper than expected from standard reionization history.
Clark et al. [305] were the first to derive the PBH constraints in the
$10^{15}-10^{17}\,$g mass range. They took only primary photons and electrons
into account and discussed prospective data from PRIZM, HERA, LEDA and SKA.
Yang [306] recomputed the $21\,$cm constraints on PBHs with the code CAMB,
based on work done for particle DM decay [307]. The primary and secondary EM
particles were taken into account and the lower mass range
$10^{13}-10^{15}\,$g was explored.

Spinning PBH constraints have been obtained by [302, 311] using BlackHawk,
with typically more stringent limits due to the enhanced HR rates.

Saha & Laha [312] noted that the SARAS 3 experiment seems to mitigate the
EDGES data. Therefore, they only considered putative $21\,$cm measurements and
conservatively produced prospective constraints, taking the PBH spin into
account and using BlackHawk to compute the EM energy injection and DarkHistory
to compute the $21\,$cm signal for the first time. Full EM cascades in the IGM
plasma and PBH evolution were finally considered by Cai et al. [302]. Indeed,
using the full redshift dependency of the $21\,$cm trough rather than just the
$21\,$cm temperature allows to place more stringent constraints [310].

The $21\,$cm constraints of Refs. [312, 302, 310] are shown in Fig. 7. It is
difficult to assess which limit is the more robust as different models and
external effects are considered in those papers. Furthermore, since the
$21\,$cm data have not yet been confirmed by other surveys, we leave these
limits as prospective but note that they would constitute the strongest to
date in the $M\sim 10^{14}-10^{18}\,$g mass range. The model of Cang et al.
[302] seems to embed all the relevant PBH ingredients.

Figure 7: $21\,$cm constraints on PBHs. All constraints are shown as dashed
lines as the EDGES results have not yet been confirmed by other instruments.
The limit of Saha & Laha [312] is shown as a blue dashed line (SL21), with the
EDGES error bars echoed as a blue shaded area. The limits of Cang et al. [302]
(C22) and Mittal et al. [310] (M22) are represented as resp. green and red
dashed lines.

##  8 Cosmological backgrounds

The cosmological backgrounds are conceptually the simplest limits on the PBH
abundance from HR; hence they were the first ones to be computed. Assume that
there exists a cosmological (or present) density of PBHs, then they evaporate
by generating a continuous flux of (B)SM particles. If those particles are
cosmologically stable (photons, neutrinos, gravitons, DM particles), they
travel towards Earth in straight lines. If not, they decay into stable
particles that contribute as secondary components to the stable backgrounds.
The background is then obtained by stacking and redshifting the instantaneous
emission with formula (23) with correct $t_{1}$ and $t_{2}$. The limiting
times are obtained simply as follows. $t_{1}$ is the minimum time, it should
start from the PBH formation, but interaction of the particles introduces an
optical depth factor $e^{-\tau(t,E)}$ that effectively cuts off the integral
at $t\sim 1\,$s for neutrinos and $t\sim t_{\rm CMB}\sim 380\,000\,$yr for
photons. The limiting time for DM particles depends on their
(self-)interactions while gravitons propagate freely in any era of the
universe. $t_{2}$ is the maximum time, it is simply the minimum between the
age of the universe $t_{0}$ and the lifetime of PBHs $\tau$. These backgrounds
are described in Sections 8.1, 8.2, 8.3 and 11.4.

The quasi-instantaneous contribution from particles currently produced by PBHs
in the Galaxy or in defined targets (M31, dSphs, etc.) may be added to
describe more accurately the flux along a particular line of sight (see
Section 8.4).

###  8.1 Photons

The photon constraint was the first ever constraint on PBHs linked to their
HR. Hereafter, we denote as “EGXB” the extragalactic $\gamma$/X-ray
background. Chapline [313] depicted the PBH evaporation contribution to the
photon background as:

> […] a continuation of the 3 K blackbody spectrum and [rising] slowly to peak
> at an energy determined by the smallest black hole mass now existing. [313]

The strongest PHL limit would come from PBHs evaporating today with $M\sim
M_{*}$ initially, as they dominate the integral (23). This corresponds to a BH
temperature $T\sim 100\,$MeV, an energy range that was poorly explored.

This background is the simplest as the photon detection techniques have been
mastered for a long time. Furthermore, photons are reasonably products of any
interacting (particle) DM candidate, unless they interact only
gravitationally. The characteristics of the diffuse photon background could
then discriminate between the different DM models [314, 315], and the
refinement of the “Grand Unified photon spectrum” [316] over the years [317]
has progressively tightened the constraints.

Following the first numerical calculation of the HR rate by Page [98, 106,
107], Page & Hawking [318] estimated precisely the constraint that the
$100\,$MeV measurement of the diffuse background by Fichtel et al. [319]
imposes on the cosmologically stacked abundance of PBHs. They obtained that
the stacked flux follows a power-law ${\rm d}F/{\rm d}E\sim E^{-3}$, which
when compared to $\gamma$-ray measurements raised a constraint $\Omega_{\rm
PBH}\lesssim 10^{-8}$ around $M_{*}$. This computation could be applied to
other fields, with somewhat different optical depths, but:

> It would be very difficult to detect the gravitons or neutrinos because they
> have such small interaction cross-section. [318]

(see however below). Concerning CRs, one further difficulty is that:

> The charged particles would be deflected by magnetic fields and so would not
> propagate freely to Earth. [318]

thus requiring sophisticated propagation models. This first constraint was so
tight that it immediately implied that PBHs could not be very abundant in the
universe, and the resulting limit on their (original) power-law spectrum
resulted in strong constraints on the amplitude of primordial fluctuations. It
was denoted as the “Page Hawking limit” (PHL). No observations of the diffuse
photon backgrounds were performed until the 1990’s. Hence, the PHL was not re-
assessed until then.

Halzen et al. [320] returned to the diffuse photon bound in their 1991 review.
This review is very interesting in that it is the first to use the newly
published MG&W model [110, 111]. The low-energy tail of the photon emission
spectrum is greatly enhanced by taking into account the decay of QCD jets, but
this also corresponds to regions of more abundant background. Hence, the DM
fraction into PBHs is not constrained above $M\sim 10^{17}\,$g. A more precise
constraint $\Omega_{\rm PBH}=7.6\times 10^{-9}$ was obtained using the same
old data from Ref. [319]. MacGibbon & Carr [322] obtained the same order of
magnitude with the lower energy [321] data. The interactions between the PBH
photon flux and the IGM are discussed, with the conclusion that only photons
emitted by PBHs with $M\gtrsim 10^{13}\,$g after the CMB can reach detectors.

Great improvement came from the COMPTEL [323] (1996) and EGRET [324] (1998)
measurements which observed the MeV$-$GeV energy range, first used by Kim et
al. [325] in the PBH context. Following the discovery of the “critical
collapse” mechanism for PBH formation, Kribs et al. [326] refined the PHL for
this extended mass function numerically, which can be compared to the
analytical extrapolation of Yokoyama [52]. The first estimation of the limits
for “bumpy” mass functions was given in [327]. Early accretion impact was
considered in [328].

The cosmological PBH flux slope ($E^{-3}$) is not the same as the diffuse
background one, thus PBHs cannot be the sole source of this background
[329].151515See however Refs. [125, 330] in which models challenging the MG&W
one are used to comply more to the data. Barrau et al. [331] were the first,
to our knowledge, to include a background component to the PHL. They argue
that the contribution of blazars and active galactic nuclei (AGNs) is
“guaranteed”, translating in a reduced parameter space for PBHs $\Omega_{\rm
PBH}\lesssim 3.3\times 10^{-9}$.

In their 2010 review, Carr et al. [232] used for the first time the Fermi-LAT
data [332] above the GeV energy scale to obtain the limit on PBHs taking the
secondary spectra into account (see also [9] with more recent Fermi-LAT data
[333]). The secondary spectra calculation relies on PYTHIA, but the transition
from primary only to primary and secondary photon spectrum is abrupt because
of the computation method. The strongest limit $\Omega_{\rm PBH}\sim 10^{-10}$
is obtained at $M_{*}$.

Arbey et al. were the first to compute the PHL limit with numerical GFs thanks
to BlackHawk [334]. This was in fact the first application of the code to PBH
constraints, with a complete study of the impact of extended mass
distributions and non-zero PBH spin (see also [335, 336]). The latest Fermi-
LAT data were used [333]. The enhanced HR rates of KBHs results in more
stringent constraints, and extended mass functions broadens the excluded
region towards higher masses (up to $10^{18}\,$g for broad functions of near-
extremal PBHs). This opened a new period of intense work on the PHL. A low
energy X-ray background was assumed in Ref. [337] to obtain more stringent
constraints. Ref. [338] in particular considered the latest AGN and star
forming galaxy models as well as PBH electron-positron contribution to the
photon flux. Low-energy e± annihilation was also studied in Refs. [339] and
multiple extended PBH distributions were constrained.

Despite repetition of the EGRET survey [340], the lower MeV range is, to this
date, still unexplored with a high resolution instrument (better than
COMPTEL). Several propositions have been made, including AdEPT, AMEGO,
ASTROGRAM, GECCO, GRAMS, MAST, PANGU and XGIS-THESEUS. On the one hand, the
possibility to reach very low energies with high accuracy and effective
detection area makes the prospective limits from XGIS-THESEUS among the most
stringent ones with the $21\,$cm signal [336]. On the other hand, the
necessity to account precisely for the MeV PBH physics requires dedicated sub-
GeV particle physics codes, such as Hazma. The latter was proposed as an
amelioration to the original low-energy extrapolation of the PYTHIA results
inside BlackHawk v1 [341], and was implemented recently in BlackHawk v2 [142].

However, low-energy hadronization suffers from large uncertainties, inner
Bremsstrahlung effects may dominate at keV energies [481], the design of
prospective instruments is not definitive, statistical treatments and
background models vary from one study to the other. All these aspects make it
very difficult to compare the PBH constraints from one study to the other.
This was clearly demonstrated by the author in Fig. 7 of [342]: letting aside
the instrumental design, the previous uncertainties result in $10^{5}$
relative error bars on the PHL. Ref. [342] advocates for clear assessment of
the assumptions in order to compare limits with unified method, using e.g. the
BlackHawk tool Isatis. Furthermore, even with high-resolution surveys down to
eV energies, the overwhelming background and the dilution of the PBH density
at $M\gtrsim 10^{20}\,$g cause the PHL to stop being a competitive bound on
the PBH abundance (see also [343]).

Finally, even if PBHs contribute to some extent to the diffuse background,
there is a degeneracy with unresolved sources and other unknown background
components. Multimessenger correlations with e.g. stochastic GWs from PBH
formation could remove the degeneracy and provide positive PBH detection [344,
345]. We show the PHLs from Refs. [9, 338, 342] in Fig. 8, taking Ref. [342]
as the most conservative one.

Figure 8: EGXB constraints on PBHs. The solid lines are the limits set by
Refs. [9] (C21, blue) and [342] (A22, green). The (dot-)dashed lines are the
ones obtained in Ref. [338] (C22) with background models of [346, 347] (model
1) and [346, 348] (model 2). The grey shaded area is robustly excluded by Ref.
[342] while the green hatched area shows the $\sim 10^{5}$ uncertainty region
of the PHL.

###  8.2 Neutrinos

As stated above, neutrinos were at first too elusive to be detected directly
inside detectors [349]. But the neutrino background is particularly
interesting as it would give a privileged access to the universe conditions
before the CMB and BBN epochs, as neutrinos decoupled from the plasma as early
as $t\sim 1\,$s. For a long time, only upper limits on the PBH neutrino
background existed [108, 350].

The expected isotropic neutrino background was computed in the MG&W model by
Refs. [322, 351], with mere prospects for future detection. Halzen et al.
[351] pointed at the differences in the flux from the different neutrino
flavors, originating in the secondary PBH spectrum calculations: $\tau$
neutrinos for example are expected to have a lower flux because their
production involves heavier primary particles. Halzen et al. compared the PBH
neutrino flux to the known astrophysical backgrounds: while the Solar
background dominates by far at lower energies, the PBH intensity could be
above the atmospheric and supernova contributions at $E\gtrsim 10\,$MeV [352].
Furthermore, the Solar background is composed of electronic (anti)neutrinos
whereas PBHs produce them in all flavors.

Bugaev & Konishchev [353, 354] were the first to set a prospective constraint
on the PBH abundance with the expected Super-Kamiokande sensitivity. They
considered a scale-invariant PBH mass function with the critical collapse
mechanism and secondary particles from analytical fragmentation functions.
They also took into account precise optical depth of neutrinos.161616Just like
for photons, a photosphere model of PBH evaporation gives very different
constraints [330]. The first effective constraints were finally set with LSD
and KamLAND data [355]. A more robust neutrino constraint was computed by Carr
et al. [232], where they estimated the secondary neutrino fluxes from PBHs and
used the first Super-Kamiokande upper limits on the electronic antineutrino
[356].

Lunardini & Perez-Gonzalez [357] were the first to compute the neutrino
constraints from PBHs in the light of their Dirac/Majorana nature, “under the
hypothesis that black holes emit neutrino mass eigenstates” rather than flavor
ones [358]. Indeed, the oscillation mechanism [359, 360] proves that neutrinos
are massive: BHs “might be the only emitters of neutrino mass eigenstates”.
They predicted that only the flux of RHNs from early evaporated PBHs, once
redshifted, could lie above the Solar background at very low energy, in the
range of the PTOLEMY detector.

Dasgupta et al. [361] recomputed the constraints from $\overline{\nu}_{\rm e}$
from Super-Kamiokande using BlackHawk, hence numerical GFs for the first time
(with only tiny difference relative to [232]). They represented the HR of
neutrinos as mass eigenstates, so that the they don’t oscillate during
propagation. However, the secondary neutrinos oscillate as they are generated
as flavor eigenstates by EW interactions, but taking the oscillation matrix
into account results only in $\sim 2\%$ modification of the constraints [362].
They also considered log-normal PBH mass functions.

Until now, no direct detection of a neutrino background has been reported at
$E>10\,$MeV, with two consequences: first, the neutrino background from
astrophysical sources is still subject to debate; second, the PBH constraints
are only upper bounds. Wang et al. [362] computed the background neutrino
constraints at the prospective JUNO observatory, Calabrese et al. [363]
considered the XENONnT and DARWIN designs and de Romeri et al. [364] computed
the prospective DUNE ($\nu_{\rm e}$) and THEIA ($\overline{\nu}_{\rm e}$)
constraints. All these prospective constraints have been recomputed recently
by Bernal et al. [12] in a unified framework with the release of the public
numerical tool nuHawkHunter. They are also the first to estimate the
constraints down to neutrino decoupling epoch from $M\sim 10^{12}\,$g PBHs. In
Fig. 9, we show the old constraint from SK [232], with the more robust to date
from Ref. [12] as well as some prospective limits (HK, JUNO, DUNE, ARGO). We
see that the most stringent limits will be set by Super-Kamiokande.

Figure 9: Neutrino constraints on PBHs. The blue lines correspond to the
Super-Kamiokande limit from Ref. [232] (C10, dot-dashed) and [12] (B22,
solid). The dashed lines are derived from the prospective instruments Hyper-
Kamiokande (green), JUNO (red), DUNE (cyan) and ARGO (purple) [12]. The grey
shaded area is robustly excluded by Super-Kamiokande [12].

(RH)Neutrino emission has other consequences. If emitted before BBN/CMB, they
could participate in the effective number of neutrino species $\Delta N_{\rm
eff}$ [357]. Second, they act as an irreducible background in particle DM
direct searches, denoted as the “neutrino floor”, because they interact with
the nuclei and electrons composing DM detectors [363]. Third, they can scatter
on and boost DM particles [365]; this would modify the constraints in direct
detection experiments because boosted DM would interact more strongly with the
detectors (the resulting constraints on $f_{\rm PBH}$ completely depend on the
interaction cross-section between neutrinos and particle DM).

###  8.3 Gravitons

If there exists a gauge boson mediating the gravitational interaction, namely
the massless spin 2 graviton, then it should be emitted together with the rest
of the SM particles by PBH evaporation. In fact, the first PBH HR papers all
considered gravitons besides photons and neutrinos, as the massless set of
particles emitted by heavy BHs. This particle has in fact never been formally
identified and a coherent quantum theory of gravitation is yet to be found.
Thus, most of the later studies have deprecated the graviton emission by PBH
HR.

Vainer & Nasel’skii [237] seem to have been the first to predict a background
of gravitons generated by cosmological PBH evaporation. This background would
have very short wavelength, so that:

> […] the discovery of a additional peak in the spectrum of the primordial
> gravitational radiation [at $\sim 10\,$Å] would serve as an important
> argument favoring the existence of low-mass PBHs. [237]

This particular signature, unmodified since the dawn of the universe but
redshift effects (the gravitons have no optical depth), would be the only
reasonable observable on PBHs evaporated before BBN ($M\lesssim 10^{9}\,$g).
As a corollary, the graviton background could be the only signature of PBH
scenarios of baryogenesis (see e.g. [217]).

Alas, those gravitons would have had very high frequency $f\propto T(M\sim
10^{4}\,\text{g})$ at emission and the subsequent redshift leaves them with
still very high frequency $f\sim 10^{14}-10^{16}\,$Hz today [237, 366, 367].
Detection of these high-frequency GWs is a technical challenge, and some
proposals have been made in that direction [368]. Most of them rely on the
Gertsenshtein effect that converts high-frequency gravitons into photons in
the presence of a magnetic field [369] (GRAPH mixing). The magnetic field
could either be of astrophysical origin or monitored in a laboratory. The most
promising PBHs are those with $M\sim 10^{9}\,$g, just at the edge of the BBN
constraint, because their radiation is less redshifted than the smaller ones,
that evaporate earlier. The spectrum of GWs shows features that depend on the
expansion history of the universe.

The graviton background from PBHs was computed by several authors, including
Refs. [366, 370, 367], with numerical GFs used for the first time in [371]. In
fact, PBHs generate several GW backgrounds, e.g. from scalar-induced effects,
on top of HR gravitons. Inomata et al. [181] computed all these contributions
and point to errors in the previous studies in the redshift calculation. The
first experimental upper limits on high-frequency GWs were derived by Ejlli et
al. [372] in 2019 using data from photon cavities. While the bounds are far
from the expected backgrounds, this gives a proof of principle.

Gravitons injected in the very early universe have another consequence: as
relativistic particles, they participate in the effective number of neutrinos
$\Delta N_{\rm eff}$ as dark radiation (DR). Hence, their abundance is
constrained by BBN and CMB observations and this in turn would constrain the
abundance of PBHs in the $M\lesssim 10^{9}\,$g mass range. Hooper et al. [72]
gave the first calculation of the graviton impact on $\Delta N_{\rm eff}$ and
concluded that they would be constrained only by the future CMB-S4 experiments
(see also [373, 374]). Arbey et al. [73] showed that precision calculation of
the graviton density at the time of matter-radiation equality is necessary as
the constraints depend strongly on the assumptions regarding the reheating
temperature associated to PBH domination as well as the model for the
effective number of relativistic dofs $g_{*}(T)$ as a function of temperature.

###  8.4 Present backgrounds

PBHs radiate particles continuously, meaning that on top of the cosmologically
integrated background, there is also a contemporary background coming from
zero-redshift sources such as the Milky Way, nearby galaxies such as M31 or
DM-dominated objects like dSphs (Reticulum II, Leo T). The emission from
distant sources is easy to assess as they are observed as a whole, thus the
only free parameter is their DM density which enters as a normalization factor
in the PBH flux. The emission from the Milky Way is trickier as one does not
have a clear account of the DM density and spatial distribution in the Galaxy.
This is a pity as the Galactic center (GC) contribution to the global
background is two orders of magnitude larger than the “darkest” known objects
[341], and it forms an irreducible background to the isotropic surveys. Many
DM halo profiles have been proposed in the past decades, with the most popular
being the isothermal and Einasto profiles [375], and obviously the
Navarro–Frenck–White profile [376] (NFW). Some generalized version of this
profile is usually used in particle DM/PBH studies to obtain the column
density of DM along one particular line of sight (los). The flux of some
stable particle $i$ from PBHs in the Galaxy per unit solid angle is given by

$\dfrac{{\rm d}\Phi_{i}^{\rm gal}}{{\rm d}E}=\dfrac{1}{A_{\rm
gal}}\dfrac{J_{\rm gal}}{4\pi}\dfrac{{\rm d}^{2}N_{i}}{{\rm d}t{\rm d}E}\,,$
(28)

where

$J_{\rm gal}\equiv\dfrac{1}{\Delta\Omega}\int_{\Delta\Omega}{\rm
d}\Omega\int_{\rm los}\rho_{\rm gal}(r(l,\Omega))\,{\rm d}l\,,$ (29)

where $\Delta\Omega$ is the field of view and $A_{\rm gal}$ is the
normalization of the local PBH density. We see that uncertainties in the halo
profile propagate to the particle flux, with nearly 2 orders of magnitude
spread in the constraints for the $68\%$ CL parameters of the NFW profile
obtained with the Isatis tool for BlackHawk [342].

To our knowledge, the first evaluation of the Galactic PBH contribution to the
$\gamma$-ray background was given by Wright [377], claiming that the Galactic
EGRET survey [378] highlighted a PBH halo. Wright noted further that the
intensity of the EGXB and Galactic contributions were coincidentally of the
same order of magnitude in the MeV $-$ GeV range of EGRET, hence the Galactic
contribution cannot be neglected. It would appear as an anisotropic feature in
the diffuse background [379]. The Galactic contribution is usually taken into
account in the latest PBH EGXB studies [335, 339, 336, 338, 380], and some
papers even focused on this particular background [381, 382, 383, 341, 384],
e.g. with the high-resolution INTEGRAL/SPI data [385] that gives access to the
spatial distribution of $\gamma$-ray sources in the GC [386, 339, 380]. The
Galactic limit is usually the strongest one in the high-energy range as the
high-energy contribution from cosmological PBHs is diluted by the redshift.

The isotropic neutrino background also embeds a Galactic contribution. As the
observation of neutrinos is an all-sky average, only 1% difference is found
when tweaking the galactic profile [362]. Furthermore, contrarily to photons,
the isotropic component dominates because only low energy neutrinos $E\sim
10-50\,$MeV are surveyed [12].

Reticulum II is for now the only DM-dominated source that benefited from a
dedicated $\gamma$-ray survey in the context of particle DM/PBH studies [387],
while globular clusters could have captured numerous PBHs at the time of their
formation, resulting in bright sources inside the Galaxy [388].

##  9 Cosmic rays

Cosmic rays (CRs) are energetic particles that interact either in the upper
atmosphere, leading to cascades of secondary particles, or that are detected
by satellites. Within this review, CRs will be used as a denomination for
charged particles such as electrons/positrons, or protons/antiprotons and
heavier nuclei. The huge difficulty regarding CRs, as stated above, is that
they interact with the magnetic field of the Galaxy and therefore do not
propagate in straight lines. A second difficulty comes from the Solar winds
that hinder low-energy CRs from reaching the Earth, which acts as an effective
energy cut-off at $E\lesssim$ GeV. Both Galactic propagation [389] and Solar
winds suffer from large uncertainties: the Galactic magnetic field is poorly
known and re-acceleration mechanisms could exist [390, 391], while the Solar
winds vary periodically following two cycles of 11 and $22\,$yr. Combination
of the two effects is described by respectively a diffusion model and a
“force-field” approximation.171717Online material is available on the
PPPC4DMID website for propagation calculation. In the context of antiproton
detection, Wells et al. [392] claimed that should a probe escape the
heliopause, it would be ridden of the Solar modulation and have a direct
access to the local ISM flux of charged particles—this does not alleviate the
Galactic propagation uncertainties.

The complicated propagation of charged particles in the Galaxy results in a
segregation of the origin of the particles. Only CRs emitted inside the Galaxy
in recent epochs could contribute to the CR background in a sizeable
way.181818Carr et al. [232] propose a different constraint from antiprotons
that were produced before galaxy formation by evaporated $M\sim 10^{13}\,$g
PBHs. The CR limits further apply to local PBHs as the continuous energy
losses prevent charged particles from travelling efficiently from distant
regions of the Galaxy. The factors that can affect the CR spectrum are
annihilation on background particles, inverse Compton scattering on CMB
photons, Bremsstrahlung effects, e± pair production off nuclei, synchrotron
radiation due to ambient magnetic fields and ionization [322].

CRs provide a mean of direct detection of local PBHs, and the fact that the CR
limits were of the same order of magnitude as the $\gamma$-ray bounds in the
1980’s was interpreted as a correlative evidence for some PBH population in
the Galactic halo [393]. This was merely coincidental, and the constraints
have evolved a lot since then.

###  9.1 Electrons & positrons

Electrons and positrons are the simplest CRs as they are fundamental
particles. They could be released by atoms ionized in the ISM, or they could
originate in the evaporation of PBHs. The fact that CR measurements show that
the ratio $n_{{\rm e}^{+}}/n_{{\rm e}^{-}}\ll 1$, while PBHs should produce
both in equal quantity, advocates for use of the e+ flux as a limiting
constraint [394]. Electrons and positrons can be detected either directly by
satellite experiments, or they can produce radiation (e.g. synchrotron) that
allows for indirect detection in radio/X-rays. A particularly interesting
feature is the $511\,$keV line corresponding to e+-e- annihilation. Only PBHs
with mass $M\lesssim 10^{18}\,$g can produce sizeable amounts of e±, thus
limiting the scope of such constraints.

####  9.1.1 Direct detection

The first direct detection constraint on electrons and positrons is due to
Carr [108] using the $100\,$MeV background of charged particles [395, 396],
with the conclusion that this limit was not competitive with the strong
$\gamma$-ray one. The positron flux was recomputed within the MG&W model in
Ref. [322], as the PBH scenario was backed by evidence that $n_{{\rm
e}^{+}}/n_{{\rm e}^{-}}\longrightarrow 0.5$ at low energy. No more direct
detection constraints were set until Boudaud & Cirelli [397] used the fact
that Voyager 1 had crossed the heliopause and measured the ISM
electron/positron flux [398] to place new cleaner constraints on the local
density of PBHs, without the complicated Solar modulation, based on the DM
study [399]. The GO limit is used and only primary electrons are considered,
with uncertainties coming from the Galactic halo profile and the propagation
model (with or without re-acceleration). The constraints were also computed
for the first time with a log-normal mass function using the Carr method.

####  9.1.2 Indirect detection

Indirect detection of e± from PBHs can proceed from different mechanisms.
Charged particles loose energy when propagating in a magnetic field due to
synchrotron radiation, which is typically in the radio wavelengths. They can
also scatter off the CMB photons, which results in $\gamma$-rays in the
visible to near ultraviolet wavelengths [400]. Finally, the EM interactions
between e± and the ISM may heat up the gas clouds which would then radiate in
unexpected wavelengths [401]. Synchrotron radiation in the Galaxy due to PBH
evaporation was first predicted by Nasel’skii & Pelikhov [394], but was
completely left aside until recently. Chan & Lee [402] used archival radio
data of the GC to constrain the rate of injection of charged particles by
PBHs, relying on BlackHawk for the secondary electron spectrum; they
considered log-normal and power-law mass distributions. Dutta et al. [403]
studied the indirect detection of energetic electrons and positron emission by
PBHs through their inverse Compton scattering on the CMB photons, which would
produce a radio signal detectable at SKA observation of nearby dSphs (Segue I,
Ursa Major II, $\omega$-Centauri) and the Coma cluster. They considered only
primary e±. The constraints depend on the propagation model but could be
competitive for some targets. Mukhopadhyay et al. [400] computed the
synchrotron and inverse Compton scattering from e± injected by PBHs in the GC
with only the GO limit for HR used. Future radio instruments like SKA would
have enough sensitivity to constrain the PBH abundance very tightly by this
mean. Lee & Chan [404] extended this analysis to constrain the abundance of
PBHs in cool-core galaxy clusters, but this is not competitive.

Kim [401] were the first to compute the heating rate of the ISM by
electrons/positrons injected by PBHs in the DM-dominated object Leo T, using
BlackHawk. They advocated for the search for other DM dominated dSphs to
constrain further the PBH abundance. The same target was used by Laha et al.in
Ref. [405] where they extended the study to spinning PBHs and accounted for
photon heating. They found very different results compared to Ref.
[401]—constraints are not competitive.

####  9.1.3 511 keV line

A special consequence of PBH radiation of e± in a dense medium is the chance
that positrons scatter off astrophysical electrons. This can result in two
signals: a peaked signal at $511\,$keV if annihilation is immediate, or a low-
energy tail at $E<511\,$keV if positronium is formed. A $511\,$keV feature was
indeed detected some 50 years ago in the direction of the GC [406, 407, 408].
A PBH origin was first suggested by [409] and numerically estimated in [410,
411], while there may be other possible origins [412]. The $511\,$keV feature