Investigating the Impact of Primordial Black Hole Evaporation on the Lithium-7 Discrepancy in Big Bang Nucleosynthesis

By
Vaibhav Parashar
-
0
103
Read the full article as a PDFDownload

Abstract

Big Bang Nucleosynthesis (BBN) offers a powerful way to study the early universe, successfully predicting the original amounts of primordial light elements, however, it consistently overestimates the abundance of lithium-7 (^7Li) by a factor of 3-4, a persistent puzzle known as the cosmic lithium problem. In this study, we investigate whether the evaporation of primordial black holes (PBHs) via Hawking radiation, which can affect the ratio of neutrons to protons before they freeze out during BBN, might offer a solution. To test this scenario, we modify the publicly available AlterBBN code, which is frequently used for modelling Big Bang Nucleosynthesis under non-standard cosmologies, by implementing a treatment of energy injection from PBHs. This lets us model how nucleosynthesis proceeds across different PBH masses and initial abundances \beta_{\mathrm{PBH}}. We validate our results by comparing outputs with those from the Parthenope code for consistency. Our goal is to determine whether PBH driven energy injection can reduce lithium-7 production while preserving accurate predictions for deuterium and helium-4 abundances. However, our results indicate that for M_{\rm PBH} = 5 \times 10^{8}~\mathrm{g}, PBH-induced energy injection worsens the lithium problem, increasing the predicted ^7\mathrm{Li}/\mathrm{H} abundance by more than a factor of three relative to the standard BBN value for \beta_{\rm PBH} \gtrsim 10^{-15}. At the same time, deuterium becomes underproduced and helium-4 overproduced, violating observational bounds and placing an upper constraint on the allowed PBH abundance for this mass range.

Introduction

Big Bang Nucleosynthesis (BBN) is a cornerstone of modern cosmology, describing how the light elements were synthesized within the first few minutes of the universe. The predicted primordial abundances depend sensitively on the baryon-to-photon ratio, the number of neutrino species, and the expansion rate of the early universe1. Observations of deuterium (D) and helium (^4He) match these predictions with high precision, providing strong support for Standard BBN model.

In contrast, the predicted 7Li abundance 7Li/H \sim 5 \times 10^{-10} exceeds the measured value in old, metal-poor halo stars \sim 1.5 \times 10^{-10} by a factor of 3-4.This persistent inconsistency known as the cosmic lithiunom problem1,2 has resisted resolution for decades.

Proposed explanations include uncertainties in nuclear reaction rates3, stellar depletion4, or exotic new physics such as decaying dark matter5 and primordial black holes (PBHs)6.

In this study, we investigate whether energy injection from evaporating PBHs in the early universe could modify BBN and alleviate the lithium problem. We adapt the open-source AlterBBN code7 to incorporate PBH evaporation effects and simulate how different PBH masses and abundances influence primordial light element yields.

Theoretical Background

Big Bang Nucleosynthesis

Big Bang Nucleosynthesis (BBN) occurred during the first \sim 20 minutes after the Big Bang, as the universe cooled from temperatures of order MeV to keV. During this period, neutrons and protons combined through nuclear reactions to form light nuclei such as deuterium (D), helium-3 (^3\mathrm{He}), helium-4 (^4\mathrm{He}), and lithium-7 ^7\mathrm{Li}8,9,10. The predicted abundances of these elements depend strongly on the baryon-to-photon ratio (\eta), which is now precisely determined by measurements of the cosmic microwave background (CMB)11. The ^4\mathrm{He} mass fraction (Y_p) and the deuterium-to-hydrogen ratio (D/H) show excellent agreement between theoretical predictions and observations, providing a remarkable success of the Standard Big Bang Nucleosynthesis (SBBN) model9,12.

The Lithium Problem

Despite the overall success of SBBN, the predicted abundance of ^{7}Li remains a major unresolved issue.
Theoretical models yield (^{7}\mathrm{Li}/\mathrm{H}){\mathrm{SBBN}} \approx 5 \times 10^{-10}, while observations of the oldest, metal-poor halo stars, the so-called Spite plateau , give (^{7}\mathrm{Li}/\mathrm{H}){\mathrm{obs}} \approx 1.5 \times 10^{-10}, a discrepancy of about a factor of 3-413,2.

Attempts to resolve this problem by revising nuclear reaction rates, such as ^{3}\mathrm{He}(\alpha,\gamma)^{7}\mathrm{Be}, reduce the predicted ^{7}Li yield only marginally14.
Stellar depletion mechanisms, where ^{7}Li is destroyed in stellar interiors, have also been challenged by the observed flatness of the Spite plateau and the lack of sufficient depletion without also destroying D or ^{4}He15.

As a result, alternative possibilities involving physics beyond the Standard Model have been considered, including energy injection from exotic particles or primordial black holes (PBHs)5,6,7.

Primordial Black Holes

Primordial black holes (PBHs) are hypothetical black holes that could have formed in the very early universe, long before the first stars or galaxies. Unlike stellar black holes, which form from the gravitational collapse of massive stars at the end of their lifetimes, PBHs can form from the collapse of high-density regions during the radiation-dominated era, when the universe was only a fraction of a second old16,17.

Such overdensities could be seeded by quantum fluctuations during inflation or by phase transitions that temporarily enhanced the local density. If the density contrast \delta \rho / \rho in a region exceeded a critical threshold (typically \delta_c \sim 0.3—0.7), gravity would overcome pressure forces, and that region would collapse into a black hole6,18.

One unique feature of PBHs is their wide range of possible masses. The mass of a PBH is approximately equal to the mass contained within the cosmological horizon at the time it formed:

(1)   \begin{equation*}M_{\mathrm{PBH}} \approx 10^{15} \mathrm{g}\,\left(\frac{t}{10^{-23}\,\mathrm{s}}\right)\end{equation*}

where t is the cosmic time at formation19,18
This implies that PBHs could have formed with masses ranging from the Planck scale \sim 10^{-5} \mathrm{g} up to several solar masses, depending on the formation time.

For PBHs lighter than about 5 \times 10^{14} \mathrm{g}, Hawking radiation causes them to completely evaporate within the current age of the universe16. In particular, PBHs in the mass range 10^{8}—10^{11} \mathrm{g} have lifetimes on the order of 10^{2}—10^{3} seconds, which is comparable to the duration of BBN, hence, they could inject significant energy into the cosmic plasma during this critical period and modify the abundances of light elements20,21.

As PBHs evaporate, they emit high-energy particles including photons, electrons, neutrinos, and hadrons via Hawking radiation. This energy injection can modify the neutron-to-proton ratio, induce photodissociation of light nuclei, and alter the expansion rate of the universe. These effects could potentially reduce the overproduction of ^7Li while preserving the consistency of D and ^4He abundances with observations, making PBHs an intriguing candidate for resolving the cosmic lithium problem22,11.

Hawking Radiation and Energy Injection

Hawking radiation predicts that black holes emit particles thermally, with a temperature inversely proportional to their mass:

(2)   \begin{equation*}T_{\mathrm{BH}} = \frac{\hbar c^3}{8\pi G M k_B} \approx 1.06~\mathrm{TeV} \left( \frac{10^{10}~\mathrm{g}}{M} \right)\end{equation*}

Smaller PBHs evaporate more quickly, emitting photons, neutrinos, electrons, and quarks/gluons20,19. The lifetime of a PBH scales roughly as M^3, so PBHs of 10^8 - 10^{11}~\mathrm{g} would have evaporated within the first 10^3 seconds just like in Eq. 11. This energy injection can alter the neutron-to-proton ratio, modify nuclear reaction rates, or photodissociate light nuclei, potentially reducing the overproduction of ^7\mathrm{Li}5,20.

Previous studies have considered the impact of PBHs on BBN and derived constraints on their abundance6,23. However, a systematic investigation of whether PBH evaporation can resolve the lithium problem while preserving the D and ^4He predictions remains incomplete. Our study aims to fill this gap by modelling PBH energy injection in the AlterBBN24 code and exploring this possibility. The simplicity of the modified codebase would allow further, deeper research if the results are promising.

Methodology

Codebase and modifications

We use the publicly available AlterBBN code24 as the foundation for our simulations.
AlterBBN is a robust, open-source package that computes the evolution of light element abundances in the early universe by numerically solving the coupled Boltzmann and Friedmann equations under standard cosmological assumptions.
It tracks the neutron-to-proton ratio, nuclear reaction rates, expansion rate, and the entropy content of the plasma to produce the predicted abundances of primordial deuterium, ^4He, D, ^3He, and ^7Li.

To incorporate the physics of evaporating PBHs into BBN, we modified AlterBBN in several significant ways:

We implemented a new module that calculates the time-dependent PBH energy density, \rho_{\mathrm{PBH}}(t), based on the evolving PBH mass and number density. We included this in the total energy density budget driving the Hubble expansion.

  • We added additional terms to the temperature evolution and entropy equations to account for the energy and entropy injection from Hawking radiation. These terms are computed dynamically at each time step, ensuring consistency with the evolving PBH parameters.
  • We compute the entropy source from PBH evaporation and adjust the radiation and baryon components accordingly. While the photon density and total entropy are updated at each step, the baryon-to-photon ratio \eta is assumed to remain close to its initial value and is not explicitly evolved as a separate dynamical variable.
  • We incorporated our own routines for the PBH mass loss rate due to Hawking radiation, explicitly including approximate greybody factors for all relevant species (photons, neutrinos, electrons, quarks, gluons), species-dependent thresholds, and degeneracy factors to compute the total power output and mass loss rate accurately25.
  • We exposed M_{\mathrm{PBH}} and \beta_{\mathrm{PBH}} as user-specified free parameters, configurable at runtime via an input file or command-line arguments, allowing systematic scans over the PBH parameter space.
  • We set up safeguards to prevent the PBH energy density from exceeding 1\% of the radiation density (capped at \rho_{\mathrm{PBH}}/ \rho_{\mathrm{rad}} \leq 10^{-2}) to maintain radiation domination, as required by cosmological observations11.

To validate our implementation, we ran the modified code with \beta_{\mathrm{PBH}} \to 0, confirming that it reproduced standard BBN predictions, consistent with published AlterBBN results24.

We also developed diagnostic outputs to aid debugging and verification, including real-time console logs of PBH mass evolution, energy density contributions, entropy source terms, and the simulation time and temperature at each integration step.
These outputs allow us to monitor the PBH effects on the cosmic plasma throughout the calculation.

At a high level, our modifications to AlterBBN consist of incorporating the effects of PBH evaporation into the cosmic energy density, entropy evolution, and Hubble expansion rate.
These additions allow us to self-consistently track how PBH energy injection alters the thermal history and nuclear reaction network during Big Bang Nucleosynthesis.

We discuss the detailed physics model implemented in our code in the next subsection.

Theory model

We model the influence of evaporating PBHs on Big Bang Nucleosynthesis (BBN) by calculating their energy injection into the radiation-dominated plasma during the first few minutes of the universe. Primordial black holes (PBHs) are assumed to form from the collapse of horizon-scale overdensities at very early times26,27.

We consider PBHs with a single (monochromatic) mass M_{\mathrm{PBH}} and a fractional abundance \beta_{\mathrm{PBH}}, defined as the ratio of the PBH energy density to the radiation energy density at the time of formation18.

At formation, the cosmic time t_{\mathrm{form}} is approximately

(3)   \begin{equation*}t_{\mathrm{form}} \simeq \gamma \frac{G M_{\mathrm{PBH}}}{c^3}\simeq \gamma \left(2.48 \times 10^{-24} \mathrm{s}\right)\left(\frac{M_{\mathrm{PBH}}}{10^{15} \mathrm{g}}\right),\end{equation*}


where we assume that M_{\mathrm{PBH}} equals the horizon mass at that epoch18.
Using M_0 = 5\times10^8 \mathrm{g}=5\times10^5 \mathrm{kg} and \gamma=0.2,

(4)   \begin{equation*}t_{\mathrm{form}} = \gamma \frac{G M_0}{c^3}\approx 2.48 \times 10^{-31}\,\mathrm{s},\end{equation*}

At this time, the radiation energy density \rho_r(t_{\mathrm{form}} is given by the Friedmann equation:

(5)   \begin{equation*}\rho_{r}(t_{\mathrm{form}}) = \frac{3}{32 \pi G t_{\mathrm{form}}^{2}},\end{equation*}

We define the PBH abundance as:

(6)   \begin{equation*}\beta_{\mathrm{PBH}} =\frac{\rho_{\mathrm{PBH}}(t_{\mathrm{form}})}{\rho_{r}(t_{\mathrm{form}})},\end{equation*}

Since the energy density of PBHs at formation is simply their number density n_{\mathrm{PBH}} multiplied by M_{\mathrm{PBH}}, the initial PBH number density is given by:

(7)   \begin{equation*}n_{\mathrm{PBH}}(t_{\mathrm{form}}) =\frac{\beta_{\mathrm{PBH}} \rho_{r}(t_{\mathrm{form}})}{M_{\mathrm{PBH}}},\end{equation*}

After the formation, PBHs behave as pressureless, non-relativistic matter, and their number density evolves as n_{\text{PBH}}(t) \propto a(t)^{-3}, where the scale factor evolves as a(t) \propto t^{1/2} during the radiation-dominated era28. At the same time, each PBH evaporates via 29 radiation26, with its mass decreasing according to:

(8)   \begin{equation*}\frac{dM}{dt} = - \frac{\hbar c^{4}}{15360 \pi G^{2} M^{2}},\end{equation*}

The evaporation rate is faster for lighter PBHs because the Hawking temperature is inversely proportional to mass26:

(9)   \begin{equation*}T_{\mathrm{PBH}} = \frac{\hbar c^{3}}{8 \pi G k_{B} M}\approx 1.06 \times 10^{13}\,\mathrm{GeV}\left(\frac{1\,\mathrm{g}}{M}\right),\end{equation*}

Accordingly, the lifetime of a PBH can be estimated by integrating the mass-loss rate. Starting from

(10)   \begin{equation*}\begin{aligned}\frac{dM}{dt} &= - K \frac{1}{M^2}, \quad K \equiv \frac{\hbar c^4}{15360 \pi G^2}, \quad \int_{M_0}^{0} M^2 \, dM = - K \int_0^{\tau} dt \\&\Rightarrow \left[ \frac{M^3}{3} \right]_{M_0}^{0} = - K \tau \quad \Rightarrow \quad \tau = \frac{M_0^3}{3 K} = \frac{5120 \pi G^2}{\hbar c^4} M_0^3 \propto M_0^3.\end{aligned}\end{equation*}

This explicitly shows that the PBH lifetime scales as \tau \propto M^3.
Numerically, using SI units, the prefactor evaluates to 5120 \pi G^2/(\hbar c^4) \approx 8.41 \times 10^{-17}\ \mathrm{s\,kg^{-3}}.
For example, a PBH with M_0 = 5 \times 10^8~\mathrm{g} has a lifetime

(11)   \begin{equation*}\begin{aligned}\tau \approx 8.41 \times 10^{-26} \, (5 \times 10^8)^3 \ \mathrm{s} \approx 10.5~\mathrm{s}.\end{aligned}\end{equation*}

Including greybody factors and the full set of Standard Model degrees of freedom modifies the numerical prefactor but preserves the \tau \propto M^3 scaling25.

To account for the fact that blackholes aren’t perfect blackbody radiators due to spacetime curvature near their horizon, we include approximate greybody factors . These modify the spectrum for each particle and depend on the dimensionless energy x = E/kT. Each species has its own greybody suppression factor \Gamma_i(x), degeneracy g_i and threshold. The total Hawking power is computed by integrating over x and summing over all species:

(12)   \begin{equation*}\begin{aligned}P_{\text{Hawking}} =\sum_{i} g_i \int_0^\infty \frac{\Gamma_i(x) \, x^3}{\exp(x) \pm 1} \, dx.\end{aligned}\end{equation*}

In our code, greybody factors \Gamma_i(x) are implemented using simple analytic approximations that match the leading-order spectral behavior. While more precise tabulated values exist, our fits are sufficient to capture the qualitative impact of Hawking radiation on the thermal plasma.

The PBH energy density at any time is given by

(13)   \begin{equation*}\begin{aligned}\rho_{\text{PBH}}(t) = M(t) \cdot n_{\text{PBH}}(t)\end{aligned}\end{equation*}

and it contributes to the total energy density driving the cosmic expansion.
The entropy injection and energy release from PBHs are modeled using the Hawking mass-loss rate.
For a monochromatic PBH population, the energy injection rate at time t is

(14)   \begin{equation*}\begin{aligned}\dot{\rho}{\rm PBH}(t) = -n{\rm PBH}(t)\, \frac{dM}{dt}\, c^2= n_{\rm PBH,0} \left(\frac{a_0}{a(t)}\right)^3 \frac{\alpha\, c^2}{\left[M_0^3 - 3 \alpha t\right]^{2/3}},\end{aligned}\end{equation*}

where \alpha = \hbar c^4/(15360 \pi G^2) includes greybody and species factors.
This captures the correct physical behavior: the injection rate is small at early times and rises sharply near the end of the PBH lifetime.
The entropy source term is then computed as

(15)   \begin{equation*}\begin{aligned}\dot{S}{\rm PBH}(t) = \frac{\dot{\rho}{\rm PBH}(t)}{T(t)},\end{aligned}\end{equation*}

assuming instantaneous thermalization of the injected energy.
Both temperature and time evolve during BBN, and in our implementation, T and t are tracked separately, consistent with the radiation-domination relation T \propto t^{-1/2}.

During its life, PBHs emit high-energy particles, including photons, electrons, neutrinos, and hadrons that interact with the plasma29. These injections modify the expansion rate of the universe, alter the neutron-to-proton ratio before freeze-out, and photodissociate light nuclei after formation, all of which affect the primordial abundance of light elements.

The full nuclear reaction network is solved using AlterBBN code, which computes the evolution of D, ^4He, and ^7Li abundances. The reaction rates are taken from experimental data11 and include weak interaction, fusion reactions, and their inverse processes. These rates carry known uncertainties at the \sim 5–10\% level2.

The QCD equation of state, determining g_*(T) and h_*(T), is included through tabulated values30 and governs the energy and the entropy densities of the plasma:

(16)   \begin{equation*}\begin{aligned}\rho_{\text{rad}}(T) = \frac{\pi^2}{30} g_(T) T^4, \quad s_{\text{rad}}(T) = \frac{2\pi^2}{45} h_(T) T^3.\end{aligned}\end{equation*}

To ensure radiation domination during BBN and maintain numerical stability, we impose the following upper bound on PBH energy density to remain consistent with CMB and light-element observational bounds:

(17)   \begin{equation*}\begin{aligned}\frac{\rho_{\text{PBH}}}{\rho_{\text{rad}}} \leq 10^{-2}\end{aligned}\end{equation*}

Note on \beta_{\text{PBH}} Convention

In this work, we define the initial PBH abundance at the time of formation as

(18)   \begin{equation*}\begin{aligned}\beta_{\text{PBH}} = \left. \frac{\rho_{\text{PBH}}}{\rho_r} \right|{t{\text{form}}},\end{aligned}\end{equation*}

where \rho_{\text{PBH}} is the energy density in primordial black holes and \rho_r is the radiation energy density at the same time. This ratio directly represents the energy fraction in PBHs when they form and is the quantity implemented in our modified AlterBBN code.

However, this definition differs from that used in much of the literature, including31, where the PBH abundance is defined in a rescaled form that more directly reflects the impact of PBHs on later cosmological evolution. Specifically, many studies define \beta as the fraction of the total energy density in PBHs at horizon re-entry, including numerical factors to account for collapse efficiency, entropy dilution, and the difference in redshifting behavior between matter and radiation.

To understand why this leads to a difference of orders of magnitude, we consider how the energy densities evolve. After formation, PBHs behave like non-relativistic matter, scaling as \rho_{\text{PBH}} \propto a^{-3} , while radiation scales as \rho_r \propto a^{-4}. Therefore, their ratio increases with the scale factor:

(19)   \begin{equation*}\begin{aligned}\frac{\rho_{\text{PBH}}(t)}{\rho_r(t)} = \beta_{\text{PBH}} \cdot \frac{a(t)}{a(t_{\text{form}})}.\end{aligned}\end{equation*}

During radiation domination, the product ( aT ) remains approximately constant, so

(20)   \begin{equation*}\begin{aligned}\frac{a(t)}{a(t_{\text{form}})} = \frac{T_{\text{form}}}{T(t)}.\end{aligned}\end{equation*}

If the PBHs form at T_{\text{form}} \sim 10^9~\text{GeV} and evaporate around T_{\text{evap}} \sim 1~\text{MeV} = 10^{-3}~\text{GeV}, then the ratio of scale factors becomes

(21)   \begin{equation*}\begin{aligned}\frac{T_{\text{form}}}{T_{\text{evap}}} \sim 10^{12}.\end{aligned}\end{equation*}

Thus, even a tiny initial abundance \beta_{\text{PBH}} \sim 10^{-25} leads to a much larger ratio \rho_{\text{PBH}}/ \rho_r \sim 10^{-13} at the time of evaporation. This time-dependent amplification is implicitly absorbed into the literature’s definition of \beta, which often represents the effective impact of PBHs near their evaporation epoch rather than their initial fraction.

Therefore, to compare our results with values from the literature, we can roughly estimate that

(22)   \begin{equation*}\begin{aligned}\beta_{\text{literature}} \sim \left( \frac{T_{\text{form}}}{T_{\text{evap}}} \right) \cdot \beta_{\text{PBH}} \sim 10^{12} \cdot \beta_{\text{PBH}},\end{aligned}\end{equation*}

depending on the specific formation and evaporation temperatures. For example, a value \beta_{\text{PBH}} = 10^{-25} in our convention corresponds to approximately \beta_{\text{literature}} \sim 10^{-13}, which is a typical value cited in PBH constraint plots.

Throughout this paper, we consistently use our definition of \beta_{\text{PBH}} as the initial energy fraction at formation, which directly corresponds to the parameter implemented in the numerical simulation.

Parameters and Setup

We explored the impact of PBHs on BBN by scanning a grid of physically motivated parameter values. The two key free parameters in our study are the PBH mass at formation, M_{\text{PBH}} , and their fractional abundance at formation, \beta_{\text{PBH}}.

Motivated by previous works20,23,21, we fixed the PBH mass to M_{\text{PBH}} = 5 \times 10^8\,\mathrm{g} , corresponding to PBHs that evaporate during the BBN epoch. This mass was chosen because PBHs in this range inject energy at timescales relevant for altering light-element abundances without completely dominating the expansion or surviving to late times as dark matter. Lighter PBHs evaporate too early to have any significant effects, while heavier PBHs survive past recombination.

For this fixed mass, we varied the initial abundance \beta_{\text{PBH}} over the range 10^{-18} to 10^{-12}31 (convention) , subject to the imposed upper bound \rho_{\text{PBH}}/\rho_{\text{rad}} \leq 10^{-2} at all times11. The lower bound allows us to probe even extremely small PBH contributions to BBN in detail.

ParameterValue
η6.1 x 10-10
MPBH5 x 108  g
ΒPBH range10-18 – 10-12
ΡPBH/ρr   cap10-2
Table 1 | Simulation parameters used throughout this work

For each choice of \beta_{\text{PBH}}, we ran the modified code and computed the resulting primordial abundances of the deuterium-to-hydrogen ratio (D/H), helium-4 mass fraction Y_p, and lithium-7-to-hydrogen ratio ^7\mathrm{Li}/\mathrm{H}. These outputs were compared to observational constraints to identify the parameter regions where PBH evaporation could help alleviate the lithium-7 discrepancy without violating constraints from deuterium and helium-4 abundances2,32,1. We adopted a baryon-to-photon ratio \eta_0 = 6.1 \times 10^{-10}, consistent with Planck CMB measurements12. Standard nuclear reaction rates and cross sections, as implemented in AlterBBN24, were used throughout.

All simulations were performed on a Linux-based desktop workstation using a GCC-compiled version of the code. Computations typically completed in seconds to minutes per run, allowing efficient exploration of \beta_{\text{PBH}} over many orders of magnitude. Results were post-processed and visualized using Python scripts, producing plots of predicted abundances as a function of \beta_{\text{PBH}}.

To ensure reproducibility, we maintained a version-controlled repository of our modified code, input configurations, and output data.

Observational Abundance Constraints

To determine whether a given PBH parameter point is consistent with astrophysical observations, we impose the current best–fit primordial abundances for deuterium, helium, and lithium. The following 1\sigma observational ranges are adopted throughout this work:
Deuterium:

(23)   \begin{equation*}{\mathrm{D/H}}_{\mathrm{obs}} = 2.547 \pm 0.025 \times 10^{-5},\end{equation*}


based on high-redshift damped Lyman-\alpha absorber measurements.

Helium-4 mass fraction:

(24)   \begin{equation*}Y_{p} = 0.2449 \pm 0.0040,\end{equation*}

inferred from metal-poor extragalactic H, {ii} regions.

Lithium-7:

(25)   \begin{equation*}{}^{7}\mathrm{Li}/\mathrm{H})_{\mathrm{obs}} = (1.6 \pm 0.3)\times 10^{-10},\end{equation*}

determined from Population~II halo stars.

These observational values are used as consistency criteria when scanning the PBH mass–abundance parameter space. Parameter points that produce primordial abundances outside these ranges are marked as observationally excluded in the Results section.
In our analysis, a parameter point is classified as observationally allowed only if all three predicted abundances (D/H, Y_p, and ^7Li/H) lie within these 1\sigma ranges simultaneously. These observational values follow the compilation in Fields et al. (2020)33.

Validation Against PArthENoPE

To validate our PBH evaporation implementation, we compare the AlterBBN predictions with PArthENoPE for the same PBH mass M_{\rm PBH} = 5\times10^{8}\,g over several orders of magnitude in PBH abundance. Because the two codes use different conventions for \beta_{\rm PBH} (formation-time in AlterBBN versus horizon-entry in PArthENoPE), we compare elemental
abundances directly rather than the numerical values of \beta_{\rm PBH}.

Helium-4 mass fraction Y_p

βPBH (PArth.)YP (PArth. )βPBH(ABBN)YP(ABBN)Percent diff.
5 x 10-190.2478 x 10-290.2471+.04
5 x 10-180.2488 x 10-280.2471-0.36
5 x 10-170.2498 x 10-270.2475-0.60
5 x 10-160.25258 x 10-260.2520-0.20
5 x 10-150.2608 x 10-250.3365+29.4
Table 2 | Comparison of Y_p between PArthENoPE and our modified AlterBBN for M_{\rm PBH}=5\times10^{8} g. Percent differences are (Y_p^{\rm ABBN} - Y_p^{\rm PArth}) / Y_p^{\rm PArth} \times 100\%.

Deuterium abundance D/H

βPBH (PArth.)D/H(PArth. )βPBH(ABBN)D/H(ABBN)Percent diff.
5 x 10-192.508 x 10-292.456-1.76
5 x 10-182.488 x 10-282.455-1.01
5 x 10-172.308 x 10-272.432+5.73
5 x 10-161.008 x 10-262.218+121.8
5 x 10-15~08 x 10-250.2485n.a.
Table 3 | Comparison of D/H (in units of 10^{-5}) between PArthENoPE and AlterBBN for M_{\rm PBH}=5\times10^{8} g.

Across the first three decades in PBH abundance, AlterBBN reproduces PArthENoPE predictions at the \lesssim 1\% level for both Y_p and D/H. Larger discrepancies at high \beta_{\rm PBH} arise from extreme energy injection, which lies in a region already observationally excluded.

Standard BBN Validation

Before introducing PBHs, we verify that our modified AlterBBN code reproduces standard BBN results when \beta_{\rm PBH}=0. Table 1 compares our abundances to the reference PArthENoPE predictions1. The agreement at the sub-percent level for Y_p and D/H, and at the few-percent level for ^7Li/H, confirms that the PBH extensions do not alter the baseline BBN evolution.

AbundanceThis workReferencePercent diff
YP0.24710.24709+0.004
D/H (x10-5)2.4562.45+0.26
7Li/H (x10-10)5.3825.62-4.24
Table 4 | Standard BBN validation with \beta_{\rm PBH}=0. “Reference” values are taken from PArthENoPE (Fields et al.~2020). Percent differences are computed as (X_{\rm this} - X_{\rm ref})/X_{\rm ref} \times 100\%.

Comparison with Observational Constraints

To assess the viability of PBH-induced energy injection during BBN, we compare our simulated elemental abundances against current observational limits for light elements.

Our results indicate that for \beta_{\text{PBH}} \gtrsim 10^{-15}, all three key abundances, \mathrm{^7Li}, D, and Y_p begin to deviate significantly from observationally inferred values:

  • \mathrm{^7Li} rises well above the already overpredicted SBBN level.
  • Deuterium falls below its lower observational bound.
  • Helium-4 rises above its upper observational limit, with Y_p \gtrsim 0.27 at high \beta_{\text{PBH}}.

At \beta_{\text{PBH} \sim 10^{-14}, all the abundances are in strong conflict with observational data, clearly ruling out such a high PBH abundance for this mass scale. These findings suggest that, within the M_{\text{PBH}} = 5 \times 10^8~\mathrm{g} range, PBH evaporation does not resolve the lithium problem, in fact it worsens it – while simultaneously violating deuterium and helium-4 abundance constraints.

To validate the credibility of our bounds, we compare our results with earlier literature, particularly the study by Boccia et al.31, which analyzed the BBN effects of PBH evaporation and reported an upper limit of \beta_{\text{PBH}} \lesssim 10^{-16} for M_{\text{PBH}} = 5 \times 10^8~\mathrm{g}.

When translated to our convention where \beta_{\text{PBH}} is defined as the initial PBH energy fraction relative to the total energy density at formation, without horizon scaling31 limit corresponds to \beta_{\text{PBH}} \lesssim 10^{-27}, which closely aligns with our derived upper bound.

This agreement reinforces the accuracy of our simulation and supports the claim that our modified AlterBBN code captures the essential physics of PBH evaporation and its effects on light element synthesis. Minor discrepancies within an order of magnitude are expected, given differences in nuclear reaction rates, greybody factor approximations, and thermalization assumptions.

D/H from our modified
AlterBBN code
Yp from our modified AlterBBN
code
D/H from Boccia et al. (using
Parthenope)
Yp from Boccia et al. (using
Parthenope)
Figure 1 | Comparison of deuterium (D/H, left) and helium-4 mass fraction (Yp, right) between our modified AlterBBN results (top row) and34 obtained with the Parthenope code (bottom row). The qualitative agreement validates our implementation.

As shown in Fig.1, our simulated abundance trends for D/H, Y_p, and ^7Li/H closely match those obtained by Boccia et al.31, once our \beta_{\mathrm{PBH}}^{\mathrm{Boccia}} values are converted to their \beta' convention via Eq. 18.

βPBHYPD/H (x 10-5)7Li/H (x 10-10)
8 x 10-170.24712.4565.382
8 x 10-160.24712.4565.382
8 x 10-150.24712.4555.388
8 x 10-140.24752.4325.464
8 x 10-130.25202.2186.260
8 x 10-120.33650.248542.46
Table 5 | AlterBBN abundances across the PBH abundance scan for M_{\rm PBH}=5\times10^{8} g.

The qualitative behaviour of each element is consistent across both works, validating our numerical setup. Small quantitative differences at high \beta_{\mathrm{PBH}} can be traced to differences in the nuclear reaction network implementation, choice of thermal history treatment, and interpolation of cross-section data. This agreement provides an important cross-check and builds confidence in the robustness of our results.

Our simulations explored the impact of PBH evaporation on the primordial abundances of light elements by injecting energy into the plasma during BBN. Here we discuss our results in detail, interpret the observed trends, compare them with previous studies, and highlight the implications, limitations, and prospects for future work.

Results and Discussion

In this section, we present the results of our simulation investigating the effect of evaporating PBHs on BBN. We have fixed the PBH mass at M_{\mathrm{PBH}} = 5 \times 10^8~\mathrm{g} and varied the initial abundance \beta_{\mathrm{PBH}} over the range 10^{-18} - 10^{-12}. The abundances of key light elements  lithium-7 ^7Li/H, deuterium D/H, and the ^4He mass fraction Y_p , were computed as functions of \beta_{\mathrm{PBH}}.

Allowed Parameter Range

From our simulations and comparisons with observational constraints on light element abundances, we find that values of \beta_{\mathrm{PBH}} \gtrsim 10^{-15} (for M_{\mathrm{PBH}} = 5 \times 10^8~\mathrm{g}) begin to produce significant deviations from the standard BBN (SBBN) predictions. In particular, we observe that for \beta_{\mathrm{PBH}} \gtrsim 10^{-14}, the helium-4 mass fraction Y_p and the lithium-7 abundance (^7Li/H) exceed their respective observational upper limits, while the deuterium abundance \mathrm{D}/\mathrm{H} drops below its lower bound.

These deviations become increasingly pronounced with rising \beta_{\mathrm{PBH}}, indicating that excessive energy injection from PBH evaporation disrupts the thermal and nuclear conditions required for accurate light-element synthesis.

Based on these results, we constrain the initial PBH abundance to \beta_{\mathrm{PBH}} \lesssim 10^{-14} for M_{\mathrm{PBH}} = 5 \times 10^8~\mathrm{g}. This upper bound is consistent with earlier constraints reported in the literature31, after accounting for differences in normalization conventions. Our findings support the conclusion that PBHs on this mass scale must have been extremely rare to avoid spoiling the predictions of SBBN.

Effects on Light Element abundances

(i) Origin of the Lithium-7 Enhancement

Lithium-7 is particularly sensitive to modifications of the thermal and nuclear reaction history during Big Bang Nucleosynthesis (BBN). Our results show that for \beta_{\rm PBH} \lesssim 10^{-15}, the predicted ^7Li/H abundance remains close to the Standard BBN (SBBN) value. However, as \beta_{\rm PBH} increases beyond this threshold, the lithium abundance rises sharply.

For \beta_{\rm PBH} \gtrsim 10^{-15}, the ^7Li/H ratio increases by more than a factor of three compared to the SBBN prediction. This enhancement is driven by two related effects. First, PBH evaporation alters the nuclear reaction network in a way that enhances the production of ^7Be, which subsequently decays into ^7Li via electron capture. Second, energy injection from PBH evaporation modifies reaction rates and effectively prolongs the time available for lithium synthesis.

Figure~2 illustrates this behavior. Panel (A) shows the continuous dependence of ^7Li/H on \beta_{\rm PBH} obtained from our modified AlterBBN code, while panel (B) presents discrete data points with associated numerical uncertainties. In both cases, PBH evaporation leads to lithium production well above the SBBN value, thereby worsening the lithium problem rather than alleviating it.

The physical origin of this behavior can be understood using a simple order-of-magnitude estimate. Although PBHs constitute only a tiny fraction of the total energy density during BBN (for example, \rho_{\rm PBH}/\rho_{\rm rad} \sim 10^{-6} for \beta_{\rm PBH} = 8 \times 10^{-26}), the energy injected per baryon is substantial. At temperatures T \sim 1\,\mathrm{MeV}, the photon number density is n_\gamma \simeq 3.2 \times 10^{37}\,\mathrm{m}^{-3}, implying a baryon number density n_b = \eta n_\gamma \simeq 2 \times 10^{28}\,\mathrm{m}^{-3}.

The PBH energy density measured in our simulations, \rho_{\rm PBH} \simeq 1.5 \times 10^{3}\,\mathrm{kg\,m^{-3}}, corresponds to an injected energy per baryon of

(26)   \begin{equation*} E_{\rm per\,baryon} \simeq \frac{\rho_{\rm PBH} c^2}{n_b} \sim 4 \times 10^{10}\,\mathrm{eV}, \end{equation*}

i.e. tens of GeV per baryon. Such energy injection is sufficient to generate energetic non-thermal photons and hadrons that significantly affect the nuclear reaction network.

These non-thermal particles efficiently destroy deuterium, which has a low binding energy, but do not efficiently destroy ^7Be or ^7Li due to their higher photodisintegration thresholds and smaller cross-sections. At the same time, the production of extra neutrons and ^3He in the non-thermal cascades enhances the reaction

^3\mathrm{He}(\alpha,\gamma)^7\mathrm{Be},

thereby increasing the ^7Be abundance. Since ^7Be subsequently decays into ^7Li, the net effect of PBH evaporation in this mass and abundance range is an increase in the final lithium-7 abundance.

Predicted 7Li/Hvs. βPBH (contin-
uous curve) from our modified
AlterBBN code
Predicted 7Li/H vs. βPBH (discrete
points with error bars) from our
modified AlterBBN code.
FIGURE 2 | Panel (A) shows the continuous trend from our modified AlterBBN code, while panel (B) includes
discrete data points with error bars for comparison. PBH evaporation enhances lithium production beyond the
standard BBN value, worsening the lithium problem.

(ii) Evolution of D/H

In contrast to lithium-7, the deuterium abundance decreases with increasing PBH abundance. For small values of \beta_{\text{PBH}}, the deuterium-to-hydrogen ratio (D/H) remains nearly constant and agrees with the SBBN prediction. However, once \beta_{\text{PBH}} \gtrsim 10^{-15}, the deuterium abundance begins to gradually decline.

This depletion becomes rapid for \beta_{\text{PBH}} \gtrsim 10^{-14}. The primary cause is photodissociation: high-energy photons emitted during PBH evaporation can efficiently destroy deuterium nuclei, which have a relatively low binding energy of \sim 2.2~\mathrm{MeV}.

Predicted (D/H) vs. βPBH (continuous curve) from our modified AlterBBN code
Predicted (D/H) vs. βPBH (discrete points with error bars) from our modified AlterBBN code.
FIGURE 3 | Predicted deuterium abundance (D/H) as a function of PBH abundance βPBH, computed with our
modified AlterBBN code. Panel (A) shows the continuous curve, while panel (B) displays the same data with
discrete points and error bars. At large PBH abundances, photodissociation processes drive D below its observationally
allowed range.

(iii) Evolution of Y_p

The helium-4 mass fraction Y_p shows a steady increase with rising \beta_{\text{PBH}} For \beta_{\text{PBH}} \lesssim 10^{-15}, Y_p \approx 0.25, consistent with the SBBN prediction. As \beta_{\text{PBH}} increases, Y_p begins to rise due to two effects: delayed neutron-proton freeze-out and partial reheating of the plasma.

At \beta_{\text{PBH}} \sim 10^{-13}, Y_p exceeds 0.29, significantly above observational constraints. This enhancement arises from the elevated neutron-to-proton ratio at freeze-out, since nearly all remaining neutrons eventually form helium-4.

Predicted Y p vs. βPBH (continuous curve) from our
modified AlterBBN code.
Predicted Yp vs. βPBH (discrete points with error bars)
from our modified AlterBBN code.
FIGURE 4 | Predicted helium-4 mass fraction (Yp) as a function of PBH abundance βPBH, computed with our modified AlterBBN code. Panel (A) shows the continuous curve, while panel (B) displays the same data with discrete points and error bars. Energy injection delays neutron–proton freeze-out, raising Yp above observational
limits for βPBH ≳ 10−15.

Global Parameter Scan: Joint Dependence on \{beta\_{PBH}} and M\_{PBH}

To explore the broader PBH parameter space, we perform a two-dimensional scan over the initial PBH abundance \beta_{\rm PBH} and the PBH mass M_{\rm PBH}, and visualize the resulting primordial abundances using the three-dimensional surface plots shown in Fig.~5. This scan reveals a narrow region in the (M_{\rm PBH}, \beta_{\rm PBH}) plane where the predicted ^7\mathrm{Li}/\mathrm{H} abundance is mildly suppressed relative to the standard Big Bang Nucleosynthesis (BBN) value. However, even within this region, the lithium abundance remains significantly above the observationally inferred value and is therefore phenomenologically irrelevant. The detailed physical origin of this suppression is not fully clear within our simplified treatment; it likely reflects a subtle interplay between the timing of PBH evaporation, neutron injection, and the onset of ^7\mathrm{Be} synthesis. A definitive explanation would require a more detailed non-equilibrium treatment of hadronic and electromagnetic cascades, which is beyond the scope of this work.

Outside this narrow band, PBH evaporation either occurs too early or too late to meaningfully influence nucleosynthesis, or the associated energy injection produces substantial distortions in deuterium and ^4He abundances, which are tightly constrained observationally. The joint scan therefore demonstrates that, although PBH evaporation can transiently suppress lithium production within a finely tuned region of parameter space, this effect is insufficient to resolve the primordial lithium problem.

FIGURE 5 | Global parameter scan of primordial abundances in the (βPBH,MPBH) plane, generated with our modified
AlterBBN code. A narrow valley appears where 7Li/H is suppressed, but abundances remain inconsistent with observations, showing that PBH evaporation cannot resolve the lithium problem.

Physical Interpretation of the trends in abundance

At small values of \beta_{\text{PBH}} \lesssim 10^{-18}, the impact of PBH evaporation is negligible. The abundances of light elements remain consistent with Standard BBN (SBBN) predictions, as expected. In this regime, the energy density of PBHs is so small compared to the radiation background that their contribution to the plasma’s thermal or dynamical evolution is insignificant35.

As \beta_{\text{PBH}} increases beyond \sim 10^{-15}, the injected energy from PBH evaporation begins to affect key processes in nucleosynthesis. We observe the following deviations from standard predictions:

  • Helium-4 Mass Fraction Y_p Increases:

The early injection of energy, particularly from high-energy particles delays the neutron-to-proton freeze-out by maintaining thermal equilibrium slightly longer. As a result, more neutrons survive until nucleosynthesis begins, leading to a higher synthesis of helium-4, since nearly all available neutrons end up in ^4\text{He} nuclei5.

  • Lithium-7 Abundance ^7\text{Li}/ \text{H} Increases:

Contrary to some expectations, we find that lithium-7 abundances increase rather than decrease. This is primarily because non-thermal hadrons and photons emitted during PBH evaporation modify the nuclear reaction flow. In particular, reactions producing ^7\text{Be} (which later decays to ^7\text{Li}) become more efficient, and destruction channels are insufficient to counteract this enhancement36. Thus, PBHs in the considered mass range aggravate rather than alleviate the lithium problem.

  • Deuterium Abundance \text{D}/ \text{H} Decreases:

Deuterium is the most weakly bound of the light nuclei and is especially sensitive to photodissociation. The injection of high-energy photons during or shortly after nucleosynthesis dissociates existing deuterium nuclei, reducing the final \text{D}/ \text{H} ratio. This is a key constraint on PBH models, as current observations of deuterium are highly precise.

These trends are physically intuitive when viewed through the lens of thermal and non-thermal processes during BBN. PBH evaporation injects entropy, disturbs the baryon-to-photon ratio \eta(t), and alters the expansion rate via increased energy density. Additionally, energetic secondary particles open new or accelerated reaction channels.

These results collectively indicate that for the PBH mass and abundance range considered, PBH evaporation worsens the lithium-7 discrepancy instead of resolving it. While such models provide a rich framework to study non-standard cosmology, the specific parameter space we explored suggests PBHs are unlikely to be the primary solution to the lithium problem without introducing additional mechanisms or constraints.

Physical intuition and implications

The results reinforce the idea that PBHs with mass M_{\text{PBH}} = 5 \times 10^8~\mathrm{g} inject energy too late, after weak interaction freeze-out (which determines the neutron-to-proton ratio) but during the nuclear burning phase when light nuclei form. This timing is crucial: the injected energy heats the plasma and alters reaction rates in a way that enhances the production of ^7\mathrm{Be}, which later decays into ^7\mathrm{Li}, thereby worsening the lithium problem rather than alleviating it.

For PBH evaporation to reduce the lithium-7 abundance, the injected energy would need to arrive earlier, ideally before or during neutron-proton freeze-out, to shift the neutron-to-proton ratio or suppress the formation of beryllium-7. However, such early evaporation would require lighter PBHs (e.g., M_{\text{PBH}} \lesssim 10^8~\mathrm{g}), which evaporate too rapidly to inject significant energy during the nucleosynthesis window18. This presents a fine-tuning challenge: PBHs must be massive enough to influence BBN but light enough to evaporate early enough to affect the relevant nuclear processes.

Our constraint \beta_{\text{PBH}} \lesssim 10^{-14} implies that PBHs must remain strongly subdominant to the radiation density during BBN to avoid disrupting the observed abundances of light elements. This reinforces earlier findings22 and places tight limits on PBH populations in this mass range. In summary, PBHs in the \sim10^8~\mathrm{g} range are unlikely to resolve the lithium problem, and any viable PBH-based mechanism would require more careful tuning of mass and abundance to achieve the desired effect without violating other observational constraints.

Limitations and future work

While our study provides a detailed and quantitative exploration of the impact of PBHs on BBN, it is subject to several simplifications, assumptions, and modeling limitations that could be addressed in future works. While the refinements below could shift the exact numerical upper bound on \beta_{\rm PBH}, for example, bringing our value even closer to34, they are not expected to produce any qualitative change.

  • Monochromatic PBH Mass Spectrum

We assume a monochromatic mass spectrum for PBHs where all black holes have the same mass M_{\text{PBH}} = 5 \times 10^8\,\mathrm{g}. While this choice allows for clear isolation of effects near critical evaporation timescales, it neglects the possible influence of broader or extended mass functions, which are generally expected in realistic early-universe collapse scenarios. Such spectra could lead to energy injection over a wider time window, with cumulative effects on nuclear reaction rates, freeze-out, and photodisintegration processes. Future work should consider log-normal, power-law, or critical collapse-inspired distributions to assess whether such effects could shift or even alleviate the lithium problem6.

  • Instantaneous Thermalization

Our modeling assumes that all energy injected by evaporating PBHs is instantly and uniformly thermalized into the background plasma. In reality, high-energy photons initiate electromagnetic cascades, which may delay full thermalization and change the dissociation rates of nuclei. These particles can initiate cascades or delay energy deposition, leading to time-dependent modifications in the photon spectrum and altered nucleosynthesis dynamics. Incorporating energy deposition kernels and transport modeling could improve the accuracy of the entropy and temperature evolution.

  •  Approximate Greybody Factor

While we included greybody factors in our Hawking radiation calculations, these were implemented via analytical approximations for tractability. The full numerical computation of greybody factors, especially for spin-1 and spin-2 particles, exhibits energy- and species-dependent deviations from blackbody spectra. Using tabulated greybody factors or solving the wave equation numerically for each species would allow more accurate modeling of the energy distribution and evaporation rate25.

  •  Neglect of Hadronic Showers

Our code tracks the effect of total energy injection on the BBN plasma but does not separately model hadronic showers, which arise when emitted quarks and gluons hadronize and interact with light nuclei. These non-thermal processes can be especially important for altering neutron-to-proton ratios and dissociating helium or lithium nuclei. A full treatment of hadronic cascades using transfer functions or Monte Carlo shower simulations would be necessary to rigorously quantify these effects22.

Future extensions of this work can focus on relaxing current assumptions and improving the PBH–BBN interaction model. Key directions include:

  • Incorporating Extended Mass Spectra

Studying a range of PBH masses with varying abundances could reveal whether certain energy injection patterns yield more favorable outcomes for the lithium-7 problem while remaining within the observational bounds for deuterium and helium. Exploring log-normal or power-law spectra may lead to qualitatively different nucleosynthesis effects due to the spread in evaporation times.

  • Lower Mass PBHs

PBHs with masses M_{\text{PBH}} < 10^8\,\mathrm{g} evaporate earlier in the thermal history, potentially before deuterium and lithium formation. Their early evaporation could affect the neutron-to-proton ratio or the baryon-to-photon ratio in ways not captured by heavier PBH models. Testing whether earlier evaporation improves lithium predictions without violating constraints on light element yields or entropy density is a promising direction.

  • Detailed Non-Equilibrium Modeling

Implementing time and energy-dependent deposition models for high-energy particles, such as photon cascades, hadronic showers, or delayed energy injection functions would provide a more realistic estimate of how PBH energy interacts with the thermal plasma and light nuclei. This would help refine our understanding of non-instantaneous thermalization effects23.

  • Dynamic Evolution of \eta(t)

Our analysis assumes a constant baryon-to-photon ratio during BBN, following standard BBN calculations. In principle, entropy injection from PBH evaporation can induce a time-dependent evolution of \eta, since baryon number is conserved while the photon number density may increase due to energy injection into the radiation bath. A fully self-consistent treatment would therefore require solving for the evolution of \eta(t) in the presence of entropy sources.

Schematically, the baryon-to-photon ratio can be written as

(27)   \begin{equation*}\eta(t) = \frac{n_b}{n_\gamma(t)}\end{equation*}

where n_b remains fixed while n_\gamma may increase due to PBH-induced entropy injection. A decrease in \eta modifies the balance between nuclear reaction rates and the expansion of the Universe, primarily affecting deuterium synthesis and, indirectly, the downstream production of heavier nuclei.

In the parameter space explored in this work, however, PBHs remain a subdominant component of the total energy density during BBN, and the associated entropy injection is correspondingly small. As a result, the evolution of \eta(t) is negligible for \beta_{\rm PBH} \lesssim 10^{-15}, and the baryon-to-photon ratio can be treated as effectively constant. Even at larger PBH abundances, variations in \eta constitute a subleading correction and are not the dominant driver of the lithium enhancement observed in our simulations.

We therefore expect that including a self-consistent evolution of \eta(t) would not qualitatively modify the behavior of the ^7Li abundance as a function of \beta_{\rm PBH}, nor alter the main conclusions of this study. Nonetheless, incorporating entropy-induced \eta evolution would be a valuable extension for future work aimed at higher PBH abundances.

  • Fixed Neutrino-Photon Temperature Ratio

In our modelling, the neutrino temperature T_\nu(t) is assumed to evolve in the standard way, maintaining a fixed ratio T_\nu/T_\gamma = (4/11)^{1/3} after e^\pm annihilation. This neglects the slight differential reheating that would occur if PBH evaporation injected entropy primarily into the photon-electron plasma. For the parameter space explored here (M_{\rm PBH} = 5\times 10^8\,\mathrm{g}, \beta_{\rm PBH} \lesssim 10^{-28}), the resulting change in T_\nu/T_\gamma and the corresponding shift in light-element abundances are expected to be at the sub-percent level, well below our current numerical uncertainties. However, for larger \beta_{\rm PBH} or different PBH masses with higher energy injection, evolving T_\nu(t) self-consistently could become important.

Overall, our simplified model already shows that PBHs at M_{\text{PBH}} = 5 \times 10^8~\mathrm{g} do not improve the lithium problem without violating other constraints. Hence, we chose to keep the model minimal and transparent rather than pursuing a more complex implementation, as the initial results were not promising enough to justify further complication. We intentionally kept the model minimal, since our results already show that PBH evaporation in this mass range worsens the lithium problem, and further complication would not qualitatively change this conclusion

Conclusion

In this work, we investigated whether energy injection from evaporating primordial black holes (PBHs) with mass M_{\rm PBH} = 5 \times 10^{8}~\mathrm{g} could alleviate the long-standing lithium-7 discrepancy in Big Bang Nucleosynthesis (BBN). By modifying the AlterBBN code to consistently include PBH evaporation effects, we computed the primordial abundances of deuterium, helium-4, and lithium-7 over a wide range of PBH initial abundances.

Our results demonstrate that PBH evaporation in this mass range does not resolve the lithium problem. Instead, increasing \beta_{\rm PBH} leads to an enhancement of the final ^7Li/H abundance, exacerbating the discrepancy with observations. At the same time, deuterium becomes underproduced and helium-4 overproduced for sufficiently large PBH abundances, yielding stringent constraints on the allowed PBH contribution during BBN. For M_{\rm PBH} = 5 \times 10^{8}~\mathrm{g}, we derive an approximate upper bound \beta_{\rm PBH} \lesssim 8 \times 10^{-14}, consistent with existing limits in the literature after accounting for differing \beta_{\rm PBH} conventions.

This PBH mass was chosen deliberately, as it corresponds to evaporation occurring near the peak of BBN, when modifications to nuclear reaction rates and freeze-out dynamics would be expected to have the greatest potential impact on lithium production. The fact that PBHs in this optimally timed mass window worsen rather than alleviate the lithium problem therefore provides a robust negative result, effectively closing off this region of PBH parameter space as a viable solution.

While our analysis focuses on a monochromatic PBH mass spectrum and assumes standard thermalization of injected particles, these assumptions are unlikely to alter the qualitative outcome. Even within a narrow region of parameter space where lithium is mildly suppressed, the reduction is insufficient to reconcile theory with observations and requires fine-tuning that renders such scenarios implausible. More elaborate PBH models, such as extended mass functions or detailed hadronic shower treatments, may modify quantitative bounds but would still need to overcome the fundamental tendency of PBH evaporation to enhance lithium production.

In conclusion, our findings indicate that evaporating PBHs during BBN do not offer a viable resolution to the primordial lithium problem under standard assumptions, and instead impose tighter constraints on early-Universe black hole scenarios. This result helps narrow the range of proposed solutions and motivates further investigation into alternative explanations, potentially involving new nuclear physics, stellar depletion mechanisms, or physics beyond the Standard Model.

Code Availability

The modified AlterBBN code and representative input files used in this work are available via a public archive at: https://drive.google.com/file/d/1JgYkSK4b6Eo1Un7eOV_IhE5eSxenBze_/view?usp=drive_link.

Acknowledgements

We thank Isaiah James Richardson for guidance and feedback, and the developers of AlterBBN for making their code publicly available.

References

  1. B. D. Fields, K. A. Olive, T.-H. Yeh, C. Young. Big-bang nucleosynthesis after planck. Journal of Cosmology and Astroparticle Physics. Vol. 2020, pg. 010–010, 2020 https://doi.org/10.1088/1475-7516/2020/03/010 [↩] [↩] [↩] [↩]
  2. R. H. Cyburt, B. D. Fields, K. A. Olive, T.-H. Yeh. Big bang nucleosynthesis: present status. Reviews of Modern Physics. Vol. 88, pg. 015004, 2016 https://doi.org/10.1103/RevModPhys.88.015004 [↩] [↩] [↩] [↩]
  3. A. Coc, J.-P. Uzan, E. Vangioni. Standard big bang nucleosynthesis and primordial cno abundances after planck. Journal of Cosmology and Astroparticle Physics. Vol. 2014, pg. 050–050, 2014 https://doi.org/10.1088/1475-7516/2014/10/050 [↩]
  4. O. Richard, G. Michaud, J. Richer. Implications of wmap observations on li abundance and stellar evolution models. The Astrophysical Journal. Vol. 619, pg. 538–548, 2005 https://doi.org/10.1086/426470 [↩]
  5. K. Jedamzik. Did something decay, evaporate, or annihilate during big bang nucleosynthesis? Physical Review D. Vol. 70, pg. 063524, 2004 https://doi.org/10.1103/PhysRevD.70.063524 [↩] [↩] [↩] [↩]
  6. B. J. Carr, K. Kohri, Y. Sendouda, J. Yokoyama. New cosmological constraints on primordial black holes. Physical Review D. Vol. 81, pg. 104019, 2010 https://doi.org/10.1103/PhysRevD.81.104019 [↩] [↩] [↩] [↩] [↩]
  7. A. Arbey. AlterBBN: a program for calculating the bbn abundances of the elements in alternative cosmologies. Computer Physics Communications. Vol. 183, pg. 1822–1831, 2012 https://doi.org/10.1016/j.cpc.2012.03.018 [↩] [↩]
  8. R. H. Cyburt, B. D. Fields, K. A. Olive, T.-H. Yeh. Big bang nucleosynthesis: present status. Reviews of Modern Physics. Vol. 88, pg. 015004, 2016 https://doi.org/10.1103/RevModPhys.88.015004 [↩]
  9. K. A. Olive. Review of particle physics. Chinese Physics C. Vol. 38, pg. 090001, 2014 https://doi.org/10.1088/1674-1137/38/9/090001 [↩] [↩]
  10. D. N. Schramm, M. S. Turner. Big-bang nucleosynthesis enters the precision era. Reviews of Modern Physics. Vol. 70, pg. 303–318, 1998 https://doi.org/10.1103/RevModPhys.70.303 [↩]
  11. P. Collaboration, N. Aghanim, Y. Akrami, M. Ashdown, J. Aumont, C. Baccigalupi, M. Ballardini, A. J. Banday, R. B. Barreiro, N. Bartolo, S. Basak, R. Battye, K. Benabed, J.-P. Bernard, M. Bersanelli, P. Bielewicz, J. J. Bock, J. R. Bond, J. Borrill, F. R. Bouchet, F. Boulanger, M. Bucher, C. Burigana, R. C. Butler, E. Calabrese, J.-F. Cardoso, J. Carron, A. Challinor, H. C. Chiang, J. Chluba, L. P. L. Colombo, C. Combet, D. Contreras, B. P. Crill, F. Cuttaia, P. de Bernardis, G. de Zotti, J. Delabrouille, J.-M. Delouis, E. D. Valentino, J. M. Diego, O. Doré, M. Douspis, A. Ducout, X. Dupac, S. Dusini, G. Efstathiou, F. Elsner, T. A. Enßlin, H. K. Eriksen, Y. Fantaye, M. Farhang, J. Fergusson, R. Fernandez-Cobos, F. Finelli, F. Forastieri, M. Frailis, A. A. Fraisse, E. Franceschi, A. Frolov, S. Galeotta, S. Galli, K. Ganga, R. T. Génova-Santos, M. Gerbino, T. Ghosh, J. González-Nuevo, K. M. Górski, S. Gratton, A. Gruppuso, J. E. Gudmundsson, J. Hamann, W. Handley, F. K. Hansen, D. Herranz, S. R. Hildebrandt, E. Hivon, Z. Huang, A. H. Jaffe, W. C. Jones, A. Karakci, E. Keihänen, R. Keskitalo, K. Kiiveri, J. Kim, T. S. Kisner, L. Knox, N. Krachmalnicoff, M. Kunz, H. Kurki-Suonio, G. Lagache, J.-M. Lamarre, A. Lasenby, M. Lattanzi, C. R. Lawrence, M. L. Jeune, P. Lemos, J. Lesgourgues, F. Levrier, A. Lewis, M. Liguori, P. B. Lilje, M. Lilley, V. Lindholm, M. López-Caniego, P. M. Lubin, Y.-Z. Ma, J. F. Macías-Pérez, G. Maggio, D. Maino, N. Mandolesi, A. Mangilli, A. Marcos-Caballero, M. Maris, P. G. Martin, M. Martinelli, E. Martínez-González, S. Matarrese, N. Mauri, J. D. McEwen, P. R. Meinhold, A. Melchiorri, A. Mennella, M. Migliaccio, M. Millea, S. Mitra, M.-A. Miville-Deschênes, D. Molinari, L. Montier, G. Morgante, A. Moss, P. Natoli, H. U. Nørgaard-Nielsen, L. Pagano, D. Paoletti, B. Partridge, G. Patanchon, H. V. Peiris, F. Perrotta, V. Pettorino, F. Piacentini, L. Polastri, G. Polenta, J.-L. Puget, J. P. Rachen, M. Reinecke, M. Remazeilles, A. Renzi, G. Rocha, C. Rosset, G. Roudier, J. A. Rubiño-Martín, B. Ruiz-Granados, L. Salvati, M. Sandri, M. Savelainen, D. Scott, E. P. S. Shellard, C. Sirignano, G. Sirri, L. D. Spencer, R. Sunyaev, A.-S. Suur-Uski, J. A. Tauber, D. Tavagnacco, M. Tenti, L. Toffolatti, M. Tomasi, T. Trombetti, L. Valenziano, J. Valiviita, B. V. Tent, L. Vibert, P. Vielva, F. Villa, N. Vittorio, B. D. Wandelt, I. K. Wehus, M. White, S. D. M. White, A. Zacchei, A. Zonca. Planck 2018 results. vi. cosmological parameters. Astronomy & Astrophysics. Vol. 641, pg. A6, 2020 https://doi.org/10.1051/0004-6361/201833910 [↩] [↩] [↩] [↩] [↩]
  12. C. Pitrou, A. Coc, J.-P. Uzan, E. Vangioni. Precision big bang nucleosynthesis with improved helium-4 predictions. Physics Reports. Vol. 754, pg. 1–66, 2018 https://doi.org/10.1016/j.physrep.2018.04.005 [↩] [↩]
  13. B. D. Fields. The primordial lithium problem. Annual Review of Nuclear and Particle Science. Vol. 61, pg. 47–68, 2011 https://doi.org/10.1146/annurev-nucl-102010-130445 [↩]
  14. A. Coc, S. Goriely, Y. Xu, M. Saimpert, E. Vangioni. STANDARD big bang nucleosynthesis up to cno with an improved extended nuclear network. The Astrophysical Journal. Vol. 744, pg. 158, 2011 https://doi.org/10.1088/0004-637X/744/2/158 [↩]
  15. O. Richard, G. Michaud, J. Richer. Implications of wmap observations on li abundance and stellar evolution models. The Astrophysical Journal. Vol. 619, pg. 538–548, 2005 https://doi.org/10.1086/426470 [↩]
  16. B. J. Carr, S. W. Hawking. Black holes in the early universe. Monthly Notices of the Royal Astronomical Society. Vol. 168, pg. 399–416, 1974 https://doi.org/10.1093/mnras/168.2.399 [↩] [↩]
  17. S. W. Hawking. Black hole explosions? Nature. Vol. 248, pg. 30–31, 1974 https://doi.org/10.1038/248030a0 [↩]
  18. A. M. Green, B. J. Kavanagh. Primordial black holes as a dark matter candidate. Journal of Physics G: Nuclear and Particle Physics. Vol. 48, pg. 043001, 2021 https://doi.org/10.1088/1361-6471/abc534 [↩] [↩] [↩] [↩] [↩]
  19. B. J. Carr, S. W. Hawking. Black holes in the early universe. Monthly Notices of the Royal Astronomical Society. Vol. 168, pg. 399–416, 1974 https://doi.org/10.1093/mnras/168.2.399 [↩] [↩]
  20. B. J. Carr, K. Kohri, Y. Sendouda, J. Yokoyama. New cosmological constraints on primordial black holes. Physical Review D. Vol. 81, pg. 104019, 2010 https://doi.org/10.1103/PhysRevD.81.104019 [↩] [↩] [↩] [↩]
  21. K. Kohri, J. Yokoyama. Primordial black holes and primordial nucleosynthesis i: effects of hadron injection from low mass holes. Physical Review D. Vol. 61, pg. 023501, 1999 https://doi.org/10.1103/PhysRevD.61.023501 [↩] [↩]
  22. K. Jedamzik. Did something decay, evaporate, or annihilate during big bang nucleosynthesis? Physical Review D. Vol. 70, pg. 063524, 2004 https://doi.org/10.1103/PhysRevD.70.063524 [↩] [↩] [↩]
  23. A. M. Green, B. J. Kavanagh. Primordial black holes as a dark matter candidate. Journal of Physics G: Nuclear and Particle Physics. Vol. 48, pg. 043001, 2021 https://doi.org/10.1088/1361-6471/abc534 [↩] [↩] [↩]
  24. A. Arbey, J. Auffinger, K. P. Hickerson, E. S. Jenssen. AlterBBN v2: a public code for calculating big-bang nucleosynthesis constraints in alternative cosmologies. Preprint at http://arxiv.org/abs/1806.11095 2019 https://doi.org/10.48550/arXiv.1806.11095 [↩] [↩] [↩] [↩]
  25. Particle emission rates from a black hole: massless particles from an uncharged, nonrotating hole | phys. rev. d. https://journals.aps.org/prd/abstract/10.1103/PhysRevD.13.198 [↩] [↩] [↩]
  26. S. W. Hawking. Black hole explosions? Nature. Vol. 248, pg. 30–31, 1974 https://doi.org/10.1038/248030a0 [↩] [↩] [↩]
  27. B. J. Carr. Primordial black holes as a probe of cosmology and high energy physics. in Vol. 631 pg. 301–321, 2003 https://doi.org/10.1007/978-3-540-45230-0_7 [↩]
  28. E. W. Kolb, M. S. Turner. The Early Universe. Vol. 69 Taylor and Francis, 2019 https://doi.org/10.1201/9780429492860 [↩]
  29. M. Kawasaki, K. Kohri, T. Moroi. Big-bang nucleosynthesis and hadronic decay of long-lived massive particles. Physical Review D. Vol. 71, pg. 083502, 2005 https://doi.org/10.1103/PhysRevD.71.083502 [↩]
  30. M. Laine, Y. Schroder. Quark mass thresholds in qcd thermodynamics. Physical Review D. Vol. 73, pg. 085009, 2006 https://doi.org/10.1103/PhysRevD.73.085009 [↩]
  31. A. Boccia, F. Iocco, L. Visinelli. Constraining the primordial black hole abundance through big-bang nucleosynthesis. Physical Review D. Vol. 111, pg. 063508, 2025 https://doi.org/10.1103/PhysRevD.111.063508 [↩] [↩] [↩] [↩] [↩] [↩]
  32. C. Pitrou, A. Coc, J.-P. Uzan, E. Vangioni. Precision big bang nucleosynthesis with improved helium-4 predictions. Physics Reports. Vol. 754, pg. 1–66, 2018 https://doi.org/10.1016/j.physrep.2018.04.005 [↩]
  33. B. D. Fields, K. A. Olive, T.-H. Yeh, C. Young. Big-bang nucleosynthesis after planck. Journal of Cosmology and Astroparticle Physics. Vol. 2020, pg. 010–010, 2020 https://doi.org/10.1088/1475-7516/2020/03/010 [↩]
  34. A. Boccia, F. Iocco, L. Visinelli. Constraining the primordial black hole abundance through big-bang nucleosynthesis. Physical Review D. Vol. 111, pg. 063508, 2025 https://doi.org/10.1103/PhysRevD.111.063508 [↩] [↩]
  35. V. Poulin, P. D. Serpico, F. Calore, S. Clesse, K. Kohri. CMB bounds on disk-accreting massive primordial black holes. Physical Review D. Vol. 96, pg. 083524, 2017 https://doi.org/10.1103/PhysRevD.96.083524 [↩]
  36. M. Kawasaki, K. Kohri, T. Moroi. Big-bang nucleosynthesis and hadronic decay of long-lived massive particles. Physical Review D. Vol. 71, pg. 083502, 2005 https://doi.org/10.1103/PhysRevD.71.083502 [↩]

RELATED ARTICLESMORE FROM AUTHOR

LEAVE A REPLY

Please enter your comment!
Please enter your name here