# Constraining the primordial black hole abundance through Big-Bang
nucleosynthesis

Andrea Boccia<EMAIL_ADDRESS>Scuola Superiore Meridionale, Largo
S. Marcellino 10, I-80138 Napoli, Italy  Università degli Studi di Napoli
“Federico II”, Via Cintia 26, 80126 Napoli, Italy  Fabio Iocco
<EMAIL_ADDRESS>Università degli Studi di Napoli “Federico II”,
Via Cintia 26, 80126 Napoli, Italy  INFN sezione di Napoli, Via Cintia 1,
80126 Napoli, Italy  Luca Visinelli<EMAIL_ADDRESS>Tsung-Dao Lee
Institute (TDLI), No. 1 Lisuo Road, 201210 Shanghai, China  School of Physics
and Astronomy, Shanghai Jiao Tong University,
Dongchuan Road 800, 201240 Shanghai, China

###### Abstract

We investigate the scenario in which primordial black holes (PBHs) with masses
$\displaystyle M_{\rm PBH}\lesssim 10^{9}$ g undergo Hawking evaporation,
around the Big-Bang nucleosynthesis (BBN) epoch. The evaporation process
modifies the Universe’s expansion rate and the baryon-to-photon ratio, leading
to an alteration of the primordial abundance of light nuclei. We present
numerical solutions for the set of equations describing this physics,
considering different values of PBH masses and abundances at their formation,
showing how their evaporation impacts the abundances of light nuclei, obtained
by incorporating the non-standard Hubble rate and baryon-to-photon ratio into
the BBN code PArthENoPE. The results are then used to place upper bounds for
the PBH relative abundance at formation in the range $\displaystyle
10^{8}{\rm\,g}\lesssim M_{\rm PBH}\lesssim 10^{9}$ g, providing the strongest
constraints existing to-date in this mass range.

##  I Introduction

Primordial black holes (PBHs) are hypothetical objects formed in the very
early Universe via gravitational collapse of large density fluctuations in the
hot plasma Zel’dovich and Novikov (1967); Hawking (1971); Carr and Hawking
(1974). This occurs when an overdense region with an amplitude greater than a
critical value re-enters the horizon Carr (1975); Shibata and Sasaki (1999);
Niemeyer and Jedamzik (1999); Musco  _et al._ (2005). These large
inhomogeneities of order unity, produced through different mechanisms, allow
for the formation of PBHs with a mass of the same order as the total horizon
mass at that time. For this, more massive PBHs form at later times Carr  _et
al._ (2021).

Owing to the mechanism described above, PBHs can in principle cover a wide
mass spectrum, being suitable both as dark matter (DM) candidates Carr  _et
al._ (2016); Green and Kavanagh (2021) and as seeds for supermassive black
holes Duechting (2004). This last statement, however, applies only to those
PBHs which are stable over cosmological times, as it is commonly thought that
black holes should emit Hawking radiation and eventually evaporate due to
quantum effects Hawking (1975). In this framework, it is expected that PBHs
lighter than $\displaystyle M_{\rm PBH}\lesssim 10^{15}$ g would have
evaporated by now. Even though these light PBHs cannot solve the late-time
cosmology puzzles, it is still worth studying them as they could have had an
important role in the thermal history of the Universe, leaving signatures of
their existence onto early Universe-related cosmological observables.

For instance, the energy density of light PBHs formed very early on would
dilute at a lower rate than the surrounding relativistic plasma as the
Universe expands, triggering an epoch of early matter domination that would
last down to PBH evaporation time $\displaystyle t_{\rm evap}$. The condition
for this to happen is that the energy density of PBHs overcomes that of
radiation before $\displaystyle t_{\rm evap}$, condition which can be
expressed in terms of the abundance at formation of PBHs $\displaystyle\beta$
and the temperature of the Universe at their formation and evaporation,
respectively $\displaystyle T_{f}$ and $\displaystyle T_{\rm evap}$. An early
matter domination is achieved if $\displaystyle\beta$ exceeds the critical
value $\displaystyle\beta_{\rm crit}\equiv T_{\rm evap}/T_{f}$.

In this work, we focus on a particular class of light PBHs of mass
$\displaystyle M_{\rm PBH}\lesssim 10^{9}$ g that would have evaporated by the
time of Big-Bang nucleosynthesis (BBN). The non-thermal injection of high
energy particles from evaporating PBHs and their effects on nuclear reactions
have been widely discussed in previous literature Kohri and Yokoyama (2000);
Carr  _et al._ (2010), see Ref. Auffinger (2023) for a recent review. Here, we
consider a different aspect, namely the consequences of an expansion rate and
baryon–to–photon ratio, altered due to the presence of a primordial fraction
of PBHs. Even for values of the PBH abundance at formation below the threshold
value $\displaystyle\beta\leq\beta_{\rm crit}$, PBHs could have had an effect
on the expansion rate at the onset of nuclear reactions. In addition, the
spectrum of highly relativistic particles injected by the PBH evaporation into
the cosmological plasma also alters its entropy, thus modifying the baryon-to-
photon ratio. These two physical effects combine inducing an alteration of the
reaction rates, thus potentially modifying the primordial abundance of light
nuclei yields. We exploit this physics to derive new bounds on the PBH
parameter space $\displaystyle(\beta,M_{\rm PBH})$, which have been overlooked
in previous literature.

The numerical resolution of the set of coupled BBN equations is provided by
employing the latest version (v3.0) of the code Public Algorithm Evaluating
the Nucleosynthesis of Primordial Elements or PArthENoPE Pisanti  _et al._
(2008); Consiglio  _et al._ (2018); Gariazzo  _et al._ (2022).111The code is
publicly available at http://parthenope.na.infn.it. While we do not modify the
BBN reactions described in PArthENoPE, we input the background evolution of
the Hubble rate and the value of the baryon-to-photon ratio at $\displaystyle
T=10$ MeV, here $\displaystyle\eta_{i}$, into the code to account for the
modified dynamics in the presence of a PBH population.

This paper is organized as follows. In Sec. II we discuss the basic equations
and the methods used throughout the analysis. Results are presented in Sec.
III and discussed in Sec. IV. Finally, conclusions are drawn in Sec. V. We
work with units $\displaystyle\hbar=c=k_{B}=1$.

##  II Methods

###  II.1 Background cosmology

We consider an expanding Universe described by a flat Friedmann-Lemaître-
Robertson-Walker metric of scale factor $\displaystyle a=a(t)$, where the
cosmic time $\displaystyle t$ enters the definition of the Hubble rate
$\displaystyle H\equiv\dot{a}/a$ as a dot over any time-dependent quantity.
The Friedmann equation relating the Hubble rate with the total energy content
$\displaystyle\rho_{\rm tot}$ at time $\displaystyle t$ is $\displaystyle
3H^{2}=8\pi G\rho_{\rm tot}$. In the standard cosmological scenario, the
evolution of the Universe in its early stages is dominated by the energy
density of the relativistic component,

$\rho_{r}(T)=\frac{\pi^{2}}{30}g_{*}(T)T^{4}\,,$ (1)

with $\displaystyle g_{*}(T)$ being the number of relativistic degrees of
freedom at temperature $\displaystyle T$. As the Universe expands and the
scale factor $\displaystyle a$ increases, the relativistic component cools off
due to the relation $\displaystyle T\propto 1/a$ that holds for the adiabatic
regime.

The model described so far through Eq. (1) does not include the contribution
from non-relativistic matter, which becomes relevant at around matter-
radiation equality (“eq”) characterized by the redshift $\displaystyle
1+z_{\rm eq}\equiv\Omega_{m}/\Omega_{r}\sim 3400$. Here,
$\displaystyle\Omega_{m}$ and $\displaystyle\Omega_{r}$ are the present
abundances of matter and radiation, respectively. For redshifts $\displaystyle
z<z_{\rm eq}$, pressure-less matter overcomes the share in radiation in
governing the Hubble rate.

The picture above potentially modifies once a population of PBHs is included,
adding up to non-relativistic matter component that is present since the
earliest stages. The PBH fractional abundance at formation with respect to the
total energy density is

$\beta\equiv\frac{\rho_{\rm PBH}(t_{f})}{\rho_{\rm
tot}(T_{f})}\simeq\frac{\rho_{\rm PBH}(t_{f})}{\rho_{r}(T_{f})}\,,$ (2)

where $\displaystyle T_{f}$ is the temperature corresponding to time
$\displaystyle t_{f}$ and the second equality holds for a radiation-dominated
cosmology. Assuming that the PBHs have a monochromatic mass distribution
centered at the value $\displaystyle M_{\rm PBH}$, their energy density at any
time $\displaystyle t>t_{f}$ is found as $\displaystyle\rho_{\rm
PBH}(t)=M_{\rm PBH}n_{\rm PBH}(t)$, where $\displaystyle n_{\rm PBH}(t)$ is
the number density. Moreover, if the primordial power spectrum has a Gaussian
shape, the initial PBH abundance is

$\beta\sim\operatorname{erfc}\left[\frac{\delta_{c}}{\sqrt{2}\delta_{\rm
rms}}\right]\,,$ (3)

where $\displaystyle\delta_{\rm rms}$ the root mean square amplitude of the
primordial power spectrum. The Hubble rate at formation is related to the PBH
mass as Carr  _et al._ (2021)

$H_{f}=\gamma/(2GM_{\rm PBH})\,,$ (4)

when the horizon mass collapses into the PBH by a fraction
$\displaystyle\gamma$, which is $\displaystyle\gamma\sim 0.2$ for a radiation-
dominated scenario Carr (1975). The PBH energy density at formation is then

$\rho_{\rm PBH}(t_{f})=\frac{3\beta\gamma^{2}}{32\pi\,G^{3}\,M_{\rm
PBH}^{2}}\,.$ (5)

Given the uncertainty in $\displaystyle\gamma$ and the number of degrees of
freedom at PBH formation, $\displaystyle g_{*f}$, the rescaled abundance is
introduced as

$\beta^{\prime}\equiv\gamma^{1/2}\,\left(\frac{g_{*f}}{106.75}\right)^{-1/4}\,\beta\,,$
(6)

so that the results in previous literature can be compared Carr  _et al._
(2010). Strictly speaking, the scaling in Eq. (6) allows to express a
combination of $\displaystyle\beta$ and $\displaystyle\gamma$ as a unique
parameter only for models of PBHs that have not yet evaporated, and it is
adopted for models with lighter black holes to draw a comparison of the
constraints. For this, we also adopt it for this work, setting the reduced
Hubble constant $\displaystyle h=0.674$. Nevertheless, our results explicitly
depend on the choice of the parameters $\displaystyle\beta$ and
$\displaystyle\gamma$ which enter the initial condition for solving the
coupled Boltzmann equation, as discussed below. If the cosmic expansion after
PBH formation is adiabatic, the PBH number density scales with the entropy
density and it is given by

$n_{\rm PBH}(t)=\beta\,\frac{T_{f}}{T}\,\frac{\rho_{r}(T)}{M_{\rm PBH}}\,.$
(7)

In this scenario, the PBH component modifies the Hubble expansion rate and
might even trigger an early matter-dominated epoch if the energy density in
PBHs prevails over that in radiation. Using Eq. (7), this occurs after the
hypothetical early matter-radiation equality (“Eeq”), corresponding to
temperatures below $\displaystyle T\lesssim T_{\rm Eeq}\equiv\beta T_{f}$.

An early matter-dominated period is avoided if the evaporation of PBHs occurs
prior the hypothetical early matter-radiation equality. The spectrum of
particles radiated by an evaporating PBH peaks at the Hawking temperature
Hawking (1975)

$T_{H}=\frac{1}{8\pi GM_{\rm PBH}}\,.$ (8)

For a non rotating and charge-less PBH Mirbabayi  _et al._ (2020); De Luca
_et al._ (2019), the rate of evaporation is Page and Hawking (1976)

$\dot{M}_{\rm PBH}=-\frac{\mathcal{G}\,g_{H}}{30720\pi G^{2}M_{\rm
PBH}^{2}}\,.$ (9)

Here, $\displaystyle\mathcal{G}\sim 3.8$ is a graybody factor that accounts
for deviations from a black body spectrum due to the scattering of the emitted
radiation on the gravitational field of the PBH, and $\displaystyle
g_{H}=g_{H}(T_{H})$ is the spin-weighted number of degrees of freedom at
temperature $\displaystyle T_{H}$. Summing all contributions from SM particles
gives $\displaystyle g_{H}\simeq 102.6$ Mazde and Visinelli (2023), including
the effect of mutual repulsion between charged particles Page (1977). The
corresponding evaporation time obtained from integrating Eq. (9) is

$t_{\rm evap}=10240\pi G^{2}M_{\rm PBH}^{3}/(\mathcal{G}\,g_{H})\,.$ (10)

The condition over which an early matter-dominated epoch is achieved for some
period in terms of $\displaystyle\beta$ reads

$\beta\gtrsim\beta_{\rm crit}\equiv T_{\rm evap}/T_{f}\,,$ (11)

with $\displaystyle T_{\rm evap}$ being the temperature of the plasma at PBH
evaporation. In the parameter space we consider, a matter-dominated era never
occurs, the Universe always being radiation dominated, although the actual
amount of radiation present is altered with respect to the standard case by
the PBH evaporation.

The time evolution of the energy densities in radiation and PBHs is described
by a set of coupled Boltzmann equations in an expanding Universe. For the PBH
mass range considered here, the Hawking temperature is $\displaystyle
T_{H}\sim\mathcal{O}(10\,\rm{TeV})$ so that all SM products of the evaporation
are in the form of radiation. For this, the continuity equations for
$\displaystyle\rho_{r}$ and $\displaystyle\rho_{\rm PBH}$ have the bulk
expressions

$\displaystyle\displaystyle\dot{\rho}_{\rm PBH}+3H\rho_{\rm PBH}$
$\displaystyle\displaystyle=$ $\displaystyle\displaystyle\frac{\dot{M}_{\rm
PBH}}{M_{\rm PBH}}\,\rho_{\rm PBH}\,,$ (12)
$\displaystyle\displaystyle\dot{\rho}_{r}+4H\rho_{r}$
$\displaystyle\displaystyle=$ $\displaystyle\displaystyle-\frac{\dot{M}_{\rm
PBH}}{M_{\rm PBH}}\,\rho_{\rm PBH}\,.$ (13)

Here, the term proportional to the Hubble rate gives the cosmic dilution,
while the mass-loss rate is given in Eq. (9). Both the energy densities
contribute to the Hubble rate as

$H^{2}=\frac{8\pi G}{3}\left(\rho_{r}+\rho_{\rm PBH}\right)\,.$ (14)

For a given temperature above $\displaystyle T_{\rm evap}$, the expansion rate
in the modified history is then higher than in the standard scenario because
of the presence of the PBH component. Combining Eqs. (12) and (13) with the
Hubble rate in Eq. (14), gives the evolution of the PBH density,

$\frac{{\rm d}\rho_{\rm PBH}}{{\rm d}T}-\frac{3}{T}\frac{H}{H_{\rm
eff}}\rho_{\rm PBH}=-\frac{\dot{M}_{\rm PBH}}{M_{\rm PBH}}\frac{\rho_{\rm
PBH}}{TH_{\rm eff}}\,,$ (15)

where the effective Hubble rate is

$H_{\rm eff}=H+\frac{\dot{M}_{\rm PBH}}{M_{\rm PBH}}\frac{\rho_{\rm
PBH}}{4\rho_{r}}\,.$ (16)

Eq. (15) accounts for the evaporation of the PBH population and the energy
injection into the radiation component through the function $\displaystyle
H_{\rm eff}$. When considering the effects of the black holes onto the
expansion rate, we solve the equation with the initial condition in Eq. (2).
One key expression related to Eq. (13) and used in the derivations below
concerns the evolution of the entropy density $\displaystyle s(T)$, which is
here given as

$\dot{s}+3H_{\rm eff}\,s=0\,.$ (17)

###  II.2 Abundances of light nuclei

The origin of light elements during BBN and the existence of the CMB are two
successful predictions of the cosmological model, which require the Universe
to be hot enough for nuclear reactions to occur in its first few minutes. In
this, BBN has been an exceptionally successful theory, describing the
formation of light nuclei up to 7Be in the hot early Universe. It can be
considered today as a parameter free theory, insofar the baryon content at
present time $\displaystyle\Omega_{b}h^{2}=0.0224\pm 0.0001$ and the effective
number of relativistic species $\displaystyle N_{\rm eff}=2.99\pm 0.17$ –the
most relevant cosmological quantities in BBN– can in principle be determined
through the CMB observations in combination with baryon acoustic oscillation
measurements Aghanim  _et al._ (2020), and used to make predictions about the
primordial light nuclei yields in the specific cosmology considered.

Phenomenologically, BBN can be considered to begin when the weak reactions

$\displaystyle\displaystyle p+e^{-}$
$\displaystyle\displaystyle\rightleftharpoons n+\nu_{e}$ (18)
$\displaystyle\displaystyle n+e^{+}$
$\displaystyle\displaystyle\rightleftharpoons p+\bar{\nu}_{e}$ (19)

fall out of equilibrium at $\displaystyle T\sim 0.7$ MeV. The neutron-to-
proton ratio at the onset of nuclear reactions $\displaystyle n/p$ freezes out
except for the decaying free neutrons. At this point the temperature is low
enough to allow for the out–of–equilibrium formation of light elements such as
4He, D, 3He, up to 7Be and 7Li, the latter synthesized at the very cold end of
BBN, through a buildup of nuclei incapable of assembling heavier elements
Iocco  _et al._ (2007) before nuclear reactions finally freeze. Predictions of
the primordial abundances of the light elements are in overall good agreement
with those inferred from observational data Schramm and Turner (1998);
Steigman (2007); Iocco  _et al._ (2009); Cyburt  _et al._ (2016). For this,
BBN has been invoked to firmly constrain deviations from the standard
cosmological model below the MeV scale Pospelov and Pradler (2010).

In particular, after the chain of nuclear reactions has started a large
fraction of the free neutrons available end up forming helium-4 (4He). This is
expected since 4He is the nucleus with the highest binding energy between the
lightest elements formed during BBN. The primordial mass fraction of 4He is
then expected to amount to around

$Y_{p}=\frac{2(n/p)}{1+n/p}\simeq 0.25\,.$ (20)

On the observational side, the helium-4 mass fraction $\displaystyle Y_{p}$ is
inferred through the emission lines of H and He atoms in low-metallicity
extragalactic HII regions. According to the most recent review by the Particle
Data Group (PDG), the value for the 4He mass fraction deduced from these
measurements is $\displaystyle Y_{p}=0.245\pm 0.003$ Workman  _et al._ (2022).
A smaller fraction of free neutrons eventually ends up into other elements and
isotopes, while the fractions of D and 3He relative to H stand at the level of
$\displaystyle\mathcal{O}(10^{-5})$, and 7Li at the level of about
$\displaystyle\mathcal{O}(10^{-10})$. Theoretical arguments for these
fractions are found in Refs. Esmailzadeh  _et al._ (1991); Mukhanov (2004).
The fractional abundance of deuterium is measured in gas clouds seen in
absorption against an unrelated background source such as a quasar Adams
(1976). The weighted mean of the 11 most precise measurements is
$\displaystyle X_{D}=(2.547\pm 0.025)\times 10^{-5}$ Workman  _et al._ (2022)

The experimental values are reported in Table 1 in comparison with the results
from PArthENoPE v3.0 Gariazzo  _et al._ (2022). For the choice of main
cosmological parameters mentioned in Sec. III the code produces results
consistent with observations in the standard scenario without PBHs.

| $\displaystyle X_{D}\times 10^{5}$ | $\displaystyle Y_{p}$
---|---|---
Observations | $\displaystyle 2.547\pm 0.025$ | $\displaystyle 0.245\pm 0.003$
PArthENoPE | $\displaystyle 2.511\pm 0.1$ | $\displaystyle 0.246\pm 0.002$
Table 1: The weighted value for the observational abundances of deuterium and
4He from the PDG Workman  _et al._ (2022), compared with the values computed
from PArthENoPE v3.0 Gariazzo  _et al._ (2022).

In an expanding Universe, the freeze-out of chemical reactions occurs when the
reaction rate drops below the expansion rate and equilibrium can no longer be
achieved. For this, even a slight modification of the Hubble rate inevitably
affects the equilibrium of the ongoing reactions. In the presence of a sizable
PBH population, the expansion rate of the Universe can be increased, leading
to an earlier freeze-out of the chemical species. This results in an
overproduction of neutrons which are the most energetically disfavored species
in the reactions in Eqs. (18) and (19). Moreover, as the temperature rapidly
drops, less neutrons will have decayed at the onset of nuclear reactions.
Consequently, the $\displaystyle n/p$ ratio and the mass fraction of 4He are
expected to increase if PBHs in the mass range considered evaporates at the
onset of BBN.

The quantity most affecting the reactions of BBN is the baryon-to-photon ratio
$\displaystyle\eta=n_{b}/n_{\gamma}$, where $\displaystyle n_{b}$ is the
number of baryons and $\displaystyle n_{\gamma}$ is the number density of
photons, per comoving volume. This ratio is constant in the standard cosmology
in absence of interactions that change the number of particles, which is the
case at the temperatures considered, at around and below weak freeze-out. The
baryon-to-photon ratio is assessed at recombination to the value determined by
the CMB analysis reported by the Planck collaboration Aghanim  _et al._ (2020)

$\eta_{b}=6.138\times 10^{-10}\,,$ (21)

which is the default value used in the analysis by the code Parthenope since
there is no change in this quantity between BBN and recombination in standard
cosmology, nor in the case considered, since the PBHs evaporate before BBN,
and the material injected becomes part of the thermalized cosmological plasma
as described.

In the standard scenario, the value of the baryon-to-photon ratio changes
between the beginning and the end of BBN due to the transfer of entropy to
photons from $\displaystyle e^{+}e^{-}$ annihilation at $\displaystyle T\sim
0.5$ MeV. The value of $\displaystyle\eta$ before the onset of BBN,
$\displaystyle\eta_{i}$, is then obtained from its value at recombination as
in Ref. Serpico  _et al._ (2004)

$\eta_{i}=\frac{g_{s}(T_{i})}{g_{s}(T_{\rm rec})}\,\eta_{b}\simeq
2.75\,\eta_{b}\,,$ (22)

where $\displaystyle T_{\rm rec}\sim\,$eV and $\displaystyle T_{i}$ is
conventionally set in the code PArthENoPE to 10 MeV.

In the scenario we explore, PBHs evaporating around the time of $\displaystyle
e^{+}e^{-}$ annihilation constitute an additional source of entropy in the
plasma, and it is indeed crucial to estimate their effect on the alteration of
the baryon to photon ratio. To do so it is useful to express
$\displaystyle\eta$ in terms of the entropy density,

$\eta\propto g_{s}(T)\frac{n_{b}}{s}\,,$ (23)

we neglect the baryon contribution from PBHs evaporation, so that the
evolution of $\displaystyle n_{b}$ is given by

$\dot{n_{b}}+3Hn_{b}=0\,.$ (24)

Combining Eqs. (23) and (24) with Eq. (17) we find

$\dot{\eta}=\frac{3}{4}\frac{\dot{M}_{\rm PBH}}{M_{\rm PBH}}\frac{\rho_{\rm
PBH}}{\rho_{r}}\,\eta\,,$ (25)

which can be numerically solved together with Eq. (15) to obtain the evolution
of $\displaystyle\eta$ for different values of the parameters
$\displaystyle(\beta,M_{\rm PBH})$. The deuterium abundance is directly
related to the baryon-to-photon ratio $\displaystyle\eta$ so that a change in
$\displaystyle\eta$ greatly affects the fraction $\displaystyle X_{D}$.
Deuterium formation begins at temperatures much lower than its binding energy
$\displaystyle\Delta_{D}\simeq 2.2$ MeV due to dissociation by high energy
photons. The condition for deuterium to form is that the number of photons
with an energy above $\displaystyle\Delta_{D}$, per baryon, which is roughly
$\displaystyle\eta^{-1}e^{-\Delta_{D}/T}$, falls below unity. In the standard
scenario this happens at $\displaystyle T\sim 0.1$ MeV. In the presence of a
higher baryon-to-photon ratio, deuterium production is anticipated compared to
the standard case and so is the production of all the other nuclei, being the
deuterium necessary to start the chain of reactions. An earlier production of
deuterium leads to its more efficient conversion into heavier elements,
ultimately resulting in a lower primordial abundance with respect to the
standard case. Conversely, a lower baryon-to-photon ratio induces a later
onset of deuterium synthesis, a less efficient build-up of heavier elements
and consequently a higher deuterium abundance at the end of BBN.

Given the same initial conditions, and ab initio baryon-to-photon ratio, a
Universe filled with PBH will experience a reduced baryon-to-photon ratio
(with respect to the standard case) after their evaporation. This is because
the higher entropy per baryon of the particles injected lowers
$\displaystyle\eta$. In this sense, it can be said that PBH evaporation would
increase the entropy of the Universe and hence induce a higher deuterium
abundance. If one adopts the opposite point of view and imposes instead that
the value of $\displaystyle\eta$ at the end of BBN is the same as what is
measured through the CMB –which must be the case if no exotic physics takes
place in the Universe after PBH evaporation– then the baryon-to-photon ratio
of the Universe pre–existing the PBH evaporation must be higher than in the
standard case (the intervening PBH evaporation lowers it to the CMB value).
This in turn forces an earlier start of the deuterium nucleosynthesis, which
in turn leads to a lower final abundance. This is the approach we adopt here,
and in Figs. 1 and 2 we show the BBN abundance of deuterium and helium-4 as a
function of increasing PBH abundance at formation. The effect on deuterium
described above is what drives the constraints we derive222We have also tested
the effect of modification of the Hubble rate, described at the beginning of
this section, by artificially (and unphysically) switching on only one of the
two at the time: as intuition would have it we do see that the Hubble rate
modification alone has a leading effect on helium, whereas the baryon-to-
photon ratio affects both helium and deuterium, with the latter undergoing the
most dramatic modification and thus leading the resulting constraints..

##  III Results

BBN is a strongly non-linear process in which multiple reactions take place at
the same time. For this, numerical methods for the integration of coupled
differential equations are required in order to obtain an estimate of its
final outcome. For such calculations we employ here the numerical code
PArthENoPE v3.0 Gariazzo  _et al._ (2022), which allows us to consider up to
26 nuclides and 100 reactions to solve for the whole system of differential
equations and compute reliable theoretical predictions for the light nuclei
abundances. While the set of BBN reactions solved by the code PArthENoPE is
not modified, we account for the altered expansion history of the Universe by
solving for Eqs. (15) and (25) numerically for given values of the parameters
$\displaystyle\beta$ and $\displaystyle M_{\rm PBH}$ and with the initial
condition in Eq. (2), following the steps detailed in previous work Masina
(2020); Mazde and Visinelli (2023). The resulting Hubble rate $\displaystyle
H(T)$ as given in Eq. (14) is then treated as an input for the background
evolution over which the coupled BBN equations are solved with the FORTRAN77
source code.

The value of the baryon-to-photon ratio at the onset of BBN is related to its
value at recombination through Eq. (22) in PArthENoPE. We account for the
modified evolution of $\displaystyle\eta$ due to the entropy injection by PBHs
by changing the default value in the source code to the value of
$\displaystyle\eta_{i}$ obtained from the solution of Eq. (25) with the
condition $\displaystyle\eta=\eta_{b}$ at recombination. The distribution of
the PBH masses is assumed to be monochromatic and centered around the value
$\displaystyle M_{\rm PBH}$.

Figure 1 shows the primordial mass fraction of 4He –$\displaystyle Y_{p}$ as
customary in literature– as a function of the parameter
$\displaystyle\beta^{\prime}$ in Eq. (6). The result is obtained assuming a
population of PBHs with a monochromatic mass distribution centered around
M${}_{\rm PBH}=5\times 10^{8}$ g. As the relative population of PBHs increases
with increasing values of $\displaystyle\beta^{\prime}$, the mass fraction
$\displaystyle Y_{p}$ starts growing significantly due to the modified
cosmology impacting around the onset of BBN. We have explored the region of
parameters around the critical value $\displaystyle\beta^{\prime}_{\rm crit}$
related to the quantity in Eq. (11) and marked by the vertical dotted line.
Once the value of $\displaystyle Y_{p}$ for given
$\displaystyle\beta^{\prime}$ is determined through PArthENoPE, the error bar
is associated assuming a fractional error of 1% Workman  _et al._ (2022). In
Fig. 1, the red band gives the 1$\displaystyle\sigma$ range associated with
the mass fraction from PDG, see Table 1. Similarly, the fraction of deuterium
$\displaystyle X_{D}$ is shown in Fig. 2 as a function of the PBH primordial
abundance $\displaystyle\beta^{\prime}$. As explained below Eq. (25), an
increase in the PBH abundance leads to a decrease in $\displaystyle X_{D}$ due
to the anticipated production of deuterium for any given PBH mass.

As described in the previous section, the deuterium fraction is more sensible
to the increasing $\displaystyle\beta^{\prime}$ with respect to helium, so it
provides more stringent bounds.

Using the methods described above, we compare the results for the primordial
abundance of light elements, obtained for different values of
$\displaystyle(\beta,M_{\rm PBH})$, with the observational constraints on the
abundances of the elements 4He and D. As we discussed, helium and deuterium
fractions are extremely sensible to the physical alterations introduced at
BBN. The code PArthENoPE v3.0 is run across the range of values $\displaystyle
10^{8}{\rm\,g}\leq M_{\rm PBH}\leq 10^{9}{\rm\,g}$ and $\displaystyle
10^{-18}\leq\beta\leq 10^{-14}$. For a given value of the PBH mass population,
the bound on $\displaystyle\beta^{\prime}$ is obtained by comparing the value
for $\displaystyle Y_{p}$ and $\displaystyle X_{D}$ obtained experimentally,
red band in Fig. 1 and Fig. 2, with those determined through PArthENoPE v3.0.
We request that our predictions from the model does not deviate from what is
reported by PDG by more than 2 sigma or 95% confidence level or CL, where the
uncertainties for the values of $\displaystyle Y_{p}$ and $\displaystyle
X_{D}$ are inferred from PArthENoPE, and reported by observations are added in
squares.

Figure 1: Primordial mass fraction of 4He for different values of the
parameter $\displaystyle\beta^{\prime}$, for a population of PBHs with
M${}_{\rm PBH}=5\times 10^{8}$ g and a monochromatic mass distribution. Also
shown is the measurement for $\displaystyle Y_{p}$ reported by the PDG Workman
_et al._ (2022), see the red band giving the 1$\displaystyle\sigma$ range.
Figure 2: Primordial abundance of deuterium for different values of the
parameter $\displaystyle\beta^{\prime}$, for a population of PBHs with
M${}_{\rm PBH}=5\times 10^{8}$ g. Also shown is the measurement for
$\displaystyle X_{D}$ reported by the PDG Workman  _et al._ (2022), see the
red band giving the 1$\displaystyle\sigma$ range.

We do not include other light elements produced during BBN in our analysis, as
their primordial abundances are suppressed with respect to deuterium and
helium, and they are further affected by the subsequent stellar evolution,
which contributes to additional uncertainties in the determination of their
primordial abundances and are thus not to be used for cosmology purposes Iocco
_et al._ (2009); Iocco (2012). The default value for the baryon-to-photon
ratio at recombination used in the analysis is given in Eq. (21), while the
neutron lifetime is set at the value $\displaystyle\tau_{n}=(879.4\pm 0.6)$ s
Zyla  _et al._ (2020).

Figure 3 returns the upper bounds on $\displaystyle\beta^{\prime}(M)$ at 95%
CL obtained from our methods of considering the effects of a modified
expansion rate by the presence of PBHs at the onset of BBN. This is shown by
the area in orange shade labeled “BIV24 (This work)”. We have considered the
effects resulting from the modifications in the helium abundance alone (orange
shaded area) and deuterium (okra shaded area). PBHs of mass $\displaystyle
M_{\rm PBH}\lesssim 10^{8}$ g evaporate too early to affect the evolution of
BBN.

On the opposite end, PBHs heavier than about $\displaystyle 10^{9}$ g would
evaporate well into the BBN epoch, so that the by-products of their
evaporation could affect the nuclear reactions. While we do not consider these
effects in our analysis, much tighter constraints can be drawn from including
these by-products as discussed in Ref. Carr  _et al._ (2010), see the gray
shaded area labeled “Carr+09”.

Figure 3: The bounds at 95% confidence level on PBHs in the mass range
$\displaystyle 10^{8}{\rm\,g}\leq M_{\rm PBH}\leq 10^{9}{\rm\,g}$ for the
helium abundance (orange shaded area) and deuterium (okra shaded area). Also
shown are the results from Ref. Carr  _et al._ (2010) (gray shaded area).

##  IV Discussion

By accounting for the physics described above, we can constrain the abundance
of PBHs at their formation $\displaystyle\beta$ in the mass range
$\displaystyle 10^{8}{\rm\,g}\lesssim M_{\rm PBH}\lesssim 10^{9}{\rm\,g}$. In
the mass region of interest, our results improve the bounds on
$\displaystyle\beta$ by several orders of magnitude with respect to previous
analyses present in the literature, which account for the secondary
gravitational wave produced after PBH evaporation and their effects on the
cosmic microwave background Domènech  _et al._ (2021). Given that for PBH
masses considered here the evaporation process ends well before the formation
of nuclei, we have deliberately omitted the non-thermal effects of high energy
particles injection on nucleosynthesis. Indeed, we have estimated that all
high energy particles emitted by the evaporation thermalize “on–the–spot”,
namely within a short time enough to be considered instantaneous at the moment
(and therefore cosmic temperature) of evaporation, which for the masses
considered is before the very onset of BBN, as recalled above. This avoids any
dissociation of nuclei by highly energetic photons and hadrons, and allows us
to remain within the well-based approximation of thermal nucleosynthesis.

In the analysis, we have neglected the effects brought in by the flux of
neutrinos emitted by PBHs before or during neutrino decoupling, which leads to
a modification of the abundance of deuterium Vainer and Naselskii (1978). We
also expect a change in the effective number of relativistic neutrinos
$\displaystyle N_{\rm eff}$ and indirectly on $\displaystyle Y_{p}$, leading
to another complementary effect in addition to our results Mastrototaro  _et
al._ (2021). In such a scenario, our results should be regarded as
conservative. These aspects will be treated in future work.

Furthermore, our results are derived under the assumption of a monochromatic
mass spectrum for the PBH distribution. However, the mass distribution could
span various decades of mass as it is expected in models where PBH formation
possesses a low mass tail from critical collapse Niemeyer and Jedamzik (1998).
Accounting for an extended mass distribution is generally non trivial and
requires both an improvement of the methods presented in Sec. II and further
investigations on the distribution to be used. Additional phenomenology, such
as accretion and merging of PBHs, may arise leading to the necessity of a
larger set of parameters. In any case with a non-monochromatic mass
distribution we expect the impact on the Hubble rate to be weaker but to last
longer due to the higher masses tail of the distribution. This may affect our
constraints loosening them, simultaneously broadening them across a wider mass
spectrum, particularly towards lower masses. The last statement applies in the
case of a narrow mass distribution, for which the definition of
$\displaystyle\beta$ can be obtained by integrating the explicit mass
distribution over the whole mass range. For a wider mass distribution, the
time dependence of the width of the distribution should be taken into account
due to the different evaporation times between the lower and higher masses.
The approach usually adopted for this case is the binning of the distribution
into intervals of size $\displaystyle\Delta M_{\rm PBH}\sim M_{\rm PBH}$ which
are treated separately.

##  V Conclusions

We examined the effects of primordial black holes as they undergo evaporation
during the weak freeze-out phase at the onset of Primordial Nucleosynthesis
(BBN). By solving the differential equations governing PBHs and radiation
density evolution alongside the Friedmann equations, we were able to model the
evolution of the background cosmology in the presence of PBHs of various
masses and abundance at formation, assuming a monochromatic mass spectrum. We
have accounted for the change in the baryon-to-photon ratio and Hubble rate
that are brought forth by the PBH evaporation. We have implemented these
results into the BBN numerical code PArthENoPE, enabling us to compute the
primordial abundance of light nuclei as a function of the parameters
$\displaystyle(\beta,M_{\rm PBH})$.

Our analysis shows significant deviations in the abundance of deuterium and
4He from the standard, no–PBH case, for values of $\displaystyle\beta$
approaching the critical threshold $\displaystyle\beta_{\rm crit}\equiv T_{\rm
evap}/T_{f}$, for PBHs masses within the range $\displaystyle
10^{8}{\rm\,g}\leq M_{\rm PBH}\leq 10^{9}{\rm\,g}$. By comparing the
abundances of deuterium and 4He theoretically predicted in each PBH scenario
against the empirically observed primordial abundances, we have derived the
constraints in the PBH parameters. Our results request the parameter
$\displaystyle\beta^{\prime}$ to be lower than the critical value
$\displaystyle\beta_{\rm crit}$ for the whole range of masses considered,
thereby ruling out the possibility of early PBH domination. The strongest
bound are obtained for the heaviest PBHs considered, with
$\displaystyle\beta^{\prime}\lesssim 10^{-16}$ for $\displaystyle M_{\rm
PBH}=10^{9}$ g.

This analysis yields novel constraints within the parameter space of PBHs,
more stringent with respect to previous results over the same mass range. The
bounds obtained with the method presented here could be strengthened by
considering the effects of high energy neutrinos emission before neutrino
decoupling that could impact $\displaystyle N_{\rm eff}$. Future work should
be devoted to a more complete treatment of these effects of evaporation
products and should extend the analysis to other mass distributions and
models.

###### Acknowledgements.

We thank O. Pisanti for help with the PArthENoPE code, P. D. Serpico for
useful comments on the manuscript, and M. Zantedeschi for the insight provided
on the phenomenology of PBHs evaporation. F.I. is partially supported by the
research grant No. 2022E2J4RK “PANTHEON: Perspectives in Astroparticle and
Neutrino THEory with Old and New messengers” under the program PRIN 2022
funded by the Italian Ministero dell’Università e della Ricerca (MUR) and by
the European Union – Next Generation EU. A.B. and F.I. acknowledge the
hospitality of the Tsung-Dao Lee Institute (TDLI) in Shanghai (China) and
Università di Ferrara (Italy) during the final phases of this work, funded by
the TAsP Iniziativa Specifica of INFN. L.V. acknowledges support by the
National Natural Science Foundation of China (NSFC) through the grant No.
12350610240 “Astrophysical Axion Laboratories”, as well as hospitality by the
Istituto Nazionale di Fisica Nucleare (INFN) Frascati National Laboratories,
the Galileo Galilei Institute for Theoretical Physics in Firenze (Italy), the
INFN section of Napoli (Italy), the INFN section of Trento (Italy), the INFN
section of Ferrara (Italy), and the University of Texas at Austin (USA)
throughout the completion of this work. This publication is based upon work
from the COST Actions “COSMIC WISPers” (CA21106) and “Addressing observational
tensions in cosmology with systematics and fundamental physics (CosmoVerse)”
(CA21136), both supported by COST (European Cooperation in Science and
Technology).

## References

  * Zel’dovich and Novikov (1967) Y. B. Zel’dovich and I. D. Novikov, Sov. Astron. 10, 602 (1967).
  * Hawking (1971) S. Hawking, Mon. Not. Roy. Astron. Soc. 152, 75 (1971).
  * Carr and Hawking (1974) B. J. Carr and S. W. Hawking, Mon. Not. Roy. Astron. Soc. 168, 399 (1974).
  * Carr (1975) B. J. Carr, Astrophys. J. 201, 1 (1975).
  * Shibata and Sasaki (1999) M. Shibata and M. Sasaki, Phys. Rev. D 60, 084002 (1999), arXiv:gr-qc/9905064 .
  * Niemeyer and Jedamzik (1999) J. C. Niemeyer and K. Jedamzik, Phys. Rev. D 59, 124013 (1999), arXiv:astro-ph/9901292 .
  * Musco  _et al._ (2005) I. Musco, J. C. Miller, and L. Rezzolla, Class. Quant. Grav. 22, 1405 (2005), arXiv:gr-qc/0412063 .
  * Carr  _et al._ (2021) B. Carr, F. Kuhnel, and L. Visinelli, Mon. Not. Roy. Astron. Soc. 501, 2029 (2021), arXiv:2008.08077 [astro-ph.CO] .
  * Carr  _et al._ (2016) B. Carr, F. Kuhnel, and M. Sandstad, Phys. Rev. D 94, 083504 (2016), arXiv:1607.06077 [astro-ph.CO] .
  * Green and Kavanagh (2021) A. M. Green and B. J. Kavanagh, J. Phys. G 48, 043001 (2021), arXiv:2007.10722 [astro-ph.CO] .
  * Duechting (2004) N. Duechting, Phys. Rev. D 70, 064015 (2004), arXiv:astro-ph/0406260 .
  * Hawking (1975) S. W. Hawking, Commun. Math. Phys. 43, 199 (1975), [Erratum: Commun.Math.Phys. 46, 206 (1976)].
  * Kohri and Yokoyama (2000) K. Kohri and J. Yokoyama, Phys. Rev. D 61, 023501 (2000), arXiv:astro-ph/9908160 .
  * Carr  _et al._ (2010) B. J. Carr, K. Kohri, Y. Sendouda, and J. Yokoyama, Phys. Rev. D 81, 104019 (2010), arXiv:0912.5297 [astro-ph.CO] .
  * Auffinger (2023) J. Auffinger, Prog. Part. Nucl. Phys. 131, 104040 (2023), arXiv:2206.02672 [astro-ph.CO] .
  * Pisanti  _et al._ (2008) O. Pisanti, A. Cirillo, S. Esposito, F. Iocco, G. Mangano, G. Miele, and P. D. Serpico, Comput. Phys. Commun. 178, 956 (2008), arXiv:0705.0290 [astro-ph] .
  * Consiglio  _et al._ (2018) R. Consiglio, P. F. de Salas, G. Mangano, G. Miele, S. Pastor, and O. Pisanti, Comput. Phys. Commun. 233, 237 (2018), arXiv:1712.04378 [astro-ph.CO] .
  * Gariazzo  _et al._ (2022) S. Gariazzo, P. F. de Salas, O. Pisanti, and R. Consiglio, Computer Physics Communications 271, 108205 (2022), arXiv:2103.05027 [astro-ph.IM] .
  * Mirbabayi  _et al._ (2020) M. Mirbabayi, A. Gruzinov, and J. Noreña, JCAP 03, 017 (2020), arXiv:1901.05963 [astro-ph.CO] .
  * De Luca  _et al._ (2019) V. De Luca, V. Desjacques, G. Franciolini, A. Malhotra, and A. Riotto, JCAP 05, 018 (2019), arXiv:1903.01179 [astro-ph.CO] .
  * Page and Hawking (1976) D. N. Page and S. W. Hawking, Astrophys. J. 206, 1 (1976).
  * Mazde and Visinelli (2023) K. Mazde and L. Visinelli, JCAP 01, 021 (2023), arXiv:2209.14307 [astro-ph.CO] .
  * Page (1977) D. N. Page, Phys. Rev. D 16, 2402 (1977).
  * Aghanim  _et al._ (2020) N. Aghanim _et al._ (Planck), Astron. Astrophys. 641, A6 (2020), [Erratum: Astron.Astrophys. 652, C4 (2021)], arXiv:1807.06209 [astro-ph.CO] .
  * Iocco  _et al._ (2007) F. Iocco, G. Mangano, G. Miele, O. Pisanti, and P. D. Serpico, Phys. Rev. D 75, 087304 (2007), arXiv:astro-ph/0702090 .
  * Schramm and Turner (1998) D. N. Schramm and M. S. Turner, Rev. Mod. Phys. 70, 303 (1998), arXiv:astro-ph/9706069 .
  * Steigman (2007) G. Steigman, Ann. Rev. Nucl. Part. Sci. 57, 463 (2007), arXiv:0712.1100 [astro-ph] .
  * Iocco  _et al._ (2009) F. Iocco, G. Mangano, G. Miele, O. Pisanti, and P. D. Serpico, Phys. Rept. 472, 1 (2009), arXiv:0809.0631 [astro-ph] .
  * Cyburt  _et al._ (2016) R. H. Cyburt, B. D. Fields, K. A. Olive, and T.-H. Yeh, Rev. Mod. Phys. 88, 015004 (2016), arXiv:1505.01076 [astro-ph.CO] .
  * Pospelov and Pradler (2010) M. Pospelov and J. Pradler, Ann. Rev. Nucl. Part. Sci. 60, 539 (2010), arXiv:1011.1054 [hep-ph] .
  * Workman  _et al._ (2022) R. L. Workman _et al._ (Particle Data Group), PTEP 2022, 083C01 (2022).
  * Esmailzadeh  _et al._ (1991) R. Esmailzadeh, G. D. Starkman, and S. Dimopoulos, Astrophys. J. 378, 504 (1991).
  * Mukhanov (2004) V. F. Mukhanov, Int. J. Theor. Phys. 43, 669 (2004), arXiv:astro-ph/0303073 .
  * Adams (1976) T. F. Adams, Astron. Astrophys. 50, 461 (1976).
  * Serpico  _et al._ (2004) P. D. Serpico, S. Esposito, F. Iocco, G. Mangano, G. Miele, and O. Pisanti, JCAP 12, 010 (2004), arXiv:astro-ph/0408076 .
  * Masina (2020) I. Masina, Eur. Phys. J. Plus 135, 552 (2020), arXiv:2004.04740 [hep-ph] .
  * Iocco (2012) F. Iocco, Mem. Soc. Astron. Ital. Suppl. 22, 19 (2012), arXiv:1206.2396 [astro-ph.GA] .
  * Zyla  _et al._ (2020) P. A. Zyla _et al._ (Particle Data Group), Progress of Theoretical and Experimental Physics 2020, 083C01 (2020), https://academic.oup.com/ptep/article-pdf/2020/8/083C01/34673722/ptaa104.pdf .
  * Domènech  _et al._ (2021) G. Domènech, V. Takhistov, and M. Sasaki, Phys. Lett. B 823, 136722 (2021), arXiv:2105.06816 [astro-ph.CO] .
  * Vainer and Naselskii (1978) B. V. Vainer and P. D. Naselskii, Soviet Astron. 22, 138 (1978).
  * Mastrototaro  _et al._ (2021) L. Mastrototaro, P. D. Serpico, A. Mirizzi, and N. Saviano, Phys. Rev. D 104, 016026 (2021), arXiv:2104.11752 [hep-ph] .
  * Niemeyer and Jedamzik (1998) J. C. Niemeyer and K. Jedamzik, Phys. Rev. Lett. 80, 5481 (1998), arXiv:astro-ph/9709072 .