Primordial Black Hole Evaporation and Dark Matter Production: II. Interplay with the Freeze-In/Out Mechanism

We study how the evaporation of primordial black holes (PBHs) can affect the production of dark matter (DM) particles through thermal processes. We consider fermionic DM interacting with Standard Model particles via a spin-1 mediator in the context of a Freeze-Out or Freeze-In mechanism. We show that when PBHs evaporate after dominating the Universe's energy density, PBHs act as a source of DM and continuously inject entropy into the visible sector that can affect the thermal production in three qualitatively different ways. We compute the annihilation cross-sections which account for the interactions between and within the PBH produced and thermally produced DM populations, and establish a set of Boltzmann equations which we solve to obtain the correct relic abundance in those different regimes and confront the results with a set of different cosmological constraints. We provide analytic formulae to calculate the relic abundance for the Freeze-Out and Freeze-In mechanism in a PBH dominated early Universe. We identify regions of the parameter space where the PBHs dilute the relic density and thermalization occurs. Furthermore, we have made our code that numerically solves the Boltzmann equations publicly available.


I. INTRODUCTION
The observation of the cosmic microwave background (CMB) radiation, and the study of its perturbations, has shed light on the composition of our Universe at the time of recombination [1,2]. As observed by the Planck satellite, the CMB spectrum's inhomogeneities have revealed that our Universe is surprisingly flat and homogeneous. Inflation is a compelling theory that explains this astonishing homogeneity and flatness. It proposes a mechanism to generate cosmological perturbations from the quantum fluctuations of a single scalar field. Furthermore, depending on the specificity of the inflation model considered, such perturbations might have been sufficiently strong that they would have gravitationally collapsed and seeded a population of primordial black holes (PBHs). This hypothetical population of PBHs would evaporate via Hawking radiation [3,4] and if their initial mass M i 10 9 g then they would have evaporated before Big Bang Nucleosynthesis (BBN). Interestingly, such a scenario may lead to the production of dark radiation in the form of gravitons and affect the spectrum of gravitational waves, in a way which could be tested experimentally in the near future [5][6][7]. Nevertheless, this possibility remains up to now unconstrained and it is pos- * andrew.cheek@uclouvain.be † lucien.heurtier@durham.ac.uk ‡ yfperezg@northwestern.edu § jessica.turner@durham.ac.uk sible that PBHs dominated the early Universe's energy budget.
In addition to the wealth of data garnered by the CMB spectrum measurements, there is a large body of evidence that dark matter (DM) constitutes approximately 26% of the Universe's energy budget. All such evidence comes from the gravitational interactions of DM. Nonetheless, hypothesizing interactions of the DM with the Standard Model (SM) provides possible production mechanisms and insights into its nature and origin. The most testable production mechanism is thermal Freeze-Out (FO) [8], as it assumes interactions between dark matter and SM particles were, at some point, frequent enough to couple the two sectors. Alternatively, Freeze-In (FI) [9,10] makes no such assumption but requires some small coupling to the SM such that the thermal bath slowly produces the required abundance of DM.
If a population of PBHs existed in the early Universe, they would affect the production of DM. In the scenario that DM is purely gravitationally interacting, and its mass was smaller than the Hawking temperature of the PBHs, PBH evaporation would be the sole source of DM production. Such a scenario has been studied extensively [5,6,[11][12][13][14][15][16][17][18]. Recently it has been shown that purely gravitationally interacting light, fermionic DM (m DM 2 MeV) would be produced relativistically by PBHs. In this mass regime, redshifting does not sufficiently cool the DM, and it would disrupt small-scale structure formation [19]. The authors of [20] came to a similar conclusion for higher spin DM candidates, and there is considerable tension on the scenario of light DM purely produced from PBH evaporation. In those differ-ent analyses, the DM mass is often neglected as compared to the PBH temperature, the geometrical-optics limit is used in order to obtain analytical results, and the code Blackhawk is used to compute numerically greybody factors. In Ref. [21], we have derived both semi-analytically and numerically the phase-space distribution of DM particles of arbitrary spins and masses, using the full greybody factors. We have also studied in thorough detail the effect of the PBH spin on DM production, enlarging the study of Ref. [22] by solving numerically Boltzmann equations and providing helpful semi-analytical expressions for the number of particles produced during evaporation including the full greybody factors.
However, the evaporation of PBHs may not account for the production of the whole DM relic abundance. This is the case if the PBH energy fraction is particularly low or if the PBH mass is large enough. In that case, it is necessary to consider other ways to produce DM in the early Universe. In Ref. [23][24][25], it was proposed that DM particles could be produced in different manners and, in particular, through the FI and FO mechanisms. It was shown that the two sources of production could conspire non-trivially to produced the DM relic abundance that is measured today and that the mechanism of entropy dilution could play an important role in the case of the FI scenario. A future detection of particle dark matter may tell us one day what is the mass of DM, how strong its interaction with SM particles are, and therefore to which extent the DM relic abundance has been produced through a FI or FO mechanism. The detection of DM could therefore impose constraints on the density fraction, mass and spin of PBHs where the PBHs would either overpopulate the DM relic density or dilute it through an entropy injection. Understanding these implications requires a greater understanding of the possible interplay between interacting DM and PBHs. This work is a step in that direction by improving on the more approximate arguments found in the literature.
This paper aims to go one step further and study how the presence of evaporating PBHs in the early Universe can affect the dynamics of the FI and FO mechanisms in addition to producing an extra contribution to the relic abundance. We highlight that two major aspects of such an interplay remain, to our knowledge, unexplored in the literature in the context of PBH evaporation: 1. A period of PBH domination leads to a non-trivial modification of the Universe's evolution, which cannot simply be reduced to an instantaneous entropy injection and therefore affects both the FI and FO mechanism dynamics in a sizeable region of the parameter space.
2. The existence of a mediator particle X between the dark and the visible sector can lead to an enhancement of the DM interactions, in particular, when m DM m X , either because additional annihilation channels can open up at large energy, or simply because boosted DM particles can accidentally hit s-channel resonances while their momentum redshifts with the Universe's expansion.
Whereas the first point was partially addressed in a more general context of early matter domination in the case of the WIMP [26] and ultraviolet FI [27] it has never been extensively studied, to our knowledge, either in the case of PBH evaporation or in the case where the FI mechanism operates at an intermediate temperature 12 . The question of the thermalization of the PBH evaporation products was raised in Ref. [23] but only in the case of cross-sections scaling like E −2 . We go beyond such an approximation and include the full energy dependence of the DM model that we consider.
For this purpose, we focus on the simple case where a fermionic DM particle interacts with SM states by exchanging a spin-1 mediator. We derive the appropriate momentum-averaged Boltzmann equations, including PBH evaporation, that are solved numerically. For this we use the infrastructure of ULYSSES [28], a publicly available python package that has been typically used to solve Boltzmann equations associated with leptogenesis. Our code, used throughout this work, is publicly available 3 . The findings of Ref. [21] provide us with full knowledge of the DM particle energy spectrum after they are produced from PBH evaporation. Using this spectrum, we question whether such particles can constitute a sizeable fraction of DM today and explore in which regime they can be expected to thermalize, either with the SM bath or with the pre-existing relic abundance of cold DM particles. We make the simplifying assumption that the population of PBHs have a monochromatic mass spectrum, which typically is the case when PBHs are produced from inflationary fluctuations or bubble collisions. More complex distributions exist depending on the scenario considered (see e.g. Ref. [29] for a review).
Our paper is organized as follows: in Sec. II we present the generic dark matter model we study and its thermal production mechanism. Then, in Sec. III, we summarize the pertinent features of PBHs relevant for this paper and in Sec. IV we discuss the interplay between thermal and PBH produced DM as well as how the PBHs can affect the evolution of the Universe which indirectly impacts DM production. The Boltzmann equations we solve, which encapsulate this interplay, are presented in Sec. V and we discuss our approach in order to solve them in the Freeze-In and Freeze-Out cases in Sec. VIII B and Sec. VIII A respectively. In Sec. VI, we derive analytical estimates of the DM relic density produced through 1 The term intermediate referring here to the fact that the FI mechanism is neither UV nor IR dominated but takes place on a resonance, at a temperature T RH > T FI > m DM . 2 Note that the authors of Ref. [24] did consider the possibility that the FO takes place during PBH domination but could only treat that case partially since they were not solving Boltzmann's equation. Freeze-In and Freeze-Out production in a background modified by the presence of PBHs and in Sec. VII we discuss the constraints on the parameter space that we use regarding the presence of warm dark matter. Finally, in Sec. VIII we present the effects of PBH evaporation on the Freeze-Out and Freeze-In mechanism as well as discussing the impact of PBH spin on DM production.

II. THERMAL PRODUCTION OF DARK MATTER PARTICLES
This Section reviews the generic particle physics model that we will study throughout this work and the two main processes of dark matter production that we will consider, aside from the evaporation of primordial black holes.

A. The Model
We consider the dark matter particle, ψ, to be a massive Dirac fermion which is a singlet under the SM gauge group but charged under a dark Abelian symmetry U (1) X . The gauge boson associated with this new symmetry is denoted as X µ . For simplicity, we assume that the mass of the latter originates from a Stückelberg mechanism [30], such that we do not need to assume the existence of any other Beyond Standard Model particles other than the dark matter particle, ψ, and the mediator, X, throughout this paper 4 . The Lagrangian can be written where f can be any fermion in equilibrium with the SM bath throughout DM production. For the sake of generality, we assume that m DM m f and so we can safely 4 Note that in principle, the existence of a UV sector may affect the following discussion as heavy particles can be produced towards the end of PBH evaporation. However, because such particles are heavy, they will only be produced in small proportions compared to DM particles and their mediator. Therefore such a contribution to the final relic abundance is expected to be subdominant.
neglect SM masses in our analysis. As is illustrated in Fig. 1, different processes can lead to the production or thermalization of DM particles from or with the SM bath, including the annihilation of SM fermions into DM particles and the annihilation of vector mediators.

B. Thermal Production
Depending on the value of the dark and visible couplings, g D and g V respectively, the DM particles may or may not thermalize with the SM bath. In the absence of PBHs, a DM relic density is expected to be produced either through the thermal decoupling of DM particles (the FO scenario), or through their out-of-equilibrium production from the annihilation of SM particles in the plasma (FI scenario). Tracking the evolution of the relic density can therefore be achieved by solving the Boltzmann equatioṅ where n DM denotes the number density of DM particles ψ, n DM,eq = g m 2 DM T /(2π 2 )K 2 (m DM /T ) is the temperature-dependent number density of DM at thermal equilibrium, and σv th is the thermally averaged cross-section of DM annihilation. Note that we use the subscript "th" to emphasize that such a thermal average is performed over the phase-space of SM particles that follow the usual Fermi-Dirac thermal distribution. Although this might sound obvious at this stage, we will see in what follows that keeping track of which thermal distribution one is averaging over will be important when PBHs are introduced. Henceforth, we will assume that the mediator mass satisfies In the case of the Freeze-Out mechanism, this condition ensures that the DM is driven out-of-equilibrium by the s-channel annihilation involving the exchange of an offshell X boson as the t-channel annihilation of DM particles into a pair of mediators is kinematically forbidden. While in the Freeze-In case, this guarantees that DM is mainly produced on the resonance when the mediator X goes on-shell at T ≈ m X . In that case, we will also assume (and check a posteriori) that the mediator couples feebly enough to the SM such that it does not reach thermal equilibrium at any time with the SM bath and that the major production mechanism is via 2 → 2 annihilation processes. Later we will assume that those two mechanisms of DM production will co-exist with a population of primordial black holes that evaporate before Big-Bang Nucleosynthesis.

C. Non-relativistic parameterization
In this paper, we will not scan over the values of the dark and visible couplings but will parametrize instead the model in terms of the SM⇐⇒DM annihilation crosssection, σv , and branching fraction of the decay of X into DM particles, Br(X → DM). In the non-relativistic limit, the former is given by where the expression of the mediator decay width, Γ X , is given in Appendix A. After it is produced by PBH evaporation, the mediator, X, eventually decays back into DM particles, according to its decay branching fraction which, in the limit, m DM , m f m X , reads When scanning over the parameter space, these relations will allow us, for any value of the non-relativistic cross-section σv, and dark branching fraction Br(X → DM), to obtain the values of the couplings g V and g D as a function of the particle masses m X and m DM .

III. PRIMORDIAL BLACK HOLE EVAPORATION
As they evaporate, PBHs can produce a population of boosted DM and mediator particles that can contribute to the final relic abundance or thermalize with the SM bath. In Ref. [21] we studied this mechanism of production extensively and described the phase-space distribution of such particles. In this Section, we recapitulate the essential elements and notations that will be used throughout the paper.
The initial PBH mass is taken to scale with the particle horizon mass at the time of PBH formation following [31] The initial fraction of the PBH energy density ρ in PBH when they are formed is related to the radiation energy density ρ in rad via the choice of the parameter β ≡ ρ in PBH /ρ in rad . However, it is common to express this fraction using the rescaling where T in denotes the temperature of the thermal plasma when PBHs are formed. As Hawking showed in [3,4], when the black holes do not have any angular momentum, they typically emit particles with a thermal spectrum whose temperature is inversely proportional to the BH mass: While they evaporate, PBHs lose mass, and therefore, according to Eq. (8), their Hawking temperature increases. Interestingly this implies that the evaporation of PBHs can always produce particles whose masses are as heavy as the Planck mass, in proportions dictated by the evaporation dynamics. Such dynamics have been extensively studied in the literature and can be described by the following equation [32,33] where M p denotes the Planck mass. The rate of mass loss will be modified if the PBH has a non-zero angular momentum. Such cases of Kerr PBHs are discussed in detail in our companion paper [21]. In Eq. 9, d 2 N i /dpdt denotes the particle emission rate per PBH for a particle species i of mass m i , spin s i and number of degrees of freedom g i . This rate can be computed in full generality for any particle spin and mass and was used to derive the spectrum of DM and mediator particles f DM produced through evaporation in Ref. [21]. The so-called mass evaporation function ε(M BH ) = i ε i (M BH ) are defined in [33,34]. It is particularly useful in order to estimate the amount of energy that is injected by one PBH into each species i at the time of evaporation, which we write as We refer the reader to Ref. [21] for a detailed study of the emission rates that are used throughout this work. Aside from SM states, PBHs emit the mediator and DM particles during their evaporation. In Ref. [21] those modes of production have been studied in depth. In particular, the emission of a spin-1 mediator acts as a secondary source for the production of DM particles, which can be enhanced if the PBHs have a non-zero spin. In Fig. 2 we reproduced the results of Ref. [21] and show the energy fraction and PBH mass that is necessary to obtain the correct relic abundance from PBH evaporation only. For a given choice of the DM mass, the region above the contours leads to an overabundance of dark matter, whereas the region below the contours leads to an underproduction of DM. Another interesting feature of the DM production from PBH evaporation is that in the region where PBHs can dominate the energy density of the Universe before evaporating (β > β c ), the relic density of DM is independent of the energy fraction β. This is because the energy densities of both the SM and  the DM sectors after evaporation are linear in β. For heavy DM candidates (m DM 10 9 GeV) PBHs of large masses need to evaporate significantly before their temperature exceeds the DM mass, which is the case after the relic density contours cross the T BH = m DM line.

IV. INTERPLAY BETWEEN PBH EVAPORATION AND THERMAL PRODUCTION OF DM PARTICLES
Throughout the Universe's evolution, the SM sector, the DM particle, ψ, and its mediator, X, are assumed to share the universe energy density with a population of primordial black holes. Therefore, the Hubble constant can be written as a function of four elementary contributions where m p ≡ M p / √ 8π denotes the reduced Planck mass. Through Hawking evaporation, PBHs are expected to radiate energy and act as a source term for all particle species in proportions dictated by their spin and mass. Consequently, PBHs can efficiently produce many particles, including those associated with a dark sector, regardless of their interaction with or their belonging to the SM sector.
As we have seen in the previous Section, the evaporation of PBHs constitutes a natural source for the production of DM particles. In addition, DM particles can be produced in our model either from FI or FO mechanism. Although in some instances, those two contributions to the final relic abundance add up without affecting each other, there are different situations in which they interfere, leading to a non-trivial evolution of the relic abundance throughout the Universe history. In Refs. [23][24][25] some of those cases were studied using a geometricaloptics approach and avoiding solving Boltzmann's equation by focusing on cases that can be studied analytically. In particular, the evaporation of PBHs was treated as an instantaneous process. Moreover, the findings of such studies strongly depend on the choice of models that was considered: In particular, DM particles were assumed to annihilate into SM states through contact operators (no mediator exchange), and the only Freeze-In scenario considered were either IR or UV dominated. For those reasons, the questions of thermalization of evaporated particles were easily avoided, and the possible evaporation of PBH after FI or FO was reduced to a simple entropy dilution factor in the computation of the relic abundance. In this work, we go one step beyond by studying in detail the temperature dependence of the different DM interactions before and after evaporation. We also treat the evaporation as a continuous process and show that its effect on the cosmological background affects the FI and FO mechanisms non-trivially. We provide a detailed description of those possible effects both analytically and numerically after solving Boltzmann's equations.

A. Modification of the Cosmological Background
Depending on the value of the fraction β, the presence of PBHs evaporating in the early Universe may or may not modify the evolution of the cosmological background. When β < β c , PBHs never dominate the energy density of the Universe and their evaporation does not affect its evolution. However, if their energy fraction is such that β > β c , they can dominate the energy density of the Universe before they evaporate. Because their mass is approximately constant before they evaporate, PBHs behave as a matter component of the Universe with the equation of state parameter ω ≈ 0 for T T ev . In this case, the Universe goes successively from the usual era of radiation domination (regime I in Fig. 3, in which ρ tot ∼ ρ SM ∝ a −4 ) to an early matterdominated era (regime II) during which the SM bath behaves as radiation (ρ tot ∝ a −3 and ρ SM ∝ a −4 ), followed by a period of entropy injection (regime III) into the SM (ρ tot ∝ a −3 and ρ SM ∝ a −3/2 ). Finally, when evaporation ends, the Universe becomes radiation dominated again (regime IV).
These four regimes lead to qualitatively different scenarios, depending on when the thermal production of DM particles through FI or FO occurs. There are three scenarios: (i) If DM's thermal production occurs during regime I, there is an preexisting relic density of DM particles in the Universe when PBHs evaporate. In that case, the FI or FO mechanism dynamics are not affected by the evaporation, but the remaining relic density gets diluted when the evaporation injects energy into the SM bath. This situation was described in detail in Ref. [24] in the context of Higgs-portal DM and [25] in the context of ultraviolet FI.
(ii) If the thermal production of DM takes place in regimes II and III, the dynamics of the thermal production are modified, and the contribution to the relic density of the FI and FO mechanisms differs from the results derived in a radiation-dominated Universe [26,35].
(iii) If thermal processes produce DM in regime IV after PBH have already evaporated, the dynamics of the FI and FO production would, of course, be unaffected by the evaporation of PBHs. The two contributions to the DM relic density may add up together in some instances, but the evaporated DM particle may also thermalize with the SM bath when produced.
When β > β c we have seen that the relic density of DM particles produced through evaporation does not depend on the fraction β. As one can see in Fig. 2, DM masses that are sufficiently low (m DM 1GeV) or sufficiently large (m DM 10 9 GeV) do not always lead to an overdensity of DM or violate the BBN bound in the region β > β c . In those regions, PBHs may dominate the energy density and evaporate later on without overclosing the Universe, while the FI or FO mechanism can occur during any of the phases described above. If PBHs never dominate the Universe's energy density, then the four regimes reduce to one single regime where the Universe evolution is unaffected by the presence of PBHs, and results for the FI and FO production are similar to those obtained in regime IV.
Denoting by M max BH (m DM ) the value of the PBH mass such that their evaporation produces the correct relic abundance of DM, the condition for the existence of such four regimes can be expressed as Equivalently, using the results of Ref. [23], this existence bound can be expressed in terms of the temperature of the plasma after evaporation In Fig. 4 we show the values of M max BH and T min ev as a function of the temperature (plain brown line). In the grey-shaded region, PBHs can only evaporate without overclosing the Universe during the radiation-dominated era, which means that the four regimes described above cannot be present in this region of the parameter space. In the green-shaded regions, PBHs can significantly reheat the Universe while producing only a small fraction of the DM relic abundance, opening the possibility that the four regimes exist in that region 5 . For this reason, if the thermal production of DM particles takes place in the window the regimes I, II and III may take place (depending on the value of β), and the FI and FO production mechanisms can be significantly affected by the evaporation of PBHs, as described above.
In the FO mechanism, such a production typically occurs at T prod m DM . Therefore, we represented by a blue dotted line the region where T ev ≈ m DM in the right panel of Fig. 4. As one can see, this line only crosses the green-shaded region around DM masses of order m DM ∼ O(10 −2 − 0) GeV meaning that the FO mechanism may only be sensitive to the PBH evaporation in this range of DM masses. We discuss this in further detail in Sec. VIII A.
However, in the FI mechanism, most of the dark matter production is achieved around T prod ≈ m X . Therefore, the dynamics of the FI mechanism may be affected by the evaporation of PBHs only if T ev ∼ T prod ≈ m X lies within the green-shaded area.
It is important to note that those considerations restrict the region of interest in terms of entropy injection during the DM production mechanism to the low DM mass region. The second green-shaded area, laying over heavier dark matter masses, is a region where the DM relic abundance produced through FI or FO has already been achieved when PBH evaporate. In this region, as it was already described in Refs. [23,24], the entropy injection into the SM bath will lead to a dilution of the preexisting DM relic abundance.

B. Thermalization of Evaporation Products
When they evaporate, PBHs may produce a large population of relativistic DM particles. The interaction rate of such DM particles with SM particles or with DM particles produced through thermal processes strongly depends on the centre-of-mass energy, E com , of the processes involved. A large boost factor for evaporated DM particles can open new channels of its annihilation. This is the case when E com > 2m X since DM particles can annihilate into a pair of mediators through a t-channel annihilation (ψψ → XX). The s-channel annihilation ψψ → X →f f can also be enhanced as soon as E com ∼ m X . Hence, it is crucial to estimate the corresponding interaction rates to consider the possible variation of the DM number density due to such processes when PBHs evaporate.
The thermalization of the evaporated products can have various effects on the DM relic density. In the FI case, because the interaction of DM particles with SM particles is typically much smaller than in the FO case, the thermalization of DM would eventually lead to an overabundance of dark matter and thus be excluded experimentally. In the FO case, the thermalization of evaporated DM particles before the FO has the effect of washing out the contribution of PBHs to the DM relic abundance. If thermalization happens after DM particles have decoupled from the SM bath, the evaporation of PBHs may destroy the predictions of the FO mechanism or recreate initial conditions for a new FO mechanism to take place. In that case, a careful treatment of the Boltzmann equation is required to track the evolution of the DM phase-space distribution after evaporation and the evolution of the relic density. We leave such calculations for future work but nevertheless estimate in the regions of parameter space where this level of calculation will be required.

V. BOLTZMANN EQUATIONS
This Section details the Boltzmann equations that we will use to track the evolution of the DM relic abundance. In addition, we include terms that track the exchange of particles between the visible and the dark sectors and terms that account for secondary production sourced by the evaporation of primordial black holes.
As mentioned above, the evaporation of PBHs transforms mass and rotational energy into a particle number. Therefore PBH evaporation acts as a source for the phase-space density distribution of the different particle species. This source term has to be defined such that energy is conserved during the evaporation process. Denoting the number density of PBHs of mass M BH evaporating in the Universe by n BH , let us define the contribution to the phase-space density distribution time derivative of the species i as where g i is the number of degrees of freedom of the species i. With such a definition, the amount of energy created in the form of particles by PBH evaporation is which exactly compensates the mass-energy loss in the PBH sector. At the level of the phase-space density distribution for the species i, it is therefore consistent to write a Boltzmann equation incorporating such a transfer of energy: in which the collision kernel C[f i ] contains the different number-changing interaction rates involving the species i.
Integrating over the phase-space, one obtains the Boltzmann equation in terms of number densitieṡ where we have defined Together with Eq. (11) and Eq. (9) we can writė The difficulty of numerically solving such equations consists of tracking the time evolution of the complete phasespace density distributions of the different species such that the collision term can be evaluated accurately at each time step. Doing so requires considerable computational resources in order to scan over the parameter space. We will thus not attempt to solve the full Boltzmann equations in this current paper. Instead, we will solve the momentum-averaged Boltzmann equations that track the number density of dark matter as a function of time. We estimate when such an approach is valid and work in regimes where either: (i) the contributions of the evaporation to the different number densities can be secluded from the particle production from the plasma (which is typically the case in the FI scenario), such that their evolution can be traced independently, or (ii) the particles produced from evaporation quickly thermalize with the plasma (which is typically the case in the FO scenario when the evaporation takes place before the dark matter freezes out). In the final Section, we will question the validity of such approximations and probe the parameter space regions in which a better treatment of the phase-space distribution evolution has to be employed. Let us now establish the Boltzmann equations that we are to use to describe the FI and FO cases.

A. Freeze-In Case
In the Freeze-In scenario, the DM particles are only very feebly coupled to the SM bath. For that reason, it is reasonable to assume -and we will check the validity of that assumption in the next section -that once DM is produced either through thermal processes or evaporation, neither the DM nor mediator thermalizes at any time throughout the universe history. Such a regime can be easily obtained by taking the limit Br(X → SM) → 0 (see e.g. Ref. [36] for a recent review). The calculation of the final dark matter relic density corresponds to summing up the three main contributions while tracking the evolution of the SM energy density in order to take into account any effects related to the injection of entropy from PBH evaporation as described in the previous section. Because we work in a regime where the particles produced from evaporation -processes (i) and (ii) -never thermalize, neither with SM particles, nor within the dark sector, it is convenient to track the separate evolutions of the DM abundance produced from evaporation (processes (i) and (ii)) and the abundance produced from thermal processes (process (iii)). Denoting the corresponding number densities by n ev DM and n th DM , we can rewrite Eq. (20) aṡ Note that in the above equations, the distinction between the thermal average over the Boltzmann distribution of SM particles in the plasma denoted by . th and the average over the phase-space distribution of the mediator X denoted by . ev is crucial since the mediator particles produced from evaporation can have an average energy E X ev T . Denoting the phase-space density distribution of the mediator particles produced through evaporation as f ev , the averaged decay width used in Eq. (21) can be expressed as When discussing the possible thermalization of the evaporation products, such a distinction between the momentum averaged performed over the thermal or evaporated distribution will also be of great importance, as we will see in the next Section.

B. Freeze-Out Case
In the FO case, the interactions between the visible and the dark sectors are strong enough to establish thermal equilibrium between the two sectors in the early Universe. After the temperature decreases, the small number density of DM particles together with a decrease of the plasma temperature enforces the decoupling of DM particles from the thermal bath. Similarly to the Freeze-In case, the PBH evaporation acts as a secondary source of DM particles. However, in the FO case, the latter particles may or may not thermalize with the SM bath when they are produced. As it was discussed in the previous Section, we will only consider the two simple cases described below:

a. Evaporation Before FO with Thermalization
If the evaporation occurs before the FO of DM particles from the plasma, we only consider the regime in which the DM particles produced through evaporation instantaneously thermalize with the SM bath. In that case, the contribution of the PBHs to the relic density is washed out, and after evaporation, the normal FO mechanism takes place. Because thermal processes between DM and SM particles are active during evaporation, it is straightforward to realize that the mediator is also thermalized at the time of evaporation. In that case, there is no need to treat the evolution of the DM and mediator number density produced through evaporation separately, and the Boltzmann equations reduce to the usual FO Boltzmann equations with a source term for radiatioṅ n DM + 3Hn DM = σv th (n 2 DM,eq − n 2 DM ) , (23a)

b. Evaporation After FO without Thermalization
If the FO mechanism occurs before the evaporation, the DM particles produced from evaporation may constitute a significant fraction of the DM relative abundance at the time of evaporation as those particles are expected to be significantly boosted. If those DM particles interact efficiently with the DM particles produced earlier via FO, it may significantly alter the FO results and lead to a non-trivial evolution of the DM relative abundance later. In order to avoid such complication, we will evaluate the capacity of the evaporation products to interact with thermally produced particles after evaporation and ensure that we consider only the case where the DM particles produced from evaporation are not able to interact efficiently after they are produced (see details in Sec. VII). In that case, we do not need to solve the full Boltzmann equations for the DM and mediator phase-space density distributions, and the different contributions to the relic density add up, similarly to the FI case. Therefore, the Boltzmann equation to consider is the same as in Eq. (21).

VI. ANALYTICAL TREATMENT
As we have seen in Sec. IV, the evaporation of PBHs whose energy fraction verifies β > β c leads to a modification of the cosmological background evolution. In this Section, we derive analytical estimations of the DM relic density produced through the Freeze-In and Freeze-Out mechanisms in the different regimes that we have exhibited in Sec. IV. The results regarding the Freeze-Out can be found in Ref. [26] in a more general context as we have adapted them to the case of PBH evaporation. It is important to note that the following expressions stand for the relic density produced exclusively through thermal processes and not from evaporation. If PBHs would produce a sizeable contribution to the DM relic abundance when they evaporate, their contribution will have to be added to our present results. An analytic derivation of the DM relic abundance produced solely from PBH evaporation, consistent with our numerical calculations up to an O(1) factor with that given in [21].

A. Cosmological Background Evolution
Before we proceed with the calculation of the relic density, it is necessary to determine the values of the scale factor at the times when PBHs start dominating the energy density (a eq ), at the critical time when energy injection starts modifying the behavior of the SM bath (a c ) and at evaporation (a ev ). The evolution of the PBH and SM energy densities will then be obtained using appropriate scaling relations between those reference points. At first, SM radiation and PBHs can be described by the scaling relations where ρ in PBH = βρ in rad . As we will see, the scaling relation of the radiation bath in regime III differs from one of regimes I and II as the entropy injection into the SM becomes relevant at this point. At matter-radiation equality, one has a eq = a in /β .
According to Ref. [37], the energy loss of PBHs is given by which allows us to writė where the effective decay constant is approximately constant during most of the evaporation process. For that reason, PBHs nearly behave like matter (w ≈ 0) during most of the evaporation phase, and the SM energy density scales like ∝ a −3/2 during that period of time [26]. Using that approximation, PBHs are expected to evaporate when Γ PBH ≡ νH ev where ν is an O(1) parameter that has to be conveniently chosen.
During PBH domination, this occurs when where m p stands for the reduced Planck mass. From this estimation, one obtains the evaporation temperature by assuming instantaneous thermalization of the evaporation products, and we have During the phase of significant entropy injection into the SM (regime III), the energy density of the SM bath can be written The critical temperature at which the regime III starts can therefore be expressed by demanding that ρ III rad (a c ) ≡ ρ rad (a c ) giving The temperatures T eq and T c can easily be obtained by using the entropy conservation in the SM bath.

B. Freeze-In Mechanism
In order to estimate the amount of DM produced through the Freeze-In mechanism, we solve the following Boltzmann equation (33) Because most of the DM production takes place on the resonance, and because couplings of the mediator to DM and SM particles are typically small in the Freeze-In mechanism, it is reasonable to use the narrow-width approximation in order to compute the integral on the center of mass energy √ s. In the limit where m f → 0, and rewriting the time derivative in terms of the scale factor, we obtain d da Solving the above equation in the four different regimes exhibited in Fig. 3 thus allows us to express H and T as a function of the scale factor in each regime and integrating from T m X to T m X . The following scaling relations are used in what follows for the calculation of the relic density: Regime I: H 2 ∝ a −4 and T ∝ a −1 , Regime II: H 2 ∝ a −3 and T ∝ a −1 , Regime III: H 2 ∝ a −3 and T ∝ a −3/8 , Regime IV: H 2 ∝ a −4 and T ∝ a −1 .
Integrating Eq. (34), we obtain for those different regimes the following estimations: where G m,n p,q are the Meijer's G functions and ρ c = 3H 2 0 m 2 p is the critical density. In practice, the argument of the Meijer's functions is larger than unity in each of the different regimes and the following limits can be used Finally, it is important to note that in Eq. (30), the choice of the parameter ν can lead to significant variations of the evaporation temperature, which can have a substantial impact on the estimation of the relic density. Moreover, in reality, the rate Γ PBH is not constant and increases during the evaporation process. This typically leads to a slight increase in the SM temperature towards the end of the evaporation process. Therefore, the energy density of the SM bath in Regime III can be overestimated if the scaling relation T ∝ a −3/8 is assumed using the correct value of the evaporation temperature. To take this deviation into account, we first fix ν = 0.45 in order to compute the effective value of a ev and T ev and check that this choice reproduces well the numerical results for our choice of parameters. Then, in Regime III, while we keep fixed the value of a ev , we re-adjust our choice to ν = 0.9 in order to calculate the value of T ev used in the expression of Ω III h 2 . In this way, the temperature of the SM bath during the FI production is comparatively lower than it would be if the scaling relation T ∝ a −3/8 would remain true up to the real evaporation temperature.

C. Freeze-Out Mechanism
The case of Freeze-Out was treated in the context where an arbitrary sector of equation-of-state parameter w reheats the universe [26]. In our case, PBHs behave mainly like matter during most of the Universe's history, corresponding to the case w = 0. However, as mentioned above, special care needs to be adopted when discussing Regime III as the evolution of T (a) might not behave like a power law at the end of the evaporation process.
Following Ref. [26], the FO production can be described by using the non-relativistic cross-section σv defined in Eq. (4) and defining N DM = a 3 n DM Defining one obtains • Regime I and IV: x FO = ln 3 2 • Regime II: • Regime III: With those definitions of x FO , one can infer the value of the relic density in each of the regimes as where s 0 = 2π 2 45 g (T 0 )T 3 0 and s eq = 2π 2 45 g (T eq )T 3 eq .

VII. CONSTRAINTS
In this Section, we review the different constraints that we consider when solving the Boltzmann equations numerically.

A. Thermalization of Evaporation Products
As described in details in Ref. [21], the DM and mediator particles produced via PBH evaporation follows a non-trivial phase-space density distribution that is dictated by the dynamics of the evaporation. In particular, for PBHs with a large Hawking temperature, those evaporation products can be significantly boosted when emitted. For an arbitrary species, i, let us introduce the phase-space distribution of the particle of that species right after evaporation as f ev (p, t ev ) ≡ dN si dp (p) . After evaporation, unless it is affected by subsequent collision processes, this distribution is simply redshifted as follows: When they are emitted, DM or mediator particles may not interact immediately with the SM plasma. However, because their phase-space density distribution evolves as the Universe expands, the corresponding interaction rates might vary and eventually lead these particles to interact efficiently with the thermal bath or with the preexisting relic abundance of DM particles produced via FO or FI. Let us consider a particle species "1" (DM or mediator) produced from evaporation and study its scattering with particle species "2" which is thermalized with the SM bath. At a given time, we denote by f ev (p 1 ) and f th (p 2 ) their respective phase-space density distributions in momentum space. We also denote the amplitude of the scattering 1+2 → 3+4 by M 1+2→3+4 , regardless of what the final states are. Immediately after PBH evaporation, the interaction rate between the two populations of particles leading to a depletion or momentum transfer of the evaporated species can be evaluated as In order to estimate the efficiency of such a scattering process after evaporation, it is convenient to define an annihilation cross-section that is averaged over both the thermal and evaporated distributions as follows Given the definition in Eq. 47, one can estimate the scattering efficiency of an evaporated particle scattering on a thermal particle by evaluating the interaction rate Similarly, one can estimate the ability of evaporated particles to self-scatter by evaluating rates of the form In our calculation, in order to ensure that thermalization does not take place -both in the FI mechanism and in the FO mechanism when evaporation takes place after the FO -we follow the following procedure: First, we compute the contribution to the DM energy density produced through evaporation and compare it to the value of the DM energy density produced via FO at the time of PBH evaporation. If this contribution does not constitute more than 10% of the total DM energy budget at the time of evaporation, we consider that later interactions will not affect the value of the final relic abundance by more than a few percents and our derivation remains valid in this case. Suppose this contribution exceeds 10% at the time of evaporation. In that case, we demand that the interaction rate of DM and mediator particles with any other particles is smaller than unity to ensure that no thermalization of the evaporation products happens, all the way from the time of evaporation to present time.
In practice, the full numerical integration of the different amplitudes over the evaporation phase-space distribution is time and resource-consuming. Therefore, it is convenient to approximate this distribution by a Boltzmann distribution to perform at least part of the integration analytically. This treatment is described in detail in Appendix B where we provide the annihilation crosssection averaged over two Boltzmann distributions with different temperatures. However, as shown in Ref. [21], approximating the distribution of evaporated particles by a Boltzmann distribution can be erroneous, especially for particles of momentum p T BH , m DM . In Fig. 5 we illustrate such an error by comparing the value of the cross-section of DM annihilation into SM model fermions σv ev+ev computed with the evaporation distribution in the geometrical-optics limit (blue dashed line) [19,21] and computed with a Boltzmann approximated distribution (green dashed line). As one can see from this figure, for T m X , the Boltzmann approximation turns out to be a good approximation. However, at a temperature T m X the s-channel resonance, the Boltzmann approximation leads to a slight underestimation of the interaction rate. This can be understood by inspecting the shape of the evaporation distribution as compared to the Boltzmann distribution. Indeed, as it was shown in Refs. [19,21], although the Boltzmann approximation is a good approximation around the peak of the distribution, the tail of the evaporation distribution at large momenta is significantly larger than in the Boltzmann approximation. Therefore, at a temperature T m X , the resonance of the distribution stands at larger momenta p T for which the evaporation distribution is larger, explaining the excess. For the sake of evaluating whether thermalization occurs, using the Boltzmann approximation is sufficient since the peak of the thermally-averaged cross-section is unaffected. For a t-channel annihilation, the calculation is not affected by the presence of any resonance, and the Boltzmann approximation turns out to be very acceptable.

B. Warm/Non-Cold Dark Matter
The dark matter particles which are produced through evaporation carry a significant momentum at the time of evaporation. After that, although they lose part of their kinetic energy via gravitational redshift, DM particles may still be quasi-relativistic at the time of structure formation and erase structures on scales below its free-streaming length. The study of the Lyman-α forest provides one of the strongest constraints on warm or noncold DM candidates (see e.g. [38,39]). In Ref. [40,41] such a constraint was derived assuming that only a fraction of the DM relic abundance is warm. In order to exploit such a constraint in our scenario, we calculated the fraction of the DM relic abundance produced via the evaporation of PBHs and calculated their average momentum p ev at the time of evaporation using the results derived in Ref. [21]. The velocity of the DM particles today is obtained by simply redshifting the value of their average momentum at the time of evaporation as Defining the fraction of evaporated particles as we can then use the constraints derived in Ref. [40,41] to determine whether a given point of the parameter space is excluded from Lyman-α measurements 6 . In particular, we note that such constraints allow the possibility that a fraction η ev 0.02 is warm today. Thus, regions of the parameter space in which PBHs produce less than 2% of the relic density are not affected by the warm dark matter constraint derived from Lyman-α measurements. Note that a more refined analysis of the WDM constraint has been achieved in Refs. [19,20] by studying the matter distribution "transfer function" using the CLASS code. The results sketched in those studies in the context of mixed warm-cold DM scenarios agree qualitatively with our findings. We leave the adaptation of this technique to our scenario for future work.

A. Freeze-Out Regime
As discussed in Sec. IV B, thermalization between the thermally produced DM and the evaporation products of the PBHs tends to occur due to the large annihilation cross-sections typical of the FO mechanism. Nonetheless, in this short Section, we study the effect PBH evaporation can have on FO and identify regions of the parameter space where dilution of the DM relic abundance can occur  due to PBH evaporation. We also highlight a regime of thermalization where a more careful treatment involving evolving the phase-space distributions of the DM should be used.
As we have seen in Sec. IV, the evaporation of PBHs may play an important role in modifying the FO of DM particles for a DM mass m DM 1GeV. In that regime of parameter space, an annihilation cross-section of order 10 −9 GeV −2 together with couplings g V , g D ∼ O(0.1) require a mediator of mass m X ∼ O(10)GeV. Such a possibility clearly has problems avoiding other experimental constraints. Any boson of this mass interacting with SM quarks or leptons would be present in collider and fixed target experiments. A potential way around is if the mediator couples only to light (sterile) neutrinos or to some other neutral fermions that remain in equilibrium with the SM when T m DM . Furthermore, indirect detection constraints could be avoided if σv is velocity suppressed, this can be achieve if our DM is Majorana. Since this work is primarily focused on the interplay between thermal DM processes and PBHs, we do not consider these constraints further.
In Fig. 6, we show how the relic abundance changes as a function of the initial PBH mass for the point β = 10 −10 , m X = 10 GeV, Br(X → DM) = 0.5 and m DM = 1.0 GeV for log 10 (σv/[GeV −2 ]) = −8.5, −9.0, −9.5 as indicated in red, blue and green respectively. The region where post-FO DM scattering processes remain active after evaporation and our calculation is no longer valid are shown by dashed lines. From Fig. 6, we observe that for log 10 (M in BH /[g]) 8.0, the PBHs have evaporated before FO occurs (corresponding to regime IV), around T ∼ m DM , and the relic abundance is unchanged by the presence of the PBHs. For log 10 (M in BH /[g]) ∼ 8.0, evaporation occurs during FO (corresponding to regime III) and the entropy injection from the PBH evaporation dilutes the relic density by approximately 80% independent of the annihilation cross-section. For heavier initial PBH masses, the Universe undergoes a stage of matter domination corresponding to regime II, but in this regime DM particles emitted by PBHs can actively scatter with other SM and/or DM particles, which invalidates our no-thermalization assumption. Naturally, the initial mass of the PBH, which triggers thermalization will change depending on the mass of the DM freezing out: the heavier the DM, the earlier the FO occurs and therefore, the lighter the PBHs would have to be to provide the dilutionary effect and thermalization. We found that for larger annihilation cross-sections, log 10 (σv/[GeV −2 ]) −8.0, thermalization always occurs for this given point.
In Fig. 7 we consider the same benchmark point as in Fig. 6 and scan over the PBH mass, energy fraction β, and annihilation cross-section, using our analytical expressions in order to estimate the FO contribution to the relic abundance. The color map stands for the value of the cross-section leading to the correct abundance. The green-shaded area labelled Thermalization stands for the region where the evaporation takes place after FO and where the DM particles produced by PBH evaporation are expected to efficiently interact with the other DM/SM particles. As we described in Sec. VII solving the Boltzmann equation for the full DM phase-space distribution is required in this region of parameter space which goes beyond our present treatment. As one can see, there exists a narrow region of the parameter space where PBHs contribute to a small fraction of the relic density, therefore requiring a smaller annihilation crosssection in order to match with observations. At smaller values of M BH the evaporation takes place before FO and thermalization processes wash out the contribution of PBHs to the relic density. We insist on the fact that this structure of the parameter space is generic to DM masses lighter than O(GeV). For that reason, there exists practically no regime in which an early PBH domination era can significantly affect the FO dynamics while escaping the post-FO thermalization constraint. Indeed, when the PBH domination is able to deplete the contribution of the FO to the relic density, the evaporation products are always sufficiently boosted in order to interact with the SM bath after evaporation. Finally, regimes with larger DM masses in which PBHs do not dominate the energy density of the universe could in principle lead to a sizeable contribution of PBH evaporation to relic abundance. However, in that case, this contribution can only be negligible in order to avoid the warm DM constraint as was shown in Ref. [5,22].

B. Freeze-In Regime
Due to the small cross-sections typical of the FI mechanism, thermalization between the population of thermally produced and PBH produced DM is not likely to occur and we expect that dilution may play an important role in certain regions of the parameter space. However, to ensure that thermalization does not occur, we check that the thermal cross-section (computed using the narrow-width approximation) multiplied by the number density of thermal DM particles does not exceed the Hubble scale at any point in the temperature evolution.
We have identified a region of the parameter space where the FI mechanism can be affected by the late time evaporation of the PBH population. In particular, this region lies in the left green-shaded triangle of the right plot of Fig. 4. The effect of PBH evaporation on FI is shown in Fig. 8 where the relic density is plotted as a function of the initial PBH mass for the two generic points which exhibits PBH domination: β = 10 −10 , m X = 10 GeV, Br(X → DM) = 0.5 and m DM = 10 −3.0 GeV (m DM = 0.1 GeV) for the left (right) plot. The solid coloured lines indicate when all relevant tests are passed, i.e. the DM is not hot and does not thermalize at or after evaporation. While the dotted coloured lines indicate regions where the DM is hot and would disrupt structure formation [40].
As we would expect, in both cases, the relic density decreases as the cross-section is reduced from 10 −41 to 10 −44 GeV −2 . Further, for log 10 (M in BH /[g]) 7, the late time entropy injection (regime III) from large mass PBHs dilutes the relic density. This dilutionary effect is more significant for smaller cross-sections as there is less DM to dilute to begin with. Interestingly, the late-time dilution offers the possibility for certain points in the parameter space, which would otherwise overproduce DM, to be consistent with the observed DM abundance. For smaller cross-sections (shown in yellow and green of the left plot of Fig. 8) and larger masses of PBHs, the relative contribution of DM produced by the PBHs compared with DM produced from the Freeze-In is larger, and therefore the hot DM constraints exclude these scenarios. We find that for larger DM masses i.e. m DM = 10 −1.0 GeV, as shown in the right plot of Fig. 8, larger initial PBH log 10 (M in BH /[g]) 7 provide the same dilution but tend to produce DM which is too hot regardless of the annihilation cross-section. The logic follows from Eq. (4) where we observe that for a fixed annihilation cross-section, branching ratio and mediator mass, the values of the visible (g V ) and dark sector couplings (g D ) increase as a function of the DM mass. Therefore, the larger the dark matter mass, the larger the proportion of DM production by the PBHs (since Γ X→DM increases quadratically in g D , see Eq. (A2)) and therefore, the warm dark matter constraint is more easily violated for heavier DM, for this set of parameters.
As smaller DM masses are less susceptible to the hot DM constraint, for these choices of parameters, we focus on the case m DM = 10 −3.0 GeV and study it in further depth. In particular, for log 10 (σv/[GeV −2 ]) = −41.0, we plot the evolution of the radiation (blue), PBH (black), thermally (red), PBH produced DM (green) and mediator (purple) energy density for log 10 (M in BH /[g]) = 5.0 (log 10 (M in BH /[g]) = 8.0) in the top left (bottom) plot of Fig. 9. The plots on the right of Fig. 9, demonstrate that these points of the parameter space never experience thermalization. Hence the two populations of DM (thermally and PBH produced) do not interact and can be treated separately. For log 10 (M in BH /[g]) = 5.0 (top plot), we observe that there is a short period of PBH domination and evaporation before freeze at T ∼ M X . Prior to evaporation, the PBH-produced DM (shown in green) is increasing in abundance. We note that the contribution of the thermally produced DM far exceeds the PBH contribution. We can contrast this scenario with the log 10 (M in BH /[g]) = 8.0 case (bottom plot) where we observe that PBH evaporation occurs after Freeze-In. As the initial PBH mass is larger and the value of β is fixed, there is a greater PBH domination than in the previous case. The effect of the entropy dump (as seen in the change in gradient of the radiation energy density shown in blue) is significant enough to dilute the relic density. Comparing the top and bottom right plots of Fig. 9, we find that the temperature evolution of σvN T H /H is also affected. While the cross-section is unaffected, the number of thermally produced dark matter in the plasma is less in the case of late-time evaporation (comparing the pink lines of the top and bottom plot of Fig. 9), and due to the reheating of the Universe, the Hubble rate is larger. Hence the suppression of σvN T H /H in the latetime evaporation case.
For the same point, with fixed β = 10 −10 , m X = 10 GeV, m DM = 10 −3.0 GeV and Br(X → DM) = 0.5, a * = 0 we would like to explore the parameter space p which we define to be: 5 ≤ log 10 (M in BH /[g]) ≤ 9, −10 ≤ log 10 (β ) ≤ −6.5 and −43.0 ≤ log 10 (σv/[GeV −2 ]) ≤ −40.0. We would like to determine which regions of p are consistent with the observed dark relic abundance at the two sigma level [42]. To perform this task, we use ULYSSES in conjunction with Multinest [43][44][45] (more precisely, pyMultiNest [46], a wrapper around Multinest written in Python). The Multinest algorithm provides a nested sampling algorithm that calculates Bayesian posterior distributions which we will utilise in order to define regions of confidence. We place flat priors on the parameter we scan in and Multinest use the log-likelihood as an objective function: where Ωh 2 ( p) is the calculated relic density for a point in the model parameter space, Ωh 2 PDG is the best-fit value of the relic density and ∆Ωh 2 is the one-sigma range of the relic abundance [42]. Once a Multinest run is finished, we use SuperPlot [47] to visualise the posterior projected onto a two-dimensional plane as shown in Fig. 10. We note that the darker blue regions shows higher posterior probabilities and conversely the lighter blue regions show regions with lower posterior probability. The regions of the parameter space, p, consistent at the one (two) sigma level with the observed relic abundance are within the solid (dashed blue) contours.
From the left plot of Fig. 10, we find that there is an expected anti-correlation between β and M in BH : for larger initial PBH masses, smaller initial fractional energy density of PBHs are required to provide the necessary dilution. Moreover, log 10 (β ) −9 provides too large a dilutionary effect for this point in the parameter space. From the right plot of Fig. 10, we observe that large cross-sections require larger values of the initial fractional energy density to recover the observed relic density. This is because a larger cross-section produces a larger relic density, requiring a larger entropy dump to dilute it sufficiently.
In Fig. 11 we used our analytical results in order to scan over the DM annihilation cross-section for a given choice of DM and mediator masses to find the value leading to the correct relic abundance. We indicate the region where the evaporation of PBHs leads to an overdensity of DM particles in brown. The orange region is excluded because PBHs produce a significant fraction of warm DM particles, which is excluded by Lyman-α measurements. The region which transition from dark blue (smaller cross-sections) to red (larger cross-section) is where the continuous dilution during and after FI from the late time evaporation of the PBHs is effective. Finally, the green area shows where the couplings are large enough to enforce the thermalization of DM and mediator particles which is inconsistent with our Freeze-In scenario. One can see that the region in which PBH dominates the Universe's energy (β > β c ) can lead to a smooth variation of the annihilation cross-section over orders of magnitude as compared to a Freeze-In scenario that would occur in a radiation-dominated background.

C. The effect of varying the mediator mass
In Sec. VIII B, we studied the effect of varying the initial PBH mass on the FI parameter space where we allowed the DM mass and cross-section to vary. We found that the late-time evaporation of heavy PBHs could sufficiently dilute the relic density such that it is consistent with observation and that for small cross-sections log 10 (σv/[GeV −2 ]) −41.0, lighter DM candidates are favoured, m DM 10 −2 GeV, for that choice of parameters, due to the hot DM constraint.
In the previous section we fixed the mediator mass at m X = 1 GeV. In this Section, we study the effect of varying the mediator mass which is important as this determines the temperature at which the FI occurs. In Fig. 12, we study the same generic point, β = 10 −10 , m DM = 10 −3.0 GeV and Br(X → DM) = 0.5, a * = 0 and log 10 (σv/[GeV −2 ]) = −41.0, while varying the mediator mass for fixed initial PBH masses, log 10 (M in BH /g) = 5.0, 6.0, 7.0, 8.0 shown in red, cyan, green and purple respectively. The dashed lines indicate where thermalization takes place and our treatment is no longer valid. We see that for smaller PBH masses, the relic abundance is not significantly affected as the mediator mass is varied. This can be seen from the solid red and blue lines of Fig. 12. This occurs because PBH evaporation occurs before Freeze-In for this scenario. However, for all shown values of initial PBH mass, DM thermalization occurs for larger mediator masses, and when thermalization occurs, the relic density plateaus. This is because the Freeze-In mechanism has become the Freeze-out scenario where the mediator mass no longer affects the relic density. Thermalization is reached in these regions of the parameter space as g 2 V g 2 D ∼ m 4 X σv/m 2 DM and hence for a fixed crosssection and DM mass, both the visible and dark sector couplings grow as m X is increased. While we have chosen a point in the model parameter space that does not suffer from the hot dark matter constraint, we find that for sufficiently large m X 10 2.5 GeV mediator, thermalization tends to occur. The effect of choosing a smaller annihilation cross-section would mean that thermalization occurs for larger mediator masses, however as we observed from the previous section, smaller cross-sections tend to be more vulnerable to the warm DM constraint.

D. The effect of PBH spin
In the previous section, we considered the effect of Schwarzschild PBHs on the FI mechanism and identified regions where DM remains cold and does not thermalize: 10 −3 m DM (GeV) 10 −2 and 10 m X (GeV) 10 1.5 . However, Hawking radiation rates are sensitive to the angular momentum of the PBH [49], and in this Section, we briefly summarize the effect of PBH spin on the FI mechanism. A more detailed discussion of the effect of spin on Hawking evaporation can be found in our companion paper [21].
In Fig. 13, we show how the relic density varies as a function of the PBH spin, a * , for various initial PBH masses with β = 10 −10 , m DM = 10 −3.0 GeV and Br(X → DM) = 0.5 and log 10 (σv/[GeV −2 ]) = −42.0. With all parameter fixed, apart from the spin of the PBH, the relic density increases as a function of a * for log 10 (M in BH /[g]) 10 6.17 . This is because spinning PBHs   evaporate faster than non-spinning PBHs and if the evaporation temperature is close to the mass of the mediator, then the thermal plasma can be reheated when Freeze-In is most effective and this contributes to an enhancement of the relic density. For instance, the largest enhance- ment, for this given point, occurs for log 10 (M in BH /[g]) = 7.2 where the effect of spin can enhance the relic density by a factor ∼ 15. For this point in the parameter space, the evaporation temperature without spin is m X /T ev,a * =0 ∼ 5, while with spin, m X /T ev,a * =0.99 ∼ 3. For larger initial PBH masses, log 10 (M in BH /[g]) > 7.2, the evaporation temperature is much lower than the mediator mass and thus the presence of spinning PBHs still increase the relic abundance, relative to non-spinning PBHs, but the effect is not as significant.
Interestingly this implies that the dilution effect dis-cussed in Sec. VIII B, can be mitigated if the PBHs have a significant spin.

IX. DISCUSSION
The evaporation of primordial black holes constitutes a natural source of gravitational production for any dark sector of particle physics. This suggests that, although the FO or FI mechanisms can lead the SM bath to produce DM particles thermally, they might not be the only contribution to the final DM relic abundance.
Two main effects of this co-production of DM from PBH evaporation and thermal processes have been studied in the literature. First, the extra contribution of PBHs to the final relic abundance leads to readjust the amount of DM particles expected from the usual thermal production to match with observations. Second, it was known that PBHs could transiently dominate the energy density of the Universe. Their subsequent evaporation into visible entropy can lead, in that case, to a dilution of a preexisting DM relic abundance. However, this evaporation was always treated as instantaneous in the computation of the dilution in previous works. In this paper, we addressed a few caveats inherent to such an approach by studying a vanilla model including a fermionic DM, ψ, which couples to SM particles through the exchange of a vector mediator, X. In particular, we focused on the case where the mediator is heavier than DM but has a mass smaller than the reheating temperature.
In the case of the Freeze-Out mechanism, we studied three different scenarios: when PBHs evaporate before FO; we checked that the DM particles produced from evaporation thermalize with the SM and therefore do not affect the dynamics of the FO mechanism nor contribute to the final relic abundance. In the case where the evaporation takes place after the FO, there exist two possible situations: (i) PBHs may dominate the energy density of the Universe before evaporating. If the FO takes place before or during a PBH dominated era (which can typically take place for m DM 1GeV), we derived useful analytical results (based on the findings of Ref. [26]) in order to calculate the relic density produced through FO and showed that the evaporation of PBHs can affect the prediction of the standard FO scenario during PBH domination. However, we showed that it is not possible that PBHs can both affect the FO production dynamics and contribute significantly to the relic density at the same time if they evaporate after the Freeze-Out time unless they destroy the FO predictions by interacting efficiently with the thermal DM remnants. This is due to the mediator's presence, which both provides additional annihilation channels for DM and leads to resonant interactions between DM and SM particles. Indeed, when they are produced through evaporation, DM particles carry a momentum p DM ∼ T BH m X , which after redshifting long enough can hit the resonance at p DM ∼ m X . Therefore, we find that the interaction of evaporated DM particles with the preexisting relic abundance or with SM particles have to be efficient if PBHs dominate the universe energy density before they evaporate. In that case, the nontrivial phase-space distribution of the evaporation products requires solving the Boltzmann equation at the level of phase-space distributions, which we let for future work.
(ii) In the case where PBHs never dominate the energy density, the contributions of the FO mechanism and PBH evaporation may sum up without affecting each other.
In the FI case, the particle interactions are small enough to avoid the thermalization constraint easily. When the temperature reaches T ∼ m X , most of the FI production takes place. In the case where PBHs dominate the energy density of the Universe before evaporating, the effect of PBH evaporation on the FI production all depends on whether this happens before, during or after PBH evaporation. We went beyond the instantaneous evaporation approximation and derived analytical results for the FI production. We scanned over the parameter space and showed that the effect of the PBH evaporation could lead to requiring an annihilation cross-section which is larger by several orders of magnitude than in the standard FI case. Thanks to the results of Ref. [21], we also derived the precise averaged momentum of the evaporated DM particles and constrained the model using Lyman-alpha constraints on warm DM.
We studied the influence of the mediator mass, as well as the possible BH spin on the value of the relic density in the Freeze-In scenario. We showed that those additional degrees can play an important role in the choice of the cross-section which is necessary in order to obtain the observed relic abundance.
Finally, we comment on the general effect of PBH evaporation on the possible detection of DM in the future. In the FO case our findings suggest that the main effect of the evaporation of PBHs on the thermal production can take place when PBHs dominate the energy density and evaporate during and after the FO. As we have seen, the evaporation products are expected to scatter efficiently with the FO relic density and with the SM bath. Our guess is that the main effect of this evaporation would therefore be to wash-out the relic density produced by FO. Although this would have to be confirmed by an appropriate solving of the Boltzmann equation, this could allow the FO mechanism to take place with smaller values of the annihilation cross-section and therefore tend to escape detection.
In the case of the FI, we have thoroughly studied how the entropy injection during and after FI can lead to significantly larger cross-sections in order to obtain the correct relic abundance. This effect enhances by orders of magnitude the detectability of our FI scenario in certain regions of the parameter space.
(A4) In the narrow width approximation, this expression becomes

(A5)
Appendix B: Thermally-Averaged Cross-Sections with Different Temperatures Significant thermalization of DM particles produced by PBH evaporation is important in this work and in-dicates when our momentum-averaged Boltzmann equations are not reliable. To determine if thermalization has occurred, we calculate the interaction rate of the PBH produced DM particles with the DM particles produced from the thermal bath. As the temperatures of these two populations can differ, we must go beyond the standard thermal-averaging calculation [50]. In this Appendix, we present the derivation of the thermally-averaged crosssection of two particles with differing temperatures: where the temperature and phase-space distribution of ψ 1 (ψ 2 ) is given by T 1 (T 2 ) and f 1 (f 2 ) respectively. We approximate the phase-space distribution of the particles to be Maxwellian i.e. f ∝ exp{−E/T } The phase-space terms can be written as a function of particle energies, three momenta and angular separation between the three momenta (θ): Applying a change of variables: assuming ψ 1 and ψ 2 have the same mass, m, the phasespace can be written as Therefore, the numerator of Eq. (B1) can be written as The integration region (E 1 , E 2 > m and |cos θ| ≤ 1) we find that this numerator is where The integration over x + yields Substituting Eq. (B9) into Eq. (B6), we find that the numerator of Eq. (B1) is We note that in the limit T 1 = T 2 , the numerator expression above simplifies to that found in [50] as expected. The denominator of Eq. (B1) is more straightforward and is given by Combining Eq. (B11) and Eq. (B13), we find that the thermally-averaged cross-section is given by Finally, performing a simple transformation to integrate with respect to z instead of s, we obtain being z min = m(T 1 + T 2 )/(T 1 T 2 ).