Limits on Dark Matter Annihilation from the Shape of Radio Emission in M31

Well-motivated models of dark matter often result in a population of electrons and positrons within galaxies produced through dark matter annihilation -- usually in association with gamma rays. As they diffuse through galactic magnetic fields, these $e^\pm$ produce synchrotron radio emission. The intensity and morphology of this signal depends on the properties of the interstellar medium through which the $e^\pm$ propagate. Using observations of the Andromeda Galaxy (M31) to construct a model of the gas, magnetic fields, and starlight, we set constraints on dark matter annihilation to $b\bar{b}$ using the morphology of 3.6 cm radio emission. As the emission signal at the center of M31 is very sensitive to the diffusion coefficient and dark matter profile, we base our limits on the differential flux in the region between $0.9-6.9$ kpc from the center. We exclude annihilation cross sections $\gtrsim 3 \times 10^{-25}$ cm$^3$/s in the mass range $10-500$ GeV, with a maximum sensitivity of $7\times 10^{-26}$ cm$^3$/s at $20-40$ GeV. Though these limits are weaker than those found in previous studies of M31, they are robust to variations of the diffusion coefficient.


I. INTRODUCTION
To date, all evidence for dark matter comes from its gravitational influence on the visible matter in the Universe.However, the majority of successful models for the production of dark matter require some level of nongravitational interactions between the visible and dark sectors.Perhaps the best-known such scenario is that of thermally-produced dark matter, where a small interaction between dark matter and the Standard Model results in a relic population of non-relativistic particles due to thermal freeze-out during the early Universe.Weakly-Interacting Massive Particles (WIMPs) are the most wellknown implementation of this class of dark matter models.In such models, the observed density of dark matter is obtained if the velocity-averaged annihilation cross section is ⟨σv⟩ ∼ 3 × 10 −26 cm 3 /s, for dark matter in the mass range ∼ 1 − 10 3 GeV [1]. 1  Thermal relics, as well as other models with dark matter-Standard Model interactions of similar magnitude, result in a number of possible experimental signatures.Of particular interest to this work is indirect detection, where present-day residual annihilation or decay of dark matter into Standard Model particles gives visible signatures that can be detected here on Earth.Annihilation to Standard Model particles will generically result in cascade decays terminating in stable e ± , photons, neutrinos, and p/p, evidence of which can reach Earth-based detectors from their astronomical point of origin.As the strength of these signals increases with dark matter density squared and decreases with the distance to target squared, the indirect detection targets with the greatest signal rate are the largest and closest conglomerations of dark matter.The highest intensity signals are therefore expected to be seen from our own Milky Way Galactic Center, but other nearby galaxies - 1 Though model-specific details can easily change these numbers by O(1) factors or more.
High-energy prompt photons, which can either come directly from the annihilation of dark matter or the cascade decays of annihilation products, travel largely unimpeded from where they were produced to Earth.Such photons are therefore the most straightforward indirect detection signal, with a morphology that is set only by the dark matter distribution.Interestingly, many groups have identified an excess of gamma rays in the energy range 1 − 3 GeV from data collected by the Fermi Large Area Telescope (LAT) [14] in the Milky Way [15][16][17][18][19], with morphology compatible with the dark matter expectations.Possible signals consistent with this gamma ray excess have been reported in M31 [2], and the LMC [4], though with less significance.These excesses can be wellfit by dark matter models with m χ ∼ O(10 − 100 GeV) annihilating to either b b or τ + τ − followed by cascade decays which result in the observed photons, with a thermally averaged cross section of ⟨σv⟩ ∼ 2 × 10 −26 cm 3 /s [4,[16][17][18][19][20].
However, the ultimate origin of this gamma ray excess remains unclear.An unresolved population of millisecond pulsars (MSPs) in the center of the Milky Way has been suggested as an alternate source of this signal [21][22][23][24][25][26].Ref. [27] has argued that the distribution of gamma rays in the Galactic Center excess (GCE) in the Fermi-LAT data contains non-Poissonian statistics, suggestive of a MSP origin.At this time, debate appears to be far from settled, with questions about the spectrum of MSPs [28], morphology and background emission modelling [29][30][31], and the non-Poissonian statistics interpretation [32,33] all remaining open.A recent analysis [34] suggests that the observed excess is best fit by a combination of point sources and a diffuse source, but uncertainties are large enough that either origin could dominate.
In this context, searches for indirect detection signals beyond prompt photons are especially interesting.Dark matter annihilation into Standard Model final states which decay into gamma rays will necessarily also have significant branching ratios into electrons and positrons.These e ± will interact with galactic magnetic fields, ambient photons (from starlight, dust and the Cosmic Microwave Background (CMB)) and interstellar gas, losing energy and emitting a range of secondary photons ranging in energy from radio up to X-rays.These signals depend on the properties of the target beyond the dark matter distribution, introducing uncertainties that do not exist in prompt photon searches; however the systematics and backgrounds are largely distinct as well.
In this work, we set constraints on dark matter annihilation by analyzing a 3.6 cm radio map of the Andromeda galaxy (M31) by the Effelsberg telescope [35].M31 has been a common target for dark matter indirect detection searches using radio emission from the center of the galaxy [36][37][38][39]; though the resulting constraints are sensitive to assumptions made about the astrophysical characteristics of the galaxy and the dark matter distribution.
The e ± injection rate (and the associated radio signal) at the center of M31 is dependent on the slope of the dark matter density distribution, which has considerable uncertainties [40].The galactic center radio signal is also dependent on assumptions of the diffusion coefficient.For example, Refs.[36][37][38] assume electrons and positrons lose all of their energy before diffusing a measurable distance, predicting larger signals in this region than analyses which assume greater diffusion (as considered in as Ref. [20]).
In order to set robust limits which are less sensitive to reasonable variations of the astrophysical parameters, we consider the morphology and intensity of the radio emission from the region of M31 between 0.9-6.9kpc from the galactic center.We compute the expected synchrotron emission from electrons and positrons in M31 produced through dark matter in the mass range 6−500 GeV annihilating to b b while varying the diffusion coefficient over the range of experimentally allowed values [35,[41][42][43].Commonly used tools for modeling the transport of e ± (such as galprop [44] and rx-dmfit [45]) use a uniform diffusion coefficient.For modeling the relatively large region of interest for our analysis, this assumption is insufficient.We develop a numerical solution that allows for radial dependence in all transport coefficients -including the diffusion coefficient.Using our numerical method, we solve for the spherically averaged electron-positron phase space density and compute the radio emission from this phase space density and an axisymmetric model of the magnetic field.We then set exclusion limits on the annihilation cross section of dark matter, while varying the diffusion coefficient normalization.
The remainder of this paper is organized as follows.In Section II, we describe the radio observations of M31 used in the analysis.Section III describes the spectrum and morphology of electrons and positrons injected into M31 through the annihilation of dark matter.We construct our models of the magnetic fields, interstellar radiation field (ISRF), and thermal gas density using a variety of relevant data in Section IV.The transport of electrons and positrons within M31 is described in Section V.In this section, we include a discussion of the physics of charged particle transport in a galaxy and our numerical method for solving the transport equation for systems with position-dependent energy loss and diffusion.In Section VI, we calculate the intensity and morphology of the resulting synchrotron emission.Our statistical method for determining a data-driven background model and setting exclusion limits is described in Section VII.In Section VIII, we present our results.Finally, we conclude in Section IX.

II. RADIO OBSERVATIONS OF M31
To constrain dark matter annihilation into high-energy electrons and positrons in M31, we use the non-thermal radio flux per unit frequency per beam dS/dν from a survey of ν = 8.35 GHz emission in M31 [35] using the Effelsberg 100-m telescope.Our data has the thermal emission subtracted, along with 38 point sources which are not associated with M31.This is the highest frequency, and therefore highest resolution, intensity map of M31 measured by the Effelsberg telescope.The frequency bandwidth is ∆ν = 1.1 GHz, while the halfpower beam-width (HPBW) is 1.5 ′ (corresponding to a physical size of 0.34 kpc at the distance of M31).The root-mean-squared (rms) noise of the data is σ rms = 0.25 mJy/beam in the inner 9.13 kpc × 9.13 kpc region and σ rms = 0.3 mJy/beam elsewhere.
In Figure 1, we show the intensity map of the data. 2he reported intensity at each pixel is the radio emission measured by the Effelsberg telescope from that location on the sky; this corresponds to the true differential flux convolved with the frequency band and the angular beam centered on that location.Our x and y coordinates are oriented so that the x axis is aligned with the semimajor axis of M31, converting angular coordinates to lengths assuming a distance to M31 of 785 kpc [46].
We note that the observations of M31 have significant negative values, well in excess of statistical expectations given the rms noise.Most notably, the data has a large negative excursion located near the center of M31, at (x, y) ∼ (−2, 1) kpc.This may be due to over-subtraction of one of the point sources identified by Ref. [35] These negative values suggest that pixels labeled as having a flux of 0 mJy/beam actually may have a significant (unknown) positive flux.This in part motivates our choice to set limits using morphology of the expected dark matter-induced radio signal, rather than overall intensity.Like the Milky Way, M31 is a spiral galaxy with approximate axisymmetry around a rotating stellar disk.We adopt cylindrical coordinates with the origin at the center of M31, the cylindrical radius R and the height away from the midplane of the disk z.The assumption of axisymmetry implies there is no dependence on the angle around the disk ϕ.We will refer to the spherical radius from the center of M31 as r.Note that we observe M31 at an angle of inclination given by β = 77.5 • [47], and so the cylindrical (R, z, ϕ) coordinates are projected on to the x − y coordinate system of Figure 1.

III. DARK MATTER PRODUCTION OF e ± IN M31
Dark matter annihilation to unstable Standard Model particles such as b-quarks, τ leptons, or W bosons will result in cascade decays involving large numbers of leptons and QCD bound states, many of which themselves will decay into prompt photons and e ± .Though the exact spectrum of stable final-state particles that result from these decay showers of course depends on the initial Standard Model pair produced in the annihilation, there are also broad similarities regardless of the progenitors.As explicit calculation of the showers show (using a particle hadronization and decay package such as pythia8 [48]), the final state particles will have energies in the O(0.1 − 10 GeV) range for dark matter with the weak-scale masses typically expected for thermal relics.In this section, we calculate the injection morphology and spectrum of electrons and positrons produced in M31 due to dark matter annihilation, assuming weak-scale masses and annihilation into b b pairs.
While the flux on Earth of prompt photons from dark matter annihilation involves the integration of the dark matter density squared along the line of sight, electrons and positrons generated far from the Earth do not propagate to detectors here.Instead, we must track the evolution of the e ± phase space density as the particles diffuse and energy is lost -a task we will take up in Section V.For now, we will quantify the rate of production of the e ± with a source term, which depends on the local dark matter density, ρ χ (x), at every location within M31 and the particle physics model of the dark matter candidate.The source term or injection density rate of e ± due to dark matter self-annihilation is given by where ⟨σv⟩ is the thermally averaged annihilation cross section, m χ is the dark matter particle mass, dN e /dE is the injection spectrum of e ± per annihilation in terms of the total energy 3 E = (m 2 e + p 2 ) 1/2 of the e ± .Here and throughout this work, we have assumed that dark matter is its own antiparticle.If it is not and the abundance of dark matter particles and anti-particles is symmetric, there is an additional factor of 1/2 in Eq. (1).
The energy spectrum of the e ± source is determined by dN e /dE which is influenced by the dark matter mass and the annihilation channel.Our choices for these parameters are motivated by fits of dark matter annihilation to the gamma ray excess in the Milky Way's Galactic center [16][17][18][19].In the Milky Way, the dark matter candidates that best fit the gamma ray excesses have m χ ∼ 30 − 50 GeV annihilating to b b or m χ ∼ 10 GeV annihilating to τ + τ − [4,[16][17][18][19].A similar signal in M31 has a best fit of dark matter with mass m χ ∼ 10 GeV annihilating to b b, or m χ ∼ 5 GeV annihilating to b b and τ + τ − democratically [20].
In part motivated by these fits, in this work we consider dark matter annihilation into b b with m χ in the range 6 − 500 GeV.We scan in dark matter mass from 6 − 500 GeV, and use pythia8 to decay, shower, and hadronize the b b annihilations to calculate dN e /dE for each choice of m χ . 4We show in Figure 2 the resulting e ± spectra for a sample of representative dark matter masses.All spectra are cut off at low energy by E = m e .
The morphology of the source is determined by the dark matter distribution of M31.This can be fit by a modified Navarro-Frenk-White (NFW) profile [49][50][51] where γ NFW is the logarithmic inner-slope, ρ 0 is the scale density and r s is the scale radius.In the Milky Way, the gamma ray excess favors an inner slope of γ NFW ∼ 1.25 [16][17][18][19].This slope is adopted by Ref. [20] for the dark matter distribution of M31.However, there is considerable uncertainty in the M31 dark matter distribution.
In keeping with available kinematic fits to the rotation curve, we use a standard NFW with γ NFW = 1.For the scale density and the scale radius, we use the best-fit values from analysis of kinematic data [43]:

IV. ASTROPHYSICAL MODEL OF M31
The intensity and morphology of a radio signal which originates from electrons and positrons in a galaxy depends greatly on how the charged particles propagate.In M31, the most important propagation effects are diffusion and energy loss [20].To calculate the effects of 4 The pythia shower was modified to allow decays of Standard Model particles which are meta-stable on detector timescales, but whose decays would be astrophysically relevant, e.g., µ, π 0 , π ± , and neutrons.these, we must first model the properties of the interstellar medium.
In this section, we present our models for the magnetic field, the interstellar radiation field (ISRF), and the various components of thermal gas within M31.As relativistic e ± propagate, the interstellar magnetic field causes them undergo synchrotron energy loss, emitting photons at radio frequencies.Random fluctuations in the magnetic field also diffuse the charged particles through space.The ISRF causes the e ± to lose energy through inverse-Compton scattering.Lastly, the various components of gas cause energy loss through bremsstrahlung and Coulomb scattering.
The measured distance to M31 has varied considerably in the literature.Early measurements found a value of 690 kpc [52], but modern techniques and measurements prefer 760 − 785 kpc [46,53,54].As a result, the references we used to construct our astrophysical models assume a variety of distances to the galaxy.Throughout this work, we use a distance to M31 of 785 kpc [46,53] obtained using measurements of tip of the red giant branch stars.When necessary, we scale the results of previous works used in our model of M31.

A. Magnetic Fields of M31
The magnetic field B in M31 is turbulent, with fluctuations on many length-scales, ranging from the size of the galaxy down to far below the resolution limit of our experimental probes.The details of small-scale fluctuations are not observationally accessible, but the expectation value of B 2 and the power spectrum of fluctuations of the magnetic field at different locations in the galaxy can be measured though radio polarization and cosmic ray propagation observations, respectively.These observationally accessible features of the galactic magnetic field suffice for our purposes.
We write the magnetic field as the product of an axisymmetric field magnitude and a random dimensionless vector-field that depends on location x: The vector b contains the local fluctuations in the field, and is the expectation value of the magnitude of the field squared, which we assume is independent of ϕ.Formally, this is an ensemble average with respect to the probability distribution that the magnetic field is sampled from.
As is conventional [20,[55][56][57][58], we characterize the magnetic field fluctuations in terms of a power spectrum normalized as where k 0 is the minimum wavenumber for which the power spectrum applies.The length-scale 1/k 0 is typically assumed to be a factor of O(10 − 100) smaller than of the characteristic length-scale of the galaxy [55].For M31, this implies O(0.1 kpc) ≲ 1/k 0 ≲ O(1 kpc) (though when setting limits we will vary this parameter over a much more conservative range).Observations of the propagation of cosmic rays in the Galaxy [59] find a diffusion coefficient of the form D ∝ E δ with δ ≃ 1/3.This is consistent with magnetic field fluctuations that follow a Kolmogrov spectrum [60], P b ∝ k −5/3 .The diffusion of charged particles moving in a magnetic field with fluctuations obeying a Kolmogrov spectrum will be discussed in detail in Section V A. We determine the R dependence of B from measurements of the M31 magnetic field (taken to be the rootmean-squared (RMS) field strength) in the disk in three regions: within the inner 1 kpc [61], in the range 6−14 kpc [62], and in intergalactic space [63].The measurements from Refs.[61,62] are shown in Figure 3.The intergalactic magnetic field has been measured to be at most 0.3 µG [63], which we take to be the 1σ upper bound of the field strength outside of M31.To require our fit to the M31 field strength to agree with this upper bound, we include B = 0.15 ± 0.15 µG at R = 300 kpc (though it is not shown in Figure 3) in addition to the measurements from Refs.[61,62].
We fit the RMS magnetic field measurements to the functional form which has sufficient flexibility to fit the available data.Our best-fit parameters of Eq. ( 7) are shown in Table I.
Assuming equipartition between the cosmic ray and magnetic field energy density, the vertical scale height of the magnetic field is approximately four times the scale height of the disk of synchrotron emission [64], which itself is approximated by the scale height of the HI disk [62].The HI disk scale height as a function of R was measured by Ref. [65]; rescaling from the distance to M31 assumed in that work (690 kpc), the magnetic field scale height is h B (R) =4h syn = 4h HI =(0.83 ± 0.17 kpc) + (0.064 ± 0.012) × R.  I) is shown in blue.

B. Interstellar Radiation Fields of M31
The ISRF of M31 has contributions from the CMB, starlight, and infrared emission from dust.We model each component, and sum the results to obtain the total energy density of radiation within M31.Most straightforward is the energy density of the CMB, which is (for our purposes) uniform and given by where k B T CMB = 2.3 × 10 −4 eV [66].We determine the energy density from stars and dust from the measured luminosity distribution of M31.The energy density of stellar radiation ρ * in M31 is related to the bolometric luminosity density Q * by In Ref. [67], the extinction-corrected stellar luminosity distribution of M31 was modeled for five different structural components (bulge, disk, nucleus, young disk and stellar halo).As these distributions are extinction corrected, they describe the luminosities that would be observed without dust absorbing stellar emission.Assuming that the dust is in equilibrium with the starlight, the luminosity it emits in IR is equal to the luminosity that it absorbs.Therefore, the extinction-corrected stellar luminosity distribution integrated over all wavelengths approximates the bolometric luminosity from stars and dust combined.
Ref. [67] models the luminosity density Q j as where j = (bulge, disk, nucleus, young disk, stellar halo) indexes the various components of the luminosity distribution and (Q 0,j , q j , A 0,j , N j ) are parameters determined separately for each component.As Ref. [67] finds that the M31 luminosity is dominated by the bulge and disk components, we only include those in our model of the ISRF.
Ref. [67] fits the parameters q j , A 0,j , N j to data, and provides the extinction corrected total luminosity of each M31 structural component in each of the ugriz filter bands.These luminosities are defined as where a = (u, g, r, i, z) represents the spectroscopic band of the measurement, and λ a is the central wavelength of the relevant band.The Q 0,j depend on the total bolometric luminosity of each component L bol,j : The bolometric luminosity of the j component can be found by integrating Therefore, we need the luminosity per unit wavelengthor spectral energy distribution (SED) -of each component.
The SED is only available in the inner 1 kpc of M31 [68], but M31's bulge is concentrated in the inner 1 − 2 kpc [67].Therefore, we use the SED from the inner 1 kpc of M31 [68] as a template for the SED of the bulge, renormalized so that it accounts for the luminosity of the whole bulge.To chose the appropriate normalization we digitize and smoothly interpolate the extinction corrected SED of the inner 1 kpc of M31 from Figure 1 of Ref. [68], and fit the proportionality constant between the inner kpc and the entire bulge by minimizing the χ 2 between the re-normalized SED and the bulge luminosity in each band, derived from the fits of Ref. [67].Our best fit rescaled SED and the measured luminosities in the ugriz bands for the bulge are shown in Figure 4(a).Using the best fit rescaled SED, we then integrate dL bulge /dλ over wavelengths λ to obtain the bolometric luminosity of the bulge.
For the functional form of the SED for the disk, we subtract our best-fit SED for the bulge from the extinction corrected SED for the whole galaxy (given in Ref. [69]).Using this functional form, we repeat the procedure that we used for the bulge: we perform a χ 2 minimization to find the proportionality constant that gives the best agreement between the rescaled SED and the luminosity values for the disk, and integrate dL disk /dλ over λ to get the bolometric luminosity of the disk.Our best fit rescaled SED and the measured luminosities in the ugriz bands for the disk are shown in Figure 4(b).
We use our bolometric luminosity values and Eq. ( 13) to calculate Q 0,j for each component.We show the values of the structural parameters (Q 0,j , q j , A 0,j , N j ) and luminosities (L u,j , L g,j , L r,j , L i,j , L z,j , L bol,j ) in Table II (for j = bulge, disk).We add our results for the luminosity density of the disk and bulge to get Q * (x) and numerically integrate Eq. ( 10) to obtain ρ * .As we derived our luminosity density distributions from extinction corrected luminosities [67] and SEDs [68,69], our result for ρ * contains contributions from starlight and dust.Our model for the ISRF ρ tot = ρ * + ρ CMB as well as the radiation density from individual components in the plane of the disk are given in Figure 5.
It is important to note that our model for the starlight luminosity of the innermost regions of M31 is significantly larger than the equivalent for the Milky Way [70] (shown with the dashed curve in Figure 5).Previous dark matter studies of radio emission from the center of M31 used a starlight model scaled with an O(1) factor from the Milky Way [20,36,38,39].The higher luminosity in the center of the galaxy that we find in our model leads to greater energy losses into X-rays from e ± inverse Compton scattering with the starlight photons.As a consequence of this increased energy-loss mechanism, the radio signature of dark matter-produced electrons and positrons in the galactic center is reduced.
starlight model more closely matches those previously assumed for M31.As we will describe in detail, our dark matter constraints are obtained from radio emission in this region rather then from the center itself.Consequentially, we are less sensitive to differences in the starlight in the core of M31.  9).The bulge and disk components come from replacing Q * with Q bulge and Q disk , respectively in Eq. (10).The Milky Way ISRF is digitized from Figure 1 of Ref. [70].11).Bottom: observed extinction-corrected luminosities in ugriz filter bands followed by our derived bolometric luminosities.The bulge values are in the second column while the disk values are in the third column.All values except for L bol and Q0 are taken from Ref. [67], see text for details of our calculations of L bol and Q0.

C. Gas in M31
The thermal gas of M31 can be split into ionized gas and neutral gas, each of which play different roles in the energy loss of relativistic electrons and positrons.Elastic collisions between the ionized gas and the e ± result in Coulomb losses in the e ± .This leads to a net transfer of energy out of the e ± and into the gas.Interactions between e ± and ionized gas also result in bremsstrahlung losses due to inelastic collisions.Neutral gas only causes bremsstrahlung losses.As the rate of energy loss depends on the properties of the ionized and neutral gas, we must model HI, H 2 , and 4 He gas separately.

Ionized Gas
Due to the difficulty of observing M31 as compared to our own Galaxy, the gas model of M31 is motivated by that of the Milky Way.Following Ref. [71], we model the ionized gas density as where ⟨n ion ⟩ is the ion density averaged over scales that are small compared to the galaxy but large compared to fluctuations in the ion density.The parameters nion , R ion , and z ion will be fit to M31 measurements of the ion density from Hα emission [72] and Faraday rotation [62], as we will discuss below.
We first extract the ion density at the mid-plane of M31 at R ≃ 9 kpc from measurements of Hα emission [72].Here, the observable is the emission measure (E) along the line of sight, which is related to the ion density by Ref. [72] measures the value of E along a line of sight at R = 9 kpc and then projects to the result that would be observed if M31 were viewed face-on, finding E = 6 − 15 pc cm −6 .The ionized gas inhabits the galaxy in clumps that are small compared to astrophysical scales.These clumps make up a fraction of the volume of the galaxy given by the fill factor ϕ, which Ref. [72] assumes to be the Milky Way value of 0.2.Under this assumption (and taking the maximal value of E = 15 pc cm −6 as well as their median scale height of ẑion = 500 pc), Ref. [72] calculates a midplane ion density of nion = 0.39 cm −3 [72].This is the density in a gas clump corresponding to a density field in the neighborhood of R = 9 kpc, z = 0 kpc of where F is the filling function which varies over short length-scales and is defined to be 1 inside a clump of ionized gas and 0 outside.Its spatial average over astrophysical length-scales is equal to the fill factor ϕ. The result obtained by Ref. [72] for nion is likely an overestimate since they obtained it using their maximal value of E. Therefore, we modify it to the result they would have obtained if they had used their central value of E. Given Eqs. ( 16)&( 17), the mid-plane ion density scales with the emission measure and the scale height as nion ∝ E ẑion .
We can use this dependence to re-scale the results of Ref. [72] for nion , as well as propagate uncertainties based on the measured value of E ≃ (10 ± 5) pc cm −6 and the range of assumed values for the scale height ẑion ≃ (500 ± 300) pc.After doing so, our re-scaled result for the mid-plane ionized gas density becomes nion = (0.32 ± 0.20) cm −3 .Using Eq. ( 17), and averaging over a small neighborhood around R = 9 kpc, z = 0 kpc leads to our inferred value of nion in this neighborhood: By observing rotation measures (RM) from Faraday rotation and assuming magnetic field equipartition, Ref. [62] determines nion in the upper layers of the thermal disk (between 0.3 − 1 kpc from the mid-plane) at three different values of R between 8 − 14 kpc.We take the distance of these measurements from the midplane to be the midpoint of the upper layers of the thermal disk (z = 0.65 kpc).Errors were not reported for these results; we make the conservative choice to use errors of 50% of the measured value.TABLE III.Values of the ionized gas density derived from observations of Hα emission [72] and Faraday rotation [62].These derived values are used to fit our model of the ionized gas density given in Eq. ( 15).
We summarize the ionized gas observations from Faraday rotation [62] and our derived value at the midplane from the measurement of emission measure [72] in Table III.Using these observations and our derived value, we perform a χ 2 fit to Eq. ( 15); we show our results in Table IV.

HI Gas
For the HI gas density, we digitize the radial column density distribution provided by Ref. [73] in the range R ∈ [1,23] kpc.We then interpolate and extrapolate the digitized column density to get a function valid for all R ≥ 0. For the extrapolation, we use a growing exponential at small R and a decaying exponential at large R, such that the value and slope of the function are continuous at each boundary.Our interpolation and extrapolation of the fit to the column density of HI within the

Gas Density Parameter
Value Ionized nion,0 (cm −3 ) 0. The next two panels list parameters for the HI and H2 distributions along the z coordinate given by Eqs. ( 20) and ( 22), respectively.
disk of M31 [73] is shown in Figure 6(a) (see the vertical axis on the right).
We assume that the gas density on the disk is proportional to the column density.For the z-dependence of the HI density we assume a decaying exponential with a scale-height of h HI (R) taken from Ref. [65], as previously discussed in Section IV A: where h HI,0 and S are reported in Table IV.We normalize the HI density to the total HI mass reported in Table 1 of [73].The resulting gas density on the disk is shown in Figure 6(a) (see the vertical axis on the left).

H2 Gas
For the H 2 density, we digitize the radial column density distribution provided by Ref. [73] in the range R ∈ [2.5,18] kpc.We then repeat our interpolation/extrapolation procedure to obtain a model of the column density of H 2 on the disk.The results are shown in Figure 6(b).We again assume that the density of H 2 on the disk is proportional to the column density.
We assume a decaying exponential for the zdependence, however there is little observational data available for the scale height of H 2 in M31.We therefore use the H 2 scale height in the Milky Way (derived from the data in Ref. [74]), and re-scale to M31 using a comparison of the HI scale heights in the Milky Way and M31.We digitize the fit to the H 2 scale height data in Figure 10 of Ref. [74], and smooth-out the fluctuations by fitting the digitized version of the fit to the functional form We compare our model for the scale height in the Milky Way to that of Ref. [74] in Figure 7.
To obtain a scale height for H 2 gas in M31 from the scale height of H 2 in the Milky Way, we take the ratio of the average HI scale height for M31 (Eq.( 20)) to that of the Milky Way (see Figure 6 of Ref. [74]) in the region R ∈ [0, 7] kpc.We find this ratio is 1.55, and assume this R (kpc) The parameterized fit of Eq. ( 21) is shown in blue.
ratio holds for the H 2 gas: The resulting values of h H2,0 and R H2 are given in Table IV.The errors in these two parameters are dominated by the systematic errors of converting from their values in the Milky Way so we conservatively set their errors to 50% of their values.Again, we normalize our distribution for the H 2 density based on the measured value for the total H 2 gas mass reported in Table 1 of Ref. [73].The resulting H 2 density on the disk is shown in Figure 6(b).

4 He Gas
The last significant component of neutral gas that we have to model is 4 He.The mass fraction of 4 He in simulated spiral galaxies was found to vary spatially over the range Y He ≃ 0.25 − 0.3 [75].The lower limit is the primordial value set by Big Bang Nucleosynthesis, and higher values are due to stellar production.For simplicity, we approximate that 4 He production from stars is negligible, so the ratio of 4 He to H is set by Big Bang Nucleosynthesis: The errors in our model for the Helium density from using this approximation are < 25%.As the effects of bremsstrahlung from 4 He are subdominant compared to the other components in our propagation model, the errors in our final results from this approximation will be negligible as well.
Further, we assume that the 4 He density has the same morphology as the total hydrogen density.The local density for 4 He is then derived from the HI and H 2 gas: The production of electrons and positrons by dark matter annihilation provides a source Q e for the phase spacedensity f e within a galaxy.From their initial locations, the e ± will diffuse in turbulent magnetic fields and undergo energy loss from synchrotron radiation as well as inverse Compton, bremsstrahlung and Coulomb scattering.The synchrotron losses lead to the radio signal we will use to constrain dark matter annihilation, but all forms of energy loss must be tracked to determine the evolution of f e in energy and position.
The evolution of the phase space density f e is controlled by the diffusion-loss equation: where f e (x, E) = dn e /dE is the phase space density of electrons at position x and energy E, D ij (x, E) is the diffusion matrix, and b(x, E) is the energy loss parameter.The position-dependent diffusion matrix depends on both the RMS magnetic field, B and the turbulent fluctuations of the field at small scales.The loss parameter depends on B2 , the ISRF, and the densities of the various gas components, which have been modeled in Section IV.

A. Diffusion Matrix
Fluctuations in the magnetic field cause the relativistic e ± to exhibit diffusive motion.In a uniform magnetic field with strength B, the motion of e ± is helical, with a Larmor radius of (26) Here, E is the particle energy and α is the pitch angle, which is defined as the angle of the velocity with respect to the direction of the magnetic field.For a changing magnetic field, the Larmor radius can still be defined, provided the field variations are small over the distance traversed by the particle in the time it takes for its phase to change by O(1).For the O(10 µG) magnetic fields in the disk of M31 (see Section IV A), electrons and positrons of the energies expected from the annihilation of dark matter have Larmor radii of ≲ O(10 −7 pc).As the fluctuations of the field follow a power law distribution on length-scales smaller than 1/k 0 ∼ O(1 kpc), the magnetic field is dominated by Fourier modes that are much larger than the Larmor radius.Modes that are of order the Larmor radius and smaller can be treated perturbatively.
The motion of an electron or positron under the influence of the large-scale magnetic field fluctuations is well-described by the adiabatic approximation: the particle exhibits helical motion about the local field with an axis that gradually changes as the particle moves along the slowly changing magnetic field [76].The magnetic field fluctuations with length-scales below the Larmor ra-dius perturb this adiabatic motion, causing pitch angle scattering which leads to diffusion along the axis of the local magnetic field [77].Since the direction of the largescale magnetic field slowly changes over space and time, the diffusion is shared evenly in all directions leading to isotropic diffusion [77].
With these approximations and assuming magnetic field fluctuations characterized by the Kolmogrov spectrum, introduced in Section IV A, the diffusion matrix is [55,78] where d 0 = 1/k 0 is largest length-scale over which the Kolmogrov spectrum of magnetic field fluctuations is valid.Since the diffusion matrix is isotropic, it can be written as It is conventional to refer to D the diffusion coefficient.We absorb the uncertainties in the prefactor of Eq. ( 27) and d 0 into the constant D 0 .A range of possible values for D 0 (in both the Milky Way and M31) have been suggested in the literature.We review these briefly here.
Ref. [79] infers the diffusion coefficient for e ± near star-forming regions in M31 from measurements of nonthermal radio emission at ν = 1.4 GHz and one higher frequency using two methods.The first method infers the diffusion coefficient from the difference in morphology between the two frequencies.The second method uses the difference between non-thermal emission and thermal emission at each of the frequencies, assuming that the thermal emission has a similar morphology to the source distribution of cosmic ray electrons.These methods allow Ref. [79] to extract the diffusion coefficient at two electron and positron energies: 4.1 and 7.5 GeV.
Rescaling Ref. [79]'s results to the magnetic field parameters of M31, we obtain D 0 ≃ 1.1 × 10 28 cm 2 /s using the first method and D 0 ≃ 3.5 × 10 27 cm 2 /s using the second.Neither method fully models the propagation of cosmic rays, and the variation between the two results makes it difficult to identify either value as our default D 0 value.
For further guidance about the value of D 0 in M31, we review studies of propagation in the Milky Way.The galprop cosmic ray propagation model [56] uses observations of cosmic rays in the Milky Way [80,81] to determine its best-fit diffusion coefficients [82].The assumed propagation model includes a uniform diffusion coefficient of the form D = D0 (E/4 GeV) δ for which the best-fit parameters were found to be D0 = (8.3± 1.5) × 10 28 cm 2 /s and δ = 0.31 ± 0.02 [82].Comparing the model of Ref. [82] to our form for the diffusion co-efficient, Eq. ( 27), and assuming that the cosmic rays studied were subject to a constant 10 µG magnetic field, these measurements imply D 0 = (5.2± 0.9) × 10 28 cm 2 /s.
Ref. [83] constructed MIN, MED and MAX propagation models for 1 GeV cosmic ray energies in the Milky Way, which can be interpreted as a range of D 0 = 5.9 × 10 27 − 2.0 × 10 28 cm 2 /s, assuming the relevant magnetic field is a constant 10 µG.
Underestimating the diffusion coefficient would underpredict how far particles will move before emitting most of their energy, leading to an over-prediction of the signal in regions of high dark matter density.A large value of D 0 will likewise result in a smaller signal flux for a given cross section.Therefore, to set conservative limits on the cross section of dark matter annihilation, we must avoid assuming too small a value for D 0 .To set our conservative upper limit on D 0 , we select a maximum value of d 0 ≡ 1/k 0 in Eq. ( 27) by setting k 0 equal to the wavenumber of a fluctuation with wavelength of R B,1 = 77.6 kpc (the longest scale-length in our magnetic field model).This leads to d 0 ≲ R B,1 /(2π) ≃ 12.5 kpc, implying D 0 ≲ 8 × 10 28 cm 2 /s.
Given the variation in diffusion coefficients in the Milky Way and M31, as well as our conservative upper bound, we consider D 0 in the range We select as a default value D 0 = 1 × 10 28 cm 2 /s.Though D 0 is position-independent, the diffusion coefficient D explicitly depends on the magnetic field, and our model for the magnetic field is position dependent (see Section II).As a result, the diffusion coefficient also depends on location within M31.We show the dependence on location within M31 in Figure 8 for our default diffusion coefficient normalization (D 0 = 1 × 10 28 cm 2 /s) and E = 1 GeV.The diffusion coefficient varies more rapidly with z when R is small, as a result of the magnetic field scale height increasing with R.
Prior studies of radio emission from dark matter annihilation in M31 make the approximation that the dif- fusion coefficient is zero [36][37][38] or position independent [20,39].This latter assumption is sufficient when the region of interest is small and the diffusion coefficient is nearly constant over the region.Over the length scales of interest, the variations in the diffusion coefficient must be taken into account, and we develop a numerical method for calculating the evolution of phase space density of the charged particles which can accommodate a positiondependent diffusion coefficient.We describe our numerical solution in Section V C, after introducing the energyloss terms that enter into Eq.( 25) in the next subsection.To our knowledge this is the first study that uses a position-dependent diffusion coefficient to set limits on dark matter annihilation in M31 via radio emission.

B. Energy Loss due to Radiative Processes
As the electrons and positrons diffuse through the ISM they lose energy through radiative processes.This energy loss is encoded in Eq. ( 25) by the loss parameter b.In M31, the relevant losses are inverse Compton (IC), synchrotron, bremsstrahlung, and Coloumb interactions: We treat each of these terms in turn.
The inverse Compton scattering between e ± and the ambient starlight, rescattered light from dust, and CMB emission will transfer energy from the charged particles into the photons, at a rate [76] where b (0) IC = 1.0 × 10 −15 GeV/s and ρ γ is the total radiation energy density, derived in Section IV B.
Synchrotron emission occurs due to the acceleration of charged particles in galactic magnetic fields.As described in Section IV A, the magnetic fields of M31 do not change appreciably over the Larmor radius of the relevant e ± .Additionally, due to pitch angle scattering, the pitch angles are approximately uniformly occupied.Therefore, the energy loss due to synchrotron emission can be determined by assuming a locally constant magnetic field and averaging the energy loss over all pitch angles.The expression for the loss due to synchrotron is given by [76] where b (0) sync = 2.5 × 10 −16 GeV/s.The third term in Eq. ( 30) is the contribution to the loss from bremsstrahlung emission due to e ± scattering with neutral hydrogen, neutral helium, and ionized gas: The expressions for these three components of the bremsstrahlung loss are given by [70,84] b , where b Lastly, the fourth contribution to the loss parameter is from Coulomb interactions with ionized gas and is given by [85,86] where b (0) C = 6.2 × 10 −16 GeV/s.Note that Coulomb losses are not radiative processes involving the loss of energy from the charged particles into photons, but rather are due to an energy transfer from the relativistic e ± to non-relativistic ions in the interstellar plasma.
We show the resulting loss coefficient as a function of energy in Figure 9. Figure 9(a) shows the total loss coefficient given by Eq. (30) for various values of R on the disk.

C. Solving the Diffusion Loss Equation
We now turn to the numerical solution to the diffusion loss equation (Eq.( 25)) in M31, assuming the electron and positron injection from dark matter from Section III and the astrophysical model of M31 from Section IV.
To motivate our approach, it is useful to consider the two dynamic time scales which characterize the diffusion (τ D ) and energy loss (τ b ), defined implicitly through rewriting Eq. (25) as These timescales depend on R, z, E and derivatives of f e over f e .As a result, τ b and τ D are independent of the overall magnitude of f e .The timescale for diffusion is In the approximation that f e depends on position only as a power-law in r, where L(x) is a length-scale that determines the rate that diffusion causes the phase space density to change.
Assuming f e is a power law in E, the characteristic timescale for energy loss can be approximated by The propagation is dominated by diffusion when τ D ≪ τ b and dominated by energy loss when τ D ≫ τ b .Both diffusion and loss will dominate at different values of E and x.In Figure 10 we plot the inverse timescales for diffusion and loss over a range of R and z.Along the disk, loss tends to dominate (Figure 10 FIG. 10.The inverse timescales for diffusion (solid lines) and loss (dashed lines), Eqs. ( 38) and (39).We show for comparison the inverse of the age of M31 (solid black), TM31 = 10 10 years.The shaded regions around each solid line shows the variation of the inverse timescales for diffusion as D0 is varied within the range given in Eq. ( 29).As can be seen in Figure 11, within R < 25 kpc and |z| < 15 kpc, we find τ < T M31 .Near the center of the galaxy and for higher e ± energies, τ decreases.Though some regions at large R have timescales comparable to the age of M31, these regions are far from the inner part of the galaxy where the Effelsberg radio data will be used to set limits.We are therefore justified in following the general approach of the literature [20,[36][37][38] by approximating the phase space density f e of e ± in M31 today as the equilibrium density.
If b and D do not depend on x, a semi-analytic solution exists for the equilibrium density (see e.g., Ref. [45]).When the region of interest is small, homogeneous coefficients can be obtained by averaging the diffusion and loss coefficients over the relevant volume [20,45].However, our goal in this paper is to compute the synchrotron distribution over the field of view of the radio data in Figure 1, that is, most of the galactic disk of M31.Based on the astrophysical models (described in Section IV), the diffusion and loss coefficients will vary significantly over this region.We must therefore solve Eq. ( 25) in the case of non-homogeneous coefficients.
While the source term is spherically symmetric, the diffusion and loss coefficients are axially symmetric, implying that the solution to Eq. ( 25) depends on R, z and E. However, a fully axially symmetric numeric solution is intractable given our numeric approach.To overcome this problem, we average Eq.( 25) over solid angle Ω: where (for an arbitrary function g(E, x)), Spherically averaging Eq. ( 40), we find where The averaged coefficients D and b required to solve for ⟨f e ⟩ in Eq. ( 42) themselves depend on f e .To calculate D and b, we use approximate solutions for f e , then use these averaged coefficients to numerically solve the spherically averaged diffusion loss equation for ⟨f e ⟩.
We calculate these approximate solutions for D and b in two different ways: first by assuming that f e is approximately spherically symmetric, and in the second approach taking into account approximate deviations from spherical symmetry.For the region of M31 of interest for our analysis of the radio data (namely, the region inside ∼ 10 kpc), the resulting solutions for ⟨f e ⟩ (and the resulting synchrotron emission) are similar regardless of our assumptions.
Our first approximate solution -the "unweighted" solution -assumes that deviations from spherical symmetry for f e are small.If this is the case, and we can numerically average over solid angles our models for D and b given in Sections V A and V B.
For our second solution for the diffusion and loss coefficients, we calculate the spherically-averaged D and b parameters by substituting into Eq.( 43) the approximate equilibrium solution for f e obtained from solving Eq. (36) with ∂f e /∂t = 0: where We refer to this as the "weighted" solution, as D and b are obtained as averages of D and b weighted by ∂ r f e and f e , respectively.
Figure 12 shows the spherically averaged diffusion and loss coefficients for the unweighted (dashed) and weighted (solid) averaging schemes.The two methods give very similar results for the loss coefficient across all of phase space and galactic radii.For the diffusion coefficient, the two calculations agree for the inner part of M31, r ≲ 10 kpc.As we will show, the disagreement at large radii in the diffusion coefficients does not result in significant differences in the predicted synchrotron emission from the region of M31 that we will use to set limits.As a result, the constraints derived from radio observations are robust across these different solutions.
Using either the weighted or unweighted solutions for D and b, we must solve the diffusion loss equation for ⟨f e ⟩.Defining u ≡ r⟨f e ⟩, Eq. ( 42) becomes ) Under this redefinition, the boundary condition at r = 0 can be easily written as u(0, E) = 0, as long as f e does not diverge faster than 1/r.This is satisfied if the inner slope of the dark matter density has a power law index γ NFW < 1.5.The other required boundary conditions are u(49.9kpc, E) = 0 and u(r, m χ ) = 0. To solve Eq. ( 46), we discretize r, E and t and use finite differences to approximate the derivatives.This leads to a recursive equation for u at the next time-step in t, given its value at the current t.
Forward difference schemes for solving Eq. ( 46) are only stable if the time-step satisfies ∆t ≲ (∆r) 2 /D over the whole domain [87], where ∆r is the grid-spacing.Given the approximate age of M31 and our grid-spacing ∆r = 62 pc, O(10 7 ) time-steps would be needed to reach the equilibrium solution using a forward difference FIG.
12. Spherically averaged diffusion coefficient (left) and loss coefficient (right) using the unweighted average (dashed) and the weighted average (solid) for a range of energies.We use our default value of D0 = 1 × 10 28 cm 2 /s.method.Backward differences, on the other hand, are unconditionally stable [87] for any size of time-step.We therefore use backward differences to approximate the derivatives on the right-hand-side, leading to an implicit equation for u at the next time-step, which can be solved with a sparse matrix method.We choose the time-step to be much larger than the maximum timescale in the problem to minimize the number of iterations required.Further details about our numerical method for solving the diffusion-loss equation are provided in Appendix A.
The results for the equilibrium solutions of ⟨f e ⟩ are shown in Figure 13. Figure 13(a) shows the energy dependence of ⟨f e ⟩ at r = 5 kpc for a representative set of values of m χ and our default value of D 0 .Figure 13(b) shows the dependence on r for each value of D 0 and E = 1 GeV.The results become more sensitive to changes in D 0 for D 0 ≳ 3×10 28 cm 2 /s.In both panels, the solid curves represent the weighted solution while the dashed curves represent the unweighted solution.

VI. SYNCHROTRON SPECTRUM AND MORPHOLOGY
Relativistic electrons and positrons in M31 accelerate in the galactic magnetic field, leading to synchrotron emission.The power emitted per unit frequency from an electron or positron at pitch angle α and energy E is [88] where ν 0 ≡ e B/(2πγm e ), x = 2ν/(3γ 3 ν 0 ), γ is the Lorentz factor and K n is the n th -order modified Bessel function of the second kind.The differential flux can therefore be obtained by averaging Eq. ( 47) over uniformly distributed pitch angles, and convolving with the spherically-averaged phase space density of electrons leading to: where Ω = (θ, ϕ) is the location on the sky.For sin α ∼ O(1) and x ≫ 1, dP/dν is exponentially suppressed at low energies [88], so most of the power is radiated by e ± with energies satisfying As the Effelsberg radio telescope data used in this study is at frequencies around 8.35 GHz, we are most interested in the e ± produced through dark matter annihilation with energies of ∼ 10 GeV and higher.This is shown in Figure 14 where we plot the dependence of ⟨dP/dν⟩ α on E for a variety of fixed values of B.
In Figure 15, we show the 8.35 GHz radio emission resulting from dark matter of mass m χ = 39 GeV annihilating with a cross section of ⟨σv⟩ = 2.2 × 10 −25 cm 3 /s, assuming D 0 = 1 × 10 28 cm 2 /s.In Figure 16, we show the signal along the semi-major axis (a) for a variety of values of m χ , holding D 0 constant and (b) for a variety of values of D 0 holding m χ constant.

VII. STATISTICAL METHODOLOGY
Having developed a numeric method to calculate the radio emission induced by dark matter annihilation in M31, we can now compare our predicted signal with data to set limits on the annihilation cross section to bb as a function of dark matter mass.Though the dark matter annihilation will be brightest in the center of the galaxy, this region also has significant baryonic sources whose intensities cannot be easily modelled.In addition, the flux near the center of M31 is sensitive to the value of D 0 for D 0 ≳ 1 × 10 28 cm 2 /s (as seen in Figure 16(b)).For these reasons we set limits using the expected morphology of the dark matter signal outside of the center, rather than the total intensity.This also makes our constraints insensitive to possible mismeasurement of the overall zero-level of the radio data.
This approach requires data-driven modeling of the backgrounds within the galaxy.The background emission in M31 is complicated, with numerous point sources and a prominent ring feature (see Figure 1).None of these features are morphologically consistent with the expectations of dark matter annihilation, and can safely be attributed to baryonic physics.Even so, a multi-step process is required to define a search region and construct a background model within that region that does not risk fitting-away any potential signal.
In Section VII A, we first describe how we mask the point sources, the ring of radio emission in the disk, and the bright center of M31.This will allow us to define a search region interior to the ring, where the background can be approximated as the residual emission from the ring plus a constant.In the end, the radio emission from this search region will be used to set limits on the dark matter model.Next in Section VII B, we describe how we determine the background model within the search region using the data itself -without absorbing dark matter emission (potentially present in the data) into the model.We introduce a background model with five free parameters: three morphological parameters (µ 1 , µ 2 , µ 3 ) controlling the shape of the residual background from the elliptical ring, and two coefficients (w 1 , w 2 ) which determine the intensity of each component of the background.The morphological parameters are fixed based on the data independent of the signal hypothesis, while the intensity coefficients are adjusted to their most likely values for each hypothesis.
Fixing the morphological parameters must be done carefully to avoid absorbing any signal present in the data into the background model.We leverage the fact that the signal peaks toward the center of M31, while the emission from the ring is dominant away from the center.We therefore can fix the morphological parameters by using the data away from the center of the intensity map (exterior to the dark matter-rich "signal region").The size of this signal region is determined by comparing fits of the morphological parameters of the background model assuming the presence or absence of a dark matter signal.
After defining the search region and fixing the morphology of our background model within this region, we set statistical limits on specific signal models.We use a CL s test, described in Section VII C. CL s works by building distributions of test-statistics from synthetic observations generated from background-only and signal-plus background hypotheses.This test-statistic is sensitive to the morphology of the signal in addition to the amplitude, making this ideal for the distributed signal of dark matter in M31.Our full set of limits varying over astrophysical model parameters and using the methodology we describe here will be shown in Section VIII.

A. Background Masks
The baryonic sources of radio emission in M31 are complicated and difficult to model from first principles.Overall, we expect relatively uniform background emission across the interior of the galaxy, overlaid with significant emission from the galactic center due to baryonic processes, as well as point sources throughout the galaxy.In addition, M31 contains a prominent elliptical ring-shaped structure in radio with a semi-major axis of approximately 10 kpc, due to significant star formation in this region [89,90].All of these features can clearly be seen in the radio map of Figure 1.The location of the ring correlates with the highest densities of gas in our astrophysical models from Section IV.
Notably, other than the emission at the center of the galaxy, the spatial distribution of all of these sources of radio emission is inconsistent with emission sourced by dark matter.Rather than attempting to model these baryonic sources from first principles, we mask and remove them from our statistical analysis.As the emission at the center and the point sources are localized, we are able to completely remove them using masks.The ring of bright emission is broad enough that it cannot be removed completely.Instead, we model it as a Gaussian ring and mask its brightest emission.

Point Source Masks
Ref. [35] has removed 38 point sources unrelated to M31.However, many point sources within the galaxy remain in the data.We locate point sources algorithmically by identifying circular regions (with a diameter of 0.75 times the HPBW of the beam) that are over-bright compared to the concentric annulus with inner and outer diameter of 2.25 and 2.75 times the HPBW, respectively.We classify a circular region centered on pixel i a pointsource at high (∼ 4σ) confidence if where ⟨d⟩ (cir) i and ⟨d⟩ (ann) i are the flux per beam averaged over the circle and annulus centered at the pixel i, respectively (the noise σ rms is defined in Section II).For each pixel in the radio map that passes this criteria, we mask a circular region (of diameter 0.75 times the HPBW) centered on the pixel.
In addition to these conventional point sources, there is a feature (located near x = 4 kpc, y = 0 kpc in Figure 1) that is likely an artifact of the imaging process.As this feature does not have the intensity distribution of a point source, it was not identified by our point source algorithm, and we mask it by hand. 5

Center Mask
The center of M31 is the brightest source of radio emission in the galaxy.While dark matter-induced emission would also peak in this region, much of the observed emission is likely due to difficult-to-model baryonic processes.Limits on annihilation can be set by using only this central emission [36,38,39], but the intensity of the dark matter signal here is sensitive to the diffusion parameter for D 0 ≳ 1 × 10 28 cm 2 /s (as shown in Figure 16(b)).For these reasons, we also mask the center of M31 in our analysis and set limits on dark matter using the region outside the center.Here, the lower signal rate is off-set by the lower background, and the differing morphologies of the signal and background can be used to set limits less sensitive to uncertainties in the diffusion coefficient.Our point-source masking technique also identifies a source at the center of M31, but the default point-source mask is too small to cover the entire bright center region.To determine the size of the central circular mask, we plot the intensity, averaged over concentric circular annuli (excluding pixels in point-source masks), as a function of 2D radius ρ = (x 2 + y 2 ) 1/2 in Figure 17.We mask the central region out to the minimum of this averaged flux, at ρ = 0.93 kpc.The intensity map with point sources and the center masked is shown in Figure 18(a).

Ring and Outside Masks
Finally, we must construct a mask for the elliptical ring of bright emission in the star forming region of the M31 disk [89,90].Given the morphology of this feature, it cannot be due to dark matter annihilation.Masking it therefore does not risk removing a potential signal and setting overly-strong constraints.
To construct the mask, we first fit the data to a sum of a uniform template and a Gaussian elliptical ring template which have the forms The Global Fit has the parameter values fit to the data with the center and point sources masked, while the Signal-Region Masked fit is over data with the additional mask over the central signal-rich region applied.We separately show the parameters after fitting to the entire M31 data set (labeled "Full Map Analysis"), and the data in the x > 0 right-hand side of Figure 1 (labeled "Right-Only Analysis," see Section VIII). where is the elliptical radius, µ = (µ 1 , µ 2 , µ 3 ) are free parameters of the model that control the shape and size of the ring, and w = (w 1 , w 2 ) control the intensity of each component of the background.For the remainder of the paper, we will refer to µ as the background morphological parameters and w as the background coefficients.The total background model is As the dark matter-induced annihilation signal is expected to be small at the radius of the ring, we can fit our model to the ring independent of the signal model.We minimize the χ 2 statistic between the observed flux (d i in pixel i at location x i ) and the ring plus uniform background model with respect to all components of w and µ.The resulting best fit values for the morphological parameters are listed in Table V in the first section (labeled "Full Map Analysis") and second column (labeled "Global Fit").
Figure 19 shows the radio data (with point sources and galactic center masked) and the globally fit background model averaged over concentric elliptical annuli (with the same eccentricity as the globally fit ring model) as a function of R e .We show a heatmap of the globally-fit background model superimposed with simulated errors in Figure 18(b).Our method for simulating random errors correlated over the beam size (which is much larger than the pixel size) is explained in Appendix B.
It is clear from Figure 19 that the emission outside of the ring (R e ≳ 15 kpc) is significantly brighter than the emission inside (R e ≲ 5 kpc).As the dark matter signal is expected to drop with distance from the center, it cannot be responsible for this excess emission outside the ring.Therefore, we mask exterior to the ring.The width of the best-fit Gaussian is too broad to completely mask the emission out to the level of statistical noise in the region interior to the ring.We instead mask the ring inward to 1 standard deviation from the peak of the Gaussian ring model.That is, we mask all pixels satisfying using the globally-fit values for the morphological parameters µ.where µ1 is taken to be the globally fit value.
The inner boundary of the ring mask is shown in black contours in each panel of Figure 18.The interior of this contour (minus the center and pixels masked as part of point sources) is the search region that will be used to constrain dark matter annihilation.

B. Background Model of the Search Region
Having selected our search region, we must now define our background model that we will use to construct our background-only and signal plus background hypotheses.The model has the functional form of Eq. ( 53) (used to define the ring and outside-region mask in Section VII A 3), but as the intensity map has the potential to be signal-rich, we must fix the morphological parameters when calculating our limits more carefully than when we initially defined the search region.
If there was a (known) dark matter signal in the data, the parameters of the background model would be most accurately found by subtracting that signal from the data and then fitting our background model to the result.Alternatively, if there is no dark matter signal in the data, the parameters of the background model would be most accurately found by fitting the background model to the data itself.With only the point-source and circular center masks applied, the best fit parameters for the background model are sensitive to the (unknown) presence of signal in the data.We avoid the risk of the background model absorbing any signal present in the data by leveraging the morphology of the signal maps, which peak towards the center of M31.
Unlike the procedure for constructing the globally-fit background model, if we fit the background model only using data outside the center of M31 (where dark matter contributes less to the radio flux), then fits with and without signal subtracted will be more in agreement.The level of statistical agreement will increase as we mask more of an assumed signal.To maximize the amount of signal masked for a given area masked, the mask should have the shape of a region bounded by a contour of constant signal intensity.We will call the inner signal-rich region the "signal region," and the mask that covers it the "signal-region mask." Our strategy then is to mask the signal region, and fit the parameters of Φ b to the data exterior to the mask (including data outside the search region).This fit will allow us to define where μ are the morphological parameters, fit outside the signal region and fixed for the rest of the analysis, while w are free parameters which set the amplitude of the various components of the background.These free parameters will be fit to the data (or pseudo-data) in the search region when we set limits on the presence of a dark matter signal.The region exterior to the signal-region mask, used to find μ, must be sufficiently signal-poor so that statistical tests that distinguish between signal plus background and background-only hypotheses obtain the same results regardless of whether μ is determined by assuming the presence or absence of signal in the data.
To identify this signal-poor region, we use as a benchmark signal the flux from dark matter with m χ = 38.6GeV, ⟨σv⟩ = 2.2 × 10 −25 cm 3 /s and a diffusion normalization of D 0 = 1 × 10 28 cm 2 /s.As the leakage out of the masked signal region is minimal and the signal morphology depends only weakly on the choice of mass and diffusion parameters, the resulting fits can be applied to signals with other values of m χ and D 0 .The cross-section is chosen to be approximately an order of magnitude larger than the best fit value from the GCE [16][17][18][19] and existing limits from dwarf galaxies [6][7][8][9][10]12].Our fitting procedure will ensure that the background model is not significantly influenced by the presence of dark matter signals of this intensity and weaker in the data.
We make a series of candidate signal-region masks that intersect the semi-major axis at x values between [1.0 − 8.8] kpc.For each of these masks, we fit our background model (Eq.( 53)) to the remaining unmasked data with signal subtracted (defined as "Fit A") or without signal subtracted ("Fit B") by minimizing Eq. ( 54) with respect to all components of w and µ.We take the sum in Eq. ( 54) to be over pixels not covered by the candidate signal region mask, the center mask, or point source masks.
For each candidate signal-region mask, we determine if Fits A and B of the morphological parameters lead to statistically indistinguishable results when testing for the presence of signal.To compare the background-only and signal plus background hypotheses, we introduce a test statistic, defined as This statistic will also be used in our methodology for setting limits on dark matter (see Section VIII).In Eq. ( 57), {d i } are the differential flux values in each pixel of the intensity map, and the sum runs over pixels i that are in the search region (defined in Section VII A). w b (w s+b ) are the most-likely values of the background coefficients, w = (w 1 , w 2 ), under the background-only (signal plus background) hypothesis and are determined analytically for each intensity map.The test statistic is constructed such that higher test statistic values imply that the intensity map is more background-like.
The test statistic depends on the signal hypothesis being tested through the signal plus background model where the signal flux at the pixel centered at solid angle Ω i is given by using the differential flux calculated in Section VI.The signal hypothesis is parameterized by the cross section ⟨σv⟩ and a vector θ, containing m χ , D 0 and all other default astrophysical parameters, given in Section IV.For a given candidate signal-region mask, we calculate distributions of test statistics from ensembles of pseudodata using background models with morphological parameters fixed by Fits A and B. These ensembles are drawn using the methods described in Appendix B, using the Φb i appropriate for each fit, summing over pixels in the search region for each calculation of λ ⟨σv⟩,θ .As the test statistic requires a choice of signal parameters, we use as our reference signal model m χ = 38.6GeV, D 0 = 1 × 10 28 cm 2 /s, and a cross-section for which the signal plus background hypothesis is easily distinguished from the background hypothesis: ⟨σv⟩ = 1.1 × 10 −25 cm 3 /s.If -for a given signal mask -the distribution of teststatistics is indistinguishable between background Fits A and B, then the same will be true for the result of a statistical test for distinguishing signal plus background from background.This means the morphological parameters can be fit to the data using that signal region mask without absorbing potential signal from the data.Figure 20 shows the mean test statistic using candidate Fits A and B as a function of the distance from the origin that the signal region mask intersects the semi-major axis.For a mask which intersects the semi-major axis at x = 3.65 kpc, the means of the distributions of the test statistic from Fits A and B agree within statistical noise.Selecting this signal region mask (shown with a red contour in Figure 18), we set μ to the best-fit morphological parameters of Fit B. The values of the morphological parameters from this fit are shown in Table V (the "Signal-Region Masked" column of the "Full Map Analysis" section).To construct our signal plus background and background-only hypotheses, we will use the background model with the morphological parameters fixed to μ and the coefficients w = (w 1 , w 2 ) free to float to their most likely values for each hypothesis.

C. Limits on a Signal Model
Having fixed our background model morphology, we now describe our statistical approach to setting limits on dark matter annihilation in M31, using the data in the entire search region.In order to maximize the statistical power in the morphology of the signal when setting limits, we use the CL s method [91] with pixel-level radio data and templates.
As in Section VII B, we parameterize our signals using the cross section ⟨σv⟩, and the parameter vector θ which includes the dark matter mass m χ and diffusion coefficient D 0 .We will use the test statistic λ ⟨σv⟩,θ from Eq. (57) to distinguish between background-like and signal plus background-like intensity maps.
Statistical inference for a signal parameterized by ⟨σv⟩ and θ requires the probability distributions of λ ⟨σv⟩,θ under our background-only and signal plus background hypotheses.We construct these probability distributions by generating an ensemble of simulated observations of M31 under each hypothesis and calculating λ ⟨σv⟩,θ for each simulated observation.The simulated observations are generated under the background (signal plus background) hypothesis by superimposing Φb ( Φs+b ) with randomly drawn noise maps.The probability distribution from which the noise maps are drawn has correlations between nearby pixels, as expected due to the Gaussian beam of the observations (described in Section II).More details on our procedure for producing pseudo-data are described in Appendix B.
In Figure 21, we show sample distributions of λ ⟨σv⟩,θ for m χ = 38.6GeV, D 0 = 1×10 28 cm 2 /s and two choices of ⟨σv⟩ (1.1 × 10 −26 cm 3 /s and 4.6 × 10 −26 cm 3 /s).In these examples, the blue and red histograms are the distributions of the test statistic assuming background and signal plus background (respectively) in arbitrary units.The solid curves are the Gaussian approximations of each distribution.The vertical green line in each plot is the value of λ For an arbitrary signal hypothesis parameterized by ⟨σv⟩ and θ, our distributions of λ ⟨σv⟩,θ can be used to approximate the probability distribution of λ ⟨σv⟩,θ assuming background: and signal plus background: A given observation with test-statistic λ (obs) ⟨σv⟩,θ then has a CL b value given by the probability of seeing data more background-like than observed, assuming the background-only hypothesis is correct: dλ p ⟨σv⟩,θ (λ|b).(62) Similarly, the CL s+b value for the observation is the probability of seeing a more background-like intensity map than that observed, assuming that the signal plus back-ground hypothesis is correct: dλ p ⟨σv⟩,θ (λ|s+b, ⟨σv⟩, θ).

(63)
The ratio CL s λ (obs) ⟨σv⟩,θ , ⟨σv⟩, θ ≡ CL s+b /CL b can then be interpreted as the probability of signal parameters greater than ⟨σv⟩, given data.A 95% confidence level exclusion therefore corresponds to a signal for which The expected 95% confidence limits correspond to signal parameters for which the median test statistic under the background hypothesis (CL b = 0.5) leads to CL s = 0.05.The 1 and 2σ errors of this expected limit are calculated using the corresponding percentiles of the background distribution.
Example CL curves are shown in Figures 21(c) and 21(d) (corresponding to the distributions of the test statistics in Figures 21(a) and 21(b), respectively).The CL s curve of each plot (shown in grey) is dominated by statistical noise in the simulated intensity maps when CL b and CL s are small.As seen in this example, we generically find that λ (obs) ⟨σv⟩,θ is 5 − 6σ larger than the mean of the background-only distribution for our search region, 6 implying that CL b of the observed test statistic is ∼ 2 × 10 −8 .
The fact that the observed test statistics are located on the tail of the CL b distributions is a consequence of observed radio intensities that are much less signal-like than the signal-free background model.We confirmed that this issue exists for all values of m χ and D 0 studied in this work.This means that while the background model does not describe the observations well, the observed deviations away from the background-only model are not compatible with the morphology of any signal.This will lead to limits that are much stronger than expected.
This issue requires further investigation.In The value of λ ⟨σv⟩,θ for the data is shown with the vertical green lines and the median expected test statistics from the background only hypothesis are shown with the dashed blue vertical lines.Each distribution is constructed from N = 20000 independent simulated maps.The smooth curves are Gaussian approximations of the distributions.Bottom row: example CL curves for the signal models shown above.The black curve is an approximation of CLs derived from the Gaussian approximations of the distributions of the test statistic.The test statistic that is excluded at 95% confidence according to the Gaussian approximation of CLs is shown with a vertical cyan line.
of signal parameters is excluded at 95% confidence). 7As can be seen, the residuals of the data with respect to the the signal plus background model predicts a larger flux than the background-only model in this region.
In Figures 22(b) and 22(c) we show elliptically averaged radio emission as well as best fit models in the left search region (masking the x > 0 data) and right search region (masking x < 0), respectively.As can be seen, the large negative excursion in the left search region is the source of the negative residuals we identified in Fig. 22(a) and can be traced specifically to the large negative observed flux located around (x, y) ∼ (−2, 1) kpc.This negative flux is clearly visible in Figures 1 and 18(a).The low flux measurement (well in excess of a 2σ deviation given the expected measurement errors) may be due to the over-subtraction of a point source external to M31.
Critically, given our understanding of the dark matter distribution within M31, dark matter emission cannot create such a region of low emission close to the center of M31, even if the overall baseline of zero radio flux was mismeasured.Thus, considering only the part of the data without this region of anomalously low emission will not set overly-optimistic limits on dark matter annihilation.Indeed, we will see in the next section that it sets more conservative and weaker bounds.

VIII. CONSTRAINTS ON ANNIHILATING DARK MATTER IN M31
In light of the findings from last section, we first calculate limits for the original search region (within the black contour in Fig. 18(a)) and compare them to the limits from the right and left search regions.We set 95% CL s exclusion limits for dark matter with mass in the range [ assuming the default diffusion parameter normalization D 0 = 1 × 10 28 cm 2 /s.This includes the observed limits for the original, right and left search regions with black curves of various styles, alongside the 1 and 2σ variation around the expected limits for the original search region.As expected, the observed limit for the original search region is far stronger than the 2σ variation assuming the background hypothesis.The limits from the left search region are nearly as strong as those from the original search region while the limits from the right search region are approximately as strong as expected, confirming that the stronger-than-expected limits come entirely from the left-hand side of the M31 emission, where the region of negative flux is located.
To get the most accurate limits using only the right side of the map, we recalculate the search region and refix the background morphological parameters with the left side of the data masked, using the steps described in Sections VII A and VII B. The resulting morphological parameter values used to recalculate the search region are shown in the "Right-Only Analysis" section of Table V in the second column (labeled "Global Fit").The new search region is bounded by the black contour in Figure 24.We show the resulting test statistics as a function of candidate signal region mask size for background Fits A and B in Figure 25.Based on this, we select the signal region mask for the right-only analysis which intersects the semi-major axis at x = 3.65 kpc, where the mean test statistics from Fits A and B agree within statistical error.This turns out to be the same signal region mask that we found for our full map analysis.This signal region mask is bounded by the red contour in Figure 24.We fix the morphological parameters of the background model µ to the best fit values from background model B with this signal region mask.The resulting values for the morphological parameters are shown in the third column (labeled "Signal-Region Masked") of the second half (labeled "Right-Only Analysis") of Table V.
Using this background model for the right-only analysis, we show in Figures 26(a In Figure 27, we show our limits on ⟨σv⟩ as a function of m χ , using the right-only analysis.These results our most conservative and robust limits on dark matter annihilation using the radio observations of M31.In both panels, the green and yellow bands quantify the 1 and 2σ statistical error of our expected limits for our default value of D 0 and for the weighted averaging scheme, introduced in Section V C (and otherwise default parameters).Each panel quantifies the effects of different systematics on our results.In Figure 27(a) we show the observed 95% confidence exclusion limit from each spherically averaging procedure for our default value of D 0 .As can be seen, the limits do not depend on the averaging scheme used.In Figure 27(b), we show the observed limits as D 0 is varied.The limits are relatively insensitive to variations of the diffusion coefficient in the range 3 × 10 27 cm 2 /s ≤ D 0 ≤ 3 × 10 28 cm 2 /s.The limits become about a factor of three weaker when D 0 changes from 3 × 10 28 cm 2 /s to 8 × 10 28 cm 2 /s for m χ ≳ 10 GeV.We also show best fit contours to the GCE from previous analyses [16,17,19], as well as limits from Milky Way dwarfs using Fermi Pass 8 data [12].Our limits do not exclude parameters that fit the GCE.Although our limits are weaker than the dwarf limits, they are robust to astrophysical uncertainties.FIG. 26.Same as Figure 21 but for the right-only analysis.The cross-section shown here is close to the expected and actual 95% confidence limit.The expected limit is almost the same as the actual limit since the test statistic from the data is very close to the 50 th percentile test statistic from background pseudo-data.FIG.27.Expected and actual 95% confidence limits from the right-only analysis, using the data from the search region shown in Figure 24.The two panels show the variation of the observed limits due to (a) changes in the averaging procedure (introduced in Section V C) for our default diffusion coefficient normalization and (b) changes in the diffusion coefficient normalization, D0 for our weighed averaging scheme.Both panels have the expected limits obtained using the default value of D0 (1 × 10 28 cm 2 /s) and the weighted averaging scheme.The contours show best fits to the GCE [16,17,19] and the solid blue lines show limits from dwarfs using Fermi Pass 8 data [12].

IX. CONCLUSION
In this work, we have set robust and conservative limits on dark matter annihilation to b b using the 8.35 GHz Effelsberg radio map of M31, which is sensitive to the predicted synchrotron emission of the e ± produced in the cascade b decays.These limits are based on a numeric solution to the diffusion-loss equation that accommodates non-uniform parameters, and an astrophysical model that uses observations of the gas, dust, starlight, and magnetic fields of M31.Our ISRF model for the starlight is derived directly from ugriz luminosity data, which led to notably larger values for ρ * in the center of M31 compared to previous works.
Unlike previous studies, our numerical solution to the diffusion-loss equation allows for position dependent diffusion and loss coefficients.Our method still requires spherical averaging of the background model.Though we have shown that our final limits are insensitive to the averaging procedure, additional work is needed to develop a numeric solution which is adapted to the axisymmetry of M31.
Our limits are based on the morphology of the observed flux, and our results are independent of the true zero-flux level of the intensity map.The limits are based on the radio flux only interior to the bright ring of radio emission in M31, allowing us to use data-driven models of the background based on signal-poor regions of the observed intensity map.Due to a localized anomaly of low radio flux in the search region, we choose conservatively to select only the half of the dataset without this anomaly with which to set our limits.
Our limits on dark matter annihilation in M31 do not exclude best-fit models to the GCE (shown as contours in Figures 27(a)& 27(b)) and are weaker than those found in previous radio studies [36][37][38][39].These weaker limits are due in part to the fact that in our analysis we mask the center of the galaxy, where the signal intensity is maximum -however, this choice minimizes the sensitivity to unknown astrophysical parameters at the galactic center.The weaker limits are also likely due to the differences in our astrophysical model of M31 compared to previous work.In particular, the core of M31 is much more luminous in starlight than a simple scaling of the comparable region of the Milky Way would suggest.This increased starlight flux results in increased energy losses of e ± into X-rays through inverse Compton scattering, reducing the flux of dark matter-induced radio waves.Though beyond the scope of this work, this suggests that an analysis of constraints from X-ray emission in the center of M31 from dark matter annihilation may set interesting limits.
The sensitivity of the limits to the astrophysical conditions within M31 are notable; though in this work we have taken care to construct an accurate model of M31 based on observations, future measurements and astronomical input would likely improve the model and the resulting limits.Similar analysis is likely necessary for constraints on dark matter annihilation via radio waves in other systems beyond M31.Starting from the first row, we reduce M to upperechelon form.The components of U n+1 can then be solved for, starting from the final component and recursively solving for the other components in reverse order.For the most general matrix, this procedure would require O(n E × n r ) 3 operations, which would be computationally prohibitive.In our case, the matrix M is tridiagonal with a fringe, which requires only O(n 2 r × n E ) operations.
Initially, we set the time-step to an approximation of the maximum timescale of the problem: where r s is the scale radius of the dark matter distribution, given in Section III.Using ∆t = ∆t 0 , we iteratively solve for U n+1 from U n until each component of the two vectors is different by less than 1 part in 10 3 .We then reduce the timestep by a factor of 2 and repeat, starting with the final result from the last time-step and iterating until the same convergence criteria is met.We repeat this procedurereducing the time-step by a factor of 2, and achieving convergence of the solution -until there have been at least 5 different values of ∆t and U n converges in one step for 3 values of ∆t in a row.We find that these convergence criteria are conservative as convergence is achieved after 5 values of ∆t for all solutions that we examined.

FIG. 1 .
FIG.1.Smoothed non-thermal radio intensity map of M31 from Ref.[35], showing the flux per unit frequency per beam averaged over a frequency bandwidth of 1.1 GHz.The HPBW projected into the plane of M31 is 0.340 kpc and the rms noise is given by σrms = 0.25 mJy/beam in the inner 9.13 kpc × 9.13 kpc region and σrms = 0.3 mJy/beam in the rest of the map.Digitized data for this figure was provided by the authors of Ref.[35].

1 )
FIG. 2. The number of e ± in final states per unit energy per annihilation of dark matter into b b for a representative sample of dark matter masses mχ.

FIG. 6 .
FIG.6.Number density (left axis) and surface density (right axis) of (a) HI and (b) H2 gas in the plane of the disk.The digitized and interpolated distributions from Ref.[73] are within the two vertical red lines.Outside these regions, we fit exponential extrapolations, matching the function values and first derivatives at the boundaries.

FIG. 8 .
FIG. 8. Diffusion coefficient as a function of R for E = 1 GeV and various values of z, with D0 = 10 28 cm 2 /s.

He = 3
.61 × 10 −16 GeV/s and b (0) ion = 1.74 × 10 −16 GeV/s.Our models for the density of each gas component were presented in Section IV C.

10 − 3 2 EFIG. 9 .
FIG. 9. (a) The energy dependence of the loss coefficient for various values of R at z = 0. (b) The energy dependence of the total loss coefficient (solid line) and its subcomponents at R = z = 0. Inverse Compton and synchrotron losses, which have the same energy dependence, are shown as the dashed line, bremsstrahlung as dot-dashed, and Coulomb losses as the dotted line.

Figure 9 (
Figure 9(b) shows the total loss coefficient at the origin and the contributions to it from the individual processes discussed in this subsection.Coulomb losses dominate at low energy, inverse Compton and synchrotron losses dominate at high energy, and bremsstrahlung only becomes marginally important at intermediate energies for R ≃ 10 kpc due to the large concentration of interstellar gas in the ring-like structure (discussed in Section VII B).

FIG. 11 .
FIG. 11.Dynamic timescale τ of M31 as a function of R and z for E = 0.5 GeV (the minimum energy contributing significantly to the 8.35 GHz synchrotron signal) and D0 = 3 × 10 27 cm 2 /s (the lower bound on the diffusion coefficient).

FIG. 14 .
FIG. 14.Energy response of e ± producing synchrotron emission of frequency ν = 8.35 GHz for a variety of magnetic field strength values.

FIG. 17 .
FIG. 17. Observed flux averaged over concentric circular annuli of radius ρ in the plane of the sky, not including pixels that are in the point source mask.The errors in the annulus averaged flux in a particular bin are found by averaging the rms noise over the bin and dividing by the square root of the number of beams in the bin.The radius of the center circular mask is shown with a red vertical line.

FIG. 18 .
FIG.18.Intensity maps of the (a) radio data and (b) simulated pseudo-data using the globally-fit background model, with point source and center masks (described in Sections VII A 1 and VII A 2, respectively) applied.The method of simulating the pseudodata is described in Appendix B. The search region (used to set limits on dark matter annihilation, see Section VII A 3) consists of the unmasked pixels within the black contour.The signal region, masked when defining signal-independent background templates (see Section VII B), is interior to the red contour.

FIG. 19 .
FIG.19.Synchrotron data and globally fit background model (parameters given in the "Full Map Analysis" section of Table V) averaged over elliptical annuli as a function of Re(x, µ1) where µ1 is taken to be the globally fit value.

FIG. 20 .
FIG.20.Mean test statistics from background-only pseudodata for a series of candidate fits of the morphological parameters of the background model.Each fit comes from minimizing Eq. (54) outside of the candidate signal region mask that intersects the semi-major axis of M31 at x mask , assuming the presence (green) or absence (black) of dark matter signal in the data.The distributions of test statistics are constructed for a signal from dark matter with mχ = 38.6GeV and ⟨σv⟩ = 1.1 × 10 −25 cm 3 /s and default diffusion normalization of D0 = 1 × 10 28 cm 2 /s.The size of the signal-region mask that we select is shown with the red vertical line.
FIG. 21.Top row: example histograms of λ ⟨σv⟩,θ for background and signal plus background hypotheses for a signal parameterized by mχ = 38.6GeV and D0 = 1 × 10 28 cm 2 /s.Plot (a) has ⟨σv⟩ = 1.1 × 10 −26 cm 3 /s and (b) has ⟨σv⟩ = 4.6 × 10 −26 cm 3 /s.The value of λ ⟨σv⟩,θ for the data is shown with the vertical green lines and the median expected test statistics from the background only hypothesis are shown with the dashed blue vertical lines.Each distribution is constructed from N = 20000 independent simulated maps.The smooth curves are Gaussian approximations of the distributions.Bottom row: example CL curves for the signal models shown above.The black curve is an approximation of CLs derived from the Gaussian approximations of the distributions of the test statistic.The test statistic that is excluded at 95% confidence according to the Gaussian approximation of CLs is shown with a vertical cyan line.

FIG. 22 .
FIG. 22. Flux averaged over concentric elliptical annuli as function of Re(x, µ1) (with µ1 given by the global fit of the Full Map analysis) along with the best fit background model and the excluded signal plus background model for (a) the original search region, (b) the left half search region, and (c) the right half search region.All signal plus background models shown have mχ = 38.6GeV and D0 = 1 × 10 28 cm 2 /s.For the excluded signal plus background model, we take the lowest value of ⟨σv⟩ that leads to 95% exclusion for the values of mass and diffusion normalization plotted.

1 ×FIG. 23 .
FIG.23.95% confidence limits on dark matter annihilation assuming D0 = 1 × 10 28 cm 2 /s from our Full Map analysis.The 1σ and 2σ expected limits from the original search region are shown in green and yellow, with the observed limits for this search region shown with a solid line (labeled "original").The dotted and dashed lines are the actual limits from the data in the right (x > 0) and left (x < 0) half of the search region, respectively.

FIG. 24 .FIG. 25 .
FIG. 24.As Figure 18(a) but for the right-only analysis.The search region contour is recalculated with the left side of the image masked, and thus differs slightly from the search region of the full map analysis.
) and26(b)  examples of the test statistic distributions and CL curves for pseudodata in the right search region and show sample distributions, for m χ = 38.6GeV, D 0 = 1 × 10 28 cm 2 /s, and ⟨σv⟩ = 7.4 × 10 −26 cm 3 /s.This is the cross section that is excluded at approximately 95% confidence for these values of m χ and D 0 .

TABLE V .
Best-fit morphological parameters for the ring.