Constraining Annihilating Dark Matter Using the Multifrequency Radio Flux Profiles of the M33 Galaxy

Radio data can give stringent constraints for annihilating dark matter. In general, radio observations can detect very accurate radio flux density with high resolution and different frequencies for nearby galaxies. We are able to obtain the radio flux density as a function of distance from the galactic center and frequencies S(r, ν). In this article, we demonstrate a comprehensive radio analysis of the M33 galaxy, combining the radio flux density profile S(r) and the frequency spectrum S(ν) to get the constraints of dark matter annihilation parameters. By analyzing the archival radio data obtained from the Effelsberg telescope, we show that the dark matter annihilation contributing to the radio flux density might be insignificant in the disk region of the M33 galaxy. Moreover, by including the baryonic radio contribution, we constrain the 2σ conservative upper limits of the annihilation cross section, which can be complementary to the existing constraints based on neutrino, cosmic-ray, and gamma-ray observations. Our results indicate that analyzing the galactic multifrequency radio flux profiles can give useful and authentic constraints on dark matter for the leptophilic annihilation channels.


Introduction
Recent studies using radio data of galaxies and galaxy clusters can constrain the dark matter annihilation parameters, such as the lower limits of dark matter mass and the upper limits of annihilation cross section for different annihilation channels (Colafrancesco, Profumo & Ullio 2006;Egorov & Pierpaoli 2013;Storm et al. 2013;Chan 2016;Chan et al. 2019;Regis et al. 2021;Beck & Sarkis 2023;Lavis et al. 2023).If dark matter particles can self-annihilate, they could give high-energy electrons, positrons, photons, and neutrinos.In particular, the high-energy electrons and positrons can produce synchrotron radiation in radio bands when they are moving inside the magnetic field of a galaxy or a galaxy cluster.The underlying physics is well-known and this indirect method of dark matter detection has been applied for almost two decades (Colafrancesco, Profumo & Ullio 2006;Chan et al. 2019;Beck & Sarkis 2023;Lavis et al. 2023).
Comparing with using gamma-ray data to constrain dark matter (Calore et al. 2015;Daylan et al. 2016), radio detection can provide very accurate radio flux density maps which indicate the radio flux density as a function of radius r (i.e. the distance from galactic center) and frequencies ν.Since the resolution of gamma-ray detection is relatively poor, it is very difficult to obtain the gamma-ray flux as a function of r for other galaxies or galaxy clusters.For radio analyses, especially in galaxies, some previous studies have used the radio spectrum S(ν) (Colafrancesco, Profumo & Ullio 2006;Chan & Lee 2020, 2022) or the radio flux density profile S(r) (Chan et al. 2021) to constrain dark matter parameters.These methods can provide different perspectives in handling the radio contribution of dark matter annihilation in galaxies.
In this article, by using the archival radio data of the M33 galaxy (Buczilowski & Beck 1987;Buczilowski 1988;Tabatabaei, Krause & Beck 2007), we re-construct the radio maps of the M33 galaxy for different observing frequencies.Then, we can obtain the radio flux density as a function of both radius r and frequency ν: S(r, ν).Based on this multi-variable function, we have developed a comprehensive radio analysis by combining the frequency spectral data (flux with different frequencies) and radio flux density profiles (flux at different r) to constrain the dark matter parameters.Some recent studies have applied this similar idea to analyze the multi-frequency radio flux density profiles to constrain annihilating dark matter in the Large Magellanic Cloud (Regis et al. 2021) and galaxy clusters (Beck & Sarkis 2023).This study provides the first analysis of the M33 galaxy by using this comprehensive method.

Radio flux profile construction
We have used the archival radio data of the M33 galaxy obtained by the Effelsberg radio telescope in Buczilowski & Beck (1987); Buczilowski (1988); Tabatabaei, Krause & Beck (2007).The data consist of four different central frequencies (1.42 GHz,2.70 GHz,4.85 GHz and 8.35 GHz) at different positions of the M33 galaxy.After performing the data reduction, we can plot the radio flux density maps (in mJy/beam) for different frequencies (see Fig. 1 for the map with ν = 4.85 GHz).
However, the radio flux density maps might consist of some strong foreground or background point sources which should not be included.Therefore, we subtract the bright compact point sources from the radio flux density maps to get the reduced radio flux maps.
Here, since we do not have enough point source information for the 1.42 GHz radio flux density map, we cannot get the reduced radio flux map for ν = 1.42 GHz.We will only analyse the other three radio flux density maps (2.70 GHz,4.85 GHz,and 8.35 GHz) in the followings.
Using a constant bin of r = 60 arcsec, we take the azimuthal averaging of the radio flux density in concentric bins for different frequencies from the reduced radio flux maps.The fluctuations in the radio flux density for the same bin would contribute to the 1σ error bar of the data.Following this method, we can get the radio flux density profiles S(r) for ν = 2.70 GHz, 4.85 GHz and 8.35 GHz (see Fig. 2).These three radio flux density profiles can be combined to form the multi-variable function S(r, ν).
M33 has a bulge component and a disk component.However, the radio flux density profiles only include one data point for the bulge component.To minimize the free parameters in our analysis, we will only analyze the disk region (r > 96"), but not the data points for the bulge component.For the disk region, we have total 19 data points for each of the observing frequencies in the radio flux density profiles.

The theoretical framework
A large amount of high-energy electrons and positrons would be produced via dark matter annihilation.The injection energy spectrum of these electrons and positrons dN e,inj /dE for different annihilation channels can be predicted by numerical simulations (Cirelli et al. 2011).These high-energy electrons and positrons would emit synchrotron radiation in radio bands when there is a large magnetic field strength.The magnetic field strength in M33 is constrained to be B = 8.1 ± 0.5 µG (Berkhuijsen, Beck & Tabatabaei 2013), which is relatively large among the galaxies in our Local Group.The cooling and the diffusion of the  electrons and positrons produced from dark matter annihilation is given by the diffusioncooling equation (Ginzburg & Syrovatskii 1964;Atoyan et al. 1995): where dn e /dE, D(E) and Q(E, r) are the electron/positron density spectrum, the diffusion function, and the particle-injection source, respectively.The diffusion function can be written as D(E) = D 0 (E/1 GeV) δ .We adopt the diffusion coefficient D 0 = 2.0 × 10 28 cm 2 s −1 obtained in Berkhuijsen, Beck & Tabatabaei (2013) and use the benchmark value of δ = 1/3 (Kolmogorov 1941) in our present study to model the diffusion process.In most galaxies, the cooling process is dominated by the synchrotron emission and inverse Compton scattering (ICS).The cooling function (in unit of 10 −16 GeV s −1 ) can be expressed as Here, E, B are in the units of GeV and µG, respectively, and U rad is the energy density of the Interstellar Radiation Field (ISRF) in the unit of eV cm −3 .For the ICS in M33, by adding the total radiation density ranging from infra-red to ultraviolet 0.77 eV cm −3 (Thirlwall et al. 2020) to the cosmic microwave background energy density 0.26 eV cm −3 , we get U rad ≈ 1.03 eV cm −3 .
The particle-injection source term in Eq. ( 1) is given by (Vollmann 2021) where ρ DM (r) is the dark matter density in spherically symmetric profile and σv is the annihilation cross section.By setting ∂ ∂t ( dne dE ) = 0 in Eq. ( 1) and satisfying the boundary condition dne(r h ,E) dE = 0 with diffusion halo radius r h , the general solution of the equilibrium electron density spectrum is obtained in terms of the Fourier-series representation of the Green's function as follows (Vollmann 2021): where the dimensionless variable η(E) is given by In our analysis, we take the halo size as the size of the M33 galaxy r h = 5.04 kpc.
The average power at frequency ν under magnetic field B for synchrotron emission induced by the dark matter annihilation is given by where ν g = eB/(2πm e c), r e is the classical electron radius, and F syn (x/ sin θ) = x/ sin θ ∞ x/ sin θ K 5/3 (y)dy with the quantity x defined as in which γ is the Lorentz factor of the electron/positrons, and ν p = 8890[n(r)/1 cm −3 ] 1/2 Hz is the plasma frequency with the number density of the thermal electrons (n(r) ≈ 1 cm −3 ).
The radio flux density emitted by a galaxy within a solid angle ∆Ω due to dark matter annihilation as observed from Earth is: Here, s is the line-of-sight distance and the factor 2 indicates the contributions of both high-energy electrons and positrons in the synchrotron radiation emission.
For a typical dark matter annihilation model, m DM and σv are the only free parameters.Nevertheless, standard cosmology predicts that σv = 2.2×10 −26 cm 3 /s if dark matter particles were produced thermally at the early epoch (Steigman, Dasgupta & Beacom 2012).This standard scenario will be particularly examined.
Apart from the dark matter contribution, the baryonic contribution is also very important for the radio flux density.We assume that the radio flux density contributed by the baryonic component is proportional to the baryonic density.As mentioned above, we will not analyze the bulge component as the data points are too few.We will only focus on the disk region.Observations show that the disk surface brightness profile can be best-fitted with an exponential function in −r (Seigar 2011).This can be transformed to the similar functional form of the line-of-sight baryonic contribution: where k(ν) is a constant depending on frequency ν only, and h = 1.70 ± 0.12 kpc is the disk scale length (Seigar 2011).The total radio contribution would be given by

Data analysis
Let's consider the dark matter-only model first (i.e.there is no baryonic contribution to the radio signal).For a particular frequency ν j , the goodness of fits can be determined by the reduced χ 2 value: where S obs (r i , ν j ) and σ obs (r i , ν j ) are the data of the observed flux density and uncertainties of the observed flux density at different angular radii r i respectively, N = 19 is the number of data points for each ν j , and p is the number of free parameters in the fitting.We can also combine the reduced χ 2 values for different frequencies to get the total χ 2 value: First of all, we examine the standard cosmological scenario (i.e.σv = 2.2 × 10 −26 cm 3 /s) for both dark matter density profiles.In Fig. 3, we show the χ 2 red,j and the total χ 2 as a function of dark matter mass for m DM ≥ 5 GeV.We can see that there exist some minimum χ 2 values which indicate the best fits, especially for the e and µ channels.However, the total χ 2 values for these scenarios are larger than 700, which means that the dark matteronly model cannot explain the radio data of M33.In Fig. 4, we fit the corresponding flux profiles with the best-fit scenarios and we can see that the fits are very poor indeed for both NFW and Burkert profiles.Therefore, the dark matter-only model following the standard cosmological annihilation cross section is unlikely to account for the radio flux profiles.
Then we consider the baryon-only model (i.e.there is no dark matter contribution to the radio signal).For each frequency, we have a unique proportional constant k.By minimizing the reduced χ 2 values, we can obtain the best-fit value of k for each frequency.In fact, many radio studies usually assume that the radio flux and the frequencies have a power-law relation (Tabatabaei, Krause & Beck 2007).This implies that the proportional constant k for different frequencies can be theoretically connected with a power law in frequencies as well.However, there are some cases in which the radio flux spectrum deviates significantly from a power-law description (Lisenfeld et al. 2004).Therefore, we assume that the constant k is a completely free parameter for each frequency to minimize the possible systematic error involved.In Fig. 4, we show the best-fit flux profiles with the best-fit values of k.The reduced χ 2 values are 0.47, 0.11 and 0.04 for ν = 2.70 GHz, 4.85 GHz and 8.35 GHz respectively, which are much smaller than the reduced χ 2 values for the dark matter-only model.Therefore, the baryon-only model can give much better fits compared with the dark matter-only model (see Fig. 4).In other words, the dark matter contribution may be very small compared with the baryonic contribution.
Now we test a more realistic combined model: dark matter plus baryonic contribution.
For the standard cosmological scenario (i.e.σv = 2.2 × 10 −26 cm 3 /s), adding dark matter component does not have a large impact on the total χ 2 values (see Fig. 5).There exist some minimum values of χ 2 for the e channel (m DM = 5 GeV with Burkert profile) and µ channel (m DM = 10 GeV with Burkert profile).However, the statistical significance of these two best-fit cases is very small (< 1.3σ).Generally speaking, as the dark matter contribution is relatively small, m DM ≥ 5 GeV is allowed for all annihilation channels with the thermal annihilation cross section.
Beside assuming the thermal annihilation cross section, we can set the annihilation cross section as a free parameter and obtain the 2σ upper bounds of the annihilation cross section for different m DM and annihilation channels.The annihilation cross section might be larger or smaller than the standard value if dark matter was not thermally produced, or the annihilation cross section is velocity-dependent (Yang et al. 2014).For the dark matter-only model, we can find the best-fit annihilation cross section for different m DM by minimizing the values of χ 2 red .In Fig. 6, we show the best-fit annihilation cross section for each frequency and the case of combining three frequencies.The values of the best-fit cross section are larger than the standard thermal annihilation cross section.In the χ 2 plots shown in Fig. 7, we can see that there exist two cases where the χ 2 is minimum: m DM = 10 GeV for the e channel (NFW profile with σv = 1.3 × 10 −24 cm 3 /s) and m DM = 40 GeV for the µ channel (Burkert profile with σv = 1.3 × 10 −23 cm 3 /s).Although the total χ 2 values for these two cases are slightly larger than the value in the baryon-only model, they represent very good fits for the data.We show the corresponding fits of the radio data in Fig. 8. Nevertheless, when comparing the best-fit results for this dark matter-only model with the AMS-02 positron and gamma-ray analyses, they are excluded by the 2σ upper limit bands on σv (Fermi-LAT Collaboration 2015; Cavasonza et al. 2017).For instance, the 2σ upper limits derived from the gamma-ray data of the Milky Way dwarf galaxies are ∼ 2 × 10 −26 cm 3 /s (m DM = 10 GeV) and ∼ 9 × 10 −25 cm 3 /s (m DM = 40 GeV) respectively for the e channel and µ channel (Fermi-LAT Collaboration 2015).Note that the unknown functional form of the dark matter density assumed in the gamma-ray analysis (i.e. the so-called Jfactor) and whether the annihilation cross section is velocity-dependent would contribute large systematic uncertainties in the derived bounds (Chiappo et al. 2019;Boucher et al. 2022).Now we consider the combined model (dark matter plus baryons) again.In principle, we can also obtain the best-fit annihilation cross section by tuning the value of k(ν) for each frequency.However, as the baryon-only model can provide excellent fits to the radio data, adding a dark matter component does not improve the fits significantly.There exist multiple degenerate best-fit annihilation cross sections for different values of m DM .This suggests that the contribution of dark matter annihilation in the disk region must be small compared with the baryon emission.Nevertheless, we can obtain the 2σ upper bounds of the annihilation cross section for different m DM and annihilation channels.By adopting the baryon-only model described above, we investigate on how large of the dark matter contribution would exceed the 2σ margins.In Fig. 9, following the combined model, we show the overall 2σ upper limits of the annihilation cross section against m DM for different channels by combining the radio data of the three frequencies.These are the strongest limits for dark matter annihilation cross section in the disk region of M33 galaxy as we have included baryon contribution in the radio emission.

Discussion
In this article, we present an analysis combining the radio flux profile and multiwavelength approach to get more comprehensive constraints on dark matter annihilation.Here, we use the archival data of the M33 galaxy to get the radio flux profiles for ν = 2.70 GHz,4.85 GHz and 8.35 GHz.We find that the baryon-only model can give much better fits compared with the dark matter-only model.This means that the overall radio emission in the M33 disk is likely dominated by baryon-related emission (i.e.thermal and non-thermal emissions), but not dark matter contribution.Nevertheless, if annihilating dark matter dominates the radio emission, we find that dark matter with m DM = 10 GeV annihilating via the e channel and m DM = 40 GeV annihilating via the µ channel can give the best fits for the radio data.However, the best-fit annihilation cross sections are excluded by the 2σ upper limit bands on σv (Fermi-LAT Collaboration 2015; Cavasonza et al. 2017).
Furthermore, we apply the combined model to include both contributions of dark matter and baryons.Although there is no best-fit scenario for dark matter due to the degeneracy problem, we can constrain the 2σ conservative upper bounds of the annihilating cross section for different annihilation channels.Compared with the ANTARES neutrino bound (ANTARES collaboration 2020), our upper limits are tighter for the b and W channels, and the µ channels with m DM ≤ 100 GeV.Also, our limit for the µ channel with m DM ∼ 100 GeV is generally more stringent than that in the AMS-02 antiproton analysis (Calore et al. 2022).However, our limits for the b, τ and W channels are less stringent than that in gamma-ray analysis (H.E.S.S. Collaboration 2021; McDaniel et al. 2023).Overall speaking, analyzing the radio flux profile data with different frequencies can give good limits for leptophilic channels, especially for the µ channel.The results for the µ channel are particularly important as recent studies have paid more attention to the µ-related annihilation (Ghosh et al. 2021;Abdughani et al. 2022) due to the muon (g − 2) anomaly in the measurement (Abi et al. 2021).For the b and W channels, gamma-ray detection can generally provide more stringent constraints.Note that for all analyses of dark matter annihilation cross section upper limits, there are different assumptions and uncertainties involved such as dark matter density profile, cosmic-ray propagation parameters, magnetic fields, etc.Therefore, there is no single model-independent robust upper bound for the annihilation cross section.Nevertheless, all of the bounds can be complementary to each other to give more comprehensive constraints for dark matter annihilation.Our analysis using radio data of M33 can provide one of the important contributions.Moreover, compared with the previous studies of the M33 galaxy using a single frequency radio data (Chan 2017), the constraints in this study are more realistic and authentic for annihilating dark matter.
Note that we did not consider any boost factor due to dark matter substructures in our analysis (Moliné et al. 2017;Chan 2017).Besides, our 2σ upper limits of the annihilation cross section are based on the combined model, but not the dark matter-only model.Therefore, in principle, our constraints presented here are stronger compared with that without considering baryonic contributions, provided that the baryonic model assumed is physical.For the case of assuming the annihilation cross section as a free parameter, our results are useful to constrain the non-thermal annihilating dark matter and the velocity-dependent annihilation cross section.For instance, the radio 2σ upper limits of the annihilation cross section shown in our study are generally weaker than that in Beck & Sarkis (2023).However, the constraints in current study are derived from a small galaxy while the constraints in Beck & Sarkis (2023) are derived from galaxy clusters.If the annihilation cross section is velocity-dependent, such as Sommerfeld enhanced (Sommerfeld 1931;Yang et al. 2014), the constraints derived in small galaxies and galaxy clusters would have different implications because the dark matter velocity dispersion in a small galaxy is much smaller than that in a galaxy cluster.
In this study, we demonstrate a comprehensive approach in analysing the radio flux profiles with different frequencies.The overall uncertainties of the radio flux density observed for M33 galaxy are quite significant so that we cannot get very stringent constraints for dark matter.Since there is a resolution limit (i.e. the minimum beam size) for radio observations, we can only obtain good quality radio flux profile data (with small uncertainties) for nearby galaxies (e.g.Local group galaxies).If one can obtain good quality radio map with different frequencies for nearby galaxies, we can follow this comprehensive analysis to get more stringent constraints for annihilating dark matter.

Fig. 1 .
Fig. 1.-The radio flux density map of the M33 galaxy with ν = 4.85 GHz.The region inside the white circle defines our interested region of the M33 galaxy.The colours in the map indicate different radio flux density in mJy/beam.

Fig. 2 .
Fig. 2.-The radio flux density (per beam size) profiles of the M33 galaxy for ν = 2.70 GHz, 4.85 GHz and 8.35 GHz.The red data points indicate the radio flux density in the bulge region while the black data points indicate the radio flux density in the disk region.
Fig. 3.-The graph at the left upper corner indicates the total χ 2 as a function of m DM for different annihilation channels in the dark matter-only model with the thermal annihilation cross section.The other three graphs are the reduced χ 2 for different frequencies.The colored solid and dashed lines indicate the dark matter density distributions following the NFW and Burkert profiles respectively.

Fig. 9 .
Fig. 9.-The 2σ upper bounds of the annihilation cross section for different annihilation channels in the combined model.The colored solid and dashed lines indicate the dark matter density distributions following the NFW and Burkert profiles respectively.