Search for $\gamma$-ray emission from dark matter particle interactions from Andromeda and Triangulum Galaxies with the Fermi Large Area Telescope

TL;DR

This study analyzes eight years of Fermi-LAT gamma-ray data from the Andromeda and Triangulum galaxies to search for signals of dark matter interactions, setting limits on dark matter properties and exploring astrophysical contributions.

Contribution

It provides new constraints on dark matter annihilation cross sections and decay times using gamma-ray observations of M31 and M33, considering astrophysical backgrounds.

Findings

No significant dark matter gamma-ray signal detected.

Limits probe the thermal relic cross section for dark matter masses up to tens of GeV.

Constraints are relevant for dark matter interpretations of the Galactic center excess.

Abstract

The Andromeda (M31) and Triangulum (M33) galaxies are the closest Local Group galaxies to the Milky Way, being only 785 and 870 kpc away. These two galaxies provide an independent view of high-energy processes that are often obscured in our own Galaxy, including possible signals of dark matter (DM) particle interactions. The Fermi Large Area Telescope (Fermi-LAT) preliminary eight year list of sources includes both M31, which is detected as extended with a size of about 0.4$^\circ$, and M33, which is detected as a point-like source. The spatial morphology of M31 $\gamma$-ray emission could trace a population of unresolved sources and energetic particles originating in sources not related to massive star formation. Alternatively, the $\gamma$-ray emission could also be an indication of annihilation or decay of DM particles. We investigate these two possibilities using almost 10 years of data from the Fermi LAT. An interpretation that involves only a DM $\gamma$-ray emission is in tension with the current limits from other searches, such as those targeting Milky Way dwarf spheroidal galaxies. When we include a template of astrophysical emission, tuned on $\gamma$-ray data or from observations of these galaxies in other wavelengths, we do not find any significant evidence for a DM contribution and we set limits for the annihilation cross section that probe the thermal cross section for DM masses up to a few tens of GeV in the $b\bar{b}$ and $\tau^+\tau^-$ channels. For models where the DM substructures have masses above $10^{-6}$ solar masses our limits probe the DM interpretation of the Fermi LAT Galactic center excess. We provide also the lower limit for the DM decay time assuming the same spatial models of the DM distribution in M31 and M33.

Figures

Tables

Table 1. Table 1: Summary of the most important parameters of CLUMPY to model the substructure contribution to the total J − limit-from 𝐽 J- / D − limit-from 𝐷 D- factor in M31 and M33. Note that the virial mass of M33 assuming either an NFW or Burkert profile has been derived with respect to Δ = 97.2 Δ 97.2 \Delta=97.2 times the critical density of the universe instead of Δ = 200 Δ 200 \Delta=200 as for all other profiles. Moreover, r s subscript 𝑟 𝑠 r_{s} and ρ s subscript 𝜌 𝑠 \rho_{s} of M33’s NFW profile have been derived with respect to Δ = 97.2 Δ 97.2 \Delta=97.2 , too. The choice of Δ = 97.2 Δ 97.2 \Delta=97.2 has been made by the authors of Corbelli  et al. ( 2014 ); Lopez Fune  et al. ( 2017 ) .

	smooth profile M31	smooth profile M33	${\alpha }_{\text{sub}}$	$f_{\text{sub}}$	$M_{\text{min}}\left[M_{\odot }\right]$	$c_{200,\text{sub}}$
min	Burkert Tamm et al. (2012)	Burkert Lopez Fune et al. (2017)	$1.9$	$0.12$	${10}^6$	Moliné et al. (2017)
	$r_s=9.06\text{kpc,}$	$r_s=9.6\text{kpc,}$
	${\rho }_s=3.68\times {10}^7{M}_{\odot }{\text{/kpc}}^3$	${\rho }_s=1.23\times {10}^7{M}_{\odot }{\text{/kpc}}^3$
	$M_{200}=7.9\times {10}^{11}{M}_{\odot }$	$M_{97.2}=3.0\times {10}^{11}{M}_{\odot }$
med	Einasto Tamm et al. (2012)		$1.9$	$0.19$	${10}^{-6}$
	$r_s=178\text{kpc,}$
	${\rho }_s=8.12\times {10}^3{M}_{\odot }{\text{/kpc}}^3$	NFW Lopez Fune et al. (2017)
	$M_{200}=1.13\times {10}^{12}{M}_{\odot }$	$r_s=22.41\text{kpc,}$
max	adiabatically contracted NFW Seigar et al. (2006)	${\rho }_s=2.64\times {10}^6{M}_{\odot }{\text{/kpc}}^3$	$2.0$	$0.45$	${10}^{-12}$
	$r_s=25.65\text{kpc,}$	$M_{97.2}=5.4\times {10}^{11}{M}_{\odot }$
	${\rho }_s=4.44\times {10}^7{M}_{\odot }{\text{/kpc}}^3$
	$M_{200}=5.7\times {10}^{11}{M}_{\odot }$

Table 2. Table 2: Summary table of Fermi -LAT data selection criteria used for this paper’s DM analysis.

Selection	Criteria
Observation Period	2008 August 4 to 2018 March 4
Mission Elapsed Time (s)121212$Fermi$ Mission Elapsed Time is defined as seconds since 2001 January 1, 00:00:00 UTC	239557417 to 541779795
Energy Range (GeV)	0.3–1000
Fit Region (M31)	16${}_{}{}^{\circ }\times$16${}^{\circ }$ centered on $(\alpha ,\delta )=(10{.}^{\circ }685,41{.}^{\circ }269)$
Fit Region (M33)	16${}_{}{}^{\circ }\times$16${}^{\circ }$ centered on $(\alpha ,\delta )=(23{.}^{\circ }462,30{.}^{\circ }660)$
Zenith Range	90${}^{\circ }$ and PSF2 and PSF3 for $E\in [0.1,0.3]$ GeV
	100${}^{\circ }$ and PSF1, PSF2 and PSF3 for $E\in [0.3,1.0]$ GeV
	105${}^{\circ }$ all PSF types for $E>1$ GeV
Data Quality Cut131313Standard data quality selection: DATA_QUAL$>$1 LAT_CONFIG==1 roicut==yes with the $gtmktime$ Science Tool	yes

Table 3. Table 3: Summary table for the T ​ S 𝑇 𝑆 TS of detection ( T ​ S 𝑇 𝑆 TS ) and extension ( T ​ S E ​ X ​ T 𝑇 subscript 𝑆 𝐸 𝑋 𝑇 TS_{EXT} ) in our analysis of Fermi -LAT data in the M31 ROI for the disk (left side) and Gaussian templates (right side) and using in the fit the Off IEM.

$E$	$TS$	$T{S}_{EXT}$	${\theta }_{\mathrm{EXT}}$ [deg]	$TS$	$T{S}_{EXT}$	${\theta }_{\mathrm{EXT}}$ [deg]
	Disk template			Gaussian template
$>0.1$ GeV	110	15.3	$0{.}^{\circ }33\pm 0{.}^{\circ }03$	109	13.9	$0{.}^{\circ }41\pm 0{.}^{\circ }09$
$>0.3$ GeV	98	13.6	$0{.}^{\circ }33\pm 0{.}^{\circ }04$	97	12.8	$0{.}^{\circ }42\pm 0{.}^{\circ }10$
$>0.5$ GeV	82	9.6	$0{.}^{\circ }32\pm 0{.}^{\circ }04$	81	8.6	$0{.}^{\circ }37\pm 0{.}^{\circ }09$
$>1.0$ GeV	58	9.3	$0{.}^{\circ }31\pm 0{.}^{\circ }05$	58	8.2	$0{.}^{\circ }31\pm 0{.}^{\circ }09$

Table 4. Table 4: Summary table for the significance of detection and extension in our analysis of Fermi -LAT data in the M31 ROI and M33 ROIs with the Off and Alt IEMs. We assume here a uniform disk spatial template.

IEM	$TS$	$T{S}_{EXT}$	${\theta }_{\mathrm{EXT}}\left[deg\right]$	$TS$
	M31			M33
Off	110	15.3	$0{.}^{\circ }33\pm 0{.}^{\circ }03$	$39$
Alt 1	$90.1$	$12.4$	$0{.}^{\circ }32\pm 0{.}^{\circ }04$	$42$
Alt 2	$100$	$10.2$	$0{.}^{\circ }37\pm 0{.}^{\circ }06$	$36$
Alt 3	$89$	$12.3$	$0{.}^{\circ }32\pm 0{.}^{\circ }04$	$42$
Alt 4	$85$	$9.8$	$0{.}^{\circ }42\pm 0{.}^{\circ }10$	$37$
Alt 5	$86$	$12.1$	$0{.}^{\circ }32\pm 0{.}^{\circ }04$	$39$
Alt 6	$94$	$9.6$	$0{.}^{\circ }36\pm 0{.}^{\circ }07$	$34$
Alt 7	$106$	$11.1$	$0{.}^{\circ }43\pm 0{.}^{\circ }07$	$40$
Alt 8	$84$	$10.0$	$0{.}^{\circ }32\pm 0{.}^{\circ }04$	$35$

Table 5. Table 5: Value of the T ​ S 𝑇 𝑆 TS and difference of likelihood with respect to the disk model ( Δ ​ log ⁡ L Δ 𝐿 \Delta\log{L} ) for M31 and M33 spatial models considered in our analysis. We show the results for the point source, disk and Gaussian geometrical models and for the templates taken from other wavelengths.

M31			M33
Model	$TS$	$\Delta \log L$	Model	$TS$	$\Delta \log L$
Disk	98	0	Point	41	0
Gauss	97	1	Gauss	41	0
Point	80	5	Optical ($3552$ A)	41	0
Spitzer (3.6$\mu$m)	94	3	Spitzer (24$\mu$m)	45	-2
Herschel (160$\mu$m)	75	9	Spitzer (160$\mu$m)	42	0
$N_H$	65	18	2MASS (2.2$\mu$m)	41	0

Table 6. Table 6: Summary table for the T ​ S 𝑇 𝑆 TS of detection ( T ​ S 𝑇 𝑆 TS ) and flux in our analysis of Fermi -LAT data in the M31 and M33 ROI (in the left and right side) with 3 uniform annuli as explained in the text. The top (bottom) part of the table is with the isotropic template free to vary (fixed to the control region).

ISO free	annulus 1	annulus 1	annulus 3	annulus 1	annulus 1	annulus 3
TS	0	0	8	0	0	0
Flux [ph/cm2/s]	$9.45\cdot {10}^{-10}$	$4.53\cdot {10}^{-10}$	$4.62\cdot {10}^{-9}$	$1.09\cdot {10}^{-9}$	$9.12\cdot {10}^{-10}$	$1.91\cdot {10}^{-9}$
ISO fixed	annulus 1	annulus 1	annulus 3	annulus 1	annulus 1	annulus 3
TS	5	0	20	0	0	0
Flux [ph/cm2/s]	$1.82\cdot {10}^{-9}$	$8.05\cdot {10}^{-10}$	$6.56\cdot {10}^{-9}$	$1.09\cdot {10}^{-9}$	$1.05\cdot {10}^{-9}$	$2.41\cdot {10}^{-9}$

Table 7. Table 7: Summary table for the T ​ S 𝑇 𝑆 TS , m χ subscript 𝑚 𝜒 m_{\chi} and ⟨ σ ​ v ⟩ delimited-⟨⟩ 𝜎 𝑣 \langle\sigma v\rangle of DM in M31 ROI considering SHS (left side) or SH case (right side). We assume here a b ​ b ¯ 𝑏 ¯ 𝑏 b\bar{b} annihilation channel.

SHS				SH
DM model	$m_{\chi }$	$⟨\sigma v⟩$	$TS$	DM model	$m_{\chi }$	$⟨\sigma v⟩$	$TS$
MIN	${25}_{-15}^{+40}$	${1.3}_{-0.6}^{+2.2}\cdot {10}^{-25}$	27	Burkert	${20}_{-5}^{+10}$	${1.4}_{-0.3}^{+1.0}\cdot {10}^{-25}$	43
MED	${30}_{-15}^{+35}$	${3.2}_{-1.7}^{+4.1}\cdot {10}^{-26}$	27	Einasto	${20}_{-5}^{+10}$	${5.2}_{-1.1}^{+2.9}\cdot {10}^{-26}$	55
MAX	${40}_{-20}^{+60}$	${1.1}_{-0.6}^{+2.4}\cdot {10}^{-26}$	20	Adiabatic	${20}_{-5}^{+10}$	${2.0}_{-0.5}^{+1.2}\cdot {10}^{-26}$	56

Table 8. Table 8: Same as Table 7 but for M33.

SHS				SH
DM model	$m_{\chi }$	$⟨\sigma v⟩$	$TS$	DM model	$m_{\chi }$	$⟨\sigma v⟩$	$TS$
MIN	${20}_{-10}^{+50}$	$5_{-3}^{+8}\cdot {10}^{-24}$	23	Burkert	${50}_{-30}^{+60}$	${2.6}_{-1.8}^{+3.3}\cdot {10}^{-24}$	25
MED	${30}_{-20}^{+60}$	$9_{-7}^{+16}\cdot {10}^{-25}$	18	Einasto	${90}_{-60}^{+110}$	${2.3}_{-1.6}^{+2.8}\cdot {10}^{-24}$	31
MAX	${25}_{-20}^{+45}$	$7_{-5}^{+20}\cdot {10}^{-27}$	13

Equations

$$\text{decay}:$$

$$\text{annihilation}:$$

$${\rho }_{\mathrm{ZHAO}}\left(r\right)=\frac{{\rho }_s}{{\left(\frac{r}{r_s}\right)}^{\gamma }{\left[1+{\left(\frac{r}{r_s}\right)}^{\alpha }\right]}^{\frac{\beta -\gamma }{\alpha }}}\text{,}$$

$${\rho }_{\mathrm{Ein}}\left(r\right)={\rho }_s\exp \left(-\frac{2}{\alpha }\left[{\left(\frac{r}{r_s}\right)}^{\alpha }-1\right]\right),$$

$${\rho }_{\mathrm{Burkert}}\left(r\right)=\frac{{\rho }_s}{\left(1+\frac{r}{r_s}\right)\left[1+{\left(\frac{r}{r_s}\right)}^2\right]}.$$

$$\ln \mathcal{L}(\mu ,\theta |\mathcal{D})=\sum _j\ln {\mathcal{L}}_j(\mu ,{\theta }_j|{\mathcal{D}}_j),$$

$$TS=2\Delta \log \mathcal{L}=2\ln \frac{\mathcal{L}(\mu ,\theta |\mathcal{D})}{{\mathcal{L}}_{\mathrm{null}}(\theta |\mathcal{D})},$$

Full text

I Introduction

It has been clear for many decades that the Universe includes a significant component of matter which is not made of any known particles of the Standard Model of particle physics. This form of matter is called dark matter (DM) because, though solid observational evidence exists for its gravitational influence from the earliest moments of the Universe’s history Zwicky (1933); Rubin et al. (1980); Olive (2003); Ade et al. (2014), no direct measurements have been made so far of its particle nature.

Though by no means the only possibility, a theoretically well motivated class of DM models include interactions between DM and Standard Model particles that are approximately as strong as the weak nuclear force, and a mass of similar scale ($\sim$10–1000 GeV). Such weakly interacting massive particles (WIMPs) would generically attain the observed DM density after thermal freeze-out in the early Universe. The canonical “WIMP” is a $\sim$100 GeV particle interacting through the $SU(2)_{L}$ weak force, although many other candidates have been proposed with a wide range of masses and interaction strengths Feng and Kumar (2008); Steigman et al. (2012).

The WIMP models provide a useful benchmark for DM searches designed to look for the present-day pair annihilation (or decay) of DM particles in regions of high density of DM. A thermally-averaged annihilation cross section of $\langle\sigma v\rangle\sim 3\times 10^{-26}$ cm3/s would provide approximately the correct WIMP relic density at present, and so experiments capable of seeing the present-day annihilation of DM with cross sections near this value have the sensitivity to either confirm or exclude a large number of theoretically interesting models.

In terms of couplings to Standard Model particles, there are many possibilities for dominant annihilation (or decaying) channels. Of particular interest is annihilation or decay resulting in monochromatic $\gamma$ rays, as this signature is more easily distinguished from other astrophysical sources. However, as this is a loop interaction (DM does not couple to photons directly), this channel is suppressed; thus, searches for this signature are challenging. In addition to the direct annihilation to pairs of $\gamma$-ray photons, if DM annihilates into pairs of other Standard Model particles, the resulting hadronization and/or decay will result in a continuum of $\gamma$ rays observable from Earth with an energy distribution that extends up to the rest mass of the DM particle. $\gamma$ rays produced in the local Universe are relatively unaffected by their propagation in the interstellar medium; thus they arrive at the Earth unscattered and unattenuated. This allows us to trace the $\gamma$-ray emission spatial distribution and spectral information back to its original source. Thus, $\gamma$-ray observations together with separate information or assumptions about the distribution of DM in the regions under study and models for the hadronization allow for measurement of, or determination of upper limits for, the annihilation cross section.

With this motivation in mind, $\gamma$-ray data measured by the Large Area Telescope (LAT) carried by the Fermi Gamma-ray Space Telescope (Fermi-LAT) are of great interest. The Fermi-LAT is a pair-conversion telescope. Incoming $\gamma$ rays pass through the anti-coincidence detector and convert in a Silicon strip tracker to $e^{+}/e^{-}$ pairs. The charged particle direction is reconstructed using the information in the tracker, and the energy is estimated from depositions in a CsI calorimeter. Detailed descriptions of the LAT and its performance can be found in dedicated papers Atwood et al. (2009, 2012). At the present time, the Fermi-LAT is one of the most sensitive instruments to DM particles with weak-scale mass and producing $\gamma$ rays. Analysis of the LAT $\gamma$-ray data can place strong limits on—or discover—DM annihilation with cross sections near the canonical thermal value into a wide variety of Standard Model particles Charles et al. (2016).

A large number of DM searches have been performed using Fermi-LAT data. Since annihilation rates are proportional to the square of the DM density, lower annihilation cross sections can be probed by targeting regions in the local Universe with the greatest densities of DM, such as the center of the Milky Way (MW) Goodenough and Hooper (2009); Hooper and Goodenough (2011); Boyarsky et al. (2011); Abazajian and Kaplinghat (2012); Huang et al. (2013); Gordon and Macias (2013); Daylan et al. (2016), satellite dwarf spheroidal galaxies of the MW Ackermann et al. (2011, 2014); Geringer-Sameth and Koushiappas (2011); Geringer-Sameth et al. (2015a); Ackermann et al. (2015a); Geringer-Sameth et al. (2015b); Drlica-Wagner et al. (2015), unresolved halo substructures Ackermann et al. (2012); Belikov et al. (2012); Berlin and Hooper (2014); Buckley and Hooper (2010), galaxy clusters Ackermann et al. (2010); Dugger et al. (2010), and the Large and Small Magellanic Clouds (LMC and SMC) Buckley et al. (2015); Caputo et al. (2016).

The direction with the predicted brightest $\gamma$-ray emission from DM is towards the Galactic center. The LAT observations of the Galactic center indicate that the region is brighter than expected from standard models for Galactic diffuse emission at GeV energies, and the spatial distribution is broadly consistent with our expectations for a DM signal (see, e.g., Ackermann et al. (2017a)). However, previously unconsidered astrophysical backgrounds could match the observed morphology and spectrum and the true source of the $\gamma$ rays remains a subject of much debate, with an unresolved population of millisecond pulsars (MSPs) or a past transient event at the Galactic center being some of the most popular interpretations Abazajian (2011); Wharton et al. (2012); Hooper et al. (2013); Carlson and Profumo (2014); Petrovic et al. (2015, 2014); Mirabal (2013). Considering both the broad interest in indirect searches for DM, and the current questions raised by the Galactic center excess, it is important to identify new high-density targets for DM annihilation indirect searches.

Nearby galaxies, such as Andromeda (M31) and Triangulum (M33) offer a great opportunity to test the origin of the Galactic center excess and to look for signals of DM particles. They are close enough (approximately 785 kpc and 870 kpc away respectively111https://ned.ipac.caltech.edu) so that their stellar disks and bulges can be resolved as two separate components between radio and X-ray energies which is not possible in our Galaxy, since our bulge is obscured by the bright disk emission. It is worth noting that Andromeda in particular was one of the astrophysical objects where compelling evidence for the existence of DM was first brought forward, and shaped our understanding of the Universe (see e.g. Bertone and Hooper (2018)).

M31 was first detected in $\gamma$ rays by the LAT at a significance of 5.3 standard deviations ($\sigma$) with a marginal spatial extension (significance of $\sim$1.8$\sigma$) Abdo et al. (2010). M31, with its disk, has an angular size of over 3*∘, and is therefore one of the rare nearby galaxies that can be spatially resolved by Fermi-LAT. A more recent analysis Ackermann et al. (2017b) revisited M31 with more than 7 years of Pass 8 observations detected the galaxy at a significance of nearly 10$\sigma$ and confirmed the spatial extension at 4$\sigma$ significance with a size of about 0\mbox{~{}.!!^{\circ}}4. M31 has been detected with a similar size of extension also in Ackermann et al. (2018). Its spectrum is consistent with a power law with an index of $2.4\pm 0.1$. The spatial distribution of the emission is consistent with a uniform-brightness disk over the plane of sky with a radius of 0∘*.4 and no offset from the center of M31. The flux from M31 appears confined to the inner regions of the galaxy and does not fill the disk of the galaxy or extend far from it. Since the spatial morphology of the $\gamma$-ray signal is not compatible with the M31 disk, which is a region rich of gas and star formation activity, the emission probably is not interstellar in origin, unless the energetic particles radiating in $\gamma$ rays do not originate in recent star formation activities. Alternative interpretations are that the emission results from a population of MSPs located in the bulge of M31 Eckner et al. (2018) or from the decay or annihilation of DM particles, similar to what has been proposed to account for the Galactic center excess found in Fermi-LAT observations of the MW.

On the other hand, M33 was not detected significantly by the LAT team in Ackermann et al. (2017b) in an analysis using only $\gamma$ rays with energies above 1 GeV, and has been detected as point-like in the Fermi-LAT 8 year source list222https://fermi.gsfc.nasa.gov/ssc/data/access/lat/fl8y/gll_psc_8year_v5.fit (FL8Y) with a significance of about $4.1\sigma$.

M31 and M33 are thus natural targets for DM indirect detection searches. As a general strategy of this work, we will apply the techniques developed in the search for DM in the LMC and SMC Buckley et al. (2015); Caputo et al. (2016), which are both extended $\gamma$-ray sources as is M31.

In this paper we analyze almost 10 years of Fermi-LAT data, about two more years than previous analysis, and we dedicate our search to any evidence of a possible DM contribution.

In Sec. II, we describe the DM distribution in M31 and M33 and how it relates to searches for indirect signals of DM annihilation or decay. In Sec. III, we discuss the analysis setup and technique that we apply, we show how the DM signal from these galaxies would be detected by the LAT and we describe a search for such a signal in the M31 and M33 regions. Finally in Sec. IV we report the results for the annihilation and decay of DM particles and we conclude in Sec. V.

II M31 and M33 dark matter flux and spatial distribution

The flux spectrum $\textrm{d}\phi/\textrm{d}E_{\gamma}$ of $\gamma$ rays originating in DM decay/annihilation processes can be factored into astrophysics- and particle physics-dependent terms Ullio et al. (2002):

[TABLE]

The quantities in the first parentheses are the DM decay rate $\tau_{\chi}$ or the thermally-averaged annihilation cross section with respect to the velocity distribution of DM particles $\langle\sigma v\rangle$, respectively. Moreover, there is the differential yield of $\gamma$ rays from a single DM annihilation $\textrm{d}N_{\gamma}/\textrm{d}E_{\gamma}$, the mass of the DM particle $m_{\chi}$, and a normalization factor $x$ which is unity if the DM is its own antiparticle and $1/2$ otherwise. All of these depend on the unknown particle physics model that includes DM particles. The typical approach for DM indirect detection searches, as we will follow here, is to set an upper (lower) bound, if no excess is observed, on $\langle\sigma v\rangle$ ($\tau_{\chi}$) as a function of the DM mass $m_{\chi}$ while assuming a particular annihilation (decay) channel and its associated spectrum $\textrm{d}N_{\gamma}/\textrm{d}E_{\gamma}$.

In this paper, we assume $x=1$ and consider the final states $b\bar{b}$ and $\tau^{+}\tau^{-}$, which have been of particular interest given the Galactic center excess. Other sets of Standard Model final states are possible, but have sufficient similarity to the channels selected that bounds can be extrapolated reasonably. In this work, we calculate the spectrum $\textrm{d}N_{\gamma}/\textrm{d}E_{\gamma}$ for each chosen final state and DM mass using a code available as part of the Fermi-LAT ScienceTools.333The DMFitFunction spectral model described at http://fermi.gsfc.nasa.gov/ssc/data/analysis/documentation/Cicerone/Cicerone_Likelihood/Model_Selection.html, see also Ref. Jeltema and Profumo (2008). Note that our implementation does not include electroweak corrections Kachelrieß and Serpico (2007); Berezinsky et al. (2002); Kachelrieß et al. (2009); Ciafaloni and Urbano (2010); Ciafaloni et al. (2011). Such corrections can be important for heavy DM ($m_{\chi}\gtrsim 1$ TeV); in any case, they would increase the resulting flux and, thus, strengthen the resulting bounds Ciafaloni et al. (2011); Cirelli et al. (2011); Ackermann et al. (2015b).

In order to describe experimental results in terms of the particle physics parameters, the astrophysical quantities in the second set of parentheses in Eq. 2 must be known. The integral of the DM density along the line of sight and over a solid angle $\Delta\Omega$ corresponding to the region under study, is known as the $J-$factor (or $D-$factor in the case of decaying DM), and encapsulates the dependence of an indirect detection search on the distribution of DM in the search target. Of particular interest is the case of annihilating DM where the $J-$factor depends on the DM density squared and also implicitly on the inverse distance squared. Hence, targeting nearby overdensities of DM yields larger values of the $J-$factor. Such targets are, thus, very well suited to probe smaller annihilation cross sections $\langle\sigma v\rangle$.

To apply the indirect search pipeline that has been developed to study the DM content of the LMC and SMC, we must first determine the DM density distribution of M31 and M33; that is to calculate their expected $J-$factor ($D-$factor). This task is of a complex nature as $N-$body simulations of the formation and evolution of MW-sized galaxies predict a hierarchical formation scenario. The DM halo of spiral galaxies, like the one of M31 or M33, is expected to form by mergers of small overdensities which are referred to as subhalos. Depending on the particle resolution of the respective $N-$body simulation (Springel et al. (2008); Kuhlen et al. (2008a); Griffen et al. (2016)), around 10 - 20 $\%$ of the mass of a MW-sized galaxy’s DM halo has been found to be present in the form of substructure. An extrapolation of these results to less-massive, and yet unresolved, subhalos seems to predict that in the most extreme scenario about 50$\%$ of a DM halo’s mass stems from substructure. This strongly affects $\gamma$-ray indirect DM searches because a population of DM subhalos can boost the $J-$factor of the parent halo substantially Kuhlen et al. (2008b); Sanchez-Conde and Prada (2014). The $D-$factor of a DM halo, on the other side, is mostly unaffected by the presence of substructure since it grows linear in the DM density. Nonetheless, the region exhibiting the largest $J-$/$D-$factor in a galaxy is its center where the DM density is dominated by the profile of the smooth parent halo. The observed origin of the extended gamma-ray emission from M31 is coinciding with its central region so that we are required to carefully select smooth DM halos for M31 (and M33) that, on one side, cover the full variety of existing DM profiles types and, on the other side, are in accordance with the available stellar data.

As a matter of fact, M31 seems to be the only well-studied galaxy which was argued to require the effect of adiabatic contraction around its central region Seigar et al. (2006); Gnedin et al. (2011). Adiabatic contraction is caused by baryonic physics and gravitational interactions between baryons and DM in galaxies. During the formation of a galaxy, typical processes like gas dissipation, supernova feedback and star formation lead to substantial energy losses of a sizable fraction of galactic baryons which hence fall into the central region of their host galaxy. As first reported in Zel’dovich et al. (1980), these particles deepen the gravitational potential of the galactic center so that the surrounding DM follows the baryonic pull creating a compressed DM halo in the central region. Subsequent hydrodynamical simulations of galaxy formation (Gustafsson et al. (2006); Colin et al. (2006); Gnedin et al. (2011); Tissera et al. (2010); Zemp et al. (2012); Sommer-Larsen and Limousin (2010)) confirm the prediction of adiabatic contraction of DM halo profiles obtained from DM-only simulations. In fact, an adiabatically contracted Navarro-Frenk-White (NFW) profile Navarro et al. (1996, 1997) enhances the $J-$factor in the center of M31 which is remarkable since, following a DM interpretation, a large $J-$factor could be the source of the observed extended gamma-ray emission from M31 in this region. The NFW profile is a particular instance (with $\left(\alpha,\beta,\gamma\right)=\left(1,3,1\right)$) of the general Hernquist-Zhao profile Hernquist (1990); Zhao (1996):

[TABLE]

where $\rho_{s}$ is a density normalization, $r_{s}$ refers to the profiles’ scale radius and $\alpha,\beta$ and $\gamma$ determine the inner and outer slope of $\rho_{{\rm ZHAO}}$ as well as the transition between both regimes. These parameter labels are also used in the definitions of other DM density profiles considered in this analysis.

We adopt the smooth adiabatically contracted NFW profile from Seigar et al. (2006) where it is called “M1 B86”. In detail, we read off the “Halo” mass-to-radius curve in their Fig. 6. Afterwards, we convert it into a radial density profile via $\rho(r)=\left(4\pi r^{2}\right)^{-1}\textrm{d}M/\textrm{d}r$ and interpolate the obtained data points. In fact, the resulting density profile cannot be described by a single set of parameters of the Hernquist-Zhao profile. However, for $r<25.65$ kpc (the profile’s scale radius) it provides an acceptable – albeit far from good – approximation with the parameters $\left(\alpha,\beta,\gamma\right)=\left(0.38,3.84,1.54\right)$.

Alongside this non-standard DM density distribution, we consider two distinct but frequently-used profiles, namely an Einasto profile that provides a better fit to the profile of DM halos derived from $N-$body simulations Einasto (1968); Navarro et al. (2004) (representing the family of cuspy profiles)

[TABLE]

and a Burkert profile Burkert (1996) (representing the family of cored profiles)

[TABLE]

These three DM profiles bracket the range of cosmologically and astrophysically viable DM halo morphologies according to the current understanding of large structure formation and baryonic feedback. As a matter of fact, baryonic physics has been identified not only to be the driving force of an adiabatic contraction of the innermost region of a DM halo but it can also have the opposite effect leading to the formation of a DM core Macciò et al. (2012). The infall of baryons into the center of a DM halo is described to trigger a large number of enhanced star formation periods which each time create a massive outflow of baryons from the central region. The DM follows the baryonic flow as this flow causes a shallower gravitational potential so that the inner cusp is successively washed out by the cycles of baryonic infall and a burst in star formation Mashchenko et al. (2006, 2008); Pontzen and Governato (2012). There have been attempts to implement this kind of baryonic feedback in cosmological simulations of structure formation. They confirm a flattening of the inner DM cusp observed in DM-only simulations Governato et al. (2010); Macciò et al. (2012). However, as for example pointed out in Di Cintio et al. (2014), the impact of baryons on the evolution of a DM halo depends on its properties like the stellar-to-halo mass ratio. Hence, the ultimate effect of baryonic feedback can vary from galaxy to galaxy.

Regarding the Einasto and Burkert profile of M31, we adopt the set of parameters that has been derived in Tamm et al. (2012) to match its available kinematical stellar data. As concerns M33, we make the same distinction between cored and cuspy profiles by selecting a Burkert profile444The admissible Burkert core radii of M33 seem to span a rather wide range Corbelli et al. (2014); Lopez Fune et al. (2017). The difference is most likely attributed to the chosen analysis approach, rotation curve data selection and fitting scheme. and an NFW profile where the respective parameters are taken from Lopez Fune et al. (2017). A summary of the adopted profile parameters is found in Table 1. In case of M33’s NFW profile and M31’s adiabatically contracted NFW profile, we had to calculate $r_{s}$ and $\rho_{s}$ from the given quantities in the respective reference.

Since Eq. 2 involves a line of sight integral, we also have to consider the DM distribution in the MW when pointing towards M31 or M33. Moreover, it is not possible to exclude a priori that the DM halo of the MW has left no traces in the $\gamma$-ray data of Fermi-LAT. Classically, the MW DM halo is fit by either an NFW or a Burkert profile to cover profiles featuring either a central cusp or core. However, the morphological difference between both profiles in the Galactic center is a marginal aspect with respect to the sky position of M31 ((\textrm{RA},\textrm{DEC})=(10\mbox{{}.!!^{\circ}}685,41\mbox{{}.!!^{\circ}}269)) or M33 ((\textrm{RA},\textrm{DEC})=(23\mbox{{}.!!^{\circ}}462,30\mbox{{}.!!^{\circ}}660)) so that we choose the NFW profile from Nesti and Salucci (2013) as the smooth MW DM halo. Comparing the MW $J-$factor within a circular region of interest (ROI) of radius $1^{\circ}$ centered on M31 or M33, the relative difference between a Burkert and an NFW profile is about $15\%$ in both cases. For instance, while the average $J-$factor of M31 from an adiabatically contracted NFW profile inside a circular ROI of $0.4^{\circ}$ (corresponding to the spatial extension of M31’s $\gamma$-ray emission) is about a factor of 100 larger than the respective $J-$factor from the MW NFW halo, the situation reverses not far outside this particular ROI due to the almost perfect isotropy of MW’s $J-$factor (cf. Fig. 1(a)). In another extreme case, namely choosing the same circular ROI of $0.4^{\circ}$ and centering it on M33, the $J-$factor from M33’s DM halo following a Burkert profile is only about $90\%$ of the corresponding $J-$factor due to the MW’s DM halo.

We aim to go beyond the zeroth order $J-$factor estimates by accounting for the effect of substructure inside the smooth halos of the MW, M31 and M33. However, quantitative statements about the net effect of the DM substructure boost are hard to formulate as a precise prediction of the present-day subhalo population properties – e.g., their radial distribution in the host halo or the subhalo survival probability until present time – of MW-sized galaxies remains an objective of ongoing research. To account for those uncertainties, we define two limiting substructure scenarios (MAX and MIN) that model the expected substructure of the galaxies under study according to the extreme cases still allowed by $N-$body simulations. Moreover, we also create a benchmark scenario (MED) that features a DM substructure distribution based on the best-fit parameters of recent observations and numerical simulations of structure formation. To this end, we use the public code CLUMPY Hütten et al. (2019); Bonnivard et al. (2016); Charbonnier et al. (2012) to generate 2D $J-$/$D-$factor sky maps of M31 and M33 as source model for our analysis pipeline as well as MW DM templates as an additional background source.

The parameters governing the substructure distribution in a galaxy that have the largest impact on the expected $J-$factors are:

The index $\alpha_{\textrm{sub}}$ of the subhalo mass function $\textrm{d}n/\textrm{d}M$ which was found to follow a power-law Kuhlen et al. (2008a); Springel et al. (2008),

• -

the fraction of the DM halo mass which is stored in substructure $f_{\textrm{sub}}$,

• -

the minimal mass of DM subhalos $M_{\textrm{min}}$ and

• -

the subhalo concentration parameter $c_{200,\,\textrm{sub}}$ Bullock et al. (2001).

We rely on the most recent model of the concentration parameter of subhalos Moliné et al. (2017). We make use of a developer’s version of CLUMPYv3 which features this concentration model for extended extragalactic objects. This model reports a flattening of the concentration of subhalos towards the low mass tail of the relation and, furthermore, it includes a dependence on the position of the subhalo within its host halo555Neither the release version nor the developer’s version of CLUMPYv3 implement the spatial dependence of the subhalo concentration parameter so that this model can only be considered as a slight improvement with respect to the previous parametrization given in Sanchez-Conde and Prada (2014). Using this description of the subhalo concentration relation implicitely assumes that the DM profiles of subhalos follow an NFW profile which is a fair assumption given the large uncertainty on this quantity as obtained from $N-$body simulations Springel et al. (2008); Kuhlen et al. (2008a); Griffen et al. (2016).

The remaining mutually dependent parameters are chosen as to be consistent with the findings of DM-only $N-$body simulations of MW-sized DM halos. In fact, it has been theoretically established that the minimal subhalo mass $M_{\textrm{min}}$ depends on the particle physics nature of DM so that it might cover orders of magnitude down to values like $10^{-12}\,M_{\odot}$ Binder et al. (2017); Bringmann (2009). In fact, $M_{\textrm{min}}$ cannot be constrained very well even with astronomical data from current-generation instruments. A natural upper limit on $M_{\textrm{min}}$ are the masses of dSphs that have already been detected and are resolved in the MW halo. The faintest of them possess masses as low as $\approx 10^{7}\,M_{\odot}$ Strigari et al. (2007); Walker et al. (2007). From this range, we fix $M_{\textrm{min}}=10^{6}\,M_{\odot}$ in the context of the MIN scenario (being very conservative) and $M_{\textrm{min}}=10^{-12}\,M_{\odot}$ with respect to the MAX scenario (being overly optimistic)666CLUMPY restricts the user to at most $M_{\textrm{min}}=10^{6}\,M_{\odot}$ for extragalactic objects.. Our benchmark case MED assumes a typical value of $M_{\textrm{min}}=10^{-6}\,M_{\odot}$ which is often used in the context of a $\Lambda$CDM cosmology with cold, thermal WIMP DM Green et al. (2004). From these values of the minimal subhalo mass we can infer the expected fraction of subhalos in a galaxy for a given subhalo mass function index $\alpha_{\textrm{sub}}$ using the reported behavior found, e.g., in the Aquarius project or the Via Lactea simulation Kuhlen et al. (2008a); Springel et al. (2008). We summarize the definitions of our three substructure scenarios and their assignments to the smooth DM halos in Table 1. Those assignments were made in order to bracket the theoretically and observationally allowed $J-$factors from M31 and M33 according to the combination of smooth halo and substructure parameter set.

To complete our model definition, we assume that the radial distribution of subhalos follows the smooth DM density profile of their host halo. As a consequence of this choice, subhalos are expected to populate even the central region of the host DM halo despite the impact of tidal forces and other effects, like baryonic physics, in this part of the galaxy. Nonetheless, DM-only $N-$body simulations, e.g. the aforementioned Aquarius and Via Lactea II simulations Kuhlen et al. (2008a); Springel et al. (2008), poorly constrain the subhalo population in this particular region Pieri et al. (2011). In fact, while certain simulations including baryonic feedback during the galaxy formation process reveal a strong depletion of subhalos in a large volume around the center of a galaxy Bullock and Boylan-Kolchin (2017), there are opposing works like van den Bosch and Ogiya (2018); Peñarrubia et al. (2010) arguing that the observed tidal disruption of subhalos is a numerical artifact due to the particle resolution of the simulations. In the context of this ongoing debate, our choice seems as justified as any other.

Having established a full description of our models and substructure parameter sets, we show in Fig. 1 the generated radial profiles of $J-$/$D-$factors. To stress it again, in the case of decaying DM, we only use the here discussed smooth DM profiles since the $D-$factor is proportional to the DM density and does not feel the boost due to subhalos.

III Data analysis

The analysis technique we apply in this paper closely follows those used in the DM searches from dwarf spheroidal galaxies Ackermann et al. (2015a); Drlica-Wagner et al. (2015), Smith High-Velocity Cloud Drlica-Wagner et al. (2014), and LMC and SMC Buckley et al. (2015); Caputo et al. (2016). In the next sections we describe the data selection and the different steps of the pipeline.

III.1 Data selection

We use 115 months of Fermi-LAT data (from 2008 August 4 to 2018 March 4) selecting Pass 8 SOURCE class events and using the corresponding instrument response functions (IRFs) P8R3_SOURCE_V2. We choose an energy range from 300 MeV to 1 TeV and select events with reconstructed directions within a 16\mbox{~!!{}^{\circ}}\times16$~\!\!{}^{\circ}$ region of interest (ROI) centered at the infrared position of M31 and M33. We test also different lower bound values (e.g., 100 MeV) for the energy range and different sizes of the ROI to see how the results are affected by these parameters (see Sec. IV.2). We bin the data using 8 energy bins per decade in energy and 0$~{}.\!\!^{\circ}$08 pixel size.

The Pass8 event reconstruction and selection called Pass 8 introduces a generalization of the conversion type classification in the form of event types. PSF event types are event-level quantities indicating the quality of the reconstructed direction. The data are divided into quartiles, from the lowest quality quartile (PSF0) to the best quality quartile (PSF3) Atwood et al. (2013).

We apply zenith angle cuts to the data in order to reduce the contamination from the low-energy Earth limb emission. We select for $E=[0.1,0.3]$ GeV PSF2 and PSF3 event types with zenith angles $\theta_{z}>90^{\circ}$, for $E=[0.3,1.0]$ PSF1, PSF2 and PSF3 event types with zenith angles larger than $\theta_{z}<100^{\circ}$ and finally above 1 GeV we keep all PSF types with zenith angles $\theta_{z}>105^{\circ}$. We apply the same cuts used in the construction of FL8Y (and 4FGL, see later in the text) source list; these reduce the contribution of the Earth limb at that zenith angles to less than $10\%$ of the total background. See Table 2 for a summary of the analysis setup.

We construct a background model of each region that includes the FL8Y sources in the region, as well as an interstellar emission model and an isotropic emission template. Very recently, the Fermi-LAT 8-year Point Source Catalog (4FGL) has been created using a new IEM and isotropic templates Fer (2019). The FL8Y and 4FGL catalogs have been created with the same years of data and no significant difference in the characteristics of the sources present in M31 and M33 ROIs are present thus we do not find any relevant difference in our result by using the 4FGL instead of FL8Y (see Sec. IV.2). Specifically, we use the interstellar emission model (IEM) released with Pass 8 data (Acero et al., 2016a) (i.e., gll_iem_v06.fits) since this is the model routinely used in Pass 8 analyses. We will label this IEM as Official (Off). This model is derived by performing a template fitting to Fermi-LAT $\gamma$-ray data. It is thus based on the spatial correlations between $\gamma$-ray data and a linear combination of gas and inverse Compton scattering maps. This model contains patches to account for extended excess emissions of unknown origin. However, the M31 and M33 regions do not contain any of these patch components. We use for the isotropic emission the template associated to this IEM (iso_P8R3_SOURCE_V2.txt)777For descriptions of these templates, see http://fermi.gsfc.nasa.gov/ssc/data/access/lat/BackgroundModels.html..

To approximately study the systematic uncertainties from the mismodelling of the diffuse emission, we also run our analysis using the 8 alternative IEM models and corresponding isotropic templates used in the first Fermi-LAT supernova remnant (SNR) catalog Acero et al. (2016b). These models were generated by varying the cosmic-ray (CR) source distribution, height of the CR propagation halo, and HI spin temperature in order to test the effect of the choice of the IEM in the flux and spatial distribution of SNRs. These 8 models, which are all based on the GALPROP888http://galprop.stanford.edu/ CR propagation and interaction code, have been used in the SNR catalog to explore the systematic effects on SNRs fitted properties, including the size and morphology of the extension, caused by IEM modeling. We will label these models as alternatives (Alt).

It is important to stress that the Off and Alt IEMs have been designed to model the diffuse background for analysis of point and small extended sources. Because they are fit to the data, they are not suited for studies of very extended sources and/or large-scale diffuse emissions. Since both M31 and M33 signal are extended at $\leq$ degree scales, and are not correlated with any other diffuse template, these diffuse models are applicable for our analysis.

In fact, the extended sources studied in the Fermi-LAT SNR catalog Acero et al. (2016b) have similar spatial extension as M31 and M33. Finally, we note that the Off IEM and isotropic templates have been routinely used in previous DM analysis from Fermi-LAT. We also stress that if the excess signal would have been found in our ROI with any of these models, then a dedicated diffuse analysis would be required to determine its properties, but it will be shown below that this is not the case here999This approach is complementary to the recent work Karwin et al. (2019), in which an almost isotropic emission on the 10 degree scales from M31, degenerate with the isotropic component of the MW was studied. In that case a careful modeling of the diffuse emission was necessary and was indeed undertaken in that work..

We employ FermiPy to perform our analysis of Fermi-LAT data. FermiPy is a Python package that automates analyses with the Fermi Science Tools (Wood et al., 2017)101010See http://fermipy.readthedocs.io/en/latest/.111111We use version 17.3 of FermiPy and 11-07-00 of the ScienceTools.. We will explain in detail the analysis pipeline in the next sections.

III.2 Baseline fit

We first perform a broadband fit to our ROIs using the sources from the catalog, the Off IEM and isotropic template. The results of this baseline fit (without the DM template) will be used as an input to the DM dedicated analysis described below. The size of the ROI has been taken to be much larger than the DM contribution. Indeed, we see from Fig. 1 that for an angular distance $>5^{\circ}$ from the center of the ROIs the contribution to the $J$ and $D$ factors are negligible. During the broadband fit the spectral energy distribution (SED) parameters of all the point sources in the ROI, the normalization and spectral index of the IEM and the normalization of the isotropic templates are free to vary. At this stage in the analysis M31 and M33 are modeled as point-like sources.

Then, we relocalize all the sources in the ROIs, including M31 and M33. Since we are using more years of data than FL8Y, we identify new sources detected with a test statistic141414The $TS$ is defined as twice the difference in maximum log-likelihood between the null hypothesis (i.e., no source present) and the test hypothesis: $TS=2(\log\mathcal{L}_{\rm test}-\log\mathcal{L}_{\rm null})$ Mattox et al. (1996). $TS>25$. In the last step of the procedure to define the baseline model for the ROI we refine the astrophysical model for M31 and M33. M31 has been detected in Ackermann et al. (2017b), using about 7 years of data above 100 MeV, as extended with $TS_{\rm{EXT}}=16$ and size 0\mbox{{}.!!^{\circ}}38\pm 0\mbox{{}.!!^{\circ}}05. We re-run the extension analysis at $E>0.3$ GeV finding that $TS_{\rm{EXT}}=13.6$ and the size (i.e., the 68% containment radius) is 0\mbox{{}.!!^{\circ}}33\pm 0\mbox{{}.!!^{\circ}}04 for a disk template and $TS_{\rm{EXT}}=12.8$ with a size of 0\mbox{{}.!!^{\circ}}42\pm 0\mbox{{}.!!^{\circ}}10 for a Gaussian template. We also run this analysis for $E>0.1/0.5/1$ GeV to see if there is an energy dependence of the spatial extension. We report the results for the detection and the spatial morphology of M31 in Table 3 for the disk and Gaussian spatial templates. The sizes of extension for M31 in the different energy ranges are all compatible within $1\sigma$. We significantly detect M31 as extended also at $E>0.1$ GeV with a similar size of extension and $TS$. This justifies our choice of using an extended template that is energy independent. Since the disk morphology is slightly preferred we use this geometry for M31 as the benchmark case in the rest of our analysis (see Sec. IV.5, IV.4 and IV.3).

M33 is detected as a point-like source ($TS_{\rm{EXT}}\approx 0$) with a $TS=41.9$ ($TS=39.4$) with the analysis performed in the energy range $E\in[0.1,1000]$ GeV ($E\in[0.3,1000]$ GeV).

Fig. 2 shows the $TS$ map of the ROI of M31 and M33 without these galaxies included in the model. It is clear from these plots that M31 is much brighter than M33 and that it has an extension of the order of the size detected in our analysis. In the same figure we also show the $TS$ map for the baseline model, i.e. with all sources included and with their positions refined. There are no significant residuals in these maps meaning that the baseline model is an appropriate fit to the two ROIs. The highest $TS$ peaks are of the order of $2-3\sigma$ significance. We use this model as a baseline for the search of a DM signal.

We also run the baseline analysis for M31 and M33 using the 8 Alt IEM models and corresponding isotropic templates used in the first Fermi-LAT SNR catalog Acero et al. (2016b). Depending on the Alt IEMs, M31 is detected with $TS\in[84-110]$, $TS_{\rm{EXT}}\in[9-15]$, disk size \theta_{\rm{EXT}}=[0\mbox{{}.!!^{\circ}}30,0\mbox{{}.!!^{\circ}}45] and Gaussian template size \theta_{\rm{EXT}}\in[0\mbox{{}.!!^{\circ}}32,0\mbox{{}.!!^{\circ}}43] (see Table 4). In short, no significant changes are found by using different IEMs. Similarly, M33 is detected as a point source with all the IEMs ($TS_{\rm{EXT}}\approx 0$) with a $TS$ in the range $34-42$ (see Table 4).

We test a power-law shape for the SED and also a log-parabola (LP) or power-law with an exponential cutoff (PLEC). There is no significant preference for the LP or PLEC over the simple power-law so we decide to assume this SED shape in the rest of the paper when we consider the models used in this section.

We test a power-law shape for the SED and also a log-parabola (LP) or power-law with an exponential cutoff (PLEC). There is no significant preference for the LP or PLEC over the simple power-law so we decide to assume this SED shape in the rest of the paper when we consider the models used in this section.

III.3 Baseline fit using astrophysical models from the observations of M31 and M33 in other wavelengths

The baseline fit reported in the previous section, with M31 modeled with a disk template and M33 with a point-like source, is tuned directly on $\gamma$-ray data. It is thus a phenomenological way to explain the $\gamma$-ray emission from M31 and M33 without any direct relation with what is observed from these galaxies in other wavelengths. $\gamma$ rays are produced predominantly by the so-called interstellar emission that is traced by atomic gas density, radio and infrared emissions. Therefore, we used maps derived from observations in other wavelengths as templates for the astrophysical components of the $\gamma$-ray emission from M31 and M33. For M31 we use the Herschel/PACS map at 160 $\mu$m (which traces the star formation), Spitzer/IRAC map at 3.6$\mu$m (that traces old stellar population), and an atomic gas column density $N_{H}$ map from Braun et al. (2009) (which traces gas densities). On the other hand, for M33 we consider the 2MASS infrared map at 2.2$\mu$m, the Spitzer/IRAC infrared maps at 24$\mu$m and at 160$\mu$m, and the Mayall optical map at $3552$ A. The analysis used here is the same as in Sec. III.2 and we assume for the SED shape of these templates a simple power-law. We report the values of the $TS$ for the different components in comparison with the result of the point source, disk and Gaussian templates in Tab. 5.

Among the templates considered for M31, the $N_{H}$ map yields the fit with the lowest likelihood and $TS$. Indeed, we see in Fig. 3 that the $TS$ map for the case with $N_{H}$ template has more residuals than the one with Spitzer/IRAC template. The reason is that the $N_{H}$ map traces the disk of M31 that is extended about $3^{\circ}$ across the sky while the $\gamma$-ray emission detected by Fermi-LAT is concentrated within about $0.^{\circ}4$ with a spherical symmetry.

For the same number of degrees of freedom, the Herschel/PACS provides a slightly better fit that however it is not favored compared to a simple point source at the center of M31. Finally the Spitzer/IRAC map provides the best fit to the data among the templates considered from other wavelengths. The Spitzer/IRAC infrared map is the most similar to those obtained with geometrical models because it is dominated by the bulge component and it is thus similar with what is observed in $\gamma$-ray data. Therefore, the $\gamma$-ray emission that we detect with Fermi-LAT from the M31 direction is dominated by the bulge emission. Indeed, the stellar bulge of M31 has a total mass of about $3.1\cdot 10^{10}$ $M_{\odot}$ and a size of about 5 kpc Tamm et al. (2012). The size is perfectly compatible with the extension that we measure in $\gamma$ rays. Moreover, the M31 stellar bulge is a factor of about 5 more massive than the MW one. These results are compatible with what has been presented in Ackermann et al. (2017b).

The astrophysical templates we try for M33 give almost all the same significance which is the of the same order of the point source scenario. Indeed, the infrared and optical emission observed from this galaxy have an extension of about 0\mbox{~{}.!!^{\circ}}2 that is of the order of the point source emission. Since the simple point source emission gives the same $TS$ as the templates derived from observations in other wavelength we use directly this model in the rest of the analysis.

III.4 M31 and M33 DM morphology and correlations with other background sources

As we have seen in the previous sections we detect with Fermi-LAT a $\gamma$-ray signal that is well fit with a radial disk of size 0\mbox{~{}.!!^{\circ}}33 for M31 and with a point source for M33. We have demonstrated that the $\gamma$-ray flux in the direction of M31 is compatible with the emission from the stellar bulge while M33 with infrared and optical maps.

A possible contribution might also come from DM so it is important to understand the Fermi-LAT sensitivity to such a signal. In particular we are interested in determining the spatial morphology of a putative DM emission (i.e., determine whether it would be detected as point like or extended source and in the latter case estimate the size of extension) and calculating the correlations between this component and the other background sources. We use Fermipy to simulate a DM signal from M31 and M33 DM templates and we derive its size and spatial morphology. We take the baseline model from Sec. III.2, we remove M31 and M33 sources from the model and we add a DM contribution using the templates reported in Sec. II. We simulate a DM signal that would give a detection at $9\sigma$ significance with a power-law SED with index $2.0$. We use such a simple power-law SED, as we are interested here only in finding the spatial morphology of the DM signal. We show in Fig. 4 the $TS$ map for the $\gamma$-ray emission from this DM signal for M31 and M33 for the MED DM model. A very similar $TS$ map is found when using the MIN or the MAX DM model. Then, we fit this excess with an extended source finding that it is well fitted with a Gaussian template with size of about 0\mbox{{}.!!^{\circ}}5 and $TS_{\rm{EXT}}=30$ for M31 and 0\mbox{{}.!!^{\circ}}9 and $TS_{\rm{EXT}}=25$ for M33. The Gaussian template has about the same likelihood value as the disk one. The correlation coefficients between the M31 DM template normalization (spectral index) and the isotropic and IEM template normalizations are -0.07 (-0.14) and -0.23 (-0.12), respectively. Instead, the correlation coefficients between the M33 DM template normalization (spectral index) and the isotropic and IEM template normalizations are -0.04 (-0.19) and -0.16 (-0.15), respectively. Therefore, no significant correlations are present between the DM template of M31 and M33 and IEM and isotropic components. The correlation coefficients between the DM component normalization and the M31 disk normalization is -0.87; instead with the M33 point source SED normalization it is -0.61. Therefore, SED parameters of the M31 and M33 DM templates are correlated with the SED of the disk template for M31 and M33, respectively. We consider these correlations in our analysis as we will explain in Sec. III.5. We find similar results using the MIN or MAX DM models.

In the line of sight of M31 and M33 a contribution of $\gamma$ rays could also come from DM present in the MW. As shown in Fig. 1, the $J$-factor for the MW is much smaller than the M31 and M33 component in the inner few degrees from M31. Moreover, the DM MW signal is almost isotropic with a variation of about a factor of $10\%$ across the M31 and M33 ROIs. We test a possible effect of the presence of this additional DM component by taking the simulation done before, which includes the DM M31 component, and adding also the contribution of DM from the MW. We then ran a fit and found that the MW contribution is almost completely absorbed by the isotropic template. Moreover, there is no correlation between the MW and the M31 DM templates. Therefore, we decide in the rest of our analysis to not add this component to the model.

III.5 DM SED and significance

At this stage of the analysis we add to the baseline model, which includes the template for astrophysical emission from M31/M33 (i.e., a disk template for M31 and a point like source for M33), the DM template. We run a fit and then compute the likelihood profile as a function of energy and energy flux of DM. We scan in each energy bin the likelihood as a function of the flux normalization for the assumed DM signal which is specified by the choice of decay or annihilation, by the channel and the DM mass. In the rest of the paper we will consider the decay and annihilation into $b\bar{b}$ and $\tau^{+}\tau^{-}$. The case of $b\bar{b}$ is representative of hadronic channels such as quarks and gauge bosons while $\tau^{+}\tau^{-}$ represents the leptonic channels $\mu^{+}\mu^{-}$ and $e^{+}e^{-}$. For this bin-by-bin scan, we fix the SED parameters of the sources that have an absolute value of the correlation parameter smaller than 0.10. By analyzing each energy bin separately, we avoid selecting a single spectral shape to span the entire energy range at the expense of introducing additional degrees of freedom into the fit. For the fit in any given bin, the only free parameter describing the DM component is the normalization. In Fig. 5 we report the SED of the DM template for M31 and M33 with the MED model.

In the final step of the fitting procedure we convert the bin-by-bin likelihood curve in flux into a likelihood curve in $\langle\sigma v\rangle$ for each spatial profile and annihilation or decay channel, which determines the spectrum. We scan DM masses ($m_{\chi}$) for 5–10000 GeV (when kinematically allowed in the annihilation or decay channel under consideration), and the pair-annihilation or decay final states $\tau^{+}\tau^{-}$ and $b\bar{b}$. For each DM spectrum, we extract the expected flux, $F_{j}$, in each energy bin and calculate the likelihood of observing that flux value. The log-likelihood in each energy bin is summed to get the log-likelihood curve, defined as:

[TABLE]

where $\mathcal{L}$ is the likelihood, $j$ runs over the energy bins of Fermi-LAT data ($\mathcal{D}$), $\mu$ are the DM parameters ($\langle\sigma v\rangle$ or $\tau$ and $m_{\chi}$) and $\theta$ are all the other parameters in the background model, i.e. the nuisance parameters.

Therefore, the DM SED has all the information needed to find if the presence of DM is significant or not for any possible DM annihilation channel and mass. Indeed, we can choose a DM annihilation or decay channel and convert the SED into the likelihood profile as a function of the DM $\langle\sigma v\rangle$ and $m_{\chi}$.

We can also evaluate the significance of the DM hypothesis using the $TS$ defined as:

[TABLE]

where $\mathcal{L}_{\mathrm{null}}$ is the likelihood for the null signal of DM, $\mathcal{L}$ is the likelihood for the presence of DM. For the energy bins up to about 10 GeV the statistics are large enough that Chernoff’s theorem applies, and we expect the $TS$-distribution to follow a $\chi^{2}$ distribution Chernoff (1938). At higher energies, the counts per bin are in the Poisson regime and the $\chi^{2}$ distribution moderately over-predicts the number of high $TS$ trials observed in simulated data.

From $\mathrm{ln}~{}\mathcal{L}(\mu,\theta|\mathcal{D})$ we can evaluate the one-sided 95% confidence level (CL) exclusion limit on the flux as the point at which the p-value for a $\chi^{2}$ distribution with 1 degree of freedom is 0.05 when we take the maximum likelihood estimate as the null hypothesis. That is, the 95% CL upper limit on the flux assigned to DM is the value at which the log-likelihood decreases by 1.35 with respect to its maximum value. An example of this is shown in Fig. 6 for $b\bar{b}$ annihilation channel for M31 and M33.

IV Results

In Sec. IV.1, we first show a model-independent analysis where we search for an excess in M31 and M33 ROIs without using any DM model. In Sec. IV.2 we show how the limits on the DM annihilation cross section change by considering a different energy range for the analysis, size of the ROI, assuming different templates for the astrophysical emission of M31 and M33. Then, in Sec. IV.3 we use 100 simulations to calculate the expected limits on DM in the null signal hypothesis. In Sec. IV.4 we assume that all signal comes from DM and we derive best-fit contours in cross section and mass parameters space. As there is no compelling evidence that the emission from these galaxies is solely due to DM and in fact the best fit region is in tension with the DM search from other targets (e.g., dSphs or the Galactic Center), in Sec. IV.5 we set DM limits. We do this for the $b\bar{b}$ and $\tau^{+}\tau^{-}$ channels and for the MAX, MED and the MIN DM distributions considered in Section II. These channels were previously considered in the dwarf spheroidal analyses Ackermann et al. (2015a); Drlica-Wagner et al. (2015).

IV.1 DM model independent search for an excess in M31 and M33 ROIs

As discussed in Sec. II the DM spatial distribution can vary significantly by assuming a MAX, MED and the MIN model or considering annihilation or decay of DM particles. In addition to this, the DM SED is also uncertain and can vary for the different annihilation or decay channels.

In this section we describe a search for a radial dependent excess by adding to the M31 and M33 ROIs three uniform annuli with radial shapes: $r\in[0.4^{\circ},3.5^{\circ}]$, $r\in[3.5^{\circ},6.0^{\circ}]$ and $r\in[6.0^{\circ},8.0^{\circ}]$. For this analysis we use a ROI width of $20^{\circ}$ to avoid edge effects in the farthest annulus and we select an energy range between $0.3-1000$ GeV. The analysis pipeline is the same presented for the baseline fit in Sec. III.2 where we leave free to vary in the fit the SED parameters of point sources, the normalization of the isotropic template and the normalization and slope of the IEM. We model the emission from M31/M33 with templates described in Sec. III.2.

We show in Fig. 7 the count spectrum with the residuals and the $TS$ map for this fit to M31 and M33 ROIs using the Off IEM and isotropic templates. The change in normalization of the isotropic template has a best fit value of 1.00 (0.914) while the change of normalization and slope of the IEM are 1.010 (1.048) and $-0.03$ (-0.02) for M31 (M33). Therefore, for both M31 and M33 the deviation of the isotropic and IEM SED parameters from their input values are minimal.

Now we add to the model the three annuli and we re-do the fit. The results are reported in Table 6 where we can see that the $TS$ for all the three annuli is very small. We run the same analysis for M31 also with the 8 Alt IEMs finding even tighter constraints. The $TS$ is $\approx 0$ for all the three annuli with upper limits for the flux of $3-4\cdot 10^{-10}$ ph/cm2/s, $2.5-3.5\cdot 10^{-10}$ ph/cm2/s and $2.5-3.5\cdot 10^{-9}$ ph/cm2/s for the first, second and third annulus, respectively, and depending on the Alt IEM considered.

In Fig. 8 we show the intensity (the flux divided by the solid angle) upper limits found above together with the intensity of the DM contribution for MAX, MED and MIN DM models. We report the case where the DM flux fits the astrophysical contribution from M31 taken as a disk template and M33 as a point source. These DM intensity profiles are in tension with the upper limits for the annuli derived in this section. This analysis shows that if DM contributes entirely to the $\gamma$-ray emission from these two galaxies we should be able to detect a signal also in the outer region of M31 and M33, i.e. for $r\geq 2^{\circ}$, which would be absorbed by the annuli. These results might change considering different assumptions for the DM distribution in M31 and M33. We also display the DM intensity profiles for the upper limits of the flux reported in Sec. IV.5. These contributions are compatible with the upper limits found in this section since they have to be considered as upper limits.

A possible weak point of this analysis is that we are leaving free to vary also the normalization of the isotropic template. If residuals are present in the ROI and their spatial distribution is mostly isotropic they can be absorbed by the isotropic template. Therefore, we decide to change the analysis made before by fixing the isotropic normalization in the fit to the value found in a control region. We choose the center of the control region to be at the same longitude of the M31 and M33 ROIs but at a latitude $20^{\circ}$ below. In this way the contribution of the isotropic is higher and we are able to constrain effectively its normalization.

We find a best fit value for the isotropic normalization of 0.866 (0.889) for M31 (M33) in the control regions and we apply these values to the M31 and M33 ROIs fixing this parameter in the fit. Even if we fix the isotropic in the fit, the residuals in the count map are still basically compatible to 0 in the entire energy range, similar to what is shown in Fig. 7. Also the ROI $TS$ maps remain almost unchanged with no significant larger-scale residuals.

The results with the three uniform annuli are reported in Table 6. The $TS$ for the annuli is 5 (0) for the inner, 0 (0) for the second and 20 (0) for the outermost. The upper limits for the flux are slightly larger than in the previous case where we leave the isotropic template free to vary. We run the same analysis for M31 also with the 8 Alt IEMs finding even tighter constraints. The $TS$ is $\approx 0$ for the inner two annuli and between $5-10$ for the outermost with upper limits for the flux of $5-10\cdot 10^{-10}$ ph/cm2/s, $4-7\cdot 10^{-10}$ ph/cm2/s and $3-5\cdot 10^{-9}$ ph/cm2/s for the first, second and third annulus, respectively.

IV.2 Systematics in the DM results

Different choices in the details of the data analysis can modify at different extent the results for our DM search. In this section we report the study of the change in the results for the upper limits on $\langle\sigma v\rangle$ with respect to the lower bound of the energy range and ROI width considered in our analysis, size and model of the astrophysical emission of M31 and M33, prior on $J$ and spatial distribution for the DM component. Other assumptions might change the results by negligible factors. The results are collected in Fig. 9 for M31 (for M33 we find very similar results). The differences are calculated with respect to the baseline setup which is given by the following choices: energy range $E\in[0.3,1000]$ GeV, ROI width $16^{\circ}\times 16^{\circ}$, disk template with size 0\mbox{~{}.!!^{\circ}}33, uncertainty for the $J$ factor $\log_{10}\sigma_{J}=0.2$ and Off IEM.

The choice of the lower bound of the energy range considered in the analysis affects the limits for $m_{\chi}>100$ GeV as an overall normalization. For these masses, the results for $E>0.1$ and $>0.3$ GeV differ only by a factor of about 15% with the limits for the former which are better. For $m_{\chi}<100$ GeV the limits found with $E>0.1$ GeV become larger than the one for $E>0.3$ GeV because very small residuals at low energies are absorbed by the DM template. On the other hand the limits for $\langle\sigma v\rangle$ derived assuming $E>0.5$ and $>1$ GeV are worse by about a factor 2 and 3, respectively. Since the limits found with $E>0.3$ GeV are the tightest and the PSF and acceptance are much better at 0.3 GeV than at 0.1 GeV we decide to use this energy range in the rest of the analysis. This choice is also motivated by the fact that the significance for the detection of M31 and M33 and the $TS_{\rm{EXT}}$ is about the same order for $E>0.1$ and $>0.3$ GeV (see Table 3).

We found in Sec. III.2 that in the baseline model M31 is extended with a size for a uniform disk template of 0\mbox{{}.!!^{\circ}}33\pm 0\mbox{{}.!!^{\circ}}04. During the search of a DM contribution in our pipeline we do not vary the size of the disk component. We run our search for DM with a size of the disk template modified by $\pm 1\sigma_{\rm{EXT}}$ or decreasing it by $3\sigma_{\rm{EXT}}$ from the best-fit and assuming that $\sigma_{\rm{EXT}}$ is the $1\sigma$ error for the extension size. These changes in the size of the astrophysical emission of M31 affects by a negligible contribution the upper limits for $\langle\sigma v\rangle$.

The disk template is a phenomenological model created to fit the $\gamma$-ray emission from M31 and is tuned directly on Fermi-LAT data. Employing this geometrical template could hide part of the DM emission. We perform the analysis by substituting the disk template with the following templates motivated by observations of M31 in other wavelengths: the Herschel/PACS map at 160 $\mu$m, Spitzer/IRAC map at 3.6$\mu$m, and the atomic gas column density $N_{H}$ map from Braun et al. (2009). The results for the DM search is that we find no evidence for a DM contribution in any annihilation or decay channel and the upper limits on $\langle\sigma v\rangle$ are larger for the case with the $N_{H}$ (Herschel/PACS) map by a factor of about $2-3$ ($1.5-2$) for $m_{\chi}\in[10,1000]$ GeV. This is due to the fact that by using the $N_{H}$ or the Herschel/PACS templates the $TS$ of the astrophysical component decreases and more residuals remain in the M31 region that are partially absorbed by the DM component (see Fig. 3). On the other hand the results found with the Spitzer/IRAC template are very similar to the baseline model. Indeed, the $TS$ derived with this infrared map is close to the one obtained with the disk template (see Table 5).

Our benchmark DM spatial templates include a spatially extended map with an uncertainty for the $J$ factor that is 0.20 in $\log_{10}$ units. We test how much the limits for DM change by assuming a point-like DM template and with $\log_{10}\sigma_{J}=0$ and $\log_{10}\sigma_{J}=0.35$. These cases embed the uncertainty in the $J$ factor of the main halo of M31 and M33 found by assuming different functions for the DM distribution, which are $\log_{10}\sigma_{J}=0.35$ for M31 and $\log_{10}\sigma_{J}=0.25$ for M33 Tamm et al. (2012); Lopez Fune et al. (2017). The upper limits decrease (increase) by a factor of about 15% (25%) using $\log_{10}\sigma_{J}=0$ ($\log_{10}\sigma_{J}=0.35$). These differences are subdominant with respect to the ones reported above (e.g. the spatial template for M31). Changing the DM template into a point-like morphology strengthens the limits by a maximum of about a factor of 8 for $m_{\chi}>300$ GeV and of about a factor of 4 for $m_{\chi}\sim 10$ GeV. However, we know that a point-like DM template is not physically motivated for such close galaxies. Choosing $\log_{10}\sigma_{J}=0$ decreases the limits by a normalization factor of the order of 15%.

Moreover, changing the ROI width to $14^{\circ}\times 14^{\circ}$ or $20^{\circ}\times 20^{\circ}$ changes by a negligible amount the limits for $m_{\chi}>100$ GeV while at smaller masses the choice of a smaller ROI width can give significantly larger limits because of residuals generated by edging effects that are absorbed as DM signal at low energy.

Finally, we run the analysis with the Off, the Alt IEMs and using the newest 4FGL catalog and IEM and isotropic templates Fer (2019). The results change by at most $30\%$ between the Off and Alt IEMs while the 4FGL gives a difference of at most 50% at $m_{\chi}=200$ GeV. Even if the limits with the 4FGL catalog are higher than the other cases the significance for the presence of DM is still negligible.

Similar conclusions are also valid for the $\tau^{+}\tau^{-}$ annihilation channel and for the decay case.

IV.3 Null and injected signal simulations

The pipeline that we employ in this paper can also be used to perform simulations. In particular simulations are generated by Fermipy that take the source model and randomize it with Poisson statistics. This method is much faster than the tool gtobssim which is included in the fermitools151515https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/help/gtobssim.txt and it is usually used for the same scope. We consider here two types of simulations to validate our analysis pipeline: null and injected signal simulations.

Null simulations are made by taking the model from the baseline fit, i.e. without the DM contribution, and simulating the ROI. Then we run the search for DM on the simulated data. The goal of these simulations is to calculate the expected limits in the absence of any DM signal. We run the null signal simulations on both M31 and M33, for the $b\bar{b}$ and $\tau^{+}\tau^{-}$ channels and for MIN, MED and MAX DM models. We show in Fig. 10 the upper limits for $\langle\sigma v\rangle$ that we derive for 100 simulations of the null signal for M31 and M33. This is done for the $b\bar{b}$ annihilation channel and for the MED model. The plots show the median upper limits and the $95\%$ and $68\%$ containment bands over the 100 simulations. As expected, the search for DM with these simulations gives a $TS\sim 0$. Therefore, we show the results in form of upper limits for $\langle\sigma v\rangle$. The median is well contained in the $68\%$ containment bands and the limits rule out the thermal cross section for $m_{\chi}<50$ GeV for M31 and $m_{\chi}<20$ GeV for M33.

On the other hand, for the injected signal simulations a DM signal for a specific annihilation channel, cross section and mass is added to the model. Then, Fermipy generates the simulated data, which is analyzed in the same way as the actual data. We perform these simulations to verify that our pipeline is able to recover an injected signal. We choose to inject in the M31 ROI a DM signal with $m_{\chi}=100$ GeV, $\langle\sigma v\rangle=10^{-25}$ cm2/s and for an annihilation channel $b\bar{b}$. We run 100 simulations and we use the MED DM distribution. This signal can be detected at most with a $TS$ of 15. The contour plot for the cross section and DM mass is reported for this simulation in Fig. 10. The best fit cross section and DM mass is perfectly compatible with the characteristics of the injected signal. Since the $TS$ for detection is below 25 we decide to calculate upper limits for each simulation. In Fig. 11 we show the median and the $95\%$ and $68\%$ containment bands over the 100 simulations. These limits are consistent with the cross section of the injected signal demonstrating once again that our pipeline is able to recover an injected DM signal.

IV.4 DM-only interpretation

In this section we make the assumption that the $\gamma$-ray emission in the direction of M31 and M33 is entirely given by DM particle interactions. In order to do so, we remove from the source model the astrophysical model for the $\gamma$-ray emission from M31 and M33, we include the DM template and we run the pipeline described in Sec. III.

Table 7 and 8 contain the best fit and $1\sigma$ errors for the DM mass and cross section for $b\bar{b}$ annihilation channel and for M31 and M33 respectively. We report these results using only the main DM halo (SH) or including also DM substructures (SHS) and for each of these cases we try the MIN, MED and MAX models for the DM distribution.

We first focus on the results found for M31. The $TS$ for the presence of DM for the $b\bar{b}$ annihilation channel, is 27 for the MED DM model in the SHS case and 55 for the Einasto profile in case of SH. This DM candidate has a mass of 30 (20) GeV and a cross section of $3.2\cdot 10^{-26}$ cm2/s ($5.2\cdot 10^{-26}$ cm2/s) for SHS (SH) case. The MIN and MAX DM distributions for SHS and the Burkert and Adiabatic for the SH case provide similar significances and DM masses but larger and lower values for the best-fit cross sections. On the other hand, for the $\tau^{+}\tau^{-}$ annihilation channel the best fit mass is at the lower limit of the DM mass considered in our analysis, i.e. 5 GeV, and a cross section of about $1.5\cdot 10^{-26}$ cm3/s for the MED model. The decay scenario provides a very low significance for the Burkert and Einasto DM profiles and both $b\bar{b}$ and $\tau^{+}\tau^{-}$ channels with $TS$ values of the order of $5-8$. The only case that provides a large significance is with the Adiabatic DM profile for which in the $b\bar{b}$ channel the best fit is $m_{\chi}=30$ GeV and $\tau=6.5\cdot 10^{24}$ s for a $TS=21$ while for the $\tau^{+}\tau^{-}$ channel $m_{\chi}\leq 5$ GeV and $\tau\leq 9.9\cdot 10^{24}$ s for a $TS=20$. Only the Adiabatic DM profile provides a high significance because its spatial profile is similar to the $\gamma$-ray signal. On the other hand the Burkert and Einasto models have a spatial distribution that is much broader (see Fig. 1).

The results for M33 for the DM annihilating into $b\bar{b}$ for the MIN, MED and MAX models in case of SHS and for Einasto and Burkert DM profiles for the SH scenario are reported in Table 8. The significance for the presence of DM for the $b\bar{b}$ annihilation channel, considering the MED DM model, is 18 for the SHS case and 31 for the Einasto profile in the SH scenario. This DM candidate has a mass of 30 (20) GeV and a cross section of $9\cdot 10^{-25}$ cm2/s ($2.3\cdot 10^{-24}$ cm2/s) for SHS (SH) case. The MIN and MAX DM distributions provide similar significances and DM masses but larger and lower values for the best-fit cross sections. The decay scenario with $b\bar{b}$ channel provides $TS$ values of the order of 12 for Einasto and 16 for the Burkert profile with $m_{\chi}=15$ GeV and $\tau=1.3\cdot 10^{25}$. On the other hand the $\tau^{+}\tau^{-}$ channel gives $TS$ values of the order of 11 for Einasto and 14 for the Burkert profile with $m_{\chi}\leq 5$ GeV and $\tau\leq 2.1\cdot 10^{25}$.

The DM candidates for either $b\bar{b}$ and $\tau^{+}\tau^{-}$ annihilation channels are in strong tension with limits found from dwarf spheroidal galaxies Albert et al. (2017) (see Fig. 12). This tension is present for the MIN and MED DM models while it is alleviated for the MAX DM distribution.

IV.5 DM plus astrophysical emission interpretation

In this section we make the more realistic assumption that in addition to a putative DM signal there is also an astrophysical contribution from the galaxy itself. We use the astrophysical emission from M31 and M33 as given in the model found with the baseline fit, thus using an extended source with a disk template with size 0\mbox{~{}.!!^{\circ}}33 for M31 and a point source for M33. Moreover, we include in the analysis the correlations between the SED parameters of DM with the ones of M31 and M33 and the other background sources (see, e.g., Sec. III).

We do not find any significant emission when we include the astrophysical contribution of M31 and M33. Indeed, the $TS$ for the presence of DM is very close to 0 for both annihilation and decay and $b\bar{b}$ and $\tau^{+}\tau^{-}$ channels.

The disk template for M31 is a phenomenological model tuned directly on $\gamma$-ray data and can hide a possible DM contribution. Therefore, we calculate the DM $TS$ also with the following templates: the Herschel/PACS map at 160 $\mu$m, the Spitzer/IRAC map at 3.6$\mu$m, and the atomic gas column density $N_{H}$ map from Braun et al. (2009).

We find no evidence for a DM contribution since the $TS$ is at most of a few considering all the annihilation or decay channels.

We therefore set upper limits on the annihilation cross section $\langle\sigma v\rangle$ or lower limits for the decay time $\tau$.

In Fig. 13 and 14 we show the upper limits for the annihilation cross section $\langle\sigma v\rangle$ and in Fig. 15 we show the lower limit for the decay time $\tau$. Together with the observed limits we also report the expected limits in case there is no signal (see Sec IV.3). First of all we note that the observed limits for the DM annihilation scenarios are in all cases included in the $95\%$ containment bands. For M33 the observed limits are systematically larger than the median expected limits from the null simulations but since for all DM scenarios they are included in the $95\%$ containment band this difference is not significant. Second, the limits derived with the MAX DM distribution model are the strongest for both M31 and M33. This is expected because the $J$ factor for this model is higher than the MED and MIN models. In case of the MED DM model the limits found for M31 constrain the thermal cross section up to about 50 GeV while in case of M33 only the MAX DM model is able to reach the thermal cross section. In the case of decay of DM particles all the limits are well included in the $95\%$ containment bands except for M31 and the case with $b\bar{b}$ channel. This is probably due to a local fluctuation present in the observed limits that is not seen in the simulations but the discrepancy is not significant.

These are the limits found using our benchmark case for the data analysis (energy range and ROI width), astrophysical template for M31 and M33 (disk template for M31 and point source for M33), IEM and isotropic templates. The limits slightly change assuming a different choice for the above cited parameters and models and we reported in Sec. IV.2 the magnitude of these changes.

V Discussion and Conclusions

We performed a systematic study of the $\gamma$-ray emission from M31 and M33 galaxies with a particular focus on a possible DM contribution. We first used our analysis to find the best geometrical model that explains the $\gamma$-ray flux from these sources. For M31 the best model is a uniform disk with a size 0\mbox{~{}.!!^{\circ}}33 and $TS_{\rm{EXT}}=13$ while for M33 is a point-like source.

We also fit $\gamma$-ray emission with templates derived from other wavelengths: (far-)infrared and gas column density maps. These templates provide worse fits for M31 and M33, with respect to data-driven templates, leaving more residuals in the model. The templates that best explain the data are infrared maps that trace the emission from the stellar bulge, meaning that most of the $\gamma$-ray flux that we observe comes probably from this component.

If we interpret the flux from M31 and M33 using only DM we have DM candidates with a mass around 20-50 GeV and a cross section that is close to the thermal one for M31 and is around $10^{-25}$ cm2/s for M33 and the MED DM model. All the DM candidates found with the MED and MIN DM models and for $b\bar{b}$ and $\tau^{+}\tau^{-}$ annihilation and decay channels are ruled out by the current limits found from the MW dwarf spheroidal galaxies Albert et al. (2017).

Finally we made the more realistic assumption that the flux from M31 and M33 comes at least partially from the galaxy. We use in this case the disk template for M31 and the point source morphology for M33 or the templates from infrared or Hydrogen gas column densities for M31. We do not find any excess for the presence of DM for all these cases so we put limits on the annihilation cross section or the decay time that for the MAX and MED DM models constrain the thermal cross section up to 200 GeV and 70 GeV (50 GeV and 10 GeV) for M31 (M33), respectively.

In Figure 16, we compare our results with different limits set by other studies for the $b\bar{b}$ channel. We see that the MED DM model is able to constrain the DM interpretation of the GC excess and that our limits for M31 are similar to the ones derived with Dwarf Spheroidal Galaxies up to about 1 TeV.

Acknowledgements.

The Fermi LAT Collaboration acknowledges generous ongoing support from a number of agencies and institutes that have supported both the development and the operation of the LAT as well as scientific data analysis. These include the National Aeronautics and Space Administration and the Department of Energy in the United States, the Commissariat à l’Energie Atomique and the Centre National de la Recherche Scientifique / Institut National de Physique Nucléaire et de Physique des Particules in France, the Agenzia Spaziale Italiana and the Istituto Nazionale di Fisica Nucleare in Italy, the Ministry of Education, Culture, Sports, Science and Technology (MEXT), High Energy Accelerator Research Organization (KEK) and Japan Aerospace Exploration Agency (JAXA) in Japan, and the K. A. Wallenberg Foundation, the Swedish Research Council and the Swedish National Space Board in Sweden. Additional support for science analysis during the operations phase is gratefully acknowledged from the Istituto Nazionale di Astrofisica in Italy and the Centre National d’Études Spatiales in France. This work performed in part under DOE Contract DE-AC02-76SF00515. X.H. acknowledges supports by the National Natural Science Foundation of China (NSFC-11503078, NSFC-11661161010 and NSFC-11673060) and the “Light of West China” Program of the Chinese Academy of Science. Resources supporting this work were provided by the NASA High-End Computing (HEC) Program through the NASA Advanced Supercomputing (NAS) Division at Ames Research Center.