Gamma-ray image reconstruction of the Andromeda galaxy

We analyze about 12 years of Fermi-LAT data in the direction of the Andromeda galaxy (M31). We robustly characterize its spectral and morphological properties against systematic uncertainties related to the modeling of the Galactic diffuse emission. We perform this work by adapting and exploiting the potential of the SkyFACT adaptive template fitting algorithm. We reconstruct the γ-ray image of M31 in a template-independent way, and we show that flat spatial models are preferred by data, indicating an extension of the γ-ray emission of about 0.3 — 0.4° for the bulge of M31. This study also suggests that a second component, extending to at least 1°, contributes to the observed total emission. We quantify systematic uncertainties related to mis-modeling of Galactic foreground emission at the level of 2.9%.


Introduction
Over the past decade, Fermi-LAT has detected several star-forming galaxies among which we find the Large Magellanic Cloud, the Small Magellanic Cloud, Andromeda galaxy, M82, NGC 253, NGC 2146, and Arp 220. The γ-ray emission coming from these distant galaxies is due to the interactions between the large-scale population of cosmic rays and the interstellar medium and is also produced by high-energy sources such as supernova remnants and pulsars therein [1]. Studying external Milky Way-like galaxies can help provide an independent and outside perspective of the γ-ray astrophysical processes which also occur in our own Galaxy [2]. Therefore, it would bring out additional and complementary knowledge to our understanding of the Milky Way, especially if we can discriminate individual contributions to the γ-ray emission of such extragalactic systems. In this work, we focus on the Andromeda spiral galaxy (M31), our nearest galaxy neighbor, which spans 3.2 • × 1 • on the sky [2,3] and is located at a distance of about 785 ± 25 kpc [4] at Galactic coordinates (l = 121.285 • , b = −21.604 • ) [5]. Thanks to its close proximity, we are able to optically resolve its stellar disk and bulge as two distinct components. This distinction is not possible in our Milky Way Galaxy as the bright emission of the disk obscures that of the bulge [1]. As for γ-ray studies, M31 has the advantage of lying at high Galactic latitudes, away from the Galactic plane which makes it less polluted by the diffuse γ-ray foreground emission of the Milky Way [6]. Several γ-ray analyses have been performed in order to characterize the spectral and morphological shape of M31 where they only find an evidence for a point-like source detection at ∼ 10σ and a 4σ preference for extended emission at best [7]. In addition, a weak evidence for the presence of "Fermi bubble"like structures perpendicularly to M31's galactic plane was found by [8] while [9,10] showed the existence of an extended excess up to 120−200 kpc (8−15 • ) away from the center of M31. In spite of the growing evidence for M31 γ-ray spatial extension, the evaluation of the systematic uncertainties related to the contamination of the Galactic diffuse emission remains the main limitation of the current analyses [9,2]. Thus, the goal of this work is to use the SkyFACT adaptive template fitting algorithm [11] to robustly detect M31 as an extended source against foreground model systematic uncertainties, to characterize its spectral and spatial distributions, and to test the presence of multiple γ-ray emission components. We then aim to reconstruct M31's morphology in a fully 2. Data selection and fitting algorithm We use 624 weeks (12 years) of Fermi-LAT data from August 4th, 2008 to July 16th, 2020, collected from a region of interest (ROI) of a 10 • × 10 • centered at Galactic coordinates (l = 121 • , b = −21 • ). We select Pass 8 SOURCE class events with Point Spread Function (PSF) PSF0+PSF1+PSF2+PSF3 type and use the associated instrument response functions (IRFs) P8R3 SOURCE V2. We apply a cut on maximum apparent zenith angle (z max = 105 • ), and recommended data-quality filters (DATA_QUAL>0) & (LAT_CONFIG==1). We analyze data in the energy range 300 MeV to 100 GeV, logarithmically spaced into 22 bins. We use the Fermi ScienceTools v11r5p3 software package * to prepare the data that we spatially split into 200 × 200 angular bins of size 0.05 • in cartesian sky projection. To perform the γ-ray analysis of M31 we use SkyFACT, a code which combines a hybrid approach template fitting and image reconstruction for studying and decomposing the γ-ray emission. SkyFACT offers the advantage of accounting for expected spatial and spectral uncertainties in each model emission component of the fit, contrary to standard template fitting methods where the spatial distributions of model components are fixed. More details on the analysis code can be found in [11]. The detection of the source and the evidence for its extension are probed through a likelihood ratio (LR) statistical test, denoted TS, between our various models including M31 and our model without the source. The TS results are then interpreted in terms of standard deviations σ using the Chernoff theorem, describing a mixture between a χ 2 and a Dirac function. A detailed description of the statistical analysis framework can be found in [12,13,14,15].

Modeling the γ-ray emission
In our model, the γ-ray emission in the ROI is made up by the following contributions: the Galactic Interstellar Emission (IEM), the Isotropic Gamma Ray Background (IGRB), 13 Pointlike Sources (PS) as reported in the 4FGL [5], and M31. No extended source is present in the chosen ROI. We make use of the Galactic IEM, gll_iem_v06.fits, which includes the Inverse Compton, the π 0 decay, and the Bremstrahlung γ-ray production mechanisms, the spectral shape P8R3_SOURCE_V2.txt for the IGRB †, and the 4FGL catalog, gll_psc_v22.fits, to model the PS [5]. As for M31 galaxy, we model its γ-ray emission using various templates motivated either by previous works [3,2,6] or by studies of the stellar distribution of the galaxy [16], with the final goal of understanding which M31 morphology is preferred by data. We create three kinds of spatial models: (i) 2D gaussian profile whose width σ defines the 68% containment extension angle, (ii) uniform disks, and (iii) Einasto-based models for the nucleus, the bulge, and disk where we integrate the Einasto profile along the line of sight of each pixel and whose profile parameters can be found in [12,16]. The spectral template is given by the 4FGL catalog gll_psc_v22.fits [5].

Standard Template Fitting
Although SkyFACT was developed to simultaneously perform spatial and spectral template fits, it can also be used for traditional or standard template fitting (STF). This analysis setup will allow the comparison with the results of previous literature, as well as with the other SkyFACT runs we will explore in the next sections. Thus, we set the spatial parameters to be fixed to the input models for all components while the spectral parameters are left completely free to vary. We first perform an STF to probe the evidence for M31 in Fermi-LAT data where we compare the fit of a model that only contains IEM, IGRB, and the 13 PS, i.e. our STF baseline model, to the fit of a second model including an additional component modeling M31. Using our PS-limit Gaussian spatial template of 0.001 • width, we compute the significance of an additional source in a nested model with bin-by-bin free spectral parameters. M31 is significantly detected in the 300 MeV-100 GeV energy range with a significance of 7.6σ. Figure 1 (left) shows the flux of each component of the model fitted to the data. We found a reconstructed spectral flux for M31 consistent with the one reported in the 4FGL catalog and parametrized by a power law with exponential cut-off.  3 We note that the lack of statistics generates very large uncertainties on the reconstruction of the M31 flux in the highest energy bins. We then perform a scan over the 68% extension angle θ of the Gaussian spatial profile, from 0.001 • to 1 • , that aims to demonstrate the preference for an extended emission morphology and measure the extension of the signal. We obtain a log-likelihood (−2 ln L) profile as a function of θ, as shown in Fig. 1 (right). The log-likelihood profile is minimized by non-zero extension angles, indicating a preference for an extended γ-ray emission for the M31 model component. The highest detection significance of the source is σ = 9.2, and it is associated to the Gaussian template of width θ = 0.5 • . We fit the log-likelihood profile with a cubic spline (green line in Fig. 1 (right)) and find the minimum of the curve at θ = 0.48 • ± 0.07 which corresponds to the best-fit extension angle. In addition, this determination of the extension is very robust against systematic uncertainties related to mis-modeling of the PSF. We quantify the systematic uncertainties on θ at the level of 1.5% due to PSF mis-modeling. The details of the computation are available in [12]. Our result for extension is compatible with the value of [2] who found a size of 0.42 • ± 0.10 using a Gaussian template.

Gaussian profiles
We then use SkyFACT to test systematic uncertainties of the extension related to the mis-modeling of the foreground emission template. We allow more freedom on the spatial part of the IEM diffuse component so that we can test (i) whether or not the extension of the source is still preferred against uncertainties on the IEM template, and (ii) the robustness of the reconstructed spectra of M31 and IEM components against the freedom given to the baseline parameters. We first define three configurations of baseline models allowing a 30%, 50%, and 100% variation of the IEM spatial parameters. For each of these baseline models, we perform a scan over θ using a Gaussian profile following the same procedure as in the STF configuration. We found compatible values for the best-fit angular extension, θ = 0.48 • ± 0.07 (IEM 30%), θ = 0.48 •+0.07 −0.08 (IEM 50%), and θ = 0.45 •+0.09 −0.10 (IEM 100%), from which we evaluate the systematic uncertainty on θ due to the IEM mis-modeling to be 2.9%. In addition, the reconstruction of the spectra for both the IEM and M31 model components is robust against the variation allowed on the IEM spatial parameters, and is compatible with the STF results.

Uniform profiles
We then test a uniform disk template of various sizes from 0.025 • to 1 • . We interestingly find two minima of the log-likelihood profile: the first one at an angular extension of 0.4 • ±0.04 associated to M31 bulge and compatible with the value of 0.33 • ± 0.04 found by [2] for a uniform disk template. We find a second minimum at 0.95 +0.05 −0.18 • indicating the presence of a more extended structure around M31, such as galactic bubbles [8]. Both extensions have been detected with significance of 8.5σ. We tested the presence of bubble-like structures lying above and below the M31 galactic plane on top of a centrally concentrated emission. We only find marginal evidence (< 3σ) for this additional emission component.

Einasto profiles
We finally test Einasto spatial templates as motivated by recent stellar mass models of M31. Here, we do not perform any scan over the Einasto template parameters, since we want to test specific models tracing M31 stellar components including a nucleus, a bulge, and a disk. For the nucleus and the bulge, we keep the same spectral template of M31 extracted from the 4FGL. As for the disk spectral model, we use a power law of index -2.4 [8]. We separately test the evidence for each of these three components on top of our baseline model (IEM 30%) and also run fits for all possible combinations of these Einasto spatial templates. We observe that the data prefers to complement centrally concentrated components (nucleus-or bulge-only models) with an extended disk-like component (∼ 3σ evidence), which has an angular extension of about 1.08 • (semi-major axis). This result is also compatible with the second minimum 0.95 • of the uniform disk scan. This corresponds to a further indication that the M31 γ-ray emission is indeed more extended than what was previously found in most of the literature showing an evidence for the bulge only component [2,7,3]. On the other hand, we cannot claim evidence for a centrally concentrated component (nucleus or bulge) on top of a disk-only model, as we found a marginal improvement of 0.8σ only.

Model comparison
We properly compare spatial template models through the evaluation of the AIC score of each model. The AIC score takes into account the number of effective parameters in order to assess the goodness of the models from one to another. More details can be found in [17]. Based on ∆AIC values, we find the data rather prefers flat templates to describe the γ-ray emission of M31. As a matter of fact, the uniform disk templates come first in our ranking, followed by the Einasto nucleus+disk model, then the Einasto disk and bulge+disk models, and finally the Gaussian templates. This highlights once again the preference for a more extended emission, which goes beyond M31's bulge (See [12] for more details on these sATF fits and their comparison).

Image reconstruction of M31 morphology
We finally take full advantage of SkyFACT's flexibility by performing an adaptive template fit allowing a 100% freedom on the spatial and spectral parameters of the M31 component. With this analysis mode, we do not assume any spatial input morphology, but we rather define a large circular region of 1.5 • radius around the source position, within which the M31 component is adjusted to the data in each pixel and each energy bin. We keep the same spectral shape from the 4FGL for the M31 component. The reconstructed morphology of M31 is presented in Fig. 2 (left).   We build the intensity profile of M31 as shown in Fig. 2 (right) by summing the contributions of all pixels in concentric annuli from the reconstructed M31 image. The annuli are spaced of 0.1 • from 0 • to 1.5 • and are centered around M31 4FGL position. The intensity of each annulus is averaged over the solid angle ∆Ω annulus = N ∆Ω pixel , where N is the number of pixels in a given annulus.
The intensity profile appears rather flat and does not follow a simple known function. We fit the 5 intensity using a combination of a Gaussian function and a power law given by: where the best-fit coefficients are found to be a = 1.48 × 10 −6 , µ = −1.04 × 10 −1 , σ = 4.35, b = 4.72 × 10 −8 , and c = 2.37. This function Ξ(x) suggests a superposition of two distinct components of M31 emission where the Gaussian function would model M31 central emission as it is peaked at the center of the source and decreases at larger radii, while the power law could correspond to the extended structures around the galaxy, such as galactic bubbles, the outer part of the disk, and/or a larger cosmic-ray halo. We find an inflection point of the intensity profile at 0.3 • by taking its second derivative. This corresponds to the angle at which the power law starts kicking in and may delimit the apparent size of the bulge in γ rays. This value is compatible with the result of [2] who found an extension angle of 0.33 • ± 0.04 using a uniform disk template to model M31's bulge, and also marginally consistent with our first minimum of the uniform disk scan 0.40 • ± 0.04.