MNRAS 000, 1–10 (2021) Preprint 4 January 2022 Compiled using MNRAS LATEX style file v3.0
An INTEGRAL/SPI View of Reticulum II: Particle Dark Matter and
Primordial Black Holes Limits in the MeV Range
Thomas Siegert,1,4? Celine Boehm,2 Francesca Calore,3 Roland Diehl,4 Martin G. H. Krause,5
Pasquale D. Serpico,3 and Aaron C. Vincent6
1Institut für Theoretische Physik und Astrophysik, Universität Würzburg, Campus Hubland Nord, Emil-Fischer-Str. 31, 97074 Würzburg, Germany
2Sydney Consortium for Particle Physics and Cosmology, School of Physics, The University of Sydney, NSW 2006, Australia
3Univ. Grenoble Alpes, Univ. Savoie Mont Blanc, CNRS, LAPTh, F-74940 Annecy, France
4Max-Planck-Institut für extraterrestrische Physik, Gießenbachstraße, 85741 Garching b. München, Germany
5Centre for Astrophysics Research, Department of Physics, Astronomy and Mathematics, University of Hertfordshire, College Lane, Hatfield, Hertfordshire AL10 9AB, UK
6Department of Physics, Enigneering Physics and Astronomy, Queen’s University, Kingston, ON, K7L 3N6, Canada
Accepted XXX. Received 4 January 2022; in original form ZZZ
ABSTRACT
Reticulum II (Ret II) is a satellite galaxy of the Milky Way and presents a prime target to investigate the nature of dark matter
(DM) because of its high mass-to-light ratio. We evaluate a dedicated INTEGRAL observation campaign data set to obtain
γ-ray fluxes from Ret II and compare those with expectations from DM. Ret II is not detected in the γ-ray band 25–8000 keV,
and we derive a flux limit of . 10−8 erg cm−2 s−1. The previously reported 511 keV line is not seen, and we find a flux limit of
. 1.7 × 10−4 ph cm−2 s−1. We construct spectral models for primordial black hole (PBH) evaporation and annihilation/decay of
particle DM, and subsequent annihilation of e+s produced in these processes. We exclude that the totality of DM in Ret II is made
of a monochromatic distribution of PBHs of masses . 8 × 1015 g. Our limits on the velocity-averaged DM annihilation cross
section into e+e− are 〈σv〉 . 5 × 10−28 (mDM/MeV)2.5 cm3 s−1. We conclude that analysing isolated targets in the MeV γ-ray
band can set strong bounds on DM properties without multi-year data sets of the entire Milky Way, and encourage follow-up
observations of Ret II and other dwarf galaxies.
Key words: galaxies: individual: Reticulum II – gamma-rays: galaxies – dark matter – black hole physics
1 INTRODUCTION
Dark matter (DM) is the most abundant component of the matter
content of the Universe. Evidence for DM exists at different scales,
from galaxy to cosmological background radiation, however the very
nature of this elusive component remains a mystery (see Bertone &
Hooper 2018, and references therein for a historical review).Owing to
the vast parameter space for DMcandidates and possible interactions,
different astroparticle observables can be used to tackle the nature of
DM, from cosmic surveys to near field cosmology, from high-energy
radiation to gravitational waves (see Alves Batista et al. 2021, for a
state-of-the-art review). In particular, signals of DM in high-energy
radiation from radio frequencies to γ-rays are expected from different
targets. The choice of the ‘best’ target of interest is determined, on
the one hand, by maximising the expected DM signal, and, on the
other hand, by reducing the usually strong astrophysical fore- and
background.
Despite their abundance, dwarf galaxies are dim objects, some-
times including only a few thousand stars, thus resulting among the
most pristine, chemically unevolved, and possible DM dominated
self-gravitating structures known (Simon 2019). As a result, they
are a primary target for near-field cosmology and indirect searches
? E-mail: tho.siegert@gmail.com
for DM, benefiting of a high (putative) exotic signal to astrophys-
ical background ratio (e.g., Winter et al. 2016). High-energy pho-
tons in the X- and γ-ray band from such objects may come from
the annihilation of conventional DM candidates like weakly inter-
acting massive particles (WIMPs), from decays of alternative DM
candidates like sterile neutrinos, or from the Hawking evaporation
process (Hawking 1971) of sufficiently light black holes, produced
in the early universe. Interestingly, Siegert et al. (2016c) reported
a 3.1σ signal of a 511 keV line from the direction of the dwarf
galaxy Reticulum II (Ret II; distance d = 30 ± 2 kpc; Koposov et al.
2015). The line flux from this previousworkwas extraordinarily high,
F511 = (1.7± 0.5) × 10−3 ph cm−2 s−1, which would make Ret II the
most luminous 511 keV source in the sky. The signal has been inter-
preted in terms of a recent binary neutron star (NS) merger (Fuller
et al. 2018), possibly outshining the entire galaxy. Ret II has also
been reported to show a 2.3–3.7σ excess at GeV energies (Geringer-
Sameth et al. 2015) which might point to a dark matter (DM) origin
of the signal. Follow-up analyses found no flux enhancement above
the diffuse GeV emission (Albert et al. 2017). Intrigued by these
findings, the INTEGRAL (Winkler et al. 2003) satellite performed
a dedicated observation campaign to follow up the signal. In this
work, we analyse the new 1Ms data together with archival data from
the coded-mask γ-ray spectrometer onboard INTEGRAL, SPI (Ve-
drenne et al. 2003). We use our findings to set limits on possible DM
© 2021 The Authors
2 Thomas Siegert et al.
models, including primordial black holes (PBHs) and annihilation or
decay of (light) DM particles.
This paper is structured as follows: In Sec. 2, we describe the new
observations, and how SPI data are analysed in general. We build
photon-emission models expected from the evaporation of PBHs and
annihilating/decaying particle DM in Sec. 3. Here, we also include
the possible annihilation of positrons (e+s) as a DM product into
the galaxy emission spectra. Our results are presented in Sec. 4. We
conclude in Sec. 5.
2 NEW OBSERVATION CAMPAIGN
The previous analysis by Siegert et al. (2016c) focussed on the search
of 511 keV signals from the 39 satellite galaxies of the Milky Way
(MW) known at the time. In their data set, the entire sky was mod-
elled and the diffuse 511 keV emission from the MW itself was
included. Before the new observation campaign, Ret II had never
been targeted for observations by SPI because 1) the galaxy was only
discovered in 2015 (Koposov et al. 2015), and 2) other known X-
or γ-ray sources are separated from Ret II by at least 18◦. However,
thanks to SPI’s wide fully-coded field of view (FCFOV) of 16◦×16◦
(30◦ × 30◦ partially coded), Ret II data were collected also before its
discovery. Outside the full coding region, the effective area of the
instrument sharply drops which placed previous Ret II observations
right onto exposure edges with a cumulative on-target observation
time of 550 ks, split over ten years. While these effect were taken into
account in the analysis, dedicated observations, pointed directly to
Ret II, would be suitable to validate earlier findings.
2.1 Data Set
Our new data set includes previous observations between 2003 and
2013 and adds ∼ 1Ms during the Ret II campaign in 2018. We
selected all pointed observations (pointings) whose distance to Ret II
is less than 16◦ so that our target is always in the FCFOV. This
amounts to 603 pointings with an average observation time of 3300 s
each. The mean dead time of a working detector is about 18.2%.
Since the launch of the INTEGRAL mission in October 2002, four
of the 19 SPI detectors have failed, reducing the effective area by
≈ 20%. The total dead-time- and efficiency-corrected exposure time
then amounts to 1.38Ms. We use SPI’s bandpass between 25 and
8000 keV, and perform our analysis in 29 logarithmic energy bins,
plus a single narrow bin for the 511 keV line from 508 to 514 keV.
2.2 General Analysis Method
SPI data analysis relies on the comparison of expected detector il-
luminations, that is, patterns in the 19-detector data space, with the
measured count rate per pointing and energy. Fluxes are estimated
by maximising the Poisson likelihood,
L (D |M) =
∏
p
m
dp
p exp(−mp)
dp!
, (1)
where dp are the measured counts in pointing p, and mp is the
model expectation. MeV γ-ray measurements are described by a two
component model, one representing the instrumental background
mBG, and one representing the celestial emissionmSKY. Instrumental
background originates from cosmic-ray interactions with the satellite
material and leads to continuum (C) and line (L) backgrounds (Diehl
et al. 2018). Spatial emission models are convolved with the imaging
response function of SPI, that is, the mask coding, leading to a 19-
element vector of expected source counts for each pointing. The total
model reads
mp = mBGp + m
SKY
p =
=
∑
t,t′

∑
i∈{C,L }
βi,tRBGp,i +
NS∑
j=1
αj,t′RSKYp;lbM
lb
j (φ j )
 , (2)
where RBG
p,i
is the backgroundmodel (response) appropriate for point-
ing p, scaled by the amplitudes βi,t for background components
i ∈ {C, L} varying on timescales t, RSKY
p;lb is the coded-mask response
in pointing p, converting the j = 1 . . . NS sky models Mlbj (φ) from
image space lb to the data space.
The source models may have additional parameters φ j , for ex-
ample defining the position of point sources, and are scaled by the
amplitudes αj,t′ , possibly varying on another time scale t ′. Given
SPI’s angular resolution of 2.7◦, and the expected maximum extent
of Ret II’s DM halo of 1◦ (Evans et al. 2016), we model the galaxy
as point source,
MRet II(l, b) = δ(l − l0)δ(b − b0), (3)
with the Galactic coordinates of Ret II (l0, b0) = (266.30,−49.74).
All other sources that can be expected in our data set, EXO0748-
676, ESO33-2, SMCX-1, and LMCX-4 (c.f. Bouchet et al. 2011) are
modelled as point sources as well. The fitted parameters in Eq. (2) are
βi,t and αj,t′ . The fit is performed with spimodfit (Halloin 2009),
which creates a spectrum by looping over the energy bins defined in
our data set. Details about the procedure are found in Siegert et al.
(2019).
2.3 Spectral Fitting
The spectrum created in this way is subject to the instrument disper-
sion, that is, the probability that a photon with initial energy Ei is
measured at a final energy E f owing to scattering in the instrument. In
order to correct for the intrinsic dispersion, expected spectral models
dF
dE are convolved with the energy-redistribution matrix H(Ei, E f )
that forward-folds the models into the appropriate data space. Thus,
the spectral models read
dF
(
E f ;ψ
)
dE f
=
∫
dEiH(Ei, E f ) dF (Ei ;ψ)dEi , (4)
with dFdE f being the folded model which is compared to the extracted
flux values from Sec. 2.2, and dFdEi is the intrinsic source model that
depends on a set of spectral parametersψ. We use theMulti-Mission-
Maximum-Likelihood (3ML, Vianello et al. 2015) framework to
perform spectral fits. Details about the choices of prior probabilities
and parameter ranges are given in AppendixA.
We note that, ideally, the steps explained in Secs. 2.2 and 2.3
should be performed in one single step to enhance the sensitivity of
the parameter estimation, which is, however, not yet implemented in
the current software release.
3 EXPECTATIONS IN THE MEV BAND
No MeV emission has ever been detected from a (dwarf) satellite
galaxy. Therefore, the expected signals have a large variety, ranging
from γ-ray line emission due to radioactive decays and (subsequent)
e+ emission. Population synthesis models estimate a 1.8MeV γ-ray
MNRAS 000, 1–10 (2021)
MeV Dark Matter Limits in Reticulum II 3
line flux from the decay of 26Al of the order of 10−6 ph cm−2 s−1
from the Large Magellanic Cloud (LMC), about one order of magni-
tude below SPI’s sensitivity Diehl et al. (2018). Because the expected
line flux is related to the star formation and supernova rate, any other
satellite galaxy, and in particular Ret II, can be expected to show
no significant excess at 1.8MeV. Cordier et al. (2004) performed a
search for a 511 keV line owing to DM annihilation/decay from the
Sgr dwarf galaxy, but found no excess. Siegert et al. (2016c) ex-
tended this search to all satellite galaxies of the MW, which resulted
in one 3σ signal from the direction of Ret II out of 39 tested galax-
ies. Where this signal, if true, comes from is unknown, and might
point to different origins in terms of NS mergers (Fuller et al. 2018),
accreting X-ray binaries (Siegert et al. 2016a), or DM (Boehm et al.
2004). All these scenarios would show additional emission features
in the MeV band besides a sole 511 keV line. In particular additional
broad γ-ray lines from other nucleosynthesis products inNSmergers,
continuum emission from a population of X-ray binaries, or radia-
tive corrections in the final products from DM decay/annihilation or
PBH decay. Therefore, we extract the spectrum of Ret II and nearby
sources using the logarithmic energy binning defined in Sec. 2.1,
plus a single bin for the 511 keV line. In this way, we can assess the
previously found signal from Ret II in a model-independent fashion,
and determine model dependent parameters using the additionally-
expected components in Sec. 3.1.
3.1 Modelling the MeV Band in Ret II
We focus on the γ-ray spectra expected from the evaporation of
PBHs in the mass range 1014–1018 g (see Sec. 3.1.1) and annihi-
lating/decaying particle DM in the mass range 0.05–300MeV (see
Sec. 3.1.2). These black hole and particle mass ranges would result
in soft γ-ray (MeV) emission with characteristic spectra that INTE-
GRAL/SPI is sensitive to. The absolute fluxes as well as detailed
spectral shapes of such DM signals are strongly model dependent
and are based on unknown astro- and particle-physics input such as
the size of the galaxy’s halo, particle cross sections, or the magnitude
of secondary particles involved. Using literature values (see below),
a DM signal from Ret II that imprints in the MeV photon band would
bewithin reach in our renewed observation data set. Other signals, for
example the γ-ray emission from old NS mergers remnants, are not
part of our predictions because the expected emission of one 10 kyr
old remnant in Ret II is about onemillion times lower (Korobkin et al.
2020) than the instrument sensitivity. While emission from GeV DM
would also imprint in the MeV band, SPI is most sensitive to annihi-
lation/decay to e+e− and subsequent emission processes considering
the final state radiation (FSR) of the pairs as well as the annihilation
of e+s during propagation and after thermalisation. Therefore, we
restrict our models to the two channels e+e− and γγ, and present one
µ+µ− example for comparison. We follow the formalism of Fortin
et al. (2009) throughout the paper.
If the DM microscopic emission process (e.g., annihilation cross-
section 〈σv〉) does not depend on kinetic properties (notably veloc-
ity distribution) the flux from DM haloes factorises into a ‘particle
physics’ (owing the underlying DM model) and an ‘astrophysics’
component (which depends on the distribution of DM in astrophysi-
cal systems) (Bergström et al. 1998) as
dFγ
(
Eγ; κ,M,ψ, n
)
dEγ
=
κ
4piM
dNγ
(
Eγ;ψ
)
dEγ︸                  ︷︷                  ︸
particle physics: PBHorDM
∫
los
drdΩρnhalo(r,Ω)︸                     ︷︷                     ︸
astrophysics: D or J
,
(5)
Model κ M n ψ
PBH fPBH MBH 1→ D
{
fPs, E
max
kin
}
DM Annihilation 〈σv〉 2mDM 2→ J
{
fPs, E
max
kin
}
DM Decay 1/τ mDM 1→ D
{
fPs, E
max
kin
}
Table 1. Parameters for considered DM models in Eq. (5).
where κ, M , and ψ describe the parameters of interest for the specific
cases of PBH evaporation, self-conjugated DM annihilation or decay
(see Tab. 1), and the last term is the D- (n = 1) or J-factor (n = 2)
of the galaxy. They describe the extent and normalisation of the
emission, calculated by an integration of the galaxy’s halo profile
ρhalo(r,Ω) over the line of sight (los). The different scaling for decay
(evaporation) vs. annihiliation originates from the fact that while
for the former process only one particle (PBH) induces the γ-ray
emission, for the latter two particles are needed for the annihilation
to take place, and therefore the emissivity scales as ρ2halo(r,Ω). The
photon spectrum dNγ (Eγ ;ψ)dEγ is related to the photon density of states
as 1Nγ
dNγ
dEγ
= 1κ
dκ
dEγ
so that the number of emitted photons per
process is Nγ =
∫
dEγ
dNγ
dEγ
. Nγ is thus normalising the source
spectrum, Eq. (5), in terms of the totally-emitted photons.
We use the range of D- and J-factors of Ret II from the lit-
erature (Bonnivard et al. 2015; Evans et al. 2016; Albert et al.
2017), in particular D ≈ (1–4) × 1018 GeV cm−2, and J ≈ (0.2–
3.7) × 1019 GeV2 cm−5. The conversion of the D- and J-factor into
other often-used units are D = (0.9–3.5) × 10−2M pc−2 = (1.9–
7.3) × 10−6 g cm−2, and J = (0.6–8.4) × 10−3M2 pc−5 = (0.8–
11.8) × 10−29 g2 cm−5. We use these ranges as uniform priors in our
spectral fits to include uncertainties in the DM halo as well as the
distance to Ret II. Depending on the DM model, the parameters of
interest change and are listed in Tab. 1. Here, fPBH is the fraction
of DM made of PBHs (unitless), 〈σv〉 is the velocity-averaged DM
annihilation cross section in units of cm3 s−1, τ = Γ−1 is the decay
time of a DM particle in units of s, MBH is the mass of PBHs in
units of g, and mDM is the DM particle mass in units of keV. The
additional spectral model parameters ψ =
{
fPs, Emaxkin
}
appear when
e+ annihilation is included (see below).
3.1.1 Primordial Black Hole Evaporation
The spectrum of an evaporating black hole (BH) of mass MBH due
to Hawking radiation (Hawking 1975; Page & Hawking 1976) is(
dNi
dEi
)
BH
=
1
2pi
Γi(Ei,MBH)
exp
(
Ei
TBH
)
− (−1)2si
. (6)
In Eq. (6), Ei is the energy of particle i and si its spin. Γi(Ei,MBH) is
the ‘greybody’ factor that alters the expected blackbody distribution
of possible emitted particles, given the BH temperature
TBH =
~c3
8piGkBMBH
= 6 × 10−8
(
M
MBH
)
K = 1.06
(
1016 g
MBH
)
MeV,
(7)
where ~ is reduced Planck’s constant,G is the gravitational constant,
c is the speed of light, and kB is Boltzmann’s constant.
We use BlackHawk (Arbey & Auffinger 2019) to calculate the
spectra of all relevant particles in the energy range between 1 keV and
1GeV. BlackHawk allows to include secondary particle production
due to hadronization, fragmentation, decay, and other processes as a
result of BH evaporation, which we add to our spectra.
MNRAS 000, 1–10 (2021)
4 Thomas Siegert et al.
101 102 103 104
E  [keV]
10 7
10 6
10 5
10 4
10 3
10 2
10 1
100
E
×
F
 [k
eV
ph
cm
2
s
1
ke
V
1 ]
PBH Evaporation (primary)
PBH Evaporation (secondary)
Positronium Decay
Annihilation in Flight
Annihilation in Flight (weighted)
Internal Bremsstrahlung
Internal Bremsstrahlung (weighted)
Total Spectrum 
(MBH = 1016 g,
 D = 2.5 × 1018 GeV cm 2, 
fPs = 1.0, Ekin, maxe = 4.0 MeV)
102 104
Ekine  [keV]
10 1
101
Ek
in e
×
F e
 [1
0
3 ]
101 102 103 104
E  [keV]
10 7
10 6
10 5
10 4
10 3
10 2
10 1
100
E
×
F
 [k
eV
ph
cm
2
s
1
ke
V
1 ]
PBH Evaporation (primary)
PBH Evaporation (secondary)
Positronium Decay
Annihilation in Flight
Annihilation in Flight (weighted)
Internal Bremsstrahlung
Internal Bremsstrahlung (weighted)
Total Spectrum 
(MBH = 3 × 1015 g,
 D = 2.5 × 1018 GeV cm 2, 
fPs = 0.1, Ekin, maxe = 10.0 MeV)
102 104
Ekine  [keV]
10 1
101
Ek
in e
×
F e
 [1
0
3 ]
Figure 1. Expected γ-ray spectra from PBH evaporation and subsequent an-
nihilation of positron-electron pairs, plus internal bremsstrahlung of emitted
pairs. The two panels show different model parameters (see legends) for a D-
factor of 2.5×1018 GeV cm−2, consistent with Reticulum II. The small insets
show the distribution of electrons and positrons from evaporation, together
with the maximum kinetic energy Emaxkin considered by the shaded area, and
the mean kinetic energy of the electron/positron population (dashed line). The
shaded area of the electron/positron distribution corresponds to the shaded
areas in the γ-ray spectra.
The production of e+s from BH evaporation has no kinematic
threshold based on the mass of the BH. However, the number of e+s
produced that may annihilate subsequently can lead to a 20–8000 keV
flux stronger than the evaporation signal itself which is given by a
BH mass of. 1.1× 1017 g. Given the particle spectrum of e+s from
BlackHawk, the differential e+ flux in a galaxy with D-factor D is
given by
dFe
dEe
=
1
4piMBH
dNe
dEe
D, (8)
so that the total possible e+e−-annihilation flux FAnn in that galaxy
is
FAnn =
∫ Emaxkin
0
dEkin
dFe
dEe,kin
. (9)
Here, Emaxkin is the maximum kinetic energy of a e
+ that annihilates
within the galaxy and does not escape into the intergalactic medium.
For the MW, Emaxkin has been estimated to be around 3–7MeV (Bea-
com & Yüksel 2006; Sizun et al. 2006). We caution, however, that
these estimates may be inaccurate because the statistical methods to
compare the expected spectrumwith the measurements of two differ-
ent instruments are not rigorously correct. Therefore, we leave Emaxkin
as a free parameter in our spectral fits. Because the transport and
environmental conditions in Ret II are largely unknown, this proce-
dure effectively takes into account these uncertainties. For reference,
we discuss the two extreme cases of no and maximal e+ annihila-
tion in our results, Sec. 4, together with the intermediate solution of
marginalising over Emaxkin .
The total annihilation flux is composed of three components,
FAnn = F511 + FoPs + FIA, (10)
where F511 is the flux in the 511 keV line due to direct annihilation
with e−s and intermediate formation of Positronium (Ps), FoPs is the
three-photon decay flux of ortho-Ps, and FIA is the total flux of e+s
annihilating in flight (IA) before stopping/thermalising.
The ortho-Ps and the line flux are related by quantum statistics as
FoPs =
9 fPs
8 − 6 fPs F511 =: r32F511, (11)
with fPs being the Ps fraction, that is, the number of e+s forming
Ps before annihilating (Leventhal et al. 1978). The scaling factor r32
ranges between 0 ( fPs = 0) and 4.5 ( fPs = 1). In general fPs depends
on the annihilation conditions, for example the gas in which e+s anni-
hilate, which is a function of the temperature, density, and ionisation
fraction, among others (e.g., Jean et al. 2006; Churazov et al. 2005,
2011; Siegert et al. 2016b). Since the true conditions in Ret II are
unknown, we leave fPs as free parameters in our spectral fits. The
ortho-Ps spectrum
(
dNγ
dEγ
)
oPs
has been calculated by Ore & Powell
(1949), and we model the 511 keV line,
(
dNγ
dEγ
)
511
as a Gaussian at
511 keV with a width according to SPI’s energy resolution (Diehl
et al. 2018; Attié et al. 2003).
Beacom & Yüksel (2006) derived a relation between FIA and the
511 keV line flux,
FIA =
1
1 − 34 fPs
1 − P
P
F511 =: rIAF511, (12)
where P = P(Ekin,me) is the probability of a e+ with initial kinetic
energy Ekin to annihilate before stopping/thermalising. Therefore,
the total annihilation flux can be expressed as
FAnn = (1 + r32 + rIA)F511. (13)
The survival probability in Eq. (12) is given in general by
P(Ekin,me) = exp
©­­«−nX
∫ Ekin
me
dE
σ(E) dE(E,nX )dx 
ª®®¬ , (14)
where σ(E) is the total annihilation cross section, and
 dE(E,nX )dx  is
the stopping power (energy loss rate, cooling function) of a e+ propa-
gating in a mediumwith density nX . For simplicity, we only consider
Coulomb losses which are expected to dominate up to 100MeV for
MW-like conditions (Jean et al. 2009),
 dE(E,nX )dx Coulomb ∝ nX , so
that nX cancels. The total IA spectrum for mono-energetic injection
of e+s thus reads(
dFγ
dEγ
)mono
IA
= F511
1
1 − 34 fPs
1
P
(
dNγ
dEγ
)
IA
, (15)
with(
dNγ
dEγ
)mono
IA
=
nX
2me
∫ Ekin
E1
dE
P(Ekin, E) dσ(E,Eγ)dE dE(E,nX )dx  . (16)
where the integration limit E1 = max
(
(Eγ−me )2+E2γ
(Eγ−me )+Eγ me, Eγ
)
is
MNRAS 000, 1–10 (2021)
MeV Dark Matter Limits in Reticulum II 5
bound to obey energy conservation (Svensson 1982; Beacom&Yük-
sel 2006). We note that
∫
dEγ
(
dNγ
dEγ
)
IA
= 1 − P. Since in the case
of PBHs, the injection energy is not mono-energetic, we calculate a
weighted photon spectrum given the distribution of kinetic energies
from BlackHawk as(
dNγ
dEγ
)
IA
=
∫ Emaxkin
0 dEkin
dNe
dEe
(
dNγ (Ekin)
dEγ
)mono
IA∫ Emaxkin
0 dEkin
dNe
dEe
. (17)
Radiative corrections in the production of e+e−-pairs can be de-
scribed as internal bremsstrahlung (IB)which can be calculatedwith-
out detailed knowledge of the previous process. This final state radia-
tion (FSR) arises as a result of the production itself and is independent
of the astrophysical environment. The photon source spectrum of IB
per 511 keV line photon is given by Beacom et al. (2005),(
dFγ
dEγ
)mono
IB
=
1
2
1
1 − 34 fPs
(
dNγ
dEγ
)mono
IB
, (18)
with(
dNγ
dEγ
)mono
IB
=
α
2pi
[
M2 + (M − 2Eγ)2
M2Eγ
ln
(
M(M − 2Eγ)
m2e
)]
, (19)
where here, M = Ekin + me is the particles’ injection energy from
PBH evaporation. Because the e+-injection spectrum is again not
mono-energetic, we weight over expected energy distribution such
that(
dNγ
dEγ
)
IB
=
∫ Emaxkin
0 dEkin
dNe
dEe
(
dNγ (Ekin)
dEγ
)mono
IB∫ Emaxkin
0 dEkin
dNe
dEe
. (20)
Finally, the total photon source spectrum of PBH evaporation in-
cluding subsequent e+-annihilation and IB in a galaxy is the sum the
components:(
dFγ
dEγ
) tot
PBH
=
(
dFγ
dEγ
)
PBH
+
(
dFγ
dEγ
)
IB
+
(
dFγ
dEγ
)
Ann
=
=
(
dFγ
dEγ
)
PBH
+
(
dFγ
dEγ
)
IB
+
∑
i∈A
(
dFγ
dEγ
)
i
, (21)
with A = {511, oPs, IA}.
(
dFγ
dEγ
)
PBH
is a function of MBH, D and
fPBH, and
(
dFγ
dEγ
)
Ann
is a function of fPs and Emaxkin . We show two
examples of the total emission spectrum from PBH evaporation in
Ret II in Fig. 1.
3.1.2 Particle Dark Matter Annihilation and Decay
We consider light (mDM . 300MeV) DM particles that either anni-
hilate or decay directly into standard model particles, leading to an
expected photon emission spectrum. Therefore we evaluate the cases
(i) DM + DM→ e+ + e− + (γ)FSR,
(ii) DM + DM→ γ + γ,
(iii) DM→ e+ + e− + (γ)FSR,
(iv) (a) DM→ γ + γ, and
(b) DM→ Y + γ,
whereY is any relativistic particle, for example a sterile neutrino that
decays into γ + ν. The photon spectra of the FSR in cases 1. and 3.
are identical to the IB spectrum,(
dNγ
dEγ
)
FSR
≡
(
dNγ
dEγ
)
IB
, (22)
101 102 103 104
E  [keV]
10 7
10 6
10 5
10 4
10 3
10 2
10 1
100
E
×
F
 [k
eV
ph
cm
2
s
1
ke
V
1 ]
DM Annihilation (FSR)
Positronium Decay
Annihilation in Flight
Total Spectrum 
(MDM = 3.8 MeV, v = 3 × 10 27 cm3 s 1
 J = 2.0 × 1019 GeV2 cm 5, 
fPs = 1.0, Ekin, maxe = 10.0 MeV)
102 104
Ekine  [keV]
10 2
10 1
100
Ek
in e
×
F e
101 102 103 104
E  [keV]
10 7
10 6
10 5
10 4
10 3
10 2
10 1
100
E
×
F
 [k
eV
ph
cm
2
s
1
ke
V
1 ]
DM Annihilation (FSR, e + e )
DM Annihilation (FSR, + )
DM Annihilation ( ) ×10 5
Positronium Decay
Annihilation in Flight
Total Spectrum 
(MDM = 10.3 MeV, v = 3 × 10 24 cm3 s 1
 J = 2.0 × 1019 GeV2 cm 5, 
fPs = N/A, Ekin, maxe = 1 MeV)
102 104
Ekine  [keV]
10 1
101
Ek
in e
×
F e
Figure 2. Expected γ-ray spectra from light dark matter annihilation and
subsequent annihilation of positron-electron pairs, similar to Fig. 1. The two
panels show different model parameters (see legends) for a J-factor of 2 ×
1019 GeV2 cm−5. The lower panel also shows model expectations from final
state radiation of the µ+µ−-channel and direct decay into two photons. No Ps
decay or in flight annihilation would be expected in the latter case.
however now,M = 2mDM for annihilation andM = 1mDM for decay.
Another factor of 12 is required in Eq. (19) for Dirac DM annihila-
tion/decay to account for distinct particle-antiparticle DMwith equal
densities. The expected photon spectrum of direct annihilation/decay
into two photons (cases 2. and 4.) is a delta-function,(
dNγ
dEγ
)
γγ
= 2δ
(
Eγ − M2
)
. (23)
Case 4. (b) differs from case 4. (a) only by an additional factor of 2.
We model Eq. (23) as Gaussians with the energy resolution of SPI.
Similar to the PBH case, the emitted e+s from DM particles in a
galaxy’s halo might annihilate during (IA) and after (Ps) propagation
in the galaxy, or escape into the intergalacticmedium if their injection
energy is larger than a threshold Emaxkin . The only difference is now
that the energy distribution of produced pairs is mono-energetic, so
that either all or none of the pairs lead to additional e+-annihilation
related spectra. Given the equivalence in Eq. (22), one can estimate
the total 511 keV line flux from DM particle annihilation as
FAnn511 =
(
1 − 3
4
fPs
)
J〈σv〉
2pim2DM
, (24)
and from DM particle decay as
FDec511 =
(
1 − 3
4
fPs
)
D
2pimDMτ
. (25)
The remaining spectral components from ortho-Ps and annihilation
MNRAS 000, 1–10 (2021)
6 Thomas Siegert et al.
45
60
75
90
105
120
135
90
       
75
60
45
30
15
LMC X-4
SMC X-1
ESO 33-2
EXO 0748-676
Reticulum II
25-68 keV
45
60
75
90
105
120
135
90
       
75
60
45
30
15
LMC X-4
SMC X-1
ESO 33-2
EXO 0748-676
Reticulum II
504-514 keV
45
60
75
90
105
120
135
90
       
75
60
45
30
15
LMC X-4
SMC X-1
ESO 33-2
EXO 0748-676
Reticulum II
940-2000 keV
Figure 3. Count residuals projected back to the sky by applying the imaging response backwards. The radial and azimuthal coordinates are the Galactic latitude
and longitude, respectively. Shown are the bands 25–68 keV, 508–514 keV, and 940–2000 keV as examples. Only LMC X-4 is significantly detected in the first
energy band.
in flight are the same as in Sec. 3.1.1, except that now the mono-
energetic cases are considered. The total source photon spectrum is
again the sum of the components, here FSR, 511 keV line, ortho-Ps,
and in-flight annihilation. We show examples of the expected DM
annihilation spectra in Fig. 2.
4 RESULTS
We examine our model fits to measured data through a careful anal-
ysis of the residuals in different data space dimensions. Herein, it
proves valuable to apply a back-projection onto the sky of the count
residuals. For cases of background model imperfections, we would
discover irregular or large-scale patterns on the sky, whereas for im-
perfections of the sky model (e.g., missing a point source) would
show up as a single peak at the source’s location. In Fig. 3, we
show the residual count map of an instrumental-background-only fit
in three example energy bands. Only the high-mass X-ray binary
LMCX-4 is detected with a significance of 21σ up to an energy of
∼ 200 keV. The other sources, and in particular Ret II, are not seen
in any energy band. The count residuals show no strong patterns as
expected from an otherwise empty field. In the 511 keV band, the
residuals appear more structured, which is a result of the reduced
count rate compared to other bands. For reference, the particularly
bright spot between ESO33-2 and LMCX-4 at (l, b) = (−78◦,−37◦)
has a significance of ∼ 2.5σ, and is thus considered a background
fluctuation. Details about the 511 keV line from Ret II are given in
Sec. 4.2. If not otherwise stated, upper bounds on fluxes are given as
99.85th percentile (3σ).
4.1 Total Spectrum
We show the spectrum measured by SPI from the direction of Ret II,
spatially modelled as a point source, Eq. (3), in Fig. 4. Because the
source is not detected in any of the energy bins, we plot 2σ upper
limits. For comparison, we show several models that are excluded,
given these fluxes. Detailed exclusion plots are found in the following
sections. From a spectral fit with a generic power-law, ∝
(
Eγ
1MeV
)ν
,
we derive an upper bound on the total flux in the band 20–8000 keV
of 10−8 erg cm−2 s−1.
4.2 511 keV Line
Wefind a 511 keV line flux limit of 1.7×10−4 ph cm−2 s−1. This value
is about one order of magnitude smaller than the 3σ signal reported
101 102 103 104
Energy [keV]
10 2
10 1
100
101
102
103
E2
×
F E
 [k
eV
2
cm
2
s
1
ke
V
1 ]
Powerlaw ( = 0.65)
DM Annihilation (FSR; MDM = 4 MeV,
v = 4 × 10 25 cm3 s 1)
PBH (primaries, MBH = 2 × 1015 g)
PBH (secondaries)
NS merger (×106; @10 kyr)
SPI ( < 2 )
Figure 4. SPI spectrum from the position of Reticulum II. Shown are 2σ
upper limits on the flux as a function of photon energy. A selection of excluded
models is shown (see Secs. 4.3 and 4.4).
by Siegert et al. (2016c) of (1.7 ± 0.5) × 10−3 ph cm−2 s−1. Given
the enhanced exposure time and the source being inside the FCFOV
all the time thanks to the new dedicated observation campaign, the
improvement in sensitivity is plausible.
The previous detection of Ret II might be due to an intrinsic vari-
ability in time, as could be expected from the outburst of a micro-
quasar (Siegert et al. 2016a) or a NS merger (Fuller et al. 2018),
for example. We therefore split the data set into different time bins
according to chance observations of nearby targets as well as the new
observation campaign. The resulting 511 keV light curve is shown
in Fig. 5. Considering only observations in which Ret II is inside the
FCFOV results in upper limits consistent with the previous measure-
ment by Siegert et al. (2016c). We find that the significance of the
511 keV line rises up to ∼ 2σ until the end of the older data set.
Including more data after 2013 and in particular the additional 1Ms
in 2018 (IJD ∼ 6700), leads to a greatly reduced upper flux bound as
well as a lower significance. Assuming the source to be variable in
time, the upper limit on the 511 keV line is ∼ 6.1×10−4 ph cm−2 s−1
– barely consistent with the previous estimate. Another factor here is
the analysed spatial region: While Siegert et al. (2016c) included the
entire sky with the expected diffuse 511 keV emission (four compo-
nents) plus 39 dwarf galaxies, in this work we narrowed the region
of interest and number of components included (Ret II plus 4 nearby
MNRAS 000, 1–10 (2021)
MeV Dark Matter Limits in Reticulum II 7
2000 3000 4000 5000 6000 7000
IJD [MJD-51544]
10 4
10 3
F 5
11
 [p
hc
m
2
s
1 ]
new observation campaign
previous data set
Siegert et al. (2016c)
this work (variable source)
this work (constant source)
SPI data (3  upper limits)
Figure 5. Lightcurve of the 511 keV line from Ret II. No positive excess is
found at any time. Therefore, shown are 3σ upper limits in the band 508–
514 keV. The inferred flux from Siegert et al. (2016c) is shown in comparison
to the new INTEGRAL observation campaign. Additional serendipitous ob-
servations between IJD 4800 and 6600 are also included.
sources). This reduces the covariances between the fitted components
which can result into a slightly changed significance.
We conclude that the 3σ signal fromRet II by Siegert et al. (2016c)
wasmost likely due to an instrumental background fluctuation, paired
with short exposures outside the FCFOV. We can not, however, ex-
clude that the previously reported signal was of astrophysical origin.
4.3 Primordial Black Holes
We perform fits to the extracted spectrum, Fig. 4, with Eq. (21), and
sub-cases thereof. Here we emphasise that we do not solve for the
interesting parameters until a certain figure of merit is above a thresh-
old, but determine the joint posterior distributions of all the free pa-
rameters in Eq. (21). That means we include the uncertainties on D
and estimate fPBH conditioned on all the other unknown parameters,
and illustrate the limits on fPBH as a function ofMBH, as it is typically
done.However herewe consider the 95th percentile in fPBH-direction
of the marginalised posterior pi( fPBH,MBH |yi, σi) with yi and σi as
the flux values and corresponding uncertainties of energy bin i in the
spectrum. See AppendixA for details. In Fig. 6, we show our upper
bounds on fPBH, and compare them to estimates from the recent
literature, including the MW (Laha 2019; Laha et al. 2020), cosmic
rays (Boudaud & Cirelli 2019), the cosmic microwave background
(CMB, Clark et al. 2017), and the cosmic γ-ray background (CGB
Iguaz et al. 2021).
The primary component of PBH evaporation alone cannot con-
strain fPBH in Ret II. Including the secondary particles, the upper
bound on fPBH = 1 excludes monochromatic PBH distributions with
masses of . 0.4 × 1016 g. Assuming that all e+s created by PBHs
in Ret II annihilate strengthens the bound on fPBH, then excluding
masses . 2.9 × 1016 g. Since this assumption is too optimistic, we
included the limiting factor for e+ annihilation, that all e+s above
the threshold Emaxkin escape from the galaxy and do not contribute to
the annihilation spectrum. This describes our most realistic bound
on fPBH = 1, which excludes masses . 0.8 × 1016 g.
Laha (2019) performed a similar estimate for the flux of the
511 keV line in the MW, however assumed that the entire annihi-
lation flux, from integrating over the e+-distribution of PBHs, ends
up in the 511 keV line. Since the Ps fraction in the MW is within
a range 0.97–1.00, the 511 keV line flux from PBHs is at least four
times smaller, because most of the annihilation flux goes into the
three-photon decay of ortho-Ps. In return this means that the bounds
on fPBH from Laha (2019) should be shifted by about a factor of four
1015 1016 1017
MBH [g]
10 1
100
f P
BH
 (9
5t
h 
pe
rc
en
til
e)
PBH (prim.)
PBH (sec.)
PBH (sec.) + Ann.
PBH (sec.) + lim. Ann. + AiF + IB
Laha+2020; MW; SPI
Boudaud+2019; MW; Voyager 1
Iguaz+2021; CGB
Clark+2017; CMB
Figure 6. Fraction of DM that can be constituted of a monochromatic PBH
distribution with mass MBH in Ret II (solid lines), compared to literature
estimates (dashed). The circles mark the masses when the 95th percentile
of the posterior pi( fPBH, MBH |yi, σi ) in each case is reached. Stars mark
the literature limits of 100% PBH DM. The most realistic bound is found
for the case of PBH evaporation, including secondary particles, limited Ps
annihilation, annihilation in flight, and internal bremsstrahlung (purple).
‘vertically’,which decreases themass at fPBH = 1 from∼ 8.9×1016 g
to ∼ 6.4 × 1016 g. We test this by making the same assumptions
with our Ret II spectrum, and find fPBH = 1 is excluded for masses
MBH . 4.0 × 1016 g.
Our conservative and realistic bounds are about one order of mag-
nitude in PBHmass below previous estimates. Nevertheless, our data
set is considerably shorter than all the other data sets, often compris-
ing more than a decade of observation time. Thanks to the much
smaller region of interest in the case of Ret II, basically constituting
a point source of SPI, our systematic uncertainties can be consid-
ered smaller compared to studies of the MW. The entire Galactic
(DM) profile is rather uncertain, and also the Galactic centre itself
bears problems when DM models are applied, as is often the case
for diffuse emission in the GeV band. While our bounds are not yet
as constraining, we show that reasonable progress can be made with
even a small data set. In particular, we show that when all uncertain-
ties on the spectral parameters are included also more conservative
and realistic bounds can result.
4.4 Particle Dark Matter
Similar to the PBH case, we perform a fit to the extracted Ret II spec-
trum and marginalise over the uncertain parameters to determine the
posteriors pi(〈σv〉,mDM |yi, σi) and pi(τ,mDM |yi, σi) for DM annihi-
lation and decay, respectively. In what follows, the excluded regions
consider the 95th percentile of the posteriors in 〈σv〉- and τ-direction,
respectively.
4.4.1 Annihilation
In Fig. 7, we show the excluded region of the velocity-averaged DM
annihilation cross section into e+e−, followed by FSR and subsequent
e+ annihilation. Without additional e+ annihilation, the bounds from
Ret II are barely comparable to literature values from other γ-ray
experiments (Essig et al. 2013) or the CMB (Slatyer 2016). When
limited e+ annihilation is included, we can improve the limits in the
MNRAS 000, 1–10 (2021)
8 Thomas Siegert et al.
100 101 102
mDM [MeV]
10 31
10 29
10 27
10 25
10 23
10 21
10 19
v
 [c
m
3
s
1 ]
 (9
5t
h 
pe
rc
en
til
e)
+  (FSR, this work)
e + e  (FSR, this work)
e + e  (FSR + Ps, this work)
e + e  (FSR, Essig+2013; MW; SPI)
e + e  (FSR, Essig+2013; MW; COMPTEL)
e + e  (CMB, Slatyer+2016; Planck)
Figure 7. DM annihilation cross section into e+e− with subsequent FSR and
possible annihilation of the pairs in Ret II (case 1. in Sec. 3.1.2; gray lines),
compared to literature limits. The shaded regions are excluded.
10 2 10 1 100 101 102
mDM [MeV]
10 32
10 30
10 28
10 26
10 24
10 22
10 20
v
 [c
m
3
s
1 ]
 (9
5t
h 
pe
rc
en
til
e)
 (this work)
 (Laha et al. 2020; SPI, MW)
 (Ng et al. 2019; NuSTAR, M31)
 (Slatyer et al. 2016; Planck, CMB)
Figure 8. Similar to Fig. 7 but for DM + DM→ γ + γ.
DM mass range 0.5–1.0MeV by up to two orders of magnitude. We
find a general trend of the annihilation cross section at tree level of
〈σv〉 . 5×10−28 (mDM/MeV)2.5 cm3 s−1, based on the shape of the
measured limits as a function of energy in Fig. 7. For completeness,
we show the limits on DM annihilation into µ+µ−, being hardly
constrained by Ret II data.
The cross section limit for the case DM + DM→ γ + γ is shown
in Fig. 8. Our limits on DM annihilation into two photons are compa-
rable to those from NuSTAR observations of M31 (Ng et al. 2019),
and about one order of magnitude worse than fromMWobservations
with SPI (Laha et al. 2020). Because our data set extends to the full
energy range of SPI, we can set however limits out tomDM ≤ 8MeV,
varying between 10−28 and 10−25 cm3 s−1. However, limits from the
CMB (Slatyer 2016) are still more stringent at there energies.
4.4.2 Decay
The properties of decaying DM particles cannot be inferred directly
from the previous case of annihilating DM. First, the kinematic
thresholds for standard model particle production is shifted by a
factor of two which leads to altered (shifted) spectral shapes in par-
100 101 102
mDM [MeV]
1020
1021
1022
1023
1024
1025
1026
1027
 [s
] (
95
th
 p
er
ce
nt
ile
)
e + e  (FSR, this work)
e + e  (FSR + Ps, this work)
e + e  (FSR, Essig+2013; MW; SPI)
e + e  (FSR, Essig+2013; MW; COMPTEL)
e + e  (CMB, Slatyer+2017; Planck)
Figure 9. Similar to Fig. 7 but for DM→ e+ + e− + (γ)FSR.
10 2 10 1 100 101 102
mDM [MeV]
1017
1019
1021
1023
1025
1027
1029
 [s
] (
95
th
 p
er
ce
nt
ile
)
 (this work)
 (Laha+2020; MW; SPI)
 (Essig+2013; MW; HEAO)
 (Essig+2013; MW; SPI)
 (Essig+2013; MW; COMPTEL)
 (Ng+2019; M31; NuSTAR)
 (CMB, Slatyer+2017; Planck)
Figure 10. Similar to Fig. 7 but forDM→ γ + γ. For the caseDM→ γ + Y,
the limits are smaller by a factor of 2.
ticular for FSR and IA. Second, the J-factor of Ret II is about twice as
uncertain as its D-factor, which allows to populate different regions
in the full posterior. We therefore perform spectral fits for the DM
decay case separately to infer bounds on the DM lifetime τ. This also
results in slightly worse bounds when considering decay compared
to annihilation (cf. NuSTAR limits).
In Fig. 9 we show the exclusion region of τ as a function of the
DM mass. Similar to the DM annihilation case, FSR alone is barely
constraining, but can be pushed by including e+ annihilation. We
find that DM particles decaying into e+e− pairs must have a lifetime
longer than 1–5× 1024 s between 1 and 2MeV. In the mass range 1–
2MeV, our bounds are as strong as the limits from the CMB Slatyer
& Wu (2017). Above this mass range, previous limits from SPI and
COMPTEL in the MW (Essig et al. 2013), and in particular those
from the CMB are more constraining.
Finally, our DM decay time limits considering the direct decay
into two photons is not improving beyond previous constraints. This
is expected from the pure narrow line sensitivity of SPI.
MNRAS 000, 1–10 (2021)
MeV Dark Matter Limits in Reticulum II 9
5 SUMMARY AND CONCLUSION
Ret II is not seen in soft γ-rays between 0.025 and 8MeV, up to
an upper flux limit of 10−8 erg cm−2 s−1 (5 × 10−10 erg cm−2 s−1
within 20–80 keV). The previously reported 511 keV line signal in
Ret II has not been consolidated in this work, and we provide an
upper limit on its flux of 1.7 × 10−4 ph cm−2 s−1. Given the distance
to the galaxy, this constrains the e+-annihilation rate (including Ps
formation and annihilation in flight) to 0.9–3.7×1043 e+ s−1, at most
the annihilation rate of the MW1.
We calculated parametrised spectral models for PBH evaporation
and DM annihilation/decay with subsequent e+ annihilation and fit-
ted these models to the total spectrum from Ret II. Our conservative
limit on the monochromatic PBH mass distribution constituting the
entirety of DM in Ret II is 0.8 × 1016 g. We improve the limits on
the velocity-averaged DM annihilation cross section of decay time
of DM particles into e+e− in the range 0.5–2MeV. For annihilation
or decay into two photons, our limits are weaker than previous es-
timates. We provide the currently only limits on the cross section
〈σv〉γγ in the range 3–8MeV based on γ-ray observations, ranging
between 10−28 and 10−25 cm3 s−1.
While our estimates are at most on the same level as previous
constraints, we arrive at our conclusions with a very small data set
(two weeks) compared to the multi-year or even decade-long obser-
vations. With our detailed spectral modelling including the effects of
possible e+ annihilation and the proper statistical treatment of known
unknowns, we show that dedicated observations of selected targets
provide valuable information about the DM conundrum. We encour-
age follow-up observations of Ret II and other dwarf galaxies with
INTEGRAL, and foresee a bright future for soft γ-ray instruments to
come, such as COSI (Tomsick et al. 2019), to contribute significantly
in the research for the nature of DM. With its more than one order
of magnitude increased sensitivity, in the 0.5–5MeV band COSI can
reach and surpass the strongest CMB limits within its nominal two-
year mission from observations of dwarf galaxies only. Thanks to
its better spatial resolution as well as its more uniform coverage of
the MeV sky, disentangling fore- and background components of the
Milky Way itself will result in a clearer separation of putative DM
signals from the still weakly constrained MeV spectrum.
SOFTWARE
OSA/spimodfit (Halloin 2009), BlackHawk (Arbey & Auffinger
2019), numpy (Oliphant 2006), matplotlib (Hunter 2007), astropy
(Collaboration et al. 2013), scipy (Virtanen et al. 2019), 3ML
(Vianello et al. 2015),MultiNest (Feroz & Hobson 2008; Feroz et al.
2009, 2019).
ACKNOWLEDGEMENTS
Thomas Siegert is supported by the German Research Foundation
(DFG-Forschungsstipendium SI 2502/3-1). We thank Ranjan Laha,
Kenny Ng, Shunsaku Horiuchi, Marco Cirelli, Sam McDermott,
Tracy Slatyer, and Joaquim Iguaz for providing limits on dark matter
properties, and Oleg Korobkin for γ-ray spectra from neutron star
1 We note that in Siegert et al. (2016c), there is a typo, reducing the estimate
for Ret II by a factor of ten – the actual annihilation rate in the previous work
should be on the order of 1044 e+ s−1.
mergers. We thank John Beacom and Hassan Yüksel for details on
the in-flight annihilation spectrum.
DATA AVAILABILITY
The data underlying this article will be shared on reasonable request
to the corresponding author.
REFERENCES
Albert A., et al., 2017, The Astrophysical Journal, 834, 110
Alves Batista R., et al., 2021, arXiv.org, p. arXiv:2110.10074
Arbey A., Auffinger J., 2019, The European Physical Journal C, 79, 693
Attié D., et al., 2003, Astronomy & Astrophysics, 411, L71
Bayes T., 1763, Proceedings of the Royal Society of London, Philosophical
Transactions of the Royal Society, 53, 370
Beacom J. F., Yüksel H., 2006, Physical Review Letters, 97, 071102
Beacom J. F., Bell N. F., Bertone G., 2005, Physical Review Letters, 94,
171301
Bergström L., Ullio P., Buckley J. H., 1998, Astroparticle Physics, 9, 137
Bertone G., Hooper D., 2018, Reviews of Modern Physics, 90, 045002
Boehm C., Hooper D., Silk J., Cassé M., Paul J., 2004, Physical Review
Letters, 92, 101301
Bonnivard V., et al., 2015, The Astrophysical Journal Letters, 808, L36
Bouchet L., Strong A.W., Porter T. A., Moskalenko I. V., Jourdain E., Roques
J.-P., 2011, The Astrophysical Journal, 739, 29
Boudaud M., Cirelli M., 2019, Physical Review Letters, 122, 041104
Churazov E., Sunyaev R., Sazonov S., Revnivtsev M., Varshalovich D., 2005,
Monthly Notices of the Royal Astronomical Society, 357, 1377
Churazov E., Sazonov S., Tsygankov S., Sunyaev R., Varshalovich D., 2011,
Monthly Notices of the Royal Astronomical Society, 411, 1727
Clark S. J., Dutta B., Gao Y., Strigari L. E., Watson S., 2017, Physical Review
D, 95, 083006
Collaboration A., et al., 2013, Astronomy & Astrophysics, 558, A33
Cordier B., et al., 2004, in Schoenfelder V., Lichti G., Winkler C., eds, 5th
INTEGRAL Workshop on the INTEGRAL Universe. p. 581
Diehl R., et al., 2018, Astronomy & Astrophysics, 611, A12
Essig R., Kuflik E., McDermott S. D., Volansky T., Zurek K. M., 2013,
Journal of High Energy Physics, 11, 193
Evans N. W., Sanders J. L., Geringer-Sameth A., 2016, arXiv.org, 93, 103512
Feroz F., Hobson M. P., 2008, Monthly Notices of the Royal Astronomical
Society, 384, 449
Feroz F., Hobson M. P., Bridges M., 2009, Monthly Notices of the Royal
Astronomical Society, 398, 1601
Feroz F., Hobson M. P., Cameron E., Pettitt A. N., 2019, The Open Journal
of Astrophysics, 2, 10
Fortin J.-F., Shelton J., Thomas S., Zhao Y., 2009, arXiv.org, p.
arXiv:0908.2258
Fuller G. M., Kusenko A., Radice D., Takhistov V., 2018, arXiv.org
Geringer-Sameth A., Walker M. G., Koushiappas S. M., Koposov S. E.,
Belokurov V., Torrealba G., Evans N. W., 2015, Physical Review Letters,
115, 081101
Halloin H., 2009, |spimodfit| Explanatory Guide and Users Manual. Max
Planck Institut für extraterrestrische Physik Max Planck Institut für ex-
traterrestrische Physik, Giessenbachstraße 1, 85748 Garching, Germany,
version 2.9 edn
Hawking S., 1971, Monthly Notices of the Royal Astronomical Society, 152,
75
Hawking S. W., 1975, Communications In Mathematical Physics, 43, 199
Hunter J. D., 2007, Computing in Science & Engineering, 9, 90
Iguaz J., Serpico P. D., Siegert T., 2021, arXiv.org, p. arXiv:2104.03145
Jean P., Knoedlseder J., GillardW., Guessoum N., Ferrière K., Marcowith A.,
Lonjou V., Roques J. P., 2006, Astronomy & Astrophysics, 445, 579
Jean P., Gillard W., Marcowith A., Ferrière K., 2009, Astronomy & Astro-
physics, 508, 1099
MNRAS 000, 1–10 (2021)
10 Thomas Siegert et al.
Koposov S. E., Belokurov V., Torrealba G., Evans N. W., 2015, arXiv.org,
805, 130
Korobkin O., et al., 2020, The Astrophysical Journal, 889, 168
Laha R., 2019, Physical Review Letters, 123, 251101
Laha R., Muñoz J. B., Slatyer T. R., 2020, Physical Review D, 101, 123514
Leventhal M., MacCallum C. J., Stang P. D., 1978, Astrophysical Journal,
225, L11
Ng K. C. Y., Roach B. M., Perez K., Beacom J. F., Horiuchi S., Krivonos R.,
Wik D. R., 2019, Physical Review D, 99, 083005
Oliphant T. E., 2006, A guide to NumPy. Vol. 1, Trelgol Publishing USA
Ore A., Powell J. L., 1949, Physical Review, 75, 1696
Page D. N., Hawking S. W., 1976, Astrophysical Journal, 206, 1
Siegert T., et al., 2016a, Nature, 531, 341
Siegert T., Diehl R., Khachatryan G., Krause M. G. H., Guglielmetti F.,
Greiner J., Strong A. W., Zhang X., 2016b, Astronomy & Astrophysics,
586, A84
Siegert T., Diehl R., Vincent A. C., Guglielmetti F., Krause M. G. H., Boehm
C., 2016c, Astronomy & Astrophysics, 595, A25
Siegert T., Diehl R., Weinberger C., Pleintinger M. M. M., Greiner J., Zhang
X., 2019, Astronomy & Astrophysics, 626, A73
Simon J. D., 2019, Annual Review of Astronomy and Astrophysics, 57, 375
Sizun P., Cassé M., Schanne S., 2006, Phys. Rev. D, 74, 063514
Slatyer T. R., 2016, Physical Review D, 93, 023527
Slatyer T. R., Wu C.-L., 2017, Physical Review D, 95, 023010
Svensson R., 1982, ] 10.1086/160081, 258, 321
Tomsick J. A., et al., 2019, arXiv.org, p. arXiv:1908.04334
Vedrenne G., et al., 2003, Astronomy & Astrophysics, 411, L63
Vianello G., et al., 2015, arXiv.org, p. arXiv:1507.08343
Virtanen P., et al., 2019, arXiv.org, p. arXiv:1907.10121
Winkler C., et al., 2003, Astronomy & Astrophysics, 411, L1
Winter M., Zaharijaš G., Bechtol K., Vandenbroucke J., 2016, The Astro-
physical Journal Letters, 832, L6
APPENDIX A: DETAILS ON SPECTRAL FITS
The error bars in the extracted spectrum in Sec. 4.1, Fig. 4 approxi-
mately follow a normal distribution. Therefore, we use the likelihood
Lnormal(D |M(ψ)) =
30∏
i=1
1√
2piσi
exp
(
−1
2
[
yi − mi(ψ)
σi
]2)
, (A1)
with yi as themeasured flux in energy bin i = 1 . . . 30 in the spectrum,
σi the corresponding uncertainties, and mi(ψ) the forward-folded
model, Eq. (4), that depends on a set of spectral parameters ψ.
Since we want to include the astrophysical and particle physics
uncertainties of our source, we construct the joint posterior of all N
parameters ψ1 . . . ψN through Bayes’ theorem (Bayes 1763) as
pi(ψ1 . . . ψN |yi, σi) ∝ Lnormal(yi, σi |ψ1 . . . ψN )pi(ψ1 . . . ψN ).
(A2)
Here, pi(ψ1 . . . ψN ) is the joint prior distribution of the individ-
ual parameters. We use independent priors so that pi(ψ1 . . . ψN ) =
pi(ψ1) · · · pi(ψN ), that is, each parameter obtains its individual prior
probability distribution function. By integrating out the parameters
ψ3 to ψN , for example, we can construct the marginalised joint pos-
terior distribution of ψ1 and ψ2,
pi(ψ1, ψ2 |yi, σi) =
∫
· · ·
∫
dψ2 · · ·ψN pi(ψ1 . . . ψN |yi, σi). (A3)
To obtain the upper bound on ψ1 as a function of ψ2, we integrate the
marginalised posterior pi(ψ1, ψ2 |yi, σi) in ψ1-direction and equate to
103 104 105 106
mDM [keV c 2]
10 29
10 27
10 25
10 23
10 21
v
ee
 [c
m
3
s
1 ]
v ub(mDM)
MultiNest Samples
Figure A1. Marginalised posterior distribution pi(〈σv〉,mDM |yi, σi ) (red
dots), together with upper bounds of the annihilation cross section as a func-
tion of the DM mass (black line).
a certain percentile P,
P =
∫ ψub1 (ψ2)
ψmin1
dψ1pi(ψ1, ψ2 |yi, σi). (A4)
Solving Eq. (A4) for ψub1 (ψ2) provides the upper bound on ψ1 as
function of ψ2.
A written example considering the case of annihilating DM into
e+e− with subsequent FSR and e+ annihilation would calculate
the upper bound on the velocity-averaged annihilation cross sec-
tion 〈σv〉 as a function of the DM mass mDM. Thus first, the
join posterior pi(〈σv〉,mDM, J, fPs, Emaxkin |yi, σi) is estimated, then
marginalised over the J-factor, the Ps fraction and the kinetic energy
threshold to avoid escape,
pi(〈σv〉,mDM |yi, σi) =∫ Jmax
Jmin
dJ
∫ fPs,max
fPs,min
dfPs
∫ Emaxkin,max
Emaxkin,min
dEmaxkin
pi(〈σv〉,mDM, J, fPs, Emaxkin |yi, σi), (A5)
and finally integrated up to the bound 〈σv〉ub(mDM) by
P =
∫ 〈σv〉ub(mDM)
0
d〈σv〉pi(〈σv〉,mDM |yi, σi). (A6)
We use MultiNest (Feroz & Hobson 2008; Feroz et al. 2009, 2019)
in 3ML to evaluate the joint posterior distributions, marginalisations,
and upper bounds. We show the samples of this example in Fig. A1.
We list our prior distributions and parameter ranges for all considered
models in Tab. A1.
This paper has been typeset from a TEX/LATEX file prepared by the author.
MNRAS 000, 1–10 (2021)
MeV Dark Matter Limits in Reticulum II 11
Model κ M D J fPs Emaxkin
PBH (prim.) U (0, 1) logU (1014, 1018) U (1018, 4 × 1018) – – –
PBH (sec.)
.
.
.
.
.
.
.
.
. – – –
PBH (sec.) + Ann.
.
.
.
.
.
.
.
.
. – U (0, 1) –
PBH (sec.) + lim. Ann.
.
.
.
.
.
.
.
.
. –
.
.
. logU (103, 105)
2DM→ e+e−γ logU (3 × 10−34, 3 × 10−20) logU (511, 106) – U (2 × 1018, 4 × 1019) U (0, 1) logU (103, 105)
2DM→ 2γ
.
.
. logU (25, 8000) –
.
.
. – –
DM→ e+e−γ logU (1015, 1024) logU (1022, 106) U (1018, 4 × 1018) – U (0, 1) logU (103, 105)
DM→ 2γ logU (1015, 1029) logU (50, 16000)
.
.
. – – –
Table A1. Prior distributions for considered DM models in Eq. (5).
MNRAS 000, 1–10 (2021)