Baryogenesis via Asymmetric Evaporation of Primordial Black Holes

TL;DR

This paper explores how primordial black holes can generate the universe's matter-antimatter asymmetry through asymmetric Hawking radiation influenced by higher-dimensional operators, accounting for entropy effects and detailed evolution.

Contribution

It introduces a comprehensive numerical framework for PBH baryogenesis, including entropy dilution, greybody factors, and various PBH mass spectra, advancing previous models.

Findings

Baryon asymmetry can be achieved with specific PBH mass spectra.

Entropy dilution effects are significant in the evolution.

Viable parameter space exists for multiple PBH mass distributions.

Abstract

We revisit baryogenesis from the asymmetric evaporation of light primordial black holes, focusing on scenarios where gravitational effects induce a matter antimatter asymmetry. In particular, we consider a higher-dimension operator coupling the Kretschmann scalar to a baryon-number-violating current which generates an effective chemical potential at the black hole horizon and leads to asymmetric Hawking radiation. Relative to earlier studies, we account for entropy dilution from evaporation, incorporate chemical potential dependent greybody factors and numerically track the fully coupled evolution of a PBH population in an expanding universe. We show that the observed baryon asymmetry can be reproduced within a viable region of parameter space for several PBH mass spectra including log-normal, critical-collapse, and power-law distributions.

Equations

$$\frac{1}{M_{\star}^{4}}\,\partial_{\alpha}\!\bigl{(}\mathcal{R}_{\mu\nu\rho\sigma}\mathcal{R}^{\mu\nu\rho\sigma}\bigr{)}\,J^{\alpha}\,,$$

$$\mathcal{K}\equiv {\mathcal{R}}_{\mu \nu \rho \sigma }{\mathcal{R}}^{\mu \nu \rho \sigma }=\frac{3M(t{)}^2}{4{\pi }^2{M}_{\mathrm{pl}}^4{r}^6},$$

$$\mu=\frac{\partial_{0}\mathcal{K}}{M^{4}_{*}}\bigg{\rvert}_{r=r_{\rm S}}=-\frac{3}{2}(8\pi)^{6}\varepsilon(M)M_{\mathrm{pl}}\left(\frac{M_{\mathrm{pl}}}{M}\right)^{7}\left(\frac{M_{\mathrm{pl}}}{M_{\star}}\right)^{4}\,,$$

$$f_{\mathrm{PBH}}(M)=\frac{1}{n_{\mathrm{PBH}}}\frac{d{n}_{\mathrm{PBH}}}{dM},$$

$${\int }_0^{\infty }{f}_{\mathrm{PBH}}(M)dM=1.$$

$${\rho }_{\mathrm{BH}}^{\mathrm{in}}=n_{\mathrm{PBH}}^{\mathrm{in}}{\int }_0^{\infty }M{f}_{\mathrm{PBH}}(M)dM,$$

$$M_{\mathrm{scl}}\equiv M_{\mathrm{scl}}(T_0)=\frac{4\pi }{3}\kappa \frac{{\rho }_R(T_0)}{H^3(T_0)},$$

$$T_0\simeq 4.3\times 10^{15}\mathrm{GeV}{\left(\frac{106.75}{g_{\ast }}\right)}^{1/4}{\left(\frac{1\text{leavevmode}\mathrm{g}}{M_{\mathrm{scl}}}\right)}^{1/2}{\left(\frac{\kappa }{0.2}\right)}^{1/2},$$

$$\beta \equiv \frac{{\rho }_{\mathrm{BH}}^{\mathrm{in}}(T_0)}{{\rho }_R(T_0)}.$$

$$f_{\mathrm{PBH}}(M)=\delta (M-M_{\mathrm{scl}}).$$

$$f_{\mathrm{PBH}}(M)=\frac{1}{\sqrt{2\pi }\sigma M}\exp \left[-\frac{{\log }^2(M/M_{\mathrm{scl}})}{2{\sigma }^2}\right],$$

$$f_{\mathrm{PBH}}(M)\propto \{\begin{array}{ll}{M}^{-\alpha },& \text{for }{M}_{\mathrm{scl}}\le M\le M_{\mathrm{scl}}\times 10^{\sigma },\\ 0,& \text{otherwise}.\end{array}$$

$$\alpha \equiv \frac{4w+2}{w+1}.$$

$$1<\alpha \le 3.$$

$$f_{\mathrm{PBH}}(M)\propto M^{1.25}\exp \left[-{\left(\frac{M}{M_{\mathrm{scl}}}\right)}^{2.85}\right],$$

$$f_{\mathrm{PBH}}(M)=\frac{\eta }{\mu \Gamma (\kappa /\eta )}{\left(\frac{M}{\mu }\right)}^{\kappa -1}\exp \left[-{\left(\frac{M}{\mu }\right)}^{\eta }\right],$$

$$r_{\mathrm{S}}\sim 2.7\times 10^{-5}\left(\frac{M}{10^{10}\mathrm{g}}\right)\mathrm{fm},$$

$$T_{\mathrm{BH}}=\frac{M_{\mathrm{pl}}^2}{M}\simeq 1\text{leavevmode}\mathrm{TeV}\left(\frac{10^{10}\mathrm{g}}{M}\right).$$

$$\frac{d^2{n}_i}{dpdt}=\frac{g_i}{2{\pi }^2}\frac{{\sigma }_{s_i}\left(M,p,{\mu }_i\right)}{\exp \left[\left(E_i(p)+{\mu }_i\right)/T_{\mathrm{BH}}\right]-(-1{)}^{2s_i}}\frac{p^3}{E_i(p)}.$$

$${\Gamma }_{s_i}\equiv \frac{{\sigma }_{s_i}{p}^2}{\pi },$$

$$\frac{dM}{dt}=-\sum _i{\int }_{m_i}^{\infty }\frac{d^2{n}_i}{dtdE}EdE=-\frac{(8\pi {)}^2\epsilon (M)M_{\mathrm{pl}}^4}{M^2},$$

$$\tau \simeq \frac{160}{\pi g_{\ast }}\frac{M_{\mathrm{in}}^3}{M_{\mathrm{pl}}^4},$$

$$T_{\mathrm{ev}}^{(\mathrm{rad})}\simeq 1.2\times 10^{10}\mathrm{GeV}{\left(\frac{g_{\ast }}{106.75}\right)}^{1/4}{\left(\frac{1\mathrm{g}}{M_{\mathrm{in}}}\right)}^{3/2}.$$

$$T_{\mathrm{ev}}^{(\mathrm{PBH})}\simeq 4.1\times 10^{10}\mathrm{GeV}{\left(\frac{g_{\ast }}{106.75}\right)}^{5/12}{\left(\frac{10^{-3}}{\beta }\right)}^{1/3}{\left(\frac{0.2}{\kappa }\right)}^{1/6}{\left(\frac{1\mathrm{g}}{M_{\mathrm{in}}}\right)}^{11/6}.$$

$$\Gamma \equiv \int dp\left[\frac{d^2{n}_{\overline{L}}}{dpdt}-\frac{d^2{n}_L}{dpdt}\right]=\frac{8\pi M_{\mathrm{pl}}^2}{M}\int d\omega \left[\frac{d^2{n}_{\overline{L}}}{dpdt}-\frac{d^2{n}_L}{dpdt}\right]=\frac{8\pi M_{\mathrm{pl}}^2}{M}{\tilde{\Gamma }}_{\mathrm{BH}\to \mathrm{\Delta L}},$$

$$\Gamma \approx -\frac{g_L{\sigma }_0\mu T_{\mathrm{BH}}^2}{6},$$

$$\Gamma \approx -{\sigma }_0\frac{g_L{T}_{\mathrm{BH}}^3}{{\pi }^2}\left({\mathrm{Li}}_3(-e^{-\mu /T_{\mathrm{BH}}})-{\mathrm{Li}}_3(-e^{\mu /T_{\mathrm{BH}}})\right).$$

$$\frac{n_L}{s}=\frac{n_{\mathrm{PBH}}}{s}\delta N_L={\Omega }_{\mathrm{PBH}}\frac{{\rho }_{\mathrm{tot}}}{s}\frac{\delta N_L}{M_{\mathrm{ev}}}=\frac{3{\Omega }_{\mathrm{PBH}}{T}_{\mathrm{ev}}\delta N_L}{4M_{\mathrm{ev}}},$$

$$\delta N_L={\int }_{M_{\min }}^{M_{\mathrm{PBH}}^0}dM\frac{{\dot{n}}_L-{\dot{n}}_{\bar{L}}}{\int dpp\frac{d^2{n}_L}{dpdt}+\int dpp\frac{d^2{n}_{\bar{L}}}{dpdt}+\int dpp\frac{d^2{n}_{\mathrm{other}}}{dpdt}}.$$

$$\delta N_L=\frac{1}{2}{\left(\frac{M_{\mathrm{cr}}}{M_{\mathrm{pl}}}\right)}^2{\int }_{X_{\min }}^{X_0}dXg(X).$$

Full text

1 Introduction

The universe contains more matter than antimatter, with a baryon-to-entropy ratio, $n_{\mathrm{B}}/s\simeq 9\times 10^{-11}$ [1]. Explaining this asymmetry remains an open problem in fundamental physics and many mechanisms have been proposed (see Refs. [2, 3, 4, 5]). Two of the most studied are electroweak baryogenesis [6] and leptogenesis [7]. In the absence of CPT violation, such scenarios must satisfy Sakharov’s conditions [8]: baryon-number violation, C and CP violation and departure from thermal equilibrium. In electroweak baryogenesis, anomalous sphaleron processes provide baryon-number violation, while CP-violating sources, typically arising in extensions of the Standard Model, operate near the walls of bubbles during a first-order electroweak phase transition. The bubble expansion supplies the required non-equilibrium environment. Leptogenesis creates a lepton asymmetry through the out-of-equilibrium, CP-violating decays of heavy Majorana neutrinos that violate lepton number. Standard Model sphalerons then partially convert this lepton asymmetry into a baryon asymmetry. In contrast, in baryogenesis from primordial black holes (PBHs), the out-of-equilibrium condition is inherent in the Hawking evaporation and we consider the possibility that a higher-dimension operator coupling the Kretschmann scalar to a baryon-number–violating current can generate an effective chemical potential at the horizon, biasing the emitted Hawking flux.

Alternative baryogenesis scenarios involving black holes were also explored by Hawking [9] and Carr [10] in the 70s. In particular, they studied the evaporation of light PBHs into nucleons that subsequently decay asymmetrically producing more baryons than antibaryons. Baryogenesis via PBH evaporation has since attracted a lot of attention [11, 12, 13, 14, 15, 16, 17]. Recently, several authors [18, 19, 20, 21] revisited the role of PBHs in the context of baryon asymmetry. Specifically, they studied a scenario where a chemical potential around the PBH is dynamically generated due to the presence of a higher order ($d>4$) operator, biasing Hawking evaporation towards baryon over antibaryon production. Several mechanisms can induce such an effect, including spontaneous baryogenesis [22, 23], gravitational fluctuations produced during inflation [24] or the mass loss of an evaporating PBH [19]. In this work, we consider the final scenario, similar to Ref. [19], where it was pointed out that as Hawking radiation steadily reduces the black hole mass, $M=M(t),$ the curvature invariants that are a function of $M$, become time-dependent. In particular, let us consider the following dimension-eight, CP-violating operator

[TABLE]

where $M_{\star}$ is the cut-off scale of the new physics that induces the CP-violating coupling, $\mathcal{R}_{\mu\nu\rho\sigma}$ is the Riemann tensor and $J^{\alpha}$ is the baryon (or lepton) number–violating current. This operator is non–vanishing whenever the time derivative of the Kretschmann scalar, $\mathcal{K}\equiv\mathcal{R}_{\mu\nu\rho\sigma}\mathcal{R}^{\mu\nu\rho\sigma}$, is non-zero. The Kretschmann scalar outside a PBH is given by

[TABLE]

where $r$ is the Schwarzschild radius, $M_{\mathrm{pl}}\equiv(8\pi G)^{-1/2}\simeq 2.43\times 10^{18}\,\mathrm{GeV}$ is the reduced Planck mass with $G$ being Newton’s constant, and one finds that $\dot{\mathcal{K}}\neq 0$ whenever $\dot{M}\neq 0$. This non-vanishing time derivative induces the operator $(\partial_{\alpha}\mathcal{K})J^{\alpha}/M_{\star}^{4}$ to be non-zero and creates an effective CP-violating chemical potential at the black hole event horizon [19],

[TABLE]

where $\varepsilon(M_{\rm BH})$ is the mass evaporation function of the black hole and $r_{\rm S}$ is the event horizon radius. The evaporation function encapsulates the total degrees of freedom available for emission at a given black hole mass, as described in Ref. [25]. The induced chemical potential grows as the black hole evaporates and produces a baryon (or lepton) asymmetry. In Ref. [19], the authors find a viable region in parameter space where the evaporation of light PBHs can account for the observed baryon asymmetry. However, the results are obtained under a set of simplifying assumptions. First, the authors only consider the evaporation of PBHs into spin [math] and spin $1/2$ particles, and second, they neglect the entropy dilution coming from the evaporation into photons. Third, the dependence of the greybody factors on the effective chemical potential is also neglected, and finally, the baryonic yield is obtained from an approximate analytic treatment rather than a full solution of the coupled evolution equations for an evaporating PBH population in an expanding universe.

In this work, we build on the results of Ref. [19], by going beyond the approximated analytical analysis and investigating the interplay between this CP-violating operator and PBHs via numerical simulations. Our paper is laid out as follow: we review the physics of PBHs in Sec. 2 and provide some analytical understanding of the lepton asymmetry generated from a PBH in Sec. 3. In Sec. 4, we compare the analytical approach with our numerical setup and present our results and conclusions in Secs. 5 and 6, respectively. In Appendix A, we provide the details of the computation of the greybody factors including the chemical potential, which we use in our numerical simulations.

2 Primordial Black Holes

We first review four representative BH mass spectra: monochromatic, log-normal, power-law, and (generalized) critical-collapse distributions in Sec. 2.1 and follow in Sec. 2.2 by detailing the evaporation of microscopic Schwarzschild BHs.

2.1 Formation and Mass Spectrum

Primordial black holes can form in the early Universe through the collapse of overdensities generated by mechanisms such as inflation [26, 27, 28], first-order phase transitions [29, 30, 31, 32, 33, 34, 35, 36], or topological defects [37, 38, 39, 40, 41, 42, 43, 44, 45]. We consider a population of Schwarzschild PBHs formed during the radiation-dominated era, potentially characterized by an extended initial mass distribution, $dn_{\rm PBH}/dM$. We define the normalized initial number distribution as

[TABLE]

such that it satisfies the normalization condition

[TABLE]

From this distribution, the initial comoving PBH energy density is given by

[TABLE]

where $n_{\rm PBH}^{\rm in}$ denotes the PBH number density at formation. This number density can be related to the properties of the primordial plasma. We assume that the distribution $f_{\mathrm{PBH}}(M)$ is centered around a characteristic mass scale $M_{\rm scl}$, which we associate with the mass contained within the particle horizon at the time of formation, corresponding to a radiation temperature $T_{0}$

[TABLE]

where $\kappa\simeq 0.2$ parametrizes the efficiency of collapse, $\rho_{R}(T_{0})$ is the radiation energy density, and $H(T_{0})$ is the Hubble parameter at temperature $T_{0}$. Using Eq. (2.3), the formation temperature for a PBH of mass $M_{\rm scl}$ is approximately

[TABLE]

where $g_{*}$ is the effective number of relativistic degrees of freedom at the time of formation. This implies that gram-scale PBHs produced from inflationary fluctuations must have formed at extremely high temperatures, $T_{0}\gtrsim 10^{15}\,\mathrm{GeV}$.

The fractional contribution of PBHs to the total energy density at formation is defined as

[TABLE]

The exact form of the PBH mass distribution $f_{\mathrm{PBH}}(M)$ depends on the specific production mechanism. We will consider the following four well-known cases for the normalized mass distribution.

• •

Monochromatic. A simple scenario corresponds to consider that all PBH were generated with the same mass $M_{\rm scl}$ such that

[TABLE]

Although unrealistic, this monochromatic scenario is helpful to understand the main features of PBH evaporation and its generation of baryon asymmetry.

• •

Log-Normal (LN). The formation of PBHs during inflation typically requires a brief phase of ultra-slow-roll dynamics, which amplifies the primordial power spectrum of scalar curvature perturbations and generates a localized peak [46, 47, 48, 49]. This enhancement generally results in a log-normal PBH mass function [50], a feature that has been confirmed both numerically and analytically in the context of slow-roll inflation [51, 52]. The resulting mass distribution takes the form

[TABLE]

where $M_{\rm scl}$ denotes the central mass of the distribution and $\sigma$ characterizes its width.

• •

Power-Law (PL). An alternative mechanism for PBH formation involves the collapse of a nearly scale-invariant spectrum of primordial curvature perturbations in a Universe dominated by a perfect fluid with a constant equation-of-state parameter $w$ [53]. In this scenario, the resulting PBH mass distribution follows a power-law form,

[TABLE]

where the exponent $\alpha$ is determined by

[TABLE]

The mass range $[M_{\rm scl},\,M_{\rm scl}\times 10^{\sigma}]$ is set by the range of scales over which the scale-invariant spectrum extends. In physically relevant scenarios where the Universe is no longer inflating, the equation-of-state parameter typically lies in the range $-1/3<w\leq 1$, which corresponds to

[TABLE]

• •

Critical Collapse (CC). The application of critical scaling to gravitational collapse provides a more refined description of PBH formation from primordial density fluctuations [54, 55, 56, 57]. While early models assumed that PBHs form with masses equal to the horizon mass at the time of collapse, critical phenomena reveal a more complex picture: there exists an upper cutoff near the horizon mass, but the mass distribution extends to smaller values, forming a characteristic low-mass tail. This behavior has been shown to arise generically across a wide range of inflationary models [57]. The resulting PBH mass function takes the form

[TABLE]

where $M_{\rm scl}$ denotes the characteristic mass scale of the distribution. However, one can generalize the CC mass function to better reproduce the PBH mass distribution generated by several representative primordial curvature spectr, see Ref [58]. Such generalized critical collapse (GCC) mass function is parametrized as

[TABLE]

with $\eta>0,\kappa>1$ [58, 59] are the parameters of the GCC distribution, such that it peaks at $M_{\rm scl}=\mu[(\kappa-1)/\eta]^{1/\eta}$.

Having considered the initial distribution of our PBH population, we can now focus on its possible evaporation in the Early Universe.

2.2 Evaporation and Hawking Radiation

The Schwarzschild radius, $r_{\mathrm{S}}$, of a light black hole can be microscopic,

[TABLE]

and when this radius becomes comparable to the Compton wavelength of a particle species, the PBH can emit particles of that species as Hawking radiation. The temperature associated with this Hawking radiation is inversely proportional to the black hole mass:

[TABLE]

As a black hole evaporates it converts its mass into a spectrum of Hawking quanta. The differential emission rate, per unit time and per unit momentum, for a particle species $i$ of mass $m_{i}$, spin $s_{i}$, internal degrees of freedom $g_{i}$ and chemical potential $\mu_{i}$ is

[TABLE]

The cross-section $\sigma_{s_{i}}(M,p,\mu_{i})$ is often recast as the greybody factor:

[TABLE]

which quantifies how the gravitational and centrifugal barriers surrounding the black hole partially reflect outgoing quanta, thereby modifying the Hawking spectrum [9, 60, 61, 62].

In Figure 1 we show the difference of the differential emission rate between leptons and antileptons. As expected, in the limit where the chemical potential is negligible ($\mu/T_{\rm BH}\ll 1$), the differential rate for both leptons and antileptons is roughly the same, or, in other words, Hawking evaporation is symmetric. Therefore, the differential emission rate asymmetry is negligible. On the other hand, as the value of $\mu/T_{\rm BH}$ becomes significant, one can see that Hawking evaporation yields more antileptons than leptons. The large suppression at large values of $\omega$ comes from the exponential function in Eq. (2.15).

Integrating the differential rate of Eq. (2.15) over momentum and summing over all particle species $i$ (for simplicity we also neglect the chemical potential), one finds the net mass loss is [63, 64, 65]

[TABLE]

where $t$ denotes time and $M$ the instantaneous black hole mass. Let us note that in order to obtain the mass lose rate of Eq. (2.17), we have assumed the validity of the semi-classical approximation, in which quantum fields propagate on a fixed classical spacetime, that is, neglecting backreaction effects111The “memory burden” hypothesis [66] proposes late-time deviations from the standard mass loss rate, motivated by analogies with solitonic configurations in non-Abelian gauge theories and assuming the N-graviton portrait of a BH. As it remains speculative, we do not include these effects in our analysis.. This treatment is expected to fail near the Planck scale, where the adiabaticity condition $\dot{T}_{\rm BH}/T_{\rm BH}^{2}\ll 1$ breaks down [67]. To remain conservative, we will consider the effect of stopping the evolution at $M=M_{\rm pl}$ and $M\approx 100\,M_{\rm pl}$, showing the effect on the generated baryon asymmetry. Further complications arise from the information loss problem. After the Page time [68, 69, 70], the von Neumann entropy of Hawking radiation exceeds the Bekenstein-Hawking entropy, suggesting potential non-thermal corrections to the emission spectrum. As the PBHs of interest have long surpassed this point, such effects could in principle modify the spectrum [70]. However, in the absence of a well-established framework for these corrections [71, 72], we assume Eq. (2.17) as valid throughout the PBH lifetime until near the Planck scale.

Integrating the mass loss rate of Eq. (2.17) from the initial mass $M_{\mathrm{in}}$ down to $100\cdot M_{\rm pl}$ gives the total lifetime of a PBH,

[TABLE]

where we assumed that $M_{\mathrm{in}}\gg 100\cdot M_{\rm pl}$ and $\varepsilon(M)\simeq g_{*}/30720\pi$. Using this lifetime and assuming the Universe continues to be radiation dominated, we can compute the radiation temperature when the PBH population completely evaporates [13]:

[TABLE]

If, on the other hand, PBH evaporation occurs when the Universe’s energy density is PBH-dominated, the radiation temperature at evaporation now depends on the PBH’s abundance:

[TABLE]

Both Eqs. (2.19) and (2.20) are approximations that describe the Standard Model plasma just before the PBHs evaporate suddenly.

3 Lepton asymmetry generated from a PBH

An unbiased black hole emits leptons via Hawking radiation at a rate determined by the lepton’s internal degrees of freedom and the energy-dependent greybody factor (see e.g. Eq. (2.15)). Consequently, antileptons are produced at the same rate and the net lepton number vanishes. However, the presence of the curvature–induced operator of Eq. (1.1) creates a local chemical potential $\mu(M)$ that can favor the production of leptons over anti-leptons, or vice versa. To quantify the bias, we define the following lepton asymmetry rate:

[TABLE]

where we have rescaled the momentum to the dimensionless variable $\omega\equiv Mp/(8\pi M_{\mathrm{pl}}^{2})$. In our numerical analysis, we work with the dimensionless rate $\widetilde{\Gamma}_{\rm BH\to\Delta L}$. We will now consider the lepton asymmetry rate in several limits. First, we will assume that $\mu\ll T_{\mathrm{BH}}$ and that the absorption cross-section for particle and antiparticle are identical and momentum independent such that $\sigma_{s_{i}}\left(M,p,\mu\right)=\sigma_{0}$. To first order in $\mu/T_{\rm BH}$, we find that the lepton asymmetry rate is

[TABLE]

and hence the lepton asymmetry, increases linearly with $\mu$. The second analytic limit we consider is that simply the absorption cross-section is momentum and chemical potential independent and we find that the lepton asymmetry rate can be approximated as

[TABLE]

For large $\mu/T_{\rm BH}$, the polylogarithms grow approximately as $(\mu/T_{\mathrm{BH}})^{3}/6$ and therefore the asymmetry rises faster than linearly, $\Gamma\propto\mu^{3}$, at large chemical potential. Nonetheless, for a more accurate estimation of the lepton asymmetry from evaporation, we compute the greybody factors, including the dependence on the chemical potential $\sigma_{1/2}(M,p,\mu)$ and without relying on any assumption for a wide range of values for $\mu/T_{\rm BH}$ (see discussion in Appendix A).

We compare the lepton asymmetry production rate for different greybody factor prescriptions in Figure 2. Specifically, we show the dependence of the dimensionless rate $\widetilde{\Gamma}_{\rm BH\to\Delta L}$, assuming for $\sigma_{0}$ the geometric optics limit $EM/M_{\mathrm{pl}}^{2}\gg 1$ (green), the low-energy approximation $EM/M_{\mathrm{pl}}^{2}\ll 1$ (orange), and the full numerical greybody factors (purple). The bottom panel displays the ratio of the two limiting cases, geometric optics (olive) and low energy (light blue), to the exact result. All three prescriptions exhibit the same qualitative behavior: for $\mu/T\gg 1$, the rate scales as $(\mu/T_{\rm BH})^{3}$, while for $\mu/T_{\rm BH}\ll 1$, it scales linearly with $\mu/T_{\rm BH}$. As expected, the two limiting cases, which assume $\sigma_{0}$ to be independent of $\mu/T_{\rm BH}$, differ from each other by an overall normalization. However, the full rate shows a slightly different dependence on $\mu/T_{\rm BH}$. In the $\mu/T_{\rm BH}\ll 1$ regime, the geometric optics approximation overestimates the correct rate due to its significant overestimation of the greybody factors. This discrepancy diminishes near $\mu/T_{\rm BH}\sim 1$, where contributions from higher angular momentum modes become relevant. As seen in Fig. 1, for some range of momenta $p$, the condition $\mu/T_{\rm BH}>E_{i}(p)/T_{\rm BH}$ is satisfied, reducing the exponential suppression typically affecting higher-$\ell$ modes in the absence of a chemical potential. As a result, the exact rate exceeds the geometric optics and low-energy approximations by a factor of approximately 1.15 and 15, respectively. This highlights the critical importance of using accurate greybody factors when evaluating lepton asymmetry production in this context. In the following sections, we use the exact greybody factors for our study of lepton asymmetry production.

4 Friedmann and Boltzmann Equations for PBH Baryogenesis

In this section, we begin in Sec. 4.1 by reviewing the analytic treatment on how a baryon asymmetry can be generated from a population of biased PBHs and follow in Sec. 4.2 by presenting our numerical treatment of the problem.

4.1 Analytic Treatment

To compute the lepton asymmetry produced by a population of PBHs in the early Universe, given the CP- and lepton number $L$-violating operator of Eq. (1.1), one must take into account the complete set of coupled equations describing both PBH and radiation energy density evolution in an expanding Universe. This level of detail was not fully addressed in Ref. [19], which instead approximated the lepton asymmetry as:

[TABLE]

where $\Omega_{\rm PBH}$ is the ratio of the PBH energy density to the total energy density at the time of evaporation (when the Universe’s temperature is $T_{\rm ev}$), $M_{\rm ev}$ is the PBH mass at evaporation, and $\rho_{\rm tot}/s=3\,T_{\rm ev}/4$. The quantity $\delta N_{L}$ estimates the lepton asymmetry generated by a single black hole and is obtained by integrating the lepton asymmetry rate relative to the radiated energy carried per each type of emitted particle over the lifetime of the black hole:

[TABLE]

Changing the integration variable from $M$ to $X=(M/M_{\rm cr})^{2}$ (see Eq. (A.15)), we can rewrite Eq. (4.2) as

[TABLE]

Note that the value of $X_{\rm min}$ will be in general different from zero because as the PBH evaporates both the black hole temperature and chemical potential increase and will eventually reach the cutoff scale $M_{\star}$, beyond which our effective field theory breaks down. Therefore, $X_{\rm min}$ will be given by the following condition

[TABLE]

where the critical mass $M_{\mathrm{cr}}$ is set by the condition $\mu\sim T_{\rm BH}$, which corresponds to the point in the black hole’s evolution at which the lepton asymmetry reaches its maximum. This quantity is defined to be

[TABLE]

One can check that $X_{\rm min}=1$ for $M_{\star}=10^{-3}M_{\rm pl}$, for which both the temperature and chemical potential reach $M_{\star}$ at the same time. For $M_{\star}<10^{-3}M_{\rm pl}$ ($X_{\rm min}>1$), it is the PBH temperature that sets the limit, while for $M_{\star}>10^{-3}M_{\rm pl}$ ($X_{\rm min}<1$), it is the chemical potential that will reach the cutoff scale first. As discussed in Ref. [19], this provides a conservative estimate of the lepton asymmetry, since, in principle, the PBH would still evaporate asymmetrically beyond this point.

The behavior of $g(X)$ is shown in the left panel of Figure 3. We see that $g(X)$ peaks at around $X\sim 0.6$, or, equivalently, $M\sim 0.8M_{\rm cr}$. This is the point when the PBH most efficiently produces a lepton asymmetry. As explained in Ref. [19], for larger values of $X$ the chemical potential becomes small, and Hawking evaporation is essentially symmetric, whereas for smaller values, although the bias in Hawking evaporation is maximal, the emitted particles (in this case, leptons or antileptons) carry a lot of energy.

In deriving Eq. (4.1) and Eq. (4.2), it is assumed that PBHs only evaporate into massless scalars and fermions, using the low-energy limit for the greybody factors. Moreover, entropy dilution, which is important for estimating both dark matter production and the baryon asymmetry in PBH scenarios [12, 73, 74, 13], is neglected. Following the procedure in Ref. [19], we reproduce the results shown in Figure 3 right, which shows the resulting lepton yield (see colour bar) as a function of the cutoff scale $M_{\star}$ and initial PBH mass $M$ for $\Omega_{\rm PBH}=1$. The grey triangle in the lower left corner indicates the region of parameter space where PBHs have either a temperature (below the yellow line) or a chemical potential (below the pink line) above $M_{\star}$ at formation, and hence the effective field theory is immediately invalid. The grey-shaded band highlights regions of parameter space excluded by inflation222This bound only applies to PBH formation mechanisms related to inflation, such as the collapse of large overdensities upon horizon reentry. However, we note that effects such as critical collapse can lead to arbitrarily light PBHs. In any case, we do not assume any particular formation mechanism and consider that PBHs above the Planck mass form with a mass distribution of the kind discussed in Section 2.1. [75, 76]. Note that there is a clear change of behavior at $M_{\star}\simeq 10^{-3}M_{\mathrm{pl}}$. This is due to the integral cutoff $X_{\rm min}$ in Eq. (4.3). For a fixed PBH mass, as we decrease the cutoff scale $M_{\star}$, the point of maximum lepton asymmetry production is reached earlier (see Eq. (4.5)). At the same time, from Eq. (4.4), the time at which our effective theory breaks down is also reached earlier. Therefore, two opposing effects are at play as we decrease the value of $M_{\star}$. Due to the different dependence of the breakdown point on $M_{\star}$, we observe two different behaviors above and below $M_{\star}\simeq 10^{-3}M_{\mathrm{pl}}$.

We can now estimate the impact of entropy dilution by following Ref. [13]. Using energy conservation before and after PBH evaporation, we find that

[TABLE]

where $n^{0}_{\rm PBH}=n_{\rm PBH}(T_{0})$ is the initial PBH number density, $\tilde{T}$ is the temperature of the SM plasma after PBH evaporation and $\beta=\rho_{\rm PBH}(T_{0})\!/\!\rho_{\rm tot}(T_{0})$, with $\rho_{\rm PBH}(T_{0})$ and $\rho_{\rm tot}(T_{0})$ denoting the energy densities of the PBHs and the total energy of the Universe at formation, respectively. We also define all relevant quantities at PBH formation time since the PBH mass is constant for most of its evolution and $n_{\rm PBH}(T_{\rm ev})/s(T_{\rm ev})\simeq n_{\rm PBH}(T_{0})/s(T_{0})$. Applying the entropy dilution factor, the estimated lepton asymmetry is

[TABLE]

where we used $\rho_{\rm tot}/s=3\,T_{0}/4$. We present the resulting leptonic yield in Figure 4. Although the inclusion of entropy dilution suppresses the lepton asymmetry relative to Eq. (4.1), there is an additional factor of $T_{0}/T_{\rm ev}$ that partially compensates for this suppression. We observe a similar behavior than in Figure 3.

4.2 Numerical Solution

Having gained this insight, we now move on to a more detailed numerical treatment of the interplay between PBHs and the radiation component, closely following the approach in Refs. [12, 77]. First, we write down the mass loss rate for the PBH population,

[TABLE]

where $\varrho_{\mathrm{BH}}\equiv a^{3}\,\rho_{\mathrm{BH}}$ is the comoving energy density of the PBHs, $a$ is the scale factor (initialized to $a_{0}=1$), $M_{\rm in}^{t}\equiv M_{\rm in}(M,t)$ is the initial mass that evolves to $M$ at a time $t$. Finally, $dM_{\rm in}/dM$ is the Jacobian that transforms the initial mass to the mass at a time $t$. It can be readily computed (see Refs. [78, 77]):

[TABLE]

where we have changed the variable from time to scale factor according to, $d/dt\equiv aHd/da$ as $H=\dot{a}/a$. Since PBHs behave like matter, their energy density is given by the product of the PBH number density and the mass, with the latter evolving in time. Additionally, we include equations for the comoving radiation energy density as well as the Hubble parameter

[TABLE]

where $\epsilon_{\rm SM}(M)$ includes only the Standard Model contribution to the evaporation and

[TABLE]

with $g_{*}(T)$ and $g_{*S}(T)$ denoting the effective and entropic numbers of relativistic degrees of freedom, respectively.

Since PBH evaporation injects additional entropy into the Universe, the overall entropy is not conserved. Therefore we track the evolution of the ambient temperature $T$ following Refs. [79, 80, 81, 77]:

[TABLE]

We solve Eqs. (4.10) and (4.12) as a function of the scale factor, typically stopping at $M_{\rm BH}\approx 100\cdot M_{\rm Pl}$. In practice, the black hole will stop generating a lepton asymmetry before reaching this limit as the effective field theory breaks down beyond a certain scale determined by the condition given in Eq. (4.4). We adopt a conservative prescription by assuming symmetric Hawking evaporation from when the effective field description becomes invalid. This choice is conservative because, in reality, the chemical potential would continue increasing until the black hole fully evaporates, thereby amplifying the evaporation bias toward leptons (or antileptons) and yielding a larger final lepton asymmetry.

The rate of change of the PBH mass and plasma temperature, $T$, will be unchanged from previous calculations 333In principle one should take into account consistently the evolution of the chemical potential as well. However, it will be negligible for most of the PBH evolution time and only become significant at the very end before reaching the breakdown point. We have verified numerically that it only induces a negligible correction. . The evolution of the lepton asymmetry is

[TABLE]

where $\widetilde{\Gamma}_{\rm BH\to\Delta L}$ (see Eq. (3.1)) is shown in Figure 2, $n_{\overline{L}-L}$ denotes the leptonic asymmetry number density. As PBH evaporation occurs significantly before the electroweak phase transition, electroweak sphalerons will convert the lepton asymmetry to a baryonic one.

One can ask if the lepton asymmetry gets washed out or removed via competing processes. If the lepton asymmetry is generated at energies below the scale of the new physics responsible for neutrino masses, then the dominant washout processes are $\Delta L=2$ interactions. Their cross-section can be approximated as

[TABLE]

resulting in a rate

[TABLE]

Comparing this rate to the Hubble expansion parameter, one finds that the washout processes fail to reach equilibrium whenever

[TABLE]

assuming $m_{\nu}\simeq 0.1\,\text{eV}$. Hence, as long as the Universe’s temperature remains below this temperature during the bulk of the lepton-asymmetry production, the washout never becomes efficient, and it can safely be ignored in the Boltzmann equation of Eq. (4.13).

5 Results

To obtain the total lepton asymmetry, we solve the complete set of Eqs. (4.10) and (4.13) accounting for the full evolution of the Universe with the exact greybody factors for different PBH mass distributions. Let us discuss the monochromatic scenario first. In Figure 5, we scan over the parameter space given by the cutoff scale $M_{\star}$ and the PBH mass at formation. We show the results for a value of the initial abundance $\beta=10^{-3}$ (left) and $\beta=10^{-5}$ (right) where the black line indicates the central value of the baryon asymmetry. We observe the same qualitative features as in the results in Figure 4 with the analytic treatment. Namely, a constant leptonic yield for a fixed initial PBH mass and cutoff $M_{\star}<10^{-3}M_{\rm pl}$, and a decrease at larger values of the cutoff scale. Also, larger initial PBH masses lead to a lower final leptonic yield and vice versa. Indeed, for a fixed abundance, a population of heavier PBHs corresponds to a lower number density, hence there will be less black holes contributing to the lepton asymmetry. Also, there will be a larger entropy dilution, as heavier PBHs will convert a larger fraction of their mass into radiation not contributing to the asymmetry. There is however an overall noticeable difference. We observe a very significant increase of the leptonic yield below $M\sim 10^{6.8}M_{\rm pl}$ ($M\sim 10^{8.8}M_{\rm pl}$) for $\beta=10^{-3}$ ($\beta=10^{-5}$), and up to a factor $\sim 4$ ($\gtrsim 10$) for the lowest PBH initial masses. We have checked that for small initial PBH masses the entropy dilution becomes irrelevant and the difference is due to the use of the exact greybody factors. On the other hand, the full numerical analysis results in a gradual reduction of the leptonic yield as we move to larger PBH initial masses, and up to 70$\%$ less compared to the analytical values. For these values of the PBH mass, even if the enhancement from the exact greybody factor still applies, the entropy dilution becomes dominant. The difference is then due to the different prescriptions used to compute the entropy dilution, namely Eq. (4.6) in the analytical treatment versus the full numerical setup.

In Figure 5 we evolve the PBH masses down to the Planck mass, after which Hawking evaporation is assumed to halt. We show the results when changing this prescription and evolving down to a final mass $M^{f}_{\rm PBH}=100M_{\rm pl}$ for $\beta=10^{-5}$ in Figure 6 left. Since we impose a cutoff at a fixed final mass $M^{f}$, the numerical evolution shows that, for some initial PBH masses, the evolution stops before reaching $M_{\rm min}$, while for others the cutoff is reached first. Our algorithm iteratively solves the PBH system in steps, reducing the mass from its initial value to 10% of that value. As a result, the point at which either $M_{\rm min}$ or $M^{f}_{\rm PBH}$ is reached varies periodically, producing the zigzag pattern observed in this figure. Therefore, apart from this numerical artifact, no significant differences are seen between Figure 5 and Figure 6.

We also scan over several values of $\beta$ for a given choice of the cutoff scale $M_{\star}=10^{-4}M_{\rm pl}$ in Figure 6 right. Interestingly, the resulting leptonic yield becomes independent of the PBH abundance $\beta$ at large initial PBH masses as well as for $\beta\gtrsim 10^{-5}$ (above the dashed line in Figure 6 right). In this region $1\ll\beta T_{0}/T_{\rm ev}$, and the universe goes through an epoch of early PBH domination before PBH decay. Therefore, from the analytical treatment, in particular Eq. (4.7) and Eq. (2.20), one can see that for a fixed initial PBH mass the resulting leptonic yield indeed becomes independent since $T_{\rm ev}\propto\beta^{-1/3}$. On the other hand, when $1\gg\beta T_{0}/T_{\rm ev}$, the entropy dilution is negligible and the leptonic yield increases linearly with $\beta$.

Moving on to extended mass distributions, in Figures. 7 and 8, we show the results for two example scenarios with a log-normal and a generalised critical collapse mass distribution, respectively, and compare them with the monochromatic case (dashed line). The right panel shows the energy density of radiation and PBH in blue and black, respectively. In the left panel, we show the leptonic yield $n_{\text{L}}/s$. In both cases the Universe presents an early phase of PBH domination starting around $\log(a)\sim 5$, as observed in the right panel of Figure 8. Note that this phase is also present for a monochromatic mass function for large values of $M_{\rm scl}$. The duration of this epoch is shorter for the critical collapse case, due to the less extended high-mass tail compared to the log-normal mass function, which corresponds to heavier PBHs that evaporate later. On the other hand, the enhanced low-mass tail leads to a larger leptonic yield already at very early times compared to the log-normal scenario. In both cases, the leptonic yield increases initially at a similar pace as in the monochromatic case, essentially given by the linear increase of the lepton asymmetry rate discussed in Figure 2. As the lighter PBHs complete their evaporation, the leptonic yield gradually stabilises at the observed value $n_{L}/s\sim 10^{-10}$. The resulting curve is smoother compared to the monochromatic case as PBHs completely evaporate at different moments over some extended period of time rather than having a population of PBHs evaporating at the same instant, where the dashed line suddenly increases by several orders of magnitude. Note that there is a slight decrease in the very end due to the effect of heavier PBHs producing a large amount of radiation-like particles that do not contribute to the lepton asymmetry. We also observe the same qualitative features with a power-law mass distribution, for which the observed baryon asymmetry can also be reproduced.

6 Conclusion

In this paper we have revisited the role of light PBHs in the context of baryogenesis (and leptogenesis). In particular, we studied the evaporation of a population of PBHs in the framework of an effective field theory with a dimension-eight, CP-violating operator coupling the derivative of the Kretschmann scalar with a baryon (or lepton) number-violating current, see (1.1). The evolution of the PBH mass generates a non-zero time derivative that induces an effective chemical potential $\mu$ on the PBH horizon, biasing Hawking evaporation towards more baryons (or leptons) over antibaryons (or antileptons).

We have improved upon previous results on several fronts. Firstly, we performed our computations using $\mu$-dependent greybody factors. As shown in Figure 2, the resulting lepton asymmetry rate of a single PBH with exact greybody factors can differ substantially from results obtained by approximating the greybody factors in the high- or low-energy limits. Secondly, we considered the full PBH evaporation spectrum and, in particular, the contribution from evaporation into photons, which had previously been neglected. We find this component to play a crucial role through entropy dilution, which constrains the parameter space capable of reproducing the observed baryon asymmetry, as illustrated in Figure 6 (right). Finally, we consistently solved the complete set of equations describing the full evolution of the Universe, including radiation and a PBH population characterised by an extended mass distribution. We presented two examples in Figure 7 and Figure 8, corresponding to a critical collapse and a log-normal mass function, respectively. In both cases, the observed baryon asymmetry is successfully reproduced.

In conclusion, our findings demonstrate that the observed baryon asymmetry can be explained by a biased evaporating population of light PBHs. Our results hold for the three types of mass distributions we considered: critical collapse, log-normal, and power-law. Moreover, based on our study of the monochromatic mass function, we expect our conclusions to remain valid across a wide range of initial PBH abundances $\beta$ and effective field theory cutoffs $M_{\star}$, with entropy dilution representing the main physical limitation. At the same time, as discussed in Ref. [19], our results should be regarded as conservative, since PBHs are expected to continue generating baryon asymmetry beyond the scale $M_{\star}$, where our effective description ceases to be applicable. A complete UV theory could therefore open additional regions of parameter space capable of accounting for the observed baryon asymmetry. Further work is needed to explore the possible origin of the CP-violating higher-dimensional operators.

Acknowledgments

JT would like to thank the Quantum Field Theory Centre at the University of Southern Denmark for their hospitality during the completion of this work. YFPG has been supported by the Consolidación Investigadora grant CNS2023-144536 from the Spanish Ministerio de Ciencia e Innovación (MCIN) and by the Spanish Research Agency (Agencia Estatal de Investigación) through the grant IFT Centro de Excelencia Severo Ochoa No CEX2020-001007-S.

Appendix A Greybody Factors for fermions including chemical potential

Let us consider the Dirac equation in the Schwarzschild background, including a new interaction in the form of a chemical potential. The action, following Ref. [19] is

[TABLE]

where ${\cal\mathcal{K}}=R_{\alpha\beta\gamma\delta}R^{\alpha\beta\gamma\delta}$ is the Kretschmann scalar for the Schwarzschild metric, given by

[TABLE]

where $M$ is the instantaneous BH mass. The Dirac equation in the two-component spinor notation is given by [82]

[TABLE]

where $\nabla_{BB^{\prime}}$ is the covariant differentiation in te pseudo-Riemannian geometry related to the Schwarzschild spacetime, $A_{BB^{\prime}}$ is the vector field associated to the chemical potential, $\mu_{e}$ is the fermion mass, $P^{B},Q^{B}$ are the two component spinors, such that

[TABLE]

We use the Newman-Penrose formalism with a null tetrad ($\mathbf{l,n,m,M_{\star}}$), we get the system of equations, see Ref. [82]

[TABLE]

For the Schwarzschild spacetime, one has the contravariant four vectors

[TABLE]

where $\Delta=r(r-M/(4\pi M_{\mathrm{pl}}^{2}))$. The components relevant for the Dirac equation are

[TABLE]

with $\omega$ represents the energy. Taking $A_{\mu}=-(A_{0},0,0,0)$, $A_{0}$ is the time component of the chemical potential four vector, we have

[TABLE]

Separating $r$ and $\theta$ variables by defining

[TABLE]

Here $R_{\pm 1/2}(r)$ and $S_{\pm 1/2}(\theta)$ are functions depending only on the radial and angular variables. Since the new interaction assumed here is spherically symmetric, we have that the angular functions $S_{\pm 1/2}(\theta)$ are the same as in the standard Schwarzschild case, see e.g. Ref. [83, 84]. To compute the greybody factors, we need to solve the equations for $R_{\pm 1/2}(r)$. Letting $R_{-1/2}\equiv R_{1}$ and $R_{+1/2}=(2/\Delta)^{1/2}R_{2}$

[TABLE]

where $K=(\omega-A_{0})r^{2}$, and $\lambda=j+1/2$, $j\in[1/2,3/2,5/2,\ldots]$, the angular eigenvalue. Performing the same transformation as in Ref. [62], $R_{1}=G+iF$, $R_{2}=G-iF$, we get

[TABLE]

Since the shape of these equations is the same as the ones given in Ref. [62], to obtain the greybody factors, we follow the procedure detailed in the Appendix of Ref. [62], with the appropriate changes due to the additional potential $A_{0}$, given by

[TABLE]

where

[TABLE]

being $r_{S}=M/(4\pi M_{\mathrm{pl}}^{2})$ the Schwarzschild radius and $\mu\left.\mu\right|_{r=r_{S}}$ is the chemical potential evaluated at the horizon, as given in Eq. (21) of Ref. [19]. When performing the variable transformation given in Ref. [62], we find that we need to work with the dimensionless variable $\widetilde{\alpha}$,

[TABLE]

which is proportional to $\left.\mu\right|_{r=r_{S}}/T_{\rm BH}$ given in the same paper, apart from a factor of $4\pi$.

The maximum value for $\left.\mu\right|_{r=r_{S}}/T_{\rm BH}$ can in principle go way beyond 1 as the black hole evaporates and its mass falls below $M_{\rm cr}$. An absolute upper bound can be derived by assuming the maximum value for the cutoff scale $M_{\star}$ and the minimum value for the black hole mass $M_{\rm BH}$ for which our current physical theories still apply. This corresponds to $M_{\star}=M_{\rm p}$ and $M_{\rm BH}=M_{\rm p}$, which leads to $\left.\mu\right|_{r=r_{S}}/T_{\rm BH}\sim 10^{6}$, where we used Eq.21 and 22 given in Ref. [19]. However, this value will not be reached because, as the black hole evaporates, the temperature and the chemical potential increase until $T_{\rm BH}=M_{\star}$ and/or $\left.\mu\right|_{r=r_{S}}=M_{\star}$ is satisfied, beyond which the effective theory is no longer valid. For a given value of the cut-off $M_{\star}$, we can only track the black hole evolution down to a mass given by Eq. 28 in Ref. [19]. This sets the maximum value for $\left.\mu\right|_{r=r_{S}}/T_{\rm BH}$ for that particular value of $M_{\star}$. In Figure 9 we show the maximum values in the range of interest for $M_{\star}$.