Primordial Black Hole Evaporation and Dark Matter Production: I. Solely Hawking radiation

Hawking evaporation of black holes in the early Universe is expected to copiously produce all kinds of particles, regardless of their charges under the Standard Model gauge group. For this reason, any fundamental particle, known or otherwise, could be produced during the black hole lifetime. This certainly includes dark matter (DM) particles. This paper improves upon previous calculations of DM production from primordial black holes (PBH) by consistently including the greybody factors, and by meticulously tracking a system of coupled Boltzmann equations. We show that the initial PBH densities required to produce the observed relic abundance depend strongly on the DM spin, varying in about $\sim 2$ orders of magnitude between a spin-2 and a scalar DM in the case of non-rotating PBHs. For Kerr PBHs, we have found that the expected enhancement in the production of bosons reduces the initial fraction needed to explain the measurements. We further consider indirect production of DM by assuming the existence of additional and unstable degrees of freedom emitted by the evaporation, which later decay into the DM. For a minimal setup where there is only one heavy particle, we find that the final relic abundance can be increased by at most a factor of $\sim 4$ for a scalar heavy state and a Schwarzschild PBH, or by a factor of $\sim 4.3$ for a spin-2 particle in the case of a Kerr PBH.


I. INTRODUCTION
The entire catalogue of experimental evidence for dark matter (DM) comes only from its gravitational effects. Despite this, the particle physics community pins many of its hopes on discovering a DM candidate that has additional interactions with the Standard Model (SM). The three main reasons for this are simple: many wellmotivated extensions of the SM include DM candidates with such interactions; there are plausible mechanisms that require interactions to provide the correct DM abundance, and importantly, many such mechanisms are testable by experiments. However, the possibility remains that DM only interacts with the SM gravitationally. If this were the case, the production of DM in the early Universe still requires an explanation. One such explanation is the focus of this paper, namely that some population of primordial black holes (PBHs) were abundant and energetic enough to evaporate and produce the relic dark matter we observe today. Notably, such a scenario relies upon particle production via Hawking radiation [1,2], a phenomenon that does not rely on the existence of additional and unobserved interactions. Instead, it arises due to the ambiguity of the definition of * andrew.cheek@uclouvain.be † lucien.heurtier@durham.ac.uk ‡ yfperezg@northwestern.edu § jessica.turner@durham.ac.uk the vacuum state in curved spacetime. The disruption of the spacetime resulting from the collapse of some matter generates a thermal flux of particles. Crucially, a black hole (BH) will emanate all existing degrees of freedom in nature, without regard to their interactions, and thus constitutes a compelling source of a purely gravitationally interacting DM. One of the earliest probes of the Universe's history comes from the cosmic microwave background (CMB) [3,4]. Perhaps the most profound lesson from the CMB is that the observable Universe is remarkably homogeneous. The current scientific consensus is that this is achieved by some early period of cosmic inflation, which also provides the seeds for small matter perturbations that eventually form galaxies. The true model of inflation is far from determined and many of which predict the existence of PBHs. This topic has surged in popularity recently because of the gravitational wave measurements of solar mass black hole binaries. It has been argued that PBHs themselves constitute DM, where their masses are constrained by a large and varied set of experimental searches [5,6]. The minimum value for the PBH mass is set by the requirement that they have not evaporated already, determined by Hawking radiation, M PBH ≥ 5 × 10 14 g [7].
Even without the requirement that PBHs constitute DM, Big Bang Nucleosynthesis (BBN) provides serious restrictions on how many PBHs existed in the early Universe for masses 10 9 g ≤ M PBH ≤ 10 14 g [7][8][9], below which PBHs have evaporated before BBN. A lower limit on the PBH mass comes from constraints on inflation; the Hubble scale during inflation has an upper bound from CMB [10], which in turn imposes the smallest possible mass to be M PBH 0.1 g [7]. Let us stress that such a value is model dependent, specifically on the details of the gravitational collapse and on the features of inflation. One obtains such a minimal value for the PBH mass by assuming a standard slow-roll scenario. For simplicity, we assume such minimal case, and take M PBH 0.1 g [7] as the lower limit. This window keeps alive the possibility that PBHs dominated the early Universe and played an important role in its evolution. The consequences of this have been well studied since the discovery of the Hawking radiation [11], and span many different and important aspects, for instance, the generation of Dark Radiation [12][13][14][15][16][17][18], matter-antimatter asymmetry production [19][20][21][22][23][24][25][26][27], and the implications for the production of DM through evaporation [13,16,17,20,[28][29][30][31][32][33][34][35][36][37][38][39][40]. Generally, DM particles produced in this way can be very light. However, if they are too light, such DM particles are expected to be relativistic and their free-streaming length will be constrained by observations regarding structure formation [16,17,41,42]. This is the first paper of a two-part series, where we return to the calculation of DM emission from PBH evaporation to improve existing treatments. We do so by ameliorating the analysis in two different aspects: solving, in detail, the momentum-averaged Boltzmann equations and including consistently the greybody factors, quantities essential for an accurate description of the Hawking evaporation. The code we use for this purpose has been made publicly available 1 . In addition, we also provide a semi-analytic solution that is consistent with our numerical analysis. Furthermore, we address the possibility of having baroque Dark Sectors, consistent with a large number of degrees of freedom. Since PBH evaporation would produce significant quantities of particles belonging to such sector, one could imagine that, in the scenario, all but one particles are unstable, the generation of the stable DM would be enhanced by such indirect production. In this paper, we assume that this Dark Sector is disconnected from the SM, avoiding thermal production mechanisms such as Freeze-In (FI) or Freeze-Out (FO). In the companion paper [43], we will consider the situation where there are interactions with the SM sector. We use the infrastructure of ULYSSES [44], a publicly available python package that has been typically used to solve Boltzmann equations associated with leptogenesis, to solve the relevant Friedmann and Boltzmann equations.
This paper is organized as follows. First, we describe the emission properties of non-rotating (Schwarzschild) and rotating (Kerr) Black Holes in Sec. II. In each case, we consider the mass and angular momentum loss rate from the BH, the rate of particle emission, and, when 1 https://github.com/earlyuniverse/ulysses possible, the total number of emitted particles. These characteristics will be crucial for the analysis in the subsequent section. Also, we consider the phase-space distribution of emitted particles, which will be helpful to address free-streaming constraints on DM. In Sec. III, we first establish the Friedmann and Boltzmann equations that we solve in the presence of evaporating PBHs. Then, we describe our results for the cases in which the PBHsboth for Schwarzschild and Kerr -are the only source of DM. We then focus on the next-to-minimal case which consists of a dark sector containing only DM together with one heavy metastable state. Finally, we make our concluding remarks in Sec. IV. We have included two appendices: App. A provides useful formulae related to the BH evaporation properties and derive some specific quantities used in the main text, and App. B, which contains the decay width of scalars, vectors and massive tensors into a fermion-antifermion pair. We use natural units where = c = k B = 1 throughout this manuscript.

II. BLACK HOLE EVAPORATION
Black holes were initially thought to be eternal and were expected to increase their mass by accreating additional matter or even other black holes. Nevertheless, when the BH quantum properties were inspected, it was shown that they also emit particles with a thermal spectrum related to BH surface gravity [1,2], making the BHs lose mass and angular momentum in the process. Hence, the properties of the emitted particles depend only on the specific characteristics of the BH, which, according to the no-hair conjecture, are its mass, angular momentum, and charge. We focus here on two distinct cases, Schwarzschild (non-rotating) and Kerr (rotating) PBHs. Next, we discuss the emission properties and the BH evaporation rates for each case separately.

A. Schwarzschild Black Holes
Schwarzschild BHs correspond to the simplest scenario, where the BHs are solely described by their mass, M BH . As Hawking demonstrated in his seminal papers [1,2], the emitted particles from the evaporation process have a thermal spectrum with temperature related to the mass as (G the gravitational constant) The emission rate of any particle species i of mass µ i , spin s i , and number of degrees of freedom g i from the evaporation of a BH, within time dt and momentum [p, p + dp] interval, is given by where E i (p) = µ 2 i + p 2 , and σ si stands for the absorption cross-section. From this emission rate, we will be able to obtain the time evolution of the BH mass and the phase-space distribution of the different particles evaporated. The absorption cross-section σ si -or the related greybody factor, Γ si ≡ σ si p 2 /π -is a crucial characteristic of the Hawking emission rate as it describes the possible back-scattering of particles due to the gravitational or centrifugal potentials [1,2,45,46]. We note that in the literature this factor is sometimes neglected. However, recent works such as Refs. [17,42] provide the most comprehensive inclusion of these greybody factors. Here, in a similar fashion, we include these factors as consistently as possible, given the results in the literature. For instance, we incorporate the the absorption cross-section for massive fermions emitted from Schwarzschild BHs, obtained in Refs. [47,48]. For massive bosons, we will only include the cross-section obtained by assuming a massless field [45]. Since particle emission is only possible when E i ≥ µ i , while the correction to the greybody factors due to the finite mass is not large for such values of energy [49], we do not expect a significant effect from such an approximation. For values GM BH p 1, and independently of the particle's spin, the greybody factors tend to the geometrical-optics limit, σ si (E, µ)| GO = 27πG 2 M 2 BH [45,46,49,50]. Hence, it is convenient to define the ratio of the full greybody factors to the geometrical-optics limit 2 [51] ( In Fig. 1 we present the reduced greybody factors, ψ si (E), for the case of massless particles and different spins, s i = 0 (emerald), s i = 1/2 (purple), s i = 1 (orange), s i = 2 (light blue), as function of E/T BH . The oscillations present in such quantities are related to the different contributions of the partial waves, each having a different value of the total angular momentum quantum number. Moreover, we observe that the low energy contributions are suppressed from higher particle spin values. This crucial characteristic will play an important role in the accurate determination of the relic abundance.
BHs lose their mass over time because of the evaporation process. The reduction in mass can be obtained by summing Eq. (2) over the different species and integrating over the phase space, to obtain [52,53] Fig. 1. Ratio of the greybody factors to the geometric optics limit for massless particles and different spins, si = 0 (emerald), si = 1/2 (purple), si = 1 (orange), si = 2 (light blue), as function of E/TBH.

1.2
where M p = G −1/2 denotes the Planck mass. Here, we have defined ε(M BH ) as the evaporation function which is dependent on the BH instantaneous mass, with the functions ε i (z i ) given by where the integration is performed over the dimensionless parameter x = E i /T BH , and z i = µ i /T BH . The spin-dependent expressions of ε i (z i ) for massless particles, in the geometrical-optics limit, and a fitted form obtained after integrating over the full greybody factors are explicitly given in the App. A. In Fig. 2 we present the different contributions to the evaporation function for particles with different spins, together with the results in the geometrical-optics limit as function of z i . As we observe in this figure, the geometrical-optics limits closely resembles the expected evaporation function for scalars where for bosons with non-zero spin, the approximated forms overestimate the mass loss rate, while for fermions there is a underestimation when z i 4.
Let us determine the momentum-integrated emission rate, Γ BH→i , and the total number of emitted particles per BH, N i . Integrating the Hawking rate, Eq. (2), over the momentum, we obtain Γ BH→i ≡ dp d 2 N i dp dt , where In the massless case µ i = 0, Ψ simply takes a numerical value which depends on the particle's spin [51] We provide useful analytic expressions for Ψ i (z i ) in App. A. The total number of emitted particles of the species i during the BH existence is simply computed by integrating the total rate over time, where τ is the BH lifetime, and with z in i = µ i /T in BH the ratio of the particle's mass to the initial BH temperature, and m j ≡ µ j /µ i the ratio of each existing particle mass to the mass of the species i. The derivation of η i (z in i ) is presented in App. A. Differently from what has been previously done in the literature, we have not assumed any relation between the particle mass and the BH temperature. Instead, N i is general: the Boltzmann suppression present when T BH < µ i is automatically included in it. Let us compare the total number of emitted particles including the greybody factors to the geometric optics limit, we therefore observe that by not including correctly the greybody factors, there is a significant overestimation of the number of produced particles by a BH.

B. Kerr Black Holes
Another possibility is that the evaporating BHs have some non-zero angular momentum. Such rotating BHs, also known as Kerr BHs, could have formed with some  initial spin or acquired their angular momenta via some specific mechanisms, such as mergers [54][55][56]. The BH temperature for the Kerr scenario is modified due to the presence of the angular momentum, where the dimensionless parameter a is related to J, the BH angular momentum, as a = JM 2 p /M 2 . Such parameter can have values a ∈ [0, 1], so that for the case of close-to-maximally rotating BHs, the temperature tends to be zero.
The spectra of emitted particles have an additional dependence on the BH angular momentum, with where Ω = (a /2GM BH )(1/(1 + 1 − a 2 )) is the angular velocity of the horizon and l, m the total and axial angular momentum quantum numbers, respectively. From the emission rate in Eq. (13) it is clear that the absorption cross-section also depends on a . In what follows we will use the procedure established in Refs. [57][58][59] in order to compute the cross-sections σ lm si appearing in Eq. (14) in the case of scalar, fermion, and vector particles in the Kerr scenario 3 . For the spin-2 case [61] we use the greybody factors from BlackHawk as a numerical input when using Eq. (14). Interestingly, the emission of higher-spin particles is enhanced for BHs with a non-zero angular momentum. Thus, we could expect an enhanced emission of spin-2 DM particles, such that it would be possible to increase the relic density. This will be explored in more detail in the next section.
Similarly to the mass depletion in the Schwarzschild case, for Kerr black holes the angular momentum decreases in time because of particle emission. The equation for the angular momentum is obtained by integrating the rate multiplied by the axial angular momentum number in Eq. (14) [45], with γ(M BH , a ) = i γ i (M BH , a ) the angular momentum evaporation function. Substituting the definition of a in Eq. (14), one finds the evolution equations as function of time for both mass and spin, The functions, γ i (M BH , a ) and ε i (M BH , a ), for the different spins can be parametrized in a similar fashion as in the Schwarzschild case, The previous definitions were chosen in order to have a smooth transition to the Schwarzschild case when a → 0. We have determined fitted forms for these factors from explicit integration of the greybody factors in the Kerr similar to those contained in the code BlackHawk [60], finding an agreement at the levels of ∼ 1.37% (∼ 0.44%) for massless scalars, ∼ 1.39% (∼ 10%) for massless fermions, and ∼ 0.55% (∼ 1.8%) for massless vectors in the case of a = 0 (a = 0.99).
case, see App. A. We parametrize the emission rate for spinning BHs similarly to the Schwarzschild case, where, analogously, we have Obtaining a closed form for the total number of particles in the Kerr case is not straightforward. It is not possible to take as an independent variable the BH mass, as done in the Schwarzschild case since the angular momentum also changes with time.
Finally, note that in the limit of an initial a = 0, one readily recovers the Schwarzschild functions. Thus, in our simulations, we solve the Eq. (16) in the cosmological context and impose a = 0 as an initial condition when analyzing the specific scenario of Schwarzschild BHs.

C. Phase-space Distribution of Evaporated Particles
The phase-space distribution of particles emitted from BHs has a significant impact on the evolution of the Universe. For the simple setup explored in this study, the mean free path of DM is the quantity of most consequence, limiting the formation of small-scale structures.
The mean free path of the emitted particles strongly depends on the evolution of their respective phase-space distributions. In the usual FO and FI cases, such distributions are dictated by the Boltzmann distributions already present in the thermal bath. In the presence of BH evaporation, such phase-space distributions may be significantly distorted. Indeed, when they evaporate, BHs produce particles with a typical momentum p(t) ∼ T BH (t). Because T BH is an increasing function of time when BHs evaporate, the momentum of the particles they produce is directly related to the dynamics of the Hawking evaporation. For a particle of mass µ i , this typically leads to two major production regimes: • µ i T in BH : most of the particles produced via evaporation are relativistic, as they carry a momentum p T in BH .
• µ i T in BH : the production is statistically suppressed until the BH temperature increases above the particle mass. Therefore, most of the production occurs when T BH ∼ µ i producing a population of non-relativistic evaporated products.
Given the expression of evaporation rate per unit of time and momentum in Eq. (2), we can derive the phase-space distribution of the different particles produced through 5 10 15 20 25 p/TBH and T BH = m DM /10 (right). We compare the distribution from the full calculation (violet), with the Boltzmann distribution (green dashed) and the Geometrical-Optics limit (blue dashed). We also indicate the average momentum, p (grey dashed) as calculated in (24).
BHs evaporation. In Ref. [41] such a distribution was derived in the geometrical-optics limit in the case where the DM mass, m DM , verifies m DM T BH . Note, however, that in Refs. [17,42], the phase-space distribution was first computed including the greybody factors, showing the crucial impact of incorporating such factors correctly. We also go beyond the geometrical-optics limit and solve those phase-space distributions using our expressions for the greybody factors by simply integrating Eq. (2) over time 4 Extending the results of Ref. [41] to the massive DM case, we can compare our results to the geometricaloptics limit of such a formula where the function f si is an integral that can be com- 4 In principle, taking into account the expansion of the Universe during the evaporation process may slightly alter this result. However, it was shown in Ref. [41] that such an effect is negligible.
puted analytically and Li n are the polylog functions of order n and i = (−1) 2si . In Fig. 3 we depict the phase-space distribution of a fermionic DM particle produced by evaporation in two representative cases where m DM T BH (left panel) and m DM T BH (right panel). We indicate in violet the phase-space distribution of DM particles that we obtain using the full greybody factors in Eq. (2). As expected, such a distribution is peaked around the BH temperature, similarly to what was obtained in Ref. [41]. We compare our results with the distribution of Eq. (21) obtained in the geometrical-optics limit and find that our distribution is slightly shifted towards larger values of the momenta. Such a shift is related to the suppression of the low momenta present in the greybody factors, similar to what was observed in Ref. [42]. We also indicate (dashed green line) the corresponding Boltzmann distribution evaluated at the temperature T BH as well as the value of the typical momentum of evaporated particles (grey dashed line). In the right panel of Fig. 3 one can see that the DM phase-space distribution instead peaks at p ∼ m DM , since BHs mainly produce DM particles after their temperature rises above m DM . Again we can notice a significant shift between our findings and the geometrical-optics limit obtained using the prescription of Ref. [41]. Finally, the authors of Ref. [41] evaluated the Boltzmann distribution at ∼ 3T BH to make the distribution peaks match. We can see that such a prescription must be modified to match a Boltzmann distribution with the full distribution we obtained because of the aforementioned shift towards larger momenta.
An important constraint that the purely gravitational production via Hawking evaporation is subject to corresponds to the warm DM bound. From our discussion above, we have found that the emitted particles could have a large average momenta depending on their masses. In such a case, the redshift resulting from the expansion of the Universe might not be large enough to make the DM non-relativistic at the moment of structure formation, hence contradicting observations. Following previous treatments [16,17,41], we compute the average DM velocity today v 0 from the expected average momentum, with a ev (a 0 ) the scale factors at evaporation (today). We will impose that such a velocity should be smaller than the maximum value allowed from Lyman-α constraints, assuming all DM coming from PBH evaporation, to have a sufficiently cold DM [41,[62][63][64].
The average momentum of an emitted particle will be computed for spinning BHs in a simple and general manner. Reversing the integration order, that is, integrating the Hawking rate first over the momentum and using the definitions of the evaporation functions, Eq. (16), and the momentum integrated Hawking rates, Eq. (18), and then integrating over time, we have This complementary approach will be used in our numerical procedure to enforce the warm DM constrain in our results. It it worth noting that more accurate determinations of the WDM constraint have been undertaken by the authors of Refs. [17,41,42] where the DM phase space has been used as an input to the cosmic linear perturbation solver CLASS [65][66][67].

III. PRODUCTION OF DARK MATTER VIA PRIMORDIAL BLACK HOLE EVAPORATION
Several mechanisms lead to the formation of PBHs in the early Universe after inflation [7,12,68]. For simplicity, we assume that a population of PBHs was formed with a monochromatic mass spectrum. Let us stress that assuming such a simple spectrum allows us to give more generic statements, since it decouples the particle production from the details of the PBH formation. Clearly, the PBH mass spectrum obtained from a given mechanism will depend on specific parameters related to the model. For instance, PBHs formed from collapse of inhomogeneties relies upon the critical value of the overdensities that enter the horizon [7,12]. Other specific models, such as collapses from multi-field inflatons, cosmic strings, bubble collisions, domain walls, or even the PBH formation in an early matter dominated era, will produce distinct mass spectra. Note, however, that the assumption of a monochromatic spectrum is not totally unrealistic, as the PBHs could have formed at very specific time, thus having a rather narrow spectrum. We leave the extension of our results to more realistic mass distributions for future work. Moreover, we consider that the initial PBH mass is proportional to the particle horizon mass at the moment of formation in a radiationdominated era [12] where γ is a factor related to the gravitational collapse, assumed here to be equal to ( The initial PBH population is characterized by the initial fraction of the PBH energy density, ρ in PBH , with respect to the total energy density ρ in , which can be expressed through the parameter β ≡ ρ in PBH /ρ in , or, more commonly, using the definition where T in is the plasma temperature at the time of the PBH formation, and the additional factors are included as the initial PBH fraction always appears corrected by them [12]. Since the PBH energy density scales as a −3 , it is possible to have a PBH-dominated era depending on the initial value of β . Such a possibility will play an important role when we consider the effects of the evaporation on the DM production. Furthermore, for the case of Kerr PBHs, we assume a monochromatic angular momentum distribution, similarly to the mass, such that all BHs had the same initial value of the angular momentum. As with the mass spectrum, this simplification allows us to remain moderately independent of the PBHs formation mechanisms. For the specific case of Kerr BHs, PBHs can acquire a non zero spin via different models, such as mergers, accretion, or even the formation mechanism mentioned before could produce BHs with some spin (see, e.g. [69,70], and Ref. [31] for a first computation of ∆N eff from Dark Radiation considering an extended spin distribution). Therefore, the early Universe will be comprised of three different energy density components, the PBH population plus radiation related to the SM and, possibly, a Dark Sector (DS). The Hubble parameter, therefore, should take into account these three elementary contributions, By means of Hawking evaporation, PBHs will not only change the evolution of the Universe but also emit a large number of particles, regardless of their possible interactions. The set of produced particles will affect the Universe's energy budget and, as we have mentioned before, could lead to the generation of the observed DM. The capacity of PBHs to produce DM particles when they evaporate strongly depends on two factors: (i) whether the temperature of the black holes is smaller or larger than the DM mass, and (ii) whether PBHs evaporate in a matter or radiation dominated era [13, 16, 17, 34-36, 38, 39]. In order to track effectively the number of DM particles produced by a PBH population in the early Universe, we must specify how the phase-space distribution of such states changes over time. We define for the species i 5 where n BH is the PBH number density. Hence, it is possible to write a Boltzmann equation for such a species in a FLRW Universe, where we have included possible interactions via a collision term C[f i ]. In the following, however, we assume that the DM does not interact with the SM thermal plasma, so that such a collision term will be absent. We can obtain the usual equation for number densities after integrating over the phase space, The Friedmann equations for the ρ PBH , ρ SM PBH and SM radiation energy densities, respectively, are given bẏ where the energy produced by the evaporation depends on the mass loss rate since 5 We include in the definition the factor of p 2 /(2π 2 ) because the integration of the Hawking rate over momentum and time directly gives the total number of emitted particles.
where we used that ρ PBH = M BH n BH . The set of Friedmann equations includes two different effects related to the presence of a PBH population. First, PBHs behave as matter, ρ PBH ∝ a −3 , enabling the possibility of early matter domination, as mentioned above. Second, the evaporation produces SM particles that reheat the Universe. Thus, to determine the DM generation consistently, we solve the system of equations, Eq. (31), together with the mass and angular momentum PBH loss rates, Eq. (16), and an equation for the DM number density in the lines of Eq. (30). The solution is found using the ULYSSES python package [44], which allows for a rapid determination of the resulting DM relic abundance, including the PBH evaporation. Some words are in order about the numerical procedure. In general, it is not possible to naïvely apply a differential equation solver to the full system of equations, especially when the DM mass is much larger than the initial PBH temperature because of the stiffness present in the mass loss rate. Such stiffness is a consequence of the explosive nature of the particle emission in the final stages of the BH lifetime. Starting with a relatively large PBH mass, M in BH 1 g, it is not possible to reach M p , a value which we aim to attain when we solve the equations, with direct use of a numerical solver. Instead, we use a zoom-in procedure: We iteratively solve the Boltzmann equation on smaller and smaller time scales until the PBH mass reaches the Planck mass, M p . We have checked that the solutions are stable and correctly account for the case when the particle emission only occurs during the final moments of the BH existence.
Once our coupled equations have reached a stable point, where the Universe is radiation dominated and there is no longer any production of DM, we can use the temperature at which the evaporation occurs, T ev and entropy conservation to obtain today's dark matter density parameter (T 0 is the present temperature), We present in Fig 4 prototypical solutions of the Friedmann-Boltzmann equations for m DM = 0.1 GeV, β = 10 −7 and M in BH = 10 6 g Schwarzschild (left) and Kerr (right) PBHs. The time evolution of ρ i a 3 is displayed for the SM, PBH and DM energy densities. The value of the relic abundance is also shown. We observe in both cases, PBHs modify the evolution of the Universe and generate DM. After a radiation-dominated phase, the PBH density, in this case, leads to an early matter dominated era, which ends when the PBHs evaporate. During the final states of the evaporation, a large entropy injection into the SM takes place, while DM production is accelerated. Such entropy injection is modified if the PBHs had a non-zero a . We will return to these solutions in more detail in the next subsections.
Similarly to our semi-analytic expression in Eq. (9), for the total number of DM particles produced per Schwarzchild BH, N DM , we can obtain the same parameter, where n ev BH is the BH number density at the evaporation, which for a monochromatic mass spectrum can be related to the initial number density n in BH by n ev BH (a ev ) 3 = n in BH (a in ) 3 and thus In general, it is difficult to get a good approximation for all the above values at evaporation (see, however, [43]). Nevertheless, in the case where the populations of PBHs remain a negligible component of the Universe's energy density, entropy conservation can be assumed, leading to the simpler form of the relic density which is fully calculable using Eq. (9) and the initial conditions, Eq. (25) -(26), leading to Where we apply this method, we find agreement with the fully numerical method to the level below 1%. In the case where PBHs play a much greater role in the cosmological evolution, we use the approximations in [43] and obtain values that agree up to some O(1) multiplicative factor. This gives us a high degree of confidence in the accuracy of our calculation. By numerically solving the Boltzmann equations and including the greybody factors as accu-rately as possible, we believe that this work constitutes a step forward in the work connecting DM production and PBHs. As previously mentioned, Refs. [17,42] take great care in consistently using the greybody factors but use approximate analytic solutions to obtain Ω DM h 2 . These approximations are most appropriate when the PBH population does not affect the thermal history of the universe as seen in our validation of the code. Moving to the nu-merical framework for solving these systems allows for a more sophisticated analysis to be performed, where dark sectors for have non-gravitational interactions with the Standard Model. Next, we describe our results regarding the DM production from Schwarzschild and Kerr PBHs, and then we analyze the effects of having a baroque dark sector composed of a large number of particles, whose lightest particle is stable and thus constitutes the perfect candidate to be the DM present in the Universe.

A. Direct Production
In the case where PBHs are the only source of DM, the values of β and M in BH leading to the correct relic abundance are indicated in Fig. 5 for various values of the DM mass. For any of those masses, a point above the corresponding coloured contours leads to an overproduction of DM (Ωh 2 > 0.11) while DM is underproduced in points below the coloured contour.
In the limit where M in BH → 0, the Hawking temperature T in BH ∝ (M in BH ) −1 is always larger than the DM mass. Therefore PBHs produce DM particles during the entire evaporation process. In that limit, the relic density of DM particles produced from of evaporation is linearly related to the fraction of PBHs β. A too-large value of this fraction leads to an overabundance of DM, which sets an upper bound on β. For larger PBH masses, T in BH might be smaller than the DM mass while PBHs still evaporate during a radiation-dominated era (this is typically the case for DM masses above 10 9 GeV). In that case, the larger M in BH , the fewer DM particles are being produced during evaporation, which explains why the relic density contours go up after crossing the T in BH = m DM line in Fig. 5. For even larger M in BH , PBHs dominate the universe energy density before they evaporate and reheat the SM bath at a temperature T ev . This is the case if their energy fraction β at the time of PBH formation T = T in is larger than β c ≡ T ev /T in . In that case, the relic abundance of DM particles does not depend on the PBH fraction anymore but rather only on the PBH mass, this is reflected by the contours being vertical past the line β = β c . Interestingly, on the right of those vertical lines, PBHs can significantly reheat the Universe, and therefore modify the evolution of the SM thermal bath while not overproducing DM particles. Note that in most of the previous works, the contours depicted in Fig. 5 were derived analytically, ignoring the greybody factors [36] and/or fully tracking the Boltzmann equations, we indicate such a result with dashed lines. Our studies used the evaporation rates, including the full greybody factors, leading to significantly shifted contours towards larger PBH masses (plain coloured lines), assuming the DM to be fermionic. Since the Hawking rate departs from being a full blackbody spectrum because of the absorption probabilities, the number of emitted particles is larger than expected in the approximated purely-Planckian form. Moreover, Ref. [36], whereas the plain lines were derived in this work, including the full greybody factors. We assume the DM to be a fermion.
the evaporation temperature is greater when including the greybody factors. Thus, we observe that smaller values of β are required to give the correct relic abundance. Let us notice that our results coincide qualitatively with those from Ref. [42] for the case of a PBH dominated Universe, when fully including the greybody factors for the different spins. For a Universe where there was not PBH domination, but the BHs constitute an important contribution to the energy budget of the Universe, our results differ from those in Ref. [42]. Such a difference arises because in our code we always include the PBH contribution to the evolution of the Universe, which can alter the final relic density.
Yet another effect of including the greybody factors correctly is that the relic abundance depends on the spin of the DM particle. In the top panel of Fig. 6, we show such dependence for several values of the DM mass and for spin s i = 0, 1/2, 1, 2. For lighter DM masses relative to the initial BH mass, T in BH m DM , we observe that larger values of β are required to produce the correct Ωh 2 for larger spins. Similar conclusions were found in Refs. [34,42] and it is a direct result of the suppression of the number of emitted particles for higher spins due to the greybody factors, see Eq. (8) and Fig. 1. On the other hand, when T in BH m DM , the cutoff induced by the non-zero mass affects the scalar case especially, such that the number of emitted particles is reduced in comparison to the case where the PBH temperature is higher than the DM mass. Hence, the required value of β necessary to obtain the observed Ωh 2 is larger and becomes similar to the values needed in the case of a Vector DM. As specified in the previous section, we have applied the warm dark matter constraint in the same figure. We indicate where the DM would violate the cold dark matter condition by marking the line with yellow. From this, we observe that for light masses, m DM 10 GeV for M in BH 10 4 g, the DM particles emitted from the evaporation are too hot, thus in tension with the observations [16,36,41]. We find that our method proves to be more conservative than is reported in Refs. [16,36,41,42]. It is true that our procedure described in Sec. II C is certainly more approximate than that of Refs. [41,42] and we leave the implementation of such methods in our code to future work. The warm dark matter constraint becomes less relevant for heavier masses, and for m DM 100 GeV, the full parameter space would obey such a limit.

B. Effect of the BH spin
As mentioned in Section II B, Kerr PBHs could have a unique impact on the DM generation given the peculiar features present in such a case. Specifically, the enhanced emission of spin-2 particles that can compensate for the large initial fractions required to account for all the DM. In Fig. 6, bottom panel, we present the energy fraction β as a function of the initial PBH mass for the Kerr case, assuming a value of a = 0.99. Interestingly, we observe that, when T in BH m DM is valid in all the parameter space, the values of β that give the correct relic abundance coincide for scalars, fermions and vectors. Such agreement is related to the increase of high-spin emission reflected in the greybody factors. Moreover, the initial PBH fraction that gives the correct relic abundance for spin-2 DM is reduced by ∼ 2 orders of magnitude with respect to the Schwarzschild case. For the case where T in BH m DM , a similar behavior to the non-rotating case is present; the emission cutoff due to the DM mass diminishes the overall particle production, specifically for scalars and fermions. Nevertheless, for tensor DM, there is an interesting effect when T in BH m DM . Even though a large Boltzmann suppression is still present, the enhanced emission of tensor particles [71] is a significant countervailing effect, which leads to enhanced particle production. In Fig. 6, such amplification is responsible for the structure observable for m DM = 10 7 GeV (10 15 GeV) for BH masses of M in BH ∼ 1 − 10 g (10 8 − 10 9 g). Such an effect is also present for scalars, fermions, and vectors, although much less conspicuously 6 . Ref. [17] also investigated the affect of Kerr BHs on DM production, focusing predominantly on light m DM , however one can observe the amplification of Tensor particle production at high a in their Fig. 10. Finally, regarding the warm DM constraint, we observe that the BH spin increases the parameter space that is excluded by such a limit in comparison to the non-rotating case, this is also in agreement with Ref. [17]. Still, we have demonstrated that the DM production from Kerr BHs has many compelling features not encountered before.

C. Indirect Production: Presence of additional dark sector particles
The DM could be part of a much larger dark sector, containing a large quantity of particles. Such a baroque scenario should not be inconceivable from what we have learnt about the SM sector. Indeed, supersymmetric (SUSY) models constitute the perfect example of UV complete scenarios that are expected to contain many additional degrees of freedom [30]. Let us assume that the DM particle belongs to an extended sector that does not interact with the SM. Moreover, for simplicity, let us consider that just one particle is stable, just like the lightest superpartner in SUSY with some R-parity. Suppose that there are i copies of X particles, where X = {S, F, V, G} indicates whether the particles are scalars, fermions, vectors or tensors, respectively. The total number of final DM particles produced via PBH evaporation N tot DM will be the sum of all emitted particles being n Xi→DM the number of DM particles resulting from the decay of X i , such that n Xi→DM ≥ 2. Following our previous analytical estimation of the final relic abundance, we can examine the enhancement of Ωh 2 with respect to the case where there is just the DM particles, where T t ev , a t ev are the Universe temperature and scale factor at evaporation in the extended dark sector case. From this, we observe that the effect of having additional dark sector particles is twofold. First, the increase of DM particles evidently enlarges the final relic density. Second, since the emission of the additional particles affects the BH lifetime, the Universe properties when the PBHs evaporate are changed, and thus Ωh 2 .
Let us be more specific and consider the situation in which the dark sector is only composed by the DM particle and a heavier state X, assuming for simplicity one decay channel, X → DM + DM, such that n X→DM = 2. In order to be consistent in our treatment, we solve the same set of Eqs. (31), plus the following equations for the number density of X and DM including the exchange termṡ with the thermally averaged decay width of X given by where "ev" indicates that the average is taken with respect to the BH temperature, and Γ X→DM the decay width of X in vacuum 7 (for further details, see the com- panion paper [43]).
In Fig. 7, we present Ωh 2 for a Fermionic DM as a function of the mass of X. We show the result for different types of X, scalar (emerald), vector (light blue), and massive tensor (orange), for m DM = 10 5 GeV, β = 10 −17.75 , and M in BH = 10 6 g. Without accounting for greybody factors, one could expect that Ωh 2 should increase by a factor of 3 since the X would decay into two DM particles. However, the more accurate calculation leads to enhancements of ∼ 3.7, ∼ 1.9 and ∼ 1.3 for a scalar, vector or tensor X respectively. Once more, we are seeing the greybody factors affect the emission of higher spin particles more significantly, reducing the contribution of X to the total. We note that M in BH = 10 6 g means that T in BH ∼ 10 7 GeV, suggesting that the suppression of X particle emission should occur when m X 10 7 GeV. This is corroborated by Fig. 7 and we see that by m X ∼ 2 × 10 8 GeV enhancements of Ωh 2 from X decay is negligible. Of course this suppression is independent of the particle's spin, so we see the behaviour across the three cases in the figure.
Interestingly, once there is a large enough separation of scales between X and DM, the warm DM bounds need to be considered once more. Unlike previously, we now can have a mix of cold and warm DM from Hawking emission and dark sector decay respectively. According to Ref. [42,64] when the fraction of warm/cold DM is less than ∼ 0.2, constraints from structure formation do not apply. Despite this, even when the parent dark sector particles do not contribute much, such as the X = G, Schwarzschild PBHs case, the fraction of DM particles which could be warm is sufficiently large ∼ 0.3. This has implications for the coexistence of dark sectors that contain a dark matter candidate as its lightest member and PBHs. For the decay products, the average momentum is given by where p X is calculated following Eq. (24) but for the parent particle. Plugging Eq. (42) into Eq. (23) gives the velocity today which can be compared with Lyman-α constraints [41,64]. By taking the setup in Fig. 7, we find that m X would have to be above 10 10 GeV to produce warm dark matter with mass 10 5 GeV. At which point, the constraint would be irrelevant because the contribution decayed DM has on the relic abundance is negligible. The scale separation that leads to a heavy warm dark matter component is highly dependent on the time of BH evaporation, the later evaporation occurs the less the particles will red-shift. For example, taking the Fig.  7 setup but M in PBH = 10 8 g, now m X ≥ 10 8 GeV produces warm DM when m DM = 10 5 GeV.
If the PBH population had an initial non-zero angular momentum, we find that the spin of X plays a crucial role in the final Ωh 2 . We present the relic abundance as a function of m X in Fig. 8 for three different values of a = {0, 0.5, 0.99999} corresponding to full, dashed and dotted lines, respectively, assuming the heavier state to be a vector and considering the same parameters as in Fig. 7. We find two different effects at play here; when the particle X is kinematically accessible by the evaporation, the indirect DM production is largely enhanced because of the PBH spin. Meanwhile, the relic density is decreased when m X T in BH in comparison to the Schwarzschild case since Kerr PBHs inject much more entropy to the early Universe due to the amplified production of SM boson states. From such effects, we have that the increase in the final relic density is ∼ {1.9, 2.0, 2.2} for a = {0, 0.5, 0.99999} with reference to the Schwarzschild value without any additional state, respectively. Such an augmentation is more stringent if the X particle has a spin of 2, reaching a value of ∼ 4.3 for a = 0.99999. Thus, one can see how having a rich dark sector at high masses can quickly overclose the Universe even if there is a tiny number of PBHs in the early Universe.

IV. CONCLUSIONS
Black holes are one of the most fascinating objects predicted by General Relativity. Initially thought to be everlasting, we have learnt that they instead evaporate by emitting a thermal flux of particles, losing simultaneously their mass and angular momentum. Such evaporation in the early Universe could have critical consequences on our understanding of how the Universe came to be what we observe. In particular, since the Hawking radiation is democratic in nature, i. e., BHs emit all existing degrees of freedom in nature, the observed relic abundance could be the result of the evaporation of PBHs, even in the case that the DM only interacts gravitationally. In this paper, we have addressed distinct effects that impact the DM production by PBHs in the purely gravitationally interacting scenario thoroughly. We have solved the system of Friedmann -Boltzmann equations, and investigated systematically the distinct features present in this scenario for both Schwarzschild and Kerr PBHs. Especially, after including consistently the greybody factors in the description, we have demonstrated how the DM relic abundance can depend on the particle's spin, in such a way that the initial PBH fraction necessary to obtain the observed values is ∼ 2 orders of magnitude larger for massive tensors than for scalar DM. Besides, by correctly including the mass cutoff due to Boltzmann suppression, we have identified the modifications of the required fractions when the initial PBH temperature is smaller than the DM mass. In such a case, the emission only occurs in the last stages of the BH lifetime. Regarding the warm DM bounds that affect this scenario, we have computed the average momenta of the emitted particles, finding it to be larger than estimated before because of the energy dependence present in the greybody factors. Light DM masses are thus in tension with small scale structure measurements, similarly to the previous results present in the literature. We have also illustrated the properties of DM production in the case that PBHs had an initial non-zero angular momentum. The enhancement of the emission expected for bosons, particularly for spin-2 particles, reduces the initial fractions needed to generate the observed DM. Interestingly, we have identified the regions in the parameter space where the enhanced Hawking emission due to the BH spin plays a significant role in the particle production.
Finally, we have analysed the impact of having a large dark sector containing a unique stable particle, the DM candidate. For such models, PBHs would also emit the additional unstable particles of the dark sector during its evaporation, will would produce an additional surplus of DM particles. Such indirect production alters not only the number of DM particles during the PBH evaporation but also the PBH lifetime and impacts the Universe's evolution via entropy injection. We scrutinized a minimal scenario where there exist just one additional heavier particle that decays into the DM. In this case, we found that the increase on the relic abundance can be as large as a factor ∼ 4 in the case that the heavy particle is a scalar. For other types of spins, the factor is smaller. This dependence on the spins is simply understood as the effect due to the greybody factors. In this regard, we also investigated the indirect mechanisms for Kerr PBHs, finding, as expected, an enhancement by a factor of ∼ 4.3 for the tensor case when the PBHs initially had a closeto-maximal angular momentum. Assuredly, an extended dark sector can lead to a rich phenomenology. Moreover, if we assume the existence of interactions with the SM, there could be significant modifications to the results presented here. Such a treatment is left for the second part of this series [43]. Here we go through the analytic derivation of the emission rates and total number of particles of Schwarzschild BHs including greybody factors. The Hawking spectrum, for a given particle species, i is parametrized as where g i is the internal d.o.f, s i is the spin and ψ si (E) is the absorption cross section normalized to the geomet-ric optics limit, and G and M BH are the gravitational constant and the mass of the BH respectively. Introducing the dimensionless parameters x ≡ E/T BH and z i = µ i /T BH , the total emission rate per particle species is where for now we simply take the integral result unspecified as Ψ i (z i ). An assumption that is often made is that ψ(x, z) = 1, that is, take the greybody factors equal to the geometrical-optics limit, which allows one to perform the integral analytically, being Li n the polylog functions of order n, and i = (−1) 2si . We then have Therefore, under this assumption and taking µ i = 0 this allows one to make a comparison between the calculation with the full greybody factors in the massless limit. We can carry out a similar procedure for the evaporation function ε i (z i ) per particle species, defined by where we now have defined the function ε i (z i ) in a similar fashion to Ψ i (z). Its fairly straightforward to obtain the massless geometric optics limit for ε i (0), for Bosons in agreement with Ref. [41].
To find the total number of emitted particle, we need to integrate over the lifetime, τ of the BH. where we have chosen the time the BHs are formed to be t = t in = 0. Using the mass loss rate eq.(4) we can make a change of variables where ε ≡ i g si ε(M BH ). For z = 0, and taking the Geometric Optics limits, one recovers the results from Refs. [36,41] To keep the greybody factors in the equation we can rewrite ε(M BH ) as (A13) Writing ε(z i ) = j g sj ε j (m j z i ) and using that z i = 8πGM µ i , we obtain where (A15) In order to obtain a semi-analytic approximation all functions, Ψ i (z) and ε i (z) have been fitted to a generalized logistic form where the parameters {A s , B s , C s , ν s } depend on the spin of the particle. We give such parameters in Tab. I, and we present an example for the fit function for Γ BH→i for all particle types together with the values obtained by direct integration in Fig. 9.

Kerr case
For initially rotating BHs we can perform a similar analysis. The Hawking rate is modified because of the presence of the non-zero angular momentum, being ψ lm si (M BH , p, a ) the greybody factor dependent associated to the partial wave with quantum numbers l, m, and normalized to 27π 2 G 2 M 2 BH . The total emission rate, obtained after integration over the energy, is being where x = E i /T S BH and z i = µ i /T S BH , being T S BH = (8πGM BH ) −1 the temperature for Schwarzschild BHs. After considering the greybody factors for each spin type, and integrating numerically, we have performed a fit to Ψ i (z i , a ) in the form The fitting parameters {α s j , β s j , η s j , δ s j }, j = 1, . . . , 5 for each spin are given in Tables II-V. We present an example for the fit function for Γ BH→i for all particle types together with the values obtained by direct integration in Fig. 10.
The mass ε i (z i , a ) and angular momentum γ i (z i , a ) evaporation functions per spin defined by  being J the BH angular momentum. These evaporation functions are parametrized as with x and z i defined as before. Similar to the particle emission rate, we fit these evaporation parameters according to a general logistic form, Eq. (A20), where the parameters are also given in Tables II-V.

Appendix B: Decay Widths
In this appendix we quote the decay widths of scalars, vectors and massive tensors into a fermion-antifermion pair, used in the subsection III C. For X having a mass