{'fraction_non_alphanumeric': 0.038837137308980384, 'fraction_numerical': 0.015422283511102355, 'mean_word_length': 4.246419607264137, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The EMU view of the Large Magellanic Cloud: Troubles for sub-TeV WIMPs

We present a radio search for WIMP dark matter in the Large Magellanic Cloud (LMC). We make use of a recent deep image of the LMC obtained from observations of the Australian Square Kilometre Array Pathfinder (ASKAP), and processed as part of the Evolutionary Map of the Universe (EMU) survey. LMC is an extremely promising target for WIMP searches at radio frequencies because of the large J-factor and the presence of a substantial magnetic field. We detect no evidence for emission arising from WIMP annihilations and derive stringent bounds on the annihilation rate as a function of the WIMP mass, for different annihilation channels. This work excludes the thermal cross section for masses below 480 GeV and annihilation into quarks.


Introduction
The nature of dark matter (DM) is one of the defining mysteries of modern physics. One of the most attractive candidates proposed to solve the DM puzzle is given by hypothetical particles that are more massive than baryons and weakly interacting, the so-called WIMPs [1].
WIMPs have masses in the GeV-TeV range and can annihilate in pairs into lighter particles. In particular, electrons and positrons can be directly or indirectly injected by WIMP annihilations, and a sizable final branching ratio of annihilation into e + − e − is a rather generic feature of WIMP models (see, e.g., Fig. 4 in Ref. [2]). Emitted in an environment with a background magnetic field, high-energy electrons and positrons give rise to synchrotron radiation typically peaking at radio frequencies.
The Large Magellanic Cloud (LMC) is the most massive satellite galaxy of the Milky Way. The large dark matter mass (with the LMC virial mass being around 10 11 M [3]) and the proximity to Earth (distance around 50 kpc [4]) make the LMC one of the best targets for indirect searches of WIMPs. The so-called J-factor, namely, the integral of the density squared over the line-of-sight and solid angle, amounts to ∼ 10 20 GeV 2 /cm 5 [5], second only to the Galactic center. The presence in LMC of a µG magnetic field [6] suggests an investigation of a possible WIMP signature at radio frequencies.
The idea of deriving bounds on WIMPs from radio observations of the LMC is not new [7,8]. The improvement presented in the current analysis with respect to Refs. [7] and [8] arises from new, more sensitive, data, an updated model (with the inclusion of spatial diffusion in the computation of the DM signal), and the choice of statistical approach (comprising a morphological analysis with pixel by pixel comparison of the observed and predicted signals).
We use observations of the LMC obtained by the Australian Square Kilometre Array Pathfinder (ASKAP [9]), as part of the ASKAP commissioning and early science (ACES, project code AS033) verification and made available to the Evolutionary Map of the Universe (EMU) project [10]. These data were used to obtain a deep image of the LMC at 888 MHz [11], which will be the starting point of the analysis of this work.
The structure of the paper is as follows. Section 2 describes the ASKAP observations and how the LMC radio image has been created. The model of DM and interstellar medium in the LMC is presented in Section 3. We introduce the statistical analysis and present the results in Section 4, while Section 5 provides a comparison with previous works. Conclusions are drawn in Section 6. The Appendix is devoted to consistency checks.

Observational Maps
In this work we make use of the observations of the LMC taken as part of the ASKAP commissioning and early science at 888 MHz with a bandwidth of 288 MHz, and analysed as part of the EMU project [11]. The observations cover a total field of view of 120 deg 2 , with a total exposure time of about 12h40m. Data processing was performed using ASKAPsoft by the ASKAP operations team and the resulting images are available on the CSIRO ASKAP Science Data Archive. The beam size of the map shown in Fig. 1 (left) is 13.9 × 12.1 and the median Root Mean Squared (RMS) noise is ∼ 58 µJy/beam. For more details see Ref. [11].
Structures on scales 2 • can be recovered thanks to the short baselines of the ASKAP array (with shortest one being 22 m). We expect the DM-induced emission to be diffuse but showing variations on scales below 2 • (see next Section). We confirm the sensitivity of the image to DM diffuse emissions by performing a test detailed in the Appendix. Our search looks for a possible diffuse emission associated to the LMC halo, and all the small-scale discrete sources (in the LMC or in the background/foreground) are a contaminant that we attempt to mask. We identify discrete sources by running the publicly available tool SExtractor [12], which is also used to derive the RMS map, in the same way as described in Ref. [13]. The threshold for a source to be masked is set to 3 × RM S, and the result is shown in Fig. 1 (right). We also mask negative pixels using the same criterion, i.e., absolute value larger than 3 × RM S. They are likely due to missing short-spacing data.
Since the expected emission from DM has a size of several arcmin, we further smooth the masked image (using the task SMOOTH in Miriad [14]) to FWHM= 2 , in order to be more stable against small-scale residuals and fluctuations.  [11]. Right panel: Same as left panel but blanking masked pixels. The red cross indicates the position of the dynamical center and the blue circle encloses the region of interest for this work.

LMC description
We compute the radio emission induced by WIMP DM by combining the synchrotron power associated with the LMC magnetic field with the equilibrium distribution (n e ) of electrons and positrons injected by DM. In order to compute n e , we solve a transport equation describing the cooling and spatial diffusion experienced by the electrons and positrons after injection. We describe it in the limit of spherical symmetry and stationarity: where f (r, p) is the e + /e − distribution function at the equilibrium, 1 at a given radius r (from the LMC dynamical center) and at a given momentum p. The distribution f is related to the number density in the energy interval (E, E + dE) by: n e (r, E)dE = 4π p 2 f (r, p)dp; analogously, for the e + /e − source function, we have q e (r, E)dE = 4π p 2 s(r, p)dp. The first term on the left-hand side describes the spatial diffusion, with D(r, p) being the diffusion coefficient. The second term accounts for the energy loss due to radiative processes;ṗ(r, p) = i dp i (r, p)/dt is the sum of the rates of momentum loss associated to the radiative process i. The DM source q e scales with the number density of WIMP pairs locally in space, i.e., with ρ 2 /(2 M 2 χ ), where ρ(r) is the halo mass density profile, and M χ is the mass of the DM particle. 2 We neglect substructure contributions and assume the DM spatial distribution to be spherically symmetric and static. The source term associated to the production of e + /e − is given by: where σ a v is the velocity-averaged annihilation rate, and dN a e /dE is the number of electrons/positrons emitted per annihilation in the energy interval (E, E + dE) for a given annihilation channel.
We solve Eq. (3.1) numerically using finite-differencing Crank-Nicolson scheme, for details see Ref. [13]. Boundary conditions are set to be Neumann's one at the centre and Dirichlet's one at the farthest boundary, the latter chosen to be ten times the radius of our region of interest (RoI). A recent semi-analytical treatment of Eq. (3.1) can be found in Ref. [15].
The properties we want to constrain are the DM mass and annihilation rate, while the ingredients we need to model are the DM spatial profile, the magnetic field, the spatial diffusion coefficient, the CMB and LMC interstellar radiation fields (ISRF, for inverse Compton losses) and the gas density (for bremsstrahlung losses), that we describe in detail in the following sections. Since our goal is to derive conservative bounds on the WIMP signal, we will model the above quantities taking lower limits for the DM profile and magnetic field, while upper limits for diffusion coefficient and ISRF and gas densities.
To limit uncertainties in the model description, our RoI for the analysis will be defined by a relatively small region around the LMC center, corresponding to 1.3 kpc in radius (1.5 • in angular units). The loss in terms of J-factor, if compared to considering the full LMC halo, is limited, around a factor of two, depending on the profile. As we will describe in the following, in such RoI we have a more robust determination of the various ingredients entering the computation, such as the magnetic field, the gas and ISRF distributions, and the DM profile. Moreover, we exclude the bulk of the contamination from the 30 Doradus region (south-west in Fig. 1).

Dark matter profile
To model the radio emission of the LMC due to annihilating DM, the spatial particle distribution ρ(r) is a key ingredient, see Eq. (3.2). Previous work has explored different functional forms for the DM density profile in the LMC [5,8,[16][17][18][19][20] and analysed H i rotational [21] and carbon star data [22] to constrain the parameters of ρ(r). In this work we are interested in the inner region, so we make use of the H i data [21] that provide the most accurate rotation velocities at small distances from the LMC dynamical center. We explore four different profiles and fit ρ(r) up to a radius of ∼ 2.7 kpc from the center 3 , which corresponds to about twice our RoI. In particular, we consider two different "cuspy" DM profiles from the NFW model (Navarro-Frenk-White) [23] and Hernquist [24]: and two "cored" profiles, the isothermal sphere [25] and the Burkert profile [26]: The main reason for considering different shapes is related to our poor knowledge about DM physics at small scales, including the possible role of baryons in affecting the DM spatial profile. The range of possibilities encompassed by the above functional forms should bracket the uncertainty. The free parameters of the different profiles are the scale radius r s and the normalization ρ s . We fit these values for each profile using the H i rotation velocity data and the fact that the velocity v(r) measured up to a radius r is given by the expression v(r) = GM tot (r) r , (3.5) where M tot (r) is the total mass contained within a radius r, given by DM plus contributions from the stellar and gas components. We model the stellar potential φ (R, z) using a Plummer-Kuzmin disk [27] as a function of the disk radial distance R and vertical height z in cylindrical coordinates: where a and b are the radial scale length and vertical scale height, respectively, for which we take a = 1.7 kpc and b = 0.34 kpc [28,29]. The stellar mass M is left as a free parameter in the fit. The contribution to the mass density ρ g (R, z) from the gas follows the expression given in Ref. [29] (with radial scale length a g = a and vertical scale height b g = b ): Once we have the total mass contribution from the different components to the rotation velocity at some radius r, we proceed to fit the parameters r s , ρ s and M using a least squares method through the python package scipy.optimize. The best-fit parameters for the different DM density profiles are shown on Table 1. For simplicity, and since the gas component provides a subdominant contribution to the matter density, we set M g = 5 × 10 8 M . Nevertheless, we checked that the results for r s and ρ s reported on Table 1 are unchanged if M g is left as a free parameter in the fit.
Previous work has suggested that the LMC virial mass is around 2 × 10 11 M , but estimates can have significant uncertainties [3,18,20]. On the other hand, the virial mass is mostly related to the DM profile at larger radii than the one relevant for our analysis (the enclosed mass in our RoI is 2 × 10 10 M ), making these uncertainties of little relevance for our results. We adopt different values for the LMC virial mass, corresponding to Table 2 from Ref. [20], in order to normalize the profiles. We find that the rotation velocity within ∼ 5 kpc from the center of the LMC, as well as the profile parameters determined by our fit, do not change considerably with different choices for the mass normalization. Therefore we adopt a (low) virial mass of 10  The results are shown in Figure 2, where we report the contributions to the rotation velocity from DM (dashed lines), stellar and gas components (red dots and green crosses respectively). The orange points represent the H i rotation data from Ref. [21] and the solid lines show the contribution from the sum of all components.

Magnetic field and diffusion coefficient
The strength of the large-scale coherent component of the LMC magnetic field is found to be 1 µG, as determined via Faraday rotation measure of polarized background sources [6], and with diffuse polarized data [30] (see also [31]). The turbulent component is expected to be larger than the regular one, by a factor > 3, as generically found in galaxies (see, e.g., Ref. [32]), and confirmed by the scatter observed in the LMC rotation measures [30]. In Ref. [6], the total LMC magnetic field strength on large scales has been estimated to be B = 4.3 µG, and we take this as our reference value. We focus on a relatively small and central region of LMC, where rotation measures do not show significant dependence on the radial distance, so we can assume a uniform strength, in agreement with the model in Ref. [30]. Amplifications on small scales [6] are disregarded. A recent analysis based on the equipartition assumption and on observations of the LMC synchrotron emission at 166 MHz [33] suggests a higher value, B = 7.7 ± 1.1 µG. Since the magnetic field is (together with the DM properties) the most crucial ingredient of our analysis, we show how our results change for a range of total strength B = 2 − 8 µG. This range is in agreement with estimates stemming from the cosmic-ray density derived from γ-ray data and again applying the equipartition assumption [30].
Data on supernovae remnants [34] and on large-scale diffuse emission [35] indicate that the transport of cosmic-rays in the LMC proceeds in a similar way as in other nearby galaxies, and can be explained as diffusive propagation in a turbulent regime. In this scenario, the diffusion coefficient can be estimated as [35]: where d L is the diffusion length and τ is the cooling time. Radio observations suggest d L 1 − 2 kpc in the vertical direction (larger along the disk), consistent with findings in other galaxies where the confinement region is a few times the disk height. The discussion on radiative losses below leads to τ 10 15 s, in agreement with the limit of τ > 10 14 s estimated in Ref. [34]. Thus observations point towards a value somewhat lower than in the Galaxy. For clarity, and in the spirit of making conservative assumptions, we assume the same strength and energy dependence of the diffusion coefficient as at large scales in the Galaxy, taking the latest determination from Ref. [36] (BIG model, which provides 2 × 10 28 cm 2 /s at 4 GeV). Recall that the larger the diffusion coefficient the smaller the DM signal, since diffusion can remove electrons and positrons from the RoI before they emit synchrotron radiation at the frequency of interest.

Gas and interstellar radiation fields
We determine the central value of the gas density from Eq. (3.7) and taking M g = 5 × 10 8 M , which is the neutral hydrogen mass observed by Ref. [21]. Then we assume a flat spatial profile, normalized to the maximal value, i.e., the central value, which leads to the gas number density n g = ρ(0, 0)/m H = 0.8 cm −3 (where m H = 0.938 GeV is the hydrogen atom mass). This is clearly a conservative approach. Moreover it simplifies the computation by avoiding uncertainties related to the radial and vertical scale lengths of the gas distribution and allowing us to keep assuming spherical symmetry. We checked this description by deriving the gas density from the hydrogen column density in Fig.4 of Ref. [37] divided by a disk thickness of 350 pc [37]. The spatial profile does not show significant variations in our RoI (justifying a flat model) and the average value for the column density is 6 × 10 20 cm −2 , which translates into n g 0.5 cm −3 , confirming the above assumption for n g as an upper limit.
We assume the gas to be composed solely of neutral atomic hydrogen, since molecular hydrogen and ionised gas are subdominant components [6,38,39], negligible in this analysis.
The ISRF spectrum is taken to have the same shape as that of the Milky-Way [40]. Observationally, this is found to be a good approximation [41]. Moreover, even though we implement a full computation, the Klein-Nishina corrections are subdominant (the energy of the emitting electrons is 10 GeV), so the exact ISRF spectral shape is not critical, and the size of the inverse Compton losses is essentially set by the integral over energy.
The normalization is chosen such that the integral of the spectrum provides U ISRF = 1 eV/cm 3 , consistent with the parameter X ISRF found in Ref. [41]. We note that a somewhat lower density can be derived from the LMC SED of Ref. [42], see Ref. [43] who found U ISRF = 0.57 eV/cm 3 , and using the estimate of Ref. [35] where U ISRF 0.3 eV/cm 3 . Again, our choice is in the spirit of adopting a realistic upper limit.
As for the gas density, we conservatively take a spatially flat profile.

Results
We assume the likelihood associated with the LMC diffuse emission to be described by a Gaussian: where S i th is the theoretical estimate for the flux density in the pixel i, S i obs is the observed flux density and σ i rms is the r.m.s. error, both described in Sec. 2. N pix is the total number of pixels in the RoI (excluding masked pixels) and N F W HM pix is the number of pixels within the FWHM of the synthesized beam. We only include the DM signal coming from inside a sphere of radius of 1.3 kpc around the LMC center (thus disregarding other line-of-sight DM contributions inside the angular region of 90 ). The theoretical estimate is provided by the WIMP emission, computed from the solution of Eq. (3.1) and following Sec. 4.1 in Ref. [13], plus a disk component and a spatially flat term.

The disk is described through a Gaussian
where θ d,i is the angular distance of the pixel i from the axis of the disk, and S 0 and θ 0 are free parameters. The position angle of the LMC disk has been found to be between 122.5 • and 170.5 • , depending on the tracer (see Ref. [44] and references therein). In our analysis we assume a value determined by fitting the map without including the DM component. We find 138 • , similar to that found recently using Gaia DR2 data [45].
In Fig. 3 (left), we show the shape of the model given by disk plus DM. We use arbitrary normalization and fix the FWHM of the Gaussian describing the disk to 0.45 • (which is the best-fit value found in the fit).
On top of the disk component we add a spatially flat term S f lat . This is included in the fit to account for possible offsets or a large-scale foreground component. The parameters S f lat , S 0 and θ 0 are treated as nuisance parameters.
In Fig. 3 (right), we show the WIMP emission as a function of the angular distance from the center, for different masses, annihilation channels and DM density profiles. Note that the size of the DM source is below 2 degrees. The NFW profile is the most concentrated case (together with the Hernquist profile, which is not shown since it is nearly identical to the NFW at small distances). The Burkert and isothermal profiles provide more extended distributions. Note that high masses and leptonic channels imply a less concentrated profile than low masses and hadronic channels. This can be understood by noting that at the frequency of the observations (888 MHz) and for a magnetic field of 4.3 µG, the synchrotron emission is mostly provided by e + /e − with energy around a few GeV. High energy electrons take time to cool down to few GeV and thus can travel long distances, flattening the central overdensity. Recall that leptonic channels provide harder e + e − spectra than in the bb case.
Bounds on the parameter σ a v are computed at any given mass M χ through a profile likelihood technique [46], namely "profiling out" the nuisance parameters Π = (S f lat , S 0 , θ 0 ). We assume that follows a χ 2 -distribution with one d.o.f. and with one-sided probability given by P = ∞ √ λc dχ e −χ 2 /2 / √ 2 π, where σ a v b.f. denotes the best-fit value for the annihilation rate at that specific WIMP mass. Therefore, the 95% C.L. upper limit on σ a v at mass M χ is obtained by increasing the signal from its best-fit value until λ c = 2.71, keeping Π fixed to its best-fit value.
The possible presence of a DM signal is investigated by evaluating the difference , We always find ∆χ 2 < 1, and thus no evidence for a diffuse component associated with WIMP-induced emission. Results are shown in Fig. 4, reporting the upper limits on σ a v . The boundaries of the uncertainty band are determined by taking the weakest/strongest bound among those obtained using the four different DM profiles described in Section 3.1. More concentrated profiles provide more stringent constraints. The NFW and Hernquist cases set the lower boundary of the band, while Burkert at low masses and Isothermal at high masses set the upper boundary. The dashed black line is the so-called thermal cross section, namely the self-annihilation cross section needed in the early Universe in order to provide the DM mass density observed today [47]. A common way to see Fig. 4 is to consider "canonical" WIMPs excluded for masses where the bound is below the thermal value.
The trend of the bound is similar for the bb (blue) and W + W − (green) channels, on one hand, and for τ + τ − (red) and µ + µ − (orange) channels, on the other. The reason is related to the injection spectrum of e + e − . Let us first remind that the key quantity is the density of e + /e − induced by the specific DM model at energies of a few GeV. The injection spectrum of e + e − is harder in the leptonic cases, where WIMPs with mass of tens of GeV have therefore a significant injection of e + e − with energy around the peak of the synchrotron power. This makes the bounds in the τ + τ − and µ + µ − cases very tight at low masses. Clearly, the picture is the opposite at high masses where the injection energy is "too high" and the e + e − undergo energy losses and diffusion before emitting synchrotron radiation. In the cases of bb and W + W − , the peak in the injection of e + e − occurs at around M χ /20. Therefore, they are more efficiently constrained in the range of masses around hundreds of GeV (so that, again, the production of e + /e − is peaked around a few GeV).
Since we consider non-relativistic DM, the WIMP mass has to be larger than the mass of the annihilation products, and this is the reason of the cut in the W + W − bound.
The overall increase of the constraints with the WIMP mass, occurring for all the channels, can be understood from Eq. (3.2).
The bottom line of Fig. 4 is that the thermal cross-section is excluded for masses below (480,358,192,164) GeV for the (bb, W + W − , τ + τ − , µ + µ − ) annihilation channel.  In Fig. 5 (left), we show how the magnetic field strength affects the bound on σ a v . We consider an example with M χ = 100 GeV, annihilation into bb, and NFW DM profile. The bound approximately scales with the inverse of the square of the magnetic field for small strengths, and flattens as the strength increases. As already stated, in Fig. 4 we adopted B = 4.3 µG.
Throughout the paper, we have been assuming annihilating DM. In the right panel of Fig. 5 we show the bounds that can be obtained for decaying DM. The only difference from the above analysis consists in the replacement of Eq. (3.2) with where Γ d is the decay rate, and dN d e /dE is the number of electrons/positrons emitted per decay in the energy interval (E, E + dE).
The different behaviour of the four decaying channels can be understood in a very similar way to that already discussed above for the annihilating case. Note that the uncertainty band of the curves in Fig. 5 is smaller than in the annihilating cases of Fig. 4. This is because the source function of annihilating DM depends on ρ 2 , while the decaying scenario scales linearly with ρ, and thus uncertainties in the DM spatial profile affect the former more than the latter.

Comparison with radio analyses
An analysis similar to the one presented here was conducted in Ref. [8]. They employed ATCA+Parkes data and obtained the black curve in Fig. 6 (left). The great improvement in the constraining power of our analysis is not due to the model, for which we adopted a more conservative description than Ref. [8], but to the different statistical approach and to the more limited residuals and higher rms sensitivity in the ASKAP image compared to the ATCA+Parkes image. Concerning the statistical analysis, Ref. [8] derived the bound from individual lines of sight, while we developed a morphological analysis. The latter allows us to ascribe part of the LMC emission to a disk component and combines ∼ 10 4 lines of sight. This is important since the constraining power scales roughly as the square root of this number.
There have been a few attempts to detect WIMP-induced radio signals in dwarf spheroidal galaxies of the Local Group. We expect the signal from LMC to be stronger than from dSphs since the magnetic field strength is higher (it is actually unknown in dSph, but typically assumed to be around 1 µG [13]), the J-factor is higher than (or at the level of) that in the most promising dSphs, and the LMC is bigger (which means diffusion effects are less relevant in depleting the signal than in dSphs). In Fig. 6, we include bounds derived from the observations of different samples of dSphs with the ATCA from Ref. [57] (taking their model with B = 1 µG , D 0 = 3 × 10 28 cm 2 /s, red solid line) and Ref. [58] (AVE model, red dashed line), GMRT [59] (B = 2 µG , D 0 = 3 × 10 28 cm 2 /s, violet line), LOFAR [60] (B = 1 µG , D 0 = 10 27 cm 2 /s, orange line), MWA+GMRT [61] (B = 2 µG , D 0 = 3 × 10 28 cm 2 /s, magenta line). Other relevant campaigns have been conducted with the GBT [62][63][64] and MWA [65]. Their bounds are not in a suitable form to be shown in Fig. 6, but correspond to about σ a v 10 −24 cm 3 /s for M χ = 100 GeV and the bb channel.

Comparison with γ-ray analyses
In Fig. 6 (right), we compare the results of this work with the bounds obtained by the Fermi-LAT Collaboration from the analysis of the LMC [16] and dSphs [66]. For completeness we also show the expected LMC bounds from the Cherenkov Telescope Array [67], since they can be more constraining than those from Fermi-LAT at high WIMP masses.
One can immediately see that the LMC radio constraints are much stronger than γ-ray ones. This should not come as a surprise. Indeed, generically, in WIMP models, the luminosity associated with the induced γ-rays is comparable to or smaller than that of the injected electrons and positrons. For hadronic channels, the γ-ray emission mainly proceeds through the production and decay of neutral pions, while electron/positron injection is related to charged pions, and so the two mechanisms are tightly related. For leptonic channels, GeV electrons and positrons have a larger yield than the final state radiation of γ-ray photons. Therefore, if the cooling time/diffusion length is small enough, so that the energy of the electron/positron is radiated within the source, and the synchrotron loss is the dominant (or, at least, among the most relevant) radiation mechanism in the astrophysical object under investigation, the luminosity produced as synchrotron radiation is comparable to or higher than that from γ-rays [2]. Since radio telescopes are much more sensitive than γ-ray telescopes for all sources having related mechanisms of emission in the two bands (see, e.g., the level of detail in the ASKAP LMC image compared to the γ-ray image of the LMC [68]), radio bounds on WIMPs are significantly stronger, when above conditions are satisfied, and in particular for objects with low astrophysical diffuse radio background, such as the LMC.
The picture is different for dSphs, since there the diffusion length is typically larger than the galaxy itself and the magnetic field is thought to be rather small (and so too the synchrotron radiation), which implies a less favourable ratio between radio and γ bounds, even though still comparable (see red solid line in the left panel versus orange line in the right panel). Fig. 6 shows that the bound derived in this work is the most stringent bound on WIMPs coming from indirect searches in extragalactic objects. Comparison with other radio searches, see text for details. "DW" labels analyses of dwarf spheroidal galaxies. Right panel: Comparison with γ-ray Fermi-LAT searches on dwarf spheroidals [66] (orange) and on the LMC [16], and with the expected sensitivity of CTA for LMC [67].

Conclusions
We analysed the ASKAP radio image at 888 MHz of the LMC, in order to search for synchrotron emission induced by WIMP DM annihilations.
The large J-factor of the LMC implies it is one of the best targets for DM indirect searches. The presence of a magnetic field with strength > 1 µG makes radio searches in the LMC particularly suited for this purpose.
We detect no evidence for emission arising from WIMP annihilations and derive stringent bounds. Annihilations into leptonic channels provide the most constraining bounds at low masses with the thermal cross-section excluded for masses below 192 GeV (τ + τ − ) and 164 GeV (µ + µ − ). Annihilations into quarks and gauge bosons are the most constraining cases at intermediate and high masses with the thermal cross-section excluded below 480 GeV (bb) and 358 GeV (W + W − ).
The comparison with the state-of-the-art in Fig. 6 shows that the bounds on WIMPs derived in this work are extremely competitive. We adopted a simple and conservative approach, limiting the analysis to a relatively small region, where simplified assumptions and well-motivated data-driven descriptions can be taken for the astrophysical ingredients entering the model prediction. For the two most relevant quantities, the DM spatial profile and the magnetic field, we defined reference models according to observations [6,21]. For the other components (which are important for the computation, but to a somewhat lesser extent), such as the spatial diffusion coefficient, the interstellar radiation fields, and the gas density, we consider their upper or lower limits (all in the direction of minimising the DM signal).
Our results imply there is little hope to detect a signal from low mass thermal WIMPs in laboratories (i.e., in direct and collider searches), whilst very massive thermal WIMPs remain a viable DM candidate. They can be probed by different techniques, including observations from future radio telescopes, such as the SKA, in particular going to higher frequencies. On a shorter timescale, the addition of short spacing data coming from forthcoming observations with the Parkes telescope will provide a complete picture of the LMC at different scales. With such image at hand, a more refined 3D modeling of the synchrotron emission from the entire LMC can be attempted, with the possibility of further tightening the bounds derived in this analysis. interested in understanding the trend but not the absolute value, we performed the analysis on Stokes V and considering one of the LMC pointings (all pointings were taken at similar directions and times). The rms sensitivity normalized to one at 2 arcmin is plotted in Fig. 7 (red line) as a function of the tapering size.
In the same figure, the blue line reports the total flux of a source which is excluded at 95% C.L. by the analysis described in the main text, as a function of the angular scale of the source. Again values are normalized to one for a source with FWHM= 2 arcmin. One can quickly check that, in the case of uniform rms, no masking and no confusion, Eq. (4.1) would imply a linear scaling. The blue curve is derived by considering Gaussian sources of different widths (essentially replacing the DM component with a Gaussian emission and then repeating all the steps described in Sects. 3 and 4). The actual behaviour is close to linear scaling.
The bottom-line of Fig. 7 is that the "theoretical" degradation of the sensitivity as the angular scale of the sources increases, as assumed by the analysis described in the paper, is higher if compared to the results obtained through the tapering test. This ensures our bounds are conservative.
To test against possible systematics associated to the selected region of the image, we re-do the same analysis outlined in the main text, but now centering the DM distribution on different positions across the map. Since we are dealing with a non-detection, they should all provide similar bounds, because the RMS sensitivity is approximately uniform across the map. We compare five different positions, as listed in Fig. 8 (left). Concerning the DM model we take the same description as in Section 3, but centered in the new positions. This is clearly not realistic, but functional to our test. We find that positions that are far from the LMC disk (i.e., C and D in Fig. 8) provide slightly more constraining bounds, with the component S disc compatible with zero. Positions located on the LMC disk (i.e., B and E in Fig. 8) lead to bounds similar to the ones described in the main text (case A), with the fit requiring a disk component different from zero. Fig. 8 (left) reports the bound in the case of the NFW profile and annihilation into bb. In Fig. 8 (right) we compare our results with the maximal sensitivity that can be achieved with the image we have at hand. The latter is derived by keeping the original resolution (FWHM=13") and evaluating , which can be seen as setting to zero all pixels in the map (after RMS determination). We show the result for the Burkert profile since it is the most extended case, so where the number of pixels relevant for the χ 2 determination is largest, which implies the sensitivity difference is largest. As expected the bound derived with such ideal sensitivity is more constraining than for our reference analysis, but by a rather limited factor (between 2 and 3).
In the same Figure, we derive WIMP bounds from the LMC map at 1.4 GHz presented in Ref. [69]. Such map contains all scales above 40 being a combination of ATCA and (single-dish) Parkes data, contrary to the ASKAP image having only interferometric data. This means that the very large-scale emission would need a more careful treatment than the simple model introduced in Section 3. For this reason, we perform the comparison in the "ideal-sensitivity" case. The analysis of ATCA+Parkes data is performed in the same way as for the ASKAP map. From the ratio of the RMS sensitivity of the two maps (300 versus 58 µJy/beam), and considering the different frequency (1.4 versus 0.888 GHz) and beam (40" versus 13"), we expect the ATCA+Parkes bound to be a factor around 5 weaker than the ASKAP one. Results are along the line of expectations, see dotted versus dashed lines, providing a consistency check.