Weakly interacting massive particle cross section limits from LOFAR observations of dwarf spheroidal galaxies

Context. Weakly interacting massive particles (WIMPs) can self-annihilate, thus providing us with a way to indirectly detect 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 µ G) magnetic fields in dSph galaxies exist, the particles produced in DM annihilation emit synchrotron radiation in the radio band. Aims. We used the non-detection of 150 MHz radio continuum emission from dSph galaxies with the LOw Frequency ARray (LO-FAR) to derive constraints on the annihilation cross section of WIMPs in electron–positron pairs. Our main underlying assumption is that the transport of the cosmic rays can be described by the di ff usion approximation, which necessitates the existence of magnetic fields. Methods. We used observations of six dSph galaxies in the LOFAR Two-metre Sky Survey (LoTSS). The data were reimaged, and a radial profile was generated for each galaxy. We also used stacking to increase the sensitivity. In order to derive upper limits on the WIMP cross section, we injected fake Gaussian sources into the data, which were then detected with 2 σ 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 explored the uncertainty due to the choice of di ff usion and magnetic field parameters by constructing three di ff erent 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 ( ⟨ σ v ⟩≈ 2 . 2 × 10 − 26 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 γ -ray observations of multiple dwarf galaxies, and they are especially strong for low WIMP masses.


Introduction
Dark matter (DM) is known to interact gravitationally, and only weakly via other fundamental forces of nature, which makes it difficult to observe.Among the most promising candidates for DM are weakly interacting massive particles (WIMPs; Jungman et al. 1996), (quantum chromodynamics) axions (Peccei & Quinn 1977a,b;Weinberg 1978;Wilczek 1978) or axion-like particles (Kim 1987;Jaeckel & Ringwald 2010), massive compact halo objects (Alcock et al. 2000), sterile neutrinos (Ibarra 2015), primordial black holes (Hawking 1971), and modification of the Newtonian dynamics 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 ⟨σv⟩ ≈ 2.2 × 10 −26 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 freezein, 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), andothers. nian et al. 2008) are able to provide new stringent upper limits for the WIMP annihilation cross section that depend on the particle mass.Recent work by Delos & White (2022) shows that if DM is thermally produced WIMPs, thousands of Earth-mass prompt cusps should be present in every solar mass of DM.This is manifested in a drastic increase in the expected DM annihilation signal compared to a smooth DM distribution, substantially tightening observational cross-section limits.The most recent γ-ray results (Hoof et al. 2020), assuming a smooth DM distribution, exclude the thermal relic annihilation cross section for WIMPs with masses of ≲ 100 GeV annihilating through the quark and τ-lepton channels; and for WIMPs with masses of ≲ 20 GeV annihilating through the electron-positron channel.
Radio continuum observations provide a complementary approach to constraining the WIMP annihilation cross section.Such searches exploit the fact that highly energetic electronpositron pairs (e ± 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 e ± pairs.For WIMP masses on the order of 10 GeV and microgauss (µG) magnetic fields, the critical frequency of synchrotron radiation, ν = 3eBE 2 /(2πm 3 e 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 that traces 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(Regis et al. , 2021)).
Dwarf spheroidal galaxies are very promising systems to look for DM-related emission since 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 dSph galaxies 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 are therefore close to Earth.This allows us to resolve the spectra of their stellar light well, and the estimates of their DM content are relatively precise (Hütten & Kerszberg 2022).
The LOw Frequency ARray (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 the high angular resolution needed to identify compact background sources with the 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 was already performed by Vollmann et al. (2020).That proof-of-concept study shows that, for an individual galaxy, the limits are comparable to or even better than that from Fermi-LAT, assuming reasonable values for the magnetic field strength and diffusion coefficient.This work expands upon that study by considering six galaxies observed with LOFAR HBA.In order to improve the WIMP annihilation cross-section limits, we stacked the signal from these six galaxies in various ways.We thus present improved upper limits on the annihilation of low-mass WIMPs into e ± pairs in the GeV-TeV mass range using new radio continuum data.
This work is organized 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 via both profile and image stacking.Section 4 contains our results.We finish off in Sect.5, which contains the discussion and conclusions.

Cosmic ray injection
The injection of standard model particles by WIMP annihilation is described by for example Lisanti (2017).The injection rate of CR electrons or positrons is expressed as where m χ 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 ρ(r) and the annihilation cross section ⟨σv⟩.
We assume a Navarro-Frenk-White (NFW; Navarro et al. 1997) model for the density profile given by where ρ s is the characteristic density and r s is a scale length.The velocity-and spin-averaged cross section for WIMP annihilation into e ± pairs per unit energy are obtained as 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 is 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 Once an electron or positron is injected into the DM halo, it diffuses through the turbulent magnetic field in the halo of the dSph galaxy while emitting synchrotron radiation.

Cosmic ray diffusion
The electron mass is always much lower than the WIMP mass, m e ≪ m χ , so the e ± pairs are ultra-relativistic.Since magnetic fields in galaxies are mostly turbulent (Beck 2015), we assume the propagation of e ± 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, such that the diffusion-loss equation depends only on galactocentric radius, r, and CR energy: 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 h is on the order of the half-light radius r ⋆ .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) sync , and inverse Compton scattering, b(E, B) ICS , losses.Other energy losses, such as bremsstrahlung and ionization losses, are suppressed due to the low density of both ionized and neutral gas (Regis et al. 2015).The total energy loss of ultrarelativistic CR e ± in nearby galaxies (at redshift z = 0) is then where B is the magnetic field and 3.24 µ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 ⋆ , is quite low in comparison to that of the CMB, u CMB , with u ⋆ /u CMB ≈ 0.3 % L V /10 6 L ⊙ r ⋆ /kpc −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 −1 , frequency ν = 144 MHz, beam size 20", spectral index α = −0.8,path length l = 400 pc, and proton/electron ratio K 0 = 0, an equipartition magnetic field strength of ≈2 µG is calculated for an e ± -plasma (Beck & Krause 2005).We adopted a more conservative value of 1 µG here, but varied it by an order of magnitude to account for the large uncertainties.
For the diffusion coefficient, we assume an energy-dependent power law, where D 0 = D(1 GeV) and the power-law index, δ, describes the energy dependence, which is determined by the adopted turbulence model.The diffusion coefficient and its energy dependence in dSph galaxies are prone to uncertainties.We adopted a value of 10 27 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 varied D 0 within two orders of magnitude.
For the energy dependence, we assumed Kolmogorov-like turbulence resulting in δ = 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, if it is present 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 CRs 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 dependence for CRs with less than 10 GeV energy (Heesen 2021).This is hence another source of uncertainty, at least for low WIMP masses.

Model scenarios
In order to deal with the uncertainties in the diffusion coefficient, D 0 , and the magnetic field strength, B, we employed three different model scenarios.Our standard values (D 0 = 10 27 cm 2 s −1 , B = 1 µG) define the "intermediate" (INT) scenario.In the "optimistic" (OPT) scenario, we chose values that boost the DM signal.With a diffusion coefficient of 10 26 cm 2 s −1 and an average magnetic field strength of 10 µ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 used a high diffusion coefficient of 10 29 cm 2 s −1 in a comparably weak magnetic field of 0.1 µ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 magneticfield strengths are even weaker than O(0.1 µG).In this case, a ballistic description for the CR electron propagation would be more appropriate.Since we cannot exclude this possibility theoretically nor observationally, its study is left for future work.

Diffusion regimes
Because the full solution of the diffusion-loss equation (Eq.4) is rather complicated, Vollmann (2021) defines three regimes, which allow one to simplify the solution.These regimes depend on the ratio of the CR diffusion timescale to the energy loss timescale.The diffusion timescale is where r h is again the radius of the diffusion halo.As already mentioned, we assume that r h is on the order of the half-light radius, r ⋆ .The loss timescale from synchrotron radiation and inverse Compton scattering of CRs is When τ diff ≫ τ loss , one can assume that the CRs 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 τ diff ≈ τ 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 the sum of Fourier-like modes as a function of radius, where we consider only the leading term in this expansion.For τ diff ≪ τ 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.

Synchrotron signal occurrence
The radio continuum intensity, I ν , is the integral of the radio emissivity, j ν (r), along the line of sight (LoS): Following Vollmann ( 2021), the emissivity can be separated into a halo part, H(r), a spectral part, X(ν), and a normalizing prefactor, 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 where h B is the halo factor in units of GeV 2 cm −5 , which contains the part of H B that is independent of radius: It should be noted 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 (see Vollmann 2021, Appendix B).
For the spectral part, it is not viable to only consider one regime, as X(ν) strongly depends on the environment, for example which energy loss mechanism is dominant in the specific situation.Hence, we used all three spectral parts for the regimes A, B, and C, respectively, from Vollmann (2021): The function F(x) is described by Ghisellini et al. (1988) as where K i (x) is the modified Bessel functions of the second kind.
For monochromatic e ± injection, the CR energy is given by E(ν) = 2πm 3 e ν/(3eB) and the spectral injection function is S (E) = δ(E − m χ ).All three spectral part formulae are related, where X A and X C are the limits for X B when assuming η → 0 and η → ∞, 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 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 get an expression for the cross section in terms of N B : where j ∈ {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.

LoTSS observations
The data used for our analysis are observed as part of the LOFAR Two-metre Sky Survey (LoTSS; Shimwell et al. 2017Shimwell et al. , 2019) ) and published in the second data release (LoTSS-DR2; Shimwell et al. 2022).LoTSS is a deep low-frequency survey with LO-FAR 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 were important for the subtraction of pointlike sources, whereas the low-resolution data allowed us to detect extended emission at high signal-to-noise ratios.The rms noise at 20" resolution is 50-100 µJy beam −1 (Shimwell et al. 2022).We analyzed six dSph galaxies that 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 nondisturbed dSph galaxies observed with the LOFAR HBA at this point in time.Additional four galaxies with declinations above +20

Reimaging the LoTSS data
We use recalibrated LoTSS data, where the calibration was specially tailored to our dSph galaxies (van Weeren et al. 2021).
We reimaged the (u, v) data with WSClean v2.9 (Offringa et al. 2014;Offringa & Smirnov 2017).The points in the (u, v)-plane were weighted using Briggs robust weighting as a compromise between uniform and natural weighting.A robustness parameter of robust=−0.2was found to produce the highest signalto-noise ratio for the extended emission on the scales that we are interested in.Further imaging parameters are listed in Appendix B.
We excluded emission on large angular scales (≳1 • ) that 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 (u, v) range between 60 and 400 λ, corresponding to angular scales from to 7 ′ to 46 ′ , making sure that these scales are not smaller than the size of the galaxy.Compact sources were subtracted from the (u, v) data prior to imaging.This was 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 have a point-like nature, we additionally used 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 automasking 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 reimaging steps (as well as further steps in the cross-section limits calculation) were automated using Python.1

Calculation of the cross-section upper limits
To constrain the annihilation cross section of WIMPs, we used central amplitudes of the radial intensity profiles.To generate the radial profiles, we averaged intensities within annuli of increas-ing radius and constant width.For every galaxy, the width of an annulus was 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 ⋆ .We note that the halo function depends on the density squared so the FWHM should indeed be approximated with r ⋆ instead of 2r ⋆ , which we would expect from the density distribution.We varied the Gaussian central amplitude to best fit the observed radial profile.Since the size of the dSph galaxy 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 varied the FWHM of the Gaussians using values that are higher and lower than r ⋆ .
To mimic a DM halo, we injected fake sources directly into the point-source subtracted (u, v) data.The fake source was 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.The simplification, however, has only a negligible effect on our inferred limits.The amplitude of the Gaussian was varied until we got a 2σ detection by fitting a Gaussian to the radial profile, as for the purely observational profile.Since the Gaussian FWHM was fixed for each dSph galaxy, the only free parameter was 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 was done using the halo factor (Eq. 12): 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 ⋆ /(0.4d).Additionally, the spectral function approximation was calculated for the appropriate scenario.The upper limits on the cross section were determined as a function of WIMP mass by inserting N B into Eq.(18).

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

Stacking radial profiles
Our first approach to stacking the data was to stack the radial intensity profiles that are rescaled to the stellar radius.We used the stellar radius instead of the NFW scaling radius as it is the much more reliable observable.We set the width of the annuli in which the radio intensity is averaged to 0.05r ⋆ for each galaxy.This was larger than the beam FWHM for most, but not all, galaxies, which might cause a slight correlation between adjacent data points.Each data point in the radial profile was expressed in terms of r ⋆ and the intervals between them were equal even if the actual size on the sky is different.We note that the intensities did not have to be corrected for distance.The rescaled radial profiles were 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σ so in the stacked profile we could 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 remained the same as in the individual analysis.This factor was chosen to achieve a detection with a significance of 2σ in the stacked radial profile.
Once we had the stacked radial intensity profile for all galaxies, we repeated 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 used this value for a and averaged all the other terms in Eq. ( 18) for each WIMP mass and each scenario.This is justified because by averaging the other terms we calculated the average emission from galaxies and this is exactly what we got when stacking.

Stacking images
Our second approach to stacking was to combine the galaxies in image space.This was done in the following steps: 1.A rescaled cutout image of every galaxy was created with a size of 4r ⋆ and 1367×1367 pixel 2 .The dimension in pixels was the median of all galaxies while ensuring there were at least seven pixels per beam FWHM.After the rescaling process, a single pixel sampled a larger section of the sky for galaxies with a larger angular diameter, compared to those with a small angular diameter, but this was 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. The flux density variance σ was calculated inside an annulus with an inner radius of r ⋆ and an outer radius of 2r ⋆ .3. All images were stacked using the weighted mean.The weight was adopted as the inverse square variance (1/σ 2 ) of each image, so that galaxy images with lower noise contributed more to the stacked image.Cosmological surface brightness dimming is negligible because the galaxies are in the Local Group so we have not applied 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.

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.

Individual limits
Of the individual galaxies, we present first results of CVnI, which is 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σ significance.This is indeed the case, as the best-fitting amplitude for the profile including the injected source is (36 ± 13) µJy beam −1 whereas without fake source it is (9 ± 15) µ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.  Fluctuations that cannot 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 ⋆ = 110 ′′ ) and CVnII (r ⋆ = 95 ′′ ).For WilI, a negative Gaussian amplitude for the fit to the purely observational data is 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 ⋆ and 4r ⋆ for WilI and CVnII, respectively.
Since the size of the DM halo is uncertain, we varied the FWHM of the Gaussian fit to the radial intensity profiles.Here, we wanted 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 ⋆ , 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σ confidence intervals for almost the entire range of FWHM.On the other hand, the 2σ 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 was to calculate the diffusion and energy-loss timescales.This identified the diffusion regime that then lead to the set of equations needed to estimate limits on the annihilation cross section.For the INT scenario, the diffusion timescale for CVnI is ≈30 Myr, whereas the energy-loss timescale is ≈110 Myr.Since both timescales are on the same order of magnitude, we used diffusion regime B with Eq. ( 14) to calculate the limits on the WIMP annihilation cross section.We summarize the diffusion and energy-loss timescales together with the resulting diffusion regimes for the three model scenarios for our six dSph galaxies in Table 3.We note that both timescales depend on the CR energy.We used a benchmark-energy of E = 10 GeV; a different CR energy may change the choice of diffusion regimes and hence slightly affect the limits.
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).

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 obs, profiles = (0.3 ± 4.3) µJy beam −1 and a obs, images = (−1.8±2.9)µ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 as before for the individual galaxies.Again, the observed amplitudes are consistent with zero as shown in Fig. 5.The flux density of the injected fake sources from the individual galaxies was divided by a factor to determine how much fainter is the signal that can be detected by stacking.The factor was chosen to achieve a 2σ detection, it equals 1.8 for profile stacking and 2.8 for image stacking.We again fitted Gaussians to the resulting stacked radial profile and the Gaussian amplitudes are a inj, profiles = (9.4± 4.5) µJy beam −1 for profile stacking and a inj, images = (6.6 ± 2.8) µJy beam −1 for image stacking, showing that we have detected the fake source at 2σ confidence.
Upper limits for the WIMP annihilation cross section were calculated from the best-fitting amplitudes of the radial profiles obtained with either method.A unique value of amplitude a was used and the other terms in Eq. ( 18) were given as an average value for all galaxies.This average was 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.

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 assumed that the diffusion coefficient is 10 27 cm 2 s −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 cm 2 s −1 as found in late-type spiral galaxies (Heesen et al. 2019) or even 10 29 cm 2 s −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 cm 2 s −1 measured in dwarf irregular galaxies (Heesen et al. 2018).The magnetic field is much more uncertain but values range between 0.1 µG and 10 µ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 treated them independently and this is an additional source of uncertainty.
A further uncertainty is introduced by the adopted NFW profile parameters r s and ρ s .While including full posterior probability-distribution functions (PDFs) 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 contend ourselves by using the best-fit values for r s and ρ 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 chose one of the diffusion regimes (A, B, or C) and used the appropriate equation (Eq.13, 14, or 15).The regime was chosen by comparing the diffusion and energy-loss timescales (Table 3).The diffusion timescale depends on the diffusion coefficient, D 0 , and the energy-loss timescale depends on the average magnetic field, Fig. 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.Red, black, and blue correspond to the OPT, INT, and PES model scenarios, respectively.The thermal freeze-out cross section (Steigman et al. 2012, gray) and upper limits from Fermi-LAT γ-ray observations (Hoof et al. 2020, green) are shown for comparison.
B. For this reason, regimes were chosen independently for each model scenario.An inappropriate choice of the approximation regime for certain ratio of timescales affects the cross-section limits much less than 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 timescales (Vollmann et al. in prep.), which would completely eliminate such uncertainties.
There is some evidence for tidal disruption or nonequilibrium 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 find 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 that are small compared to the model scenario uncertainty.First, we assumed that the density distribution in dSph galaxies 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 assumed 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 dSph galaxies, the mentioned uncertainties would become more important.

Discussion and conclusions
We first compared 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 the DM density profile and distance.The best constraints on the annihilation cross section are derived from UMaII.This is also true for γ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 counterintuitive since UMaII has the largest apparent size (r ⋆ ≈ 16 ′ ), whereas CVnII is the smallest (r ⋆ ≈ 1 ′ .6).But reviewing the halo factor defined in Eq. ( 11), it becomes evident that the apparent size is not the most important parameter; the most important are instead the characteristic density, ρ s , the scale length, r s , and the half-light radius, r ⋆ .According to Eq. ( 11), the combination of these parameters is high for 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 galaxies2 , and 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.
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 by a factor of approximately two and the stacking of images by three.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 controlling the uncertainties is more difficult.We only compared 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 galaxy 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 is 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 γ-ray excess in the Milky Way and other galaxies.
To date, the results from γ-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.Our 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 galaxies 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 six years of Fermi-LAT data.Compared to these observations, we achieve 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 drawbacks 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 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).

Appendix A: Results for the individual galaxies
Fig. 1: Radial intensity profiles for Canes Venatici I. Data points show the mean intensity in 20 ′′ wide annuli, with the error bars showing the standard deviation of the mean.Blue data points are for purely observational data, and red data points are for the same data with an additional 20 mJy fake source.Dashed lines show the best-fitting Gaussians with a fixed FWHM of r ⋆ equivalent to 8 ′ .2.

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

Fig. 3 :
Fig. 3: Individual upper limits on the WIMP annihilation cross section for the OPT, INT, and PES model scenarios.Line colors represent 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).

Fig. 4 :Fig. 5 :
Fig. 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, and 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 ⋆ .The radius is expressed as an apparent angle, θ, scaled to the apparent size of the stellar radius, θ ⋆ = r ⋆ /d.(a)

Fig. 7 :
Fig.7: Individual upper limits on the WIMP annihilation cross section for the dSph galaxies in our sample for the INT model scenario.
Fig. A.1: Radial intensity profiles from each of the six dSph galaxies: Canes Venatici I (panel a), Ursa Major I (b), Ursa Major II (c), Ursa Minor (d), Willman I (e), and Canes Venatici II (f).Blue indicates purely observational data, and red indicates data with an additional flux density from a fake source.The individual flux densities are listed inTable 2. The dashed lines show the best-fitting Gaussians with a FWHM of r ⋆ .The half-light radii are listed in Table 1.
Fig. A.2: Best-fitting Gaussian amplitudes for the radial intensity profiles in Figs.A.1a-A.1f: Canes Venatici I (panel a), Ursa Major I (b), Ursa Major II (c), Ursa Minor (d), Willman I (e), and Canes Venatici II (f).Blue indicates purely observational data, and red indicates data with an additional fake source.The added individual flux densities are listed inTable 2. The shaded areas are 1σ intervals.

Table 1 :
• (Geringer-Sameth et al. 2015)could be observed in the future.LOFAR sensitivity is greatly reduced at lower declinations Properties of the galaxies in our sample from LoTSS-DR2.

Table 2 :
Gaussian amplitude a Gauss, obs for purely observed data and amplitude a Gauss, inj for data with an additional injected fake source.S FS is the fake source flux density, and FWHM is the width of the Gaussian, here assumed to be equivalent to r ⋆ .

Table 3 :
Timescales of CR diffusion and loss for every model scenario with a benchmark CR energy of 10 GeV.The resulting diffusion regime is also noted.

Table 2 .
The dashed lines show the best-fitting Gaussians with a FWHM of r ⋆ .The half-light radii are listed in Table1.

Table 2 .
The shaded areas are 1σ intervals.