# Dark matter production from two evaporating PBH distributions

Arnab Chaudhuri<EMAIL_ADDRESS>Indian Institute of Technology
Gandhinagar, Gujarat 382055, India  Baradhwaj Coleppa<EMAIL_ADDRESS>Indian Institute of Technology Gandhinagar, Gujarat 382055, India  Kousik Loho
<EMAIL_ADDRESS>Indian Institute of Technology Gandhinagar, Gujarat
382055, India

###### Abstract

Particulate Dark Matter (DM), completely isolated from the Standard Model
particle sector, can be produced in the early universe from Primordial Black
Hole (PBH) evaporation. However, Big Bang Nucleosynthesis (BBN) observations
put an upper bound on the initial mass of PBH requiring the PBH to evaporate
completely before the advent of BBN. DM particles in the mass range
$\sim(1-10^{9})$ GeV can not explain the observed relic abundance for an early
matter dominated universe due to this BBN constraint. However, this assumes
the presence of only one monochromatic PBH mass distribution in the early
universe. In this work, we explore the simple possibility of achieving the
observed relic with DM masses from the above mentioned range for an early
matter dominated era with two monochromatic evaporating PBH mass distributions
and demonstrate that the fermionic DM masses consistent with BBN change
slightly.
Keywords: Dark Matter, Primordial Black Hole, Early Matter Domination.

##  I Introduction

The Standard Model (SM) of particle physics continues to be the most
successful theory that describes the physics of elementary particles and
interactions barring the gravitational ones. However, there are certain issues
with the SM, of which a very important one is that the SM can not explain the
particulate nature of Dark Matter (DM). The presence of DM is evident in the
universe [1, 2, 3, 4] and it accounts for more than one fourth of the energy
budget of the universe [5, 6]. Some of the popular mechanisms of DM production
are the WIMP (weakly interacting massive particle) [7, 8, 9, 10] and FIMP
(feebly interacting massive particle) [11, 12] scenarios which require the DM
candidate(s) to have a portal coupling with the SM sector. However, a DM
particle, completely isolated from the SM sector, can be produced from the
Primordial Black Hole (PBH) evaporation via Hawking radiation [13, 14] and is
hence immune to various direct detection [15, 16, 17] and collider constraints
[18].

PBHs can be produced in the early radiation dominated universe from quantum
fluctuations among other production mechanisms mentioned in section II and
later on dominate the energy density of the universe only to finally evaporate
away completely by radiating DM particles along with other SM states before
the advent of the Big Bang Nucleosynthesis (BBN). The phenomenology of DM
production from PBH evaporation in the context of a single PBH initial mass
has been widely studied in the literature [19, 20, 21, 22, 23, 24, 25, 26, 27,
28, 29, 30, 31] especially since the detection of gravitational waves (GW). A
recent study [32] has looked into DM production from PBH mass and spin
distributions where the PBHs are produced simultaneously at the moment when
the PBH corresponding to the peak mass value is produced and also the
evolution of the mass function is taken into account in [33]. In the scenario
of a non-interacting DM production from the evaporation of a single
monochromatic PBH distribution, a DM in the rather wide mass range of the
order ($1-10^{9}$) GeV (a more specific range is mentioned in section III)
cannot satisfy the relic ($\Omega h^{2}=0.120\pm 0.001$ as per Planck data
[5]) in the PBH dominated region of parameter space [34] due to BBN
constraints [35, 36, 37, 38, 39]. An obvious question to ask in this context
is whether this rather strong constraint can be lifted in the presence of
multiple PBHs with different initial masses. In this work, we investigate the
simplest possibility of DM being produced in the early universe from two
separate monochromatic PBH mass distributions to check if at least a portion
of this disallowed DM mass region can be redeemed. This would, of course, mean
that a richer spectrum of PBHs can easily serve as the sole originator of the
presently measured DM relic – this observation serves as the motivation for
the present work. We have done the analysis for a fermionic DM in this work
and expect a qualitative extension of the results to the other types (i.e.
scalar DM, vector DM etc.) of DM to be true.

The paper is organized as follows: in Sec. II, we provide a quick review of
PBH evaporation before discussing the DM phenomenology in Sec. III and we
offer our concluding remarks in Sec. IV.

##  II Primordial Black Holes: A Review

PBHs have been of interest to physicists for decades. They can be created when
the density fluctuations ($\frac{\delta\rho}{\rho}$) at the horizon level is
greater than unity via what is commonly known as the Zel’dovich$-$Novikov
mechanism [40, 41, 42]. Other possible ways of formation are from topological
defects [43] like the collapse of cosmic strings from second order phase
transitions [44, 45, 46], inflationary perturbations [47], or due to bubble
collisions which arise from first order phase transitions in the early
universe [48, 49, 50]. Interestingly, PBHs which originate from the collapse
of highly over-dense regions have no upper bounds on their mass - such cases
are generally studied by fitting into a mass spectrum [51]. This, however, is
not true in the case when PBHs are formed from topological defects, because
the mass of such PBHs is defined by the correlation length of the respective
phase transition [52, 53].

The study of PBHs is particularly important from a phenomenological point of
view. PBHs with mass $M_{BH}^{in}\leq 10^{9}\rm{g}$ have evaporated well
before the onset of the Big Bang Nucleosynthesis (BBN). In spite of this,
their impact on the present day universe is quite strong. They can highly
impact the baryon asymmetry of the universe [54], the fraction of dark matter
particles [55, 56], and would lead to the rise of density perturbations at
relatively small scales [57].

The initial mass $M_{BH}^{in}$ of a PBH which is formed in the radiation
dominated (RD) stage can be written in terms of the energy density of
radiation $\rho_{R}$ [58, 59, 60]:

$M^{in}_{BH}=M_{BH}(T_{in})=\frac{4\pi}{3}\gamma\frac{\rho_{R}(T_{in})}{H^{3}(T_{in})},$
(1)

where $T_{in}$ is the temperature of the universe when the PBH was formed and
$H(T_{in})$ is the value of the Hubble parameter at $T=T_{in}$. The numerical
factor $\gamma$ depends on the gravitational collapse and is taken to be $\sim
0.2$ for PBHs which are formed in the RD stage. The energy density of
radiation - which is a function of temperature - takes the form

$\rho_{R}(T)=\frac{\pi^{2}}{30}g_{*}(T)\,T^{4},$ (2)

where $g_{*}(T)$ is the relative SM degrees of freedom and is a function of
temperature as well.

As the PBH is formed, the temperature of the PBH is given by [13, 14]:

$T_{BH}=\frac{M_{p}^{2}}{M_{BH}},$ (3)

where $M_{p}$ is the reduced Planck mass i.e. $M_{p}=\frac{1}{\sqrt{8\pi G}}$
with $G$ being the gravitational constant. Upon evaporation by Hawking
radiation the PBH emits particles and the rate of mass loss of the PBH is
governed by the following equation [61, 62, 63]:

$\frac{dM_{BH}}{dt}=-\sum_{i}\mathcal{E}_{i}(M_{BH})\frac{M_{p}^{4}}{M_{BH}^{2}}.$
(4)

Here, the evaporation function $\mathcal{E}_{i}$ depends on the mass of the
$i^{th}$ particle as well as the mass of the PBH and Eqn. 4 is summed over all
the particle species. This early matter dominance in the form of PBH can
distort the cosmological history. This in fact can modify the dark sector of
the universe by injecting new dark matter particles by evaporation of the PBH.
At the moment of PBH formation, the fraction of the PBH energy density to the
total energy density denoted by $\beta$ is given by

$\beta=\frac{\rho_{BH}^{in}}{\rho_{BH}^{in}+\rho_{R}^{in}}\approx\frac{\rho_{BH}^{in}}{\rho_{R}^{in}}.$
(5)

The evolution history of the universe is also perturbed due to the
modification of the Hubble parameter – this is also mathematically represented
by $\beta$ as shown in Eqn. 5 and an often-used rescaling of this parameter
[34] is

$\beta^{\prime}=\gamma^{\frac{1}{2}}\bigg{(}\frac{g_{\star}(T_{in})}{106.75}\bigg{)}^{-\frac{1}{4}}\beta\approx\gamma^{\frac{1}{2}}\beta.$
(6)

Dominance of either component is determined by comparing the value of $\beta$
with the critical value $\beta_{c}$, which is given by [28]:

$\beta_{c}=\gamma^{-\frac{1}{2}}\left(\frac{\mathcal{G}g_{*,H}(T_{BH})}{10640\pi}\right)^{\frac{1}{2}}\frac{M_{Pl}}{M^{in}_{BH}}.$
(7)

$\mathcal{G}$ is commonly known as the grey factor and $M_{Pl}$ is the Planck
mass. If the value of $\beta$ is greater than $\beta_{c}$ the universe goes
into the PBH (matter) domination at some stage in the early universe while
$\beta<\beta_{c}$ implies the universe is radiation dominated throughout. The
evaporation of the PBH is governed by the Boltzmann equation [60]:

$\frac{d\rho_{BH}}{dt}+3H\rho_{BH}=\frac{\rho_{BH}}{M_{BH}}\frac{dM_{BH}}{dt},$
(8)

where the rate of mass decay of the PBH $\frac{dM_{BH}}{dt}$ follows Eqn. 4.

An early matter domination might change the BBN predictions and hence, in
order not to violate the well established results like the BBN temperature,
the CMB results, and structure formation which took place at the later stage
of the expansion of the universe, the PBH must evaporate before the onset of
BBN and this necessitates the initial mass of the PBH to be bounded from
above: $M^{in}_{BH}<10^{9}\rm{g}$, as mentioned above. The lower bound on PBH
mass arises from inflationary scales: $M_{BH}^{in}>10^{-1}\rm{g}$ [64].
Moreover, constraints from GW impose an upper bound on $\beta$ [65, 66, 60]:

$\beta<\frac{10^{9}\rm{g}}{M_{BH}^{in}}10^{-4}.$ (9)

We now turn to understanding how some of the constraints might play out in the
next-to-minimal scenario of two evaporating PBH distributions presenting a
detailed phenomenological study of DM production in this case in Sec. III. We
propose a new way to relax the BBN constraint in the same scenario as well.

##  III Dark Matter Phenomenology

One of the key aspects of particulate dark matter that continues to elude us
till date is its origin. In contrast to more traditional WIMP and FIMP
scenarios, DM production from PBH evaporation does not require any interaction
or portal between the SM sector and the DM candidate whereas the WIMP and FIMP
mechanisms necessitate an interaction, however small. Thus the PBH evaporation
as a DM production mechanism can explain the measured relic abundance for a DM
particle interacting only gravitationally to the SM sector as has been assumed
in this work. The key ingredients that constitute the parameter space of such
DM production from a single non-spinning PBH distribution are
$\beta^{\prime}$, $M_{BH}^{in}$, and the DM mass ($m_{DM}$). It is often
convenient to represent the relic contours of DM in the
$M_{BH}^{in}-\beta^{\prime}$ plane with each contour representing a different
$m_{DM}$ value. While in the PBH dominated region ($\beta>\beta_{c}$), the
relic is independent of $\beta^{\prime}$, in the radiation dominated region
($\beta<\beta_{c}$) the relic contours can have two different slopes depending
upon the relative values of $m_{DM}$ and the temperature of the PBH ($T_{BH}$)
[20]. To radiate DM particles, the PBH temperature has to be higher than the
DM mass. Hence, a PBH with small initial mass can easily produce DM particles
due to its high $T_{BH}^{in}$. However, for a large $M_{BH}^{in}$, the
corresponding $T_{BH}^{in}$ may not be large enough to compete with $m_{DM}$
(if $m_{DM}$ is relatively large) and hence cannot produce DM to begin with.
But as the PBH evaporates through other SM states and loses mass, its
temperature finally achieves a value high enough to start radiating DM
particles as well.

Due to BBN constraints (discussed in Sec. II) on the PBH distribution, the DM
masses that are consistent with this scenario is limited as is evident in the
literature [34]. For a fermionic DM in a PBH dominated scenario this forbidden
DM mass region is given precisely by $(2.63-6.9\times 10^{8})$ GeV for a
monochromatic mass distribution of PBH. It must also be noted that a DM
particle of mass less than $\approx$ 1 TeV is forbidden anyway by the warm DM
constraints. A few example cases, spanned throughout the DM mass range in
discussion, have been shown in Fig. 2 where in the PBH dominated region the
relic contour for a single PBH distribution (given by the dotted green curve)
falls inside the grey shaded region disallowed by BBN constraints111Of course,
it is evident from Fig. 2 that this range is perfectly allowed in a _radiation
dominated_ universe, but our goal here is to see if the same can be true in a
_matter dominated_ universe as well. This is motivated by recent studies in
non-standard cosmology wherein one can have an early matter dominated
universe.. The area above the contour is disallowed by DM overabundance. For a
scenario with monochromatic PBHs, the GW and inflation bounds mentioned in
Sec. II are given respectively by the light yellow and light green shaded
regions in Fig. 2. There is an interesting possibility of satisfying the relic
in the PBH dominated region for a DM candidate of that forbidden mass range
mentioned above with the existence of an additional PBH distribution with a
different initial mass. The primary motivation for this is the idea that
having more black holes in the spectrum at an early matter-dominated stage can
affect the Hubble parameter, as their temporally separated evaporations
enhance the DM energy density more than once and thus we should be able to
account for the observed relic abundance of DM for lower initial PBH masses
bringing the relic contour out of the BBN constrained region. However, the
presence of two distributions comes with a few new degrees of freedom in the
parameter space which demands to be defined properly and compared to the
single monochromatic case.

Before taking forward this discussion further as a comparative study between
the single-PBH and two-PBH scenarios, the details of a two-PBH scenario need
to be formulated carefully. To start with, the timeline can be described in
the following manner: in primordial times at temperature $T_{1}$, black holes
of initial mass $M_{BH_{1}}^{in}$ are created and at a later time at
temperature $T_{2}$, another set of black holes of initial mass
$M_{BH_{2}}^{in}$ come into the picture. The obvious assumptions that go
without saying are the conditions $T_{1}>T_{2}$ and hence,
$M_{BH_{1}}^{in}<M_{BH_{2}}^{in}$ (see Eqn. 1 with $\rho_{R}$ and $H$ replaced
by their dependence on the temperature of the universe). Furthermore, now
there will be two $\beta^{\prime}$ parameters that are defined as follows:

$\beta_{1}^{\prime}=\sqrt{\gamma}\frac{\rho_{BH_{1}}(T_{1})}{\rho_{BH_{1}}(T_{1})+\rho_{R}(T_{1})},\,\,\textrm{and}$
(10)
$\beta_{2}^{\prime}=\sqrt{\gamma}\frac{\rho_{BH_{1}}(T_{2})+\rho_{BH_{2}}(T_{2})}{\rho_{BH_{1}}(T_{2})+\rho_{BH_{2}}(T_{2})+\rho_{R}(T_{2})},$
(11)

where $\rho_{BH_{i}}$ denotes the energy density of the $i$-th PBH
distribution and $\rho_{R}$ corresponds to the radiation energy density of the
universe. Fig. 1 depicts one such scenario where there are two PBH
distributions contributing to the radiation and DM relic as they evaporate via
Hawking radiation. Here, we have plotted the evolution of co-moving energy
densities of both the PBH distributions and DM along with the radiation
component that is dimensionally appropriate to compare (co-moving energy
density of radiation, $\rho_{R}a^{4}$, can not be compared with the co-moving
matter). One can notice the co-moving energy density of DM increases as the
PBHs evaporate with a sudden enhancement as the PBH-I comes close to complete
evaporation. A similar trend is evident around the evaporation of the PBH-II.
The radiation component also gets similarly elevated around those two regions.
Once both the PBH distributions are evaporated completely the co-moving DM
energy density stabilises as expected. We have used the ULYSSES code [34, 67]
to generate PBH contributions as input and then solved the Boltzmann equations
mentioned in Appendix A numerically to generate this plot.

Figure 1: Evolution of various components of the early universe for the
following benchmark values: $M^{in}_{BH_{1}}=10^{5.0}$ g,
$M^{in}_{BH_{2}}=10^{9.543}$ g, $\beta^{\prime}_{1}=10^{-5.0}$,
$\beta^{\prime}_{2}=10^{-1.0}$ and $m_{DM}=10$ GeV. The relic abundance
criteria is satisfied for these benchmark values. The evolution parameter is
chosen to be inverse of the temperature of the universe ($T$) scaled
appropriately by the DM mass.

The scenario with two monochromatic distributions of different initial mass
involves more parameters as compared to the one with a single monochromatic
PBH distribution. In order to draw a comparison between the two scenarios with
mismatching number of parameters one will have to resort to a specific way of
describing the scenario and the parameter space. Here, we have chosen to
analyze the PBH-II of the two-PBH scenario in comparison to the PBH of the
single-PBH scenario. The motivation behind this choice is to understand how
the DM production from a PBH distribution (i.e. PBH-II) is modified in
presence of another preexisting DM-producing PBH distribution (i.e. PBH-I).
For that purpose we set the $\beta^{\prime}_{2}$ of the two-PBH scenario side
by side with the $\beta^{\prime}$ of the single-PBH scenario and similarly
contrast $M_{BH_{2}}^{in}$ against $M_{BH}^{in}$. In Fig. 2, we have compared
the relic contours of both single-PBH and two-PBH scenario for the same DM
mass in the $M_{BH_{2}}^{in}(M_{BH}^{in})-\beta^{\prime}_{2}(\beta^{\prime})$
plane. We have chosen four representative DM mass values from the range of
interest: one ($4\times 10^{8}$ GeV) near the upper end, two ($3$ GeV and $6$
GeV) in the lower end and the remaining one ($10^{5}$ GeV) around the middle
of that range. We have fixed the excess degrees of freedom in the parameter
space of the two PBH scenario by choosing
$Log_{10}(\beta^{\prime}_{2})-Log_{10}(\beta^{\prime}_{1})=$ 1.0 and
$Log_{10}(M_{BH_{2}}^{in})-Log_{10}(M_{BH_{1}}^{in})=$1.0, 0.5 as benchmark
values. We find that for a fixed DM mass, having two PBH spectrums can
significantly mitigate the BBN bound compared to the single-PBH scenario. It
is evident in Fig. 2 that for the two-PBH case, the DM relic abundance
criteria can be satisfied _even in the PBH dominated region_ in some cases
which were disallowed in the single-PBH case because of the BBN bounds. For
example, DM masses around the edges of the $1-10^{9}$ GeV region can now
satisfy the relic with the addition of an extra PBH in the spectrum as has
been demonstrated with DM masses of 3 GeV, 6 GeV and $4\times 10^{8}$ GeV in
Fig. 2. However, for a similar logarithmic initial PBH mass difference and
logarithmic $\beta^{\prime}$ difference, it can been seen that certain DM
masses from the midrange are still forbidden by the BBN constraint even after
the inclusion of an extra distribution of PBHs as has been shown in Fig. 2 for
a DM mass of $10^{5}$ GeV. The precise DM mass range where DM relic can not be
satisfied in the PBH dominated region now shrinks to $(7.24-8.91\times
10^{7})$ GeV for $\Delta Log_{10}(M_{BH}^{in})=1.0$ (in grams) and
$(5.56-2.67\times 10^{7})$ GeV for $\Delta Log_{10}(M_{BH}^{in})=0.5$ (in
grams). It has to be noted that the warm DM constraints [19, 68] become very
much relevant for the lower DM mass range. The DM produced from the first
monochromatic PBH mass distribution does not contribute sufficiently enough to
the final DM relic to be constrained by the warm DM limit. However the energy
density of the first distribution can still affect the evaporation temperature
of the second one. Taking this fact into account, we have implemented the warm
DM limits on the DM produced from the second PBH mass distribution and found
that there is no significant change in the parameter space for bimodal
production compared to the production from a single monochromatic PBH mass
distribution. The whole PBH dominated region and a good portion of the
radiation dominated region thus become disallowed for the DM masses of 3 GeV
and 6 GeV regardless of whether produced from a single monochromatic
distribution or from two such distributions as indicated by the light cyan
shade in Fig. 2.

Interestingly, the effect of the initial mass differences of PBH is not the
same in different DM mass regions. In the higher DM masses, $m_{DM}$ is higher
than $T_{BH_{2}}^{in}$ (the initial PBH-II temperature) and that leads to the
opposite order of relic contours for
$Log_{10}(M_{BH_{2}}^{in})-Log_{10}(M_{BH_{1}}^{in})$ values of 1.0 and 0.5 in
the PBH dominated region compared to the low DM masses (where
$m_{DM}<T_{BH_{2}}^{in}$ ) as can be seen by comparing the $m_{DM}=4\times
10^{8}$ GeV graph of Fig. 2 with the $m_{DM}=3$ GeV graph. In the $m_{DM}=6$
GeV graph, the difference of initial mass differences between the two
monochromatic distributions plays a crucial role in determining whether or not
the relic can be satisfied in the PBH dominated region. Thus we see a pattern
emerging wherein a larger initial mass difference between the two black holes
tends to accommodate smaller DM masses better while a smaller initial mass
difference does a slightly better job of safeguarding more massive DM from the
BBN bound. There are certain DM masses (as is illustrated in the
$m_{DM}=10^{5}$ GeV plot) wherein neither of our two parameter choices help in
overcoming the BBN bound. Another crucial aspect in this context is the
addition of an extra component of PBH energy density both in the numerator and
denominator in the definition of $\beta_{2}$ in the two-PBH scenario compared
to the $\beta$ in the single-PBH case, which elevates the $\beta_{2}$ in
comparison to the $\beta$ and thus reaching the PBH dominated region for
comparatively lower initial mass of PBH. Even though we are mainly interested
in the BBN constraints (light gray shaded region in the plots), the
constraints from GW and inflation are also shown in the contour plots. It is
imperative to remember that these constraints are applicable to both the PBHs.
Since $M_{BH_{1}}^{in}$ is smaller than $M_{BH_{2}}^{in}$, PBH-I automatically
satisfies the BBN bounds if satisfied by PBH-II. The same argument can be used
for the GW bounds in the $M_{BH_{1}}^{in}$-$\beta^{\prime}_{1}$ plane which
will already be satisfied once it is respected in the
$M_{BH_{2}}^{in}$-$\beta^{\prime}_{2}$ plane222We have taken a simplified
approach in applying the GW constraints on the two monochromatic PBH
distributions separately in lieu of a detailed collective understanding of the
two distributions.. Hence, in Fig. 2 the GW and BBN bounds remain the same as
the single-PBH scenario and are marked by light yellow and light grey shades
respectively. However, PBH-II satisfying the inflation bound does not
guarantee the same for PBH-I simply because $M_{BH_{1}}^{in}<M_{BH_{2}}^{in}$
as per our prescription. Hence, the inflation bound has to be applied on PBH-
II in such a way that PBH-I automatically respects the bound as well. The
three cases and corresponding contours are described as follows:

  * •

For the single-PBH case the inflation bound demands that $M_{BH}^{in}>10^{-1}$
g. This bound is shown by the light green shade. To make it easier to
comprehend, the relic contour corresponding to the single PBH scenario is
shown by the same colour i.e. the green dotted curve.

  * •

For the two-PBH scenario with logarithmic initial mass difference of 0.5, the
inflation bound needs to be modified because $M_{BH_{2}}^{in}=10^{-1}$ g
already implies that $M_{BH_{1}}^{in}=10^{-1.5}$ g and thus $M_{BH_{1}}^{in}$
is in violation of the inflation bound. So, the inflation bound on
$M_{BH_{2}}^{in}$ is modified in such a way that PBH-I also satisfies the
bound. Hence, a region given by light blue shade is disallowed in addition to
the green shaded region and the corresponding relic contour is given by the
blue dashed curve.

  * •

Similarly for the initial logarithmic mass difference of 1.0, some more extra
region will be disallowed which is given by the light red shaded part. In this
case setting a lower bound on $M_{BH_{2}}^{in}$ at $10^{0}$ g makes sure that
$M_{BH_{1}}^{in}$ respects the inflation bound of $10^{-1}$ g. Hence, the
total region, consisting of the light green, the light blue and the light red
shaded parts, is disallowed to account for both the PBHs in this case and the
corresponding relic contour is given by the solid red curve.

It is to be noted that with the addition of an extra monochromatic PBH
distribution in the spectrum, more parameter space opens up in the radiation
dominated region compared to the single-PBH case.

Figure 2: Dark matter relic contours in the
$M_{BH_{2}}^{in}(M_{BH}^{in})$-$\beta^{\prime}_{2}(\beta^{\prime})$ parameter
plane for DM masses 3 GeV (top left), 6 GeV (top right), $10^{5}$ GeV (bottom
left) and $4\times 10^{8}$ GeV (bottom right). The relic contours have been
drawn for two choices of
$Log_{10}(M_{BH_{2}}^{in})-Log_{10}(M_{BH_{1}}^{in})=$1.0 and 0.5 with
$Log_{10}(\beta^{\prime}_{2})-Log_{10}(\beta^{\prime}_{1})=$ 1.0. The
corresponding relic contours for the single-PBH scenario are also displayed.
Exclusion bounds from BBN and GW valid for all cases are given by the light
grey and light yellow shaded regions respectively. The inflation bound for a
single-PBH scenario is given by the light green shaded region. For a two-PBH
scenario the inflation bound becomes more severe as given by the additional
light blue and light red shaded regions. The light cyan shade reflects the
portion of the relic contour disallowed by the warm DM constraints.

##  IV Conclusion

PBH evaporation has furnished a method for the production of DM that is
completely isolated from the SM sector. However, for an early PBH dominated
era, the observed DM relic abundance can not be achieved for a rather large
range of DM masses as it is forbidden by the BBN predictions. In this work, we
have demonstrated that the presence of a second PBH distribution has the
potential to mitigate this constraint - specifically, employing a few
benchmark values of DM mass from that forbidden region, we have shown that in
the presence of another preexisting PBH distribution the relic abundance can
be easily satisfied for DM masses around the edges of that forbidden region.
However, the midrange DM masses still remain disallowed as has been shown for
a specific benchmark value from midrange. This is a natural consequence of
having one more PBH distribution that produces DM particles as well as an
extra component of PBH energy density encapsulated in the definition of
$\beta^{\prime}$ as $\beta^{\prime}_{2}$ in Eqn. 11. Even though we have
displayed the results by taking just a few DM masses as benchmark values, the
results seem fairly interpolatable to a range of DM masses. Our analysis shows
that some parts of the DM mass range, forbidden by BBN constraint from
satisfying the relic in the PBH dominated region, can be redeemed with the
help of this two-PBH scenario. However, the redemption of the whole region
with just two monochromatic distributions seems difficult. A thorough study of
a scenario with many PBHs with a fixed mass distribution spanned throughout
the allowed PBH mass range produced at various time instances remains an
interesting prospect in this context.

##  Appendix A Relevant Evolution Equations

The relevant Boltzmann equations along with Eqn. 8 are the following [63]:

$\frac{d\rho_{R}}{dt}+4H\rho_{R}=-\frac{\varepsilon_{SM}(M)}{\varepsilon(M)}\frac{1}{M}\frac{dM}{dt}\rho_{BH},\,\textrm{and}$
(12)
$\frac{dT}{dt}=-\frac{T}{\Delta}\bigg{(}H+\frac{\varepsilon_{SM}(M)}{\varepsilon(M)}\frac{1}{M}\frac{dM}{dt}\frac{g_{\star}(T)}{g_{\star
s}(T)}\frac{\rho_{BH}}{4\rho_{R}}\bigg{)},$ (13)

where $\Delta=1+\frac{T}{3g_{\star s}(T)}\frac{dg_{\star s}(T)}{dT}$ and the
evolution of the number density of DM from evaporating PBH is given by

$\frac{dn_{DM}}{dt}+3Hn_{DM}=\frac{\rho_{BH}}{M_{BH}}\frac{dN_{DM}}{dt},$ (14)

where $\frac{dN_{i}}{dt}$ is the emission rate of species $i$ from the
evaporating PBH.

###### Acknowledgements.

BC acknowledges support from the Department of Science and Technology, India,
under Grant CRG/2020/004171. AC thanks IIT-Gandhinagar for support. KL thanks
Sujay Shil for helpful discussions.

## References

  * Clowe et al. [2006] D. Clowe, M. Bradac, A. H. Gonzalez, M. Markevitch, S. W. Randall, C. Jones, and D. Zaritsky, Astrophys. J. Lett. 648, L109 (2006), eprint astro-ph/0608407.
  * Sofue and Rubin [2001] Y. Sofue and V. Rubin, Ann. Rev. Astron. Astrophys. 39, 137 (2001), eprint astro-ph/0010594.
  * Metcalf et al. [2004] R. B. Metcalf, L. A. Moustakas, A. J. Bunker, and I. R. Parry, Astrophys.J. 607, 43 (2004), eprint astro-ph/0309738.
  * Bartelmann [2010] M. Bartelmann, Class. Quant. Grav. 27, 233001 (2010), eprint 1010.3829.
  * Aghanim et al. [2020] N. Aghanim et al. (Planck), Astron. Astrophys. 641, A6 (2020), [Erratum: Astron.Astrophys. 652, C4 (2021)], eprint 1807.06209.
  * Hinshaw et al. [2013] G. Hinshaw et al. (WMAP), Astrophys. J. Suppl. 208, 19 (2013), eprint 1212.5226.
  * Bertone et al. [2005] G. Bertone, D. Hooper, and J. Silk, Phys.Rept. 405, 279 (2005), eprint hep-ph/0404175.
  * Bergstrom [2009] L. Bergstrom, New J.Phys. 11, 105006 (2009), eprint 0903.4849.
  * Arcadi et al. [2018] G. Arcadi, M. Dutra, P. Ghosh, M. Lindner, Y. Mambrini, M. Pierre, S. Profumo, and F. S. Queiroz, Eur. Phys. J. C 78, 203 (2018), eprint 1703.07364.
  * Bauer and Plehn [2019] M. Bauer and T. Plehn, _Yet Another Introduction to Dark Matter: The Particle Physics Approach_ , vol. 959 of _Lecture Notes in Physics_ (Springer, 2019), eprint 1705.01987.
  * Hall et al. [2010] L. J. Hall, K. Jedamzik, J. March-Russell, and S. M. West, JHEP 03, 080 (2010), eprint 0911.1120.
  * Bernal et al. [2017] N. Bernal, M. Heikinheimo, T. Tenkanen, K. Tuominen, and V. Vaskonen, Int. J. Mod. Phys. A 32, 1730023 (2017), eprint 1706.07442.
  * Hawking [1974] S. W. Hawking, Nature 248, 30 (1974).
  * Hawking [1975] S. W. Hawking, Commun. Math. Phys. 43, 199 (1975), [Erratum: Commun.Math.Phys. 46, 206 (1976)].
  * Aprile et al. [2018] E. Aprile et al. (XENON), Phys. Rev. Lett. 121, 111302 (2018), eprint 1805.12562.
  * Akerib et al. [2017] D. S. Akerib et al. (LUX), Phys. Rev. Lett. 118, 021303 (2017), eprint 1608.07648.
  * Cui et al. [2017] X. Cui et al. (PandaX-II), Phys. Rev. Lett. 119, 181302 (2017), eprint 1708.06917.
  * Boveia and Doglioni [2018] A. Boveia and C. Doglioni, Ann. Rev. Nucl. Part. Sci. 68, 429 (2018), eprint 1810.12238.
  * Baldes et al. [2020] I. Baldes, Q. Decant, D. C. Hooper, and L. Lopez-Honorez, JCAP 08, 045 (2020), eprint 2004.14773.
  * Bernal and Zapata [2021a] N. Bernal and O. Zapata, Phys. Lett. B 815, 136129 (2021a), eprint 2011.02510.
  * Dai et al. [2009] D.-C. Dai, K. Freese, and D. Stojkovic, JCAP 06, 023 (2009), eprint 0904.3331.
  * Fujita et al. [2014] T. Fujita, M. Kawasaki, K. Harigaya, and R. Matsuda, Phys. Rev. D 89, 103501 (2014), eprint 1401.1909.
  * Morrison et al. [2019] L. Morrison, S. Profumo, and Y. Yu, JCAP 05, 005 (2019), eprint 1812.10606.
  * Hooper et al. [2019] D. Hooper, G. Krnjaic, and S. D. McDermott, JHEP 08, 001 (2019), eprint 1905.01301.
  * Bernal and Zapata [2021b] N. Bernal and O. Zapata, JCAP 03, 007 (2021b), eprint 2010.09725.
  * Jyoti Das et al. [2021] S. Jyoti Das, D. Mahanta, and D. Borah, JCAP 11, 019 (2021), eprint 2104.14496.
  * Gondolo et al. [2020] P. Gondolo, P. Sandick, and B. Shams Es Haghi, Phys. Rev. D 102, 095018 (2020), eprint 2009.02424.
  * Masina [2020] I. Masina, Eur. Phys. J. Plus 135, 552 (2020), eprint 2004.04740.
  * Borah et al. [2022] D. Borah, S. Jyoti Das, and R. Roshan (2022), eprint 2208.04965.
  * Lennon et al. [2018] O. Lennon, J. March-Russell, R. Petrossian-Byrne, and H. Tillim, JCAP 04, 009 (2018), eprint 1712.07664.
  * Coleppa et al. [2023] B. Coleppa, K. Loho, and S. Shil, JCAP 06, 027 (2023), eprint 2209.06793.
  * Cheek et al. [2023] A. Cheek, L. Heurtier, Y. F. Perez-Gonzalez, and J. Turner, Phys. Rev. D 108, 015005 (2023), eprint 2212.03878.
  * Gehrman et al. [2023] T. C. Gehrman, B. Shams Es Haghi, K. Sinha, and T. Xu, JCAP 02, 062 (2023), eprint 2211.08431.
  * Cheek et al. [2022a] A. Cheek, L. Heurtier, Y. F. Perez-Gonzalez, and J. Turner, Phys. Rev. D 105, 015023 (2022a), eprint 2107.00016.
  * Sarkar [1996] S. Sarkar, Rept. Prog. Phys. 59, 1493 (1996), eprint hep-ph/9602260.
  * Kawasaki et al. [2000] M. Kawasaki, K. Kohri, and N. Sugiyama, Phys. Rev. D 62, 023506 (2000), eprint astro-ph/0002127.
  * Hannestad [2004] S. Hannestad, Phys. Rev. D 70, 043506 (2004), eprint astro-ph/0403291.
  * de Salas and Pastor [2016] P. F. de Salas and S. Pastor, JCAP 07, 051 (2016), eprint 1606.06986.
  * De Bernardis et al. [2008] F. De Bernardis, L. Pagano, and A. Melchiorri, Astroparticle Physics 30, 192 (2008), ISSN 0927-6505, URL https://www.sciencedirect.com/science/article/pii/S0927650508001230.
  * Zel’dovich and Novikov [1967] Y. B. Zel’dovich and I. D. Novikov, Soviet Astron. AJ (Engl. Transl. ), 10, 602 (1967).
  * Harrison [1970] E. R. Harrison, Phys. Rev. D 1, 2726 (1970).
  * Zeldovich [1972] Y. B. Zeldovich, Mon. Not. Roy. Astron. Soc. 160, 1P (1972).
  * Cotner and Kusenko [2017] E. Cotner and A. Kusenko, Phys. Rev. Lett. 119, 031103 (2017), eprint 1612.02529.
  * Polnarev and Zembowicz [1991] A. Polnarev and R. Zembowicz, Phys. Rev. D 43, 1106 (1991), URL https://link.aps.org/doi/10.1103/PhysRevD.43.1106.
  * Kawana and Xie [2022] K. Kawana and K.-P. Xie, Phys. Lett. B 824, 136791 (2022), eprint 2106.00111.
  * Vilenkin et al. [2018] A. Vilenkin, Y. Levin, and A. Gruzinov, JCAP 11, 008 (2018), eprint 1808.00670.
  * Choudhury and Mazumdar [2014] S. Choudhury and A. Mazumdar, Phys. Lett. B 733, 270 (2014), eprint 1307.5119.
  * Baker et al. [2021] M. J. Baker, M. Breitbach, J. Kopp, and L. Mittnacht (2021), eprint 2105.07481.
  * Hasegawa and Kawasaki [2019] F. Hasegawa and M. Kawasaki, JCAP 01, 027 (2019), eprint 1807.00463.
  * Deng and Vilenkin [2017] H. Deng and A. Vilenkin, JCAP 12, 044 (2017), eprint 1710.02865.
  * Khlopov [2013] M. Khlopov, Int. J. Mod. Phys. A 28, 1330042 (2013), eprint 1311.2468.
  * Rubin et al. [2000] S. G. Rubin, M. Y. Khlopov, and A. S. Sakharov, Grav. Cosmol. 6, 51 (2000), eprint hep-ph/0005271.
  * Martin et al. [2020] J. Martin, T. Papanikolaou, and V. Vennin, JCAP 01, 024 (2020), eprint 1907.04236.
  * Dolgov et al. [2000] A. D. Dolgov, P. D. Naselsky, and I. D. Novikov (2000), eprint astro-ph/0009407.
  * Chaudhuri and Dolgov [2021] A. Chaudhuri and A. Dolgov, J. Exp. Theor. Phys. 133, 552 (2021), eprint 2001.11219.
  * Chattopadhyay et al. [2022] P. Chattopadhyay, A. Chaudhuri, and M. Y. Khlopov (2022), eprint 2209.11288.
  * Dolgov et al. [2020] A. D. Dolgov, A. G. Kuranov, N. A. Mitichkin, S. Porey, K. A. Postnov, O. S. Sazhina, and I. V. Simkin, JCAP 12, 017 (2020), eprint 2005.00892.
  * Carr et al. [2010] B. J. Carr, K. Kohri, Y. Sendouda, and J. Yokoyama, Phys. Rev. D 81, 104019 (2010), eprint 0912.5297.
  * Carr et al. [2021] B. Carr, K. Kohri, Y. Sendouda, and J. Yokoyama, Rept. Prog. Phys. 84, 116902 (2021), eprint 2002.12778.
  * Bernal and Zapata [2021c] N. Bernal and O. Zapata, JCAP 03, 015 (2021c), eprint 2011.12306.
  * MacGibbon and Webber [1990] J. H. MacGibbon and B. R. Webber, Phys. Rev. D 41, 3052 (1990).
  * MacGibbon [1991] J. H. MacGibbon, Phys. Rev. D 44, 376 (1991).
  * Perez-Gonzalez and Turner [2021] Y. F. Perez-Gonzalez and J. Turner, Phys. Rev. D 104, 103021 (2021), eprint 2010.03565.
  * Akrami et al. [2020] Y. Akrami et al. (Planck), Astron. Astrophys. 641, A10 (2020), eprint 1807.06211.
  * Domènech et al. [2021] G. Domènech, C. Lin, and M. Sasaki, JCAP 04, 062 (2021), [Erratum: JCAP 11, E01 (2021)], eprint 2012.08151.
  * Papanikolaou et al. [2021] T. Papanikolaou, V. Vennin, and D. Langlois, JCAP 03, 053 (2021), eprint 2010.11573.
  * Cheek et al. [2022b] A. Cheek, L. Heurtier, Y. F. Perez-Gonzalez, and J. Turner, Phys. Rev. D 105, 015022 (2022b), eprint 2107.00013.
  * Baur et al. [2017] J. Baur, N. Palanque-Delabrouille, C. Yeche, A. Boyarsky, O. Ruchayskiy, E. Armengaud, and J. Lesgourgues, JCAP 12, 013 (2017), eprint 1706.03118.