# Refining Galactic primordial black hole evaporation constraints

Pedro De la Torre Luque<EMAIL_ADDRESS>Instituto de Física Teórica,
IFT UAM-CSIC, Departamento de Física Teórica, Universidad Autónoma de Madrid,
ES-28049 Madrid, Spain  Jordan Koechler<EMAIL_ADDRESS>Laboratoire
de Physique Théorique et Hautes Energies (LPTHE), UMR 7589 CNRS & Sorbonne
Université, 4 Place Jussieu, F-75252, Paris, France  Shyam Balaji
<EMAIL_ADDRESS>Physics Department, King’s College London, Strand,
London, WC2R 2LS, United Kingdom

(, )

###### Abstract

We revisit the role of primordial black holes (PBHs) as potential dark matter
(DM) candidates, particularly focusing on light asteroid-mass PBHs. These PBHs
are expected to emit particles through Hawking evaporation that can generate
cosmic rays (CRs), eventually producing other secondary radiations through
their propagation in the Milky Way, in addition to prompt emissions. Here, we
perform a comprehensive analysis of CR signals resulting from PBH evaporation,
incorporating the full CR transport to account for reacceleration and
diffusion effects within the Milky Way. In particular, we revisit the
$e^{\pm}$ flux produced by PBHs, using Voyager 1, and study for the first time
the diffuse X-ray emission from the up-scattering of Galactic ambient photons
due to PBH-produced $e^{\pm}$ via the inverse Compton effect using Xmm-Newton
data, as well as the morphological information of the diffuse $511$ keV line
measured by Integral/Spi. In doing so, we provide leading constraints on the
fraction of DM that can be in form of PBHs in a conservative way, whilst also
testing how different assumptions on spin and mass distributions affect our
conclusions.

††preprint: KCL-2024-34

##  I Introduction

The continuous non-detection of definitive, non-gravitational signals from
dark matter (DM) candidates, especially those within the weakly interacting
massive particles category [1], has sparked renewed interest in primordial
black holes (PBHs) [2, 3, 4]. While gravitational lensing significantly
constrains the DM fraction that PBHs can account for at higher masses [2], it
is nearly impossible to detect PBHs with masses much below $10^{20}$ g due to
finite-size source effects [5]. Nevertheless, black holes around this mass
range are expected to emit intense non-thermal radiation through Hawking
evaporation, providing an alternative avenue to constrain their contribution
to DM [3].

The temperature of black holes is inversely proportional to their mass, with
$T\sim 0.1$ keV $(10^{20}\text{ g}/M_{\rm BH})$, meaning the radiation emitted
for masses below $10^{20}$ g is anticipated in the X-ray to $\gamma$-ray
spectrum. Furthermore, if $M_{\rm BH}\lesssim 10^{16}$ g, the evaporation
process also produces positrons in abundance. These positrons, upon
annihilation with ambient electrons, generate detectable $\gamma$-rays,
particularly at an energy of $E_{\gamma}=m_{e}=511$ keV [6]. This finding has
prompted efforts to determine the mass fraction of black holes in the vicinity
of $M_{\rm BH}\sim 10^{16}$ g or potentially greater, by utilizing 511 keV
observations [6, 7, 8, 9]. Notably, the recent study Ref. [10] established
rather cautious constraints based on the 511 keV emissions detected from the
center of our galaxy.

The 511 keV line is not the only signal though, with indirect searches for
charged particle injection (that does not form positronium) from PBH
evaporation like antiprotons, electrons, and positrons having been examined as
well [11, 12]. The major challenge with using these particles to constrain
PBHs is their susceptibility to the Sun’s influence at sub-GeV energies, which
significantly supresses their flux when entering the heliosphere. On the other
hand, low-energy physics is crucial because, the greater the mass of the PBH,
the lower the energy of the emitted particles. Since the Voyager 1 probe has
passed beyond the heliopause, it has detected low-energy electrons and
positrons [13, 14] that may originate from PBH evaporation [15] and is not
affected by solar screening.

Also, non-observation of Hawking radiation signatures from PBHs evaporation in
the keV-MeV energy band probing the inner regions of the Milky Way [16, 17,
18] have been considered using data from various satellites like Xmm-Newton
[19]. Data from this satellite have previously been used to constrain feebly
interacting particles [20, 21] and sub-GeV DM [22, 23, 24] as well.

For a PBH DM explanation of the 511 keV line, previous analyses, such as Ref.
[9] indicated that the black hole number density of the local halo would yield
a median distance that falls approximately within the confines of our Solar
System. Considering the expected PBH velocity in our galactic halo, it is
plausible that we could anticipate that PBHs could travel through our Solar
System relatively frequently. Despite the proximity, detecting Hawking
radiation emitted by a PBH of this mass would pose a considerable challenge.
Conversely, by analyzing the widespread MeV-scale $\gamma$-ray emissions from
the Milky Way halo, we could conclusively verify these types of scenarios
using upcoming observations.

Here, we present new constraints on the PBH fraction from the Galactic diffuse
X-ray emission, using Xmm-Newton data, and leverage the recently reported
longitude profile of the $511$ keV emission line. In addition, we revisit the
PBH constraints from the local interstellar $e^{\pm}$ flux measurements from
state-of-the-art propagation scenarios, using Voyager 1 data. The main goal of
this work consists of performing a more realistic computation of the various
cosmic ray (CR) signals that would be produced by PBH evaporation when
considering the full-fledged CR transport in the Milky Way, and evaluate the
impact of uncertainties present in the propagation and injection setup on our
limits. We utilize a fully numerical approach that does not require
approximations and uses current state-of-the-art astrophysical ingredients.

The paper is organized as follows. In Sec. II we cover the fundamentals of PBH
evaporation, in Sec. III we discuss $e^{\pm}$ propagation from particle
injections by PBHs, in Sec. IV we discuss our main results and finally in Sec.
V we conclude.

##  II PBH evaporation

Black holes are known to evaporate over time and emit particles with masses
below or comparable to the temperature $T$ of the black hole through Hawking
radiation [25]. This temperature is directly related to the mass $M$ of the
black hole and its spin parameter $a\equiv J/M$ with $J$ being its angular
momentum [26] (with $\hbar=c=k_{B}=G=1$):

$T=\frac{1}{2\pi}\left(\frac{r_{+}-M}{r_{+}^{2}+a^{2}}\right)\;,$ (1)

where $r_{+}\equiv M+\sqrt{M^{2}-a^{2}}$ is the Kerr black hole horizon
radius. For Schwarzschild black holes ($a=0$), we retrieve $T=1/(8\pi M)$.
Then, the emission spectrum of primary particles $i$ is given by

$\frac{d^{2}N_{i}}{dtdE_{i}}=\frac{1}{2\pi}\sum_{\text{d.o.f.}}\frac{\Gamma_{i}(E_{i},M,a^{\star})}{e^{E^{\prime}_{i}/T}\pm
1}\;,$ (2)

where $a^{\star}\equiv a/M<1$ (if $a^{*}=1$, then $T=0$, which is forbidden by
thermodynamics) is the reduced spin parameter, $E^{\prime}_{i}\equiv
E_{i}-m\Omega$ is the energy of the emitted particle where $\Omega\equiv
a^{\star}/(2r_{+})$ is the black hole horizon rotation velocity and
$m=\\{-l,...,l\\}$ the projection on the black hole axis of the particle
angular momentum $l$. The $\pm$ signs depend on the spin of the particles
radiated: $+$ for fermions and $-$ for bosons. Finally, $\Gamma_{i}$ are the
so-called ‘greybody factors’ which encode the deviation from black-body
physics, since emitted particles have to escape the gravitational well of the
black hole. The sum is performed over the degrees of freedom (d.o.f.) of the
emitted particles (spin, color, helicity and charge multiplicities). In order
to compute the spectra of primary particles, we use the numerical code
`BlackHawk v2.2` [27, 28].

Then, evaporated particles can hadronise, decay or emit soft radiations.
`BlackHawk` also has the possibility of dealing with such processes by
including tables from the particle physics codes `PYTHIA` [29], `Herwig` [30]
and `Hazma` [31]. However, their domains of validity differ, `PYTHIA` and
`Herwig` handle processes with particle energies above 10 GeV very well,
whereas `Hazma` excels below the QCD scale ($\sim$ 250 MeV). Since we are
interested in the production of sub-GeV $e^{\pm}$, we only use `Hazma` to
treat secondary processes, and limit ourselves to its domain of validity,
which corresponds to a black hole mass range of $M\gtrsim 10^{14.5}$ g. An
upper limit on the black hole mass can also be defined when the evaporation
into $e^{\pm}$ is not possible anymore ($T\ll m_{e}$) corresponding to
$M\gtrsim 10^{17.5}$ g. In this black hole mass range, $e^{\pm}$,
$\nu_{e,\mu,\tau}$ and $\gamma$ are emitted through evaporation, and to some
extent (for lower black hole masses) $\mu^{\pm}$ and $\pi^{0,\pm}$. `Hazma`
handles the decay of $\mu^{\pm}$ and $\pi^{0,\pm}$, as well as final state
radiation from all charged particles. From now on,
$\frac{d^{2}N_{i}}{dtdE_{i}}$ will describe the emission spectra of secondary
particles.

Although it is believed that PBHs cannot acquire a substantial spin from their
production process [32] unless formed in the matter-dominated universe [33],
it has been argued that they can do so by repeatedly merging with other black
holes. Moreover, black holes can also acquire a spin due to the accretion of
gas surrounding them. The maximum spin value a black hole can obtain through
this process is $a_{\textrm{lim}}^{*}\approx 0.998$ [34], known as the Thorne
limit, and can slightly vary depending on the considered accretion model [35].
A similar limit can be derived for black hole mergers [36]. Nevertheless, PBHs
formed during the matter-dominated universe could evade these limits,
providing a smoking gun signature of their existence. In our study, we decide
to remain agnostic on the PBH production process and consider two extreme
benchmarks to quantify the impact of the choice of the spin distribution on
our results. We therefore examine the case where all PBHs are Schwarzschild
black holes, and another one where they are all near-extremal Kerr black holes
with a spin of $a^{\star}=0.9999$.

In Figs. 9 and 9 we show the spectra of secondary $e^{\pm}$ and $\gamma$
respectively for a range of black hole masses, and for spins of $a^{\star}=0$
and $0.9999$. For increasing values of $a^{\star}$ the black hole evaporates
more particles with higher energies, due to the transfer of angular momentum
from the black hole to the emitted particles.

All of the discussion above only takes into account the physics from a single
black hole. We then have to take into account the energy density of PBHs in
our Galaxy. We investigate the possibility of PBHs constituting a fraction
$f_{\textrm{PBH}}$ of the total amount of DM in the Universe. Therefore the
number of $e^{\pm}$ injected by PBHs evaporation at the position vector
$\vec{x}$ in our Galaxy per unit of time, energy and volume is written

$Q_{e}(E_{e},\vec{x})=f_{\textrm{PBH}}\rho_{\textrm{DM}}(\vec{x})\int_{M_{\textrm{min}}}^{\infty}\frac{dM}{M}\frac{dN_{\textrm{PBH}}}{dM}\frac{d^{2}N_{e}}{dtdE_{e}}\;,$
(3)

where $\rho_{\textrm{DM}}$ is the DM energy density profile and
$\frac{dN_{\textrm{PBH}}}{dM}$ is the mass distribution of PBHs normalized to
1, since $f_{\textrm{PBH}}\,\rho_{\textrm{DM}}$ already represents the spatial
distribution of the PBH energy density in the Galaxy. $M_{\textrm{min}}\approx
7.5\times 10^{14}$ g corresponds to the minimal mass of PBHs today. As shown
in Fig. 1, PBHs formed in the early Universe with a mass below
$M_{\textrm{min}}$ should have all evaporated by now, whereas PBHs with an
initial mass of $10^{15}$ g have experienced their mass decreasing by
$\mathcal{O}(20\%)$. Thus, we follow the approximation where all PBHs with
$M_{\textrm{PBH}}\leq M_{\textrm{min}}$ do not exist today, and the remaining
ones have the same mass distribution as in their formation in the early
Universe.

Figure 1: Evolution of the mass $M$ of Schwarzschild BHs for different initial
masses $M_{0}$ at $t=0$. The $x$-axis represents the time in terms of
fractions of the age of the Universe and $y$-axis the BHs mass in terms of
fractions of its initial mass.

We consider the following PBHs mass distribution

$\frac{dN_{\textrm{PBH}}}{dM}=\frac{1}{\sqrt{2\pi}\sigma
M}\exp\left(-\frac{\log^{2}(M/M_{\textrm{PBH}})}{2\sigma^{2}}\right)\;,$ (4)

where $M_{\textrm{PBH}}$ is the peak PBH mass, and $\sigma$ is the standard
deviation of the distribution. This mass function is relevant when assuming
the formation of PBHs from a large enhancement in the inflationary power
spectrum [37, 38, 39, 40]. We explore values of $\sigma$ varying from 0 to 2,
the case $\sigma\rightarrow 0$ corresponding to a monochromatic mass
distribution ($\frac{dN_{\textrm{PBH}}}{dM}=\delta(M-M_{\textrm{PBH}})$).

##  III Electron-positron propagation and diffuse emissions

###  III.1 Electron-positron propagation

Light asteroid-mass PBHs evaporating in the Galaxy produce a continuous
injection of electrons and positrons that eventually leads to a diffuse and
steady-state flux peaking at energies above the MeV scale. Once injected, as
described in Eq. 2, these particles propagate and interact with the Galactic
gas, magnetic fields and ambient light. For particles propagating with
energies below a few tens of MeV the main process is energy losses, dominantly
from ionisation of neutral gas, although bremsstrahlung and Coulomb
interactions with the ambient plasma are relevant too [22, 41]. Diffuse
reacceleration [42, 43], that is the product of energy exchange between
relativistic charged particles and the perturbations of the plasma, also
becomes very relevant at such energies, since the low mass and energy of the
$e^{\pm}$ produced from these PBHs can be easily boosted to much higher
energies, even reaching GeV energies [22]. All of the mentioned processes are
encoded in the diffusion-advection-convection-loss equation [44, 45]

$-\nabla\cdot\left(D\vec{\nabla}f_{e}+\vec{v}_{c}f_{e}\right)-\frac{\partial}{\partial
p_{e}}\left[\dot{p}_{e}f_{e}-p_{e}^{2}D_{pp}\frac{\partial}{\partial
p_{e}}\left(\frac{f_{e}}{p_{e}^{2}}\right)-\frac{p_{e}}{3}\left(\vec{\nabla}\cdot\vec{v}_{c}f_{e}\right)\right]=Q_{e}\;,$
(5)

which has to be solved for $f_{e}\equiv\frac{dn_{e}}{dp_{e}}$, the density of
$e^{\pm}$ per unit momentum at a given position. This equation includes: i)
spatial diffusion with coefficient $D$, ii) momentum losses $\dot{p}_{i}$ due
to interactions with the galactic environment, iii) momentum diffusion (or
reacceleration) with diffusion coefficient $D_{pp}$, iv) convection due to the
galactic wind velocity $\vec{v}_{c}$. We solve Eq. (5) using the numerical
code `DRAGON2` [45, 46]. As a benchmark, we use a Navarro-Frenk-White (NFW)
[47] DM distribution. We have computed tables with the $e^{\pm}$ injection
spectra from PBH evaporation ($\frac{d^{2}N_{e}}{dtdE_{e}}$) using
`BlackHawk`. These are an input in the `DRAGON2` code and are used analogously
to the case of decaying DM. These tables are available upon request from the
authors and will be released along with a new public version of `DRAGON2`.

In our benchmark setup, the diffusion parameters are the same as prescribed in
Ref. [22], where we refer the reader for full details. These parameters are
obtained from detailed fits of existing CR data at Earth. However, given that
systematic uncertainties are difficult to asses and different CR analyses can
find slightly different results [48, 42, 43], we consider other extreme
propagation scenarios that allow us to evaluate the impact of these
uncertainties in the predicted $e^{\pm}$ spectra and the bounds on the
fraction of DM constituted by PBHs.

Figure 2: Local $e^{\pm}$ spectrum ($e^{+}$ \+ $e^{-}$) generated from the
evaporation of light PBHs under different assumptions. Top-left: Comparison of
the spectrum for Schwarzschild PBH of different masses, assuming their mass
distribution to be monochromatic. Top-right: Comparison of the expected local
$e^{\pm}$ spectrum from Schwarzschild ($a^{*}=0$) and near-extremal Kerr
($a^{*}=0.9999$) PBHs, assuming they are monochromatic in mass. Bottom-left:
Comparison of the predicted local $e^{\pm}$ spectrum for Schwarzschild PBHs
assuming a log-normal mass distribution with different values of the standard
deviation $\sigma$. Bottom-right: Comparison of the local $e^{\pm}$ spectrum
for monochromatic Schwarzschild PBHs, for different levels of reacceleration
(i.e. values of $V_{A}$).

We first repeat our calculations for ‘realistic’ variations of the propagation
parameters found in our analysis, which consist of taking the values that
maximize the difference in flux from the benchmark case at $3\sigma$. Similar
to the case of DM decay, the parameters with a greater effect on the diffuse
spectra produced from PBHs are the Alfvén velocity $V_{A}$ parameter that
controls the level of diffuse reacceleration [21, 22] and the height $H$ of
the halo, which dictates the volume where CRs are confined and where PBHs
produce particles that can reach us. In this way, to obtain a realistic
uncertainty band in our predictions, we use a conservative setup where $H=4$
kpc and $V_{A}=7$ km/s, which produces a lower (and therefore more
conservative) flux of $e^{\pm}$. In turn, the more aggressive setup is meant
to increase the flux of $e^{\pm}$ from PBHs, and uses values of $H=12$ kpc and
$V_{A}=20$ km/s. As a point of reference, we recall that our benchmark values
are $H=8$ kpc and $V_{A}=13.4$ km/s. We tested an even more ‘general’ and
extreme variation of propagation setup, that ensures that the flux of
particles must be between the two extremes: The ‘optimistic’ case, where
$H=16$ kpc and $V_{A}=40$ km/s, is much higher than the typical values111We
consider that $V_{A}=40$ km/s is the maximum realistic value for $V_{A}$,
since it already implies that most of the injected energy of CRs are coming
from the perturbations of the interstellar plasma and not from supernova
remnants, which will break the standard paradigm of CR propagation, see Ref.
[49, 50]. Then, the most ‘pessimistic’ case will be that with no
reacceleration ($V_{A}=0$ km/s) and $H=3$ kpc. These two cases are unlikely,
given the fact that propagating CRs implies energy exchange with plasma waves
and therefore non-zero reacceleration, and values below $H=3$ kpc seem to be
strongly disfavoured from CR analyses and other existing constraints [51, 52,
53, 54]. As an example, we show in the bottom-right panel of Fig. 2 the
dramatic effect of reacceleration in the $e^{\pm}$ spectrum at Earth from a
$M_{\textrm{PBH}}=10^{16}$ g PBH, for different values of $V_{A}$, comparing
our benchmark scenario with the aforementioned optimistic and pessimistic
ones.

Figure 3: Comparison of the predicted DM-induced X-ray emission with diffuse
X-ray data from Xmm-Newton in a region close to the Galactic Center. We show
the prediction is shown for different values of $M_{\textrm{PBH}}$ when the
PBH mass distribution is monochromatic (upper left panel), of $a^{\star}$
(upper right panel) and $\sigma$ when the distribution is log-normal (lower
panel).

We will obtain constraints on the fraction of DM that can be in the form of
PBHs using Voyager 1 [55, 13] measurements of the local flux. A comparison of
the total electron-positron (i.e. $e^{+}$ \+ $e^{-}$) flux measured by Voyager
1 and the predicted local electron spectrum for monochromatic PBHs of
different masses, for our benchmark propagation setup and NFW DM profile, is
shown in the top-left panel of Fig. 2. In addition, we also illustrate the
spectra predicted assuming a log-normal PBH mass distribution (see Eq. 4) with
$\sigma=1$ (green line) and $\sigma=2$ (red line) in the bottom-left panel.
This allows one to see how the $\sigma$ parameter affects our predictions,
given that physically one must expect a non-zero $\sigma$. It can be seen that
the higher $\sigma$ is, the higher is the expected flux and the higher is the
energy reached by the electrons. The main reason for this is that the
contribution from lower-mass PBHs is very important and dominates the spectra
of these particles.

Then, in the top-right panel we compare the spectra produced from PBHs with
different values of the spin parameter $a^{*}$. In particular, we show the
cases of Schwarzschild ($a^{*}=0$) and near-extremal Kerr ($a^{*}=0.9999$)
PBHs, as well as for $a^{*}=0.99$. As one can see from the figure, spin of the
PBH always leads to a higher flux produced by PBH evaporation although not
changing its spectral shape appreciably, in agreement with what was found in
Refs. [56, 28].

###  III.2 Diffuse X-ray emission

During their propagation, the population of electron-positron pairs injected
in the Galaxy produce different secondary radiations that can be used to track
their distribution and density. Especially, their interaction with the
Galactic magnetic field will generate a diffuse synchrotron emission that can
be observed at kHz-MHz. On top of that, they will interact via bremsstrahlung
with the ionized gas in the interstellar medium (ISM), leading to $\gamma$-ray
radiation at the MeV scale. Here, we focus on the X-ray diffuse emission
produced from the inverse Compton interaction of this population with the
Galactic radiation fields (mainly the cosmic microwave background, optical and
UV light from stars and infrared from the scattering of the latter on dust),
because of the high constraining power of X-ray measurements, as mentioned
above.

Figure 4: Comparison of the expected longitude profile of the $511$ keV
emission from PBH evaporation with Integral/Spi data. The left panel shows a
comparison of the predicted profile for different DM distributions, for the
case of Schwarzschild PBH monochromatically distributed with
$M_{\text{PBH}}=10^{16}$ g, while the right panel shows a comparison of the
signals expected from different PBH masses, assuming again Schwarzschild and
monochromatically distributed PBHs.

To calculate the diffuse X-ray emission generated from the diffuse (steady-
state) distribution of $e^{\pm}$ in the Galaxy that we obtain with `DRAGON2`,
we employ the `Hermes` code [57], a software designed to compute full maps of
nonthermal radiative processes such as radio, X-ray and $\gamma$-rays, as well
as neutrino emissions. The total X-ray flux also includes photons directly
emitted during PBH evaporation, as well as final state radiations produced by
evaporated $e^{\pm}$, $\mu^{\pm}$ and $\pi^{\pm}$. It turns out that this
component is sub-dominant compared to the X-ray flux emitted during the
transport of evaporated $e^{\pm}$. We compute $2\sigma$ bounds from the
diffuse Galactic X-ray emission observed by Xmm-Newton [58, 59] in the $2.5-8$
keV band, as done in Ref. [22, 23], where we refer the reader for more
details. In Fig. 3, we compare the X-ray diffuse emission expected from PBH
evaporation. In the top-left panel, we show the case of monochromatic PBHs
with masses of $10^{15}$ g, $10^{16}$ g and $10^{17}$ g, for our benchmark
propagation setup and NFW DM profile. In the top-right panel, we compare the
emission expected from a Schwarzschild PBH ($a^{*}=0$) and Kerr PBHs with
$a^{*}=0.99$ and the extreme case of $a^{*}=0.9999$, all for
$M_{\textrm{PBH}}=10^{16}$ g. Note that the fraction of PBHs comprising DM is
different for every case, as indicated in the legend. In the bottom panel, we
show results for Schwarzschild PBHs distributed log-normally with a mean mass
of $10^{16}$ g and different values of $\sigma$, ranging from $\sigma=0$
(monochromatic case) to the wider $\sigma=2$. The conclusions for the impact
of these parameters in the X-ray Galactic diffuse emission are similar to
those found for the local flux of Fig. 2. We finally remark that the
associated bremsstrahlung emission is negligible at keV energies.

###  III.3 511 keV line profile

As the injected positrons are propagating in the Galaxy, they lose their
energy until they eventually reach thermal energies of the medium that they
are travelling through. After a typical time scale of $0.1-10$ Myr [60, 61],
thermal positrons will form a positronium state with ambient electrons and
decay into a pair of $511$ keV photons with $25\%$ probability (through the
para-positronium state), that lead to a bright line emission. Assuming that
the thermal positrons follow the distribution of the steady-state diffuse
positrons, we have calculated the intensity and distribution of the $511$ keV
line from low mass DM annihilating and decaying in Ref. [62] and this case
leads to very similar signals to those expected from PBH evaporation.

In this work, we obtain constraints from the longitude profile of the $511$
keV line emission following the procedure described in Ref. [63], where the
profile of the line is directly proportional to the distribution of propagated
(steady-state) positrons. We have tested that our results are compatible with
previous evaluations applied to other exotic sources of positrons [64, 65, 20,
21, 66].

Figure 5: Comparison of the 95% confidence limits on $f_{\textrm{PBH}}$
derived in this work with other existing ones. The color of the lines
represent the different probes used to set the constraints: green for the
$e^{\pm}$ measurements from Voyager 1, blue for X-ray diffuse observations,
red for the 511 keV excess reported by Integral and grey for the 21-cm line
measurements from Edges. The two different line styles correspond to either
the bounds derived either in this work (solid) or in the literature [15, 67,
68, 69] (dashed) reported in ``PBHbounds— [70].

In the left panel of Fig. 4 we show the predicted longitude profile of the
$511$ keV emission for PBHs of masses between $10^{15}$ and $10^{17}$ g
compared to Integral/Spi data [71], assuming a NFW profile and with the PBH
fraction of DM $f_{\textrm{PBH}}$ specified in the legend for each case. It
can be seen that the most constraining data points are those obtained at high
longitudes. Given that these points are also those expected to be more
affected by systematic uncertainties (mainly background noise and the need of
templates to extract measurements) and the limited statistics of the
measurements, the bounds that we derive are conservatively calculated by
multiplying by a factor of $2$, as a proxy for the effect of systematic
uncertainties, as done in Refs. [21, 63]. We show in the right panel of Fig. 4
a comparison of the predicted line profile with the NFW DM distribution with
other popular DM profiles, namely a Moore profile ($\gamma=1.5$) [72], a
contracted NFW profile similar to the one fitting the Galactic Center excess
($\gamma=1.25$) [73], for a monochromatic $10^{16}$ g mass PBH. A Burkert
[74], or other cored DM distributions [75], will simply lead to a flatter
profile. As one can see, the predicted profiles are very similar at high
longitudes and only change significantly around the central longitudes.
Therefore, the uncertainties in the derived constraints from the DM
distribution are very small. Similarly, a different spin or adoption of
$\sigma\neq 0$ has no significant consequences on the shape of the profile and
essentially changes only the normalization of the signal.

##  IV Results and comparison with other work

In this section, we discuss the limits on PBHs we derived in this work, for
our benchmark scenario. In addition to displaying the limits using the three
probes ($e^{\pm}$, 511 keV line and diffuse $X$-rays), we also show the impact
on these limits when assuming different PBH mass distributions, spin
distributions and propagation models. Here, we set $2\sigma$ bounds on
$f_{\textrm{PBH}}$ and $M_{\textrm{PBH}}$ by applying the criterion

$\sum_{i}\left(\frac{\textrm{Max}\left[\phi_{\textrm{PBH},i}(M_{\textrm{PBH}})-\phi_{i},0\right]}{\sigma_{i}}\right)^{2}=4\,,$
(6)

where $i$ denotes the data point, $\phi_{\textrm{PBH}}$ is the predicted flux
induced by PBH evaporation, $\phi_{i}$ is the measured flux and $\sigma_{i}$
the associated standard deviation of the measurements.

In Fig. 5 we show our benchmark limits on Schwarzschild PBHs, assuming a
monochromatic mass distribution and a NFW DM profile, and compare them to
existing ones. The solid lines represent the bounds derived in this work,
while the dot-dashed lines represent some of the most stringent limits on
$f_{\textrm{PBH}}$ reported in `PBHbounds` [70] across the $10^{15}-5\times
10^{17}$ g PBH mass range.

Figure 6: Limits we derive using Voyager 1 (upper left panel), Xmm-Newton
(upper right panel) and the 511 keV line reported by Integral (lower panel)
for mass distributions with $\sigma=0$ (monochromatic, solid black line), 1
(red line) and 2 (green line), as well as for the extreme assumption where all
PBHs have a spin of $a^{\star}=0.9999$ (dashed line).

We show the Voyager 1 limits in green on Fig. 5, where the dashed line
corresponds to the limit reported in Ref. [15] without background subtraction.
The authors used a propagation model with strong reacceleration named ‘model
B’. Our Voyager 1 limit is comparable to the existing one for
$M_{\textrm{PBH}}\lesssim 10^{16}$ g and gets more stringent for higher PBH
masses. The reason is likely due to the differences in how reacceleration is
implemented in the `DRAGON2` code with respect to the the semi-analytical code
`USINE` [76]222https://dmaurin.gitlab.io/USINE/ used in their work, where
reacceleration only takes place in a thin disk, instead of adopting uniform
reacceleration across the whole Galaxy, that is important given that CR
particles spend most of their time in the Galactic halo while propagating. In
addition, to model energy losses, which are key for MeV particles, `USINE`
needs to make use of the pinching method [77].

Figure 7: Uncertainties in the limits we derive using Voyager 1 (upper left
panel), Xmm-Newton (upper right panel) and the 511 keV line reported by
Integral (lower panel). The solid blue lines correspond to the limits using
our fiducial propagation model. The blue bands show how our limits are
impacted when varying the Alfvén speed $v_{A}$ and the halo height $H$ within
their $3\sigma$ uncertainty. The grey bands correspond to a more conservative
scenario where we vary $v_{A}$ between 0 and 40 km/s and $H$ between 3 and 16
kpc.

The limits from diffuse X-ray emissions are shown in blue on Fig. 5, where the
dashed line is the limit set in Ref. [68]. The authors have computed the flux
of (prompt) X-ray emissions from the evaporation of extragalactic PBHs and
compared it to the isotropic cosmic X-ray background measurements, without
considering the secondary inverse Compton emission, to set a limit on
$f_{\textrm{PBH}}$. Remarkably, the low energy part of the X-ray measurements
are those most constraining. Therefore, X-ray diffuse measurements at lower
energies are expected to improve these limits significantly. However, X-ray
emission starts to be severely absorbed by interstellar gas, which can make it
more difficult to improve these constraints using lower energy data. Current
work in progress with eRosita Galactic diffuse data indicate that our limit
can improve by up to an order of magnitude. Furthermore, we note that the most
constraining X-ray data corresponds to the inner regions (see Refs. [22, 23]).

Our $511$ keV bound, which we weaken by a factor of two to account for
systematic uncertainties in the data (as mentioned when discussing the
calculation of the $511$ keV line in Sec. III), is shown in red in Fig. 5,
where we compare with the limit reported in Ref. [67] (red dashed line). They
used the rate of positron injection needed to explain the total $511$ keV flux
from Integral in the bulge. Given that the high-longitude measurements of the
$511$ keV line emission are the most constraining measurements, the use of
longitudinal profile lead to more stringent results compared to using the
bulge emission [62]. In addition, the authors include only the emission from a
NFW DM profile within the inner $3$ kpc from the Galactic Center and did not
model positron propagation. Therefore our 511 keV bound appears to be more
stringent than the one of Ref. [67]. In a recent paper [63], we show that
using in-flight positron annihilation emission can improve the limits on
feebly interacting particles (whose $e^{\pm}$ emission is also concentrated at
tens of MeV and follow a similar spatial morphology) with respect to the $511$
keV constraints, given that measurements of the diffuse $\gamma$-ray emission
above a few MeV have a reduced systematic uncertainty and more reliable
background models can be used.

Finally, we report the bound derived by requiring that the amount of heating
of the intergalactic medium from PBH evaporation is constrained by 21-cm
observations by the Edges experiment [69], which we show as a dotted gray line
in Fig. 5. All in all, our limit from the longitude profile of the $511$ keV
line is competitive with the Edges limit below $M_{\textrm{PBH}}\simeq
10^{16}$ g and becomes the most stringent limit to date for PBH masses between
$10^{16}$ and $2\times 10^{17}$ g, for this propagation setup.

In Fig. 5, we assumed PBHs to be Schwarzschild BHs with a monochromatic mass
distribution. The main reason is because non-rotating and monochromatic PBHs
represent the most conservative case. However, if their mass distribution were
instead log-normal, as in Eq. (4), there would be a low-mass PBH population
that contributes to most of the flux of evaporated $e^{\pm}$ and photons,
leading to a strengthening of the limits. Actually, for increasing values of
the standard deviation $\sigma$ of the distribution, the low-mass population
increases and therefore the limits become more and more stringent.
Alternatively, if PBHs were Kerr BHs, they would produce more particles at
high energies, ending up with a strengthening of the limits as well. Fig. 6
illustrates the impact of the choice of mass and spin distributions on the
limits on $f_{\textrm{PBH}}$. We also notice that in the case of a cuspier DM
distribution than the NFW (i.e. slope of $\gamma>1$) the Xmm-Newton
constraints could beat those from Edges or the $511$ keV line. In turn, for a
cored one, like a Einasto [78, 79] the limits will be slightly weaker, while
for a Burkert profile [80], the limits become much weaker.

Finally, in Fig. 7 we report the uncertainties on the limits we derive in this
work, showing the impact of the choice of the propagation model, by using the
propagation scenarios explained above. The blue bands (labeled ‘realistic’)
correspond to the variation of $V_{A}$ and $H$ up to their 3$\sigma$
uncertainties. For the lower side of the band we use $V_{A}=20$ km/s and
$H=12$ kpc, while for the upper side we use $V_{A}=7$ km/s and $H=4$ kpc. Then
the gray bands in Fig. 7 (labeled ‘general’) represent more conservative
uncertainties, where for the lower side we adopt $V_{A}=40$ km/s and $H=16$
kpc, and for the upper side $V_{A}=0$ km/s and $H=3$ kpc. In general, these
variations may affect our limits by up to an order of magnitude or more. In
the case of the $511$ keV signals, we observe that the uncertainty bands are
in general smaller. This is due to the fact that the morphology of the
predicted $511$ keV line does not change significantly for different values of
reacceleration (at high longitudes, where the main constraints come from,
reacceleration does not appreciably change the emission, which is spatially
very flat in any case). In the case of Voyager 1 and X-ray constraints, we
observe that the uncertainty band typically broadens at higher PBH masses. The
reason is that, high mass PBHs inject lower energy $e^{\pm}$, which
reacceleration affects much more. For example, in the no reacceleration case
the emission from PBHs of mass higher than a few times $10^{16}$ g lies below
the Voyager 1 data points (therefore, no constraint can be set). Additionally,
we note that uncertainties from the choice of DM distribution will have very
little effect in the constraints from Voyager 1 and the $511$ keV line, while
the limits from Xmm-Newton can be significantly affected, given that the most
constraining X-ray data is that coming from the inner regions of the Galaxy,
which is where our predictions are more affected by uncertainties in the DM
distribution.

##  V Conclusion

In this work, we have conducted a thorough analysis of signals emanating from
PBH evaporation thereby refining constraints on their role as DM candidates.
Our study leverages observations from the Integral/Spi, Voyager 1, and Xmm-
Newton, integrating a comprehensive CR transport model that encapsulates
reacceleration and diffusion within the Milky Way. By numerically solving the
diffusion equations and employing current CR propagation frameworks, we have
honed the limits on PBHs as DM, particularly for those with masses around
$10^{16}$ g.

Our findings indicate that the limits derived from the $511$ keV line,
$e^{\pm}$, and diffuse X-rays are significantly impacted by the assumptions
regarding PBH mass, spin distributions and propagation models. They complement
each other and significantly probe the parameter space available for PBHs as
DM. The most compelling result comes from the $511$ keV line emission at the
Galactic disk (specifically, the high-longitude measurements of the
longitudinal profile of the signal), where our analysis, assuming a NFW DM
profile, yields the most stringent limits to date for PBH masses between
$10^{16}$ and $2\times 10^{17}$ g. This bound is further corroborated by the
heating constraints of the intergalactic medium from PBH evaporation as
observed by the Edges experiment. We additionally remark that our limits are
conservatively derived without including backgrounds, that are expected, in
all the studied cases, to be dominant. We note, however, that the X-ray
constraints produced from inverse Compton emission, are stronger than those
from the $511$ keV line for optimistic diffusion parameters, as shown in Fig.
7, and also in the case of a cuspier DM distribution than the benchmark NFW.

Uncertainties in our limits, depicted in our analysis, underscore the
sensitivity of our results to the choice of propagation model parameters, such
as the Alfvén speed and halo height. Moreover, considering alternative mass
distributions, like a log-normal distribution, or the inclusion of Kerr BHs,
would lead to even more stringent limits due to the increased flux of
evaporated particles and photons. We note that the low mass part of the
asteroid-mass gap is currently well probed and PBHs can constitute a
significant fraction of the DM only in the gap between $10^{18}-10^{21}$ g in
the monochromatic mass case and $\sim 10^{19}-10^{21}$ g in the case where
more realistic log-normal distribution is adopted for their mass distribution.

In conclusion, our comprehensive approach to analysing CR signals from PBH
evaporation has not only refined current astrophysical constraints on PBHs as
DM candidates but also highlighted the critical influence of propagation
models and PBH distributions on these limits. Our work paves the way for
future studies to further explore the intriguing possibility of PBHs
constituting the elusive DM in our universe and also studying late-forming
evaporating PBHs.

###### Acknowledgements.

We acknowledge useful discussions with Alexandre Arbey, Joe Silk and Marco
Cirelli. SB is supported by the STFC under grant ST/X000753/1. PDL is
supported by the Juan de la Cierva JDC2022-048916-I grant, funded by
MCIU/AEI/10.13039/501100011033 European Union “NextGenerationEU”/PRTR. The
work of PDL is also supported by the grants PID2021-125331NB-I00 and
CEX2020-001007-S, both funded by MCIN/AEI/10.13039/501100011033 and by “ERDF A
way of making Europe”. PDL also acknowledges the MultiDark Network, ref.
RED2022-134411-T. This project used computing resources from the Swedish
National Infrastructure for Computing (SNIC) under project Nos. 2021/3-42,
2021/6-326, 2021-1-24 and 2022/3-27 partially funded by the Swedish Research
Council through grant no. 2018-05973.

##  Appendix A Evaporation spectra from BlackHawk

Figure 8: Total spectra of $e^{\pm}$ from the evaporation of a single BH with
spin of $a^{\star}=0$ (left panel) and $a^{\star}=0.9999$ (right panel), for
the following BH masses: $M=10^{14.5}$ g (red), $10^{15}$ g (orange),
$10^{15.5}$ g (yellow), $10^{16}$ g (lime), $10^{16.5}$ g (cyan), $10^{17}$ g
(blue), $10^{17.5}$ g (purple) and $10^{18}$ g (black). Output of
``BlackHawk+Hazma—.

Figure 9: Spectra of secondary $\gamma$ from the evaporation of a single BH
with spin of $a^{\star}=0$ (left panel) and $a^{\star}=0.9999$ (right panel),
for the following BH masses: $M=10^{14.5}$ g (red), $10^{15}$ g (orange),
$10^{15.5}$ g (yellow), $10^{16}$ g (lime), $10^{16.5}$ g (cyan), $10^{17}$ g
(blue), $10^{17.5}$ g (purple) and $10^{18}$ g (black). Output of
``BlackHawk+Hazma—. Same color scheme as in Fig. 9.

In this appendix we discuss the spectra of secondary $e^{\pm}$ and $\gamma$
from the evaporation of a single BH. These are obtained using `BlackHawk`,
which deals with the evaporation of at most $\gamma$, $\nu_{e,\mu,\tau}$,
$e^{\pm}$, $\mu^{\pm}$ and $\pi^{0,\pm}$ in the PBH mass range we consider in
our study. For $M_{\textrm{PBH}}=M_{\textrm{min}}$ all of these particles are
produced, whereas for $M_{\textrm{PBH}}=10^{18}$ g only $\gamma$,
$\nu_{e,\mu,\tau}$, $e^{\pm}$ are evaporated, recalling Eq. 1. Then,
$\mu^{\pm}$ and $\pi^{0,\pm}$ decay after being evaporated, in turn producing
more $\gamma$, $\nu_{e,\mu,\tau}$, $e^{\pm}$. Moreover, all the charged
particles we mentioned can emit final state radiations, increasing the count
of $\gamma$ that originates from PBH evaporation as well. To compute the
spectra of particles coming from these processes, we use the version of
`Hazma` integrated in `BlackHawk`. In Figs. 9 and 9 we show the total spectra
of $e^{\pm}$ and $\gamma$ from BH evaporation for different BH masses and for
$a^{\star}=0$ (0.9999) in the left (right) panel. The spectrum of $e^{\pm}$
solely come from their emission from the PBH for $M_{\textrm{PBH}}\gtrsim
10^{14}$ g, however for lower masses, $\mu^{\pm}$ and $\pi^{\pm}$ start to be
produced and their decay into $e^{\pm}$ contribute to the low-energy bump in
the $e^{\pm}$ spectrum at $M_{\textrm{PBH}}=10^{14.5}$ g. Also, in the
$\gamma$ spectra and for masses where $e^{\pm}$ start to be produced
efficiently, the low-energy ramp corresponds to final state radiations, for
which $dN_{\textrm{FSR}}/dE\propto 1/E$. Finally, we can witness that Kerr BHs
emit more particles at higher energies than Schwarzschild ones.

## References

  * Workman  _et al._ [2022] R. L. Workman _et al._ (Particle Data Group), PTEP 2022, 083C01 (2022).
  * Carr  _et al._ [2021] B. Carr, K. Kohri, Y. Sendouda, and J. Yokoyama, Rept. Prog. Phys. 84, 116902 (2021), arXiv:2002.12778 [astro-ph.CO] .
  * Auffinger [2023] J. Auffinger, Prog. Part. Nucl. Phys. 131, 104040 (2023), arXiv:2206.02672 [astro-ph.CO] .
  * Escrivà  _et al._ [2022] A. Escrivà, F. Kuhnel, and Y. Tada, (2022), arXiv:2211.05767 [astro-ph.CO] .
  * Smyth  _et al._ [2020] N. Smyth, S. Profumo, S. English, T. Jeltema, K. McKinnon, and P. Guhathakurta, Phys. Rev. D 101, 063005 (2020), arXiv:1910.01285 [astro-ph.CO] .
  * MacGibbon and Carr [1991] J. H. MacGibbon and B. J. Carr, Astrophys. J. 371, 447 (1991).
  * Bambi  _et al._ [2008] C. Bambi, A. D. Dolgov, and A. A. Petrov, Phys. Lett. B 670, 174 (2008), [Erratum: Phys.Lett.B 681, 504–504 (2009)], arXiv:0801.2786 [astro-ph] .
  * Carr  _et al._ [2010] B. J. Carr, K. Kohri, Y. Sendouda, and J. Yokoyama, Phys. Rev. D 81, 104019 (2010), arXiv:0912.5297 [astro-ph.CO] .
  * Keith and Hooper [2021] C. Keith and D. Hooper, Phys. Rev. D 104, 063033 (2021), arXiv:2103.08611 [astro-ph.CO] .
  * DeRocco and Graham [2019] W. DeRocco and P. W. Graham, Phys. Rev. Lett. 123, 251102 (2019), arXiv:1906.07740 [astro-ph.CO] .
  * Barrau [2000] A. Barrau, Astropart. Phys. 12, 269 (2000), arXiv:astro-ph/9907347 .
  * Carr  _et al._ [2016] B. J. Carr, K. Kohri, Y. Sendouda, and J. Yokoyama, Phys. Rev. D 94, 044029 (2016), arXiv:1604.05349 [astro-ph.CO] .
  * Stone  _et al._ [2013] E. Stone _et al._ , Science 341, 150 (2013).
  * Cummings  _et al._ [2016] A. C. Cummings, E. C. Stone, B. C. Heikkila, N. Lal, W. R. Webber, G. Jóhannesson, I. V. Moskalenko, E. Orlando, and T. A. Porter, Astrophys. J. 831, 18 (2016).
  * Boudaud and Cirelli [2019] M. Boudaud and M. Cirelli, Phys. Rev. Lett. 122, 041104 (2019), arXiv:1807.03075 [astro-ph.HE] .
  * Laha  _et al._ [2020] R. Laha, J. B. Muñoz, and T. R. Slatyer, Phys. Rev. D 101, 123514 (2020), arXiv:2004.00627 [astro-ph.CO] .
  * Berteaud  _et al._ [2022] J. Berteaud, F. Calore, J. Iguaz, P. D. Serpico, and T. Siegert, Phys. Rev. D 106, 023030 (2022), arXiv:2202.07483 [astro-ph.HE] .
  * Siegert  _et al._ [2022] T. Siegert, C. Boehm, F. Calore, R. Diehl, M. G. H. Krause, P. D. Serpico, and A. C. Vincent, Mon. Not. Roy. Astron. Soc. 511, 914 (2022), arXiv:2109.03791 [astro-ph.HE] .
  * Malyshev  _et al._ [2022] D. Malyshev, E. Moulin, and A. Santangelo, Phys. Rev. D 106, 123020 (2022), arXiv:2208.05705 [astro-ph.HE] .
  * De la Torre Luque  _et al._ [2024a] P. De la Torre Luque, S. Balaji, and P. Carenza, Phys. Rev. D 109, 103028 (2024a), arXiv:2307.13731 [hep-ph] .
  * De la Torre Luque  _et al._ [2024b] P. De la Torre Luque, S. Balaji, and P. Carenza, Phys. Rev. D 109, L101305 (2024b), arXiv:2307.13728 [hep-ph] .
  * De la Torre Luque  _et al._ [2024c] P. De la Torre Luque, S. Balaji, and J. Koechler, Astrophys. J. 968, 46 (2024c), arXiv:2311.04979 [hep-ph] .
  * Cirelli  _et al._ [2023] M. Cirelli, N. Fornengo, J. Koechler, E. Pinetti, and B. M. Roach, JCAP 07, 026 (2023), arXiv:2303.08854 [hep-ph] .
  * Cirelli  _et al._ [2021] M. Cirelli, N. Fornengo, B. J. Kavanagh, and E. Pinetti, Phys. Rev. D 103, 063022 (2021), arXiv:2007.11493 [hep-ph] .
  * Hawking [1975] S. W. Hawking, Commun. Math. Phys. 43, 199 (1975), [Erratum: Commun.Math.Phys. 46, 206 (1976)].
  * Arbey  _et al._ [2020a] A. Arbey, J. Auffinger, and J. Silk, Mon. Not. Roy. Astron. Soc. 494, 1257 (2020a), arXiv:1906.04196 [astro-ph.CO] .
  * Arbey and Auffinger [2019] A. Arbey and J. Auffinger, Eur. Phys. J. C 79, 693 (2019), arXiv:1905.04268 [gr-qc] .
  * Arbey and Auffinger [2021] A. Arbey and J. Auffinger, Eur. Phys. J. C 81, 910 (2021), arXiv:2108.02737 [gr-qc] .
  * Sjöstrand  _et al._ [2015] T. Sjöstrand, S. Ask, J. R. Christiansen, R. Corke, N. Desai, P. Ilten, S. Mrenna, S. Prestel, C. O. Rasmussen, and P. Z. Skands, Comput. Phys. Commun. 191, 159 (2015), arXiv:1410.3012 [hep-ph] .
  * Bellm  _et al._ [2016] J. Bellm _et al._ , Eur. Phys. J. C 76, 196 (2016), arXiv:1512.01178 [hep-ph] .
  * Coogan  _et al._ [2020] A. Coogan, L. Morrison, and S. Profumo, JCAP 01, 056 (2020), arXiv:1907.11846 [hep-ph] .
  * Harada  _et al._ [2021] T. Harada, C.-M. Yoo, K. Kohri, Y. Koga, and T. Monobe, Astrophys. J. 908, 140 (2021), arXiv:2011.00710 [astro-ph.CO] .
  * Harada  _et al._ [2017] T. Harada, C.-M. Yoo, K. Kohri, and K.-I. Nakao, Phys. Rev. D 96, 083517 (2017), [Erratum: Phys.Rev.D 99, 069904 (2019)], arXiv:1707.03595 [gr-qc] .
  * Thorne [1974] K. S. Thorne, Astrophys. J. 191, 507 (1974).
  * Sadowski  _et al._ [2011] A. Sadowski, M. Bursa, M. Abramowicz, W. Kluzniak, J.-P. Lasota, R. Moderski, and M. Safarzadeh, Astron. Astrophys. 532, A41 (2011), arXiv:1102.2456 [astro-ph.HE] .
  * Kesden  _et al._ [2010] M. Kesden, G. Lockhart, and E. S. Phinney, Phys. Rev. D 82, 124045 (2010), arXiv:1005.0627 [gr-qc] .
  * Dolgov and Silk [1993] A. Dolgov and J. Silk, Phys. Rev. D 47, 4244 (1993).
  * Carr  _et al._ [2017] B. Carr, M. Raidal, T. Tenkanen, V. Vaskonen, and H. Veermäe, Phys. Rev. D 96, 023514 (2017), arXiv:1705.05567 [astro-ph.CO] .
  * Balaji  _et al._ [2022] S. Balaji, J. Silk, and Y.-P. Wu, JCAP 06, 008 (2022), arXiv:2202.00700 [astro-ph.CO] .
  * Qin  _et al._ [2023] W. Qin, S. R. Geller, S. Balaji, E. McDonough, and D. I. Kaiser, Phys. Rev. D 108, 043508 (2023), arXiv:2303.02168 [astro-ph.CO] .
  * Boudaud  _et al._ [2017a] M. Boudaud, J. Lavalle, and P. Salati, Phys. Rev. Lett. 119, 021103 (2017a), arXiv:1612.07698 [astro-ph.HE] .
  * Osborne and Ptuskin [1987] J. Osborne and V. Ptuskin, in  _International Cosmic Ray Conference_ , Vol. 2 (1987) p. 218.
  * Seo and Ptuskin [1994] E.-S. Seo and V. S. Ptuskin, ApJ 431, 705 (1994).
  * Ginzburg and Syrovatskii [1969] V. L. Ginzburg and S. I. Syrovatskii,  _Topics in Astrophysics and Space Physics_ (Gordon & Breach Publishing Group, 1969).
  * Evoli  _et al._ [2017] C. Evoli, D. Gaggero, A. Vittino, G. Di Bernardo, M. Di Mauro, A. Ligorini, P. Ullio, and D. Grasso, JCAP 02, 015 (2017), arXiv:1607.07886 [astro-ph.HE] .
  * Evoli  _et al._ [2018] C. Evoli, D. Gaggero, A. Vittino, M. Di Mauro, D. Grasso, and M. N. Mazziotta, JCAP 07, 006 (2018), arXiv:1711.09616 [astro-ph.HE] .
  * Navarro  _et al._ [1996] J. F. Navarro, C. S. Frenk, and S. D. White, Astrophys. J. 462, 563 (1996), arXiv:astro-ph/9508025 .
  * Strong and Moskalenko [1998] A. W. Strong and I. V. Moskalenko, The Astrophysical Journal 509, 212 (1998).
  * Drury and Strong [2017] L. O. Drury and A. W. Strong, Astron. Astrophys. 597, A117 (2017), arXiv:1608.04227 [astro-ph.HE] .
  * Drury and Strong [2015] L. O. Drury and A. W. Strong, “Cosmic-ray diffusive reacceleration: a critical look,” (2015), arXiv:1508.02675 [astro-ph.HE] .
  * De La Torre Luque  _et al._ [2021] P. De La Torre Luque, M. N. Mazziotta, F. Loparco, F. Gargano, and D. Serini, JCAP 03, 099 (2021), arXiv:2101.01547 [astro-ph.HE] .
  * Evoli  _et al._ [2020] C. Evoli, G. Morlino, P. Blasi, and R. Aloisio, Physical Review D 101 (2020), 10.1103/physrevd.101.023013.
  * Weinrich  _et al._ [2020] N. Weinrich, M. Boudaud, L. Derome, Y. Génolini, J. Lavalle, D. Maurin, P. Salati, P. Serpico, and G. Weymann-Despres, Astronomy &Astrophysics 639, A74 (2020).
  * De la Torre Luque  _et al._ [2022] P. De la Torre Luque, M. N. Mazziotta, A. Ferrari, F. Loparco, P. Sala, and D. Serini, JCAP 07, 008 (2022), arXiv:2202.03559 [astro-ph.HE] .
  * Cummings  _et al._ [2016] A. Cummings _et al._ , Astrophys. J. 831, 18 (2016).
  * Arbey  _et al._ [2020b] A. Arbey, J. Auffinger, and J. Silk, “Primordial kerr black holes,” (2020b), arXiv:2012.14767 [astro-ph.CO] .
  * Dundovic  _et al._ [2021] A. Dundovic, C. Evoli, D. Gaggero, and D. Grasso, Astron. Astrophys. 653, A18 (2021), arXiv:2105.13165 [astro-ph.HE] .
  * [58] “XMM-Newton Technical Details, ESA,” http://www.cosmos.esa.int/web/xmm-newton/technical-details-epic.
  * Foster  _et al._ [2021] J. W. Foster, M. Kongsore, C. Dessert, Y. Park, N. L. Rodd, K. Cranmer, and B. R. Safdi, Phys. Rev. Lett. 127, 051101 (2021).
  * Guessoum  _et al._ [2005] N. Guessoum, P. Jean, and W. Gillard, Astron. Astrophys. 436, 171 (2005), arXiv:astro-ph/0504186 .
  * Guessoum  _et al._ [1991] N. Guessoum, R. Ramaty, and R. E. Lingenfelter, The Astrophyssical Journal 378, 170 (1991).
  * De la Torre Luque  _et al._ [2023] P. De la Torre Luque, S. Balaji, and J. Silk, (2023), arXiv:2312.04907 [hep-ph] .
  * De la Torre Luque  _et al._ [2024d] P. De la Torre Luque, S. Balaji, P. Carenza, and L. Mastrototaro, (2024d), arXiv:2405.08482 [hep-ph] .
  * Calore  _et al._ [2021] F. Calore, P. Carenza, M. Giannotti, J. Jaeckel, G. Lucente, and A. Mirizzi, Phys. Rev. D 104, 043016 (2021), arXiv:2107.02186 [hep-ph] .
  * Calore  _et al._ [2022] F. Calore, P. Carenza, M. Giannotti, J. Jaeckel, G. Lucente, L. Mastrototaro, and A. Mirizzi, Phys. Rev. D 105, 063026 (2022), arXiv:2112.08382 [hep-ph] .
  * Carenza  _et al._ [2024] P. Carenza, G. Lucente, L. Mastrototaro, A. Mirizzi, and P. D. Serpico, Phys. Rev. D 109, 063010 (2024), arXiv:2311.00033 [hep-ph] .
  * Laha [2019] R. Laha, Phys. Rev. Lett. 123, 251101 (2019), arXiv:1906.09994 [astro-ph.HE] .
  * Tan and Xia [2024] X.-h. Tan and J.-q. Xia, (2024), arXiv:2404.17119 [astro-ph.CO] .
  * Mittal  _et al._ [2022] S. Mittal, A. Ray, G. Kulkarni, and B. Dasgupta, JCAP 03, 030 (2022), arXiv:2107.02190 [astro-ph.CO] .
  * Kavanagh [2019] B. J. Kavanagh, “bradkav/pbhbounds: Release version,”  (2019).
  * Siegert  _et al._ [2016] T. Siegert, R. Diehl, G. Khachatryan, M. G. H. Krause, F. Guglielmetti, J. Greiner, A. W. Strong, and X. Zhang, Astron. Astrophys. 586, A84 (2016), arXiv:1512.00325 [astro-ph.HE] .
  * Moore  _et al._ [1999] B. Moore _et al._ , The Astrophysical Journal 524, L19–L22 (1999).
  * Ackermann  _et al._ [2017] M. Ackermann _et al._ , The Astrophysical Journal 840, 43 (2017).
  * Burkert [1995a] A. Burkert, Astrophys. J. Lett. 447, L25 (1995a), arXiv:astro-ph/9504041 .
  * Cirelli  _et al._ [2011] M. Cirelli, G. Corcella, A. Hektor, G. Hutsi, M. Kadastik, P. Panci, M. Raidal, F. Sala, and A. Strumia, JCAP 03, 051 (2011), [Erratum: JCAP 10, E01 (2012)], arXiv:1012.4515 [hep-ph] .
  * Maurin [2020] D. Maurin, Computer Physics Communications 247, 106942 (2020).
  * Boudaud  _et al._ [2017b] M. Boudaud, E. F. Bueno, S. Caroff, Y. Genolini, V. Poulin, V. Poireau, A. Putze, S. Rosier, P. Salati, and M. Vecchi, Astron. Astrophys. 605, A17 (2017b), arXiv:1612.03924 [astro-ph.HE] .
  * Merritt  _et al._ [2006] D. Merritt, A. Graham, B. Moore, J. Diemand, and B. Terzić, The Astronomical Journal 132, 2685–2700 (2006).
  * Navarro  _et al._ [2009] J. F. Navarro, A. Ludlow, V. Springel, J. Wang, M. Vogelsberger, S. D. M. White, A. Jenkins, C. S. Frenk, and A. Helmi, Monthly Notices of the Royal Astronomical Society 402, 21–34 (2009).
  * Burkert [1995b] A. Burkert, The Astrophysical Journal 447 (1995b), 10.1086/309560.