# Redshift Effects in Particle Production from Kerr Primordial Black Holes

Andrew Cheeka<EMAIL_ADDRESS>Lucien Heurtierb<EMAIL_ADDRESS>Yuber F. Perez-Gonzalezb<EMAIL_ADDRESS>Jessica Turnerb
<EMAIL_ADDRESS>a Astrocent, Nicolaus Copernicus Astronomical
Center of the Polish Academy of Sciences, ul.Rektorska 4, 00-614 Warsaw,
Poland  b Institute for Particle Physics Phenomenology, Durham University,
South Road, Durham, U.K.

###### Abstract

When rotating primordial black holes evaporate via Hawking radiation, their
rotational energy and mass are dissipated with different dynamics. We
investigate the effect of these dynamics on the production of dark radiation –
in the form of hot gravitons or vector bosons – and non-cold dark matter.
Although the production of higher-spin particles is enhanced while primordial
black holes are rotating, we show that the energy density of dark radiation
experiences an extra redshift because their emission effectively halts before
PBH evaporation completes. We find that taking this effect into account leads
to suppression by a factor of $\mathcal{O}(10)$ of $\Delta N_{\rm eff}$ for
maximally rotating black holes as compared to previous results. Using the
solution of the Friedmann and Boltzmann equations to accurately calculate the
evolution of linear perturbations, we revisit the warm dark matter constraints
for light candidates produced by evaporation and how these limits vary over
black hole spins. Due to the interplay of enhanced production and late
dilution, we obtain that higher spin particles are most affected by these
bounds. Our code FRISBHEE, FRIedmann Solver for Black Hole Evaporation in the
Early universe, developed for this work can be found at .

††preprint: IPPP/22/49

##  I Introduction

Determining the evolution of the early Universe before Big-Bang
Nucleosynthesis (BBN) era is not an easy task. Although observational data
regarding the Cosmic Microwave Background (CMB) may constrain very early
stages of the Universe evolution such as cosmic inflation, very little is
known about the post-inflationary phase of our Universe’s history. Indeed, the
thermalization of Standard Model (SM) particles at early time is expected to
erase all traces of the initial conditions which led to the hot big bang era
we are familiar with. The discovery of Gravitational Waves (GWs) Abbott  _et
al._ (2016, 2017) opened the door to new ways of probing the evolution of the
early Universe, either through the detection of recent astrophysical events or
through the measurement of a stochastic background of GWs produced in the
primordial Universe (see, e.g. Ref. Auclair  _et al._ (2022) and references
therein for a review). Historically, however, measurements of the number of
effective relativistic degrees of freedom $N_{\rm eff}$ has provided important
information regarding the Universe before BBN, for example, by confirming the
description of neutrino decoupling. In the realm of physics beyond the SM,
measurements of $N_{\rm eff}$ provide vital constraints on models with light
particles. Many different models predict a deviation of $N_{\rm eff}$ from the
minimal situation where only neutrinos and photons dominate the Standard Model
(SM) radiation energy density, see e.g. Refs. Jeong and Takahashi (2013);
Buen-Abad  _et al._ (2015); Lesgourgues  _et al._ (2016); Buen-Abad  _et al._
(2018); Blinov and Marques-Tavares (2020); Chacko  _et al._ (2004, 2005);
Berlin and Blinov (2018).

A particularly appealing scenario that has attracted renewed attention is the
formation and evaporation of Primordial Black Holes (PBHs) Zel’dovich and
Novikov (1967); Hawking (1971); Carr and Hawking (1974). PBHs could have
dominated the evolution of the Universe before the BBN epoch, leading to
observational tests such as GWs Inomata  _et al._ (2020); Domènech  _et al._
(2021) or contributions to $N_{\rm eff}$ as PBHs emit all existing degrees of
freedom in nature Hawking (1974, 1975) including light states that could
modify the number of relativistic degrees of freedom. Among the many Beyond
the Standard Model (BSM) light particles predicted in diverse scenarios, one
candidate stands out, the massless spin-2 two mediator of Gravity, the
_graviton_. As many consider the graviton’s existence to be guaranteed, one
could, in principle, constrain the presence of PBHs in the pre-BBN Early
Universe by measuring either modification to $N_{\rm eff}$, parameterized by
$\Delta N_{\rm eff}$, or GWs in future experiments Dong  _et al._ (2016).
However, let us note that if the PBHs were of the Schwarzschild type, graviton
emission is suppressed compared to other degrees of freedom. The rationale
behind such diminished production lies in the dependence of the Hawking
emission rate on absorption probabilities, or Graybody Factors, which are
reduced at low energies for larger values of the total angular momentum of the
incoming — or emitted — particle.

If, on the other hand, PBHs had an initial non-zero angular momentum, the
emission of higher spin particles becomes enhanced Page (1976a, b, 1977),
leading to an augmented number density of emitted gravitons, which in turn
could lead to observable modifications to $N_{\rm eff}$, as first studied in
Refs. Hooper  _et al._ (2019, 2020); Masina (2020, 2021); Arbey  _et al._
(2021). These works computed the contribution to $\Delta N_{\rm eff}$
considering the time evolution of both mass and spin and demonstrated that
future CMB-HD Aiola  _et al._ (2022) experiments could, in principle, probe
the existence of PBHs. In this paper, we revisit the derivation of $\Delta
N_{\rm eff}$ for hot gravitons produced by the evaporation of spinning — Kerr
— PBHs. Since the BH angular momentum is depleted much faster than its mass,
the enhanced emission of spin-2 particles only occurs in the first stages of
evaporation, such that graviton production effectively stops before the PBHs
have entirely evaporated. Thus, there is an additional and crucial redshift
that gravitons experience, which, in the past has been neglected or partially
accounted for, when calculating $\Delta N_{\rm eff}$ We thus improve on
previous treatments by carefully tracking the time evolution of PBH mass and
angular momentum, together with the Universe’s evolution. We assume a
monochromatic mass and spin distribution and consider that the PBH population
acquired a significant angular momentum via some unspecified mechanism. We
find a difference of a factor of $\mathcal{O}(10)$ to previous estimates for
close-to-maximally rotating PBHs in the determination of $\Delta N_{\rm eff}$.
Moreover, we determine the sensitivity to the initial PBH density fraction as
a function of the initial PBH mass in the range $10^{-1}\lesssim M_{\rm
BH}^{\rm in}(~{}{\rm g})\lesssim 10^{9}$ for gravitons and new light vectors
from future measurements of $\Delta N_{\rm eff}$.

The unique dynamics that Kerr PBHs possess will also produce important
modifications in the generation of additional and stable BSM states, which
could constitute the dark matter (DM) Fujita  _et al._ (2014); Lennon  _et
al._ (2018); Morrison  _et al._ (2019); Hooper  _et al._ (2019); Masina
(2020); Baldes  _et al._ (2020); Gondolo  _et al._ (2020); Bernal and Zapata
(2021a, b); Auffinger  _et al._ (2021); Bernal and Zapata (2021c); Arbey  _et
al._ (2021); Masina (2021); Jyoti Das  _et al._ (2021); Cheek  _et al._
(2022a, b); Sandick  _et al._ (2021); Bernal  _et al._ (2022). Constraints
from structure formation, especially for light DM, where $m_{\rm DM}\lesssim
10^{3}\,\rm{GeV}$, significantly impact the allowed parameter space, as
demonstrated in Baldes  _et al._ (2020); Auffinger  _et al._ (2021); Masina
(2021). In the second part of this work, we return to the derivation of such
limits on light DM by solving the complete set of Friedmann-Boltzmann
equations together with angular momentum and mass depletion rates to provide
accurate inputs for calculating the evolution of linear perturbations in the
Universe. For the latter task, we use CLASS Lesgourgues (2011); Blas  _et al._
(2011); Lesgourgues and Tram (2011) to obtain updated warm dark matter
constraints. Moreover, we perform a systematic study of how these constraints
vary over black hole spins, finding that higher spin particles are most
affected, where the same story of enhanced production and late time dilution
plays out, but with the additional effect of redshift acting to slow the
matter particles.

This manuscript is organized as follows. In Section II, we describe the time
evolution of Kerr PBHs, emphasizing the angular moment evaporation. Next, we
detail the determination of $\Delta N_{\rm eff}$ for gravitons and additional
light vectors in Section III, comparing with existing results in the
literature. Finally, Section IV is devoted to the specific scenario of light
DM emission from Kerr PBHs, where we determine the redshifted phase-space
distribution of DM and the subsequent rederivation of hot DM constraints. We
draw our conclusions in Section V. Throughout, we assume natural units where
$\hbar=c=k_{\rm B}=1$.

##  II Evaporation of Kerr PBHs

When PBHs evaporate, they emit particles of all kinds at rates that differ
depending on the particle’s quantum numbers and the mass and angular momentum
of the evaporating PBH. As consequence, their mass and angular momentum are
depleted at rates that depend on the initial PBH properties. To derive such
evaporation equations, let us consider the Hawking emission rate of any
particle species $i$, which within a time ($\differential t$) and momentum
($\differential p$) interval, given by Hawking (1974, 1975)

$\displaystyle\frac{\differential^{2}\mathcal{N}_{i}}{\differential
p\differential t}$
$\displaystyle=\frac{g_{i}}{2\pi^{2}}\sum_{l=s_{i}}\sum_{m=-l}^{l}\frac{\differential^{2}\mathcal{N}_{ilm}}{\differential
p\,\differential t}\,,$ (1)

with

$\displaystyle\frac{\differential^{2}\mathcal{N}_{ilm}}{\differential
p\differential t}$ $\displaystyle=\frac{\sigma_{s_{i}}^{lm}(M_{\rm
BH},p,a_{\star})}{\exp\left[(E_{i}-m\Omega)/T_{\rm
BH}\right]-(-1)^{2s_{i}}}\frac{p^{3}}{E_{i}}\,,$ (2)

where $M_{\rm BH}$ is the instantaneous PBH mass, $E_{i}$ ($g_{i}$) is the
energy (internal degrees of freedom) of particle $i$,
$\Omega=(a_{\star}/2GM_{\rm BH})(1/(1+\sqrt{1-a_{\star}^{2}}))$ is the angular
velocity of the horizon, $G$ is Newton’s constant, and $l,m$ are the total and
axial angular momentum quantum numbers, respectively. An essential quantity in
Eq. (2) is the absorption cross-section, $\sigma_{s_{i}}^{lm}$, which
characterizes the possible back-scattering of particles in the presence of
centrifugal and gravitational potentials Hawking (1974, 1975); Page (1976a,
1977). Such cross-sections are determined from reflection and transmission
coefficients of scatterings by the gravitational potential after solving
equations of motions in a curved spacetime around the BH. For the specific
case of spinning BHs, we have followed the method established in Refs.
Chandrasekhar and Detweiler (1975, 1976, 1977), where the problem is reduced
to solve Schrödinger-like one-dimensional wave equations whose potentials
depend on the particle’s spin.

The BH mass and spin loss rates are calculated by summing Eq. (1) over the
different species and integrating over the phase space, to obtain MacGibbon
and Webber (1990); MacGibbon (1991)

$\displaystyle\frac{dM_{\rm BH}}{dt}$ $\displaystyle=-\epsilon(M_{\rm
BH},a_{\star})\frac{M_{p}^{4}}{M_{\rm BH}^{2}}\,,$ (3)
$\displaystyle\frac{da_{\star}}{dt}$ $\displaystyle=-a_{\star}[\gamma(M_{\rm
BH},a_{\star})-2\epsilon(M_{\rm BH},a_{\star})]\frac{M_{p}^{4}}{M_{\rm
BH}^{3}}\,,$ (4)

where $M_{p}=1.22\times 10^{19}\,{\rm GeV}$ denotes the Planck mass,
$a_{\star}=JM_{\rm Pl}^{2}/M_{\rm BH}^{2}$ is the dimensionless spin parameter
of the BH, and $\epsilon(M_{\rm BH},a_{\star})$ ($\gamma(M_{\rm
BH},a_{\star})$) is the evaporation (angular momentum evaporation) function
which depends on the BH instantaneous mass and spin and the total degrees of
freedom existing in nature. For a given particle type $i$, the functions
$\gamma_{i}(M_{\rm BH},a_{\star})$ and $\varepsilon_{i}(M_{\rm BH},a_{\star})$
are

$\displaystyle\varepsilon_{i}(M_{\rm BH},a_{\star})$
$\displaystyle=\frac{g_{i}}{2\pi^{2}}\int_{0}^{\infty}\sum_{l=s_{i}}\sum_{m=-l}^{l}\frac{\differential^{2}\mathcal{N}_{ilm}}{\differential
p\differential t}\,EdE\,,$ (5) $\displaystyle\gamma_{i}(M_{\rm BH},a_{\star})$
$\displaystyle=\frac{g_{i}}{2\pi^{2}}\int_{0}^{\infty}\sum_{l=s_{i}}\sum_{m=-l}^{l}m\frac{\differential^{2}\mathcal{N}_{ilm}}{\differential
p\differential t}\,dE\,,$ (6)

where the details on how to calculate these parameters can be found in Ref.
Cheek  _et al._ (2022a).

Fig. 1: The purple curve shows the evolution of $\widetilde{M}_{\rm BH}^{\rm
max}(\xi)\equiv M_{\rm BH}(\xi)/M_{\rm BH}^{\rm in}$, the BH mass normalised
to initial mass while the green dashed curve presents the angular momentum
depletion $a_{\star}^{\rm max}(\xi)$ as a function $\xi$, the time normalised
to the BH total lifetime. The presented solutions are obtained for an
initially maximally rotating BH ($a_{\star,\,\rm in}=1$). The time evolution
for any other $a_{\star,\,\rm in}<1$ can be derived from these curves after
performing a suitable transformation, see text.

Let us now consider the time evolution of a close-to-maximally rotating BH,
having an initial spin parameter $a_{\rm\star,in}\to 1$ and an initial mass
$M_{\rm BH}^{\rm in}$. In Fig. 1, we present the numerical solution to the
system of equations where $f(\xi)$ denotes alternatively $\widetilde{M}_{\rm
BH}^{\rm max}(t)\equiv M_{\rm BH}(\xi)/M_{\rm BH}^{\rm in}$ (purple) or
$a_{\star}^{\rm max}(\xi)$ (green dashed), being $\xi\equiv t/\tau\in[0,1]$,
with $\tau$ denoting the PBH lifetime. As seen from this figure, the spin of a
Kerr BH is generally depleted via Hawking evaporation much faster than the
mass. Indeed, at around $t\sim 0.65\tau$, one can observe that the spin has
dropped over one order of magnitude ($a_{\star}\sim 0.1$), whereas the mass
has only decreased by a half ($M\sim 0.5M_{\rm in}$). As we will see in the
proceeding sections, this rapid decrease of the spin, compared to the mass,
implies that a careful treatment of Hawking evaporation, which simultaneously
takes into account the evolution of the Universe, is required when dealing
with a particle production that is strongly spin-dependent.

Interestingly, it was shown in Ref. Page (1976b) that, in the limit of light
evaporation products, the solution of Eqs. (3)-(4) can be generically obtained
for any initial mass and spin from the evolution of a close-to-maximally
rotating BH, $a_{\rm\star,in}=1$. To obtain the generic solution, the
following transformations must be performed:

$\displaystyle a_{\star}(\xi)$ $\displaystyle=a_{\star}^{\rm
max}([1-x_{i}]\xi+x_{i})\,,$ (7a) $\displaystyle\widetilde{M}_{\rm BH}(\xi)$
$\displaystyle=\frac{\widetilde{M}_{\rm BH}^{\rm
max}([1-x_{i}]\xi+x_{i})}{\widetilde{M}_{\rm BH}^{\rm max}(x_{i})}\,,$ (7b)

where $x_{i}$ is such that $a_{\star}^{\rm max}(x_{i})=a_{\star,\,\rm in}$. We
then can extend the previous discussion about the time dependence of higher
spin particle emission. For the case where $a_{\star,\,\rm in}=0.7$, average
spin value from hierarchical merger rates Fishbach  _et al._ (2017), we find
that only $\sim 40\%$ of the initial mass has been depleted when
$a_{\star}=0.1$, demonstrating that the importance of correctly tracking the
evolution of the PBH mass and spin with respect of the overall evolution of
the Universe.

##  III Determination of $\Delta N_{\rm eff}$

The number of relativistic degrees of freedom in the early Universe is
extremely sensitive to the thermal history of the Universe’s particle physics
content. All SM particles, including neutrinos, are in thermal equilibrium at
high temperatures. As the temperature of the thermal plasma decreases, SM
particles become non-relativistic and decoupled from the thermal bath,
releasing their energy to lighter species through annihilation or decay
processes. Neutrinos, however, decouple while relativistic, effectively
following a Fermi-Dirac distribution, despite their interaction rate becoming
inefficient compared to the Universe’s expansion rate. Assuming that the
neutrinos decoupled instantaneously, the subsequent decoupling of the
electrons from the thermal bath is expected to inject entropy exclusively in
the photon bath, which effectively shifts the neutrino temperature as compared
with the photon temperature:

$T_{\nu}=(4/11)^{1/3}T_{\gamma}\,.$ (8)

In the absence of any further considerations, one would expect the energy
density of radiation today to follow the relation

$\rho_{\mathrm{R}}^{\mathrm{SM}}=\rho_{\gamma}\left[1+\frac{7}{8}\left(\frac{T_{\nu}}{T_{\gamma}}\right)N_{\mathrm{eff}}^{\rm
SM}\right]\,,$ (9)

where $\rho_{\gamma}$ is the energy density of the photon bath and
$N_{\mathrm{eff}}^{\rm SM}$ is the number of relativistic neutrinos, which
would naively be $3$. However, this is in fact $N_{\mathrm{eff}}\approx 3.045$
Mangano  _et al._ (2002); de Salas and Pastor (2016). The reason for this
extra contribution is twofold: $(i)$ electrons transfer a fraction of entropy
into the neutrino sector through QED corrections, and, most importantly,
$(ii)$ the decoupling of neutrinos is not instantaneous. The amount of entropy
transferred from the electron sector to neutrinos out-of-equilibrium is
temperature-dependent. Therefore, the non-instantaneous nature of neutrino
decoupling is critical as it affects how much entropy the electrons can
transfer to the neutrino sector after decoupling.

In the event that a dark sector may contribute to the energy density of
radiation today, one usually defines this extra amount of relativistic degrees
of freedom, $\Delta N_{\mathrm{eff}}$, as

$\rho_{\mathrm{R}}\equiv\rho_{\gamma}\left[1+\frac{7}{8}\left(\frac{T_{\nu}}{T_{\gamma}}\right)(N_{\mathrm{eff}}^{\rm
SM}+\Delta N_{\mathrm{eff}})\right]\,.$ (10)

Denoting the energy density of dark radiation as $\rho_{\mathrm{DR}}$ and
using Eq. (9) together with the fact that
$\rho_{\mathrm{R}}=\rho_{\mathrm{R}}^{\mathrm{SM}}+\rho_{\mathrm{DR}}$ allows
to invert this relation to obtain

$\Delta N_{\rm
eff}\equiv\left\\{\frac{8}{7}\left(\frac{4}{11}\right)^{-\frac{4}{3}}+N_{\rm
eff}^{\rm SM}\right\\}\frac{\rho_{\rm DR}}{\rho_{\rm R}^{\rm SM}}\,,$ (11)

where both $\rho_{\mathrm{R}}^{\mathrm{SM}}$ and $\rho_{\mathrm{DR}}$ are
evaluated at present time. Computing the contribution of PBH evaporation to
$\Delta N_{\mathrm{eff}}$ therefore requires the precise evaluation of the
energy density of dark radiation, $\rho_{\mathrm{R}}$, produced during the
evaporation process.

The possibility that PBHs contribute significantly to $\Delta
N_{\mathrm{eff}}$ was thoroughly studied by different groups Masina (2020,
2021); Hooper  _et al._ (2020); Arbey  _et al._ (2021); Hooper  _et al._
(2019). Those studies showed that the Hawking production of particles with
spin $\geqslant 1$ becomes significantly enhanced in the presence of a large
PBH spin. In the context of dark matter production from PBH evaporation, this
same conclusion was reached in Refs. Cheek  _et al._ (2022a, b) in a complete
numerical analysis, including the full grey body factors and the dynamical
solving of Boltzmann and Friedman equations.

By nature, the evaporation of PBHs is not an instantaneous process and,
similarly to the decoupling of neutrinos in the SM, it is natural to question
whether the time evolution of the evaporation process for Kerr PBHs may affect
predictions for the value of $\Delta N_{\mathrm{eff}}$. As we have seen in
Section II, the dynamics of the evaporation are encoded in a two-dimensional
differential system, describing the non-trivial evolution of both the mass and
spin of the BH as a function of time (see Eq. (3) and (4)). Spinning PBHs tend
to lose most of their spin in the first phase of the evaporation before the
remaining mass evaporates in the second phase. These complicated dynamics have
two major effects: $(i)$ the total lifetime of a spinning PBH is affected by
the initial value of its spin, and $(ii)$ the production of vector and tensor
particles is significantly enhanced during the first phase of the evaporation
(when the spin is non-zero) but becomes similar to the one of a Schwarzschild
PBH towards the end of the evaporation process. In Refs. Arbey  _et al._
(2021); Hooper  _et al._ (2020, 2019) the first of those two effects was
considered via calculation of the PBH lifetime in the presence of the BH spin.
However, instantaneous evaporation was assumed when computing different
particle abundances, affecting the contribution to $\Delta N_{\mathrm{eff}}$.
To go beyond this approximation, we solve the Friedmann equations
simultaneously with the evaporation dynamical system of Eq. (3) and (4). We
will assume PBHs were formed after overdense perturbations entered the Hubble
horizon during a radiation-dominated era. This common assumption leads to the
generic one-to-one relation between the initial mass of the PBH and the
initial energy density of the Universe at the time of PBH formation

$\displaystyle M_{\rm BH}^{\rm
in}=\frac{4\pi}{3}\gamma\frac{\rho_{i}}{H_{i}^{3}}\,,$ (12)

where $\gamma=0.2$ which is a dimensionless gravitational collapse parameter
and $H_{i}$ is the Hubble parameter Carr  _et al._ (2020). For further
convenience, we will quantify the energy fraction of PBHs in the Universe, at
the time of their formation, using the parameter

$\displaystyle\beta^{\prime}\equiv\gamma^{1/2}\left(\frac{g_{\star}(T_{\rm
in})}{106.75}\right)^{-1/4}\frac{\rho^{\rm in}_{\rm PBH}}{\rho^{\rm in}}\,.$
(13)

The set of equations that are to be solved numerically consists of the
Friedmann equation

$\frac{3H^{2}M_{p}^{2}}{8\pi}=\rho_{\mathrm{R}}^{\rm SM}+\rho_{\rm
DR}+\rho_{\rm PBH}\,,$ (14)

together with the set of Boltzmann equations describing the evolution of the
SM radiation (R), dark radiation (DR) and PBH energy densities throughout the
history of the Universe are respectively given by Gutiérrez  _et al._ (2018);
Cheek  _et al._ (2022a)

$\displaystyle\dot{\rho}_{\mathrm{R}}^{\rm SM}+4H\rho_{\mathrm{R}}^{\rm SM}$
$\displaystyle=-\left.\frac{\differential\log M_{\rm BH}}{\differential
t}\right|_{\rm SM}\rho_{\rm PBH}\,,$ (15a) $\displaystyle\dot{\rho}_{\rm
DR}+4H\rho_{\rm DR}$ $\displaystyle=-\left.\frac{\differential\log M_{\rm
BH}}{\differential t}\right|_{\rm DR}\rho_{\rm PBH}\,,$ (15b)
$\displaystyle\dot{\rho}_{\rm PBH}+3H\rho_{\rm PBH}$
$\displaystyle=\frac{\differential\log M_{\rm BH}}{\differential t}\rho_{\rm
PBH}\,,$ (15c)

respectively. We determine the solutions of the Friedmann, Boltzmann and BH
evolution equations using our code called FRIedmann Solver for Black Hole
Evaporation in the Early universe, FRISBHEE, which will be publicly available
in 111A detailed instruction manual for the usage of FRISBHEE will be released
in the near future..

Fig. 2:  Evolution of the energy density of PBH (blue) and SM radiation
(black) in the case of Schwarzschild black holes (left) and Kerr black holes
with spin $a_{\star}=0.99$ (right). The bottom panels depict the evolution of
the graviton contribution to the total energy density of radiation during the
evaporation process.

Once the evaporation has completed, the SM evolves independently from DR, and
the energy densities of the two sectors at present time can easily be
extrapolated from their values after evaporation using entropy conservation.
The value of $\Delta N_{\mathrm{eff}}$ can therefore be calculated using
Hooper  _et al._ (2019)

$\displaystyle\Delta N_{\rm
eff}=\left\\{\frac{8}{7}\left(\frac{4}{11}\right)^{-\frac{4}{3}}+N_{\rm
eff}^{\rm SM}\right\\}\frac{\rho_{\rm DR}(T_{\rm ev})}{\rho_{\rm R}^{\rm
SM}(T_{\rm ev})}\left(\frac{g_{*}(T_{\rm ev})}{g_{*}(T_{\rm
eq})}\right)\left(\frac{g_{*S}(T_{\rm eq})}{g_{*S}(T_{\rm
ev})}\right)^{\frac{4}{3}}\,,$ (16)

where $N_{\rm eff}^{\rm SM}=3.045$ denotes the effective number of
relativistic neutrinos in the SM de Salas and Pastor (2016), $T_{\rm ev}$ is
the plasma temperature when the PBHs have evaporated, $T_{\rm eq}=0.75$ eV is
the temperature at which matter-radiation equality occurs.

###  III.1 Graviton contributions to $\Delta N_{{\rm eff}}$

Despite having no hopes of detecting the graviton anytime soon, its existence
is often assumed when quantum field theorists consider gravitation. As any
massless spin-2 field would give rise to gravity, this further cements the
belief that such a particle is likely to exist. Given the graviton’s status as
an ‘honorary’ member of the Standard Model, one can view the constraints on
its energy density produced by evaporating BHs as a pure test of the possible
PBH abundance in the early Universe.

To understand the effect of redshifting during the non-instantaneous
evaporation, we present in Fig. 2 an example of the evolution of the energy
densities for a PBH with $M_{\rm BH}^{\rm in}=10^{5}$ g and assuming
$\beta^{\prime}=10^{-3}$. In the left (right) panel, we show the evolution of
the energy density of the SM radiation and PBHs in blue and black,
respectively, for a Schwarzschild (Kerr with $a_{\star}=0.99$) PBHs. In the
bottom panels, the energy density of the gravitons relative to the SM
radiation is shown. In the Schwarzschild case, the PBHs come to dominate the
energy density of the Universe at $a\sim 10^{3}$ and their evaporation becomes
relevant at $a\sim 10^{12}$; however, the SM radiation energy density begins
to increase around $a\sim 10^{9}$. We see that the fraction of the radiation
energy density consisting of gravitons increases simultaneously as the visible
radiation energy density increases and then plateaus as the gravitons
contribute a constant fraction towards the total radiation energy density. For
Kerr PBHs, we observe that gravitons are amply produced while the PBHs spins
are non-zero as the fractional energy density of the gravitons to the total
radiation increases sharply (is overall larger than in the Schwarzschild case)
and plateaus. However, in the final stage of evaporation, when the rotational
energy has been depleted, the PBH remains massive and continues to evaporate.
Meanwhile, the gravitons are free-streaming thus only experiencing redshift.
Therefore, the graviton contribution to the total radiation bath is reduced as
photon production dominates (shown by the steep decrease after
$\log_{10}(a)\sim 11$ in the lower panel).

Fig. 3:  PBH contribution to $\Delta N_{\rm eff}$ assuming only the Standard
Model and a massless graviton. We assume PBH domination and monochromatic BH
distributions in mass $M_{\rm in}$ and spin $a_{\star}$. The solid lines
correspond to our calculation, solving for the number densities in a time
evolving manner. The teal, purple and orange lines correspond to
$a_{\star}=0.0,0.7$ and $0.99$ respectively. We overlay sensitivities of
current and future experimental sensitivities to $\Delta N_{\rm eff}=0.027$ at
the $2\sigma$ C.L. from CMB-HD Aiola  _et al._ (2022). We also compare our
results to Ref. Hooper  _et al._ (2020), Ref. Arbey  _et al._ (2021), and Ref.
Masina (2021) referred to as Paper A, Paper B, and Paper C, respectively. Also
included are the results to our calculation neglecting the effects of redshift
and when, in the BHs life gravitons at emitted, see text for details.

In Fig. 3, we show the predictions for $\Delta N_{\rm eff}$ that we obtain
when primordial black holes of various spins dominate the energy density of
the Universe evaporate and produce massless gravitons as they evaporate (solid
lines). The heavier the initial PBH mass, the more significant the
contribution to $\Delta N_{\rm eff}$ as heavier PBHs have longer lifetimes and
produce more gravitons. Furthermore, we observe that the higher the initial
spin of the PBH populations, the more significant the contribution to the
$\Delta N_{\rm eff}$ as the graviton production is enhanced. To visualize how
the non-instantaneous nature of the evaporation affects the results, we also
derived similar results where we account for the redshift of the evaporation
products only after the PBHs entirely evaporated (dashed lines). As one can
see, the results match in the case of Schwarzschild PBHs ($a_{\star}=0$),
whereas a discrepancy develops for increasing spins. For the case
$a_{\star}=0.99$, we observe that the effect of including redshift during the
PBH evaporation process reduces the contribution to $\Delta N_{\rm eff}$ by
almost an order of magnitude over the entire considered mass range of PBHs
when we compare our results (solid line) with the findings of Ref. Hooper  _et
al._ (2020) (dot-dashed line), Ref. Arbey  _et al._ (2021) (dotted line), and
Ref. Masina (2021) (dot-dot-dashed horizontal bar) which we refer to as Paper
A, Paper B and Paper C, respectively, in Fig. 3. In Paper B it was noted that
the dynamics of the evaporation are modified in the case of a Kerr black hole,
which affects the value of $\rho_{\rm R}(T_{\rm EV})$ at the end of
evaporation. However, it is assumed that the energy density of dark radiation
at evaporation $\rho_{\rm DR}(T_{\rm EV})$ is affected in the same proportion,
which led the authors of Paper B to claim that this does not affects the
calculation of $\Delta N_{\rm eff}$ since only the ratio of the two quantities
appears in Eq. (16). The authors of Paper A do not explicitly explain whether
their calculation of $\Delta N_{\rm eff}$ includes the effect of the dark-
radiation redshift, they obtain quantitatively similar results to Paper B,
which suggests that they used a similar approximation. Interestingly, Paper C
in the figure (horizontal dot-dot-dashed bar), does show a discrepancy with
the former papers in the case of Kerr PBHs. This reference considers the
redshift when integrating over the evaporation dynamics. However, Paper C also
assumes that the Universe is exactly matter-dominated until PBHs fully
evaporate. This does not fully account for the cosmological dynamics of the
evaporation and explains why the results of Paper C are somewhat intermediate
between our results and the results presented in Paper A and Paper B.

Unsurprisingly, our results which do not account for the redshift of the
evaporation products during the evaporation process, match the results of
Paper A and Paper B relatively well. This is easy to understand: the
production of gravitons via PBH evaporation is enhanced for larger PBH spin.
However, as we have seen in Section II, the spin of a Kerr PBH is depleted by
evaporation much faster than its mass. Gravitons, which are thus produced
copiously at the beginning of the evaporation, behave like radiation, while
the mass energy of the remaining PBHs behaves like matter until PBHs fully
evaporate mostly into SM particles. In the case of PBHs with large spin, the
time difference between the production of gravitons and SM particles is large.
In that case the energy density of gravitons redshifts significantly until the
final production of SM particles during which the production of gravitons is
negligible. This has the effect of depleting the value of $\rho_{\rm
DR}(T_{\rm EV})$ used in Eq. (16), which, in turn, explains the large
suppression of $\Delta N_{\rm eff}$ that we obtain for non-zero PBH spins.
This redshift effect can only be captured by solving both the mass and spin
loss rate of the PBHs while simultaneously tracking in time the evolution of
the Universe. We note that our results and those of Papers A and Papers B all
exhibit the characteristic bump at $M_{\rm{in}}=8\times 10^{7}\,$g
(corresponding to a reheating temperature of $\sim\,100$ MeV), which manifests
from the fact that the number of degrees of freedom is sensitive to the QCD
equation of state. In Fig. 3 we also show the current and future predicted
sensitivities to $\Delta N_{\rm eff}$ by the CMB-HD experiments on $\Delta
N_{\rm eff}\gtrsim 0.027$ at the $2\sigma$ C.L. Aiola  _et al._ (2022). For a
review of the sensitivities of other cosmic surveys including CORE Di
Valentino  _et al._ (2018), DESI Aghamousa  _et al._ (2016) and EUCLID
Laureijs  _et al._ (2011) see Ref. Lattanzi and Gerbino (2018). We observe
that including the effect of redshift bears a pessimistic result for current
or future experiments. Nonetheless, it is essential to realize that current
CMB results are not close to constraining early maximally spinning PBH
populations through the emission of gravitons.

Fig. 4:  Colour map showing PBH contribution to $\Delta N_{\rm eff}$ in the
$\beta^{\prime}$ and $M_{\rm BH}^{\rm in}$ plane. As in Fig. 3 we assume only
the Standard Model and massless graviton as well as monochromatic $M_{\rm
BH}^{\rm in}$ and $a_{\star}$ distributions. Here we show the results for
$a_{\star}=0.99$.

In Fig. 4, we show the contribution of hot gravitons to $\Delta N_{\rm eff}$
from the evaporation of Kerr PBHs according to our calculation described
above. The PBH distributions considered are monochromatic in $M_{\rm BH}^{\rm
in}$, which is varied along the x-axis, and we keep $a_{\star}=0.99$ fixed.
More realistic cases of mass and spin distributions will be investigated in
forthcoming works. We no longer assume that $\beta^{\prime}>\beta_{c}$, where
$\beta_{c}$ is the critical value of the initial PBH density from which PBHs
would eventually dominate the evolution of the Universe. We observe that, if
one assumes instantaneous evaporation (black dashed line), one would expect
CMB-HD experiments can achieve sensitivity below $\beta_{c}$. However, when
the effect of redshift during the non-instantaneous evaporation is considered
(brown dashed line), the projected sensitivity for CMB-HD is above the
$\beta_{c}$ except for when $M_{\rm BH}^{\rm in}\gtrsim 10^{7}\,\,{\rm g}$.
Finally, it can be noted that the production of gravitons from PBHs with spin
$a=0.7$ escapes the sensitivity of future experiments, as can also be seen
from Fig. 3.

###  III.2 BSM contributions to $\Delta N_{{\rm eff}}$

This section explores possible contributions to $\Delta N_{\rm eff}$ from
light particles, which are marginally more hypothetical than the graviton.
However, with the recent determination of the muon anomalous magnetic moment,
$a_{\mu}=(g-2)_{\mu}/2$, by the E989 experiment at Fermilab Abi  _et al._
(2021) and planned CERN NA64$\mu$ Gninenko (2018, 2019) and Fermilab $M^{3}$
Kahn  _et al._ (2018) experiments, the number of known light vectors may be
about to change.

As is well established, the emission of spin-1 particles is enhanced for
rotating BHs as they attempt to shed their angular momentum. However, this
enhancement is less pronounced than in the case of the graviton as they have a
lower spin. On the other hand, the absorption cross-section, $\sigma_{s_{i}}$
as seen in Eq. (1), calculated in the Schwarzchild metric for spin-1
particles, is much less suppressed than for spin-2 particles. This means that
a more significant proportion of the radiation will be produced at late times,
meaning they will not experience as much dilution as their spin-2
counterparts.

Fig. 5:  Same as Fig. 4, but for a massless vector with two different values
for BH spin $a_{\star}=0.99$ and $a_{\star}=0.7$ on the left and right panels
respectively.

In Fig. 5, we show the effect this new vector would have on $\Delta N_{\rm
eff}$ if there was a population of PBHs in the Early Universe. First, let us
focus on the left panel and compare our results for the contribution of hot
gravitons or dark vector bosons to $\Delta N_{\rm eff}$. Generally, the
picture looks less pessimistic. For both cases, taking the instantaneous
evaporation approximation overestimates the sensitivity projections for CMB-
HD, but to a much lesser extent. For $a_{\star}=0.99$, the projections when
accounting for non-instantaneous evaporation is much more constraining than
the tensor case in Fig. 4. Furthermore, as we see in the right panel of Fig.
5, CMB-HD will be sensitive to vector radiation from a population of black
holes with a more modest spin, $a_{\star}=0.7$. In fact above $M_{\rm BH}^{\rm
in}\sim 10^{8}\,\,{\rm g}$ the sensitivity we estimate is very similar for
both values of $a_{\star}$, where the lines go below $\beta_{c}$.

##  IV Phase space distribution of matter

If the emitted particles from PBHs are massive, stable and sufficiently inert,
they may constitute part or all of the dark matter required by our current
understanding of structure formation in the Universe. This mechanism of
gravitational production of dark matter does not rely on any other interaction
and is extremely efficient, particularly for heavy dark matter masses, up to
$m_{\rm DM}\sim M_{p}$, see Refs. Fujita  _et al._ (2014); Lennon  _et al._
(2018); Hooper  _et al._ (2019); Masina (2020); Baldes  _et al._ (2020);
Gondolo  _et al._ (2020); Bernal and Zapata (2021a, b); Auffinger  _et al._
(2021); Bernal and Zapata (2021c); Arbey  _et al._ (2021); Masina (2021);
Jyoti Das  _et al._ (2021); Cheek  _et al._ (2022a, b); Sandick  _et al._
(2021); Bernal  _et al._ (2022) for further details. If PBH evaporation is
solely responsible for producing massive non-interacting dark matter, then the
dark matter would cool with the expansion of the Universe and be consistent
with the cold dark matter picture. This may not necessarily be the case for
masses below $~{}10^{3}$ GeV, and if the dark matter was produced too late,
their free-streaming length might be incompatible with observations of small-
scale structures Bode  _et al._ (2001).

In this section, we improve on our previous estimates for the warm dark matter
constraints given in Ref. Cheek  _et al._ (2022a) by following the methodology
of Baldes  _et al._ (2020); Auffinger  _et al._ (2021); Masina (2021). The key
difference from our previous work is that in this paper, we track the complete
evolution of the dark matter phase-space distribution to obtain the input
parameters for the computational tool for calculating Anisotropies, CLASS
Lesgourgues (2011); Blas  _et al._ (2011); Lesgourgues and Tram (2011).
Furthermore, what separates this work from the literature Refs. Baldes  _et
al._ (2020); Auffinger  _et al._ (2021); Masina (2021), is that we have
numerically solved the Friedmann and Boltzmann equations and by using the
precise evolution of the plasma temperature $T(t)$ and scale factor $a(t)$.

Following Baldes  _et al._ (2020), one needs to determine the phase space of
the dark matter distribution, $f_{\rm DM}$ which is defined as,

$g_{\rm DM}f_{\rm DM}=\frac{{\rm d}n_{\rm DM}}{{\rm d}^{3}p},$ (17)

where $g_{\rm DM}$ is the number of degrees of freedom, $p$ is the three-
momentum, and $n_{\rm DM}$ is the dark matter number density. One provides
$f_{\rm DM}$ at the time of dark matter production to CLASS . In our case this
is $t_{\rm ev}$, which is easily defined for the monochromatic distributions.
Therefore, one can write

$f_{\rm DM}=\left.\frac{n_{\rm BH}\left(t_{\rm in}\right)}{g_{\rm
DM}}\left(\frac{a(t_{\rm in})}{a(t)}\right)^{3}\frac{1}{p^{2}}\frac{{\rm
d}\mathcal{N}_{\rm DM}}{{\rm d}p}\right|_{t=t_{\rm ev}},$ (18)

where $t_{\rm in}$ is the scale factor at formation, and we have used the fact
that the co-moving number density of BHs is constant up until they evaporate.
The particle emission rate per momentum is calculated by taking the emission
rate given in Eq. (1) and integrating over time,

$\frac{{\rm d}\mathcal{N}_{\rm DM}}{{\rm d}p}=\int_{0}^{\tau}{\rm
d}t^{\prime}\frac{a(\tau)}{a(t^{\prime})}\times\frac{{\rm
d}^{2}\mathcal{N}_{\rm DM}}{{\rm d}p^{\prime}{\rm
d}t^{\prime}}\left(p\frac{a(\tau)}{a(t^{\prime})},t^{\prime}\right)\,,$ (19)

one can see that the non-instantaneous nature of particle emission is taken
into account by including the ratio of scale factors at a given $t$. Since the
above determination of ${\rm d}\mathcal{N}_{\rm DM}/{\rm d}p$ is the norm in
the literature, we do not expect differences on the scale found in Section
III, however we note that, as pointed out in Ref. Cheek  _et al._ (2022a), the
relic density calculation is highly sensitive to the value of $a(t)$ and
$T_{\rm ev}$. In fact, $T_{\rm ev}$ also enters in the determination of the
matter power spectrum as outlined in Ref. Baldes  _et al._ (2020), when one
evaluates the time-independent NCDM temperature,

$\mathcal{T}_{\rm ncdm}=T_{\rm in}\frac{a(t_{\rm ev})}{T(t_{0})}=\frac{T_{\rm
in}}{T_{\rm ev}}\left(\frac{g_{s\star}(T_{0})}{g_{s\star}(T_{\rm
ev})}\right)^{1/3}\,,$ (20)

where $T_{\rm in}$ is the plasma temperature at the time of the PBH formation,
which we assume to take place during radiation domination222 We can simply
determine this initial temperature from the initial PBH mass in Eq. (12). ,
and $T(t_{0})$ is the temperature today.

Equipped with the inputs from Eqs. (19) and (20) one can determine the matter
power spectrum, $P(K)$, from CLASS, and compare to that of the cold dark
matter (CDM), matter power spectrum by way of the transfer function,

$P(k)=P_{\rm CDM}(k)T^{2}(k)\,,$ (21)

where $k$ is the wavenumber. Barring inconclusive or short lived anomalies,
the CDM paradigm is largely supported by astrophysical observation Fairbairn
(2022) meaning that $T^{2}(k)\approx 1$ in the regions where experimental
observations have been made. Naturally, observing structures at greater $k$ is
increasingly difficult because such small galaxies do not support star
formation. Measurements from Lyman-$\alpha$ are the probes of $P(k)$ at
highest $k$ Viel  _et al._ (2013); Palanque-Delabrouille  _et al._ (2020), see
Refs. Varma  _et al._ (2020); He  _et al._ (2022) for work on proposed studies
at higher $k$. Following Refs. Bode  _et al._ (2001); Viel  _et al._ (2005);
Baldes  _et al._ (2020) we perform parameter fits in the parameterisation

$T(k)=(1+(\alpha k)^{2\mu})^{-5/\mu}\,,$ (22)

where $\mu$ is dimensionless exponent which is fixed to $\mu=1.12$ as in Ref.
Viel  _et al._ (2005) and $\alpha$ is the breaking scale, which we take to be
saturated at $\alpha=1.3\times 10^{-2}\,{\rm Mpc}\,h^{-1}$ Viel  _et al._
(2005); Palanque-Delabrouille  _et al._ (2020); Garzilli  _et al._ (2021);
Baldes  _et al._ (2020).

The above methodology has been performed previously in the literature. For
example, in Refs. Baldes  _et al._ (2020); Auffinger  _et al._ (2021) these
constraints were calculated for Schwarzschild black holes where Ref. Auffinger
_et al._ (2021) uses emission rates shown in Eq. (1) from the BlackHawk code.
Ref. Masina (2021) estimated the effects Kerr distributions would have on the
warm dark matter constraints by using a weighted integral of the emission
rate. The primary difference in our work is that we use the solution of the
coupled Boltzmann equations such that we have fully tracked the evolution of
the scale factor and plasma temperature. For scalar, fermion and vector
particles we use emission functions calculated following the procedure in
Refs. Chandrasekhar and Detweiler (1975); Chandrasekhar (1976); Chandrasekhar
and Detweiler (1977) to compute the cross-sections $\sigma_{s_{i}}^{lm}$,
where we find good agreement with BlackHawk Arbey and Auffinger (2019). For
the spin-2 case Chandrasekhar and Detweiler (1976) we use the greybody factors
from BlackHawk Arbey and Auffinger (2019) which we checked to be consistent
with Ref. Dong  _et al._ (2016).

Using FRISBHEE, for a given $m_{\rm DM}$, $a_{\star}$ and intrinsic particle
spin, $s_{\rm DM}$, we determine the values of $\beta^{\prime}$ and $M_{\rm
BH}^{\rm in}$ such that the dark matter relic abundance matches that of the
Planck collaboration’s measured value, $\Omega_{\rm DM}\,h^{2}=0.12$ Aghanim
_et al._ (2020). We then take those parameter values to calculate $f_{\rm DM}$
and $\mathcal{T}_{\rm ncdm}$ using Eq. (18) and Eq. (20).

Fig. 6: Figure showing the $\beta^{\prime}$ and $M_{\rm BH}^{\rm in}$ values
required to fully produce the correct relic abundance (solid black lines) for
a range of dark matter masses. This is for the case where dark matter is a
Fermion and all black holes are Schwarzschild. The black $+$ symbols indicate
the precise point on the line where the transfer function, $T(k)$, becomes
inconsistent with Lyman-$\alpha$ constraints, i.e. the $\alpha$ parameter of
Eq. (22) is $\geq 1.3\times 10^{-2}\,{\rm Mpc}\,h^{-1}$. The red dashed line
corresponds to the mean value found for $\eta$ as per Eq. (23), where the
values that violate the inequality are shaded red. The dotted line corresponds
to the result reported in Ref. Baldes  _et al._ (2020). The purple shaded
regions signify the inflationary (low $M_{\rm BH}^{\rm in}$) Akrami  _et al._
(2020) and Big Bang Nucleosynthesis (high $M_{\rm BH}^{\rm in}$) constraints
on primordial black hole abundances Carr  _et al._ (2010); Hasegawa  _et al._
(2019); Carr  _et al._ (2020); Keith  _et al._ (2020).

Fig. 6 illustrates our procedure, in which we show the $\beta^{\prime}$ and
$M_{\rm BH}^{\rm in}$ values needed to produce the correct relic abundance
through evaporation only for a range of $m_{\rm DM}$ values. Above these
lines, dark matter would be overproduced and below it would be underproduced.
Along each line, a black + symbol appears, which signifies the region of
parameter space where the bound from Lyman-$\alpha$ is saturated, such that
points with greater $M_{\rm BH}^{\rm in}$ will be in tension with the observed
matter power spectrum.

In agreement with the literature before us, we find that for monochromatic
distributions, the warm dark matter constraint, applied to dark matter with
masses, $m_{\rm DM}\leq 10^{3}\,$ GeV, can be be well approximated by a simple
relation,

$\beta^{\prime}\leq\eta\left(\frac{M_{p}}{M_{\rm BH}^{\rm in}}\right)\,,$ (23)

where $\eta$ is a dimensionless parameter that varies depending on
$a_{\star}$. In Refs. Baldes  _et al._ (2020); Auffinger  _et al._ (2021) the
limit is presented as a multiple of $\beta_{c}$, the $\beta^{\prime}$ value
required such that PBHs dominate at some point. Since there exists multiple
definitions of $\beta_{c}$ in the literature Gondolo  _et al._ (2020); Baldes
_et al._ (2020) and changes depending on $a_{\star}$, we believe that
presenting our results in the form of Eq. (23) is easier to interpret. In
order to determine a computational error on the value of $\eta$, we sample
multiple $m_{\rm DM}$ values for a given $a_{\star}$ and calculate the
standard error of the mean from the set of $\eta$ values.

Returning to Fig. 6, we observe that the black + symbols, signifying where the
warm dark matter constraints become relevant, align for different $m_{\rm DM}$
values such that one can write Eq. (23). For Schwarzschild PBHs, fermion dark
matter case shown in this figure, we find $\eta=(6.26\pm 0.08)\times 10^{-4}$,
this is the red dashed line in the figure, which borders the red shaded region
representing where the full dark matter production through PBH evaporation is
inconsistent with the matter power spectrum. This allows one to interpolate
between the $m_{\rm DM}$ values shown in the figure. The dotted red line shows
the warm dark matter constraint found in Refs Baldes  _et al._ (2020);
Auffinger  _et al._ (2021). This level of disagreement, we believe, is
acceptable because of the potential for quite sizeable differences in the
ratios of scale factors in Eq. (18) and temperature in Eq. (20).

Fig. 7: The variation of $\eta$ for different $a_{\star}$ values of the pbh
distribution. Different colours correspond to different dark matter spins. The
error bars indicate the computational variation of this value by way of the
standard error of the mean calculated from the samples of different $m_{\rm
DM}$ values.

In Fig. 7, we show our determination of $\eta$ and its standard error for
scalar, fermion, vector and tensor dark matter particles with green, purple,
orange and blue lines and error bars, respectively. One can see that there is
very little variation of $\eta$ for each value of $a_{\star}$ thus validating
the expression in Eq. (23). However, as $a_{\star}$ varies, we see a
substantial variation in $\eta$, where the effect is most pronounced for the
tensor dark matter case. For the scalar, fermion and vector case, the warm
dark matter constraint is hardly changing with $a_{\star}$ except for the very
maximal case. This is qualitatively in agreement with Ref. Masina (2021),
where the author also concluded that changes in $a_{\star}$ mainly altered the
warm dark matter constraints for the tensor particle. Our analysis, which
solves both the Friedmann and Boltzmann equations, places this conclusion on
firmer ground.

When one compares the $a_{\star}=0$ case in Fig.7, we can make a comparison
with warm dark matter constraints reported in Ref. Auffinger  _et al._ (2021).
Similarly to the fermion case as mentioned above, we observe a slight
deviation with our calculation. However the general pattern is consistent, as
with Ref Auffinger  _et al._ (2021) the scalar and fermion constraints are
found to be very similar and the vector and tensor limits are less aggressive.
More precisely the order of magnitude difference we see between scalar and
tensor dark matter constraints is also reported in Ref. Auffinger  _et al._
(2021). This difference is largely attributable to the relative red-shift felt
from particles with different intrinsic spin, as seen in section III.

Throughout the analysis we additionally performed 2-dimensional fit for the
parameterisation of Eq.(22) were we fit for both $\alpha$ and $\mu$. We find
in general quite good agreement with $\mu=1.12$, with only slight variations.
Furthermore with the 1-dimensional fit, the $T(k)$ always matches well with
simulation results. We leave a more sophisticated study for future work. This
is particularly relevant for the situation where PBHs produce only a sub-
component of dark matter, where Eq.(22) is often not valid Baur  _et al._
(2017); Auffinger  _et al._ (2021).

##  V Conclusions

In previous works, it was argued that the production of dark vector or tensor
bosons (such as dark photons or hot gravitons) is enhanced for dark particles
with spin larger than one when the evaporating black holes approach
extremality. We have shown in this paper that this enhancement is mitigated
when one thoroughly considers the dynamics of evaporation. While primordial
black holes evaporate, their spin is generically dissipated first, whereas
their mass is depleted in a second phase. Since the production of particles
with spin larger than one is enhanced when primordial black holes are
rotating, the largest contribution of the evaporation to dark radiation is
emitted during the first phase of the evaporation. Consequently, these dark
particles are subject to more significant redshifting, which can only be
accounted for by solving the complete set of equations tracking the PBH’s
evolution and the Universe’s. On the contrary, the second phase of the
evaporation dilutes this primary contribution, leading to a subsequent
reduction of $\Delta N_{\rm eff}$. Considering this dilution, we showed that
the maximum contribution of PBH evaporation to $\Delta N_{\rm eff}$ is reduced
by a factor of $\mathcal{O}(10)$ for maximally rotating black holes. In the
case of lower spin PBHs, in the limit $a_{\star}\to 0$, our results align with
those in the literature.

Using FRISBHEE we could also reconstruct the full phase-space distribution of
light dark-matter particles when they are produced through PBH evaporation and
use the code CLASS to determine the corresponding matter power spectrum and
consequently the warm dark matter constraints on dark matter produced from
black hole evaporation. Our present study builds on previous works by
introducing, for the first time, the fully evolved scale factor, $a(t)$ and
plasma temperature $T$. We also provide a useful parametrization of warm dark
matter limit for scalar, fermion, vector and tensor dark matter candidates,
respectively. We show how this parametrization varies over a monochromatic
Kerr PBH distribution of different $a_{\star}$ values.

Our results have been obtained assuming monochromatic distributions for mass
and spin of the PBH. Thus, one may wonder how these conclusions would be
affected if we consider more realistic spectra. We are currently investigating
scenarios with extended distributions, and we plan to report our results in
the near future.

## Acknowledgments

The authors would like to thank Julia Stadler and Deanna Hooper for their help
with using CLASS. AC is supported by the grant “AstroCeNT: Particle
Astrophysics Science and Technology Centre” carried out within the
International Research Agendas programme of the Foundation for Polish Science
financed by the European Union under the European Regional Development Fund.
The work of LH is funded by the UK Science and Technology Facilities Council
(STFC) under grant ST/P001246/1. This work used the DiRAC@Durham facility
managed by the Institute for Computational Cosmology on behalf of the STFC
DiRAC HPC Facility (www.dirac.ac.uk). The equipment was funded by BEIS capital
funding via STFC capital grants ST/P002293/1, ST/R002371/1 and ST/S002502/1;
Durham University and STFC operations grant ST/R000832/1. DiRAC is part of the
National e-Infrastructure. This work has made use of the Hamilton HPC Service
of Durham University.

## References

  * Abbott  _et al._ (2016) B. P. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. Lett. 116, 061102 (2016), arXiv:1602.03837 [gr-qc] .
  * Abbott  _et al._ (2017) B. P. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. Lett. 119, 161101 (2017), arXiv:1710.05832 [gr-qc] .
  * Auclair  _et al._ (2022) P. Auclair _et al._ (LISA Cosmology Working Group), (2022), arXiv:2204.05434 [astro-ph.CO] .
  * Jeong and Takahashi (2013) K. S. Jeong and F. Takahashi, Phys. Lett. B 725, 134 (2013), arXiv:1305.6521 [hep-ph] .
  * Buen-Abad  _et al._ (2015) M. A. Buen-Abad, G. Marques-Tavares, and M. Schmaltz, Phys. Rev. D 92, 023531 (2015), arXiv:1505.03542 [hep-ph] .
  * Lesgourgues  _et al._ (2016) J. Lesgourgues, G. Marques-Tavares, and M. Schmaltz, JCAP 02, 037 (2016), arXiv:1507.04351 [astro-ph.CO] .
  * Buen-Abad  _et al._ (2018) M. A. Buen-Abad, M. Schmaltz, J. Lesgourgues, and T. Brinckmann, JCAP 01, 008 (2018), arXiv:1708.09406 [astro-ph.CO] .
  * Blinov and Marques-Tavares (2020) N. Blinov and G. Marques-Tavares, JCAP 09, 029 (2020), arXiv:2003.08387 [astro-ph.CO] .
  * Chacko  _et al._ (2004) Z. Chacko, L. J. Hall, T. Okui, and S. J. Oliver, Phys. Rev. D 70, 085008 (2004), arXiv:hep-ph/0312267 .
  * Chacko  _et al._ (2005) Z. Chacko, L. J. Hall, S. J. Oliver, and M. Perelstein, Phys. Rev. Lett. 94, 111801 (2005), arXiv:hep-ph/0405067 .
  * Berlin and Blinov (2018) A. Berlin and N. Blinov, Phys. Rev. Lett. 120, 021801 (2018), arXiv:1706.07046 [hep-ph] .
  * Zel’dovich and Novikov (1967) Y. B. Zel’dovich and I. D. Novikov, Soviet Astron. AJ (Engl. Transl. ), 10, 602 (1967).
  * Hawking (1971) S. Hawking, Mon. Not. Roy. Astron. Soc. 152, 75 (1971).
  * Carr and Hawking (1974) B. J. Carr and S. W. Hawking, Mon. Not. Roy. Astron. Soc. 168, 399 (1974).
  * Inomata  _et al._ (2020) K. Inomata, M. Kawasaki, K. Mukaida, T. Terada, and T. T. Yanagida, Phys. Rev. D 101, 123533 (2020), arXiv:2003.10455 [astro-ph.CO] .
  * Domènech  _et al._ (2021) G. Domènech, C. Lin, and M. Sasaki, JCAP 04, 062 (2021), [Erratum: JCAP 11, E01 (2021)], arXiv:2012.08151 [gr-qc] .
  * Hawking (1974) S. W. Hawking, Nature 248, 30 (1974).
  * Hawking (1975) S. W. Hawking,  _Euclidean quantum gravity_ , Commun. Math. Phys. 43, 199 (1975), [,167(1975)].
  * Dong  _et al._ (2016) R. Dong, W. H. Kinney, and D. Stojkovic, JCAP 10, 034 (2016), arXiv:1511.05642 [astro-ph.CO] .
  * Page (1976a) D. N. Page, Phys. Rev. D 13, 198 (1976a).
  * Page (1976b) D. N. Page, Phys. Rev. D 14, 3260 (1976b).
  * Page (1977) D. N. Page, Phys. Rev. D 16, 2402 (1977).
  * Hooper  _et al._ (2019) D. Hooper, G. Krnjaic, and S. D. McDermott, JHEP 08, 001 (2019), arXiv:1905.01301 [hep-ph] .
  * Hooper  _et al._ (2020) D. Hooper, G. Krnjaic, J. March-Russell, S. D. McDermott, and R. Petrossian-Byrne, (2020), arXiv:2004.00618 [astro-ph.CO] .
  * Masina (2020) I. Masina, Eur. Phys. J. Plus 135, 552 (2020), arXiv:2004.04740 [hep-ph] .
  * Masina (2021) I. Masina, Grav. Cosmol. 27, 315 (2021), arXiv:2103.13825 [gr-qc] .
  * Arbey  _et al._ (2021) A. Arbey, J. Auffinger, P. Sandick, B. Shams Es Haghi, and K. Sinha, (2021), arXiv:2104.04051 [astro-ph.CO] .
  * Aiola  _et al._ (2022) S. Aiola _et al._ (CMB-HD), (2022), arXiv:2203.05728 [astro-ph.CO] .
  * Fujita  _et al._ (2014) T. Fujita, M. Kawasaki, K. Harigaya, and R. Matsuda, Phys. Rev. D89, 103501 (2014), arXiv:1401.1909 [astro-ph.CO] .
  * Lennon  _et al._ (2018) O. Lennon, J. March-Russell, R. Petrossian-Byrne, and H. Tillim, JCAP 1804, 009 (2018), arXiv:1712.07664 [hep-ph] .
  * Morrison  _et al._ (2019) L. Morrison, S. Profumo, and Y. Yu, JCAP 1905, 005 (2019), arXiv:1812.10606 [astro-ph.CO] .
  * Baldes  _et al._ (2020) I. Baldes, Q. Decant, D. C. Hooper, and L. Lopez-Honorez, JCAP 08, 045 (2020), arXiv:2004.14773 [astro-ph.CO] .
  * Gondolo  _et al._ (2020) P. Gondolo, P. Sandick, and B. Shams Es Haghi, Phys. Rev. D 102, 095018 (2020), arXiv:2009.02424 [hep-ph] .
  * Bernal and Zapata (2021a) N. Bernal and O. Zapata, JCAP 03, 007 (2021a), arXiv:2010.09725 [hep-ph] .
  * Bernal and Zapata (2021b) N. Bernal and O. Zapata, Phys. Lett. B 815, 136129 (2021b), arXiv:2011.02510 [hep-ph] .
  * Auffinger  _et al._ (2021) J. Auffinger, I. Masina, and G. Orlando, Eur. Phys. J. Plus 136, 261 (2021), arXiv:2012.09867 [hep-ph] .
  * Bernal and Zapata (2021c) N. Bernal and O. Zapata, JCAP 03, 015 (2021c), arXiv:2011.12306 [astro-ph.CO] .
  * Jyoti Das  _et al._ (2021) S. Jyoti Das, D. Mahanta, and D. Borah, (2021), arXiv:2104.14496 [hep-ph] .
  * Cheek  _et al._ (2022a) A. Cheek, L. Heurtier, Y. F. Perez-Gonzalez, and J. Turner, Phys. Rev. D 105, 015022 (2022a), arXiv:2107.00013 [hep-ph] .
  * Cheek  _et al._ (2022b) A. Cheek, L. Heurtier, Y. F. Perez-Gonzalez, and J. Turner, Phys. Rev. D 105, 015023 (2022b), arXiv:2107.00016 [hep-ph] .
  * Sandick  _et al._ (2021) P. Sandick, B. S. Es Haghi, and K. Sinha, Phys. Rev. D 104, 083523 (2021), arXiv:2108.08329 [astro-ph.CO] .
  * Bernal  _et al._ (2022) N. Bernal, Y. F. Perez-Gonzalez, and Y. Xu, (2022), arXiv:2205.11522 [hep-ph] .
  * Lesgourgues (2011) J. Lesgourgues, arXiv e-prints , arXiv:1104.2932 (2011), arXiv:1104.2932 [astro-ph.IM] .
  * Blas  _et al._ (2011) D. Blas, J. Lesgourgues, and T. Tram, JCAP 2011, 034 (2011), arXiv:1104.2933 [astro-ph.CO] .
  * Lesgourgues and Tram (2011) J. Lesgourgues and T. Tram, JCAP 2011, 032 (2011), arXiv:1104.2935 [astro-ph.CO] .
  * Chandrasekhar and Detweiler (1975) S. Chandrasekhar and S. L. Detweiler, Proc. Roy. Soc. Lond. A 345, 145 (1975).
  * Chandrasekhar and Detweiler (1976) S. Chandrasekhar and S. L. Detweiler, Proc. Roy. Soc. Lond. A 350, 165 (1976).
  * Chandrasekhar and Detweiler (1977) S. Chandrasekhar and S. L. Detweiler, Proc. Roy. Soc. Lond. A 352, 325 (1977).
  * MacGibbon and Webber (1990) J. H. MacGibbon and B. R. Webber, Phys. Rev. D 41, 3052 (1990).
  * MacGibbon (1991) J. H. MacGibbon, Phys. Rev. D 44, 376 (1991).
  * Fishbach  _et al._ (2017) M. Fishbach, D. E. Holz, and B. Farr, Astrophys. J. Lett. 840, L24 (2017), arXiv:1703.06869 [astro-ph.HE] .
  * Mangano  _et al._ (2002) G. Mangano, G. Miele, S. Pastor, and M. Peloso, Phys. Lett. B 534, 8 (2002), arXiv:astro-ph/0111408 .
  * de Salas and Pastor (2016) P. F. de Salas and S. Pastor, JCAP 07, 051 (2016), arXiv:1606.06986 [hep-ph] .
  * Carr  _et al._ (2020) B. Carr, K. Kohri, Y. Sendouda, and J. Yokoyama, (2020), arXiv:2002.12778 [astro-ph.CO] .
  * Gutiérrez  _et al._ (2018) E. M. Gutiérrez, F. L. Vieyro, and G. E. Romero, Mon. Not. Roy. Astron. Soc. 473, 5385 (2018), arXiv:1710.03061 [astro-ph.CO] .
  * Di Valentino  _et al._ (2018) E. Di Valentino _et al._ (CORE), JCAP 04, 017 (2018), arXiv:1612.00021 [astro-ph.CO] .
  * Aghamousa  _et al._ (2016) A. Aghamousa _et al._ (DESI), (2016), arXiv:1611.00036 [astro-ph.IM] .
  * Laureijs  _et al._ (2011) R. Laureijs _et al._ (EUCLID), (2011), arXiv:1110.3193 [astro-ph.CO] .
  * Lattanzi and Gerbino (2018) M. Lattanzi and M. Gerbino, Front. in Phys. 5, 70 (2018), arXiv:1712.07109 [astro-ph.CO] .
  * Abi  _et al._ (2021) B. Abi _et al._ (Muon g-2), Phys. Rev. Lett. 126, 141801 (2021), arXiv:2104.03281 [hep-ex] .
  * Gninenko (2018) S. Gninenko (NA64 Collaboration),  _Addendum to the Proposal P348: Search for dark sector particles weakly coupled to muon with NA64 $\mu$_, Tech. Rep. CERN-SPSC-2018-024. SPSC-P-348-ADD-3 (CERN, Geneva, 2018).
  * Gninenko (2019) S. Gninenko (NA64 Collaboration),  _Proposal for an experiment to search for dark sector particles weakly coupled to muon at the SPS_ , Tech. Rep. CERN-SPSC-2019-002. SPSC-P-359 (CERN, Geneva, 2019).
  * Kahn  _et al._ (2018) Y. Kahn, G. Krnjaic, N. Tran, and A. Whitbeck, JHEP 09, 153 (2018), arXiv:1804.03144 [hep-ph] .
  * Bode  _et al._ (2001) P. Bode, J. P. Ostriker, and N. Turok, Astrophys. J. 556, 93 (2001), arXiv:astro-ph/0010389 .
  * Fairbairn (2022) M. Fairbairn, Symmetry 14, 812 (2022).
  * Viel  _et al._ (2013) M. Viel, G. D. Becker, J. S. Bolton, and M. G. Haehnelt, Phys. Rev. D 88, 043502 (2013).
  * Palanque-Delabrouille  _et al._ (2020) N. Palanque-Delabrouille, C. Yèche, N. Schöneberg, J. Lesgourgues, M. Walther, S. Chabanier, and E. Armengaud, JCAP 04, 038 (2020), arXiv:1911.09073 [astro-ph.CO] .
  * Varma  _et al._ (2020) S. Varma, M. Fairbairn, and J. Figueroa, (2020), arXiv:2005.05353 [astro-ph.CO] .
  * He  _et al._ (2022) Q. He, A. Robertson, J. Nightingale, S. Cole, C. S. Frenk, R. Massey, A. Amvrosiadis, R. Li, X. Cao, and A. Etherington, Mon. Not. Roy. Astron. Soc. 511, 3046 (2022), arXiv:2010.13221 [astro-ph.CO] .
  * Viel  _et al._ (2005) M. Viel, J. Lesgourgues, M. G. Haehnelt, S. Matarrese, and A. Riotto, Phys. Rev. D 71, 063534 (2005), arXiv:astro-ph/0501562 .
  * Garzilli  _et al._ (2021) A. Garzilli, A. Magalich, O. Ruchayskiy, and A. Boyarsky, Mon. Not. Roy. Astron. Soc. 502, 2356 (2021), arXiv:1912.09397 [astro-ph.CO] .
  * Chandrasekhar (1976) S. Chandrasekhar, Proc. Roy. Soc. Lond. A 348, 39 (1976).
  * Arbey and Auffinger (2019) A. Arbey and J. Auffinger, Eur. Phys. J. C 79, 693 (2019), arXiv:1905.04268 [gr-qc] .
  * Aghanim  _et al._ (2020) N. Aghanim _et al._ (Planck), Astron. Astrophys. 641, A6 (2020), [Erratum: Astron.Astrophys. 652, C4 (2021)], arXiv:1807.06209 [astro-ph.CO] .
  * Akrami  _et al._ (2020) Y. Akrami _et al._ (Planck), Astron. Astrophys. 641, A10 (2020), arXiv:1807.06211 [astro-ph.CO] .
  * Carr  _et al._ (2010) B. J. Carr, K. Kohri, Y. Sendouda, and J. Yokoyama, Phys. Rev. D 81, 104019 (2010), arXiv:0912.5297 [astro-ph.CO] .
  * Hasegawa  _et al._ (2019) T. Hasegawa, N. Hiroshima, K. Kohri, R. S. L. Hansen, T. Tram, and S. Hannestad, JCAP 12, 012 (2019), arXiv:1908.10189 [hep-ph] .
  * Keith  _et al._ (2020) C. Keith, D. Hooper, N. Blinov, and S. D. McDermott, (2020), arXiv:2006.03608 [astro-ph.CO] .
  * Baur  _et al._ (2017) J. Baur, N. Palanque-Delabrouille, C. Yeche, A. Boyarsky, O. Ruchayskiy, E. Armengaud, and J. Lesgourgues, JCAP 12, 013 (2017), arXiv:1706.03118 [astro-ph.CO] .