Bounds on warm dark matter from Schwarzschild primordial black holes

We consider light dark matter candidates originated from the evaporation of Schwarzschild primordial black holes, with masses in the range $10^{-5}-10^9$ g. These candidates are beyond Standard Model particles with negligible couplings to the other particles, so that they interact only gravitationally. Belonging to the category of warm dark matter, they nevertheless spoil structure formation, with a softer impact for increasing values of the candidate spin. Requiring such candidates to fully account for the observed dark matter, we find that the scenario of black hole domination is ruled out for all spin values up to 2. For the scenario of radiation domination, we derive upper limits on the parameter $\beta$ (the primordial black hole energy density at formation over the radiation one), which are less stringent the higher the candidate spin is.


Introduction
Dark matter (DM) candidates must be beyond Standard Model (SM) particles, neutral and stable. Having so far escaped detection, they must have tiny interactions with SM particles. It would be even possible that they interact only gravitationally.
A possible production mechanism for DM particles, taking place in the early universe, is via evaporation [1] of primordial Black Holes (BHs), with masses in the broad range 10 −5 − 10 9 g. In this case, all particles with mass below the Hawking temperature of the BH are emitted, with weights simply given by their number of degrees of freedom (dof). It has been proposed that the particles produced via the evaporation mechanism might be responsible for the excess of baryons over anti-baryons [2,3], for the observed dark matter abundance [4][5][6] and, if sufficiently light, also for dark radiation [5,[7][8][9]. Apart from the case of gravitino production [10,11], the primordial BH density at formation for the range 10 −5 − 10 9 g is at present unconstrained, as reviewed e.g. in ref. [12]. However, ref. [13] derives an upper bound on the fraction of the universe collapsed into primordial BHs in this very mass range from possible backreaction gravitational waves from primordial BHs.
Depending on the fraction of primordial BHs at formation with respect to radiation β, there is a possibility that the universe was BH dominated before the evanescence of the BHs [4,14,15]: this situation is referred to as BH domination. The case in which the BHs evaporate before they dominate the energy content of the universe is called radiation domination.
Fujita et al. [4] calculated the contribution to DM by primordial BH evaporation into new particles beyond the SM: they found that a significant contribution to DM could come from stable particles that are either superheavy or light, that is with masses in the MeV range. In the light case, DM candidates would be warm, while in the superheavy case they would be cold. Exploiting the warm DM velocity constraints available at that time [16], ref. [4] first discussed also the lower limits on the mass of the light DM candidates, using an order-of-magnitude argument essentially based on the geometrical optics approximation for Hawking radiation. This approximation ignores the low energy suppression in the greybody factors [17,18], accounting quite well for the case in which the warm DM candidate has s = 0, but missing to reproduce the case of different spins. For an up-to-date presentation of this argument see [19,20].
A more sophisticated analysis was done by Lennon et al. [5]. They also adopted the geometrical optics approximation, but included the redshift effect in the calculation of the momentum distribution of the emitted particles. Their result is an estimate of the number of particles that are still relativistic, with a spin dependence reintroduced a posteriori and based on greybody factors derived from the older literature [18,21]. As a rough-and-ready criterion for successful structure formation, they impose that when the temperature of the universe drops below 1 keV (at which stage the horizon mass is about 10 9 solar masses), less than 10% of the DM is relativistic. The result of this ingenious, but quite arbitrary, argument is that, for BH domination, warm DM candidates with s ≤ 1 are excluded, those with s = 3/2 are marginally allowed, while those with s = 2 naively survive. Summarizing, for the lower spin values (say s = 0, 1/2, 1), the order of magnitude results of ref. [4] were confirmed by ref. [5], but the latter analysis was however not fully conclusive for the higher spins (s = 3/2, 2).
The more recent analysis of Baldes et al. [20] goes some step further. As suggested in [5], they include the redshift effect in the momentum distribution of the emitted particles at evaporation, and derive the related phase-space distribution as an input for the Boltzmann code CLASS [22][23][24]. The latter allows to extract the matter power spectrum for warm DM from primordial BHs, and to compare it to the standard cold DM case thanks to the transfer function. This enables to constrain warm DM from primordial BHs using the structure formation bounds from Lyman-data already derived for the well known case of DM thermal relics. The analysis of ref. [20] however relies on the geometrical optics approximation and, in particular, provides quantitative results only for the s = 1/2 case, which agree with previous order-of-magnitude estimates [4,19], also based on the geometrical optics approximation. The case for the higher spins could thus not be quantitatively clarified (apart from a qualitative mention of the greybody effects in appendix A of ref. [20]) with respect to the results of ref. [5].
Given the present lack of robust results about the fate of warm DM candidates with high spin values, we think it would be useful and timely to make a dedicated study. The aim of this work is precisely to provide a complete and updated study on the viability of warm DM candidates from the evaporation of primordial BHs.
In order to numerically account for the greybody factors associated to the different spins, we use the recently developed and publicly available code BlackHawk [25]. We also compare the numerical results from BlackHawk with the analytical ones derived in the geometrical optics approximation. Taking into account the redshift effects as suggested in ref. [5], we study the impact on structure formation by calculating the transfer function with CLASS, as suggested in ref. [20]. We derive the transfer function for all spins values, finding that, assuming BH domination, the scenario of warm DM from primordial BHs is excluded for all spins and for all BH masses in the range 10 −5 − 10 9 g. Our results for the s = 0 case, agree with previous order-of-magnitude estimates [4,19]. For radiation domination, we derive the upper limits on β (or, equivalently, on the warm DM mass) for the various warm DM spins. For the case s = 1/2 (the only for which the comparison is possible), we find conceptual differences with respect to the results of ref. [20], but substantial numerical agreement.
In this work we consider BH evaporation as the only production mechanism. The consequences of allowing for other production mechanisms have been recently explored in refs. [26] and [27,28]. For a mixed model of DM production, ref. [26] proved that a primordial BH dominated period of DM creation by evaporation cannot explain the abundance observed today. For an updated analysis of the possibility that the matter-antimatter asymmetry is due to particles produced by primordial BHs evaporation, we refer the interested reader to ref. [29] for GUT baryogenesis, and to ref. [30] for leptogenesis. DM and baryogenesis in the case of stable remnants from thermal 2-2-holes have been studied in ref. [31].
The paper is organized as follows. In sec. 2 we introduce our notation and review basic ideas about formation and evaporation of primordial BHs. In sec. 3 we discuss the instantaneous primary spectrum for the emitted particles. In sec. 4 we discuss the dynamics of the primordial BH abundance. Sec. 5 deals with the momentum distribution at evaporation, and sec. 6 with the calculation of the DM phase space distribution. The calculation of the DM abundance and the impact on structure formation are presented in secs. 7 and 8 respectively. The discussion of the results and our conclusions are presented in sec. 9.
In order to have a better control of our formulas for dimensional analysis and numerical computations, we do not use natural units.

Preliminaries on primordial BHs
The proposal of the early existence of collapsed objects, later called primordial BHs, dates back to 1967 [32]. The formation of primordial BHs from early universe inhomogeneities was considered in refs. [1,33,34]. However, since inflation removes all pre-existing inhomogeneities, any cosmologically interesting primordial BH density has to be created after inflation. Various mechanisms have been proposed, as for instance: that they formed from large inhomogeneities arising through quantum effects during inflation or that some sort of phase transition may have enhanced primordial BHs formation from primordial inhomogeneities or triggered it. We refer to [12,35] for reviews of these proposals, with proper references to the associated literature. In the following, we present our notation for the early universe dynamics and review the primordial BH formation and evaporation mechanisms.

Radiation dominated era
According to the first Friedmann equation (neglecting the curvature and cosmological constant terms), the early universe evolution is described by where a(t) is the scale factor, H(t) is the Hubble parameter, ρ(t) is the energy density of the universe and G is the Newton gravitational constant, G 6.674×10 −11 m 3 /(kg s 2 ). Instead of G, it might be convenient to rather use the Planck mass, In the early hot and dense universe, it is appropriate to assume an equation of state corresponding to a fluid of radiation (or relativistic particles). During radiation domination, ρ ∝ a −4 , a(t) ∝ t 1/2 , and At relatively late times, non-relativistic matter eventually dominates the energy density over radiation. A pressureless fluid leads to the expected dependence ρ ∝ a −3 , a(t) ∝ t 2/3 , and The radiation energy density (at high temperatures) can be approximated by including only those particles which are in thermal equilibrium and have masses below the temperature T of the radiation bath where k B is the Boltzmann constant, k B 8.617 × 10 −5 eV/K, is the reduced Planck constant, 6.582 × 10 −16 eV s, c is the velocity of light in the vacuum, c 2.998 × 10 8 m/s, and g B(F ) is the number of degrees of freedom (dof) of each boson (fermion).
Below the electron mass, only the photon (g γ = 2) and three light left-handed neutrinos contribute, so that g * (T ) = 7.25. Below the muon mass, also the electron (and the positron) has to be included, so that g * (T ) = 10.75. For the full SM, here defined including three light left-handed neutrinos, g * (T ) = 106.75. Adding to the SM three light right-handed neutrinos (as in the case of neutrinos with Dirac nature or in the case of a low-scale seesaw mechanism), g * (T ) = 112. At higher temperatures, g * (T ) will be model-dependent. Including the massless graviton (g G = 2) has the effect of adding 2 units to the previously mentioned values of g * (T ).

Formation of primordial BHs
As reviewed for instance in ref. [12], if a primordial BH forms at the time t f during the radiation dominated era, typically its mass is close to the value enclosed by the particle horizon near the end of inflation where γ 1 is a numerical factor that depends on the details of the gravitational collapse, ρ R (t f ) and H(t f ) are respectively the radiation density and the Hubble parameter at the formation of the BH, and in the last equality we used eq. (2). Using eq. (1), we can also write where the last lower bound follows from the fact that CMB observations put an upper bound on the Hubble scale during inflation, H I 3 × 10 14 GeV at 95% C.L. [36], and H(t f ) H I . In the literature the value γ = 1/(3 √ 3) ≈ 0.2 is usually taken as reference value [12]; in this case the lower limit would become M BH 0.07 g. For values of γ smaller than 0.2, the lower bound on the BH mass would get accordingly weaker. However, this constraint applies only to conventional inflationary scenarios, i.e. the standard slow-roll models of inflation with Einstein gravity (see e.g. the review [37] and references therein). In more sophisticated scenarios, like e.g. the recent model [38], the scale of inflation can not be determined by CMB observations. In any case, the mass of primordial BHs should be larger than the Planck mass, namely M BH 10 −5 g. As is well known, there is also an upper bound on the abundance of primordial BHs of mass M BH 10 9 g (not a theoretical bound), because of their effects on BBN yields, see e.g. ref. [39] for a recent analysis. The range of primordial BH masses between these bounds is at present generically unconstrained [12].
Recalling eq. (2), the primordial BHs formation time is easily calculated from eq. (6) As for the radiation temperature at formation, combining eqs. (1), (4) and (6), we have The temperature and the time at formation of primordial BHs, as a function of the their mass at formation, are plotted e.g. in fig. 1 of ref. [19].
It is useful to introduce the parameter β defined as the BH energy density over the radiation energy density at the formation time where n BH (t f ) is the primordial BH number density at formation and the last equality holds only for a monochromatic mass distribution.

Evaporation of primordial BHs
Here we review the basic formulas describing the evaporation mechanism, in the case of a Schwarzschild (that is, uncharged and non-rotating) primordial BH.
Consider a Schwarzschild BH of mass M BH (t) (we neglect the time dependence only when we refer to the formation time). Hawking radiation mimics thermal emission from a blackbody with a temperature T BH (t), given by [1] Hereafter we denote by T BH the Hawking temperature at formation, namely T BH = T (t f ). As discussed e.g. in [40], at the time t, such a hole emits particles of type i and spin s i and total energy between (E, E + dE) at a rate, per dof, given by where E 2 = p 2 c 2 + m 2 c 4 , the greybody factor Γ s i is a dimensionless absorption probability for the emitted species (it is in general a function of E, M BH (t) and the particle's internal dof and rest mass), and g i are the internal dof of the i-th particle, which account for spin, polarization and color. For the counting of the internal dof, we follow the notation of [25] (see their table 3).
Let us consider in some detail the SM. For the Higgs boson (s = 0), g h 0 = 1. For the massless (s = 1) photon and the 8 gluons, g γ = 2 and g g = 16. For the massive (s = 1) W ± and Z bosons, As for the fermions (s = 1/2): the charged leptons, being Dirac fermions, have g e = g µ = g τ = 4; the neutrinos have g νe = g νµ = g ντ = 2(4) in the case they are Majorana (Dirac) particles, respectively; the quarks have g u = g c = g t = g d = g s = g b = 12. Finally, one might also include the graviton (s = 2), with g G = 2.
As E → +∞, each species approaches the geometrical optics limit but falls off more quickly as E → 0, with the higher spins producing the stronger cutoffs.
When a massless particle scatters off a non-rotating, uncharged hole, the low-energy (that is GM E/ c 3 1) analytic form of Γ s i , averaged over all orientations of the hole with respect to the spin-weighted spherical harmonics and angular momentum quantum numbers of the incoming field, has been computed in ref. [17] (see also [18]). The non-0 rest mass m DM of the DM particles acts as a cutoff in their emission at low energy, but the precise shape of the greybody factor around this cutoff is not relevant here since k B T BH m DM c 2 at all times.

Rate of mass loss and BH lifetime
The rate of mass loss for an evaporating BH is proportional to the total power emitted To parametrize this, Page [18] introduced the adimensional Page function f (M BH ) such that where the time dependence of M BH is understood.
Defining t ev as the time of the BH evanescence, we practically have τ = t ev .
For the SM, the function f (M BH ) is constant over the range of BH masses we are interested in, namely 10 −5 − 10 9 g, because all SM dof are already radiated by a 10 9 g BH. Its value (not including gravitons) is f (M BH ) = 4.26 × 10 −3 . This is the value that we are going to use in the following.
Assuming f (M BH ) constant over the BH lifetime, the latter is easily found to be which corresponds to the following time dependence of the BH mass, 3 Instantaneous primary spectrum In this paper we consider the possibility that, on top of the SM, a DM candidate of mass m DM is produced in the evaporation process. It is well known that there are two possible solutions, denoted as "light" and "heavy" DM, according to the fact that the particles are produced during all the BH lifetime or just in its final stages (see e.g. [19] and references therein). If m DM c 2 < k B T BH , the DM candidate belongs to the "light" category, otherwise to the "heavy" one. In this paper we focus on light DM.
Eq. (11) gives the instantaneous spectrum of the particles of type i emitted by a single BH. The maximum of the energy distribution is at E ∼ k B T BH . For sufficiently light DM, the ultra relativistic limit, E ≈ pc > m DM c 2 , is thus justified. As the BH evaporates, its temperature increases, and the relativistic limit is satisfied a fortiori. In the relativistic limit, the instantaneous distribution of emitted momentum, per dof, is where we introduced explicitly the time dependence for the sake of clarity and s is the spin of DM.
To obtain the number densities (per dof) of the emitted particles, one has to multiply the above expression by the number density of the BHs n BH (t)

Instantaneous distribution: geometrical optics approximation
Using eq. (18) with the geometrical optics limit, eq. (12), at a fixed time t, the instantaneous distribution (per particle dof) has the form Consider for definiteness the formation time t f . It is useful to introduce the quantity as it allows to get rid of the BH mass dependence. Indeed, the instantaneous distribution becomes In the left plot of fig. 1, we show the instantaneous momentum distribution, as a function of x(t f ), for bosons (B) and fermions (F) respectively, in the geometrical optics approximation. Right: the same using the numerical results from BlackHawk.

Instantaneous distribution: numerical results from BlackHawk
This has to be confronted with the corresponding numerical results obtained using BlackHawk [25], which provides the instantaneous momentum distribution at various times, from formation to evaporation. BlackHawk has been modified to obtain the additional DM spectra, a modification already highlighted in the manual and which will be incorporated in a future version of the public code. In particular, greybody factors for the spin 3/2 case have been computed. For a direct comparison with the geometrical optics approximation, the distributions in the right panel of fig. 1 are normalized per particle dof. We can see that the s = 0 case is quite similar with respect to the geometrical optics limit for the boson, while the s = 1/2 case is a bit suppressed with respect to the geometrical optics limit for the fermion. The higher the spin is, the more the distribution is suppressed at low energies, so that the mean momentum gets higher.
A more precise comparison with the geometrical optics approximation can be established by studying the ratio of the instantaneous distributions at formation for BlackHawk over the corresponding geometrical optics limit, as shown in fig. 2. As is well known (see e.g. [18] and references therein), the ratio for s = 0 and s = 1/2 at low energy is a constant, respectively equal to 16/27 and 2/27; the figure correctly reproduces this low energy behavior.

From formation to evaporation
Let us define f (t) ≡ ρ BH (t)/ρ R (t) ∝ a(t). As time increases, it is then possible that BHs come to dominate the energy content of the universe before they completely evaporate [4,14,15]: this situation is referred to as BH domination. The scenario of evaporation before this happens is referred to as radiation domination.

Radiation domination
The BHs evaporate during the radiation epoch if f (t ev ) 1. We defineβ the maximum value of β corresponding to radiation domination, namely the value of β leading to f (t ev ) 1; this value can be obtained from the following relation where we used eqs. (7) and (16) to obtain the last equality. The value ofβ as a function of the BH mass is shown in fig. 3, taking for definiteness f (M BH ) as in the SM and γ = 0.2. For all the values of β β (red solid line), the primordial BHs evaporate before they come to dominate the energy content of the universe and, also in this case, the increase in the scale factor is in the scale factor during BH domination period is smaller (mild) or greater (strong) than the change during the radiation domination period. Right: g * (t ev ) for radiation domination (clearly, this is an approximation by a step function, while the real function g * (t ev ) is a smooth continuous function).

BH domination
In the case of BH domination, we have instead to consider: first, the radiation dominated period from the formation time, t f , to the time when BHs start to dominate, t BH , such that f (t BH ) 1; and second, the matter dominated period from t BH to t ev . The first period is characterized by the following increase in the scale factor The second period is characterized by Putting together Notice also that the two periods have the same increase in the scale factor, a(tev) , if β =β 4/7 , so that t BH = t f /β 8/7 . We can define as mild and strong BH domination the regions where the increase in the scale factor in the first period is respectively larger and smaller than the second one, as shown in fig. 3. Notice that in the very strong BH domination region (close to β ∼ 1) one can neglect the first period of radiation domination, having

Radiation temperature at evaporation
For radiation domination, combining eqs. (1) and (2), we have Using also eqs. (4) and (16), we obtain The values of g * (t ev ) as a function of the BH mass are shown in the right panel of fig. 3.
For full BH domination, we can grossly assume that all the energy stored in the BH density goes, after their evaporation, into the radiation energy density of the (SM and possibly beyond SM) particles emitted by the BHs. These particles rapidly thermalize as soon as they are emitted, so that ρ R (t + ev ) ≈ ρ BH (t − ev ). Combining eqs. (1) and (3), we have We then have For BH domination, the radiation temperature after evaporation gets slightly enhanced with respect to radiation domination, the difference is about 15%.

Distribution at evaporation
The distribution of DM momentum at evaporation F (cp(t ev ), t ev ) is a superposition of all the instantaneous distributions, each redshifted appropriately from its time of emission t em (see e.g. ref. [41]) Notice that t em might be larger than t f if the initial BH temperature is smaller than the particle mass. Since we are interested in light DM, we have t em = t f .
For radiation domination from formation to evaporation, the ratio of scale factors in eq. (33) is For BH domination, the integral of eq. (33) should be split into two contributions, corresponding to a first period of radiation domination, and a second of BH domination For the second period of BH domination, starting at t BH = t f /β 2 and ending at t ev , the ratio of scale factors to be put in the integrand is For the first period of radiation domination, starting at t f and ending at t BH , the ratio of scale factors to be put in the integrand is rather For very strong BH domination the dominant contribution comes from F BH .
In order to get rid of the BH mass dependence, it is useful to define (as suggested by Lennon et al. in ref. [5]) the adimensional momentum (notice that it is not the same as in eq. (21)) with the related adimensional momentum distribution at evaporatioñ

Distribution at evaporation: geometrical optics approximation
In the geometrical optics limit of eq. (20), eq. (33) becomes We derive simplified analytical expressions in the case that f (M BH ) is constant over the BH lifetime: using eqs. (10), (16) and (17), we obtain F (cp(t ev ), t ev ) = (cp(t ev )) 2 9 2π where for full radiation domination while for full BH domination The adimensional momentum distribution at evaporation, eq. (39), becomes We show this quantity in the left panel of fig. 4 for the two scenarios of radiation domination with β =β (solid) and for full BH domination (dotted). Notice that the difference between the two scenarios is quite small. To reproduce the case of radiation domination with a different value of β, one has just to suppressF by the factor β/β. These results agree with the ones of Lennon et al. [5].

Distribution at evaporation: numerical results from BlackHawk
The quantityF derived from BlackHawk is shown in the right panel of fig. 4, assuming radiation domination with β =β (solid) and full BH domination (dotted). The suppression due to the different values of the spin is manifest. The case s = 0 is quite similar to the bosonic case of the geometrical optics limit. For higher spins, the geometrical optics approximation becomes worser. Again, we can see that the difference between radiation domination with β =β and full BH domination is quite small.

The DM phase space distribution
We now calculate the DM phase space distribution, as it is the essential ingredient to derive both the DM abundance and, using the publicly available code CLASS, the transfer function for structure formation. The DM phase space distribution (psd) per dof, f DM , at time t, is defined as where g DM is the number of DM dof, n DM the DM number density (scaling as a −3 from evaporation time), p is the DM momentum (scaling as a −1 ) and dΩ = 4π is the solid angle.
For DM produced by evaporating BHs, by using eqs. (19) and (33), we obtain that the psd at time t ev is where in the last equality we switched to the adimensional quantityF defined in eq. (39). The BH number density at the time of evaporation is related to the one at formation by Recalling the definition of β in eq. (9) we have where the last equality has been obtained by combining eq. (5) for the BH mass at the formation time and the Friedman equation (1). By exploiting the last two equations, we have Using also eqs. (24) and (28), we finally have 1 1 In order to match our expression with the results of Baldes et al. [20], we define ξ = (M P l c 2 ) 2 (kT BH ) 2 = (8π) 2 M BH M P l 2 . Using eqs. (24) and (28), we obtain that, at the evaporation time which is the same as eq. (3.8) of Baldes et al. [20], apart from the fact that they have a factor of 4 instead of 3 for BH domination, probably due to a typo. with In order to proceed further in the calculation of the DM abundance and transfer function, it is useful to establish a comparison with the well known case of a thermal DM candidate.

Psd for thermal DM
An hypothetical thermal DM decoupling a t = t d would have the following psd (per dof) where T DM (t d ) is the temperature of the DM at decoupling, which for a thermal DM candidate has to be identified with the temperature of the radiation bath from which it decouples, namely At later times, t > t d , both p(t) and T DM (t) scale as the inverse of the scale factor. It is then useful to define the time independent parameter and re-express the psd in terms of q It is also useful to express the DM temperature now, T DM (t 0 ), in units of the photon temperature now, T γ (t 0 ) 2.7 K.
If decoupling happened just before recombination, t d = t − r (or more generally when, in addition to DM, only photons, neutrinos and electrons were relativistic), clearly T DM (t − r ) = T ν (t − r ) = T γ (t − r ). Because of the successive reheating of the photons at recombination, at the present time t 0 we have T DM (t 0 ) = T ν (t 0 ) = 4 11 1/3 T γ (t 0 ).
If decoupling happened much earlier, when also other SM particles were relativistic, allowing for a possible entropy non conservation from decoupling to recombination,ᾱ(sa 3 where g * ,S (t d ) = 106.75 for the complete SM, g * ,S (t r ) = 10.75 (including photons, neutrinos and electrons). The DM temperature instead simply scales as the inverse of the scale factor Hence, the DM particles do not share the entropy release from the successive annihilations and their temperature is suppressed at recombination by a factor Since from recombination to now both T DM and T ν scale as the inverse of the scale factor, the same relation holds today Notice that CLASS actually uses the DM temperature now (or rather the temperature of the radiation bath at decoupling, rescaled to its value now) in units of the photon temperature now. For an early decoupled particle then where for the last approximation we usedᾱ = 1, g * ,S (t d ) = 106.75 and g * ,S (t r ) = 10.75.

Psd for BHs as a function of q
DM particles from BHs were never in thermal equilibrium. Nevertheless we can imagine to deal with their distribution as they were "decoupling" at t ev , so that they would have the decoupling temperature T DM (t ev ) = T R (t + ev ), which has already been derived, see eqs. (30) and (32). The time independent quantity q, defined in eq. (53), is thus where, for radiation domination while, according to eq. (32), for BH domination α BH = (9/16) 1/4 α R .
In terms of q, the psd of eq. (50) then becomes 2 with A R,BH given by eq. (51). Since the evaporation process happens at early times, eq. (59) applies also in this case, with the substitution g * ,S (t d ) → g * ,S (t ev ).
2 For comparison with previous literature, notice that Baldes et al. [20] assume TDM (tev) = TBH , so that which, in our opinion, is not consistent. In the right panel of fig. 5 we show the psd obtained from BlackHawk, taking M BH = 1 g, for radiation domination withβ (solid) and full BH domination (dashed). For radiation domination with other values of β, the psd has to be suppressed by a factor β/β. This has to be directly compared with the same situation, calculated in the geometrical optics approximation, as shown in the left panel. The suppression in the psd associated to higher spins is manifest. In the left panel we also show the thermal distribution (54) with one dof, g i = 1. It is clear that the spectrum of particles emitted by a 1 g BH is much harder than the thermal one.
The psd for other values of the BH masses can be easily reconstructed in the following way. If for M BH = 1 g the peak is at about Log 10 q ∼ 3, in general it is at Log 10 q ∼ 3 + Log 10 ( M BH 1 g ) 1/2 . In addition, the psd gets suppressed by the factor ( 1 g M BH ) 1/2 .

Dark Matter abundance
The present abundance of a stable DM particle produced by evaporation is directly related to the number-to-entropy density of such particles, Y DM (t 0 ) = n DM (t 0 )/s(t 0 ), where ρ c = 1.88 × 10 −26 h 2 kg/m 3 , h being the dimensionless Hubble parameter. Accounting for a possible entropy non-conservation from evaporation to now, namelyᾱ(sa Hence The DM number density at evaporation can be calculated by integrating over all momenta the DM spectrum at evaporation. Using the definition of the psd given in eq. (45), we have where in the second-to-last equality we changed the integration variable to q using eq. (60) and in the last one we exploited the fact that T DM (t ev ) scales as the inverse of the scale factor.
Inserting the last result in eq. (64), the DM abundance is simply given by where T DM (t 0 ) is given by eq. (59) in units of T γ (t 0 ).

DM from primordial BHs
Assuming that the DM is fully given by a stable DM candidate from the BH evaporation, we now calculate the required value for its mass. Our analysis is, to our knowledge, the first that accounts for the differences due to the spin of the DM candidate.
For comparison, we exploit both the analytical results for the psd obtained in the geometrical optics limit and the numerical ones obtained using BlackHawk. As for the dof, we consider those of the SM, with the addition of the DM candidate ones. For the DM dof we make a minimal choice, considering just those associated to polarization. For the geometrical optics approximation, we consider a boson B with g DM = 1, and a fermion F with g DM = 2. For BlackHawk we consider massive DM particles with g DM = 1, 2, 3, 4, 5 for s = 0, 1/2, 1, 3/2, 2 respectively. For M BH 10 7 g, we obtain the following results: for full BH domination m DM c 2 =ᾱ Ω DM 0.25 0.1 g * ,S (tr) g * ,S (tev) while for radiation domination where the values ofm BH,R are collected in table 1.
The scaling with the BH mass of eqs. (67) and (68) slightly breaks at M BH 10 7 g, where one obtains smaller values form BH/R , as an effect of the decrease of g * (t ev ), see fig. 3. For instance, for M BH = 10 8 (10 9 ) g, the suppression inm BH/R is by a factor of 0.74(0.56).
As expected, the result for s = 0 from BlackHawk is consistent with the geometrical optics limit for a boson. It then perfectly matches with previous literature results obtained with the alternative strategy based on the geometrical optics approximation (see e.g. [19] and references therein), which we report in fig. 6 (taken form [19] andm R are really close in the geometrical optics limit for the bosonic case and for the spin 0 numerical result, the last being slightly lower than the previous. This difference depends on the precise shape of the peak in the psd. For other spins, the suppression of the psd compared to geometrical optics limit causes a clear increase of the massesm BH/R compared to this limit. Notice that for radiation domination with β =β, the values ofm BH/R are systematically smaller by a factor of about 0.75 with respect to BH domination. Such difference could not be appreciated within the alternative strategy used in the previous literature, as the latter method does not make any difference between the two scenarios. This discrepancy is indeed due to the difference of the ratio of the scale factors in the integrands, not to the slight difference in T R (t ev ) which was already included in previous literature. with regions of BH and radiation domination. From [19].
The strong increase in the masses for increasing values of the spin allows to hope to escape bounds from structure formation. In the following section we will study in detail the impact on structure formation of these warm DM candidates, assuming entropy conservation,ᾱ = 1.
Before doing this, and in view of the comparison to be done in the next section, it is useful to recall the relevant formulas for the DM abundance for a thermal DM candidate.

DM from early decoupled thermal relics
We consider a fermion X, with g X = 2, as thermal DM candidate. It is well known that (see e.g. [16]), if decoupling happens just before reheating, so that T X (t 0 )/T ν (t 0 ) = 1, the value of the DM mass for which Ω X = 0.25, would be m X c 2 = 11 eV; for an early decoupling, when g * ,S (t d ) = 106.75, the temperature of the DM candidate now is suppressed, (T X (t 0 )/T ν (t 0 )) 3 = g * ,S (t r )/g * ,S (t d )/ᾱ = 0.1/ᾱ, so that its mass gets enhanced to m X c 2 = 110 eV forᾱ = 1.

Constraints from structure formation
Candidates of DM particles are classified according to their velocity dispersion, which defines a free-streaming length. On scales smaller than the free-streaming length, fluctuations in the DM density are erased and gravitational clustering is suppressed.
The velocity dispersion of cold DM (CDM) particles is by definition so small that the corresponding free-streaming length is irrelevant for cosmological structure formation. That of hot DM, e.g. ordinary light neutrinos, is only one or two orders of magnitude smaller than the speed of light, and smoothes out fluctuations in the total matter density even on galaxy cluster scales, which leads to strong bounds on their mass and density. Between these two limits, there exists an intermediate range of DM candidates generically called warm DM (WDM).
The matter power spectrum P (k) is very sensitive to the presence of warm DM particles with large free-streaming lengths. Due to their free-streaming velocity, warm DM particles slow down the growth of structure and suppress P (k) on scales smaller than their free-streaming scale. The effect of the free-streaming on the matter distribution can be described by a relative "transfer function" which is the square root of the ratio of the matter power spectrum in the presence of warm DM to that in the presence of purely cold DM, for fixed cosmological parameters.
For the majority of the cosmological models in which the universe contains only warm DM (in addition to the usual baryon, radiation and cosmological constant components), the transfer function can be approximated by the analytical fitting function (see e.g. [42]) where α B (labelling the scale of the break), ν and γ are free parameters sensitive to the details of the warm DM candidate. For pure warm DM models, the combined data on the CMB and the Lyman-α forest [16] provide a lower bound on the scale where the transfer function starts to fall. This lower limit was estimated in [16] to be α B 0.11 Mpc/h at the 2σ confidence level.
A well known case is the one of thermal relics 3 as warm DM. In such a case eq. (72) simplifies into (see e.g. [43]) where α B is a function of the thermal relic mass and temperature, and the index ν is fixed (ν = 1.12).
The bound on α B derived in the pioneering analysis of [16], translated into a lower bound on the thermal relic mass, gives m X c 2 0.5 keV. A more recent analysis [44] showed that the lower limit is now m X c 2 3 keV.
Standard thermal relics (those with mass m X c 2 = 11 eV) are completely ruled out, as well as early decoupled ones (those with g * (t d ) = 106.75), with mass m X c 2 = 110 eV. Only very early decoupled thermal relics could manage to be as heavy as 3 keV: this would require a huge amount of dof at decoupling. We show the transfer function for such early and very early decoupled thermal relics in fig. 7.
Using the psd's obtained with BlackHawk (taking the dof of the SM plus those of the DM candidate), for various values of M BH and of the DM particle spin, we calculated with CLASS the associated transfer function, requiring that the DM from the evaporation of primordial BHs accounts for all of the observed DM (Ω DM = 0.25). We also require that entropy is conserved, α = 1. By fixing the BH mass and the DM spin, the mass of the warm DM candidate is univoquely determined, as shown by eqs. (67) and (68).
We show the transfer functions in fig. 7. The top panels applies to the case of full BH domination (dashed lines) and the case of radiation domination with β =β (solid lines): these scenarios give quite similar results. Notice that the transfer functions of fig. 7 actually apply to all values 4 of the BH masses in the range 10 −5 − 10 7 g, for which the corresponding values of the mass providing Ω DM = 0.25 is given by eqs. (67) and (68), according to the spin. It turns out that, even for the higher spins, there is a conflict with the constraints from structure formation (at contrary with the expectations of ref. [5]).
For BHs with masses in the range 10 7 − 10 9 g, the situation is even worse, because the parameterm R/BH is smaller and the ratio T DM (t 0 )/T γ (t 0 ) is larger than in the previous case, being proportional to 1/g * ,S (t ev ) 1/3 . The only possibility left is then radiation domination with a sufficiently small value of β. In the bottom panel of fig. 7, we consider radiation domination with increasingly smaller values of β. In particular, focusing on the s = 0 case, we consider the transfer function with β =β (solid) and compare it with the ones (dot dashed) obtained taking β/β = 1/10, 1/30, 1/100. We can see that the upper limit on β/β is about 1/100: this confirms previous estimates [4,19] based on a simplified method. The region to be excluded in fig. 6 is thus the same as derived in [19].
For the different values of the spin, s = 0, 1/2, 1, 3/2, 2, by fitting the right and bottom panels of fig. 7 with eq. (72), we can derive a general formula for α B , Since for the transfer function associated to the 3 keV thermal relic one has α B 0.03 Mpc h , we can derive an upper limit on β/β. From eq. (74), for the different values of the spin, s = 0, 1/2, 1, 3/2, 2, the maximum value of β that allows to satisfy the bounds from structure formation turns out to be β β (0.013, 0.015, 0.029, 0.15, 0.16) .
One might also be interested in mixed scenarios with the simultaneous presence of both cold and warm DM, with the warm candidate coming from the evaporation of primordial BHs accounting only partially for the full DM. For the case of thermal relics it has been shown (see e.g. [16]) that in mixed models small-scale structures are not completely erased below α B and the free-streaming effect leads to a step-like transfer function, with the standard ΛCDM matter power spectrum recovered in the limit Ω W DM → 0. However, in such models the scale of the break α B becomes larger than pure warm DM models, increasing with the inverse of the mass of the DM candidate. We have verified that the same scenario arises with candidates coming from the evaporation of primordial BHs, see e.g. fig. 8. In such a case the mass of the warm candidate drops as we reduce its abundance, according to eq. (66). 9 Discussion and conclusions g. Only radiation domination is allowed, if the values β are smaller than indicated in eq. (75), according to the warm DM spin.
A couple of possibilities to evade the previous conclusions should be mentioned, which rely on the introduction of additional particles or fields. As suggested in [4], some mechanism providing entropy non conservation and taking place after the evaporation of primordial BHs (like e.g. moduli decay) might succeed in saving the warm DM candidates with BH domination. Another possibility to evade the bound would be to have a huge increase in the dof at evaporation. This is not possible to be studied within the present version of the BlackHawk code, but since analytically it is known that Ω X ∝ 1/f (M BH ) 1/2 (see e.g. [19]), the required increase in the Page function f (M BH ) would be really huge: a factor of about 10 4 for s = 0, 1/2.
While the effect of accretion does not help [19], it is plausible that the primordial BHs spin (Kerr case) might help. Primordial BHs are usually thought to form during a radiation dominated era, and thus to have really small spin. However, if formed during a transient matter dominated era, they could have sizeable to near extremal spin [45]. The Hawking radiation yield of high spin particles is grealty enhanced by the coupling with the BH spin. Hence, with the same initial density of primordial BHs, a greater fraction of their initial mass could be emitted as high spin warm DM particles.
In this work we considered non-interacting DM from primordial BHs evaporation. Interestingly enough, allowing for self-interacting DM offers the possibility to escape the structure formation bound in the light case for BH domination [46]. Thermalization in the DM sector indeed decreases the mean DM kinetic energy and, together with number-changing processes, can have a strong impact, in particular enhancing the DM relic abundance by several orders of magnitude.