© The Authors 2023

1. 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⁻²⁶ cm³ s⁻¹ 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 γ-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) 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 ≲-pagination100 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 electron–positron 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, , 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, 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.

2. Dark matter calculations

2.1. 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

(1)

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

(2)

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

(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 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.

2.2. 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:

(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_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 ultra-relativistic CR e^± in nearby galaxies (at redshift z = 0) is then

(5)

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⁶ L_⊙)(r_⋆/kpc)⁻², 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 ν = 144 MHz, beam size 20″, spectral index α = −0.8, path length l = 400 pc, and proton/electron ratio K₀ = 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,

(6)

where D₀ = 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²⁷ cm² s⁻¹ for D₀, 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₀ within two orders of magnitude. For the energy dependence, we assumed Kolmogorov-like turbulence resulting in (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.

2.3. Model scenarios

In order to deal with the uncertainties in the diffusion coefficient, D₀, and the magnetic field strength, B, we employed three different model scenarios. Our standard values (D₀ = 10²⁷ cm² s⁻¹, 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²⁶ cm² s⁻¹ 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²⁹ cm² s⁻¹ 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 magnetic-field strengths are even weaker than 𝒪(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.

2.4. 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

(7)

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

(8)

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.

2.5. Synchrotron signal occurrence

The radio continuum intensity, I_ν, is the integral of the radio emissivity, j_ν(r), along the line of sight (LoS):

(9)

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

(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

(11)

where h_B is the halo factor in units of GeV² cm⁻⁵, which contains the part of H_B that is independent of radius:

(12)

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):

(13)

(14)

(15)

with . The function F(x) is described by Ghisellini et al. (1988) as

(16)

where K_i(x) is the modified Bessel functions of the second kind. For monochromatic e^± injection, the CR energy is given by 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

(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 get an expression for the cross section in terms of N_B:

(18)

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.

3. Methodology

3.1. LoTSS observations

The data used for our analysis are 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 were important for the subtraction of point-like 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⁻¹ (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 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).

3.2. 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.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 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 & Mohan 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 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 reimaging steps (as well as further steps in the cross-section limits calculation) were automated using PYTHON¹.

3.3. 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 increasing 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)):

(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_⋆/(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).

3.4. 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.

3.4.1. 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.

3.4.2. 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². 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/σ²) 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.

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 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⁻¹ whereas without fake source it is (9 ± 15) μJy beam⁻¹ 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.

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 .

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–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.

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.

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).

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).

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_{obs, profiles} = (0.3 ± 4.3) μJy beam⁻¹ and a_{obs, images} = (−1.8 ± 2.9) μJy beam⁻¹ 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.

Fig. 4. Stacked radial intensity profiles. Panel a: 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.

Fig. 5. Best-fitting Gaussian amplitudes for the stacked radial intensity profiles. Panel a: 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σ uncertainties. The FWHM of the Gaussian is expressed as the apparent size of the stellar radius, θ⋆ = r⋆/d.

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⁻¹ for profile stacking and a_{inj, images} = (6.6 ± 2.8) μJy beam⁻¹ 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.

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.

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 assumed that the diffusion coefficient is 10²⁷ cm² s⁻¹ 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²⁸ cm² s⁻¹ as found in late-type spiral galaxies (Heesen et al. 2019) or even 10²⁹ cm² s⁻¹ 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²⁶ cm² s⁻¹ 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 (Eqs. (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₀, and the energy-loss timescale 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 timescales affects the cross-section limits much less than the previously mentioned uncertainty from the model parameters (D₀, 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.

5. 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 (). 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 galaxies², 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.

Fig. 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 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).

Acknowledgments

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

Appendix A: Results for the individual galaxies

A.1. Radial intensity profiles

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 in Table 2. The dashed lines show the best-fitting Gaussians with a FWHM of r⋆. The half-light radii are listed in Table 1.

A.2. Fitting plots

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 in Table 2. The shaded areas are 1σ intervals.

Appendix B: Reimaging parameters

Appendix C: Stacked limits excluding the galaxy Willman I

Fig. C.1. 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 (colors have the same meaning as in Fig. 6).

All Tables

All Figures

	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 .
In the text

	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.
In the text

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).

In the text

Fig. 4. Stacked radial intensity profiles. Panel a: 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.

In the text

	Fig. 5. Best-fitting Gaussian amplitudes for the stacked radial intensity profiles. Panel a: 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σ uncertainties. The FWHM of the Gaussian is expressed as the apparent size of the stellar radius, θ⋆ = r⋆/d.
In the text

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.

In the text

	Fig. 7. Individual upper limits on the WIMP annihilation cross section for the dSph galaxies in our sample for the INT model scenario.
In the text

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 in Table 2. The dashed lines show the best-fitting Gaussians with a FWHM of r⋆. The half-light radii are listed in Table 1.

In the text

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 in Table 2. The shaded areas are 1σ intervals.

In the text

	Fig. C.1. 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 (colors have the same meaning as in Fig. 6).
In the text