Primordial Black Hole Dark Matter in the Context of Extra Dimensions

Theories of large extra dimensions (LEDs) such as the Arkani-Hamed, Dimopoulos&Dvali scenario predict a"true"Planck scale $M_\star$ near the TeV scale, while the observed $M_{pl}$ is due to the geometric effect of compact extra dimensions. These theories allow for the creation of primordial black holes (PBHs) in the early Universe, from the collisional formation and subsequent accretion of black holes in the high-temperature plasma, leading to a novel cold dark matter (sub)component. Because of their existence in a higher-dimensional space, the usual relationship between mass, radius and temperature is modified, leading to distinct behaviour with respect to their 4-dimensional counterparts. Here, we derive the cosmological creation and evolution of such PBH candidates, including the greybody factors describing their evaporation, and obtain limits on LED PBHs from direct observation of evaporation products, effects on big bang nucleosynthesis, and the cosmic microwave background angular power spectrum. Our limits cover scenarios of 2 to 6 extra dimensions, and PBH masses ranging from 10 to $10^{21}$ g. We find that for two extra dimensions, LED PBHs represent a viable dark matter candidate with a range of possible black hole masses between $10^{17}$ and $10^{23}$ g depending on the Planck scale and reheating temperature. For $M_\star = 10$ TeV, this corresponds to PBH dark matter with a mass of $M \simeq 10^{21}$ g, unconstrained by current observations. We further refine and update constraints on"ordinary"four-dimension black holes.

Randall-Sundrum model [18,19] can also result in the formation of microscopic black holes; their phenomenological implications have been discussed in other works [20][21][22][23][24]. Like the PBHs produced in the ADD model, those produced in a 5D Randall-Sundrum Type II model can accrete at early times, during the high-energy regime of the braneworld cosmology. This allows the PBHs to have longer lifetimes than 4D PBHs produced at the same era and to produce evaporation radiation that can be constrained by observations at late times. However, the amount of growth is dependent on the accretion efficiency. For concreteness and simplicity we neglect these models here.
In this work, we therefore revisit the full cosmology of primordial black holes in the presence of extra dimensions, with three important results 1) we will find a full set of constraints on LED PBHs based on recent astrophysical data, 2) we will identify the region of parameter space in which LED black holes from particle collisions in the Universe could constitute a viable dark matter candidate, and 3) we will update constraints on low-mass (≲ 10 17 g) "ordinary" four-dimensional primordial black holes.
Black holes from LEDs are constrained by two important effects: first, if they are overproduced in the early Universe, they may lead to rapid absorption and loss of the primordial plasma, leading to a matter-dominated Universe incompatible with ΛCDM. BHs that do survive into observable cosmological epochs will be constrained by their evaporation products. We will compute the so-called greybody factors that describe evaporation of these BHs, along with the spectra of secondary particles, and use these to place limits on LED PBHs from their effects on BBN, the CMB, galactic and extragalactic gamma rays. The new greybody factors and constraints are packed in the CosmoLED code, which will soon be made publicly available. In all cases, the BHs produced in LED collisions are light enough that lensing and dynamical constraints do not apply.
We will find that, in the case of n = 2 extra dimensions only, PBHs can be produced which survive until today and reproduce the observed cold dark matter abundance. These dark matter candidates require a specific combination of the Planck scale M ⋆ and reheating temperature. For M ⋆ = 10 TeV, this leads to a population of PBH dark matter with a monochromatic mass M ≃ 10 22 g, which lie in the open window between evaporation and lensing constraints.
Finally, we will provide updated constraints in the low mass range on the evaporation of ordinary 4D primordial black holes. Our inclusion of secondary particles and angular information in the 511 keV flux from positron annihilation will lead to some of the strongest constraints yet from galactic gamma rays. Our updated BBN and CMB constraints also include more precise greybody and secondary particle production than prior work, leading to similar, but modified parameter space constraints.
This article is structured as follows. In Sec. II, we describe the formation of PBHs in the LED scenario and model their accretion and evaporation, including the greybody factors appropriate to 4+n-dimensional BHs, and the hadronization and decay products from primary particles. In Sec. III, we present the observational constraints we have derived from PBH evaporation's impact on: high-energy Galactic radiation (III A), isotropic photon backgrounds (III B), the rescattering of CMB photons (III C), and the relic abundances of primordial elements from Big Bang nucleosynthesis (III D). In Sec. III E, we combine the above constraints-our full results are summarized in Fig. 17. We present our conclusions and a discussion of future prospects in Sec. IV.
Throughout the text, we use units in which c = ℏ = k B = 1 and Planck 2018 cosmological parameters of BHs remain spherically symmetric in all spatial dimensions when the horizon radius is much smaller than the size of extra dimensions, i.e. r h ≪ R. As BH mass increases, the horizon approaches the boundary of extra dimensions. Larger BHs saturate the bulk and the majority horizon area will lie in the brane. For r h ≫ R, LED BHs will behave identically to classical 4D BHs, i.e., will share the same Hawking temperature, greybody spectra and lifetime, feeling the weak 4D gravitational constant rather than the true fundamental scale M ⋆ . The exact mass above which BHs behave like ordinary 4D BHs depends on the compactification scheme. We estimate it with the mass of a 4D BH whose Schwarzschild radius matches the size of extra dimensions, As displayed in Table I, at M ⋆ = 10 TeV, the maximum LED BH mass ranges from about 10 24 g to 10 14 g as the number of extra dimensions increase from n = 2 to n = 6. As with ordinary four-dimensional BHs, LED BHs also lose mass via Hawking evaporation. However, since Hawking evaporation is geometric and the horizon area of a black hole depends on the number of extra dimensions, LED black holes will have a modified Hawking temperature [26] T H = n + 1 4πr h .
The Hawking temperature of BHs in different dimensions is depicted in the left panel of Fig. 2. LED BHs in fewer extra dimensions typically radiate particles at a lower temperature than high-n BHs. They also remain considerably colder than 4D BHs, benefiting from the low bulk Planck scale. It is also worth noting that an LED BH with mass M 4D may not share precisely the same Hawking temperature with 4D BHs of the same mass, i.e. some discontinuity might be observed during extra dimension-to-4D transition. This is expected for two reasons: 1) The radius of M 4D LED BHs is not identical to the size of extra dimensions due to the different mass-radius relations for n > 0 and n = 0.
2) The LED Hawking temperature given in Eq. (7) explicitly contains n. This discontinuity is not very large: it can be seen in Fig. 2 by observing that the solid n ̸ = 0 lines do not end exactly on the blue n = 0 (4D) line. BHs may evaporate into every degree of freedom that couples to gravity so long as it is not too thermally suppressed, i.e., the Hawking temperature is not too far below the mass of the particle. Since SM particles are confined to the brane, the emission of SM particles is limited to our three dimensional space. In contrast, gravitons are free to propagate in the bulk with significantly larger emission phase space. The distribution of particles from BH evaporation resembles a black body spectrum, up to a correction due to the gravitational potential of the BH. The emission of an SM particle degree of freedom j is given by where σ j (E) is the absorption cross section, or greybody factor, which quantifies the correction. Here, the energy of a single particle is E = p 2 + m 2 j . The greybody factor can be computed via partial wave scattering theory. It is obtained by solving the wave equation of a particle near the horizon and at infinity, and by summing up the contribution from all emission modes. Because the black hole horizon behaves as a black body, the ratio of ingoing radiation at the horizon to the ingoing radiation at infinity yields the absorption coefficient A l , which is related to the absorption cross section through for brane-localized SM particles, where the sum runs over all angular momentum modes. We follow the numerical framework outlined in [30,31] and solve for the greybody spectrum for scalars [32], fermions and gauge bosons [33] in the massless limit for non-rotating higher dimensional black holes. The effect of particle mass is mainly to introduce a lower limit for the emission spectrum [34]. The greybody factors σ s for spin s = 0, 1/2 and 1 are shown in Figure 3. We note that M ⋆ does not appear in the wave equations explicitly, and the results remain valid for an arbitrary bulk Planck scale. At E → 0, the scalar greybody factor σ 0 = 4πr 2 regardless of the number of extra dimensions. In contrast to scalars and fermions, the emission of gauge bosons is suppressed at low energies. In the high energy limit Er h ≫ 1, the greybody factors for all three particle types have the asymptotic value of σ/(πr 2 h ) → 4 −1/(n+1) (n+3) (n+3)/(n+1) /(n+1). Unlike SM particles, gravitons may propagate in the bulk and thus have access to larger phase space. The emission spectrum of gravitons is more conveniently expressed by the absorption probability |A l | 2 after integrating the angular distribution over the 3+n dimensional sphere where the multiplicities of states N l are given in Ref. [37]. Graviton emission in the bulk can be decomposed into a traceless symmetric tensor, a vector and a scalar mode. We solve for the radial parts of these three components separately and sum up their absorption probabilities. The total graviton absorption probability is displayed in the last panel of Figure 3. Our numerical results agree with the exact solutions in Refs. [37,38], but differ from Ref. [39] by a constant factor. Similarly to gauge bosons, the absorption probability is suppressed in the low energy region Er h ≪ 1. At high energies, it scales asymptotically as (Er h ) n+2 . The full BH mass loss rate is obtained by integrating the particle emission spectra in Eq. (8) and (10) while accounting for the particle degree of freedom g dof . For convenience, we define ξ and α, which are related to the BH evaporation rate by where for SM particles FIG. 3. Greybody spectra for the emission of scalars, fermions and gauge bosons in the brane, and the emission of gravitons in the bulk from the evaporation of higher dimensional black holes. Different colours correspond to n = 1 to n = 6 extra dimensions. Scaled absorption cross sections are depicted for scalar, fermions and gauge bosons, and the absorption probabilities are shown for gravitons where the contributions from scalar, vector and tensor perturbations are aggregated. The n = 0 greybody spectra for all particle types are obtained from BlackHawk [35,36]. and for gravitons It is evident that the emission probability of a particle depends on the ratio between particle mass and the Hawking temperature. When m j > T H , the emission will be exponentially suppressed. This is accounted for approximately by fitting ξ j with the functional shape where ξ j,0 is evaluated at m j = 0. Numerically, we obtain b j ≃ 0.3 and c j ≃ 1.3 for SM scalars, fermions and gauge bosons. The relevant parameters for different number of extra dimensions are given in Table II. At high temperatures T H ≫ m j , we may sum over all SM particles, gravitons and their helicity states to obtain an approximately constant value for α(n, T H ) ≃ α 0 , which is also listed in Table II. In this limit, the contribution from the total emission power of gravitons in BH mass loss ranges from 0.1% to 14.4% for four dimensional (n = 0) black holes to n = 6 dimensional black holes, as also obtained in Ref. [38]. From Eqs. (2) and (7) we find the BH mass loss rate dM/dt ∝ M −2/(n+1) . As the number of extra dimensions increases, BHs tend to evaporate faster. However, they remain substantially longer lived than 4D BHs, owing to M ⋆ ≪ M pl . The right panel of Fig. 2 shows the lifetimes of BHs and Table I lists the lightest BHs that do not entirely evaporate before today. While an n = 6 BH which does not saturate the bulk does not survive until today, n = 2 BHs as light as 10 7 g may still exist now. This has striking implications which will change the BH landscape we expect: BHs in the Universe might be lighter with a larger number density, and may thus escape gravitational lensing searches but still affect astrophysical and cosmological observations through evaporation or coalescence.

B. Black hole formation in the early Universe
The Hoop Conjecture [40,41] posits that a black hole will be formed if the impact parameter b of two colliding particles is smaller than twice the horizon radius r h . Equivalently, a microscopic black hole of mass M = E CM can be created if the center of mass energy E CM is larger than M ⋆ 1 . The BH production cross section can thus be approximated by the geometric size of the scattering The high temperature primordial plasma consisted of quarks, leptons, higgs and gauge bosons. The kinetic energy of plasma particles is characterized by the reheating temperature T RH . The plasma temperature then drops due to expansion, and could also be affected by plasma loss from accretion. Given the thermal distribution of particles, T RH need not exceed M ⋆ in order for BH production to take place. During radiation domination, the BH formation rate per unit volume per unit mass is given by [14] dΓ where g ⋆ (T ) is the effective number of relativistic particle species and we have approximated the phase space distribution as a Maxwell-Boltzmann distribution. The step function Θ is added to ensure E CM ≥ M ⋆ . If the plasma temperature T ≳ 200 GeV, then g ⋆ = 106.75. To see the asymptotic behavior, we approximate the relative velocity v rel = |⃗ v 1 − ⃗ v 2 | ≃ 1 in radiation domination and carry out the integral explicitly. This yields where K ν (x) is the modified Bessel function of the second kind. In the low temperature limit T ≪ M ⋆ ≲ M , the Bessel function K ν (M/T ) ∼ T /M exp(−M/T ). This implies that there is a limited temperature window when BHs could be copiously produced. As the plasma temperature drops below M ⋆ , BH formation becomes exponentially suppressed. Without the approximation v rel ≃ 1 Eq. (16) is evaluated to be 1 We neglect the mass loss in the formation stage and assume the minimum black hole mass M min = M⋆. For discussions see Ref. [42] and the reference therein.
The difference between Eq. (17) and (18) when integrating over M is only fractional. Considering the BH production rate is extremely susceptible to the reheating temperature, the results are rather insensitive to the choice of the production formalism. To reduce computation cost, we therefore use Eq. (17) in the numerical analysis.
C. Black hole accretion and decay in an expanding universe If BHs are produced at a plasma temperature T ≲ M ⋆ , most of them acquire a mass just above the Planck scale since more massive BH production is severely limited by kinematics. However, being immersed in the radiation bath of the primordial plasma, BHs are capable of trapping any particle that crosses the horizon and become progressively more massive. The accretion rate is proportional to the horizon area and the energy density of the plasma, with an O(1) accretion efficiency f acc depending on the mean free path of the plasma particles and the peculiar velocity of the black holes [43][44][45]: with the plasma radiation density Combining the evaporation in Eq. (11) and accretion, BH mass evolves as where β = π 120 (n + 1) 2 f acc g ⋆ and α(n, T H ) is defined implicitly in Eq. (11). Depending on the sign of the bracket on the right hand side of Eq. (21), newly born BHs may either decay away or accrete and grow. Since α varies only mildly with T H , dM/dt is susceptible to the ratio T /T H . If initially dM/dt > 0, the Hawking temperature will decrease as accretion persists, further escalating the accretion rate. The accretion halts when the Universe becomes sufficiently cold to match the Hawking temperature again, at which time BHs may have accreted enough energy to appear macroscopic. For a BH created at the mass M = M ⋆ , the watershed plasma temperature between decay and accretion reads which ranges from 0.17M ⋆ to 0.62M ⋆ for n = 1 to 6 extra dimensions assuming g ⋆ = 106.75. For concreteness we have set f acc = 1. A different accretion efficiency will slightly modify T th as T th ∝ f −1/4 acc . However, the formation of massive BHs is not shut down entirely at a temperature T < T th , as BHs that are born with a mass sufficiently higher than M ⋆ may still have low enough Hawking temperature to ensure dM/dt > 0. This amounts to the production of a BH with initial mass M i where i.e., BHs that are created at a mass above M i,min may accrete rather than decay immediately after formation. On the other hand, the production rate of M i,min > M ⋆ BHs is exponentially suppressed by M/T as seen in Eq. (17). The mass evolution can be solved in a straightforward way assuming radiation dominates throughout. Relating the plasma temperature to time using the Friedmann equations in a radiation-dominated universe, Eq. (21) becomes Some examples of the BH mass evolution are given in the left panel of Fig. 4, obtained by numerically solving Eq. (25). In cases where accretion wins out, the BH mass shoots up by many orders of magnitude at the initial stage of accretion. Because of this, the accreted matter contributes nearly all of the mass, and the final BH mass is independent of the initial BH mass M i . However, the temperature dependence of the process means that the process is very sensitive to the temperature of the plasma at production, T i . As the temperature falls T (t) ≪ T i , the black hole mass grows to its asymptotic value where γ n = f acc π 3 20 g ⋆ a 2 n n−1 n+1 . It is derived when the evaporation is neglected and g ⋆ is assumed to be constant. The asymptotic BH masses are shown as a function of T i in the right panel of Fig. 4. A dotted grey line displaying T th defined in Eq. (22) is also drawn in the middle of the panel. Right of the line, BHs of any mass above M ⋆ will accrete and grow. To the left, M i has to exceed M i,min to avoid immediate decay. The production of BHs at such temperatures is more kinematically suppressed. Special attention should be paid to n = 2 BHs. For M ⋆ = 10 TeV, if the production temperature is above 6.7 TeV, BH accretion will saturate the extra dimensions at some point. After that, they behave as four dimensional BHs and continue accreting material. As the 4D Hawking temperature drops more swiftly than that of LED BHs, Eq. (25) indicates that the accretion will become much more efficient, and an asymptotic mass is missing in this scenario. These BHs may keep accreting until the plasma density is almost exhausted, indicated by the vertical line in Fig. 4.
Next, we proceed to solve for the mass and number density of BHs produced in the primordial plasma. For more precise solutions to BH evolution that do not assume radiation domination (e.g., the BH density could be large enough to affect the expansion rate of the Universe H(t)), we must solve a set of coupled integro-differential equations detailed in Appendix A. Numerical study of these equations shows that, if the BHs are able to accrete, their mass distribution function will always be very close to a monochromatic spectrum. This can be understood qualitatively, as the evolution follows two broad scenarios.
For high reheating temperatures (T RH ≳ T th ), collisional production of BHs is efficient, and the high plasma density ensures rapid accretion. As seen in Fig. 4, BH masses quickly approach M as in a radiation dominated universe until they drain a significant fraction of energy density from the radiation bath, and the rapid cooling of the plasma suppresses the subsequent production of BHs. Here, the first BHs are created approximately with an initial number density n i ≃ ρ r (T RH )/M as (T RH ). As they grow and dominate the energy budget of the Universe, the collisional production of lighter BHs is severely limited. The accreted BHs eventually decay and dump energy into the plasma. However, they must not imprint on any cosmological observations as the dominant component of the Universe. It follows that these BHs will decay before BBN and lead to an early matter domination era.
In the second scenario, the BHs initially produced at a mass M i > M i,min accrete but their energy density remains inferior to radiation density until T ≲ eV. In a radiation dominated universe all BHs are able to accrete to a mass close to M as . The second scenario usually happens at T RH < T th , otherwise BHs will be overpopulated. Similarly to the first scenario, as the expansion of the Universe cools the plasma, BH production will also cease quickly since it is kinematically suppressed by M i,min /T . The initial BH number density is therefore given by n i = t(TRH) t(T f ) dt dM dΓ/dM . Because of the suppression, the choice of final production temperature T f does not change n i so long as dΓ/dM (T f ) ≪ dΓ/dM (T RH ). The transition between these two scenarios happens at a reheating temperature T c RH which satisfies .
Below T c RH , n i is given by the left hand side of the equation, and above that n i is determined by the right. In both scenarios, the time or temperature window for BH production is extremely limited, and BHs created during that time always accrete to similar masses, leading to a distribution that is very near to a delta function. Consequently, the integro-differential equations in Appendix A can be greatly simplified to with dM/dt given by Eq. (21) and ρ • = M n • . To solve the equations, we assume the instant production of BHs with number density n i determined by the left and right of Eq. (27), contingent on the reheating temperature. We assume the all BHs are born with the minimum mass M i = M i,min (T RH ). We then evolve the BH mass and number density as a function of time, including both accretion and evaporation. Eqs. (28) to (30) reproduce the BH mass and energy density quite precisely for low and intermediate reheating temperatures, as can be seen from Figs. 19 and 20 in the appendix. At very high reheating temperature, the production of microscopic BHs becomes more efficient than BH accretion, and BHs may not reach the asymptotic mass. The precise solution of BH spectrum and mass evolution in this scenario is quite involved, which we leave for future work. Two caveats remain for this approach. First, the connection between the first and second scenarios may not be entirely smooth as we have assumed an abrupt transition. Second, Eq. (29) assumes the entropy from BH evaporation is all dumped to the radiation plasma and thermalizes instantaneously. A dedicated study, including the effects of particle decoupling and non-thermal injection, is left for future work. Examples of the solutions are shown in Figure 5, in the presence of radiation and black holes only. For reference, we include horizontal lines that indicate the density at which BBN and matter-radiation equality occur in the standard ΛCDM scenario. For n = 2 extra dimensions, if the reheating temperature T RH = T th , BHs dominate the Universe after a mere 10 −15 s, then their number density drops as the scale factor a −3 while radiation is washed away, preventing standard Big Bang cosmology from unfolding. However, if T RH = 0.375T th , the BH energy density remains subdominant until 10 12 s, when it becomes close to the radiation density near matter-radiation equality, behaving as cold dark matter should. We have not shown evolution past this time, since these illustrative models do not include a realistic treatment of baryons, dark energy, or an additional CDM component.
For n = 6 extra dimensions and T RH = T th , we still expect BHs to exhaust the radiation density promptly. However, these BHs do not survive until matter-radiation equality, their decay at about 10 −3 s replenishes the thermal bath, causing the temperature of the plasma to decrease less efficiently. BH production and decay lead to an era of early matter domination.
Early matter domination before BBN typically does not leave any detectable features. However, the decay may produce gravitational waves which do not thermalize but still contribute to N ef f , or to the stochastic gravitation wave background to be discovered at more sensitive gravitational wave observatories. Only gravitons that are localized to the brane instead of propagating in the bulk will contribute to the stochastic gravitational waves. The greybody factor of these gravitons can be obtained by solving the wave equations of gravitons on the brane, which we leave for future work.
Next, we vary the bulk Planck scale M ⋆ and solve the evolution of BHs under different reheating temperatures T RH . We derive the constraints on T RH based on two conditions: 1) As we will find more precisely in Sec. III D 2, if BHs survive until BBN terminates, the fraction of BH energy density must be less than 10 −3 in the Universe at the neutrino decoupling temperature T dec ≃ 2.33 MeV. We conservatively require these BHs not to have evaporated significantly until 1 keV, far below the temperature when all nucleosynthesis processes freeze out. In other words, if BHs do not live long enough, they are not subject to this BBN constraint. 2) If BHs survive until the plasma temperature drops to about 0.75 eV, when matter radiation equality is expected in standard cosmology, BHs must remain subdominant in order not to change the sound horizon in a significant way, which would contradict CMB observations. The results are shown in Fig. 6. The dash-dotted line corresponds to condition 1) and the solid line stems from condition 2). The regions above the lines are excluded. For n = 2 these two conditions yield very similar constraints, while the BBN constraint tends to be stronger starting from n = 3 when T RH ≳ 30 TeV. This behaviour can be understood intuitively from Fig. 2. As n rises, BH lifetime decreases sharply. The reheating temperature has to be high enough to produce massive BHs that live until matter-radiation equality, rendering weaker constraints. The same applies to the BBN condition where the constraints on T RH are weaker for larger number of extra dimensions. For n = 2, the solid line produces the right amount of BHs as 100% dark matter, which remain until today. For n ≥ 3, no reheating temperature is found such that BHs can dominate the dark matter density today for M ⋆ ≲ 10 TeV.
To investigate LED BHs as part of the dark matter today, we therefore focus on n = 2. We also add a flexible non-BH dark matter component to Eqs. (28) to (30) and evolve the energy density of dark matter over time. The non-BH dark matter energy density is adjusted such that the total cold dark matter density matches observations, Ω c h 2 = 0.120 while fixing the dark energy and baryon density to the Planck 2018 best fit [46]. We then solve for the fraction of dark matter today that is comprised of BHs, f •,0 ≡ Ω • /Ω c . We show the reheating temperature and BH mass today in Fig. 7 that corresponds to a specific f •,0 by varying M ⋆ . Since the BH production rate is exponentially suppressed when T RH ≪ M ⋆ , a minuscule change in the reheating temperature will alter f •,0 remarkably. The required reheating temperature to produce BH dark matter is roughly proportional to M ⋆ . As indicated in Eq. (26), the asymptotic BH mass, and hence the BH mass today M ∝ M −2 ⋆ with n = 2. Indeed, the fit to Fig. 7 reveals For M ⋆ ≳ 10 TeV, which evades collider constraints, the primordial BH mass today ranges from 10 17 g to 10 21 g for Planck scales below a PeV.
FIG. 6. Constraints on the reheating temperature as a function of the fundamental Planck scale M⋆. Above the solid lines, the energy density of BHs ρ• exceeds that of radiation at a photon temperature Teq = 0.75 eV, when matter-radiation equality is expected. Above the dash-dotted line, ρ• > 10 −3 ρr at the neutrino decoupling temperature T ν,de ≃ 2.33 MeV. The regions above the lines are excluded due to BH distortions to BBN or CMB. For the Planck scales considered, n = 2 BHs along the solid line will always survive until today and make the correct relic abundance, while n = 3 BHs with M⋆ < 1.4 TeV still exist. The latter n = 3 range has been excluded by collider searches.

D. Observable evaporation products
In Sec. II A we have described the primary particles from BH evaporation. If the only important observable effects of BH evaporation are the change in BH mass and injection of energy into the plasma of the early Universe, then Eqs. (8) and (10) are sufficient. However, observable stable particles (here, photons, electrons and positrons) are also produced as decay or hadronization products from heavy and coloured primary states. To correctly account for production of these secondary observable particles, we consider the contribution from several sources. As a first step, we use tabulated spectra from PPPC4DMID [47] to compute the secondary particle spectra generated from primary particles above E P = 5 GeV, which arises from the limitation of particle energy in PYTHIA [48], used in PPPC4DMID for the production of tabulated values. Below this energy, the unstable states that we include are the τ ± leptons, muons and pions. We extrapolate the τ decay spectra from PPPC4DMID down to E = m τ . The e ± and γ spectra from π and µ decay are computed and boosted to the lab frame in a similar way to Ref. [49], taking care to include the electron mass where appropriate-see Appendix B for details. These are added to the primary electrons and photons below E P produced by the evaporating BH. Overall, the secondary spectra are computed as where i = {e ± , γ}, j = {e ± , µ ± , qq, W ± , Z, g, γ, h, ν e , ν µ , ν τ }, and k = {µ ± , π ± , π 0 }. The BH primary emission spectrum is To account for QCD confinement transition, we adopt the same prescription as in Ref. [50] and include a factor where the plus sign applies for π ± and π 0 , and the minus sign for quarks and gluons. For all other species, Q j = 1. We take the confinement scale Λ QCD ≃ 300 MeV and σ = 0.1. Below Λ QCD , the emission of quarks and gluons from BH will be exponentially suppressed and the emission of hadrons is preferred. For comparison with prior work, we show the emission spectra of e ± and γ in Fig. 8, for 4D (n = 0, M ⋆ = M pl ) black holes. Our code, CosmoLED, computes the spectra of observable products from BHs, and the cosmological constraints. The dashed lines are the primary spectra obtained from Eqs. (33) and (8). The solid lines depict the total spectra of e ± and γ by considering the decay and hadronization of more energetic particles. The CosmoLED total spectra are computed using Eq. (32). For comparison, we also show the spectra obtained directly from the ExoCLASS package [50], and BlackHawk v2.1 [35,36]. Note that the ExoCLASS BH module computes the secondaries from muon and pion decay only, and BlackHawk cascades down from 5 to 10 5 GeV primary particles with the "PYTHIA" hadronization option at the present epoch. Our results agree well with BlackHawk at almost all energies, while ExoCLASS tends to underestimate the secondary spectra. As BH mass increases, the difference between CosmoLED and ExoCLASS spectra becomes less dramatic as fewer primary particles are produced above 5 GeV. However, the CosmoLED spectra remain to be larger in most of the energy range. Throughout, we assume that BHs evaporate only to standard model particles and gravitons.

III. OBSERVATIONAL CONSTRAINTS ON LED BLACK HOLES
Once produced, primordial black holes born of microscopic collisions in the early Universe will exhibit similar phenomenology to their four-dimensional cousins. In addition to affecting the energy budget of the Universe, their evaporation products will affect cosmological evolution and can interfere with Big Bang Nucleosynthesis (BBN) and the CMB, as well as produce a detectable flux of galactic and extragalactic X-rays. These constraints will not probe 4D BH masses larger than ∼ 10 18 g, and thus do not overlap with constraints from lensing and dynamical disruption of gravitational systems. In this section, we compute the dominant constraints from X-rays (Secs. III A and III B), the CMB (Sec. III C) and BBN (Sec. III D), first describing the physics, and then producing constraints from observational data. We then discuss the combined constraints (Sec. III E) as well as previous PBH constraints not studied in this work (Sec. III F).
In order to consistently compare constraints on PBHs with differing lifetimes, we define the parameter f • as where ρ • is the density of PBHs at an initial redshift, z i , before the PBHs evaporate any significant fraction of their mass and ρ DM,0 is the observed dark matter density today. With this definition, f • describes the fraction of dark matter in the early Universe comprised of PBHs. For certain observable constraints, other parameters are used to describe the abundance of PBHs. When studying galactic centre constraints we use f •,0 , the fraction of dark matter comprised of PBHs today and when studying the impact PBHs have on the expansion history near BBN we use β dec , the fraction of the total energy density comprised of PBHs at the time of neutrino decoupling.

A. Galactic constraints
For PBHs that survive until the present, the Milky Way halo is a promising source of evaporation products. Detectable sub-GeV evaporation products can consist of gamma rays, and cosmic ray electrons, positrons, protons and antiprotons. The "prompt" gamma ray flux is given by: where the D-factor is defined as an integral over the dark matter density ρ DM (⃗ x): where the integral in x is over the line of sight (l.o.s.) and ∆Ω is the solid angle of interest. ρ DM is the DM density in the Milky Way. We take it to follow an NFW profile where r is the galactocentric distance, r s is the DM halo scale radius, and the conventional the factor of 2 3−γ ensures that ρ s ≡ ρ(r s ). We employ parameters consistent with kinematic data [51] 2 , r s = 9 kpc and γ = 1.2. The DM density at the Sun's position is ρ 0 = ρ(R ⊙ ) = 0.3 GeV cm −3 , where we use recent measurements from GRAVITY [53] for the distance to the galactic centre R 0 = 8.127 kpc.
In addition to the gamma ray flux from Eq. (36), low-energy positrons produced by BH evaporation will lead to a gamma ray line signal at E γ = 511 keV from e + e − annihilation in the interstellar medium. The flux of photons from in-situ e + e − annihilation is: where f P is the positronium formation fraction and dN e + /dt is the total positron production rate per BH integrated over energy.
We employ data from INTEGRAL/SPI, the X/gamma-ray spectrometer onboard the ESA INTEGRAL satellite, launched in 2003. A full analysis of SPI data requires a template-based likelihood analysis, as there is no way to reconstruct the direction of a single photon event. Rather, SPI uses a coded mask, for which each individual photon recorded on the detector corresponds to a number of possible trajectories. This means that an image cannot be reconstructed, and one must instead compare templates using a maximum likelihood method. To sidestep this cumbersome process, we use previously-processed data reported in Ref. [54]. Although this is based on only 6 years (∼ 10 8 s) of data, it is the only published reference to include a binned reconstruction of the diffuse flux as a function of energy and galactic latitude and longitude. We follow a similar method to Ref. [55], who used this data to constrain 4D primordial black holes in the Milky Way. We employ the 5 energy bins in Figure 5 of Ref. [54] (digitized from [56]), corresponding to 27-49 keV, 49-90 keV, 100-200 keV, 200-600 keV and 600-1800 keV. These each consist of 21 latitude bins within −90 • < b < 90 • , integrated over longitudes −23.1 < ℓ < 23.1 • , with the exception of the 800-1800 keV range, which is presented in 15 bins, within −60 • < ℓ < 60 • . We do not employ the results from Figure 4, as they are drawn from the same data, but binned over latitude instead. We construct a one-sided chi-squared statistic, and obtain 95% confidence limits assuming one degree of freedom. Our limits agree with those presented by Laha et al. [55] in the n = 0 case, who instead ask that the predicted flux in every bin does not exceed the measurement by more than 2 times the reported error in that bin; using both methods, we have checked that our chi-squared approach yields identical results to the Laha et al. method, except above M ∼ 1.2 × 10 17 , where our constraints are stronger by a factor of a few. At lower masses, small differences with respect to the Laha et al. results can be attributed to a different choice of dark matter halo parameters. Our results are also similar to the very recent [57]. While their addition of Fermi and EGRET data may strengthen bounds at lower masses, they may still be superseded by the 511 keV bounds that we discuss next.
For the 511 keV signal, we may use more recent data. We have taken the binned 511 keV flux shown in Fig. 5a of Ref. [58] (black crosses). These correspond to the total 511 keV flux within galactic latitudes −10.5 • < b < 10.5 • , in 5 equally-spaced longitude bins within −30 • < ℓ < 30 • . We again produce a one-sided chi-squared, in order to establish 2σ limits on the BH fraction via Wilks' theorem. As in Ref. [59], we conservatively only consider positrons with energies less than 1 MeV, as high-energy particles may not annihilate in-situ.
Since our method slightly improves on previous results, we first show the resulting limits for the n = 0, ordinary 4D PBH case in blue, in Fig. 9. Continuum gamma-ray constraints are presented as solid lines, dash-dotted lines show the 511 keV limits from positron annihilation, and the dashed lines present the same limits, but without the E e + < 1 MeV requirement. We also show the aforementioned gamma-ray limits of Laha et al. [55] (solid yellow), as well limits based on evaporation to positrons obtained by Laha [60] and DeRocco & Graham [59]. When using the full range of positron energies, we attribute the slight improvement over DeRocco & Graham to the use of more recent data and angular information from [58]. The stronger improvement comes when comparing the E e + <MeV cases: here, our inclusion of secondary particles leads to a sizeable flux of low-energy positrons not present when only primary thermal particles are accounted for-as can be read e.g. from the left-hand panel of Fig. 8.
Constraints for n ≥ 0 are shown in Fig. 10. We arbitrarily cut the mass range to include BHs that would survive for at least 10 years, the approximate duration of the INTEGRAL mission (hence the large difference Fig. 9, which corresponds to BHs that would live for ∼ the age of the Universe or longer). Masses in the lower range are obviously "tuned" to end their lifetimes around the present day and correspond to a small sliver of initial BH masses. We will translate these constraints into cosmologically-consistent bounds in Sec. III E.
Depending on whether the Hawking temperature is high enough to produce positrons, and where the gamma-ray spectrum peaks, gamma-ray (left panel) and positron (right panel) constraints dominate for different values of M for different n. The sharp vertical jump at the right-hand side of some constraints corresponds to the transition from 4 + n-dimensional to 4-dimensional behaviour of the PBHs as they saturate the extra dimensions-i.e. masses above M = M 4D (6). We indicate with dashed lines the constraints that would be attainable in the absence of such a transition.
FIG. 9. Updated constraints on "ordinary" four-dimensional primordial black holes from INTEGRAL/SPI gamma-ray data. Solid blue: using gamma-ray continuum data from [54]; dashed blue: 511 keV line from e + e − annihilation, using data reported in [58]; dot-dashed blue: 511 keV constraints, but omitting the flux from positrons with energies higher than 1 MeV which may not annihilate in-situ. Prior results are shown from Laha et al.

B. Isotropic background light
The isotropic photon spectrum can be split into two observationally indistinguishable components. One component is the extragalactic background light (EBL) produced by extragalactic PBHs homogeneously distributed throughout the Universe. The EBL component has previously been used to constrain the abundance of extra-dimensional PBHs [39]. The other component is the isotropic part of the galactic signal, produced by PBHs within the galactic halo. Despite, the galactic halo being anisotropic (as discussed in the previous section), there is a non-zero flux in all directions. Therefore, there appears to be an isotropic component equivalent to the flux in the direction with the smallest contribution from galactic PBHs. This isotropic galactic signal has recently been used to constraint the abundance of long-lived four-dimensional PBHs [61,62].

Extragalactic photon flux
The sum of the evaporation products from all extragalactic PBHs could produce a significant isotropic flux of X-rays or gamma rays. This signal depends on the primary spectrum of photons, electrons and positrons described in Eq. (8) as well as the secondary spectrum described in Sec. II D. As the evaporation products travel from the point of evaporation to Earth, the flux changes due to the photons redshifting, being absorbed, and scattering with IGM material. By taking into account all of these processes, whose relative importance is a function of energy and redshift, we will obtain a predicted EBL flux that may be constrained by observations.
The EBL contribution to the isotropic photon flux can be found by evolving the photon spectrum over time starting at recombination. At any given redshift, z, the change in the flux of photons of energy E can be parameterized by where Φ γ,EBL is the extragalactic isotropic photon flux and i denotes the four different channels for energy injection and loss: Universe expansion, photon absorption, Compton scattering, and photon injection.
The expansion of the Universe redshifts photon energy and dilutes their number density. As shown in Appendix C, these effects may be combined into: This results in the flux per unit energy being diluted as (1 + z) 2 , as the photon number density is diluted as (1 + z) 3 while the spectral density removes a factor of (1 + z). Although Eq. (41) depends on the derivative of dΦ/dE, the discretized method that we use (Appendix D) does not actually require numerical differentiation. The processes that cause the absorption of photons are: photoionization of neutral gas, pair production from atoms and ions, photon-photon scattering, and pair production off the CMB. All of these processes either absorb a photon or remove almost all of a photon's energy. The change in photon flux due to these absorption processes is where dτ dz (E, z), as determined in [63], is the optical depth of a photon of energy E over a differential redshift step at redshift 1 + z.
Absorption of photons causes an initial flux of photons starting at redshift z i and travelling to a final redshift z f with final energy E to be suppressed by an exponential factor of e −τ (E,zi,z f ) where For photons with energies between ∼ 1 keV and ∼ 10 GeV the Universe is transparent (τ < 1) up to redshifts of order z ∼ 100. However, for photon fluxes that originate at higher redshifts, a large fraction of the photons may be absorbed.
High-energy photons can also Compton scatter with electrons, losing some amount of energy, without being entirely absorbed. The instantaneous change in photon flux due to Compton scattering is calculated as the sum of a negative loss term that accounts for the attenuation of photons of a given energy and a positive source term that accounts for all the higher-energy photons downscattered to that energy. This is given as where H(z) is the Hubble parameter, n e is the total electron density, which includes electrons bound in hydrogen and helium as the small ionization potentials do not distinguish those from free electrons (see e.g. [64,65] for more discussion), σ c is the total Compton cross section, and dσc(Ẽ) dE is the differential cross section of an incoming photon with energyẼ scattering and losing energy so that it ends up with an outgoing energy E.
Solving this integro-differential equation is computationally slow, and the effect of Compton scattering is often approximated either as an absorption process which contributes to Eq. (42) or as a process that causes all photons to continuously lose some fraction of their energy in a similar way to the expansion of the Universe. For scenarios where Compton scattering is important we utilize the full integro-differential equation. A discussion of the different computation schemes and more details on how Compton scattering was numerically calculated in this work can be found in Appendix D.
The differential Compton cross section is typically given in the rest frame of the electron in terms of the scattering angle θ by the Klein-Nishina equation whereas dσ c /dE is required to solve Eq. (44). Here, α is the fine-structure constant, m e is the electron mass, and the outgoing photon energy E is related to the incoming energyẼ and θ via The differential Compton cross section with respect to outgoing photon energy is thus The integration bounds in Eq. (44) are found by noting −1 ≤ cos θ ≤ 1 and translating that to a range ofẼ using Eq. (46). The total Compton cross section at a given energy, E, is [66] where σ T is the Thomson cross section and x = E/m e . Finally, photon injection from BH decay yields where n • is the black hole number density and dEdt is the spectrum of produced photons from a single black hole of mass M .
The photons are produced as primaries and secondaries directly from black hole evaporation, annihilation of positrons, and inverse Compton scattering (ICS) of high-energy electrons and positrons. Therefore, the rate of photon production per black hole can be split into: The photon production rate due to evaporation, d 2 N γ,evap /dEdt, is calculated as the sum of the photon greybody spectrum as expressed in Eq. (8) and the secondary photons produced by the annihilation of unstable massive particles as discussed in Sec. II D.
Sufficiently hot black holes also produce high-energy electrons and positrons. As these cool down, they yield additional X-rays by upscattering CMB photons via ICS. The production rate of photons due to ICS is given by the convolution of the electron and positron evaporation spectrum with the secondary photon spectrum produced by the cooling of a single electron or positron with a given energy. This can be expressed as where E is the photon energy, E e is the electron energy, T CMB is the CMB temperature, dEdt is the production rate of electrons from black hole evaporation, and dÑγ,ics dE is the secondary photon spectrum from a single electron or positron cooling down. The factor of 2 accounts for the fact that both electrons and positrons contribute to the ICS signal. The secondary photon spectrum from electron cooling was determined by interpolating a table calculated using DarkHistory [67].
After an energetic positron quickly loses most of its energy via ICS and other cooling processes, it will find a partner and annihilate to photons. First, positronium is formed in either the singlet or triplet state. One quarter of the positrons form positronium in the singlet (parapositronium, j = 0) state, which annihilates to two photons with E γ = m e . The remaining three quarters of the positrons form the triplet (orthopositronium, j = 1) state, which produces three photons with a spectrum first calculated in [68] and expressed in [67] as where x = E/m e and 0 ≤ x ≤ 1. Assuming 100% positronium formation, the photon yield per positron is thus Numerically, the Dirac delta function is modelled as a Gaussian with a width of 1 keV, which is a realistic approximation for the peak shape from galactic positronium annihilations [69]. Although Ref. [69] does not address extragalactic positron annihilation, cosmic expansion causes the integrated signal from all extragalactic annihilations to form a continuum below 511 keV. The resulting observed EBL flux is therefore insensitive to how the initial annihilation peak is parameterized. The production rate of photons due to positron annihilation can be found by multiplying Eq. (53) by the positron production rate dN e + /dt, including primaries and secondaries: Starting from recombination, the photon flux can be evolved forward in time using Eq. (40) to calculate the extragalactic contribution to the isotropic X-ray and gamma-ray spectrum today. Further details on how this equation was solved numerically are in Appendix D.

Galactic contribution
While the flux of evaporating black holes within the Milky Way halo would be highly anisotropic, because there is a non-zero flux in all directions, the flux in the direction that produces the smallest flux contributes an irreducible isotropic component on top of the extragalactic flux [61]. This flux can be calculated by evaluating Eq. (36) in the direction with the minimum flux, directly away from the galactic centre. Then, Eq. (36) simplifies to where f •,0 is the fraction of dark matter comprised of PBHs today, d 2 N γ /dEdt is calculated in the same way as in the EBL case except only accounting for evaporation to photons and positronium annihilation (the flux from ICS was not included in the galactic calculation), and D min is the integral of the Dark Matter density along the line of sight opposite to the galactic centre

Observational constraints
The total expected isotropic photon flux can be calculated by adding together the extragalactic contribution found by solving Eq. (40) and the galactic contribution from Eq. (55). That calculated photon flux was compared to measurements of the isotropic X-ray and gamma-ray signal compiled in [70]. The experiments included are, from lowest to highest energy: ASCA [71], RXTE [72], HEAO-1 [73] , HEAO-A4 [74], Swift/BAT [70], Nagoya [75], SMM [76], CGRO/COMPTEL [77], and CGRO/EGRET [78]. When the widths of the energy bins was not provided, it was assumed that bin widths extended to the midpoint with neighbouring bins. Although a measurement from the instruments on INTEGRAL (JEM-X, IBIS, SPI) are available [79] we do not include them, as they are less precise, and overlap with other data used here. The observed fluxes as well as a sample calculated spectrum are shown in Fig.  11.
To account for sharp features such as the 511 keV peak from Milky Way positronium annihilations, the calculated flux was averaged over each bin width to determine the expected flux for each experiment. Constraints were then set by ensuring that the expected flux does not exceed the observed flux by more than 2σ in any energy bin. This approach leads to conservative constraints on the PBH abundance because no assumptions are made about other astrophysical sources of X-rays and gamma-rays. Including models of astrophysical X-ray and gamma-ray sources can currently strengthen PBH constraints by more than an order of magnitude [61,62] and have an even larger effect when projecting the discovery potential of future X-ray telescopes [80]. Constraints from isotropic background light are shown in Fig. 12. The shapes of the n = 0 and n = 2-6 constraints are generally similar. The low mass cutoff of the constraints is given by the black hole mass that leads to evaporation before the time of recombination (taken to be z = 1100), as photons from these BHs cannot propagate freely until today. At slightly higher masses, BHs evaporate completely between recombination and today. The largest signal comes from the high-temperature emission at the end of their lives; more massive black holes evaporate closer to today such that their emitted photon spectrum has redshifted less, and the observed spectrum has a higher energy, where observed fluxes are lower. This leads to constraints strengthening with increasing initial BH mass. This trend continues until the black holes are massive enough to survive until today. Beyond this point, more massive black holes have lower temperatures and there are fewer black holes for a given energy density, causing the trend to reverse.
For n = 4, 5, and 6, as the mass increases, limits weaken sharply as the BHs Schwartzschild radius exceeds the size of the extra dimensions (as in Eq. 6), leading them to mimic the n = 0 limits. Since this transition depends on the details of the compactification, the true behaviour would not be as sharp.
We found that there are no BH masses where the inclusion of photons produced by inverse Compton scattering improves the constraints. As shown in Fig. 13, the galactic isotropic flux strengthens the constraints set on black holes that survive until today and including the flux from positron annihilations increases the strength of the constraints for black holes with temperatures close to 511 keV.

C. Cosmic microwave background
Evaporation of primordial black holes during and after recombination can lead to high-energy electrons and photons producing heating and ionization-an effect first discussed in the context of decaying heavy neutrinos [81] and later adapted to annihilating dark matter [64]. A higher ionization floor will rescatter CMB photons. During the dark ages, this has the effect of "blurring" the last scattering surface (LSS), suppressing the angular power spectrum on small scales (large ℓ). For ionization at lower redshifts, this rescattering additionally enhances power at lower multipoles in the EE polarization power spectrum, because Thomson scattering is polarized [82].
As part of the CosmoLED package, we modify the public ExoCLASS code [50], a branch of the CLASS linear anisotropy solver [83] which deals with the energy injection from WIMPs or primordial black holes. To be specific, we change the DarkAges module to incorporate LED BHs with n =1-6 and a flexible Planck scale M ⋆ . 4D BH remains a choice when n is set to 0. The electron and gamma spectrum required for the module is now computed as described in Sec. II D. We improve ExoCLASS in the following aspects: 1) We implement the complete greybody spectrum for all particles, instead of cutting the spectrum at E = 3T H and approximate the absorption cross section as σ s = 27πG 2 M 2 • . 2) We include the secondary particles from primary particles at energies above 5 GeV using the PPPC4DMID tables. 3) At low energies, we use Hazma and our own code to calculate the decay of pions and muons as a function of particle energy, instead of using the fixed decay table in ExoCLASS. A comparison of secondary particle spectra from CosmoLED and ExoCLASS can be found in Fig. 8. We have also altered the black hole mass evolution of a function of time in DarkAges module and CLASS main code. Apart from these changes, we follow the approaches in ExoCLASS to compute the energy deposition from LED BHs, which we briefly summarize below.
The injection of energy from decaying black holes with initial mass M i and initial fraction f • , relevant for CMB observation is given by where f e.m. is the fraction of BH evaporation that ends up with e ± and γ, and ρ c Ω CDM is the cold dark matter energy density today. In CosmoLED, this is computed from with the right hand side given by Eq. (33). The injected energy is then deposited at different redshift z, in the form of ionization, excitation of the Lyman-α transition and heating of the intergalactic medium. The energy deposition is therefore connected to the energy injection by and the energy deposition functions in the three channels denoted by h c can be obtained by convolving the injected electromagetic particle spectra with a transfer function that models streaming and absorption of electromagnetic products in the high-redshift IGM. We follow the treatment in ExoCLASS and employ the transfer functions precomputed in Refs. [84,85].
To constrain the initial fraction of BHs in dark matter, we use MontePython [86,87] to run a Markov Chain Monte Carlo (MCMC), which interfaces with the modified version of ExoCLASS in CosmoLED. For each PBH initial mass and n, we impose flat prior on the initial fraction of BHs, and six ΛCDM parameters {ω b , ω cdm , θ s , ln(10 10 A s ), n s , τ reio }. We adopt the Planck high-l TT,TE,EE+low l TT, EE+Planck lensing 2018 [1] likelihoods, with standard Planck nuisance parameters marginalized over. Fig. 14 shows the boundary of each 95% one-dimensional credible interval on the initial fraction of PBHs as DM, f • , as a function of the initial black hole mass M . For each n, the excluded region cuts off abruptly at low mass, where BH evaporation occurs before recombination and thus does not affect the ionization floor. The cutoff of n = 6 BHs coincides with that of 4D BHs, as from Fig. 2 M ≃ 10 14 g BHs disappear at CMB in both cases. As the mass increases, sensitivity is gradually reduced as the Hawking temperature of the PBH population falls with mass.
Even though our inclusion of secondary particles leads to a larger γ and e ± flux than in the default ExoCLASS scenario, our constraints for n = 0 are slightly weaker than those presented in Ref. [50]. These differences may be due to their implementation of a prior on τ reio , with which f • is degenerate, or the use of different Planck data sets.  testing ground for physics beyond the Standard Model. Constraints on mechanisms that modify either the expansion rate or balance of the synthesis processes during this era have been explored previously in, for example Refs. [88][89][90][91][92][93][94][95][96].
In a similar spirit, the presence and evaporation of black holes leading up to, during, and beyond BBN, can impact the resulting relic abundances in a number of ways. Weak interactions freeze out around temperatures of ∼ 1 MeV, just before the onset of BBN, setting the neutron-to-proton ratio which is critical to the eventual formation of helium. An additional black hole density component may alter the expansion history of the Universe and the subsequent freeze-out of this ratio. More specifically, an increase in the expansion rate will lead to an earlier weak interaction freeze-out, an enhanced neutron-proton ratio and eventually, a greater helium-4 abundance (see Sec. III D 2 for further discussion.) Black hole evaporation products, namely pions, may also alter the neutron-proton fraction after freeze-out via direct conversion. In addition, if the temperature of the black holes is sufficiently high, the resulting evaporation products will be able to directly contribute to the dissociation of the forming nuclei.
In order to incorporate black holes and their evaporation products correctly into the relic calculation, a complex system of reactions needs to be solved self-consistently. As most public codes do not allow for non-thermal energy injection, we will deal with these two effects separately. 3 In Sec. III D 1 we recast prior results following the method of Ref. [100]; in Sec.III D 2 we adapt the AlterBBN code to produce the light abundances from the appropriately modified expansion histories.

Photo-and hadrodissociation
If the bulk of BH evaporation occurs during or shortly after BBN, the production of high-energy particles can lead to dissociation of nuclei, affecting the relic abundance of D, He and Li. The addition of a non-thermal component to existing BBN codes is non-trivial. Kawasaki et al. [96] performed a detailed numerical analysis, deriving constraints on the lifetime of decaying dark matter during the BBN epoch as a function of its mass and density. They utilized updated reaction rates, newly implemented interconversion of energetic protons and neutrons by inelastic scattering off background nuclei, as well as the incorporation of energetic antiprotons and antineutrons. Their results use the observed relic abundance of light elements, including the primordial mass fraction of 4 He, Y p ≡ ρ( 4 He)/ρ b = 0.2449 ± 0.0040 [101], the primordial deuterium to hydrogen ratio (D/H) p = (2.53 ± 0.04) × 10 −5 [102] and the upper limit on the primordial 3 He to deuterium ratio ( 3 He/D) p < 0.83 + 0.27 [103]. Keith et al. [100] pointed out that evaporating black holes modify BBN abundances in a similar manner to decaying massive particles and recast the results of Kawasaki et al. to derive equivalent constraints for black holes. We will mostly follow the procedure outlined in Ref. [100] to recast the results in Ref. [96] for the LED BHs described in this article. The method, assumptions, limitations and results are presented below.
Ref. [100] broadly distinguishes between two phases of nuclear dissociation due to BH evaporation products: the hadrodissociation era at high plasma temperatures, and the photodissociation era at later times. Both of them lead to the dissociation of 4 He and the production of D and 3 He. We follow the same approach as Ref. [100] to account for the photodissociation of 4 He caused by BH evaporation, while for hadrodissociation, we adopt a different procedure which better captures the total number of hadrons injected by BHs. In both cases, we use the precise greybody spectrum to compute the average quark energy, instead of assuming a thermal Fermi-Dirac distribution.
If decays happen at late enough times, when the plasma temperature is lower than T ≲ 0.4 keV, all electromagnetic final states contribute to dissociation. Because a majority of SM degrees of freedom-and thus evaporation productsare in the hadronic sector, this can be mapped to previous bounds on dark matter decay to quark-antiquark pairs. Neglecting the quark masses and averaging over the quark greybody spectrum, the mean energy ⟨E q ⟩ M for a given BH mass is where d 2 N /dEdt is the radiated quark energy distribution, given by Eq. (33). Since quarks are typically produced above the QCD transition scale, the mean quark energy is obtained by averaging the emission over the lifetime of a BH when the Hawking temperature is high enough, i.e.
where dN/dM is the number of quarks produced per change in BH mass, and can be inferred from Eq. (33) (after integrating over E) and Eq. (11) considering quarks and gluons. T H,i is the initial BH Hawking temperature. The total energy injection, which is relevant for photodissociation of 4 He, of a BH with initial mass M i , should thus yield a similar effect to the decay of DM particles with mass M X ≃ 2⟨E q ⟩ into quark pairs. The step function in Eq. (61) ensures that quarks are not produced below the QCD scale. This approach is conservative, in that it ignores evaporation for Hawking temperatures below Λ QCD to other states.
At higher temperatures (T ≳ 0.4 keV), e + e − pair production from photons is efficient, and the dissociation of 4 He primarily expected to be from hadrons produced from quark and gluon jets, which builds up with the injection of more hadrons. The number of hadrons in a quark jet scales as E 0.3 q , and therefore on average, the number of quarks produced from the greybody spectrum is approximated by quarks with a single energy ⟨E h ⟩ M which satisfies Again averaging over the evaporation lifetime of a BH, the number of hadrons per unit energy, proportional to E 0.3 q /E q , is computed as The numerator gives the total number of hadrons emitted during the lifetime of a BH, and the denominator shows the total hadronic energy. This can again be mapped to dark matter which decays to quark-antiquark pairs, with the number of hadrons per unit energy given by (M X /2) −0.7 . Therefore, we have the relation The values of ⟨E q ⟩ M and ⟨E h ⟩ M , as well as the k q and k h coefficients are computed and listed in Table III. Note that these differ from values presented in Ref. [100] as we use the full greybody spectra to model the quark phase space distributions, and a different method for hadrodissociation.
To obtain the constraints on the enregy density of BHs, we find the correspondence between BHs and decaying dark matter that causes the same amount of dissociation to light elements. Conservatively we only consider the hadrons and photons produced from quarks and gluons, not other particles. If BHs initially have a Hawking temperature above the QCD transition scale, i.e. M i < M QCD (T H = Λ QCD ), the entire mass of BHs is injected to the plasma in the form of quarks (and gluons), up to an order 1 number f q which quantifies the fraction of hadronic injection. Therefore, roughly the same amount of quarks are produced in BH evaporation and dark matter decay, provided that they start from the same energy density. However, if M i > M QCD , quarks are only emitted when BH mass reduces to M QCD , and the early stage of the BH mass dump does not dissociate any nuclei. To match the number of quarks injected, the initial fraction of BHs f • that we constrain is related to the fraction of dark matter made of decaying particles, f X constrained by Kawasaki by 4D BHs always have M i < M QCD in the relevant mapping mass range. However, LED BHs can have longer lifetimes and lower Hawking temperatures, rendering the M i /M QCD factor important. The fraction of hadronic energy injection f q mildly depends on BH mass, running from 76% for 4D BHs, to 65% for n = 6 BHs, due to differences in the greybody spectra, as well as a growing fraction of graviton emission.
To complete the translation of constraints from decaying dark matter to BH evaporation, we must determine the appropriate correspondence between the lifetime of BHs τ • and dark matter decay time τ X . While we expect that τ X ≃ τ • , these processes are fundamentally different in that DM decay represents a steady injection of energetic particles, while BH evaporation products increase in energy until a dramatic spike at τ • , after which no BHs remain. As done by Keith et al. in Ref. [100], we match BHs and decaying dark matter at a time when half of the energy is injected. For decaying dark matter, this happens at a time t = τ X ln 2. For BHs with initial mass M i < M QCD , injecting half of the total energy takes the time t = f t τ • , and f t ≃ 0.5. This yields the relation τ X = f t τ • / ln 2. If however M i ≫ M QCD , the lifetime of M QCD BH is negligibly small, and we have τ X = τ • / ln 2 instead.  Our constraints are presented in Fig. 15. For BH with mass M , we find the dark matter lifetime that matches BH lifetime, and the dark matter mass M X that reproduces the dissociation effects of a BH, using the method outlined above. We then interpolate the constraint lines in Ref. [96] according to M X , using X → uū decay channel. The interpolation works well for 4D BHs. However, for LED BHs, the corresponding dark matter mass is below the smallest mass considered of 0.03 TeV in most of the parameter space due to the low Hawking temperature. Noting that the constraints on the energy density of dark matter get stronger for lighter dark matter mass, as hadrodissociation depends on the number of emitted hadrons proportional to M 0.3 X , and photodissociation is roughly determined by the total energy injection. We therefore use the M X = 0.03 TeV constraint line for any mapped dark matter mass below 0.03 TeV, to produce a conservative bound on the energy density of BHs. We present results in terms of the initial fraction of dark matter made up of black holes, f • . This can be equated to β ≡ ρ • /ρ tot and M Y , the decaying particle mass times their number density per unit entropy, used in Ref. [100] and Ref. [96] respectively, via where T 0 and T form are the CMB temperatures today and the plasma temperature at black hole formation respectively. As in previous figures, red, yellow, purple, green and light blue curves correspond to the n = 2 − 6 extra dimensional cases respectively. The rightmost dark blue curve shows the 4D results, which are well-matched to those derived in Ref. [100], though the inclusion of the relevant greybody factors and the updated method leads to some small differences at lower masses. The different mass range covered by the LED BHs also leads to a number of qualitative modifications of the 4D results. As seen in Table III, the maximum 4D BH mass translated from decaying dark matter is below M QCD . However, for any n ≥ 2 and M ⋆ = 10 TeV, in some part of the mass ranges BHs have initial Hawking temperatures that fall below the QCD transition scale. The correction due to M > M QCD is more pronounced for lower number of extra dimensions, and starts to severely restrict the parameter space that can be constrained above about 10 11 g for n = 6 BHs. For all n ≥ 2, this accounts for the f • ∝ M loss of sensitivity at higher masses. There are a number of assumptions underwriting the validity of this methodology. They mostly pertain to being able to match both the spectral and the temporal distribution of the injected energy from an evaporating BH to that of a decaying particle.
Firstly, it is assumed that the spectral shape does not significantly vary the impact on BBN, provided the average energy of the injected particles is the same. Similarly, the temporal spread of injected energy from BHs can be treated as equivalent to that of a decaying particle, as long the averaged energy is injected at approximately the same time. Keith et al. note that the spread of particle energy around the mean for the 4D case could lead to an error of around 30% for BHs evaporating after ≈ 10 7 s. The effect is larger for BHs with shorter lifetimes where errors of up to a factor of 2 are possible.

Altered expansion history
The method described above accounts for the catastrophic injection of nonthermal energy during or after nucleosynthesis leading to nuclear dissociation. In addition to this effect, the presence of extra matter in the form of black holes during BBN, as well as the smooth injection of entropy leads to an altered expansion history, baryon-to-photon ratio, and ratio of neutrino-to-plasma temperatures, which all contribute to altering the freeze-out abundances of the primordial elements. It will turn out that only the former effect has an impact on nucleosynthesis. We treat these effects separately from the dissociation discussed above, as it pertains to a slightly earlier epoch-and publicly available software allows for a more exact treatment. We modify AlterBBN [104,105] to include BHs as additional species. In the code, BHs alter BBN in two ways: 1) the energy density of BHs contributes to the expansion of the Universe, and 2) BHs dump entropy to the plasma, increasing the temperature of photons and neutrinos. As described above, only effect 1) will turn out to be constraining, though these constraints will be subdominant to those presented in Sec. III D 1. Details of the implementation in the code, and resulting constraints, are described below.
With BHs, the energy density of the Universe during BBN is given by where we include the energy density of photons, neutrinos, electron and positrons, baryons and BHs. The energy density of e ± is connected to the photon temperature, parametrized with a series of Bessel functions [105]. The baryon density is fixed by the baryon-to-photon today, and we assume η 0 = 6.1 × 10 −10 . We start evolving the code from the neutrino decoupling temperature T dec = 2.33 MeV. We assume neutrinos and photons are in thermal equilibrium separately with the temperature T γ and T ν after neutrino decoupling. The neutrino energy density ρ ν = 7π 2 120 N ν T 4 ν , where we fix N ν = 3.046, and the photon energy density ρ γ = π 2 15 T 4 γ . For each decoupled species, the continuity equation implies where ρ, P and s are the energy density, pressure and entropy of the species. BH evaporation will dump entropy to the plasma containing photon, baryons and electrons, as well as the neutrino sector. We assume the two sectors keep thermal equilibrium efficiently and separately. The net effect of the entropy dump is to raise the temperature of the plasma and neutrinos, which in turn increases the energy density of photons and neutrinos through the expressions given above. Considering BH evaporation, the neutrino entropy follows and the plasma entropy which are employed in Eq. (68) to determine the evolution of the plasma and neutrino temperatures. With AlterBBN, we compute the abundances of He and deuterium and confront them with observations. We use the most updated primordial deuterium to hydrogen abundance ratio in PDG 2020 [106] where (D/H) p = (2.547 ± 0.025) × 10 −5 , which reflects the weighted average of the most precise measurements. Similarly, the primordial 4 He abundance is determined to be Y p ≡ ρ( 4 He)/ρ b = 0.245 ± 0.003 . The numbers are slightly different from the abundances quoted in Ref. [96], but the results remain robust regardless of the subtleties. To have a sizeable effect, the BH abundance must be ∼ 10 −3 or larger at BBN, which is significantly larger than the matter density expected during BBN, i.e. f • ≫ 1. We thus present the constraints on the fraction of BH energy density in the Universe at the neutrino decoupling temperature, and show the 2σ limits on β dec Fig. 16. To the left of the lines BHs evaporate significantly before the plasma temperature drops below T de . If BHs survive through BBN, a BH fraction as low as 10 −3 will modify the expansion of the Universe substantially, resulting in He and D abundances that are inconsistent with observation. The same bound holds for higher mass BHs which barely evaporate during BBN, but becomes weaker for lighter BHs that disappear before BBN ends.

E. Combined constraints
The combined constraints on the initial fraction of dark matter in the form of PBHs are shown in Fig. 17 for M ⋆ = 10 TeV, and n = 2, 3, 4, 5 and 6 extra dimensions as well as for the regular 4D scenario, where M ⋆ ≡ M pl , denoted n = 0. Features are qualitatively similar for different n.
At low masses, rapid evaporation leads to excessive injection of high-energy hadrons and photons during and after BBN. at higher masses (≳ 10 8 -10 14 g), longer lifetimes allow BH decay to take place after recombination, leading to strong constraints from the rescattering of CMB photons on the higher ionization floor. BHs that survive longer still produce an isotropic extragalactic signal, as well as a flux of gamma rays from the Milky Way halo. As n rises from 2 to 6, the Galactic flux of gamma rays moves into, and then out of, the INTEGRAL/SPI energy window, explaining how the isotropic background light and Galactic constraints trade places as the dominant limits with varying n.
For n ≥ 5 and the 4D case of n = 0, there is a small gap in the combined constraints between the BBN constraints on low-mass PBHs and the constraints on PBHs that survive until after recombination. This gap is due to the limited range over which decaying DM BBN bounds can be recast as PBH constraints. It is expected that the BBN constraints could be extended to higher-masses, closing the gap, by directly calculating the primordial nuclei abundances in the presence of LED PBHs.
Similarly, for 2 ≤ n ≤ 4, the BBN bounds weaken significantly for black holes which evaporate completely between the time of BBN and recombination. This creates a window of PBH masses for which the constraints allow f • ≫ 1 corresponding to a scenario where in the early universe dark matter is dominated by PBHs but those PBHs completely evaporate leaving the stable dark matter density observed today. This window of weak constraints is due in part to our conservative approach to estimating the photodissociation of light nuclei expected from PBHs with temperatures greater than the QCD scale as discussed in Sec. III D 1. A dedicated calculation of the effect of PBHs on primordial abundances may be able to set stronger constraints in this mass range and to some extent close this weakly constrained window.
The galactic centre constraints in Fig. 17 have a somewhat different shape and cover a smaller mass range compared to those presented in Fig. 10. This is because Fig. 10 shows the constraints in terms of the PBH mass and abundance today, whereas Fig. 17 shows the constraints in terms of the PBH properties before evaporation occurred. For massive PBHs that have only evaporated a negligible fraction of their mass since their formation, their relative abundance is unchanged since the early Universe so f •,0 ≈ f • . Therefore, for larger masses, the galactic centre constraints in Fig.  10 and Fig. 17 look the same. However, the smaller masses presented in Fig. 10 correspond to a very narrow range of initial PBH masses in Fig. 17 for which the PBHs just happen to be in the final stages of evaporation today. This leads to the leftmost end of the galactic constraints being "compressed" in Fig. 17.
The highest mass PBHs in the range presented in Fig. 17 are constrained by microlensing of stars in the M31 galaxy [107]. These constraints, shown in grey, are the same for all values of n and therefore were not recalculated in this work because the presence of LEDs would not affect microlensing. The mass range between the grey microlensing constraints and the coloured evaporation constraints is completely unconstrained so long as PBHs are not more abundant than dark matter, corresponding to f • ≤ 1. The region of parameter space where PBHs survive until today and would be more abundant than the observed dark matter abundance is shown in Fig. 17 by the grey hatched region.

Relic abundance of LED black holes as dark matter
For n ≥ 2, the solid black lines in Fig. 17 show the predicted abundance and mass of PBHs produced by energetic collisions in the early Universe assuming that M ⋆ = 10 TeV. These lines appear vertical because they have an extremely steep slope, where each point along the line corresponds to the expected mass and fraction obtained by varying the reheating temperature, T RH . As shown in Fig. 7, a small change in T RH corresponds to a very large change in the abundance of PBHs. We solve for the energy density of BHs as described in Sec. II C, and define the initial fraction f • in Eq. (35) at the time when BH mass is largest under accretion and evaporation. Due to the strong dependence of abundance on T RH , and a comparatively weak dependence of PBH mass on T RH , the predicted abundance lines are very steep. Therefore, for a given number of LEDs, there is a narrow predicted mass range for PBHs that would be produced with a fixed M ⋆ .
In the case of n = 2, the PBHs produced are sufficiently heavy that they would only have evaporated a negligible fraction of their mass. These surviving PBHs can comprise all of the dark matter. This scenario where n = 2 and the fraction of dark matter made up of PBHs, f • = 1, would correspond to BHs with a mass of ∼ 10 21 g. These BHs are too heavy to be constrained by evaporation bounds and lighter than any lensing constraints, making them a viable unconstrained dark matter candidate.
PBHs produced in theories with n > 2 and M ⋆ = 10 TeV would have entirely evaporated before today and are therefore not dark matter candidates. However, the narrow mass window does make specific predictions about when they finish evaporating, pointing at their most promising paths to discovery. For n = 3, the PBHs finish evaporating after recombination, so that the most likely cosmological method of discovering them is from their impact on CMB anisotropies or through the isotropic X-ray and gamma ray signal they produce. PBHs with n = 4 or 5 complete their evaporation earlier, before recombination such that their most apparent cosmological imprint would stem from the destruction of primordial nuclei formed during BBN. For scenarios with n = 6, the PBHs complete evaporating so early that they would be entirely gone before BBN begins. This makes it very difficult to constrain the existence of n = 6 PBHs with M ⋆ = 10 TeV. However, the possibility of a very large abundance of PBHs forming and evaporating to all particle types in the very early Universe raises the intriguing scenario that the PBHs may have evaporated to stable dark matter particles, yielding a non-thermal (in the cosmological sense) relic abundance production mechanism. Evaporation to dark matter particles can be incorporated into any of the scenarios with n ≥ 3, but the short hot lifespan of n = 6 makes them especially interesting scenarios to explore in future work.

Comparison with prior four-dimensional results
Many previous studies of 4D PBHs (see Ref. [5] for a review) set constraints on the abundance of PBHs not in terms of f • , but instead in terms of β ′ which is defined as To make comparisons between the PBH constraints computed in this and previous work simpler, the 4D PBH constraints are shown in terms of β ′ in Fig. 18. The combined constraints in Fig. 18 also includes the BBN constraints due to changes in the expansion of the Universe as shown in Fig. 16, converted from β dec using Eq. (66). The grey shaded region in Fig. 18, shows a selection of the strongest constraints on low-mass 4D PBHs from previous work. This region combines constraints set with BBN [100], CMB anisotropies [50], isotropic photons [62], galactic centre photons [57], and galactic centre positron annihilation [59]. These are generally very similar to the strongest constraints set in this work, although there are a few differences worth noting. As discussed in Sec.III A, the 4D constraints we have set using positron annihilation in the galactic centre are stronger than those previously set in Ref. [59]. It should also be noted that Refs. [61,62] have set stronger constraints using the isotropic X-ray and gamma ray flux by modelling astrophysical sources. However, these are dependent on the astrophysical source model used, though our results are stronger than the conservative background-agnostic constraints of Ref. [62] The isotropic background light bounds set here are stronger than those in Ref. [62] for lower mass PBHs and weaker for higher mass PBHs. For lighter PBHs that would have completely evaporated this difference is driven by different approaches in calculating the secondary spectrum of photons from unstable evaporation products. For PBHs that survive until today, the difference is driven by differing assumptions for the parameterization of the Milky Way halo. Finally, the CMB constraints due to energy injection during the dark ages are a factor of ∼ 6 weaker than those presented by the authors of ExoCLASS in Ref. [50]. However, even with a fresh installation of ExoCLASS we were unable to exactly reproduce their results-our inclusion of more precise secondary spectra yields a factor of 2 improvement over constraints found using the public code as-is. This discrepancy is possibly attributable to the updates to ExoCLASS since Ref. [50] was published, or a different choice of priors or nuisance parameters.

F. Other Constraints
Previous analyses have set constraints on the existence of light 4D PBHs using more methods than we have employed in this article. In this section we discuss some of those constraints and whether they are expected to be important for the study of LED black holes.
Positrons directly injected into the interstellar medium (ISM) from BH evaporation can contribute to the local cosmic ray flux. Since these are predominantly at sub-GeV energies, they are strongly affected by solar modulation and associated uncertainties. Ref. [108] placed constraints on 4D PBH evaporation for M ∼ 10 14 -10 17 g, using data from the Voyager I spacecraft, which has recently crossed the heliopause. These are subdominant to the more recent constraints from gamma ray emission using INTEGRAL data derived by Ref. [55]. Since our INTEGRAL/SPI galactic constraints use the same dataset as Ref. [55], we anticipate that the Voyager I constraints would be similarly subdominant in the LED scenario.
Dwarf spheroidal galaxies are a prime target for gamma ray searches for dark matter decay or annihilation signatures thanks to their high mass-to-light ratio, which implies a low standard model background and a large prospective signal. Ref. [109] recently analyzed ∼ 1 Ms of observations of the galaxy Reticulum II with INTEGRAL/SPI over energies 25-8000 keV. Though this leads to improved limits on DM annihilation, the resulting limits on PBH decay are weaker than galactic centre analyses.
Radio data from the inner Galactic Centre have been used to constrain 4D PBHs [110]. Large magnetic fields cause ultrarelativistic electrons and positrons to cool via synchrotron radiation thus producing an observable radio signal. In the case of LED PBHs this is most likely to be a viable observational channel for n = 5 where PBHs that survive until today can be hot enough to produce ultrarelativistic electrons. However, constraints on 4D PBHs from radio data are always weaker than constraints based on X-ray and gamma-ray observations and therefore including radio data in this analysis is unlikely to improve the constraints we have set on LED PBHs.
Ultra-light PBHs could dominate the very early Universe and entirely evaporate before BBN evading all bounds presented in this work. However, these PBHs and associated curvature perturbations could source a measurable stochastic gravitational wave background (SGWB) [111]. Recently, that SGWB has been used to produce constraint forecasts for future space-based gravitational wave interferometers [112]. Due to the different lifetime and production mechanism of LED PBHs, these forecasts must be recomputed for the case of LEDs. Some of these constraints would apply to PBHs with masses lower than the BBN constraints presented here.
The evaporation of PBHs during the epoch of star formation and reionization could leave imprints in the highredshift 21cm signal by heating and ionizing intergalactic gas. Several studies have presented current or future limits, considering the evaporation of 4-dimensional PBHs, either motivated by the recent detection of a deep 21cm absorption trough by the EDGES [113] experiment [114][115][116][117][118][119] or looking ahead to large-scale experiments such as the Square Kilometer Array [120]. Many other studies have examined the impact of matter accretion onto macroscopic PBHs that might seed early structure formation or produce X-ray backgrounds [121][122][123][124][125]. These studies and others highlight the potential for future high-redshift 21cm observations to be highly constraining of exotic energy injection sources during the Dark Ages and the epoch of Cosmic Dawn. We expect LED PBHs may similarly have a strong impact on future 21cm observables.
A bound on PBH evaporation in the galaxy was recently placed based on measurements of the ISM temperature in the Leo T dwarf galaxy. [126]. These require careful accounting of heating and cooling effects in the ISM -based on the results of [126], which are stronger than the INTEGRAL constraints of [55] between 1 and 3 × 10 17 g, they could lead to improved limits in a small portion of the parameters space for LED BHs.  , CMB anisotropies [50], isotropic background light [62], galactic centre photons [57], and positron annihilations in the galactic centre [59].
Finally, if the compactified extra dimensions have a toroidal geometry, the production and subsequent decay of Kaluza-Klein (KK) modes during reheating sets constraints such that any reheating temperature that would result in PBHs forming would already be severely constrained [127]. However, constraints based on the production and decay of KK modes are highly dependent on the compactification geometry, the decay products and the existence of additional branes [127]. Conversely, the PBH results in this work are only sensitive to the precise compactification geometry when r h ∼ R (or alternatively stated as M ∼ M 4D ) and the results for all other values of M are insensitive to such details. This makes observational constraints based on KK mode production and PBH production complementary to each other.

IV. CONCLUSIONS
In this article we have derived the full cosmological evolution of PBHs in the presence of LEDs including their production, accretion and evaporation history. We then derived bounds on the existence of those low mass PBHs using BBN, CMB, isotropic photon flux, and galactic centre X-rays. In doing so, we have also recomputed or updated the constraints on 4D PBHs from each of those sources. The constraints on heavier PBHs from gravitational lensing, mergers, and accretion rely on physics at scales larger than the size of the LEDs and therefore will be unchanged from previous analysis.
The abundance and mass of the PBHs for a given number of extra dimensions depend on M ⋆ and T RH . We have set conservative constraints on the M ⋆ -T RH parameter space by ensuring that the PBHs are not overabundant. Stricter constraints could be set on the properties of the extra dimensions by ensuring the produced PBHs are not ruled out by the astrophysical constraints derived here. To do so would require recomputing the astrophysical constraints over a full range of M ⋆ values and has been left to future work.
We have also found that in the case of two LEDs, the PBHs produced in the early Universe would survive until today and could, with the appropriate reheating temperature, comprise the entirety of the dark matter abundance observed today. In the cases of n > 2, PBHs would still be created in the early Universe however they would be light enough such that they would have evaporated before today. In those cases the PBHs could still have interesting cosmological impacts even if they are not a dark matter candidate.
In addition to their prospect as dark matter candidates, black holes can produce all gravity-coupled degrees of freedom as they evaporate, as long as the BH temperature is high enough, and the particle mass is kinematically accessible. In the case of BHs produced at colliders, this provides a potential window into the dark sector [128]. PBHs produced in the early Universe could also evaporate to yield the relic abundance of dark matter [129][130][131][132][133][134][135]; this behaviour would change in the presence of extra dimensions.
The possibility of large extra dimensions opens a new direction in the search for primordial black holes, including the alluring possibility of producing PBH dark matter without relying on large or non-Gaussian primordial fluctuations. The distinctive evaporation rate and spectra of these BHs mean that any positive detection would point directly at the existence of higher spatial dimensions and provide tantalizing clues about the origin of the Planck scale, bringing together two of the deepest mysteries of the cosmos: dark matter, and the unification of gravity with particle physics.
If radiation dominates, dT /dt is given in Eq. (24) and After t i the right hand side of Eq. (A2) is vanishes and the mass spectrum drops as a −3 ∝ T 3 in a radiation dominated universe.
Exact solutions. Now we turn to a more rigorous treatment without assuming instantaneous accretion. At time t i , the number of microscopic BHs produced is h t ≡ dn• dt | t=ti . Since the accreted BH mass is fairly insensitive to the initial masses, we assume all BHs are created at a minimum mass required for efficient accretion M i = M i,min (T i ), defined in Eq. (23), and they evolve collectively afterwards. We use M (t i , t) to denote the mass of BHs that evolve from t i to t, and h t (t i , t) to show the evolution of the BH mass spectrum. The latter is described by Eq. (A5) can be further split into two equations, one for BH production at t i and the other for the redshift of the spectrum at t > t i ∂h t ∂t The evolution of BH mass follows Eq. (21), which reads At time t, the BH energy density is given by where t RH is the time of reheating, and the radiation density evolves as where the loss or gain of radiation is caused by the change in BH mass. The second term on the right-hand-side of the equation can usually be neglected since the accretion energy loss is supposed to be much more efficient than Planckian mass BH productions. Combining BH and radiation, the expansion of the Universe is governed by the Friedmann equation where ρ r is given in Eq. (20). Eqs. (A6) to (A11) provide a complete set of integro-differential equations to solve for the BH mass M (t i , t) and mass spectrum h t (t i , t).
Instead of solving the above equations directly, we can reduce the number of equations by switching to the temperature basis, where we find where h T (T i , T ) ≡ dn dT | T =Ti and Eq (A11) stays the same. We can substitute Eqs. (20) and (A15) into Eq. (A16) to find dT /dt Since the variation of g ⋆ is mild during accretion, we set dg ⋆ /dt = 0. Eq. (A17) can further be plugged into Eqs. (A12-A14) to obtain the equations for h T and M , where Eq. (A6) and M (T i , T = T i ) = M i,min (T i ) serve as initial conditions. However, this formalism may not work if BHs dominate the energy density of the Universe after accretion and reheat the plasma significantly when they decay. It is crucial to solving the equations on the time basis to keep track of BH evolution in this scenario.
In the above integro-differential equations, T and t always appear in the differentials while T i and t i appear in integrals. After we obtain h T and M , we can map it to the mass spectrum using the relation We solve the full integro-differential equations on time basis, and show the evolution of BH energy density in Fig. 19, as well as the mass spectrum at 10 −10 s in Fig. 20. We choose two typical scenarios. In the first scenario, n = 2, T RH = 1.09 TeV, the BH energy density remains subdominant until the plasma temperature drops below 0.75 eV. In the other scenario, BHs dominate the energy budget of the Universe at about 10 −15 s, and then evaporate away before BBN. In both cases, we find good agreement with the energy density evolution obtained from the single mass approximation described in Sec. II C. The mass distribution spreads more in the full solution. However, the peak mass of the distribution agrees with single mass approximation, up to a close to 1 factor. Appendix B: Secondary particle production from pion and muon decay We use Hazma [49] code to compute the secondary gamma from π 0 decay π 0 → γγ, and from radiative muon and charged pion decay through the processes µ − → e −ν e ν µ γ and π − → l −ν l γ where l = µ, e. For the emission of electrons from muon decay, the radiative process is subdominant and we consider the tree-level differential decay spectrum in the rest frame of a muon Given muon energy E ′ µ in the lab frame, the energy and momentum of electron in the lab frame is related to their muon rest frame values via the Lorentz boost where γ µ = E ′ µ /m µ , β µ = 1 − 1/γ 2 µ , and θ ( ′ ) is the angle between the rest (lab) frame electron momentum and muon momentum. The decay spectrum is also boosted via a Jacobian The last equality holds as the rest frame spectrum is independent of cos θ. The Jacobian can be evaluated using Eqs. (B2) and (B3). Explicitly.
where β e and γ e are defined accordingly with E ′ e . The lab frame electron spectrum is then obtained by integrating over the angular distribution, The normalized differential decay spectrum in Eq. (32) is therefore and the normalization where z ≡ m µ /m e . In the rest frame, assuming neutrinos are massless with energy E 1,2 and momentum ⃗ p 1,2 , we have This sets kinematically limits on the electron energy in the rest frame After boost we find the cutoff of electron energy at cos θ = ±1, i.e., Imposing the condition Eq. (B10) on the right hand side of Eq. (B2) we can also find the limits of the angular integral in Eq. (B6), Next we consider charged pion decay. Since π − → µ −ν µ dominates electron production, we will only consider this decay channel. The kinematics in two-body decay is rather simplified and in the pion rest frame muon obtains a single energy and the normalized decay spectrum Boost this into the lab frame and integrate over cos θ ′ we find where β π and γ π are similarly defined as before. The limits of muon energy in the lab frame are reached at E ′ µ,min / max = γ π (E µ,CM ∓ β π E 2 µ,CM − m 2 µ ) .
FIG. 21. Secondary electron spectrum from pion and muon decay. The solid, dashed and dotted lines show the decay spectrum from π ± , µ ± and π 0 respectively. The red and pink lines are obtained from CosmoLED (this work) at the primary particle energy E = 5 GeV and 0.2 GeV. The blue lines depict the secondary electron spectra computed with ExoCLASS, which are independent of primary particle energy.
The electron spectrum from π ± decay can be attained directly after integrating over the intermediate muon energy, with df µ /dE ′ e given in Eq. (B7). The secondary electrons from muon and pion decay are compared with ExoCLASS [50] spectra in Fig. 21. For direct comparison we define x ≡ E k,e /E prim , the ratio between the kinetic energy of electron and the energy of primary particles. We do not include e ± from π 0 decay, which is considered to be subdominant. The ExoCLASS spectra computed from PYTHIA v8.219 [48] are independent of energy. We show the spectra at E prim = 5 GeV, and 0.2 GeV. The high energy spectra are close to that of ExoCLASS, but the difference is more pronounced as the primary particle energy is close to their mass.
For secondary photons, ExoCLASS does not consider the contribution from muon and charmed pion decay. The ExoCLASS secondary photon spectrum π 0 decay agrees with that in Hazma at high π 0 energies. where Φ γ,EBL is the extragalactic isotropic photon flux, E i−1 and E i are the photon energies at redshifts z i−1 and z i respectively, V i is the Universe volume at redshift z i , τ is the absorption probability of a photon with energy E i travelling between redshift z i−1 and z i , and ∆z = z i − z i−1 .
The second and third term of Eq. (D1) which describe the change in flux due to Compton scattering and photon injection are determined from Eqs. (44) and (49) respectively. Calculating the change in flux due to Compton scattering in this way requires performing an integral for each energy bin in the discretized spectrum. That is computationally slow so often approximations are used to simplify this step. A more in-depth discussion about Compton scattering can be found in the next subsection.
The first term in Eq. (D1) accounts for the change of flux due to the expansion of the Universe and the absorption of photons. While the instantaneous changes in flux due to these processes are described by Eqs. (41) and (42) separately, it is convenient to combine them into one term that accounts for the total effect.
Evolving the EBL spectrum with the total effect of the Universe expanding between two redshifts also has the advantage of not needing to calculate derivatives as in Eq. (41). This is done by directly taking into account the two effects that the expansion of the Universe has on the photon flux. Firstly, the increasing volume decreases the number density of photons. This is accounted for in the Vi−1 Vi factor that contributes Secondly, the expansion causes photons to lose energy via redshifting so that As discussed in the next subsection, Compton scattering can sometimes be approximated as causing a fractional energy loss rate for all photons which would be treated as an additional term to Eq. (D3). When the fractional energy loss approximation is not used so the only difference between E i and E i−1 in Eq. (D3) comes from adiabatic expansion, the first term of Eq. (D1) can be written explicitly so that The exponent τ in the first term of Eq. (40) comes from integrating the instantaneous change due to absorption as described in Eq. (42). This exponential suppression accounts for the absorption probability over the time step due to photoionization of neutral gas, pair production from atoms and ions, photon-photon scattering, and pair production off the CMB. Depending on the treatment of Compton scattering, it may also be included in the τ calculation. Assuming that the discretized redshift steps are sufficiently small, τ can be calculated using where dτ dz is determined as in Ref. [63]. With these numerical methods, Eq. (D4) can be used to evolve the EBL spectrum and determine the expected observed flux today. The two computational bottlenecks in this method are the integrals required in solving the upscattered photon flux from ICS and the change in photons flux from Compton scattering. As discussed in Sec. III B, accounting for photons from ICS does not improve the constraints set from the isotropic X-ray and gamma ray spectrum so ICS can be safely ignored to improve computational speed. On the other hand, the treatment of Compton scattering can have an impact on the strength of the constraints so understanding which approximations can be used requires further discussion.

Compton Scattering Approximation
The instantaneous rate of change to the EBL flux due to Compton scattering is fully described by Eq. (44). However, when using this method for incorporating the effect of Compton scattering into the discretized evolution of the X-ray background as done in Eq. (D1) there are two potential issues that need to be addressed.
One potential issue is that by assuming the total change in flux due to Compton scattering is equal to dΦγ,comp dEdz ∆z as done in Eq. (D1) there is an implicit assumption that during a redshift step photons either do not scatter or scatter once. It does not allow for multiple Compton scatters of a single photon within a single step. This is valid as long as the redshift steps are sufficiently short. The maximum scattering rate is for low energy photons in the Thomson limit where σ c ≈ σ T . Therefore, the condition that must be true for this treatment of Compton scattering to be valid is σ T n e (z)∆t ≪ 1 (D6) where ∆t is the absolute time of the redshift step. The other issue with this treatment of Compton scattering is that using Eq. (44) to determine the effect of Compton scattering during a redshift step requires computing an integral to determine the change for each energy bin. This can be computationally intensive. Therefore, there are different approximations that can be used depending on the regime of interest and the accuracy needed. They are: • Attenuation -A simple approximation is to ignore the downscattered photons and assume that all photons that Compton scatter are fully absorbed. This would be implemented by treating Compton scattering as an additional component of dτ dz (E, z) in Eq. (D5) where the Compton component is given by dτ dt (E, z) compton = n e (z)σ c (E).
This is generally a conservative and computationally simple approximation to make. This approximation is able to do a good job of estimating how much the flux of high energy photons is attenuated but it breaks down with low energy photons because while they may scatter frequently, they only lose a small fraction of their energy on each scatter. Additionally, if the calculation needs to accurately calculate the shape of the low energy flux this approximation cannot be used. By ignoring the downscattered photons, the predicted flux of low energy photons will be too small.
• Fractional Energy Loss -The opposite limit of attenuation is where all photons scatter however they only lose a small fraction of their energy on each scatter. That is true in the case of photons with E ≪ m e . With the additional assumption that all photons of a given energy lose energy at the same rate which again is valid in the limit of each photon scattering many times, Compton scattering can be included as an additional form of energy loss similar to redshifting. Eq. (D3), which describes how the photon energy changes of a redshift step becomes with dE dz determined as in Ref. [63]. This does make determining the derivative dEi−1 dEi in Eq. (D1) more challenging. Therefore, when using this approximation Compton scattering and redshifting were treated sequentially. The photon spectrum was first changed accounting for redshifting and then the effect of Compton scattering was accounted for. Instead of calculating dEi−1 dEi directly, we integrated the differential flux, dΦγ dE , to determine the total flux Φ γ (E), and then took the derivative with respect to the shifted energy bins E ′ . While Ref. [63] provides an expression for dE dz for all energies, the assumptions underlying this approximation are not valid for high energy photons or when only some photons scatter during a single step. The constraints found using this approximation do match the complete calculation more closely than the attenuation approximation however due to the assumptions breaking down some accuracy is sacrificed in comparison to using Eq. (44).
• The last approximation is to use Eq. (44) to determine the proper Compton scattering effect only for black holes that have fully evaporated before today. The Universe is transparent to Compton scattering for photons originating at z < 100 and if the black holes still exist today, the signal will be dominated by photons produced recently. This is the approximation that was used to produce the final constraints in this work. For black holes evaporated before today we perform the full computationally intensive calculation and for black holes that are still around we use the fractional energy loss approximation.
A comparison of the effect the different Compton scattering approximations have on the constraints on PBH abundance with n = 2 can be seen in Fig. 22. For more massive PBHs that finish evaporating at later times Compton scattering stops being important and all approximations converge. For n > 2 the pattern is similar except the effect of Compton scattering is less and therefore the differences between the various approximations are less important.