Precision Calculation of Dark Radiation from Spinning Primordial Black Holes and Early Matter Dominated Eras

We present precision calculations of dark radiation in the form of gravitons coming from Hawking evaporation of spinning primordial black holes (PBHs) in the early Universe. Our calculation incorporates a careful treatment of extended spin distributions of a population of PBHs, the PBH reheating temperature, and the number of relativistic degrees of freedom. We compare our precision results with those existing in the literature, and show constraints on PBHs from current bounds on dark radiation from BBN and the CMB, as well as the projected sensitivity of CMB Stage 4 experiments. As an application, we consider the case of PBHs formed during an early matter-dominated era (EMDE). We calculate graviton production from various PBH spin distributions pertinent to EMDEs, and find that PBHs in the entire mass range up to $10^9\,$g will be constrained by measurements from CMB Stage 4 experiments, assuming PBHs come to dominate the Universe prior to Hawking evaporation. We also find that for PBHs with monochromatic spins $a^*>0.81$, all PBH masses in the range $10^{-1}\,{\rm g}


I. INTRODUCTION
Black hole evaporation via the emission of Hawking radiation is a well established phenomenon [1,2], with recent work towards precisely characterizing the Hawking radiation yields of relevant particles and the time evolution of the population of black holes (e.g. [3]). Primordial black holes (PBHs) are of particular interest in that their possible mass range spans many orders of magnitude and they could be relevant to the questions of dark matter and cosmological chronology, and their existence can affect observable quantities that can be probed with current (and future) cosmological experiments. Here, we undertake a precision study of Hawking evaporation of PBHs prior to Big Bang Nucleosynthesis (BBN), with particular attention to PBH spin and spin distributions, the PBH reheating temperature, and the evolution of the number of relativistic degrees of freedom, and compare our results to the current sensitivities from the cosmic microwave background (CMB) and BBN, as well as future CMB Stage 4 experiments.
PBHs may have formed in the early Universe from the collapse of primordial density inhomogeneities originating from quantum fluctuations prior to inflation or from topological defects such as cosmic strings or domain walls. Bubble collisions during a first-order phase transition can also trigger PBH formation. For a recent review of PBH formation mechanisms, we refer to [4] and references therein.
The spin of the resulting PBH population depends on the equation of state (as does the mass distribution). PBHs formed during radiation domination are believed to have negligible spin [5]. On the other hand, PBHs formed during an early matter-dominated era (EMDE) [6][7][8][9] could have sizeable to near-extremal spin [10][11][12]. PBHs can also accumulate some spin either through early accretion processes [5] or through hierarchical mergers [13]. In the last two decades, constraints have been placed on a wide range of PBH masses, assuming Schwarzschild (non-rotating) PBHs with monochromatic mass spectra 1 (for a review see e.g. [4]). Using a combination of numerical and analytical results for Hawking radiation, recent studies have started to complete the constraints on PBHs with non-zero spin [16][17][18][19][20][21][22].
Here we study the production of dark radiation in the form of gravitons coming from Hawking evaporation of populations of spinning PBHs prior to BBN. We compute the primary and secondary spectra of Standard Model (SM) particles and gravitons for realistic spin distributions of PBHs from an EMDE [11] as well as a hierarchical merger history [13]. Our calculations are performed with the public code BlackHawk [3], developed by a subset of the current authors 2 . The evolution of a given distribution of PBHs and the associated time-dependent spectrum of emitted gravitons are computed, allowing a straightforward determination of the total energy emitted in the form of 1 The distribution of PBHs can also have an extended mass function, for example if the power spectrum of primordial inhomogeneities embeds a wide peak around some spatial scale [14,15]. Extended mass functions of spinning PBHs have not yet been thoroughly studied. We leave this for future work. 2 We have implemented the possibility of adding a particle to the SM, e.g. the massless spin 2 graviton or general dark sector particles of spin 0, 1, 2 or 1 2 in BlackHawk, although in this study we focus only on massless spin 2 graviton emission. Additional dark sector particles have not yet been implemented in the public version of the code. To our knowledge, this is the first precision calculation of Hawking radiation with with non-trivial PBH spin distributions using BlackHawk. dark radiation. This affects the number of relativistic species, with the result characterized as the deviation from the SM expectation of the effective number of neutrino species, ∆N eff . We compute ∆N eff , and compare it to existing results in the literature and interpret it in the context of current limits on ∆N eff from BBN and CMB measurements. In particular, we carefully calculate BBN constraints on the dark radiation density using AlterBBN [23,24].
The main application of our results is the calculation of ∆N eff from PBHs that were formed during an EMDE and subsequently came to dominate the Universe prior to Hawking evaporation. EMDEs are highly motivated due to the ubiquity of moduli in string theory and have been extensively studied in recent years in the context of dark matter [25][26][27][28][29][30] and baryogenesis [31]. Detailed studies of PBHs formed during an EMDE have been performed by [7][8][9], with a focus on long-lived PBHs existing in the current Universe, and their interplay with dark matter physics. PBHs that evaporated before BBN are harder to constrain 3 . The authors of [11,12,35,36] have initiated much progress in this direction; of particular relevance for our work are the formation rate [11] and spin distribution [12] of PBHs formed during an EMDE. Following the spin distributions used in [12] as benchmark examples, we find that PBHs formed during an EMDE with a spin distribution due to the first-order effect are constrained by current CMB bounds on ∆N eff in the mass range 10 8 − 10 9 g; they are completely constrained in the mass range 10 −1 − 10 9 g by projections of CMB Stage 4 experiments. PBHs that formed during an EMDE with spin distribution due to the second-order effect, on the other hand, are not constrained by current BBN or CMB bounds on ∆N eff ; they too would, however, be completely constrained in the mass range 10 −1 − 10 9 g by CMB Stage 4 projections (Fig. 3).
The fact that PBHs formed during an EMDE that evaporate before BBN will be completely probed by ∆N eff measurements from CMB Stage 4 experiments is the main result of our work. Physically, this happens because PBHs formed during an EMDE are endowed with significant spin, which enhances their production of gravitons during evaporation. It should be noted that the ∆N eff constraints are only relevant if the PBHs come to dominate the Universe. Generally, this is quite restrictive on the sector that causes the EMDE. We consider a gravitationally coupled modulus that causes the EMDE and obtain conditions on the decay width (and hence the modulus mass) such that this condition holds. In terms of the modulus sector, our result is that for a variety of PBH spin distributions and fractions β of the total energy density of the Universe that is constituted by PBHs at formation time during an EMDE, moduli with masses larger than ∼ 10 8 GeV will be constrained by CMB Stage 4 experiments (Fig. 8).
We also consider the case of a spin distribution due to inspirals of PBHs under a heirarchical merger history, obtaining, for the first time, precision predictions for ∆N eff in this scenario, which will be probed by CMB Stage 4 experiment. Finally, we go on to apply our results to the case of PBHs with extremal spins regardless of origin, and find that PBHs with spin a * 0.99 and mass M BH 10 8 g are excluded by CMB stringent constraints (TT,TE,EE+low E) while those with even higher spin a * 0.999 are constrained by the CMB conservative constraints (TT+low E), but only for masses M BH 2 × 10 8 g. We further determine that the limiting value of the PBH spin that will be constrained by CMB Stage 4 experiment for all PBH masses up to 10 9 g is a * min, all 0.81 . Our paper is organized as follows. In Section II, we give an overview of the formation and evaporation of Kerr PBHs. In Section III we outline the precision calculation of the effective number of neutrino species, ∆N eff , from PBH evaporation, addressing spin distributions and the reheating temperature in Subsections III A and III B, respectively. We present the bulk of our results in Section IV. In Section IV A, we compare precision results for benchmark spins versus spin distributions, including the effects of the reheating temperature and a precision accounting of the effective degrees of freedom. In Section IV B, we explicitly focus on spin distributions relevant for an EMDE. The effect on BBN is discussed in Section V, and our conclusions are given in Section VI. Finally, we include three appendices, where we discuss the details of PBH formation and evaporation during an EMDE, PBH spin distributions from an EMDE, and PBH spin distributions from inspirals.

II. KERR PRIMORDIAL BLACK HOLES: FORMATION AND EVAPORATION
Hawking has demonstrated that black holes evaporate [1,2] by emitting quasi-thermal radiation with a temperature for the Schwarzschild solution, and 3.046, where N eff is the total number of relativistic degrees of freedom and 3.046 is the SM expectation, from PBH evaporation. The precision calculations involve two steps: taking into account the distribution of PBH spins and carefully defining the reheating temperature. We also use a precise expression for the number of accessible degrees of freedom. First we review the standard calculation of ∆N eff . Using conservation of entropy during the expansion of the Universe, one can track the evolution of the energy density of dark radiation from reheating to matter-radiation equality. For a population of PBHs with lifetime τ , the age of the Universe at formation is small relative to τ such that the evaporation time is t eva τ . Assuming instantaneous thermalization of SM particles at the end of PBH evaporation, the reheating temperature, T RH , can be obtained as where ρ PBH (τ ) is the energy density of PBHs at the time of evaporation, ρ DR (ρ SM ) is the amount of energy PBHs emit in the form of dark radiation (SM particles), and g * (T ) denotes the total number of relativistic degrees of freedom at temperature T , given by Here the sum includes all bosonic (B) and fermionic (F ) degrees of freedom with temperatures of T B and T F , respectively. The density of PBHs at evaporation is related to the density of PBHs at formation, usually expressed in terms of the fraction of the energy density of the Universe that collapsed into PBHs at PBH formation time, which is denoted by β. In this work, we assume that β is sufficiently large such that the energy density of PBHs exceeds that of radiation at some time before evaporation. A discussion of such a scenario is given in Appendix A 1 in the case of modulus decay. With this hypothesis, the density of PBHs at evaporation is fixed by the fact that SM radiation produced by PBH Hawking evaporation constitutes the main component of SM radiation at reheating. Thus, tracing the redshifted temperature of the CMB today back to reheating (from today back to the matter-radiation equality time with a(t) ∼ t 2/3 and then to the reheating time, t RH τ , with a(t) ∼ t 1/2 ), we obtain the value of T RH . The values we obtain for ∆N eff in this study should be considered as upper limits in the case of full PBH domination prior to evaporation. The constraints are generally weakened but must be recalculated if PBHs do not dominate the energy density of the Universe before evaporation.
The energy density of SM radiation (all relativistic particles) is therefore diluted as where a RH(EQ) is the scale factor at reheating (matter-radiation equality), and g * ,S (T ) counts the number of relativistic degrees of freedom contributing to the entropy, given by Similarly, the energy density of dark radiation, ρ DR , also dilutes as Therefore, the ratio of the energy density of dark radiation to the SM radiation energy density at matter-radiation equality becomes which determines the effective number of neutrino species as [42] ∆N eff = ρ DR (t EQ ) ρ R (t EQ ) N ν + 8 7 A. Extended PBH spin distributions A monochromatic distribution of non-rotating PBHs is only a convenient approximation to the more realistic extended mass distribution of rotating PBHs generated by detailed models of PBH formation, accretion, and mergers. For the purpose of this study, we focus on single-mass, rotating PBHs with a spin number distribution n(a * ) normalized to unity, Note that the assumption of a monochromatic mass distribution for PBHs is justified if the PBH production occurs at a precise time, leading to a very narrowly peaked mass distribution. The total energy that has been emitted in the form of dark radiation by the reheating time t RH can be expressed as a ratio over the SM emission, i.e. the ratio of the energy densities after evaporation is complete, where ρ DR/SM (t RH ) is the total emission integrated over the history of the Universe prior to reheating, and The emission rates for individual species, d 2 N i /dtdE, come from Eq. (5). We stress that the ratio (18) takes into account the fact that for high DR emission, which occurs for highly spinning black holes, the approximation ρ BH ρ R used in [20,22] no longer holds. This could be one of the reasons our results differ from those of [22] for high PBH spin. We recall that ρ SM = ρ R at time t RH (which occurs before matter-radiation equality), which allows to use the ratio (18) in Eq. (15) to determine ∆N eff . Furthermore, we note that the normalization of the density of PBHs ρ BH is irrelevant to the computation of ∆N eff , since it cancels out of the ratio f DR in Eq. (18).
For the purposes of this study, we have implemented in BlackHawk the possibility of including additional particles beyond those in the SM, e.g. the massless spin 2 graviton or general dark sector particles of spin 0, 1, 2 or 1 2 (although here we focus only on gravitons). We compute the evolution of a given distribution of PBHs and the associated time-dependent spectrum for this additional particle 6 . It is then straightforward to integrate over this spectrum to obtain the total energy emitted in the form of SM particles and additional dark radiation and hence the ratio f DR in Eq. (18).
The main effect of a spin distribution, relative to monochromatic spin, is to modify the rate of emission of dark radiation, and thus its ratio to SM radiation, as in Eq. (18). Indeed, it is well known that spinning black holes emit more high spin particles (s i = 1 or s i = 2) than non-spinning black holes. As we consider the emission of spin 2 massless gravitons, this effect can be quite sizeable, with the emission being enhanced by a factor of up to ∼ 10 4 [43]. The effect of this enhancement on the ratio (18) is somewhat less dramatic, since the emission of spin 0, 1, and 1 B. Reheating temperature and degrees of freedom When an extended spin distribution of PBHs is employed rather than a monochromatic spin distribution, there is some subtlety in defining the reheating temperature. As spinning black holes emit more radiation than non-spinning ones, with a continuous increase in the emission as a * increases, they evaporate faster. Although initial nonzero spin has a small effect on black hole lifetime (somewhat less than 60% diminution for extremal spin [43,47]), it does influence the way one defines the reheating time. For PBHs with lifetime τ (e.g. for PBHs with monochromatic mass and spin distributions), assuming an instantaneous reheating in Eq. (10) is justified by the fact that PBHs emit most of their Hawking radiation during a period of time that is negligibly small relative to their lifetime. However, since black holes with higher spin evaporate faster than black holes with lower spin, a distribution in initial spins causes a spread of the evaporation times and a non-instantaneous reheating scenario.
For simplicity, here we consider two possibilities for the definition of the reheating time: 1. the reheating time corresponds to the time at which the last PBHs (with the lowest spins) evaporate; and 2. the reheating time corresponds to the average PBH lifetime, weighted by the spin distribution, We believe that the second option is more physically realistic, as the averaged lifetime corresponds roughly to the peak of the emission of the Hawking radiation. We discuss both options in Section IV, where we present our results. Finally, in order to obtain the ratio (15), it is necessary to specify the quantities g * (T ) and g * ,S (T ). We stress that precise determination of these numbers of degrees of freedom are model-dependent, especially for the region of temperatures close to the QCD phase transition. Here, that corresponds to M BH ∼ 7 × 10 8 g (T ∼ 100 MeV). Refs. [22,32] use step functions which give results qualitatively similar to ours, while the model used in [20] is not made explicit and shows a significantly different behaviour. In this work, we use the tabulated values of g * (T ) and g * ,S (T ) available with the public code SuperIso Relic 7 [48,49].

IV. PRECISION RESULTS FOR ∆N eff
Here, we present precision results for ∆N eff with improvements to the calculation as described above. In Subsection IV A, we explore the effect of each of the three precision elements we have included here; spin distributions, reheating temperature, and degrees of freedom. In Subsection A 2 we present, for the first time, explicit predictions for ∆N eff from PBH spin distributions expected from an EMDE.
In all cases, we compare our results for ∆N eff to current experimental limits and projected sensitivities of future experiments. We present three relevant CMB constraints/sensitivities: two are taken from the Planck Collaboration [50] and are denoted in the plots as CMB 1 (TT+low E, conservative) and CMB 2 (TT,TE,EE+low E, more stringent). The third one is the sensitivity of the future CMB Stage 4 (CMB-S4) experiment, and represents an order of magnitude improvement over current limits (see details in [51][52][53]). Where relevant, we also include the constraint on ∆N eff from BBN, as discussed in Section V.

A. Benchmark spin scenarios -exploring precision results
In this subsection we compute ∆N eff , incorporating the precision calculations described above -spin distributions, reheating temperature, and degrees of freedom -for some benchmark PBH spin scenarios. We compare the results for ∆N eff calculated with an extended spin distribution to those obtained from monochromatic spin distributions (e.g. the central/peak value of the extended distribution), as well as ∆N eff obtained with the two reheating temperature calculations. Furthermore, we compare our results to previous calculations in the literature for a * = 0 and a * = 0.99 to demonstrate the full effects of the precision calculation.
We first make a few comments about PBH masses in the low mass regime. In our calculations, we find that changing the PBH mass in the range 10 −1 g < M BH < 10 9 g has a very small effect on the ratio ρ DR /ρ R (less than 1% over the whole mass range). This is because, for a given spin distribution, the main variation in ∆N eff as the PBH mass is varied comes from the different reheating times (and thus reheating temperatures). Below M BH 10 5 g, the reheating temperature is far above the mass of all the SM particles (T RH 10 2 GeV), so g * (T ) and g * ,S (T ) have already reached their asymptotic values. Thus, ∆N eff values for M BH 10 5 g can be safely extrapolated from their value corresponding to the case of M BH = 10 5 g. We note that our results also apply to the M BH = 10 −5 −10 −1 g mass range for PBHs. This range is sometimes excluded from analyses due to model-dependent limits on the inflationary Hubble parameter [22,54]. Below, we present results only for 10 5 g ≤ M BH ≤ 10 9 g.
To show how the prediction for ∆N eff from an extended distribution of PBH spins compares to the monochromatic approximation, we present two benchmark extended spin distributions, along with the corresponding prediction assuming a monochromatic distribution. We first consider the asymptotic spin distribution expected for multiple generation PBH inspirals [13] (see Appendix A 3 for details). The average spin in this case is a * 0.7, so we compare the results for the full spin distribution to those for the monochromatic spin distribution with a * = 0.7. As discussed above, we expect more gravitons to be emitted because there are higher spin PBHs in the extended distribution, relative to the monochromatic case. This is borne out in the results shown in the left panel of Fig. 1. We see that ∆N eff indeed does acquire greater values (by ∼ 25%) for the full distribution than for the monochromatic one. This discrepancy becomes critical for PBH masses above 7 × 10 7 g; in the case of the extended distribution, these PBHs will be probed by CMB-S4, while the average spin approximation leads to the conclusion that only PBHs with masses above 2 × 10 8 g would be accessible to CMB-S4.
In the right panel of Fig. 1, we show the results for a benchmark extended distribution from an EMDE with a * 0.64, along with a monochromatic distribution with a * = 0.64 (more details on this are discussed in Section IV B and Appendix A 2). 8 The spin distribution in the right panel of Fig. 1 due to early matter domination is significantly different from that in the left panel due to inspirals. In particular, this EMDE spin distribution is less symmetric and much more broad than the inspiral distribution. The relative discrepancy between the extended distribution and the monochromatic distribution is therefore even greater (∼ 60%) in the right panel than in the left panel of Fig. 1. For this EMDE extended spin distribution, one finds that PBHs with masses above ∼ 7 × 10 7 g will, in fact, be probed by CMB-S4. This conclusion stands in stark contrast to that inferred under the assumption of a monochromatic spin distribution at the peak or average spin. The results for the extended spin distributions in both panels of Fig. 1 are also shown for the two prescriptions for calculating the reheating temperature, as discussed in Section III; (1) instantaneous reheating at the evaporation time of the last PBH (with the lowest spin) is shown in grey, and (2) the weighted average PBH evaporation time using Eq. (21) is shown in black. In both panels, one can see that prescription (2) results in a shift in the ∆N eff curve to higher PBH mass relative to the results assuming prescription (1). This can be understood on the basis of the reheating temperature from prescription (1) being smaller than the reheating temperature from prescription (2). Indeed, higher spin PBHs evaporate faster, and are better accounted for in prescription (2). Thus, one could achieve the same reheating temperature (and therefore the same ∆N eff ) with prescription (1) by assuming a higher PBH mass.
In Fig. 2, we compare the values of ∆N eff obtained with precision calculations using BlackHawk to recent calculations in the literature. In the left panel of Fig. 2, we consider a * = 0 (Schwarzschild), and compare the ∆N eff from BlackHawk (solid) with those calculated in Refs. [20] (denoted as [H20], dashed) and Ref. [32] (denoted as [M20], dot-dashed) updated in Ref. [22] (denoted as [M21], dotted) with the use of BlackHawk. The relative discrepancies in these cases are ∼ 10%. In the right panel of Fig. 2, we consider a * = 0.99, and compare with the results of Ref. [20], where we find a ∼ 20% discrepancy. As discussed in Section III B, an important difference between our results (solid) and other calculations in the literature is that here we take the values for g * (T ) and g * ,S (T ) tabulated in the public code SuperIso Relic [48,49]. Near M BH ∼ 8×10 7 g (corresponding to T RH ∼ 100 MeV), the number of degrees of freedom is very sensitive to the QCD equation of state, and the precise behavior of ∆N eff is evident. That said, using a simple step function for g * (T ) and g * ,S (T ) gives results qualitatively similar to ours [22,32]. This precision calculation reveals that highly spinning PBHs with a * = 0.99 and masses M BH 2 × 10 8 g that dominated the Universe before BBN are, in fact, already excluded by CMB 2 constraints on ∆N eff [50].

B. Early matter domination and extremal spins
In this subsection, we present the results of precision calculation of ∆N eff for PBH spin distributions from a period of early matter domination. This is the first time a prediction for ∆N eff from PBHs produced during an EMDE has been calculated. Here, we assume that PBHs produced during the EMDE come to dominate the Universe by the time of Hawking evaporation. The validity of this assumption depends on the physics behind EMDE; as an example, we consider the conditions under which this happens when an EMDE is caused by a gravitationally coupled modulus field in Appendix A 1.
We use the spin distributions from Ref. [12] as benchmarks, the details of which are discussed in Appendix A 2. An- Upper left: ∆N eff results for PBH distributions formed during an early matter domination era due to first order effect [12] with σH = {0.01, 0.05, 1} (dotted, dot-dashed and solid respectively). Upper right: A zoom into the region where ∆N eff is constrained by current CMB limits. Lower left: ∆N eff results for PBH distributions formed during an early matter domination era due to second order effect [12] with σH = {0.01, 0.05, 1} (dotted, dot-dashed and solid respectively). Lower right: A zoom into the region where ∆N eff shows brutal change due to the step shape of g * (TRH). The black (grey) curves correspond to instantaneous reheating at the weighted average value of the black hole lifetimes (last black hole evaporation). The 95% C.L. limits on ∆N eff from CMB (shaded areas) are taken from [50] (CMB 1 : TT+low E, CMB 2 : TT,TE,EE+low E); the 95% limit from BBN is computed in Section V; the prospective CMB-S4 constraint (horizontal dashed line) is extracted from [20]. gular momentum within a comoving region of space has two components; the first-order contribution ("the first-order effect") originating from deviation of the boundary of the volume from a sphere, and the second-order contribution ("the second-order effect") sourced by density fluctuations in the comoving region (for a detailed treatment, we refer to [12]). The first-order effect usually dominates (when the initial deviation of the boundary of collapsing region from a sphere is large), but an almost spherical initial collapsing region can diminish the first-order effect and make the second-order effect the dominant one.
In Fig. 3, we present ∆N eff results for PBHs formed during an EMDE, with spin distributions due to first-and second-order effects in the upper and lower panels, respectively. In each panel, we show results for three different values of σ H , the mean variance of the density perturbations at horizon entry where δ 2 s is the variance of the density perturbations integrated over the volume of a sphere and t H is the time of horizon entry (for further details, see Appendix A 2). σ H controls the shape of the spin distribution, as well as the peak location. For both the first-and second-order effects, larger σ H leads to more broad spin distributions. Increasing σ H also shifts the peak of the second-order distribution away from a * = 1 to smaller values of a * . As mentioned in Section IV A, for σ H = 0.1, the peak average of the spin distribution from the second-order effect is located at a * = 0.64. Note that either the first-or second-order effects could dominate, as discussed in Appendix A 2.
The upper panels of Fig. 3 show ∆N eff due to spin distributions dominated by the first-order effect. We see that for σ H small enough, the largest PBH masses are already excluded by CMB 2 , and in some cases even CMB 1 . In the upper right panel, we see in detail that for σ H 0.01, M BH > 10 8 g are excluded by CMB 2 constraints. For EMDE spin distributions dominated by the first order effect, the entire PBH mass range 10 −1 g < M BH < 10 9 g will be probed by CMB Stage 4.
It is clear from the lower panels of Fig. 3 that PBH spin distributions from an EMDE are not constrained by current CMB or BBN limits on ∆N eff if the spin distribution is dominated by the second-order effect. However, these would be probed by CMB Stage 4 measurements. We can see in the lower left panel of Fig. 3 that for σ H small enough ( 0.1), all PBH masses in the range 10 −1 − 10 9 g will be probed by CMB Stage 4. For the value σ H = 0.1, only PBHs in the high mass end of this range 3 × 10 7 − 10 9 g will be accessible to CMB Stage 4.
Another noticeable feature in all panels of Fig. 3 is the shift of ∆N eff towards higher PBH masses if one takes reheating time as the average weighted lifetime τ (black curves) compared to the time of evaporation of the last PBH (grey curves). This is consistent with what was observed in Fig. 1. This shift is most sizeable for extremal spin distribution for which the average spin is a * ∼ 1, i.e. small σ H . This is especially clear in the lower right panel of Fig. 3, which zooms in to the region of strong variation of ∆N eff in the lower left panel. The difference in these results due to the different prescriptions for reheating time particularly affects ∆N eff in the mass range M ∼ 5 − 9 × 10 7 g. This is also the region where ∆N eff is most affected by the precise shape of g * (T ) and g * ,S (T ).
We next turn to an investigation of ∆N eff for near-extremal PBH spins. In Fig. 4, we present ∆N eff for monochromatic spin distributions with a * 0.9, under the assumption that the PBHs dominate the energy density of the Universe before BBN. As in Fig. 2, we see that these highly spinning PBHs are already excluded by current CMB 2 constraints for large enough PBH masses. We also see that the excluded mass range grows as the spin increases, due to the shorter lifetime of spinning PBHs. Furthermore, for the largest PBH spins we consider, the increase in ∆N eff , which is due to the enhanced emission of high spin particles (spin 2 most of all), saturates. Indeed, the Hawking emissivity of near extremal PBHs does not grow to infinity as a * → 1 but instead saturates.
To be specific, we see in Fig. 2 that the future CMB Stage 4 measurements will be sensitive to extremal values of PBH spins a * 0.9. From the right panel of Fig. 4 it is evident that PBHs with spin a * 0.99 and mass M BH 10 8 g are excluded by the CMB 2 stringent constraints. PBHs with even higher spin a * 0.999 are constrained by the CMB 1 conservative constraints, but only for masses M BH 2 × 10 8 g. This is, to our knowledge, the first constraints put on light spinning PBHs from ∆N eff from current CMB limits.
Finally, we explore the capability of the CMB Stage 4 experiment to explore PBHs with monochromatic spins, under the assumption that PBHs dominated the energy density of the Universe prior to BBN. In Fig. 5, we present the smallest monochromatic spin for which CMB Stage 4 will be sensitive to the entire mass range considered here, a * min, all , as well as the largest monochromatic spin for which CMB Stage 4 will not be sensitive to any part of the mass range, a * max, no . We find that the smallest monochromatic spin for which CMB Stage 4 will be sensitive to the whole range of masses is a * min, all 0.81 . For a monochromatic spin distribution with a * > a * min, all , CMB Stage 4 will probe all PBH masses 10 −1 g < M BH < 10 9 g. On the other hand, the smallest monochromatic spin value for which CMB Stage 4 can constrain any of the PBH masses is a * max, no 0.69. For a * a * max, no , the entire mass range would be inaccessible to CMB Stage 4, while for a * a * max, no only the heaviest PBHs (M BH ∼ 10 9 g) will be probed. While the results in Fig. 5 apply to monochromatic spin distributions, the same question can in principle be answered for various types of extended spin distributions. As discussed in Section IV A and demonstrated in Fig. 1, one can expect a ∼ 25 − 60% relative discrepancy between the ∆N eff prediction for monochromatic spin distributions relative to the extended distributions we consider here. Indeed, for a scenario such as early matter domination, which induces a particular spin distribution for PBHs, one could even explore the range of cosmological parameters that yield ∆N eff to which next generation experiments will be sensitive.

V. EFFECT ON BBN
In this Section, we outline how the dark radiation yield from light PBH evaporation can affect BBN. If PBHs evaporate before the onset of BBN (t eva 1 s or M BH 10 9 g), then the emitted SM particles thermalize to the expected plasma density and provide no measurable effect on BBN. The dark sector, which is also emitted by Hawking radiation, however, provides an additional source of density in the Friedmann equations compared to standard BBN. This sector does not interact with the SM, thus its temperature is decoupled from the plasma temperature. However, the dark radiation can be treated as an additional effective number of neutrinos N eff during BBN and up to the time of photon decoupling. Thus, the ∆N eff constraints from BBN can be used to constrain the dark radiation density before BBN.
We use the public code AlterBBN [23,24], which computes the abundances of the light chemical elements in alternative cosmological scenarios, such as with the addition of a dark radiation density. Comparison with the fiducial values for these abundances, in particular 2 H and 4 He measured in old gas clouds, provides constraints on ∆N eff .
The master parameter for BBN is the baryon-to-photon ratio η, which is related to the reduced baryon cosmological parameter Ω b h 2 via where k = 100 km/s/Mpc is the Hubble parameter scale today, M Pl the Planck mass, m b is the average baryon mass, and T 0 = 2.7255 K is the CMB temperature today. The constraints on ∆N eff are computed with Ω b h 2 as a free parameter. Its central value is Ω b h 2 = 0.0224 [50]. Inside AlterBBN, the observational values of the chemical element abundances used to obtain the updated ∆N eff constraints are The improved nuclear rate for D + p → 3 He + γ by LUNA has been implemented into the code [57]. In Fig. 6, we present the BBN constraints on ∆N eff with Ω b h 2 as a free parameter. The light and dark shaded red regions correspond to the 68% and 95% confidence level regions obtained using the 2 H and 4 He BBN constraints, as recomputed with AlterBBN for this work (these are the BBN constraints used in Figs. 3 and 4). The dot-dashed, dashed, and solid contours correspond to the 95% confidence level regions obtained by the Planck Collaboration [50] (dot-dashed: TT,TE,EE+lowE, dashed: TT,TE,EE+lowE+lensing, solid: TT,TE,EE+lowE+lensing+BAO), with the vertical and horizontal dashed grey lines marking the standard values of the parameters Ω b h 2 = 0.0224 and ∆N eff = 0, respectively. If PBHs evaporate during or after BBN, then the effects are much more complicated and require careful treatment, beyond the scope of the current study. First, the energetic hadronic emission just before BBN can trigger p ←→ n interconversion and thus modify the p/n ratio at the beginning of BBN. This ratio strongly affects the final 4 He abundance and is thus severely constrained. Second, hadronic injection (mesons) during BBN can trigger nuclear reactions through hadrodissociation and can modify the abundance of intermediate light elements. This may modify the final 2 H abundance and is also severely constrained. Third and last, the emission of energetic photons at the end of BBN can still destroy BBN products through photodissociation and thus can modify the final abundances before recombination. All these phenomena are associated with the evaporation of M 10 9 g PBHs. We refer the interested reader to [4,46,[58][59][60][61][62][63] for detailed analyses of these.

VI. CONCLUSIONS
Our purpose in this paper has been to conduct precision studies of dark radiation emanating from spinning PBHs. We have concentrated on the case of gravitons. Our precision study incorporated spin distributions of PBHs and a careful treatment of the reheating temperature and relativistic degrees of freedom. We studied the impacts of each of these three precision elements on the calculation of ∆N eff due to graviton emission from PBHs, and applied the calculation to a scenario with extended PBH spin distributions due to an early matter dominated era (EMDE).
There are two main effects related to incorporating extended PBH spin distributions relative to monochromatic spin distributions. First, since a BH's lifetime is related to its spin, a spin distribution will result in a distribution of evaporation times. The second, dominant, effect is that PBHs with high spins emit more particles with higher spins, i.e. gravitons. So a spin distribution that extends to higher spins will result in more graviton emission relative to a corresponding monochromatic spin approximation, and thus a larger prediction for ∆N eff .
In undertaking a precision study, we find that it is also important to consider a precise formulation for the number of relativistic degrees of freedom as a function of temperature. We show that different characterizations for the degrees of freedom lead to different conclusions regarding experimental sensitivity to various models. In fact, for PBHs with masses M BH few × 10 7 g that dominated the Universe before BBN, one finds very different predictions for ∆N eff . Different prescriptions for the reheating temperature due to PBH evaporation also lead to variations of ∆N eff . These are relatively small in comparison to the other effects considered, but careful attention to the reheating temperature is relevant to make a precise statement regarding experimental sensitivity for some PBH masses.
Our main application was to study gravitons coming from Hawking evaporation of PBHs created during an EMDE. If such PBHs come to dominate the Universe prior to final evaporation, the resulting dark radiation can be probed by current BBN and CMB constraints, as well as future CMB Stage 4 experiments. We have found that PBHs with spin distribution due to the first-order effect are constrained by current CMB bounds on ∆N eff in the mass range 10 8 − 10 9 g, and would be completely constrained in the mass range 10 −1 − 10 9 g by CMB Stage 4 projections. PBHs formed during an EMDE with spin distribution due to the second-order effect, while not constrained by current BBN or CMB bounds on ∆N eff , would be completely constrained in the mass range 10 −1 − 10 9 g by CMB Stage 4 experiments for all scenarios except for the largest σ H considered here. In terms of the modulus sector, we found that for a variety of PBH spin distributions and fractions β that have been considered in the literature, moduli with masses larger than ∼ 10 8 GeV will be constrained by CMB Stage 4 experiments.
We also explored ∆N eff for near-extremal PBH spins. We find that if PBHs with monochromatic spin distributions with a * 0.99 dominate the energy density of the Universe before BBN, current CMB constraints exclude PBHs with masses m BH 10 8 g. As the spin increases toward 1, ∆N eff increases until it saturates, since the Hawking emissivity of near extremal PBHs does not grow to infinity as spin approaches 1 but instead saturates. We therefore find that for increasing a * the minimal PBH mass excluded by current CMB measurements is shifted to lower PBH masses until saturation. We also find that for PBHs with monochromatic spins a * > 0.81 that dominated the energy density of the Universe prior to BBN, all PBH masses in the range 10 −1 g < M BH < 10 9 g will be probed by CMB Stage 4 experiments.
Note: Near the completion of this work, the authors became aware of the publication of Ref. [22], where the author considers Hawking radiation of light Kerr PBHs in the early Universe, in the mass range 10 −5 −10 9 g. Ref. [22] considers the emission of light dark matter particles by Kerr PBHs, as an extension of the results of [33] for Schwarzschild PBHs, as well as the effect of emission of dark radiation by light Kerr PBHs, as considered in this work. We compare the results in [22] to ours and others in the literature in subsection IV A.

Appendix A: Early matter dominated eras and PBH spin distributions
In this Appendix, we discuss possible PBH spin distributions n(a * ) that are motivated by early Universe cosmology. These distributions will then be used in Eq. (18) to obtain f DR . We will focus mainly on two benchmark scenarios: a period of early matter domination, possibly by a string modulus and scenarios in which PBHs acquire spin by inspirals. 7. Thermal history of Universe reheated at t = tRH, includes radiation (R), and a modulus field (φ). After reheating radiation dominates energy density of the Universe, so Universe experiences a radiation-dominated (RD) epoch until energy density of φ dominates at t = t φ and initiates a matter-dominated (MD) era. During this modulus-dominated era, spinning PBHs may form at t = t f . At t = t dec , φ stops to oscillate around its minimum and decays into radiation. After t = t dec Universe enters a RD epoch which may lead to a MD era at tPBH at which energy density of PBHs takes over if the lifetime of PBHs is long enough. PBHs will eventually deposit their energy content into the thermal bath at t = teva due to Hawking evaporation. This is the onset of another RD epoch which will continue until radiation-matter equality time. Time intervals in the plot are for demonstrative purposes only and not indications of actual times.

PBH formation during an early matter dominated era and subsequent evaporation
In usual studies of an early matter domination era (EMDE) phase, the scenario is the following: after inflationary reheating, the Universe is filled with radiation and a modulus field, φ. We will be agnostic about the origins of φit could be a string modulus. We will assume that it couples to other fields via gravity only. Under fairly general assumptions, it is possible that φ is displaced from the minimum of its potential and starts to oscillate. Since energy density of modulus field redshifts like energy density of matter, it eventually dominates the energy density of the Universe and causes a transition from a radiation-dominated era to a matter-dominated era. During this modulusdominated epoch, spinning PBHs can form. Modulus field will finally decay into radiation, reheat the Universe for a second time and give rise to a radiation-dominated era. Since energy density of PBHs also redshifts like matter, they can eventually dominate over radiation and lead to a matter-dominated epoch. In this case, after evaporation, their contribution to ∆N eff is not negligible.
To evaluate the initial abundance of PBHs for which a once modulus-dominated Universe may lead to a PBHdominated epoch, one needs to trace back the evolution of energy density of each component to the onset of modulusdominated era (see Fig. 7). We assume that following the reheating of Universe at t = t RH , energy density of the modulus field becomes comparable with energy density of radiation at t = t φ , i.e., ρ R (t φ ) ρ φ (t φ ). Afterward a fraction β of the total energy density of the Universe collapses into PBHs at t = t f , i.e., β ≡ ρ PBH (t f )/ρ tot (t f ). Subsequently, the modulus field decays instantaneously into radiation at t = t dec , and eventually PBHs evaporate at time t = t eva . Then, to guarantee a PBH-dominated era, we need to make sure that at some time t = t PBH , where t dec t PBH t eva , we have ρ PBH (t PBH ) ρ R (t PBH ). This leads to where a is the scale factor, and we assume that β 1. Since t φ < t f < t dec , ignoring a(t φ )/a(t dec ) and a(t φ )/a(t f ) can cause an overestimation up to a factor of 2. Demanding t PBH t eva (or equivalently a(t PBH ) a(t eva )) provides a lower bound on β given by where τ φ and τ PBH are the lifetimes of modulus field and PBHs respectively. To evaluate Eq. (A2), we use the fact that since β c corresponds to the case that PBHs dominated energy density almost at the time of their evaporation, between decay of the modulus field and evaporation time, t dec t t eva , the Universe undergoes a radiation-dominated stage. We also use t RH τ φ , t f τ PBH , and Γ φ m 3 φ /M 2 Pl . The initial abundance of PBHs, β, that formed during a matter-dominated epoch and gained angular momentum due to the first-and second-order effects (see Subsection A 2 for details), is calculated as a function of the mean variance of density perturbations at horizon entry, σ H , by Ref. [12]. A certain value of σ H can give rise to a PBH-dominated For sufficiently small black holes (M BH 10 10 g), the lifetime of a spinning black hole can be expressed as [47] τ BH = c( a * ) where c( a * ) depends on the average of the spin of the black hole, so here is a function of σ H , and is calculated by BlackHawk. By combining Eqs. (A3) and (A4), we obtain where for β(σ H ) we follow the numerically calculated curves in Fig. 5 of Ref. [12]. The authors have checked that the following semi-analytic formulae reproduce the behavior: where I = 1 and q = √ 2. For details, we refer to Ref. [12]. Fig. 8 displays the lower bound on the mass of the modulus field which can later lead to a PBH-dominated era for benchmark values of σ H that we use in this paper, for both the first-and second-order effects.

PBH spin distributions from an early matter dominated era
In the early Universe, density fluctuations, δ = δρ/ρ, grow after they enter the cosmological horizon. In a radiationdominated epoch, if density fluctuations are greater than a threshold, they can collapse into a PBH with mass bounded by the total mass within the horizon. In a matter-dominated epoch, the absence or significant reduction of the pressure gradient force enhances PBHs formation rate and it is the deviation from spherical symmetry that governs the probability of PBH formation [64].
Since in cosmological perturbation theory, the rotational mode is not growing to linear order, the effect of rotation in the formation of PBHs is naively expected to be unimportant. As a matter of fact, detailed calculation [12] shows that angular momentum plays a very important role in the formation of PBHs in the matter-dominated phase. Here we briefly review the importance of rotation in PBH formation in a matter-dominated epoch and the spin distribution of these PBHs by following the theory of angular momentum in structure formation adopted in [12].
Angular momentum within a comoving region of space has two components; the first-order contribution originating from deviation of the boundary of the volume from a sphere which can be described by an ellipsoid, and the secondorder contribution sourced by density fluctuations in the comoving region. Assuming different modes carry random phases, the variance of the angular momentum within a sphere of comoving radius r 0 can be written as where and subscripts 1 and 2 represent the first-and second-order contributions respectively. In the above expressions, M = (4π/3)ρ 0 (ar 0 ) 3 is the mass inside the spherical region of interest and R = ar 0 is the physical radius of the region, q is the dimensionless parameter of the initial reduced quadrupole moment of the mass, I is of order unity, and δ 2 s ∼ t 4/3 is the variance of δ s , the density perturbation integrated over volume of the sphere. By normalizing them at the time of horizon entry, t = t H , we have where σ H ≡ δ s (t H ) 2 1/2 . The corresponding dimensionless angular momentum can be estimated as .

(A11)
The value of angular momentum grows with time until nonlinearity becomes important. After this moment which is the time of maximum expansion, t max , linear perturbation theory is not valid any longer. The collapse of the overdense region begins and it becomes separated from the evolution of the Universe. Therefore after t max angular momentum approaches a constant value. By demanding δ s ( t max ) 2 1/2 = 1, the average value of t max can be estimated as and accordingly, the average value of the first-and second-order angular momenta are given by (1) 1/2 2 5 (2) 1/2 The dominant component is chosen as the final angular momentum; a * 2 1/2 max a * 2 (1) 1/2 . Only a minority of masses with a * 2 1/2 1 (σ H 0.1) can overcome centrifugal force and collapse directly to PBHs. Therefore, angular momentum strongly suppresses formation of PBHs and most of the PBHs are rapidly rotating at the time of formation. By comparing the first-and second-order angular momenta in Eq. (A12), we see that the magnitude of q, which quantifies initial deviations of the collapsing region from a sphere, determines dominant effect; a large q (large initial deviation from a sphere) leads to first-order dominance, on the other hand a small q (an almost spherical initial collapsing region) makes the second-order effect the dominant one.
In spite of the complicated dependence of angular momentum on different coupled modes, a hypothesis facilitates obtaining the distribution function for spins; since both L 2 (1) and δ 2 s include self-coupling of single modes while L 2 (2) consists of the coupling of two independent modes which are not parallel to each other, it is reasonable to assume that |L (1) | ∝ δ s and |L (2) | ∝ δ 2 s 1/2 δ s , or more precisely By using t max = t H δ s (t H ) −3/2 , a * can be evaluated as a * max a * (1) , a * (2) where a * The fact that a smaller δ s (t H ) leads to a larger final value for a * , can be explained by noticing that t max ∝ δ s (t H ) −3/2 . Hence for a smaller δ s (t H ), it takes a longer time to reach the nonlinear phase and consequently angular momentum has a longer time to grow. The finite duration of the early matter-dominated epoch puts a lower bound on δ s (t H ). Demanding t max < t end , where t end marks the end of the early matter-dominated era, leads to δ s (t H ) ≥ δ fd ≡ (t H /t end ) 2/3 for PBHs formation. The other lower bound on δ s (t H ) is set by requiring a * ≤ 1 or equivalently δ s (t H ) ≥ δ th(1) ≡ 3×2 2 5 3 q 2 and δ s (t H ) ≥ δ th(2) ≡ 2 5 Iσ H 2/3 . All of these conditions can be summarized as δ s (t H ) ≥ max(δ th(1) , δ th (2) , δ fd ). If δ fd < δ th (2) , the effect of finite duration is negligible, otherwise PBHs formation and the probability of formation of PBHs with large spin are severely suppressed. In this paper we assumed that δ fd < δ th (2) . It can be shown that δ th(1) (q c ) = δ th(2) (q c ) where q c ≡ 2/3(5/2) 7/6 I 1/3 σ 1/3 H . Since for δ s (t H ) = 5/3Iq −1 σ H we have a * (1) = a * (2) , there is a transition point, a * t = (2/5)(3/5) 3/4 I −1/2 q 3/2 σ −1/2 H , at which the behaviour of a * is changing or in terms of q c , a * t = (q/q c ) 3/2 ; a q > q c leads to a * t > 1 which is not acceptable. Since where N (µ, σ 2 ) represents a Gaussian distribution with mean µ and variance σ 2 . Therefore one can describe the spin distribution of PBHs with the following piecewise distribution n(a * ) = 1 N    n (1) (a * ) n (2) (a * t ) n (1) (a * t ) 0 ≤ a * < a * t , n (2) (a * ) a * t ≤ a * ≤ 1 , where n (1) (a * )da * ∝ 1 a * 3 exp − and N is the normalization factor.

PBH spin distribution from inspirals
In the early stages of the evolution of Universe, a sufficiently large ensemble of PBHs may experience mergers if the binary capture rate becomes larger than the expansion rate of the Universe and the inspiral phase ends prior to the Hawking evaporation of PBHs. Ref. [65] has studied different timescales which are relevant to mergers in a population of PBHs in early Universe. The merger rate could be enhanced if PBHs form in clusters, a hypothesis that will be testable in future experiments looking for CMB µ-distortion, as proposed recently [66]. If PBHs undergo several mergers before evaporating, the angular momentum gained during each merger causes the spin distribution of PBHs to converge to a universal distribution that is relatively independent of the mass of PBHs, the initial spin distribution of the first generation of PBHs, and the number of merger generations [13]. Although Ref. [13] considered solar mass black holes, their study is also applicable to PBHs.
The universal hierarchical merger spin distribution in [13] has been shown numerically to appear after four merger generations, and to peak at a * 0.7, with nonzero support over 0.4 a * 0.9. To understand why this universal spin distribution emerges, one needs to consider major contributions to the spin following a merger which consist of the individual spins of the two individual black holes, and the orbital angular momentum of the binary. For equal mass binary black holes the orbital angular momentum dominates over the contribution from the individual spins. Numerical simulations show that merger of non-spinning binary black holes of equal mass will result in a final black hole with a * 0.6864 [67]. The spins of the binary black holes can become important and even cancel the orbital angular momentum if they are sufficiently large and anti-aligned to the orbital angular momentum, and the mass ratio needs to be sufficiently small. This is basically why major mergers (with mass ratio ∼ 1) give rise to black holes with high spin distributions, peaked at a * 0.7.
A slightly different hierarchical merger spin distribution is reported by [68] based on the priors from LIGO/VIRGO data for mergers limited to the Milky Way. This distribution also peaks at a * 0.7.