11institutetext: Hamburger Sternwarte, University of Hamburg, Gojenbergsweg
112, D-21029 Hamburg, Germany 22institutetext: INAF - Istituto di
Radioastronomia, via P. Gobetti 101, 40129, Bologna, Italy 33institutetext:
Institut für Theoretische Physik, Auf der Morgenstelle 14, Eberhard Karls
Universität Tübingen, 72076 Tübingen, Germany 44institutetext: Thüringer
Landessternwarte, Sternwarte 5, 07778 Tautenburg, Germany 55institutetext:
Max-Planck Institut für Radioastronomie, Auf dem Hügel 69, 53121 Bonn, Germany
66institutetext: Fakultät für Physik, Universität Bielefeld, Postfach 100131,
33501 Bielefeld, Germany 77institutetext: Ruhr University Bochum, Faculty of
Physics and Astronomy, Astronomical Institute, 44780 Bochum, Germany

# WIMP cross-section limits from LOFAR observations of dwarf spheroidal
galaxies

L. Gajović 11 F. Welzmüller 11 V. Heesen 11 F. de Gasperin 2211 M. Vollmann 33
M. Brüggen 11 A. Basu 44 R. Beck 55 D. J. Schwarz 66 D. J. Bomans 77 A.
Drabent 44

(Received xxx; accepted xxx)

###### Abstract

Context. Weakly interacting massive particles (WIMPs) can self-annihilate and
thus provide us with the possibility for an indirect detection of Dark Matter
(DM). Dwarf spheroidal (dSph) galaxies are excellent places to search for
annihilation signals because they are rich in DM and background emission is
low. If O(0.1–10$\leavevmode\nobreak\ \mu$G) magnetic fields in dSph exist,
the particles produced in DM annihilation emit synchrotron radiation in the
radio band.

Aims. We use the non-detection of 150 MHz radio continuum emission from dSph
galaxies with the LOw Frequency ARray (LOFAR) to derive constraints on the
annihilation cross-section of WIMPs into electron–positron pairs. Our main
underlying assumption is that the transport of the CRs can be described by the
diffusion approximation which necessitates the existence of magnetic fields.

Methods. We use observations of six dSph galaxies in the LOFAR Two-metre Sky
Survey (LoTSS). The data are re-imaged and a radial profile is generated for
each galaxy. We also use stacking to increase the sensitivity. In order to
derive upper limits on the WIMP cross-section, we inject fake Gaussian sources
into the data which are then detected with 2$\sigma$ significance in the
radial profile. These sources represent the lowest emission we would have been
able to detect.

Results. We present limits from the observations of individual galaxies as
well as from stacking. We explore the uncertainty due to the choice of
diffusion and magnetic field parameters by constructing three different model
scenarios: optimistic (OPT), intermediate (INT), and pessimistic (PES).
Assuming monochromatic annihilation into electron-positron pairs, the limits
from the INT scenario exclude thermal WIMPs ($\langle\sigma v\rangle\\!\approx
2.2\times 10^{-26}\leavevmode\nobreak\ \rm cm^{3}\,s^{-1}$) below 20 GeV and
the limits from the OPT scenario even exclude thermal WIMPs below 70 GeV. The
INT limits can compete with limits set by Fermi–LAT using $\gamma$-ray
observations of multiple dwarf galaxies and they are especially strong for low
WIMP masses.

###### Key Words.:

astroparticle physics – dark matter – radio continuum: galaxies – galaxies:
dwarf – galaxies: magnetic fields

##  1 Introduction

Dark Matter (DM) is known to interact gravitationally and only weakly via
other fundamental forces of the nature, which makes it difficult to observe.
Among the most promising candidates for DM are weakly interacting massive
particles (WIMPs; Jungman et al., 1996), (QCD) axions (Peccei & Quinn, 1977a,
b; Weinberg, 1978; Wilczek, 1978) or axion-like particles (Kim, 1987; Jaeckel
& Ringwald, 2010), massive compact halo objects (MACHOs; Alcock et al., 2000),
sterile neutrinos (Ibarra, 2015), primordial black holes (Hawking, 1971), and
modification of the Newtonian dynamics (MOND) as an alternative to explain the
effect of DM without additional mass (Milgrom, 1983). WIMPs are particularly
appealing candidates for DM and by far the most scrutinized.

In the WIMP hypothesis, the number density of DM particles freezes out when
the expansion rate of the Universe becomes higher than their annihilation rate
so the number density of DM particles becomes constant. From this, the
theoretical thermal relic annihilation cross-section is calculated to be
$\langle\sigma v\rangle\\!\approx 2.2\times 10^{-26}\leavevmode\nobreak\ \rm
cm^{3}\,s^{-1}$ for WIMP masses above 10 GeV and predicted to increase at
lower masses (Steigman et al., 2012). There are also alternative production
scenarios such as freeze-in where DM never attains thermal equilibrium with
the primordial plasma of Standard Model particles in the early Universe (Hall
et al., 2010), as well as several scenarios where DM is in thermal
equilibrium: hidden sector freeze-out (Finkbeiner & Weiner, 2007; Cheung et
al., 2011), Zombie DM (Kramer et al., 2021), Pandemic DM (Bringmann et al.,
2021; Hryczuk & Laletin, 2021), and others.

There is ongoing research attempting to directly detect WIMPs, for example
EDELWEISS (Sanglard et al., 2005), the XENON-1T experiment (Aprile et al.,
2017), and CRESST (Abdelhameed et al., 2019). Another way to achieve a
detection is an indirect astrophysical search targeting particles that were
created by or interacted with DM. Searches using $\gamma$-ray observations of
nearby dwarf spheroidal (dSph) galaxies using the Fermi Large Area Telescope
(Fermi–LAT; Ackermann et al., 2015; Hoof et al., 2020) and the High Energy
Stereoscopic System (HESS; Aharonian et al., 2008) were able to provide new
stringent upper limits for the WIMP annihilation cross-section depending on
the particle mass. Recent work by Delos & White (2022) shows that if DM is
thermally produced WIMPs, multiple thousands of earth mass prompt cusps should
be present in every solar mass of DM. This is manifested in a drastic increase
of the expected DM annihilation signal compared to a smooth DM distribution,
substantially tightening observational cross-section limits. The most recent
$\gamma$-ray results (Hoof et al., 2020), assuming a smooth DM distribution,
exclude the thermal relic annihilation cross-section for WIMPs with masses of
$\lesssim$ 100 GeV annihilating through the quark and $\tau$-lepton channels;
and for WIMPs with masses of $\lesssim$ 20 GeV annihilating through the
electron–positron channel.

Radio continuum observations provide a complementary approach to constraining
the WIMP annihilation cross-section. Radio continuum searches exploit the fact
that highly energetic electron–positron pairs produced by the annihilation of
WIMPs give rise to synchrotron emission in the presence of magnetic fields
(Colafrancesco et al., 2007). Although other WIMP annihilation channels can be
explored, in this work we focus on WIMP self-annihilation into
electron–positron pairs. For WIMP masses on the order of 10 GeV and microgauss
($\mu$G) magnetic fields the critical frequency of synchrotron radiation
$\nu=3eBE^{2}/(2\pi m_{e}^{3}c^{5})$ is on the order of a few gigahertz
(Rybicki & Lightman, 1986), which means the WIMP signal can only be detected
at frequencies lower than that. Lowering the frequency below 1 GHz opens up
the possibility for constraining WIMPs with masses down to 1 GeV.

Cosmic ray (CR) electrons and positrons produced via DM annihilation interact
with the magnetic field and create an extended source tracing a diffusion
halo, due to the emitted synchrotron radiation. A common search strategy is to
consider radial intensity profiles of the extended emission from these halos
and compare them with modeled DM annihilation signals (Cook et al., 2020;
Vollmann et al., 2020). Vollmann (2021) presents semi-analytical formulae for
these models which were shown to be in reasonable agreement with the more
sophisticated numerical methods (Regis et al., 2015, 2021).

Dwarf spheroidal (dSph) galaxies are very promising systems to look for DM-
related emission, as in these systems the radio continuum emission is almost
uncontaminated by baryonic emission of primary CR electrons due to a low level
of star formation (Colafrancesco et al., 2007; Heesen & Brüggen, 2021). Also,
stellar-dynamical observations indicate that dSphs are DM dominated (Geringer-
Sameth et al., 2015), meaning that the signal from WIMP annihilation is
expected to be bright compared to other types of galaxies. Most known dSph
galaxies are satellites of either the Milky Way or the Andromeda galaxy and
therefore close to Earth. This allows us to well resolve the spectra of their
stellar light and the estimates of their DM content are relatively precise
(Hütten & Kerszberg, 2022).

The LOw Frequency ARay (LOFAR; van Haarlem et al., 2013) is the ideal
instrument for radio continuum searches of WIMP annihilation. LOFAR is a radio
interferometer operating at low frequencies from 10 to 240 MHz. LOFAR combines
high angular resolution needed to identify compact background sources with
high sensitivity to extended emission needed to detect the signal from WIMPs.
A LOFAR high band antenna (HBA) search for DM in the dSph galaxy Canes
Venatici I using such a strategy has already been performed by Vollmann et al.
(2020). That proof-of-concept study showed that for an individual galaxy the
limits are comparable or even better than that from Fermi–LAT, assuming
reasonable values for the magnetic field strength and diffusion coefficient.
This work expands that study by considering six galaxies observed with LOFAR
HBA. In order to improve the WIMP annihilation cross-section limits, we stack
the signal form these six galaxies in various ways. We thus present improved
upper limits on the annihilation of low-mass WIMPs into electron–positron
pairs in the GeV–TeV mass range using new radio continuum data.

This work is organised as follows. In Sect. 2, we present our theoretical
calculations, where we convert the presence of DM into a radio continuum
signal. Section 3 presents our employed methodology including how we dealt
with the LOFAR radio continuum data and how we improved our signal-to-noise
ratio by both profile and image stacking. Section 4 contains our results. We
finish off in Sect. 5 which contains the discussion and conclusions.

##  2 Dark matter calculations

###  2.1 Cosmic ray injection

The injection of Standard Model particles by WIMP annihilation was described
by for example Lisanti (2017). The injection rate of CR electrons or positrons
is expressed as

$s(r,E)=\frac{\rho^{2}(r)}{2m_{\chi}^{2}}\frac{\mathrm{d}\langle\sigma
v\rangle}{\mathrm{d}E}\Bigg{|}_{\chi\chi\rightarrow e^{+}e^{-}},$ (1)

where $m_{\chi}$ is the WIMP mass, $E$ is the CR energy and $r$ is the
galactocentric radius. The injection rate also depends on the DM density
profile $\rho(r)$ and the annihilation cross-section $\langle\sigma v\rangle$.
We assume a Navarro–Frenk–White (NFW; Navarro et al., 1997) model for the
density profile given by

$\rho_{\text{NFW}}(r)=\frac{\rho_{s}}{\frac{r}{r_{s}}\left(1+\frac{r}{r_{s}}\right)^{2}}\
,$ (2)

where $\rho_{s}$ is the characteristic density and $r_{s}$ is a scale length.
The velocity- and spin-averaged cross-section for WIMP annihilation into
$e^{\pm}$-pairs per unit energy can be obtained as

$\frac{\mathrm{d}\langle\sigma
v\rangle}{\mathrm{d}E}=\sum_{f^{+}f^{-}}\text{BR}_{f^{+}f^{-}}\langle\sigma
v\rangle\frac{\mathrm{d}N_{f^{+}f^{-}\rightarrow e^{+}e^{-}+X}}{\mathrm{d}E}\
,$ (3)

where BR${}_{f^{+}f^{-}}$ is the branching ratio that describes the weighting
of the elements for any Standard Model particle pair $f^{+}f^{-}$, into which
the WIMPs will annihilate. This can be calculated with for example the Fortran
package package DarkSUSY (Bringmann et al., 2018), but we assume annihilation
into monochromatic electron-positron pairs. Then the derivative in Eq. (3)
simplifies to $\tfrac{\mathrm{d}N_{f^{+}f^{-}\rightarrow
e^{+}e^{-}+X}}{\mathrm{d}E}=\delta(E-m_{\chi})$. Once an electron or positron
is injected into the DM halo, it diffuses through the turbulent magnetic field
in the halo of the dSph while emitting synchrotron radiation.

###  2.2 Cosmic ray diffusion

The electron mass is always much lower than the WIMP mass, $m_{e}\ll
m_{\chi}$, so the $e^{\pm}$-pairs are ultra-relativistic. Since magnetic
fields in galaxies are mostly turbulent (Beck, 2015), we assume the
propagation of $e^{\pm}$-pairs to be dominated by diffusion (Colafrancesco et
al., 2007; Regis et al., 2014). Hence, the CR propagation is described with
the stationary diffusion–loss equation, which we adopt as our transport model.
Due to the spherical geometry of dSph galaxies, we further assume an isotropic
momentum distribution, so that the diffusion–loss equation depends only on
galactocentric radius $r$ and CR energy:

$D(E)\frac{1}{r}\frac{\partial^{2}}{\partial
r^{2}}\left[n_{e}(r,E)r\right]+\frac{\partial}{\partial
E}\left[b(E,B)n_{e}(r,E)\right]+s(r,E)=0\ .$ (4)

The parameter $b(E,B)$ describes the total energy loss-rate and $D(E)$
describes the diffusion coefficient. As a boundary condition, we assume that
the CR number density (per unit volume and energy) $n_{e}$ vanishes at
diffusion-halo radius $r_{h}$, as the diffusion coefficient rises to infinity
due to vanishing magnetic fields, so that $n_{e}(r_{h},E)=0$. We assume that
$r_{\rm h}$ is on the order of the half-light radius $r_{\star}$. These
assumptions, including stationarity, are well justified and became the de-
facto standard over the years (Colafrancesco et al., 2006; McDaniel et al.,
2017).

Contributions to the total energy loss include synchrotron radiation,
$b(E,B)_{\rm sync}$, and inverse Compton scattering, $b(E,B)_{\rm ICS}$,
losses. Other energy losses, such as bremsstrahlung and ionisation losses, are
suppressed due to the density of both ionized and neutral gas (Regis et al.,
2015). The total energy loss of ultra-relativistic CR $e^{\pm}$ in nearby
galaxies (at redshift $z=0$) is then

$b(E,B)=b_{\rm sync}(E,B)+b_{\rm ICS}(E,B)\\\ \approx 2.546\times
10^{-17}\left[1+\left(\frac{B}{3.24\leavevmode\nobreak\
\mu\text{G}}\right)^{2}\right]\left[\frac{E}{1\text{GeV}}\right]^{2}\text{GeV}\,\rm
s^{-1},$ (5)

where $B$ is the magnetic field and 3.24 $\mu$G is the field with the energy
density of the cosmic microwave background (CMB) photons at z=0. Inverse
Compton scattering on radiation fields other than the CMB can be neglected.
This is because the radiation energy density due to stellar light,
$u_{\star}$, is quite low in comparison to that of the CMB, $u_{\rm CMB}$,
with $u_{\star}/u_{\text{CMB}}\approx
0.3\,\%\left(L_{V}/10^{6}L_{\odot}\right)\left(r_{\star}/\text{kpc}\right)^{-2}$
, where $L_{V}$ is the $V$-band luminosity (Vollmann, 2021).

The average magnetic field strength in dSph galaxies is generally poorly
known. In order to observe any radio continuum emission, the CR electrons have
to be confined to the plasma, which implies that the magnetic field energy
density cannot be too low in comparison with the CR energy density, even if
the exact relation is not clear. For our typical sensitivities of 0.25-0.5
mJy/beam, frequency $\nu=144$ MHz, beam size 20”, spectral index
$\alpha=-0.8$, pathlength $l=400$ pc, and proton/electron ratio $K_{0}=0$, an
equipartition magnetic field strength of $\approx$2 $\mu$G can be calculated
for an $e^{+}e^{-}$-plasma (Beck & Krause, 2005). We adopt a more conservative
value of 1 $\mu$G here, but vary it by an order of magnitude to account for
the large uncertainties.

For the diffusion coefficient, we assume an energy-dependent power law

$D(E)=D_{0}\left(\frac{E}{1\>\text{GeV}}\right)^{\delta}\ ,$ (6)

where $D_{0}=D(1\,\text{GeV})$ and the power-law index $\delta$ describes the
energy dependence which is determined by the adopted turbulence model. The
diffusion coefficient and its energy dependency in dSph galaxies is prone to
uncertainties. We adopt a value of $10^{27}\leavevmode\nobreak\ \rm
cm^{2}\,s^{-1}$ for $D_{0}$ which is in agreement with observations of nearby
dwarf irregular systems (Murphy et al., 2012; Heesen et al., 2018). As no
measurements exist for diffusion in dSph galaxies, we vary $D_{0}$ within two
orders of magnitude. For the energy dependence, we assume Kolmogorov-like
turbulence resulting in $\delta=\tfrac{1}{3}$ (Kolmogorov et al., 1991). This
model is supported by observations of the Milky Way (Korsmeier & Cuoco, 2016).
A key question is to what extent there is turbulence in the magnetic field
structure at all? Observations of the Milky Way and other galaxies usually
fall into the two categories. Either the turbulence is extrinsically
generated, where turbulence cascades down from large scales to the small
scales relevant for CR scattering; or the turbulence is intrinsically
generated by the CR themselves via Plasma instabilities; the latter scenario
is usually referred to a self-confinement (Zweibel, 2013). Because in dSph
galaxies there is presumably no external source of turbulence such as
supernova remnants, the self-confinement scenario is probably more
appropriate. In star-forming galaxies this is likely the case as well for CRs
with energy of a few GeV, where observations indicate a weaker energy
dependency for CRs with less than 10 GeV energy (Heesen, 2021). This is hence
another source of uncertainty, at least for low WIMP masses.

###  2.3 Model scenarios

In order to deal with the uncertainties in the diffusion coefficient, $D_{0}$,
and the magnetic field strength, $B$, we employ three different model
scenarios. Our standard values
($D_{0}=10^{27}\>\text{cm}^{2}\>\text{s}^{-1},B=1\>\mu\text{G}$) define the
”intermediate” (INT) scenario. In the ”optimistic” (OPT) scenario, we choose
values that boost the DM signal. With a diffusion coefficient of
$10^{26}\>\text{cm}^{2}\>\text{s}^{-1}$ and an average magnetic field strength
of $10\>\mu\text{G}$, the CRs diffuse slowly and emit more synchrotron
radiation due to the strong magnetic field. We note that these values are
highly optimistic and probably unrealistic, but they serve us as a reference
point for the maximum signal we might possibly expect. Such a high magnetic
field strength is only observed in regions of concentrated star formation in
nearby dwarf irregular galaxies (Hindson et al., 2018). In the ”pessimistic”
(PES) scenario, we use a high diffusion coefficient of
$10^{29}\>\text{cm}^{2}\>\text{s}^{-1}$ in a comparably weak magnetic field of
$0.1\leavevmode\nobreak\ \mu\rm G$. In this situation, the CRs emit less
synchrotron radiation and the DM signal is lower. For such weak fields, most
of the CR electron energy, whether primary or secondary, would be lost via
inverse-Compton radiation.

This study does not cover the case in which the magnetic-field strengths are
even weaker than $\mathcal{O}(0.1\leavevmode\nobreak\ \mu$G). In this case a
ballistic description for the CR electron propagation would be more
appropriate. Since we can not exclude this possibility neither theoretically
nor observationally, its study is left for future work.

###  2.4 Diffusion regimes

Because the full solution of the diffusion–loss equation (Eq. 4) is rather
complicated, Vollmann (2021) defined three regimes, which allows one to
simplify the solution. These regimes depend on the ratio of the CR diffusion
time-scale to the energy loss time-scale. The diffusion time-scale is

$\tau_{\text{diff}}=\frac{r_{h}^{2}}{D(E)}\ ,$ (7)

where $r_{h}$ is again the radius of the diffusion halo. As already mentioned,
we assume that $r_{\rm h}$ is on the order of the half-light radius
$r_{\star}$.

The loss time-scale from synchrotron radiation and inverse Compton scattering
of CRs is

$\tau_{\text{loss}}=\frac{E}{b(E,B)}=1.245\left(1+\left(\frac{B}{3.24\>\mu\rm
G}\right)^{2}\right)^{-1}\left(\frac{E}{1\>\text{GeV}}\right)^{-1}\mathrm{Gyr}\
.$ (8)

When $\tau_{\text{diff}}\gg\tau_{\text{loss}}$, one can assume that the CR
lose all their energy so rapidly that diffusion can be neglected and the first
term in Eq. (4) vanishes. We refer to this assumption as ”regime A” or the no-
diffusion approximation. ”Regime B” is defined such that $\tau_{\rm
diff}\approx\tau_{\rm loss}$. In this regime the we have to consider the full
solution of Eq. (4). Vollmann (2021) shows that the solution can be expressed
as sum of Fourier-like modes as function of radius, where we consider only the
leading term in this expansion. For $\tau_{\text{diff}}\ll\tau_{\text{loss}}$,
one can neglect the second term of Eq. (4), as the CRs diffuse so rapidly that
they leave the dSph galaxy without losing energy. We refer to this assumption
as as ”regime C” or the rapid-diffusion approximation.

###  2.5 Synchrotron signal occurrence

The radio continuum intensity $I_{\nu}$ is the integral of the radio
emissivity $j_{\nu}(r)$ along the line of sight (LoS)

$I_{\nu}=\int_{\text{LoS}}dl\>j_{\nu}(r(l))\ .$ (9)

Following Vollmann (2021), the emissivity can be separated into a halo part
$H(r)$, a spectral part $X(\nu)$, and a normalizing pre-factor

$j_{\nu}(r)=\frac{\langle\sigma v\rangle}{8\pi m_{\chi}^{2}}\ H(r)\ X(\nu)\ .$
(10)

Both the halo and the spectral part depend on the diffusion regime. The
regime-specific equations for the halo-factors are presented in Vollmann
(2021). For the halo part, we employ the leading-mode approximation (regime
$B$). For the NFW profile (Eq. 2) the halo function is

$H_{B}(r)=h_{B}\frac{1}{r}\sin{\left(\frac{\pi r}{r_{h}}\right)},$ (11)

where $h_{B}$ is the halo factor in units of $\rm GeV^{2}\,cm^{-5}$ which
contains the part of $H_{B}$ that is independent of radius:

$h_{B}=2\left[\text{si}\left(\pi\right)-\frac{8r_{h}}{\pi
r_{s}}+\dots\right]\frac{\rho_{s}^{2}r_{s}^{2}}{r_{h}}\ \ \text{where}\ \
\text{si}\left(x\right)=\\!\int_{0}^{x}dt\>\tfrac{\sin\left(t\right)}{t}\ .$
(12)

Note that Eq. (11) is a simplified version only valid for emissivities. In
order to compare with our measured intensities, we implemented the actual halo
factor as calculated for intensities (appendix B in Vollmann, 2021).

For the spectral part, it is not viable to only consider one regime, as
$X(\nu)$ strongly depends on the environment, for example which energy loss
mechanism is dominant in the specific situation. Hence, we use all three
spectral parts for the regimes $A$, $B$, and $C$, respectively, from Vollmann
(2021):

$X_{A}(\nu)=\frac{2\sqrt{3}e^{3}B}{m_{e}}\int
dz\>\frac{F(z)}{z}\frac{E(\nu/z)}{b(E(\nu/z))}\int_{E(\nu/z)}^{\infty}dE\>S(E)\
,$ (13) $\displaystyle X_{B}(\nu)$ $\displaystyle=$
$\displaystyle\frac{2\sqrt{3}e^{3}B}{m_{e}}\int
dz\>\frac{F(z)}{z}\frac{E(\nu/z)}{b(E(\nu/z))}e^{-\eta(E(\nu/z))}$ (14)
$\displaystyle\int_{E(\nu/z)}^{\infty}dE\>S(E)e^{\eta(E)}\ ,$
$X_{C}(\nu)=\frac{2\sqrt{3}e^{3}B}{m_{e}}\frac{r_{h}^{2}}{\pi^{2}}\int
dz\>\frac{F(z)}{z}\frac{E(\nu/z)S(E(\nu/z))}{D(E(\nu/z))}\ ;$ (15)

with
$\eta(E)=\tfrac{\pi^{2}}{r_{h}^{2}}\int_{E}^{\infty}dE^{\prime}\>\tfrac{D(E^{\prime})}{b(E^{\prime})}$.
The function $F(x)$ is described by Ghisellini et al. (1988) as

$F(x)=x^{2}\left[K_{\frac{4}{3}}(x)K_{\frac{1}{3}}(x)-x\frac{3}{5}\left(K^{2}_{\frac{4}{3}}(x)-K^{2}_{\frac{1}{3}}(x)\right)\right]\
,$ (16)

with $K_{i}(x)$ as the modified Bessel functions of the second kind. For
monochromatic $e^{+}e^{-}$ injection, the CR energy is given by
$E(\nu)=\sqrt{2\pi m_{e}^{3}\nu/(3eB)}$ and the spectral injection function is
$S(E)=\delta(E-m_{\chi})$. All three formulas are related, where $X_{A}$ and
$X_{C}$ are the limits for $X_{B}$ when assuming $\eta\rightarrow 0$ and
$\eta\rightarrow\infty$, respectively.

Assuming that the total radio emissivity is due to DM annihilation, it is
straightforward to see that the shape of the radial profile is determined only
by the halo function $H(r)$ of Eq. (10). We can therefore express the
emissivity as:

$j_{\nu}(r)=N_{B}\cdot H_{B}(r)\ ,$ (17)

where $N_{B}$ is referred to as the signal-strength parameter that contains
all the terms that do not depend on the radius (Vollmann, 2021). Now we can
get an expression for the cross-section in terms of $N_{B}$

$\langle\sigma v\rangle=\frac{8\pi m_{\chi}^{2}N_{B}}{X_{j}(\nu)},$ (18)

where $j\in\\{A,B,C\\}$ specifies the diffusion regime. A similar approach is
used in Regis et al. (2014) and Vollmann et al. (2020). The factor $N_{B}$
connects predictions to observations, as it is proportional to the intensity
of the DM signal in the radio band.

##  3 Methodology

###  3.1 LoTSS observations

The data used for our analysis were observed as part of the LOFAR Two-metre
Sky Survey (LoTSS; Shimwell et al., 2017, 2019) and published in the second
data release (LoTSS-DR2; Shimwell et al., 2022). LoTSS is a deep low-frequency
survey with LOFAR HBA at 144 MHz with 24 MHz bandwidth. LoTSS-DR2 includes
observations of 841 pointings and covers 5634 square degrees of the northern
hemisphere. LoTSS data have a maximum angular resolution of 6”, referred to as
high-resolution data, and additional low-resolution data at 20” angular
resolution. The high-resolution data are important for the subtraction of
point-like sources, whereas the low-resolution data allow us to detect
extended emission at high signal-to-noise ratios. The rms noise at 20”
resolution is 50–100 $\mu\rm Jy\,beam^{-1}$ (Shimwell et al., 2022).

We analyze six dSph galaxies which are observed in the LoTSS-DR2. These
galaxies have half-light radii between 20 and 600 pc with distances between 30
and 218 kpc (see Table 1). These galaxies are Canes Venatici I (CVnI), Ursa
Major I (UMaI), Ursa Major II (UMaII), Ursa Minor (UMi), Willman I (WilI), and
Canes Venatici II (CVnII). These are the only non-disturbed dSph galaxies
observed with the LOFAR HBA at this point in time. Additional four galaxies
with declinations above +20∘ (Geringer-Sameth et al., 2015) could be observed
in the future. LOFAR sensitivity is greatly reduced at lower declinations
(Shimwell et al., 2017), but with longer observational times it is possible to
observe additional six galaxies with declinations between +10∘ and +20∘
(Geringer-Sameth et al., 2015; Ackermann et al., 2015).

Table 1: Properties of the galaxies in our sample from LoTSS-DR2.

dSph | R.A. (J2000) | Decl. (J2000) | $r_{\star}$ | $d$ | $r_{s}$ | $\rho_{s}$ | rms-noise | References
---|---|---|---|---|---|---|---|---
| $[\rm h\>m\>s]$ | $[^{\circ}\>{}^{\prime}\>"]$ | [pc] | [kpc] | [kpc] | [GeV cm-3] | [$\mu$Jy beam-1] |
CVnI | 13:28:03.5 | +33:33:21 | 564 | 218 | 2.27 | 0.5186 | 115 | 1
UMaI | 10:34:52.8 | +51:55:12 | 319 | 97 | 3.20 | 0.5473 | 74 | 1
UMaII | 08:51:30.0 | +63:07:48 | 149 | 32 | 4.28 | 2.794 | 60 | 1
UMi | 15:09:08.5 | +67:13:21 | 181 | 76 | 0.394 | 12.10 | 103 | 1
WilI | 10:49:23.0 | +51:01:20 | 21 | 38 | 0.173 | 15.18 | 83 | 2, 3, 4
CVnII | 12:57:10.0 | +34:19:15 | 74 | 160 | 8.04 | 1.331 | 68 | 1

(1) Geringer-Sameth et al. (2015); (2) Martin et al. (2008); (3) Willman et
al. (2011); (4) Sánchez-Conde et al. (2011);

###  3.2 Re-imaging the LoTSS data

We use re-calibrated LoTSS data, where the calibration is specially tailored
to our dSph galaxies (van Weeren et al., 2021). We re-imaged the $(u,v)$ data
with WSClean v2.9 (Offringa et al., 2014; Offringa & Smirnov, 2017). The
points in the $uv$-plane are weighted using Briggs robust weighting as a
compromise between uniform and natural weighting. A robustness parameter of
robust=$-$0.2 was found to produce the highest signal-to-noise ratio for the
extended emission on the scales that we are interested in. Further imaging
parameters are listed in Appendix B.

We exclude emission on large angular scales ($\gtrsim$$1^{\circ}$) which can
be attributed to the Milky Way (Erceg et al., 2022) by excluding $(u,v)$ data
at short baselines. We used lower limits to the $uv$-range between 60 and 400
$\lambda$, corresponding to angular scales from to 7′ to 46′, making sure that
these scales are not smaller than size of the galaxy. Compact sources are
subtracted from the $(u,v)$ data prior to imaging. This is done by first
producing a source catalog with the Python Blob Detector and Source Finder
(PyBDSF; Rafferty, D. and Mohan, N., 2019). Since not all background sources
are of point-like nature, we additionally use the ”À trous wavelet
decomposition module” integrated in PyBDSF. This module decomposes the
residual maps resulting from the internal subtraction of the fitted Gaussians,
into wavelets of different scales (see Holschneider et al., 1989). We used
between two and five wavelet scales depending on the galaxy. Sources were then
subtracted as Gaussians from the $(u,v)$ data using the Default Pre-Processing
Pipeline software (DPPP; van Diepen et al., 2018).

The maps were deconvolved with the multi-scale and auto-masking options to
remove any residuals comparable to the size of the galaxies. Maps were then
restored with a Gaussian beam at 20″ angular resolution. The re-imaging steps
(as well as further steps in the cross-section limits calculation) are
automated using Python.111https://github.com/FinnWelzmueller/wimpsSoftware

###  3.3 Calculation of the cross-section upper limits

To constrain the annihilation cross-section of WIMPs, we use central
amplitudes of the radial intensity profiles. To generate the radial profiles
we average intensities within annuli of increasing radius and with constant
width. For every galaxy, the width of an annulus is set to 20”, which is equal
to the full-width-at-half-maximum (FWHM) of the restoring Gaussian beam.

The expected shape of the radial profile is described by the halo function
(Eq. 11). To analyze the observed radial profiles we approximate the shape as
a fixed-width Gaussian with a FWHM equal to $r_{\star}$. Note that the halo
function depends on the square root of the density so the FWHM should indeed
be approximated with $r_{\star}$ instead of 2$r_{\star}$, which we would
expect from the density distribution. We vary the Gaussian central amplitude
to best-fit the observed radial profile. Since the size of the dSph diffusion
radius is mostly unknown, it is important to verify the non-detection of a DM-
related signals on various scales, not only the one we assumed earlier. Hence,
we vary the FWHM of the Gaussians using values that are higher and lower than
$r_{\star}$.

To mimic a DM halo, we inject fake sources directly into the point-source
subtracted $uv$-data. The fake source is constructed as a two-dimensional
Gaussian with the FWHM equal to the stellar radius of the galaxy. This is a
simplification as the real signal may be of a different shape. It has,
however, only a negligible effect on our inferred limits. The amplitude of the
Gaussian is varied until we get a 2$\sigma$ detection by fitting a Gaussian to
the radial profile as for the purely observational profile. Since the Gaussian
FWHM is fixed for each dSph, the only free parameter is the central amplitude
$a$ which is related to the factor $N_{B}$ in Eq. (18). The transformation
between the central amplitude of the Gaussian radial intensity distribution,
$a$, and the factor $N_{B}$ is done using the halo factor (Eq. 12):

$N_{B}=\frac{\pi\left(0.4\right)^{2}}{2h_{B}}\ a\ ,$ (19)

where the numerical factors account for the different shapes of the Gaussian
source fitted to the data and assumed form of the halo function described by
Eq. (11). Specifically, the width $w$ (equivalent to the standard deviation)
of the Gaussian is equal to $w=r_{\star}/(0.4d)$. Additionally, the spectral
function approximation is calculated for the appropriate scenario. The upper
limits on the cross-section are determined as a function of WIMP mass by
inserting $N_{B}$ into Eq. (18).

###  3.4 Stacking

In addition to looking at each galaxy separately we combined the data from all
galaxies through stacking. We use two different approaches for stacking the
data. The first was to generate the radial profile for each galaxy separately,
then to re-scale to the stellar radius, and then stack the profiles. The
second approach was to re-scale and stack the images, and only then generate
the radial profile from the stacked image.

####  3.4.1 Stacking radial profiles

Our first approach to stacking the data is to stack the radial intensity
profiles that are re-scaled to the stellar radius. We use the stellar radius
instead of the NFW scaling radius at it is the much more reliable observable.
We match the width of the annuli in which the radio intensity is averaged to
be equal to 0.05$r_{\star}$ for each galaxy. For most galaxies this is larger
than the beam FWHM but not for all which might cause a slight correlation
between adjacent data points. Each data point in the radial profile is
expressed in terms of $r_{\star}$ and the intervals between them are equal
even if the actual size on the sky is different. Note that the intensities do
not have to be corrected for distance. The re-scaled radial profiles are
combined by calculating the noise-weighted mean. By fitting a Gaussian to the
stacked profiles we confirmed they are consistent with zero.

After that we needed to handle the fake sources to be able to calculate the
limits on the cross-section. We first stacked the profiles with same fake
sources as for the individual galaxies. The significance of the detection was
higher than 2$\sigma$ so in the stacked profile we can detect a fainter DM
signal. To determine exactly how much fainter, we lowered the intensity of all
the injected sources by a common factor, so the flux density ratio between the
galaxies remains the same as in the individual analysis. This factor was
chosen to achieve a detection with a significance of 2$\sigma$ in the stacked
radial profile.

Once we have the stacked radial intensity profile for all galaxies, we repeat
the fitting to determine the combined value for the Gaussian amplitude $a$. To
calculate limits on the cross-section from the combined data of our galaxies,
we use this value for $a$ and average all the other terms in Eq. (18) for each
WIMP mass and each scenario. This is justified because by averaging the other
terms we calculate the average emission from galaxies and this is exactly what
we get when stacking.

####  3.4.2 Stacking images

Our second approach to stacking is to combine the galaxies in image space.
This is done in the following steps:

  1. 1.

A re-scaled cutout image of every galaxy is created with a size of
4$r_{\star}$ and 1367$\times$1367 pixel2. The dimension in pixels is the
median of all galaxies while ensuring there are at least 7 pixels per FWHM.
After the re-scaling process, a single pixel will sample a larger section of
the sky for galaxies with a larger angular diameter compared to those with a
small angular diameter, but this is a necessary compromise that we need to
make for the image stacking. After this procedure, diffuse sources with equal
flux density would appear identical, regardless of the distance and the
stellar radius of the galaxy.

  2. 2.

The flux density variance $\sigma$ is calculated inside an annulus with an
inner radius of $r_{\star}$ and an outer radius of $2r_{\star}$.

  3. 3.

All images are stacked using the weighted mean. The weight is adopted as the
inverse square variance (1/$\sigma^{2}$) of each image, so that galaxy images
with lower noise will contribute more to the stacked image. Cosmological
surface brightness dimming is negligible because the galaxies are in the Local
Group so we do not apply any weighting with redshift.

From the final stacked image we generated the radial intensity profile (using
the same algorithm as for the individual profiles) and confirmed a non-
detection. To calculate the limits on the cross-section, we followed the same
procedure as for the profile stacking. We adjusted the multiplication factors
of injected sources and obtained the Gaussian amplitude $a$. This combined
value for $a$ was used in the calculation for the cross-section limits. Other
necessary parameters that depend on galaxy properties were averaged.

##  4 Results

Our presentation of the results is split into three parts. We start with
individual limits on the WIMP annihilation cross-section (Sect. 4.1). In Sect.
4.2 we present the combined limits from the stacking algorithm. Finally, in
Sect. 4.3 we discuss the limitations of our results.

###  4.1 Individual limits

Of the individual galaxies, we present first results of CVnI, which was
already analyzed by Vollmann et al. (2020) using the same technique but with a
slightly different implementation of other software. This galaxy serves as a
benchmark to test our data processing algorithm. The radial intensity profiles
with and without fake source are shown in Fig. 1. The amplitude of the
Gaussian fitted to the observational data should be compatible with zero
within the uncertainties to verify the non-detection of DM-related signals.
Contrary, the amplitude of the Gaussian fitted to the data with the added fake
source should be detected at 2$\sigma$ significance. This is indeed the case,
as the best-fitting amplitude for the profile including the injected source is
$\left(36\pm 13\right)\mu{\rm Jy\,beam}^{-1}$ whereas without fake source it
is $9\pm 15$ $\mu\rm Jy\,beam^{-1}$ at a FWHM of 8.2 arcmin. In Table 2, we
summarize the fitting results for each galaxy with the corresponding profiles
presented in Appendix A.

Figure 1: Radial intensity profiles for Canes Venatici I. Data points show the mean intensity in $20\arcsec$ wide annuli with the error bars showing the standard deviation of the mean. Blue data points are for purely observational data, red data points are for the same data with an additional $20\leavevmode\nobreak\ \rm mJy$ fake source. Dashed lines show the best-fitting Gaussians with a fixed FWHM of $r_{\star}$ equivalent to $8\aas@@fstack{\prime}2$. Table 2: Gaussian amplitude $a_{\text{Gauss, obs}}$ for purely observed data and amplitude $a_{\text{Gauss, inj}}$ for data with an additional injected fake source. $S_{\rm FS}$ is the fake source flux density and FWHM is the width of the Gaussian, here assumed to be equivalent to $r_{\star}.$ dSph | $S_{\text{FS, $2\sigma$}}$ | $a_{\text{Gauss, obs}}$ | $a_{\text{Gauss, inj}}$ | FWHM
---|---|---|---|---
| [mJy] | [$\mu{\rm Jy\,beam}^{-1}$] | [$\mu{\rm Jy\,beam}^{-1}$] | [arcmin]
CVnI | 20 | $9\pm 15$ | $36\pm 13$ | 8.2
UMaI | 35 | $-6\pm 10$ | $24\pm 11$ | 11.3
UMaII | 35 | $-1.1\pm 4.8$ | $12.9\pm 4.8$ | 16.0
UMi | 33 | $9.1\pm 9.6$ | $28.9\pm 9.0$ | 8.3
WilI | 5 | $-55\pm 42$ | $17\pm 39$ | 1.8
CVnII | 1 | $16\pm 28$ | $50\pm 29$ | 1.6

Fluctuations that can not be described by Gaussian statistics can affect the
fit of galaxies with small half-light radii as the number of data points is
small and any fluctuation may not average out. This is in particular the case
for WilI ($r_{\star}=110\arcsec$) and CVnII ($r_{\star}=95\arcsec$). For WilI,
a negative Gaussian amplitude for the fit to the purely observational data was
found. This is the only galaxy with a ‘signal’, but because it is negative it
can be ruled out as DM-related. For CVnII, the Gaussian amplitude is
compatible with zero albeit with a large uncertainty. We tried to mitigate the
limitation due to the small number of data points by increasing the region in
which radial profiles were measured to $3r_{\star}$ and $4r_{\star}$ for WilI
and CVnII, respectively.

Figure 2: Best-fitting Gaussian amplitudes for the radial intensity profiles
of Canes Venatici I with Gaussians of varying FWHM. The blue solid line is for
purely observational data and the red solid line the same data with an
additional $20$ mJy fake source. Shaded areas indicate $1\sigma$ confidence
intervals.

Since the size of the DM halo is uncertain, we varied the FWHM of the Gaussian
fit to the radial intensity profiles. Here, we want to investigate the
possible systematic uncertainty of the assumption that the FWHM of the signal
produced by the DM halo is equal to the stellar radius $r_{\star}$ which
itself has an uncertainty of around 10 to 15% (Geringer-Sameth et al., 2015).
For CVnI, the corresponding results are shown in Fig. 2. The data are in
agreement with zero within the 1$\sigma$ confidence intervals for almost the
entire range of FWHM. On the other hand, the $2\sigma$ detection of the
injected source over a wide range of FWHM is also evident. Only at small FWHM,
statistical fluctuations start to suppress the significance of the detection.
Signals on that scale are most likely due to fluctuations in the map and not
related to DM.

The next step is to calculate the diffusion and energy-loss time-scales. This
identifies the diffusion regime which then leads to the set of equations
needed to estimate limits on the annihilation cross-section. For the INT
scenario, the diffusion time-scale for CVnI is $\approx$30 Myr, whereas the
energy-loss time-scale is $\approx$110 Myr. Since both time-scales are of the
same order of magnitude, we use diffusion regime $B$ with Eq. (14) to
calculate the limits on the WIMP annihilation cross-section. We summarize the
diffusion and energy-loss time-scales together with the resulting diffusion
regimes for the three model scenarios for our six dSph galaxies in Table 3. We
note that both time-scales depend on the CR energy. We used a benchmark-energy
of $E=10\leavevmode\nobreak\ \text{GeV}$; a different CR energy may change the
choice of diffusion regimes and hence slightly affect the limits.

Table 3: Time-scales of CR diffusion and loss for every model scenario with a benchmark CR energy of 10 GeV. The resulting diffusion regime is also noted. dSph | Model scenario | $\tau_{\text{diff}}$ | $\tau_{\text{loss}}$ | Regime
---|---|---|---|---
|  | [Myr] | [Myr] |
CVnI | OPT | 288 | 11.4 | $A$
| INT | 28.8 | 113 | $B$
| PES | 0.288 | 124 | $C$
UMaI | OPT | 92.1 | 11.4 | $A$
| INT | 9.21 | 113 | $C$
| PES | 0.0921 | 124 | $C$
UMaII | OPT | 20.1 | 11.4 | $B$
| INT | 2.01 | 113 | $C$
| PES | 0.0201 | 124 | $C$
UMi | OPT | 29.7 | 11.4 | $B$
| INT | 2.97 | 113 | $C$
| PES | 0.0297 | 124 | $C$
WilI | OPT | 0.399 | 11.4 | $C$
| INT | 0.0399 | 113 | $C$
| PES | 0.000399 | 124 | $C$
CVnII | OPT | 4.96 | 11.4 | $B$
| INT | 0.496 | 113 | $C$
| PES | 0.00496 | 124 | $C$

In Fig. 3, we present the upper limits to the WIMP annihilation cross-section
from each individual galaxy. For comparison, we additionally plot the lower
limit on the annihilation cross-section calculated from the thermal WIMP
freeze-out mechanism by Steigman et al. (2012).

(a)

(b)

(c)

(d)

(e)

(f)

Figure 3: Individual upper limits on the WIMP annihilation cross-section for
the OPT, INT, and PES model scenarios. Line color represents the assumed model
scenario; red for OPT, black for INT, and blue for PES. Line styles represent
the assumed diffusion regime; dashed lines for regime A, solid lines for
regime B and dotted lines for regime C. Panel (a) shows CVnI, (b) UMaI, (c)
UMaII, (d) UMi, (e) WilI, and (f) CVnII. The gray line represents the lower
limit from the thermal freeze-out (Steigman et al., 2012).

###  4.2 Stacked limits

The stacked radial intensity profiles both for profile and image stacking are
shown in Fig. 4. The best-fitting Gaussian amplitudes are $a_{\text{obs,
profiles}}=0.3\pm 4.3\ \mu{\rm Jy\,beam}^{-1}$ and $a_{\text{obs,
images}}=-1.8\pm 2.9\ \mu{\rm Jy\,beam}^{-1}$ for the purely observed data
using profile and image stacking, respectively. In both stacking strategies,
the amplitudes are consistent with zero so the stacking does not reveal any
additional signal. We varied the FWHM of fitted Gaussians such as before for
the individual galaxies. Again, the observed amplitudes are consistent with
zero as shown in Fig. 5.

(a)

(b)

Figure 4: Stacked radial intensity profiles. Panel (a) shows profile stacking
and panel (b) image stacking. Data points show the mean intensity in adjacent
annuli with the error bars showing the standard deviation of the mean. Blue
data points are for purely observational data, red data points are for the
same data with an additional fake source. Dashed lines show the best-fitting
Gaussians with a fixed FWHM of $r_{\star}$. The radius is expressed as an
apparent angle $\theta$ scaled to the apparent size of the stellar radius
$\theta_{\star}=r_{\star}/d$.

(a)

(b)

Figure 5: Best-fitting Gaussian amplitudes for the stacked radial intensity
profiles. Panel (a) shows profile stacking and panel (b) image stacking. Solid
blue lines show the purely observational data, whereas solid red lines show
the same data with an additional fake source. Shaded areas indicate $1\sigma$
uncertainties. The FWHM of the Gaussian is expressed by the apparent size of
the stellar radius $\theta_{\star}=r_{\star}/d$.

The flux density of the injected fake sources from the individual galaxies is
divided by a factor to determine how much fainter is the signal that can be
detected by stacking. The factor is chosen to achieve a 2$\sigma$ detection,
it equals 1.8 for profile stacking and 2.8 for image stacking. We again fit
Gaussians to the resulting stacked radial profile and the Gaussian amplitudes
are $a_{\text{inj, profiles}}=9.4\pm 4.5\ \mu{\rm Jy\,beam}^{-1}$ for profile
stacking and $a_{\text{inj, images}}=6.6\pm 2.8\ \mu{\rm Jy\,beam}^{-1}$ for
image stacking, showing that we have detected the fake source at $2\sigma$
confidence.

Upper limits for the WIMP annihilation cross-section are calculated from the
best-fitting amplitudes of the radial profiles obtained with either method. A
unique value of amplitude $a$ is used and the other terms in Eq. (18) are
given as an average value for all galaxies. This average is calculated using
the same approximation regime as for the individual galaxies (Table 3). The
resulting upper limits for the cross-section from both approaches are shown in
Fig. 6, together with the average of the limits obtained for the individual
galaxies.

Figure 6: Upper limits on the WIMP annihilation cross-section from stacking.
Dotted lines are the averaged individual limits, dashed lines are limits
obtained by stacking radial profiles and solid lines are limits obtained by
stacking galaxy images. Colors of red, black, and blue correspond to the OPT,
INT, and PES model scenarios, respectively. The thermal freeze-out cross-
section (Steigman et al., 2012, grey) and upper limits from Fermi–LAT
$\gamma$-ray observations (Hoof et al., 2020, green) are shown for comparison.

###  4.3 Systematic uncertainties

There are several sources of inaccuracies in the obtained limits. The most
significant are the assumptions on the values of the average magnetic field
strength and diffusion coefficient. For this reason we used our three
different model scenarios in order to illustrate the influence of these
parameters. The differences between the cross-section limits in each scenario
are indeed two to three orders of magnitude showing the large uncertainty
resulting from the inaccuracy of these input parameters.

For our INT scenario we assume that the the diffusion coefficient is $10^{27}$
cm2s-1 as measured both in dwarf irregular galaxies (Murphy et al., 2012) and
in the Milky Way (Korsmeier & Cuoco, 2016). However, it is also possible that
the true value may be closer to $10^{28}$ cm2s-1 as found in late-type spiral
galaxies (Heesen et al., 2019) or even $10^{29}$ cm2s-1 such as for the halos
in edge-on spiral galaxies. On the other hand, the true value can also go down
to the value of $10^{26}$ cm2s-1 measured in other dwarf irregular galaxies
(Heesen et al., 2018). The magnetic field is much more uncertain but values
range between 0.1 $\mu$G and 10 $\mu$G assuming that the magnetic energy
density is in equipartition with the CR energy density within two orders of
magnitude. We note that this is a only a heuristic argument though as there
are no stringent physical reasons why there should be equipartition.
Generally, the diffusion coefficient depends on the magnetic field (Sigl,
2017), but to simplify our model we treat them independently and this is an
additional source of uncertainty.

A further uncertainty is introduced by the adopted NFW profile parameters
$r_{s}$ and $\rho_{s}$. While including full posterior probability-
distribution functions (PDF) of existing Bayesian fits in the literature
(e.g., Ando et al., 2020) would be more appropriate, it would also be
impractical. This is because of the overwhelming magnetic-field uncertainties
which would not be reduced by considering such PDFs. Therefore, we content
ourselves by using the best-fit values for $r_{s}$ and $\rho_{s}$ as derived
in (Geringer-Sameth et al., 2015), and direct our attention to the CR
propagation parameters (e.g., magnetic field) when discussing uncertainties.

Another source of uncertainty comes from the approximations made in the
spectral function where we choose one of the diffusion regimes ($A$, $B$, or
$C$) and use the appropriate equation (Eq. 13, 14, or 15). The regime is
chosen by comparing the diffusion and energy-loss time scales (Table 3). The
diffusion time-scale depends on the diffusion coefficient, $D_{0}$, and the
energy-loss time-scale depends on the average magnetic field, $B$. For this
reason, regimes were chosen independently for each model scenario. An
inappropriate choice of the approximation regime for certain ratio of time-
scales affects the cross-section limits much less then the the previously
mentioned uncertainty from the model parameters ($D_{0}$, $B$). There is
ongoing work to solve the transport equations regardless of the ratio of time-
scales (Vollmann et al., in prep.) which would completely eliminate such
uncertainties.

There is some evidence for tidal disruption or non-equilibrium kinematics in
WilI (Ibata et al., 1997; Willman et al., 2011; Geringer-Sameth et al., 2015)
which could cause the NFW profile parameters to be biased to higher masses.
Since that bias might artificially improve our stacked limits, we repeated the
stacked analysis without WilI but with the same injected source intensity for
other galaxies (see Appendix C). We found that the resulting limits without
WilI are stronger compared to the limits including all the galaxies. In the
end we decided to keep the weaker limits derived from all galaxies as the more
conservative estimate.

We also mention other sources of uncertainty which are small compared to the
model scenario uncertainty. First, we assumed that the density distribution in
dSph is described by a NFW density profile. This is just one of the
possibilities and alternative profiles (e.g., Einasto & Haud, 1989; Burkert,
1995) could also be used (Vollmann, 2021). Second, we assume a power-law
dependence of the diffusion coefficient on the CR energy (Eq. 6) and
Kolmogorov turbulence. If we are able to better constrain the magnetic field
and diffusion coefficient in dSphs, the mentioned uncertainties would become
more important.

##  5 Discussion and conclusions

We first compare the limits for the INT model scenario from individual
galaxies in Fig. 7. The individual limits within our sample vary due to the
different noise levels of the LoTSS maps and the different intrinsic
properties, such as DM density profile and distance. The best constraints on
the annihilation cross-section are derived from UMaII. This is also true for
$\gamma$-ray WIMP searches since this galaxy has the highest $J$-factor (Hoof
et al., 2020). The next best constraints are derived from CVnII. This is
counter-intuitive, since UMaII has the largest apparent size
($r_{\star}\approx 16\arcmin$) whereas CVnII is the smallest
($r_{\star}\approx 1\aas@@fstack{\prime}6$). But reviewing the halo-factor
defined in Eq. (11), it becomes evident that the apparent size is not the most
important but rather the characteristic density $\rho_{s}$, the scale length
$r_{s}$, and the half-light radius $r_{\star}$. According to Eq. (11), the
combination of these parameters is high for the galaxies UMaII and CVnII,
which explains the strong limits for these two galaxies. The same explanation
justifies why the limits set by WilI and CVnI are less stringent than average.
WilI has the lowest scale length within our set of galaxies222The spectral
factor $X(\nu)$ for WilI is two orders of magnitudes smaller than average
which justifies the moderate limit for this galaxy as well., CVnI has the
smallest characteristic density. Since the halo-factor depends on the square
of both values, this factor is smaller than average and hence the limits are
less stringent.

Figure 7: Individual upper limits on the WIMP annihilation cross-section for
the dSph galaxies in our sample for the INT model scenario.

We improved the limits by stacking the data for different galaxies. To test
the improvement we compared the stacked limits to the average value of all
individual limits (Fig. 6). The stacking of profiles lowered the limits
approximately two times and stacking of images three times. This is expected
because with stacking we are effectively extending the observing time. The
image stacking yields better results than the stacking of profiles but the
control of uncertainties is more difficult. We compare only the two specific
stacking strategies but in general there are many alternative ways of
combining the results of different galaxies. For example in a future study, a
more statistically rigorous approach could be to treat galaxies independently
with their own nuisance parameters while the DM model parameters are fitted
simultaneously (see e.g., Hoof et al., 2020).

The comparison of our stacking limits with other attempts to constrain the
WIMP annihilation cross-section, such as the thermal freeze-out cross-section
by Steigman et al. (2012) or dSph observations by Fermi–LAT (Hoof et al.,
2020), shows the competitiveness of our results. Our cross-section limits from
the INT scenario already exclude thermal WIMPs with masses below 20 GeV. In
the OPT scenario we can even exclude thermal WIMPs with masses below 70 GeV.
While the PES limits do not exclude thermal WIMPs in any mass range, they
still prove the validity of our concept (Vollmann et al., 2020; Regis et al.,
2014).

As customary, we assumed a smooth DM distribution. However, there has been
recent work stating that DM halos contain prompt cusps (Delos & White, 2022)
that would boost the annihilation signal, which, in turn, would lower our
limits by up to two orders of magnitude. The same effect would apply to limits
inferred from $\gamma$-ray excess in the Milky Way and other galaxies.

To date, the results from $\gamma$-ray observations of dSph galaxies presented
by Hoof et al. (2020) are among the strongest limits on the WIMP annihilation
cross-section. In the mass range below 20 GeV our INT limits are stronger than
those from Fermi–LAT. These limits are valid under the assumption of the
average magnetic field and diffusion coefficient in dSph galaxies; and
therefore have a large uncertainty. However, our results show that radio
observations of dSph can potentially constrain the WIMP annihilation cross-
section if one accepts this premise. This method is especially powerful for
WIMP masses below 10 GeV.

Considering the observational resources needed for both attempts, our method
is much more efficient. We only used LoTSS-DR2 survey data (Shimwell et al.,
2022, about 50 h in total for the six dSph galaxies) and have not performed
targeted observations. For comparison, Ackermann et al. (2015) used 6 years of
Fermi–LAT data. Compared to these observations, we achieved better limits at
lower masses in the OPT scenario and comparable limits in the INT scenario.
Hence, the advantages of our study are the sensitivity in the low-mass regime
and the efficiency in terms of observation time. The biggest drawback of using
radio continuum observations are of course the uncertainties related to the
strength of the magnetic field and the value of the diffusion coefficient. The
field strength of dSph galaxies could be measured from a grid of Faraday
rotation measures of polarized background sources with a sensitive radio
telescope, such as the Square Kilometre Array (SKA) (Johnston-Hollitt et al.,
2015).

In addition to the HBA used for LoTSS, LOFAR Low Band Antenna (LBA)
observations (de Gasperin et al., 2021) would improve the limits on the lower
mass end because the critical frequency of synchrotron radiation depends on
the square of the WIMP mass (Rybicki & Lightman, 1986). The SKA has also been
suggested as a promising future instrument for DM searches (Colafrancesco et
al., 2015).

###### Acknowledgements.

We thank the anonymous referee for a constructive report which helped to
improve the paper. This work is funded by the Deutsche Forschungsgemeinschaft
(DFG, German Research Foundation) under Germany’s Excellence Strategy –
EXC-2121 ”Quantum Universe” – 390833306. LOFAR (van Haarlem et al., 2013) is
the Low Frequency Array designed and constructed by ASTRON. It has observing,
data processing, and data storage facilities in several countries, which are
owned by various parties (each with their own funding sources), and that are
collectively operated by the ILT foundation under a joint scientific policy.
The ILT resources have benefited from the following recent major funding
sources: CNRS-INSU, Observatoire de Paris and Université d’Orléans, France;
BMBF, MIWF-NRW, MPG, Germany; Science Foundation Ireland (SFI), Department of
Business, Enterprise and Innovation (DBEI), Ireland; NWO, The Netherlands; The
Science and Technology Facilities Council, UK; Ministry of Science and Higher
Education, Poland; The Istituto Nazionale di Astrofisica (INAF), Italy. This
research made use of the Dutch national e-infrastructure with support of the
SURF Cooperative (e-infra 180169) and the LOFAR e-infra group. The Jülich
LOFAR Long Term Archive and the German LOFAR network are both coordinated and
operated by the Jülich Supercomputing Centre (JSC), and computing resources on
the supercomputer JUWELS at JSC were provided by the Gauss Centre for
Supercomputing e.V. (grant CHTB00) through the John von Neumann Institute for
Computing (NIC). This research made use of the University of Hertfordshire
high-performance computing facility and the LOFAR-UK computing facility
located at the University of Hertfordshire and supported by STFC
[ST/P000096/1], and of the Italian LOFAR IT computing infrastructure supported
and operated by INAF, and by the Physics Department of Turin university (under
an agreement with Consorzio Interuniversitario per la Fisica Spaziale) at the
C3S Supercomputing Centre, Italy. DJB acknowledges funding from the German
Science Foundation DFG, via the Collaborative Research Center SFB1491 ”Cosmic
Interacting Matters - From Source to Signal”.

## References

  * Abdelhameed et al. (2019) Abdelhameed, A. H., Angloher, G., Bauer, P., et al. 2019, Phys. Rev. D, 100, 102002
  * Ackermann et al. (2015) Ackermann, M., Albert, A., Anderson, B., et al. 2015, Phys. Rev. Lett., 115, 231301
  * Aharonian et al. (2008) Aharonian, F., Akhperjanian, A., Bazer-Bachi, A., et al. 2008, Astroparticle Physics, 29, 55
  * Alcock et al. (2000) Alcock, C., Allsman, R. A., Alves, D. R., et al. 2000, The Astrophysical Journal, 542, 281–307
  * Ando et al. (2020) Ando, S., Geringer-Sameth, A., Hiroshima, N., et al. 2020, Physical Review D, 102
  * Aprile et al. (2017) Aprile, E., Aalbers, J., Agostini, F., et al. 2017, Physical Review Letters, 119
  * Beck (2015) Beck, R. 2015, A&A Rev., 24, 4
  * Beck & Krause (2005) Beck, R. & Krause, M. 2005, Astronomische Nachrichten, 326, 414
  * Bringmann et al. (2021) Bringmann, T., Depta, P. F., Hufnagel, M., Ruderman, J. T., & Schmidt-Hoberg, K. 2021, Physical Review Letters, 127
  * Bringmann et al. (2018) Bringmann, T., Edsjö, T., Gondolo, P., Ullio, P., & Bergström, L. 2018, Journal of Cosmology and Astroparticle Physics, 2018, 033–033
  * Burkert (1995) Burkert, A. 1995, ApJ, 447, L25
  * Cheung et al. (2011) Cheung, C., Elor, G., Hall, L. J., & Kumar, P. 2011, Journal of High Energy Physics, 2011, 42
  * Colafrancesco et al. (2006) Colafrancesco, S., Profumo, S., & Ullio, P. 2006, A&A, 455, 21
  * Colafrancesco et al. (2007) Colafrancesco, S., Profumo, S., & Ullio, P. 2007, Phys. Rev. D, 75, 023513
  * Colafrancesco et al. (2015) Colafrancesco, S., Regis, M., Marchegiani, P., et al. 2015, in Advancing Astrophysics with the Square Kilometre Array (AASKA14), 100
  * Cook et al. (2020) Cook, R. H. W., Seymour, N., Spekkens, K., et al. 2020, Monthly Notices of the Royal Astronomical Society, 494, 135
  * de Gasperin et al. (2021) de Gasperin, F., Williams, W. L., Best, P., et al. 2021, A&A, 648, A104
  * Delos & White (2022) Delos, M. S. & White, S. D. M. 2022, Prompt cusps and the dark matter annihilation signal
  * Einasto & Haud (1989) Einasto, J. & Haud, U. 1989, A&A, 223, 89
  * Erceg et al. (2022) Erceg, A., Jelić , V., Haverkorn, M., et al. 2022, A&A, 663, A7
  * Finkbeiner & Weiner (2007) Finkbeiner, D. P. & Weiner, N. 2007, Phys. Rev. D, 76, 083519
  * Geringer-Sameth et al. (2015) Geringer-Sameth, A., Koushiappas, S. M., & Walker, M. 2015, The Astrophysical Journal, 801, 74
  * Ghisellini et al. (1988) Ghisellini, G., Guilbert, P. W., & Svensson, R. 1988, ApJ, 334, L5
  * Hall et al. (2010) Hall, L. J., Jedamzik, K., March-Russell, J., & West, S. M. 2010, Journal of High Energy Physics, 2010, 80
  * Hawking (1971) Hawking, S. 1971, Monthly Notices of the Royal Astronomical Society, 152, 75
  * Heesen (2021) Heesen, V. 2021, Ap&SS, 366, 117
  * Heesen & Brüggen (2021) Heesen, V. & Brüggen, M. 2021, Phys. Rev. D, 103, 048301
  * Heesen et al. (2019) Heesen, V., Buie II, E., Huff, C. J., et al. 2019, A&A, 622, A8
  * Heesen et al. (2018) Heesen, V., Rafferty, D. A., Horneffer, A., et al. 2018, Monthly Notices of the Royal Astronomical Society, 476, 1756
  * Hindson et al. (2018) Hindson, L., Kitchener, G., Brinks, E., et al. 2018, ApJS, 234, 29
  * Holschneider et al. (1989) Holschneider, M., Kronland-Martinet, R., Morlet, J., & Tchamitchian, P. 1989, Proceedings of the International Conference
  * Hoof et al. (2020) Hoof, S., Geringer-Sameth, A., & Trotta, R. 2020, Journal of Cosmology and Astroparticle Physics, 2020, 012
  * Hryczuk & Laletin (2021) Hryczuk, A. & Laletin, M. 2021, Journal of High Energy Physics, 2021, 26
  * Hütten & Kerszberg (2022) Hütten, M. & Kerszberg, D. 2022, TeV Dark Matter searches in the extragalactic gamma-ray sky
  * Ibarra (2015) Ibarra, A. 2015, Nuclear and Particle Physics Proceedings, 267-269, 323
  * Ibata et al. (1997) Ibata, R. A., Wyse, R. F. G., Gilmore, G., Irwin, M. J., & Suntzeff, N. B. 1997, The Astronomical Journal, 113, 634
  * Jaeckel & Ringwald (2010) Jaeckel, J. & Ringwald, A. 2010, Annual Review of Nuclear and Particle Science, 60, 405
  * Johnston-Hollitt et al. (2015) Johnston-Hollitt, M., Govoni, F., Beck, R., et al. 2015, in Advancing Astrophysics with the Square Kilometre Array (AASKA14), 92
  * Jungman et al. (1996) Jungman, G., Kamionkowski, M., & Griest, K. 1996, Physics Reports, 267, 195–373
  * Kim (1987) Kim, J. E. 1987, Phys. Rep, 150, 1
  * Kolmogorov et al. (1991) Kolmogorov, A. N., Levin, V., Hunt, J. C. R., Phillips, O. M., & Williams, D. 1991, Proceedings of the Royal Society of London. Series A: Mathematical and Physical Sciences, 434, 9
  * Korsmeier & Cuoco (2016) Korsmeier, M. & Cuoco, A. 2016, Phys. Rev. D, 94, 123019
  * Kramer et al. (2021) Kramer, E. D., Kuflik, E., Levi, N., Outmezguine, N. J., & Ruderman, J. T. 2021, Phys. Rev. Lett., 126, 081802
  * Lisanti (2017) Lisanti, M. 2017, in Theoretical Advanced Study Institute in Elementary Particle Physics: New Frontiers in Fields and Strings, 399–446
  * Martin et al. (2008) Martin, N. F., de Jong, J. T. A., & Rix, H.-W. 2008, The Astrophysical Journal, 684, 1075
  * McDaniel et al. (2017) McDaniel, A., Jeltema, T., Profumo, S., & Storm, E. 2017, Journal of Cosmology and Astroparticle Physics, 2017, 027
  * Milgrom (1983) Milgrom, M. 1983, ApJ, 270, 365
  * Murphy et al. (2012) Murphy, E. J., Porter, T. A., Moskalenko, I. V., Helou, G., & Strong, A. W. 2012, The Astrophysical Journal, 750, 126
  * Navarro et al. (1997) Navarro, J. F., Frenk, C. S., & White, S. D. M. 1997, The Astrophysical Journal, 490, 493–508
  * Offringa et al. (2014) Offringa, A. R., McKinley, B., Hurley-Walker, N., et al. 2014, MNRAS, 444, 606
  * Offringa & Smirnov (2017) Offringa, A. R. & Smirnov, O. 2017, MNRAS, 471, 301
  * Peccei & Quinn (1977a) Peccei, R. D. & Quinn, H. R. 1977a, Phys. Rev. D, 16, 1791
  * Peccei & Quinn (1977b) Peccei, R. D. & Quinn, H. R. 1977b, Phys. Rev. Lett., 38, 1440
  * Rafferty, D. and Mohan, N. (2019) Rafferty, D. and Mohan, N. 2019, PyBDSF 1.9.1 Documentation, https://www.astron.nl/citt/pybdsf/, Last accessed on 09.06.2021
  * Regis et al. (2014) Regis, M., Colafrancesco, S., Profumo, S., et al. 2014, Journal of Cosmology and Astroparticle Physics, 2014, 016–016
  * Regis et al. (2021) Regis, M., Reynoso-Cordova, J., Filipović, M. D., et al. 2021, Journal of Cosmology and Astroparticle Physics, 2021, 046
  * Regis et al. (2015) Regis, M., Richter, L., Colafrancesco, S., et al. 2015, Monthly Notices of the Royal Astronomical Society, 448, 3747–3765
  * Rybicki & Lightman (1986) Rybicki, G. B. & Lightman, A. P. 1986, Radiative Processes in Astrophysics, 2nd edn. (London: John Wiley & Sons, Ltd), 1–50
  * Sánchez-Conde et al. (2011) Sánchez-Conde, M. A., Cannoni, M., Zandanel, F., Gómez, M. E., & Prada, F. 2011, Journal of Cosmology and Astroparticle Physics, 2011, 011
  * Sanglard et al. (2005) Sanglard, V., Benoit, A., Bergé, L., et al. 2005, Phys. Rev. D, 71, 122002
  * Shimwell et al. (2022) Shimwell, T. W., Hardcastle, M. J., Tasse, C., et al. 2022, A&A, 659, A1
  * Shimwell et al. (2017) Shimwell, T. W., Röttgering, H. J. A., Best, P. N., et al. 2017, A&A, 598, A104
  * Shimwell et al. (2019) Shimwell, T. W., Tasse, C., Hardcastle, M. J., et al. 2019, A&A, 622, A1
  * Sigl (2017) Sigl, G. 2017, Astroparticle Physics: Theory and Phenomenology, 1st edn. (Atlantis Studies in Astroparticle Physics and Cosmology)
  * Steigman et al. (2012) Steigman, G., Dasgupta, B., & Beacom, J. F. 2012, Physical Review D, 86
  * van Diepen et al. (2018) van Diepen, G., Dijkema, T., & Offringa, A. 2018, DPPP: Default Pre-Processing Pipeline
  * van Haarlem et al. (2013) van Haarlem, M. P., Wise, M. W., Gunst, A. W., et al. 2013, A&A, 556, A2
  * van Weeren et al. (2021) van Weeren, R. J., Shimwell, T. W., Botteon, A., et al. 2021, A&A, 651, A115
  * Vollmann (2021) Vollmann, M. 2021, Journal of Cosmology and Astroparticle Physics, 2021, 068
  * Vollmann et al. (2020) Vollmann, M., Heesen, V., W. Shimwell, T., et al. 2020, Monthly Notices of the Royal Astronomical Society, 496, 2663
  * Vollmann et al. (in prep.) Vollmann, M., Welzmüller, F., & Gajović, L. in prep.
  * Weinberg (1978) Weinberg, S. 1978, Phys. Rev. Lett., 40, 223
  * Wilczek (1978) Wilczek, F. 1978, Phys. Rev. Lett., 40, 279
  * Willman et al. (2011) Willman, B., Geha, M., Strader, J., et al. 2011, The Astronomical Journal, 142, 128
  * Zweibel (2013) Zweibel, E. G. 2013, Physics of Plasmas, 20, 055501

##  Appendix A Results for the individual galaxies

###  A.1 Radial intensity profiles

Figure 8: Radial intensity profiles from each of the six dSphs. Blue color
indicates purely observational data, red color indicates data with an
additional flux density from a fake source. The individual flux densities are
listed in Table 2. The dashed lines show the best-fitting Gaussians with a
FWHM of $r_{\star}$. The half-light radii are listed in Table 1.

(a) Canes Venatici I

(b) Ursa Major I

(c) Ursa Major II

(d) Ursa Minor

(e) Willman I

(f) Canes Venatici II

###  A.2 Fitting plots

Figure 9:  Best-fitting Gaussian amplitudes for the radial intensity profiles
in Figs. 8(a)\- 8(f). Blue color indicates purely observational data, red
color indicates data with an additional fake source. The added individual flux
densities are listed in Table 2. The shaded areas are $1\sigma$ intervals.

(a) Canes Venatici I

(b) Ursa Major I

(c) Ursa Major II

(d) Ursa Minor

(e) Willman I

(f) Canes Venatici II

##  Appendix B Re-imaging parameters

Table 4: Parameters used to re-image each of the 6 galaxies. dSph | Lower | Robustness | Wavelet | Baseline-av.
---|---|---|---|---
|  $uv$-cut [$\lambda$] |  | scale | factor
CVnI | 160 | -0.2 | 564 | 8.5221
UMaI | 60 | -0.2 | 319 | 5.6113
UMaII | 60 | -0.2 | 149 | 7.0936
UMi | 160 | -0.2 | 181 | 7.0025
WilI | 400 | -0.2 | 21 | 5.2477
CVnII | 400 | -0.2 | 74 | 10.3342

##  Appendix C Stacked limits excluding the galaxy Willman I

Figure 10: Upper limits on the WIMP annihilation cross-section from stacking
without the galaxy Willman I (in orange) compared to the results of stacking
all galaxies (in standard colors form Fig. 6).