Evaporation of Primordial Black Holes in the Early Universe: Mass and Spin Distributions

Many cosmological phenomena lead to the production of primordial black holes in the early Universe. These phenomena often create a population of black holes with extended mass and spin distributions. As these black holes evaporate via Hawking radiation, they can modify various cosmological observables, lead to the production of dark matter, modify the number of effective relativistic degrees of freedom, $N_{\rm eff}$, source a stochastic gravitational wave background and alter the dynamics of baryogenesis. We consider the Hawking evaporation of primordial black holes that feature non-trivial mass and spin distributions in the early Universe. We demonstrate that the shape of such a distribution can strongly affect most of the aforementioned cosmological observables. We outline the numerical machinery we use to undertake this task. We also release a new version of FRISBHEE that handles the evaporation of primordial black holes with an arbitrary mass and spin distribution throughout cosmic history.


I. INTRODUCTION
The birth of gravitational wave astronomy [1,2] has produced a flurry of interest in Primordial Black Holes (PBHs) [3]. Unlike astrophysical black holes, which result from stellar collapse, PBHs formed in the early Universe when large over-densities collapse under their Schwarzchild radius. If proven to exist, the implications for our understanding of the post-inflationary Universe would be tremendous. This is because the primordial power spectrum implied by the cosmic microwave background (CMB) is insufficient to produce them. Evidence of PBHs would represent a concrete way to move beyond the current model of inflationary cosmology [4]. Measuring the PBH distribution would provide a direct view into their production, the early Universe [5][6][7] and inflation [8][9][10][11][12][13][14][15][16][17][18][19][20][21][22][23]. Even after their initial production, BHs can accrete the Standard Model (SM) plasma, merge with their peers, acquiring rotational momentum [24,25].
PBHs can therefore feature various mass and spin distributions depending on the cosmological scenario considered. Whereas PBHs with masses larger than 10 15 g are stable on cosmological time scales, lighter black holes may have evaporated via Hawking radiation before the current epoch [26,27]. The effect of this evaporation was studied in many different contexts (see e.g. Ref. [28] for a recent review). PBHs with intermediate masses in the range 10 8 − 10 12 g are known to evaporate after Big-Bang Nucleosynthesis (BBN) [29] and, therefore, cannot constitute a sizeable fraction of the Universe's energy density since, through Hawking radiation, they would change the neutron-to-proton ratio at the onset of BBN and therefore, the abundance of light elements that are measured today with excellent accuracy. However, a sizeable abundance of lighter PBHs may have formed in cosmic history without affecting the post-BBN era. Interestingly, such PBHs may even have dominated the energy density of the Universe before they evaporated, leading to a phase of early matter domination (EMD) and a subsequent reheating of the Universe. In both cases, several works have recently studied the imprints of such PBHs Hawking evaporation on particle and astrophysical data by studying its effect on the dark matter (DM) relic density and phase space distribution [30][31][32][33][34][35][36][37][38][39][40][41][42][43][44][45][46][47], the effective number of relativistic degrees of freedom, N eff , at the time of CMB emission [25,32,33,43,[48][49][50], the dynamics of Baryo/Leptogenesis [34,[51][52][53][54][55][56][57][58][59][60][61][62], the hydrogen 21-cm line [63], the production of gravitational waves [64][65][66], and the electroweak vacuum stability [67][68][69]. In most of these studies, the distribution of PBHs considered was monochromatic, either in mass, or in spin, see, however, Refs. [70,71] regarding the evaporation of a power-law mass distribution in the early Universe, Ref. [72] for a study of the impact of mass distributions on limits of PBH as DM, and Ref. [73] for a first attempt in considering baryogenesis in the context of extended distributions. In this paper, we consider the detailed evaporation of PBH populations that are not monochromatic, but instead spread over extended distributions in mass and/or in spin, tracking the evolution of the Universe carefully.
The paper is organised as follows: In Sec. II A we review the most well-studied mass and spin distributions that have been reported in the literature, and that we will be considering throughout this work. We then review in Sec. III the dynamics of the spin, a ⋆ , and mass, M , of a Kerr black hole during Hawking evaporation, and we extend this description in Sec. IV to the evaporation of an extended distribution in the plane (M, a ⋆ ). Details arXiv:2212.03878v2 [hep-ph] 28 Jul 2023 regarding the calculations are presented in Sec. IV and expanded in Appendix A, B and C, which are used in our numerical simulations. In the remainder of the paper, we consider the various observables that can be affected by non-monochromatic PBH distributions. Specifically, in Sec. VI A we explore the effect of extended mass distributions on the dynamics of the Universe when PBH that dominate the energy density evaporate and reheat the SM sector and follow in Sec. VI B where we study similar effects on the production of DM and explore how the DM relic density is affected by the extension of the distribution in various examples. In particular, we explore how the constraints on warm DM from measurements of the Lyman-α forest are accordingly affected. In Sec. VI C we derive constraints on PBHs that have an extended distribution in mass and spin, by looking at their contribution to ∆N eff . In Sec. VI D we comment on the possible imprints that extended distributions of PBHs could leave in the spectrum of primordial gravitational wave that could be observed in the future, and we finally conclude in Sec. VII.
Throughout, we use natural units where ℏ = c = k B = 1, and we define the Planck mass to be M p = 1/ √ G, with G the gravitational constant.

II. EXTENDED PBH DISTRIBUTIONS
Depending on their formation mechanism, the PBHs may form with a mass [74] and spin distribution. In this section, we review some well-motivated PBH mass and spin distributions that we consider throughout this paper. Further, we assume that the mass and spin distributions form at the same time in the evolution of the Universe and their distributions can be convoluted. We will discuss these assumptions later on.

A. Mass Distributions
In what follows, we will, for each of the mass distributions considered, denote by M c the value of the mass at which the corresponding mass fraction M × f PBH (M ) peaks, and by σ its width in logarithmic space.

a. Log-Normal (LN)
The production of PBHs from inflation usually requires the existence of a short period of ultra-slow-roll that produces a peak in the primordial power spectrum of scalar curvature perturbations [22,[75][76][77]. Generically, such peak is known to produce a log-normal mass function [78] and this has been numerically and analytically verified for slow-roll inflation [79,80]. The corresponding PBH mass distribution has the following form where M c is the initial peak of the distribution and σ is the width of the distribution. In the left most plot of Fig. 1, we show M × f PBH (M ) as a function of M for a central initial mass of PBH, M c = 10 6 g with varying values of the width σ.

b. Power-Law (PL)
Another possible formation mechanism of PBHs is the case where a large scale-invariant power spectrum of primordial perturbations collapses in a Universe that is dominated by a perfect fluid with constant equationof-state parameter, w. In that case, the distribution of PBHs takes the form [81] where the exponent α is given by and the mass range [M c , M c 10 σ ] depends on the domain of frequencies over which this scale-invariant power spectrum was formed. Physical situations in which the Universe is not inflating anymore typically correspond to values of w in the range −1/3 < w ⩽ 1, and thus to a scaling exponent In the central plot of Fig. 1, we show the power law mass distribution for varying α values with M c (M c 10 σ ) taking the value 10 (10 6 ) g.

c. Critical Collapse (CC)
The application of critical scaling to gravitational collapse is thought to describe the process of PBH formation from primordial fluctuations in a rigorous way [82][83][84][85]. Traditionally, overdensities were assumed to produce PBHs that had the same mass as their horizon mass. Instead, there is an upper cut-off at around the horizon mass but then a tail in the distribution for lower masses, this is a fairly generic finding over many inflationary models [85]. The resulting PBH distribution can attain a range of masses with the following form: where M c is the initial peak of the distribution. In the right most plot of Fig. 1, the brown line indicates a typical mass distribution from such a PBH formation mechanism.

d. Metric Preheating (MP)
Generically, PBHs are expected to form from the collapse of primordial perturbations that may form during or after inflation. In realistic particle physics models, inflation is usually followed by a phase of matter domination where the energy density of the Universe is dominated by the coherent oscillations of the inflaton field. In Ref. [86][87][88], it was noted that during that time, perturbations that were generated during the late inflationary era can get resonantly amplified and collapse into black holes before the Universe is reheated. Depending on the reheating temperature, the PBH mass fraction can peak at different masses. In order to obtain such a distribution, as shown by the green line in the rightmost plot of Fig. 1, one should in principle trace the collapse of the oscillating inflaton modes into PBHs numerically. In practice, such a production mechanism leads to an energy fraction of PBHs that is close to one and the PBHs formed dominate the energy density of the Universe quickly. For simplicity, we will consider in what follows the distribution exhibited in the Appendix of Ref. [88]. This distribution (let us denote it by f AV PBH ) shows a maximum around M AV c ∼ 10 5.6 g. Later, to extrapolate these results and explore different mass ranges, we assumed that the overall shape of the distribution does not change for different values of the reheating temperature, and we simply translate the distribution in log-space as follows:

B. Spin Distributions
PBHs formed from the collapse of primordial perturbations re-enetering the Hubble horizon are typically created without angular momentum, but may acquire some spin distribution via mergers or during phases of early matter domination [89][90][91][92][93]. In particular, if the rate of PBH binary capture is significant compared to the Hubble rate, the in-spiral phase may end before light PBHs start evaporating. After they merge repeatedly, the PBH spin distribution is expected to stabilise to a universal distribution which does not depend anymore on the PBH masses [94]. This spin distribution peaks strongly at ⟨a ⋆ ⟩ ∼ 0.7, with few BHs with a ⋆ < 0.4. Another spin distribution involving hierarchical mergers was proposed in Ref. [95] based on LIGO/VIRGO data regarding mergers in the Milky Way. As a matter of fact, this distribution peaks as well at a ⋆ ∼ 0.7, and in both cases, the spin distribution is almost entirely independent of the masses or the initial spin distribution of the merging binaries [94,95]. Although the black holes we consider are much lighter, as compared to the ones used in these different studies, it was shown that two non-rotating BHs with similar masses acquire a spin of order a ⋆ ≈ 0.69 after merging [96]. This universality suggests that the full distribution of PBHs before they start evaporating can be written as a product of the form In what follows, we will assume that this universality property holds for light PBHs produced in the early Universe, and consider the evaporation of distributions that take the form of Eq. (7). Note that this simplifying assumption, and our code is able to handle any general distribution f PBH (M, a ⋆ ). For concreteness, we will consider either the spin distribution obtained in Ref. [94] or simply take a Gaussian distribution centered around a ⋆ = 0.7 and having width σ a⋆ . Other mechanisms of PBH production, such as Q-balls, oscillons, or cosmic strings, may as well be created with specific spin distributions [91-93, 97, 98]. However, giving a comprehensive review of the effect of different distributions is beyond the scope of this work. In this paper, we focus on the spin distribution from hierarchical mergers as a simple case study, to exhibit the general phenomenological effects of spin distributions.

III. PRIMORDIAL BLACK HOLE EVAPORATION
In this section, we review the physics of PBH evaporation. For the sake of generality, we consider the case of a Kerr PBH, that is characterized by both its mass and angular momentum. The evaporation of a Schwarzschild PBH will simply correspond to the special case of vanishing angular momentum. In what follows, we do not consider effects from dynamical horizons on the particle production given the lack of consensus on which metric is the correct one to describe cosmological black holes, see [99] and references therein for the status of this topic. In other words, we assume that the Hawking evaporation operates in the same way as in vacuum.
We denote the mass of a PBH and its dimensionless spin parameter by M and a ⋆ = JM 2 p /M 2 respectively, where J is the angular momentum of the black hole. The Hawking emission rate for a particle species i with threemomentum p, total energy E i and g i internal degrees of freedom is given by [26,27] where In this expression, Ω = (a ⋆ /2GM )(1/(1 + 1 − a 2 ⋆ )) is the horizon's angular velocity, and l, m are the total and axial angular momentum quantum numbers, respectively. respectively. The central plot shows the power law mass distribution with Mc = 10 g and Mc10 σ = 10 6 g for α = 1, 2, 3 shown in violet, teal and yellow, respectively. The right plot shows the mass distributions generated through metric preheating- [88] (taken from Appendix C, Fig. 8) and critical collapse [82][83][84], using Mc = 10 5.6 g to align the two distributions.
In Eq. (8), it is necessary to write explicitly the sum over the angular momentum quantum numbers as the PBH spin breaks the spherical symmetry, given the explicit dependence of the Hawking rate on m. For the case of Schwarzschild PBHs, a ⋆ = 0, the emission rate is independent of m so one can replace the sum over m by a factor of (2l + 1). σ lm si appearing in Eq. (8) corresponds to the BH absorption cross section, which describes the effects of the centrifugal and gravitational potential on the particle emission [26,27,[100][101][102]. These quantities are determined by solving the spin-dependent equations of motion for a field in the Kerr spacetime. We have adopted the methodology described in Refs. [103][104][105]] to obtain such cross-sections.
To derive evolution equations for the BH mass and spin, we multiply the Hawking rate shown in Eq. (8) by the total energy of a given particle E i or by the m quantum number, and then we integrate over the phase space. Defining the evaporation functions for mass and angular momentum, ε i (M, a ⋆ ) and γ i (M, a ⋆ ) per particle i, respectively, as where g i is the particle species i's internal degree of freedom and we can sum over all existing species to obtain the following system of coupled equations [30,106,107] We refer the interested reader for further details of the derivation of these equations to Refs. [30,48]. In this equation, we introduced the total evaporation functions for mass and angular distribution, Later in the paper, we will use similar notations with additional subscripts 'SM', 'DM', and 'DR', indicating that the sum over i is restricted respectively to Standard Model, dark-matter, and dark-radiation species only. In general, obtaining the time evolution of a Kerr PBH mass and spin parameters requires using numerical methods. In our code, we solve numerically the system of equations Eqs. (12) including the dependence of the evaporation functions on the mass of the emitted particles. Thus, our approach can be readily extended to include any number of degrees-of-freedom. For instance, the production of an extended sector having a large number of degrees-of-freedom, as explored in e.g. Refs. [108,109], could be included which could lead to interesting phenomenology.
For the sake of numerical efficiency, we will, however, restrict our analysis to evaporation products with masses smaller than the initial PBH Hawking temperatures. Note that in the case where the evaporation functions ε and γ only depend on the PBH spin, that is, in the limit where the evaporation products can be considered to be exactly massless, there is a formal solution of these equations, as first described in Ref. [101]. For completeness, we briefly revisit this derivation in the App. A.

IV. EVAPORATION OF AN EXTENDED DISTRIBUTION
In this section, we derive the Boltzmann and Friedmann equations that allow tracking the evolution of PBH abundance and the Universe's energy density throughout cosmic history, for an evaporating distribution of PBH.
We start by defining f PBH (M, a ⋆ , t) as the distribution of PBH of mass M , spin a ⋆ at time t. The number density of PBHs, n BH , is then given by Similarly, one can write the comoving energy density of the Black Hole distribution as Demanding that the comoving number of PBHs is conserved throughout cosmic history 1 , we obtain the continuity equation H being the Hubble parameter. Taking the time derivative of Eq. (15) and using the relation in Eq. (16) (see App. A) one obtains the Friedmann-Boltzmann equatioṅ which has to be solved simultaneously with the equation describing the evolution of the Standard Model plasma, ρ SM , To solve these equations numerically, one needs to evaluate, at every time t, the integrals in the right-hand sides of Eqs. (17) and (18). However, the distribution f PBH (M ) is only known at the time of PBH formation t in . In order to obtain the values of these integrals, it is useful to change variables and map the mass spectrum at time t to the corresponding initial masses M in , defined such that In that case, the conservation of the infinitesimal PBH comoving number density 2 provides that where a stands for the scale factor. For future con- . The Boltzmann equations that are solved, after defining ϱ BH ≡ a 3 ρ BH and ϱ SM ≡ a 4 ρ SM , arė Note that in these expressions, dM/dt remains a function of M = M (t, M in , a in ⋆ , t in ) and a ⋆ = a ⋆ (t, M in , a in ⋆ , t in ) which can be found numerically by integrating Eqs. (12). Moreover, one should pay attention to the fact that at the end of the evaporation dM/dt diverges, and thus the mass integrals in the right-hand sides of Eqs. (12) have to be restricted to initial masses that have not yet reached 0 at time t.

V. NUMERICAL IMPLEMENTATION
Numerically, to solve the system of Eq. (21) we adopt the following procedure: We assume that the PBH distribution is formed simulatenously when the Universe is radiation dominated after inflation, and we take as formation time t in -or temperature T in -when the particle horizon mass is equal to mass of the initial peak of the distribution M c with κ ∼ 0.2. We stress that this is a simplifying assumption since it is expected that each formation mechanism will predict different initial conditions, such as t in . We note however that our code can be easily modified to specify the initial conditions and the initial time t in . Alternatively, if one had a PBH production mechanism that continued to form black holes while some smaller PBHs were evaporating, more substantial modifications of FRISBHEE would be required. The double integration is performed at every time step in the range t ∈ [t in , t fn ], t fn being the time when the whole PBH population has evaporated. We then track the evolution of the PBH population and the expansion dynamics of the Universe by solving the Friedmann and Boltzmann equations of Eq. (21), making use of the pre-calculated 2D integrals. We have implemented this numerical strategy in the publicly available code 4 FRISBHEE which follows these steps and is able to include SM and BSM particles in the Hawking evaporation spectrum considered. We have double-checked that we recover the results from Refs. [30,31,48] when we take the monochromatic limit for a given mass and/or spin distribution. To our knowledge, this is the first implementation of mass and spin distributions into the Boltzmann and Friedmann of the early Universe. Although the tool, BlackHawk [110,111] includes the facility to compute Hawking radiation from a distribution of black holes, it does not solve the cosmic evolution alongside BH evaportation. That's because BlackHawk's main focus is the accurate determination of primary and secondary spectra from black holes that are currently evaporating, M ∼ 10 14 g, which with existing constraints will simply be a background matter density in the early Universe.
Note that the relation we used in Eq. (22) does not hold in full generality. It is true to very good accuracy for evaporation products that are much lighter than the Hawking temperatures of the PBHs in the distribution. However, a PBH's lifetime starts deviating from its true value as soon as a large number of new degrees of freedom with masses larger than the Hawking temperature of the PBH are introduced in the spectrum. As a matter of fact, when adding 200 new degrees of freedom with masses m ≪ T BH , this equation is verified to 0.005% accuracy. It reduces to a 30% error for m > T BH . For the sake of numerical efficiency, we have therefore implemented this approximation in FRISBHEE to improve its calculation speed. However, we stress that employing the lifetime of PBHs for any interval [M, M +dM ]×[a ⋆ , a ⋆ +da ⋆ ], which our code is capable to perform, would not require a significant modification, although it would imply a larger run time.
It is possible to trace the evolution of the distribution at any time by computing the relevant Jacobian, As described in Appendix A-C, in the case of relativistic evaporation products, we find an explicit expression for the Jacobian We present in Fig. 2 an example of the time evolution of the mass only and mass and spin distributions as function of time. In the case of only mass distribution (top panels), we present the time evolution of a log-normal (left), power-law (middle) and critical collapse (right) scenarios, for different values of time, as indicated in the figures by the different color scale. Note that we recover, for the power-law distribution, the results derived analytically in Ref. [70]. In the case of a mass and spin distribution (lower panels) we present three different snapshots for the initial time t in (left), a time equal to 10% of the lifetime of the peak τ (middle) and t = τ (left). The shaded region represents the PBH parameters which cannot be present in the Universe since such PBHs would have a lifetime shorter than t − t in . Note that, although it looks like all these distributions get depleted with time, that the comoving number density of PBHs (including those whose masses is zero after evaporation) is conserved in such simulations. Indeed, the evaporation of a black hole is a process that accelerates with time and, therefore, the mass of PBHs that start evaporating quickly runs towards zero. Numerically, the depletion of the distributions that is visible in the figure is balanced by a delta function for mass and spin located at the origin and increasing with time. Note that in practice, FRISBHEE stops tracking the evolution of the PBH mass once it falls below 10% of the Planck Mass and sets it to zero, but keeping and tracking Planck relics in the code is possible. We have included some animations of the time evolution of mass only and mass and spin distributions for benchmark parameters. These are included in the wiki 5 page of FRISBHEE's GitHub.

VI. COSMOLOGICAL IMPRINTS OF PBH DISTRIBUTIONS
In this section, we explore the possible consequences that the evaporation of an extended distribution of PBHs can have in cosmology, and how one can attempt to probe the shape of such a distribution using cosmological observables. In Ref. [30,31,48], it was shown that the evaporation of PBHs can leave various imprints in darkmatter searches, small-scale structures, and dark radiation measurements. In the next sections, we will see how the modification of the cosmological background's evolution arising from the extension of the PBH distribution can affect such observables.
In fact, it is expected that any mechanism that involves the production of particles out of equilibrium (such as e.g. leptogenesis, see Ref. [51]) throughout cosmic history can be affected by such dynamics. Gravitational waves induced at second order in perturbation theory constitute interesting signatures of a possible PBH dominated era. Arising from the successive onset of the PBH domination and PBH evaporation, they could also reveal interesting signatures regarding the shape of the PBH distribution [66]. We will briefly discuss this possibility, and leave an extensive study of these different imprints for future works.

A. PBH Domination
In the case of a monochromatic distribution, f PBH (M ) = δ(M − M in ), PBH domination can be assimilated to an early matter-domination period that ends when the whole population of PBHs evaporates. Such sudden evaporation can have many consequences in cosmology, and possibly leave observable imprints in cosmological data. The existence of a PBH dominated era typically affects the dynamics of the Universe expansion as compared to the temperature evolution of the SM sector, as well as the time-evolution of the Universe's equation of state parameter. If primordial black holes form in the early Universe with a sufficient energy fraction, they can end up dominating the energy density of the Universe before they evaporate. When this is the case, their evaporation corresponds to an injection of entropy into the SM sector that can strongly affect its dynamic with respect to the Universe's expansion. In this section, we explore the effect that the extended PBH distributions have on the dynamics of evaporation and corresponding entropy injection into the SM bath. For the sake of simplicity, this section only considers Schwarzschild PBHs. Note that the presence of a spin distribution would however affect the following results, as spinning PBHs would have the tendency of evaporating faster. However, we believe that the qualitative dynamics of reheating that we present are not to be changed drastically. Indeed, the universality of the spin distribution considered in Eq. (7) guarantees that the lifetime of PBHs with different masses would be shortened in a relatively uniform way. Besides happening at a slightly earlier time, the qualitative dynamics of the Universe's reheating when PBHs evaporate would not be significantly affected by the shape of the spin distribution.
In Fig. 3, we present our results for the different mass distributions listed in Sec. II A. We draw the evolution of the PBH relative abundance Ω PBH ≡ ρ PBH /ρ tot as a function of the number of e-folds N ≡ log a, with a denoting the scale factor. In the case of a log-normal distribution (left panel), we vary the width σ, introduced in Eq. (1), and we fix the central mass M c to be 10 4 g. In the small-σ limit, one can observe that evaporation happens quickly, similar to the case of a monochromatic distribution. On the other hand, increasing σ corresponds to smearing the evaporation over larger time scales, and one can see that the evaporation can then happen over several e-folds of expansion. A similar observation can be made in the case of the power-law distribution introduced in Eq. (2). In the central plot of Fig. 3, we used different powers α corresponding to the formation of PBHs during matter domination (α = 2), kination (α = 0), and a nearly inflating Universe (α ≳ 1). We fixed the mass M c = 10 g and varied the width of the distribution by taking M c × 10 σ = 10 2 g, 10 3.5 g, and 10 6 g. As expected, increasing the width of the distribution leads to extend-ing the duration of the evaporation process.
It is interesting to note that for certain choices of α (namely 1 < α ⩽ 3), as shown in the middle panel of Fig. 3, the total black hole abundance initially falls but then remains essentially fixed at a non-zero value over an interval lasting many e-folds before finally dropping to zero. It turns out that this feature is a specific instance of a more general phenomenon called cosmic stasis [112], wherein the abundances of the different energy components of the Universe remain fixed over an extended interval despite cosmological expansion. In this case, the total PBH abundance is remaining fixed because its natural tendency to grow in this Universe as a result of redshifting effects is precisely counterbalanced by the loss of PBH energy density into radiation via Hawking evaporation [70,113]. It is remarkable that distributions with power 1 < α ⩽ 3 are all attracted to such a regime.
Finally, in the rightmost plot of Fig. 3, we illustrate the evaporation dynamics which is obtained when using the distributions corresponding to the critical collapse and metric preheating cases, as described in Sec. II A. Here both distributions peak towards heavy masses, and thus the evaporation process is dominated by the end tail of the distribution.

B. Dark Matter Relic Density
In this section, we investigate the impact of the evaporation of a PBH population with both mass and spin extended distributions on DM production. Besides the equations for the evolution of the SM radiation and PBH energy densities, we need to track the number density of DM particles produced from the evaporation. Here, we assume that such particles only interact gravitationally, and, although there might be other production mechanisms, we focus only on the DM from Hawking evaporation. Similarly to the Eqs. (17), we can write the equation for the DM comoving number density N DM aṡ where now the evaporation function includes the contribution of the DM particle, ε DM . At this point, we do not fix the spin of the particle, to make our discussion as general as possible. In order to compare with previous results that consider a monochromatic distribution, we parametrize the initial PBH energy density via where ρ in PBH is related to the mass and spin distribution as in Eq. (15). After the complete evaporation of the PBH population, we obtain the DM relic abundance in the standard manner. We present in Fig. 4 the PBH parameters which produce the observed relic abundance assuming a scalar DM particle with a mass of 1 GeV for a log-normal distribution (left) and power law, critical collapse and metric preheating (MP) scenarios (right) while considering a monochromatic spin distribution with a ⋆ = 0. For comparison purposes, we have included in both panels the results from a purely monochromatic distribution. In the case where the mass distribution is a log-normal, we observe a crucial dependence on the width of such distributions. For values σ ≲ 0.1 we find that the initial PBH energy density needed to produce the correct abundance is quite similar to the monochromatic case. Meanwhile, if the distributions are wider, that is, for σ > 1.0, the initial PBH fraction is significantly modified. Such modifications occur because wider distributions contain a population of much heavier PBHs than the peak value M c , these will produce much more DM particles because the number of emitted particles grows ∝ (M in ) 2 in the case where DM mass is smaller than the initial PBH temperature m DM < T in BH , cf. Ref. [30]. This implies that, for initial PBH densities which do not lead to a PBH domination, the β ′ value that leads to the observed DM abundance is reduced by a factor of ∼ O(10 4 ) for σ = 2 in comparison with a Monochromatic mass distribution centered at M = M c . Furthermore, the place where the PBH domination occurs, which leads to vertical relic abundance contours, is shifted to lower values of M c as the σ becomes larger, due to the presence of heavier PBHs. On the right panel of Fig. 4 we present the relic abundance contours for a power law distribution with α = 2.0, σ = 2 (orange), critical collapse (green) and metric preheating (red). In the case of a power law mass distribution, we observe a similar effect as in the log-normal case: the presence of heavier PBHs enhances the DM production, thus requiring a smaller initial PBH fraction. In contrast, for the critical collapse and metric preheating cases, the relic contours do not differ in a noticeable way from the monochromatic case. Since such distributions extend to lower PBH initial masses instead of higher ones, the DM production is dominated by the PBH having a mass M c , corresponding to the peak of the distribution. Furthermore, in both panels of Fig. 4 we overlay thick black dashed lines when the PBH produced DM is no longer consistent with observation and is too warm to support the small-scale structures observed in our Universe.
Of course, when DM particles emerge from the black holes, they are highly boosted. Through cosmic redshift these particles cool, and depending on their masses, can constitute warm (∼ 10 −2 − 10 3 GeV) or cold DM (≳ 10 3 GeV) [35,40,43,48,114]. Previously, the mass constraints on m DM have been applied for monochromatic distributions in mass and spin of the PBH. Since now our calculations allow for non-monochromatic distributions, we have adapted our interpretation of the warm DM constraints accordingly. Once again we make use of the computational tool CLASS [115][116][117] to determine the matter power spectrum in the CMB, where the input required for the DM phase space distribution, f DM , is where g DM is the number of degrees of freedom, p is the three-momentum, and n DM is the DM number density. With a PBH distribution, the particle emission rate per momentum is given by where similarly to above d 2 N DM /dp ′ dt ′ is a timedependent function of M = M (t, M in , a in ⋆ , t in ) and a ⋆ = a ⋆ (t, M in , a in ⋆ , t in ). The t fn here refers to the latest evaporation time we consider for a given distribution as described in section V. In practice we determine f DM by running FRISBHEE to evaluate F in and a(t), which is then used to integrate numerically Eq. (29). Also, the timeindependent NCDM temperature is needed to interface the resulting DM distribution with CLASS, where T ev and T in are the SM plasma temperatures at evaporation and PBH production respectively. Of course, because we are now working with distriutions of PBHs the evaporation temperature takes multiple values, T ev here is simply the SM plasma temperature when the entire distribution has been evaporated, T (t fn ).
To determine whether a specific DM distribution is at odds with observations of structure in the Universe, such as those of the Lyman-α forest [118,119], we use CLASS to calculate the matter power spectrum P (k), quantifying the deviation from CDM by way of the transfer function, T (k), defined in where k is the wavenumber. The scales at which T (k) can start to stray from 1 has been determined by the parameter fit where µ is dimensionless exponent which is fixed to µ = 1.12 as in Ref. [120] and α is the breaking scale, which we take to be saturated at α = 1.3 × 10 −2 Mpc h −1 [40,[119][120][121].
Returning to the results shown in Fig. 4 we see that the warm DM constraints tend to be more constraining for wider distributions such as log-normal with σ ≳ 1 and a power law. This is because the heavier PBHs evaporate later, providing less time for DM to redshift, and it is exactly these heavier PBHs that produce the lion's share of DM. For the examples which reproduce a similar result for the monochromatic, the WDM constraints are almost identical. We would like to note, that taking the 1D fit for Eq.(31) may not be the most appropriate since when we take the 2D fit, we see some variation in the µ parameter, away from the quoted µ = 1.12. We leave a more sophisticated analysis for future work, but since the α parameter really controls where the matter power spectrum diverges from CDM, we are confident that the results we show here are correct to a reasonable degree of precision.
Due to the computational expense of performing the 3D integral, we have opted to only perform the above analysis in the Schwarzchild case, and we refer readers to Refs [38,43,48] to get an impression of how spin a in ⋆ , β ′ and M c parameters will affect the relic lines. We expect that Kerr distributions of BHs only have an appreciable effect on the relic lines for spin-2 DM.
The main takeaway of this section is that broader distributions may well enable smaller β ′ values to produce the correct relic density, but this itself will introduce more aggressive warm DM constraints.
For DM having a mass larger than the initial PBH temperature, we find that the overall behaviour is similar to the monochromatic scenario, i.e., depending on the DM mass, the relic density contours are either directly or inversely proportional to M c ; the change of behaviour occurs when T in BH = m DM and it is abrupt [30]. Such change of behaviour is much more gradual for wider distributions.

C. Dark Radiation
Similar to our previous discussion on DM generation, in this section we consider the emission of massless states that will modify the number of relativistic neutrino species, N eff . To do so, we solve an equation for the Dark Radiation (DR) comoving energy densitẏ where now ε DR corresponds to the contribution to the total evaporation function coming from the DR. We track the evolution of all Universe species, and, after the full evaporation of the PBH population occurring when the plasma temperature is T ev , we determine the modification to N eff , ∆N eff as [32] ∆N eff = 8 7 where N SM eff = 3.045 are effective number of neutrinos [122], and T eq = 0.75 eV is the matter-radiation equality temperature.
Let us focus first on the specific case where the DR particle is the spin-2 massless graviton, potentially the most well-motivated undiscovered particle in fundamental physics. In Fig. 5 left, we present the modification on N eff coming from hot gravitons produced from the evaporation. Here, we assumed a close-to-monochromatic mass distribution, σ = 0.005, and considered the scenarios where the spin distribution is gaussian (colored full lines) and for the case of the fourth generation merger distribution (red dashed line), taken from Ref. [94]. For the gaussian case, the color indicates the value of the σ a⋆ ∈ [0.0, 0.2], where the value of zero corresponds to the monochromatic in spin situation. The black dashed line corresponds to the monochromatic Schwarzschild case, included here for comparison. We observe an enhancement on ∆N eff for wider gaussian distributions since such wider distributions can increase the amount of emitted DR in comparison to the monochromatic case, depending on where the mean value lies. For instance, a Gaussian distribution with σ a⋆ = 0.1, and mean value ⟨a ⋆ ⟩ = 0.71 will generate ∼ 34% more gravitons than a monochromatic distribution with the same central value. This leads to an increase of ∼ 18% in ∆N eff for σ a⋆ = 0.2. Interestingly, for broader distributions, the contribution to ∆N eff is reduced since those start to include PBHs which have a smaller angular momentum, for which the emission of gravitons is significantly decreased. For other DR particles with lower spins, this dependence is less prominent. We have found that for vectors, the enhancement only reaches ∼ 1% for gaussian spin distributions with σ s⋆ = 0.2; for fermions and scalars the modification is smaller. Such behavior is due to the mild dependence that the greybody factors have on the angular momentum parameter a ⋆ for particles having lower spins [101].
If we consider extended mass distributions instead, we observe a shift in the shape of the contribution to ∆N eff , see the right panel of Fig. 5. The black dashed lines in the same figure indicate the parameters excluded by BBN constraints. We obtain such a limit by determining whether such PBH populations reheat the Universe after evaporation with a temperature T ev smaller than the lower BBN bound of ∼ 5 MeV [123][124][125][126] Since the overall population would evaporate much later than in the monochromatic case, especially for broad distributions, the value of the evaporation temperature, T ev , is reduced. Therefore, the contribution from the dark radiation becomes more substantial. Assuming a log-normal distribution, we find that the contribution to ∆N eff is shifted by a factor of ∼ O(10 3 ) to lower masses for σ = 1.5. For other types of mass distributions, such as the power-law, the behavior is similar. Moreover, for distributions such as critical collapse and metric preheating (MP) the shift is less sizable since the evaporation temperature of such distributions is closer to the monochromatic one.

D. Gravitational Waves from Evaporation
There are several reasons why the presence of PBHs in the early Universe is expected to be accompanied by a sizeable spectrum of gravitational waves (GWs). First, it is clear that the production of PBHs, if it arises from the collapse of primordial perturbations, has to be accompanied by a large scalar power spectrum at a given scale. Such scalar perturbations are known to induce GWs at second order in perturbation theory [127][128][129][130]. However, the specific shape of the corresponding spectrum is entirely dependent on the shape of power spectrum considered, and therefore, it is not a unique prediction that can be obtained given a particular PBH distribution. Another interesting manner via which PBHs can produce gravitational waves is when they dominate the energy density of the Universe before they evaporate. In the latter case, the sudden transition from a matter to a radiation-dominated Universe, as the PBHs evaporate rapidly at the end of their lives, can cause the gravitational potential to oscillate (the so-called poltergeist mechanism) [64][65][66][131][132][133][134]. In this context, the faster the evaporation takes place, the sharper the transition between matter domination and radiation is, and the higher the peak of GWs expected to be probed by future observatories is. Therefore, it is expected that the extension of the mass and spin distributions, which mainly smears out the evaporation process (as we have seen in Sec.VI A), has the tendency to suppress the GW spectrum. This can already be seen in Ref. [64] -see also Ref. [65]-where it was shown that a width of order σ = 10 −2 for a log-normal distribution leads to a suppression of the GW spectrum of over three orders of magnitude. Another interesting aspect of the PBH formation and subsequent domination is that it can act as a source of isocurvature perturbation in the early Universe which also leads to a peak in the GW spectrum, at a different frequency than the effect described previously [50,66,131]. Similarly, we expect that such a peak would be suppressed in the presence of a broader distribution of PBHs. Indeed, PBH of different masses are expected to form at different times, and the isocurvature perturbations that would be sourced by a smeared forma-tion process are likely to lead to a broadened GW spectrum. For a fixed energy fraction of PBHs, this broadening of the spectrum has thus to be accompanied by a suppression of its overall amplitude. Such claim should in principle be verified numerically. However the study of induced gravitational waves from isocurvature perturbation requires a dedicated study that we plan to explore in a future work.
The last source of GWs arising from the evaporation of PBHs comes from the direct production of gravitons via Hawking evaporation [51,135,136]. The GW signal for a monochromatic mass spectrum with M in ∼ 1 (10 4 ) g is expected in this case to peak at frequencies of order 10 13 (10 15 ) Hz. While there are currently no technologies that can detect such high-frequency GW, there are proposed detectors that could, in principle, detect THz GWs [137][138][139][140][141][142].
In the left panel of Fig. 6, we study the effect of having an extended PBH mass distribution on these highfrequency GWs. We use M c = 10 4 g, and an initial PBH energy fraction of β ′ = 10 −7 , which are also parameters leading to a poltergeist signal detectable by LISA. We find that the GW spectrum due to direct Hawking radiation of gravitons, is not strongly affected by a nonmonochromatic PBH distribution. The general trend is that for wider PBH distributions, the peak of the stochastic gravitational wave background (SGWB) shifts to either lower or higher frequencies, depending on whether the PBH population contains more light or heavy PBHs than the monochromatic case; moreover, the amplitude of the SGWB is not affected. Of the various observables we have studied, the SGWB produced from Hawking radiation is thus the least sensitive to the population of PBHs having a non-monochromatic spectrum, at least in terms of amplitude. This is interesting as it constitutes a robust signature of the existence of PBHs evaporating in the early Universe, as opposed to the other aforementioned GW signals. Furthermore, the fact that the peak frequency can change by orders of magnitude is important because it has implications on the physics reach of proposed THz GW detectors, which will have different frequency regions where they are optimally sensitive.
For the case of GWs emitted by Kerr PBHs, Ref. [136] first computed the spectrum showing the modification coming from the enhancement due to the BH spin. We have reproduced their results for monochromatic mass and spin distributions. In the right panel of Fig. 6, we show the effect a population of highly spinning PBHs has on the SGWB signal. In order to exhibit how all effects intersect with each other we show the signals from monochromatic (black dashed), Power-Law distributed in M but monochromatic in a ⋆ (blue), monochromatic in M but Gaussian distributed in a ⋆ (orange), and Power-Law in M and Gaussian in a ⋆ (green). As was shown in Ref. [136], the amplitude of the SGWB is enhanced substantially when a ⋆ is close to maximal, comparing the left panel to the right, we observe that the peak amplitude is ∼ 10 2 larger. The peak frequency is at a similar position, but generally, the spectral shape of the SGWB is quite different thanks to the nature of the graviton emission spectrum as such high a ⋆ values. We see that introducing a Gaussian in distributed a ⋆ weakens the strength of the SGWB signal because fewer PBHs are maximally spinning. When this effect is in concert with mass distributed PBH population, we can see that the effect is similar to the Schwarzschild case, whereby the mass distribution causes the SGWB signal to shift in frequency but the amplitude remains unchanged.

VII. SUMMARY AND CONCLUSION
In this paper, we have explored the effect of PBH mass and spin distributions on various cosmological observables. We outlined our numerical method, which simultaneously evolves the mass and spin distributed PBH population in time and focused on well-studied mass (lognormal, power law, critical collapse and metric preheating) and spin distributions (spin distribution from mergers and Gaussian). While the results are unsurprising, to our knowledge, this work is the first attempt at numerically solving the evaporation of both the mass and spin of PBHs together with the Friedmann-Boltzmann equations that describe the time evolution of different Universe's components. We found that non-monochromatic mass and spin distributions reduced how suddenly the transition from PBH to radiation domination occurred. Naturally, the wider the distribution the more slowly such a transition occurs. In the case of light fermionic DM pro-duction from PBHs, we studied Schwarzschild PBHs with a non-trivial mass distribution. We found that broader mass distributions required a smaller initial number density to produce the observed relic density, as compared to a monochromatic distribution centered at M = M c . Intuitively, this occurs as a broader mass distribution includes heavier PBHs, and thus requires less to produce the same overall energy density in PBHs. Likewise, we found that the warm DM constraint was significantly more stringent for broad distributions because the heavier PBHs would evaporate later and thus diminishing the effect of redshifting the boosted DM. We studied the effect of a non-trivial mass and spin distribution on ∆N eff . We fixed the mass distribution and allowed the width of the spin distribution to vary and vice versa. Widening the spin distribution tends to produce an increase in ∆N eff as there are a higher proportion of higher spin PBHs which are efficient at producing gravitons. The effect of increasing the spin width , σ a⋆ , from 0 to 0.2 has the effect of increasing the contribution to ∆N eff by ∼ 10%. On the other hand, for a fixed central mass of PBH, for a Schwarzchild PBH, having a non-trivial mass distribution can increase the contribution to ∆N eff by a factor of a few (see the right panel of Fig. 5). Finally, we consider the effect finite mass and spin PBH widths had on the SGWB produced from direct Hawking radiation of the PBHs. We found that, unlike for other GW signals predicted from the evaporation of PBHs, the amplitude of the graviton production via Hawking emission is mainly insensitive to the extension of the PBH mass distribution. This suggests that the production of GWs from the evaporation of PBHs is amongst the most robust signals that could be searched for with observations in the future regarding the existence of light PBHs in the early Universe. Our code FRISBHEE is publicly available at the address https://github.com/yfperezg/frisbhee.