Gamma-ray detection toward the Coma cluster with Fermi-LAT: Implications for the cosmic ray content in the hadronic scenario
AAstronomy & Astrophysics manuscript no. coma˙paper © ESO 2021February 5, 2021
Gamma-ray detection toward the Coma cluster with
Fermi -LAT:Implications for the cosmic ray content in the hadronic scenario
R. Adam (cid:63) , H. Goksu , S. Brown , L. Rudnick , and C. Ferrari Laboratoire Leprince-Ringuet (LLR), CNRS, ´Ecole polytechnique, Institut Polytechnique de Paris, 91120 Palaiseau, France Department of Physics & Astronomy, The University of Iowa, Iowa City, IA, 52245 Minnesota Institute for Astrophysics, University of Minnesota, Minneapolis, MN 55455 Universit´e Cˆote d’Azur, Observatoire de la Cˆote d’Azur, CNRS, Laboratoire Lagrange, FranceReceived February 5, 2021 / Accepted –
Abstract
The presence of relativistic electrons within the di ff use gas phase of galaxy clusters is now well established, thanks to deep radioobservations obtained over the last decade, but their detailed origin remains unclear. Cosmic ray protons are also expected to accu-mulate during the formation of clusters. They may explain part of the radio signal and would lead to γ -ray emission through hadronicinteractions within the thermal gas. Recently, the detection of γ -ray emission has been reported toward the Coma cluster with Fermi -LAT. Assuming that this γ -ray emission arises essentially from pion decay produced in proton-proton collisions within the intraclustermedium (ICM), we aim at exploring the implication of this signal on the cosmic ray content of the Coma cluster and comparing itto observations at other wavelengths. We use the MINOT software to build a physical model of the Coma cluster, which includes thethermal target gas, the magnetic field strength, and the cosmic rays, to compute the corresponding expected γ -ray signal. We applythis model to the Fermi -LAT data using a binned likelihood approach, together with constraints from X-ray and Sunyaev-Zel’dovichobservations. We also consider contamination from compact sources and the impact of various systematic e ff ects on the results. Weconfirm that a significant γ -ray signal is observed within the characteristic radius θ of the Coma cluster, with a test statistic TS (cid:39) + ff use emission of the cluster itself as it is strongly degenerate with theexpected ICM emission, and extended models match the data better. Given the Fermi -LAT angular resolution and the faintness of thesignal, it is not possible to strongly constrain the shape of the cosmic ray proton spatial distribution when assuming an ICM originof the signal, but preference is found in a relatively flat distribution elongated toward the southwest, which, based on data at otherwavelengths, matches the spatial distribution of the other cluster components well. Assuming that the whole γ -ray signal is associatedwith hadronic interactions in the ICM, we constrain the cosmic ray to thermal energy ratio within R to X CRp = . + . − . % and theslope of the energy spectrum of cosmic rays to α = . + . − . ( X CRp = . + . − . % and α = . + . − . when including both the cluster and4FGL J1256.9 + + γ -ray emission is detected in the direction of the Coma cluster. Assuming that the emission is due to hadronic interactions in theintracluster gas, they provide the first quantitative measurement of the cosmic ray proton content in a galaxy cluster and its implicationfor the cosmic ray electron populations. Key words.
Galaxies: clusters: individual: Coma – Cosmic rays – γ -rays: galaxies: clusters
1. Introduction
Galaxy clusters form hierarchically via the smooth accretion ofsurrounding material and the merging of subclusters and groups(Kravtsov & Borgani 2012). During these processes, most ofthe binding gravitational energy is dissipated into the hot, ther-mal, ionized gas phase. However, shock waves propagating inthe intracluster medium (ICM) and turbulence are also expectedto accelerate both electrons and protons at relativistic energies,leading to a nonthermal population of cosmic rays (CR) that areconfined within the cluster magnetic fields (see Brunetti & Jones2014, for a review). In addition, CR may also arise from di-rect injection in the ICM, trough active galactic nucleus (AGN)outbursts (e.g., Bonafede et al. 2014a), or galactic winds asso- (cid:63)
Corresponding author: R´emi Adam, [email protected] ciated with star formation activity in cluster member galaxies(e.g., Rephaeli & Sadeh 2016). While high energy cosmic rayelectrons (CRe) should quickly lose their energy ( ∼ yearsat 10 GeV; see e.g., Sarazin 1999), cosmic ray protons (CRp)should accumulate over the cluster formation history (see, e.g.,Pfrommer et al. 2007, for the simulation of CR in clusters).Nevertheless, the CR physics and content of galaxy clusters re-main largely unknown to date.In fact, the presence of CRe and magnetic fields in galaxyclusters is now well established, thanks to the deep observa-tions of di ff use radio synchrotron emission obtained over thelast decade (see van Weeren et al. 2019, for a review). We candistinguish between radio relics, believed to trace merger shockacceleration at the periphery of clusters (e.g., van Weeren et al.2010), and radio halos, which are spatially coincident with thethermal gas in massive clusters. Radio halos are generally fur- a r X i v : . [ a s t r o - ph . H E ] F e b . Adam et al.: γ -ray emission toward Coma ther divided into mini-halos, associated with relaxed cool-coreclusters (Giacintucci et al. 2017), and giant radio halos, whichare megaparsec-size sources found in merging clusters (Cassanoet al. 2010), although the distinction is not necessarily clear(Ferrari et al. 2011; Bonafede et al. 2014b; Savini et al. 2019).Two main mechanisms have been proposed in the literature toexplain the origin of CRe that generate radio halos: the hadronicmodel (Dennison 1980; Blasi & Colafrancesco 1999; Dolag &Enßlin 2000), in which secondary CRe are the products of piondecay generated in collisions between the CRp and the ther-mal gas, and the turbulent reacceleration of a seed population ofelectrons (possibly secondary CRe; Brunetti & Lazarian 2007,2011). In both scenarios, CR might show up as a γ -ray signaldue to the inverse Compton interaction of CRe with backgroundlight. In the hadronic model, the decays of pions should also leadto additional γ -ray (and neutrino) emission, which is expected tobe the dominant component of the γ -ray flux at energies above ∼
100 MeV (Pinzke & Pfrommer 2010).Attempts to detect the γ -ray emission from galaxy clustershave been carried out over the last two decades using individ-ual targets, stacking analysis, and cross correlations to externaldatasets (see, e.g., Reimer et al. 2003; Ackermann et al. 2010;Aleksi´c et al. 2010; Dutson et al. 2013; Huber et al. 2013; Gri ffi net al. 2014; Branchini et al. 2017; Colavincenzo et al. 2020).Despite the lack of detection, these searches proved extremelyuseful in constraining the CRp content of clusters to below afew percent relative to the thermal pressure (Ackermann et al.2014).Among the relevant individual targets, the Coma cluster hasbeen one of the most promising sources to search for γ -rayemission. Indeed, it is a massive ( M (cid:39) × M (cid:12) ) andnearby system (redshift z = . γ -ray flux is expected to be large (Pinzke & Pfrommer 2010).Additionally, it is located near the galactic north pole and thusbenefits from a low galactic background. The signal is expectedto be extended even for the Fermi -LAT (characteristic radius θ ∼ .
75 deg), but not so extended that the analysis su ff erstrong systematic e ff ects in the background modeling. The re-gion around Coma is not a ff ected by the presence of bright γ -ray compact sources. Finally, the Coma cluster is a well-knownongoing merger, which has been extensively observed at vari-ous wavelengths using various probes (e.g., Kent & Gunn 1982;Briel et al. 1992; Gavazzi et al. 2009; Bonafede et al. 2010;Planck Collaboration et al. 2013; Mirakhor & Walker 2020).In particular, a well-measured giant radio halo and a radio relichave been observed, proving the presence of CRe (e.g., Willson1970; Giovannini et al. 1993; Thierbach et al. 2003; Brown &Rudnick 2011; Bonafede et al. 2020).Several analyses have focussed on the search for γ -rays to-ward the Coma cluster, using both ground-based and space-based instruments (e.g., Perkins 2008; Aharonian et al. 2009;Arlen et al. 2012; Prokhorov 2014; Zandanel & Ando 2014). Itis noteworthy that the Coma cluster was analyzed by the Fermi -LAT collaboration (Ackermann et al. 2016), with six years ofdata, who found some residual emission in the direction of thecluster, though not enough to claim a detection. Despite theseunsuccessful searches, the combination of the obtained γ -ray up-per limits with radio synchrotron data was used to show that purehadronic models cannot explain all the observed radio emissionin the case of Coma (Brunetti et al. 2012), and this allowed con-straints to be set on turbulent reacceleration models (Brunettiet al. 2017). Recently, Xi et al. (2018) claimed the first signif-icant detection of γ -ray signal toward the Coma cluster, using Fermi -LAT data. This detection was also discussed in the con- text of a dark matter interpretation of the signal in Liang et al.(2018), but only minor investigations of the consequences forthe CR physics were carried out. In 2020, the
Fermi -LAT col-laboration released an update of the 4FGL catalog (4FGL-DR2,Abdollahi et al. 2020; Ballet et al. 2020), also indicating that asource was detected in the direction of the Coma cluster, named4FGL J1256.9 + ff use emission it-self.In this paper, we explore the consequences on the CR physicsof the γ -ray emission observed in the direction of the Coma clus-ter under the assumption that this signal is associated with thedi ff use ICM. To do so, we present an analysis of the Fermi -LATdata (Atwood et al. 2009), to search for γ -ray emission within theComa cluster region, and use this measurement to constrain theproperties of the CR content of the cluster. We construct a modelfor the expected signal, in which we set the thermal gas prop-erties using thermal Sunyaev-Zel’dovich (tSZ) e ff ect and X-rayobservations. The Fermi -LAT data are analyzed, and we extractthe spectral energy distribution (SED) of the emission assumingdi ff erent spatial templates. The Fermi -LAT maps are comparedto data at other wavelengths for better interpretation of the re-sults. In particular, we compare the γ -ray signal to the radio syn-chrotron emission, to the thermal pressure traced by the tSZ sig-nal, to the galaxy distribution, and to the thermal gas distributionfrom X-ray data. We use our model to constrain the CRp popu-lation in the Coma cluster assuming that all the signal observedis due to the hadronic interactions between CRp and the ambientthermal gas, but also consider possible point source contamina-tion. Finally, we compute the amount of secondary CRe associ-ated with the γ -ray emission assuming steady state. We use thismodel to constrain the fraction of primary to secondary CRe inthe absence of reacceleration, and the boost required from reac-celeration with respect to the steady state hadronic model, toexplain the radio emission assuming that reacceleration of sec-ondary electrons is dominant, respectively.This paper is organized as follows. Section 2 describesthe multiwavelength dataset that is used through the paper. InSection 3, we present the physical modeling of the cluster andthe computation of the observables. The Fermi -LAT analysis ispresented in Section 4 and the outputs are compared to multi-wavelength data in Section 5. Section 6 discusses the implica-tions of the observed signal for the CR content of the Comacluster. In Section 7, we discuss the results presented in thispaper and compare them to the literature. Finally, a summaryand conclusions are given in Section 8. Throughout this paper,we assume a flat Λ CDM cosmology according to
Planck re-sults (Planck Collaboration et al. 2016b) with H = . − Mpc − , Ω M = . Ω Λ = . = z = . Planck cat-alog (PSZ2, Planck Collaboration et al. 2016a). Given the ref-erence cosmological model, 1 degree corresponds to 1.73 Mpcat the redshift of Coma. We use the value of R = M = . × M (cid:12) given our cosmological model.
2. Data
While this paper focuses mainly on the γ -ray analysis of theComa cluster, we also use complementary data both for mul-tiwavelength comparison of the signal and in order to buildtemplates used to model the γ -ray emission. In addition to the Fermi -LAT γ -ray data, we thus use thermal Sunyaev-Zel’dovich
2. Adam et al.: γ -ray emission toward Coma data from Planck , X-ray data from
ROSAT , radio data from theWesterbork Synthesis Radio Telescope and optical data from theSloan Digital Sky Survey. This section presents a brief descrip-tion of these data. γ -ray data The γ -ray analysis is based on the publicly available Fermi
Large Area Telescope data (
Fermi -LAT, Atwood et al. 2009),which has been collecting all sky γ -ray data at energies fromabout 20 MeV to more than 300 GeV since June 2009. In thispaper, we use almost 12 years of Pass 8 data (P8R3), collectedfrom August 4, 2008, to April 2, 2020. Events within a radius of10 degrees from the cluster center were collected for the analysis(see also Section 4 for the data selection). As a description of thecompact sources around the cluster, we use the second release ofthe 4FGL catalog (see Abdollahi et al. 2020 for the 4FGL cata-log and Ballet et al. 2020 for the 4FGL-DR2 update). The 4FGL-DR2 catalog is an incremental version of the 4FGL catalog. It isbased on 10 years of survey data in the 50 MeV - 1 TeV range(see Abdollahi et al. 2020, for the analysis improvement relativeto previous catalogs). The di ff use isotropic background and thegalactic interstellar emission are modeled using the latest avail-able templates: iso P8R3 SOURCE V2 v1 and gll iem v07 , re-spectively. The tSZ e ff ect (Sunyaev & Zeldovich 1972) is a relevant probefor comparison to the γ -ray signal as it provides a direct mea-surement of the integrated thermal gas pressure along the line-of-sight (see Birkinshaw 1999; Mroczkowski et al. 2019, for re-views). We used the MILCA (Hurier et al. 2013) Compton pa-rameter map obtained from
Planck (Planck Collaboration et al.2016c) to extract a high signal-to-noise image of the signal atan angular resolution of 10 arcmin (full width at half maximum,FWHM). In the region of Coma, the signal is not significantly af-fected by any artifacts, e.g., from compact sources, or remainingdi ff use galactic emission (see also Planck Collaboration et al.2013, for a Planck analysis of the Coma cluster). The data arepublicly available at the
Planck
Legacy Archive . The X-ray di ff use cluster emission traces the thermal gas density(see B¨ohringer & Werner 2010, for a review) and is expected tobe correlated with the γ -ray signal. Given the large size of theregion of interest (ROI) and the fact that the angular resolutionis not critical for our purpose, we used ROSAT (Truemper 1993)data as obtained from the
ROSAT
All Sky Survey (RASS) toimage the cluster at a resolution of 1.8 arcmin (FWHM). To doso, we subtracted the background and normalized the data us-ing the exposure maps. In addition, we subtracted the compactsources present in the field using both the RASS bright and faintsource catalogs (Voges et al. 1999, 2000). https://fermi.gsfc.nasa.gov/ssc/data/analysis/ http://pla.esac.esa.int/pla/ Collected from http://cade.irap.omp.eu/dokuwiki/doku.php?id=rass
The radio emission traces relativistic electrons via synchrotronemission and could be associated with the γ -ray signal. We usedthe Westerbork Synthesis Radio Telescope (WSRT) radio dataat 352 MHz (Brown & Rudnick 2011) in order to dispose of animage of the cluster halo and relic. Compact sources were sub-tracted using the method described in Rudnick (2002). The resid-ual contamination was estimated by injecting fake point sourcesinto the image, applying the filter to the image, and measuringthe residual flux. For sources that were detected by more than2 σ , the residual flux was less than 3% and as such the contami-nation is expected to be less than 3%. Nevertheless, some of theemission from the bright central radio galaxy remains blendedwith the signal, and we masked it in the following analysis us-ing a threshold of 100 mJy / beam. Additionally, the bright sourceComa A north of the radio relic causes some artifacts that shouldbe accounted for in the morphological comparison (Section 5).The image angular resolution is 4 × (FWHM).The WSRT map is also used to extract the average radialprofile of the radio halo. To do so, we first subtract the zerolevel (i.e., o ff set), and define bins of 3 arcmin width. The zerolevel arises from the standard clean deconvolution that createsa negative o ff set around bright sources such as those found inthe Coma cluster. After filtering out the point sources and con-volving the image to a lower resolution, this resulted in a lo-cal background of -20 mJy / beam. Error bars on the profile arecomputed using inverse variance Gaussian weighting with a con-stant noise root-mean-square of 5 mJy / beam. The pixels from themasked regions are excluded from the profile. We also convertfrom Jy / beam to Jy / sr assuming that the beam is Gaussian. Westress that we are considering the total (i.e., average) profile, butthe variations in profile from one direction to another (Brown &Rudnick 2011) should not have a significant e ff ect on our anal-ysis given the purpose of this paper. The profile will be used inSection 6 as a comparison to our model.In addition, we also use the integrated fluxes of the Coma ra-dio halo as compiled in Pizzo (2010). Given the di ff use nature ofthe signal, it is important to make sure that the flux we used wasmeasured in a consistent aperture for the di ff erent observations.In this respect, we used a radius of 0 . × R (i.e., 629 kpc or0.36 deg, see Brunetti et al. 2013, for details). The spatial distribution of galaxies is informative of the di ff erentgroups in the cluster, and the γ -ray emission might be associatedwith CR escaping from cluster member galaxies. We thus useSloan Digital Sky Survey (SDSS, York et al. 2000) optical data to construct a color image of the Coma region. In addition, weselect galaxies with spectroscopic information from the SDSSdatabase and use this catalog to construct a galaxy density mapof the Coma region. To to so, we bin the galaxy catalog on a grid,and we weight them with w = exp (cid:18) − ( z gal c − z Coma c ) )2 σ v (cid:19) to minimizethe background. The quantity c is the speed of light, z gal and z Coma the redshift of each galaxy and that of the cluster, respec-tively, and σ v = − (i.e., FWHM of 5000 km s − ). Themap is convolved with a Gaussian kernel with a FWHM of 10arcmin to enhance the signal-to-noise. The background is thenestimated as the mean of the map at radii above 90 arcmin fromthe cluster, and subtracted from the map. https://dr12.sdss.org/mosaics http://skyserver.sdss.org/CasJobs/
3. Adam et al.: γ -ray emission toward Coma
3. Modeling
The physical modeling of the signal, as expected in the γ -rayband at Fermi -LAT energies, is necessary for two reasons. First,we aim at constructing a template (spatial and spectral) of thedi ff use emission in order to include the cluster when fitting thedata (as done in Section 4). Then, the physical model is neededto connect the observations to the CR physics and will be usedto constrain the CR content of the Coma cluster (Section 6).Additionally, the synchrotron emission should be included inthe modeling when comparing the implication of the Fermi -LATconstraint on the radio signal.In this paper, we model the Coma cluster using the
MINOT software (Adam et al. 2020) . This allows us to have a completedescription of the cluster from the CR and thermal physics to the γ -ray signal, under the approximation of spherical symmetry.In addition, we also use spatial templates constructed usingmultiwavelength data. While this approach may provide a bettermatch to the spatial shape of the signal due to deviations fromspherical symmetry, it does not allow us to constrain the CR con-tent of the cluster because of sky projection and degeneraciesbetween the di ff erent physical components. MINOT
The
MINOT software models the thermal ICM component usingthe gas electron pressure and density profiles. The gas densityprofile is key to the γ -ray emission as it provides the target forproton-proton interactions. The thermal pressure is also particu-larly important, because it will be used to provide a normaliza-tion to the amount of CR.The electron pressure radial profile, as a function of radius r , is described by a generalized Navarro Frenk White (gNFW,Nagai et al. 2007) model according to P e ( r ) = P (cid:18) rr p (cid:19) c (cid:18) + (cid:18) rr p (cid:19) a (cid:19) b − ca , (1)with parameters taken from Planck Collaboration et al. (2013)and rescaled to our cosmological model ( P , r p , a , b , c = .
022 keV cm − , . , . , . , . β -model (Cavaliere &Fusco-Femiano 1978) according to n e ( r ) = n + (cid:32) rr c (cid:33) − β/ , (2)with parameters taken from Briel et al. (1992) andrescaled to our cosmological model ( n , r c , β = . × − cm − ,
310 kpc , . Y He = . .
3, relative to the so-lar value (see Adam et al. 2020, for more details). Finally, wenote that the radial models are truncated at a radius r > × R .This also apply for nonthermal quantities discussed below. https://github.com/remi-adam/minot The synchrotron emission, from a given CRe population, de-pends on the magnetic field. The later is thus an important com-ponent of our model when comparing the implications of ourconstraints with radio data (Section 6.3). The magnetic fieldstrength profile of the Coma cluster was inferred by Bonafedeet al. (2010) using Faraday rotation measures and we used theirresult to model it as (cid:104) B (cid:105) ( r ) = (cid:104) B (cid:105) (cid:32) n e ( r ) n (cid:33) η B , (3)where (cid:104) B (cid:105) = . µ G and η B = .
5. We note, how-ever, that uncertainties on these parameters are relativelylarge and the parameters are degenerate (see also the re-cent work by Johnson et al. 2020). They are constrained inthe range ( (cid:104) B (cid:105) = . µ G; η B = .
4) to ( (cid:104) B (cid:105) = . µ G; η B = . The CRp number density, per unit volume and energy, is mod-eled as J CRp ( r , E ) = A CRp f ( E ) f ( r ) , (4)where f ( E ) is the energy spectrum, f ( r ) the radial dependence,and A CRp is the normalization. The energy spectrum is modeledas a power law (expected from Fermi acceleration), as f ( E ) ∝ E − α CRp , (5)with α CRp the slope of the CRp energy spectrum. The radialdependence is modeled assuming that the CRp number densityscales with either the thermal density or pressure, as f ( r ) ∝ (cid:40) n e ( r ) η CRp P e ( r ) η CRp . (6)Since the radial dependence of the CRp is not well-known, wetest di ff erent values for the scaling, η CRp = (cid:16) , , (cid:17) , correspond-ing to flat, extended, and compact distributions. Finally, the nor-malization A CRp is computed relative to the thermal energy en-closed within the radius R . In the following, we use the pa-rameter X CRp ( R ) = U CRp ( R ) U th ( R ) , (7)where U CRp ( R ) is the total CRp energy enclosed within R com-puted by integrating Equation 4 and U th ( R ) the total thermal en-ergy, obtained by integrating the gas pressure profile (see Adamet al. 2020, for more details).In the end, three parameters are used to describe the CRp: thespatial scaling relative to the thermal gas η CRp , the CRp spectrumslope α CRp , and the normalization encoded in X CRp ( R ). Thevalue of these parameters will be further discussed in Section 4and 6. At the
Fermi -LAT energies, the γ -ray emission is expected tobe dominated by hadronic processes, and we neglect the in-verse Compton emission due to secondary or primary CRe in-teracting with background light (see Appendix A for more dis-cussions). However, CRe need to be included in our model
4. Adam et al.: γ -ray emission toward Coma Radius (kpc)10 n C R p ( c m ) n CRp n e n CRp P e n CRp n e n CRp = constant 10 Radius (kpc)10 X C R p n CRp n e n CRp P e n CRp n e n CRp = constant10 Radius (kpc)10 S ( c m s s r ) n CRp n e n CRp P e n CRp n e n CRp = constant 10 Energy (GeV)10 E d N d E d V ( G e V c m s ) n CRp n e , CRp = 2.8 n CRp n e , CRp = 2.3 n CRp n e , CRp = 3.3 n CRp n e Figure 1:
MINOT templates used to describe the cluster, for di ff erent assumptions regarding the CRp distributions. All models are normalized sothat X CRp ( R ) = − . Top left : Radial profile of the CRp distribution.
Top right : Enclosed CRp to thermal energy profile.
Bottom left : γ -raysurface brightness profile. Bottom right : γ -ray energy spectrum integrated within R . when computing the radio synchrotron emission in Section 6.The primary cosmic ray electrons, CRe , are modeled simi-larly as the CRp. However, CRe are a ff ected by losses andwe therefore account for this using three di ff erent models: the ExponentialCutoffPowerLaw , the
InitialInjection , andthe
ContinuouslInjection models as implemented in the
MINOT software (Adam et al. 2020). They are expressed as f ( E ) ∝ (cid:32) EE (cid:33) − α CRe1 × exp (cid:32) − EE cut , CRe (cid:33) , (8) f ( E ) ∝ (cid:32) EE (cid:33) − α CRe1 (cid:18) − EE cut , CRe1 (cid:19) α CRe1 − E < E cut , CRe E (cid:62) E cut , CRe , (9)and f ( E ) ∝ (cid:32) EE (cid:33) − ( α CRe1 + ) (cid:18) − EE cut , CRe1 (cid:19) α CRe1 − E < E cut , CRe E (cid:62) E cut , CRe , (10)respectively. Using di ff erent parametrizations will allow us toestimate the systematic e ff ects associated with the model.There are therefore four parameters used to describe the pri-mary CRe : the spatial scaling of the CRe population relativeto the thermal gas, η CRe (as in Equation 6), the CRe spectrumslope α CRe and break energy E cut , CRe , and the normalization encoded in X CRe ( R ) (as in Equation 7). The value of theseparameters will be further discussed when comparing our modelto radio data in Section 6. γ -ray and synchrotron signal Given the thermal gas, magnetic field and CRp modeling, wecompute the γ -ray emission due to hadronic collisions, the sec-ondary cosmic ray electrons assuming a steady state scenario(CRe ), and the radio synchrotron emission (from both CRe and CRe ) according to Adam et al. (2020). We do not accountfor inverse Compton emission at Fermi -LAT energies because itis expected to be negligible, as detailed in Appendix A.In brief, the production rate of γ -rays and CRe are computedas dNdEdVdt ( r , E ) = (cid:90) + ∞ E dN col dE CRp dVdt F (cid:16) E , E CRp (cid:17) dE CRp , (11)where dN col dE CRp dVdt ∝ n e ( r ) × J CRp ( r , E ) is the CRp–ICM collisionrate and the function F (cid:16) E , E CRp (cid:17) gives the number of secondaryparticles (electrons or γ -rays) produced in a collision as a func-tion of the initial energy of the CRp. This computation is basedon the parametrization by Kelner et al. (2006) and Kafexhiuet al. (2014) following the work by Zabalza (2015), as detailedin Adam et al. (2020). The distribution of CRe is then obtained
5. Adam et al.: γ -ray emission toward Coma D e c . ( d e g r ee s ) Galaxies D e c . ( d e g r ee s ) tSZ D e c . ( d e g r ee s ) X-ray D e c . ( d e g r ee s ) Radio halo D e c . ( d e g r ee s ) Radio relic Figure 2:
Template images used to describe the morphology of the γ -ray emission. From left to right and top to bottom, we present the templatesbased on the distribution of galaxies, the tSZ signal, the X-ray emission, the radio halo emission, and the radio relic emission. Images are 3 × , except for the radio relic (bottom right), which is only 1 × . Units are arbitrary, but the integral of all maps over the solid angle is setto unity. For display purpose, the scale is linear from zero to 20% of the peak value, and logarithmic above. accounting for losses under the steady state assumption. The γ -ray emission profile and spectrum are computed by integratingEquation 11 along the line-of-sight, accounting for the distanceto the Coma cluster, and integrating over the energy or the solidangle, respectively. The γ -ray attenuation due to interaction withthe extragalactic background light (EBL) is also included and in-duces a cuto ff in the spectrum above E (cid:38) GeV. These spatial(profile) and spectral templates are used in Section 4 to accountfor the cluster in the modeling of the ROI.The synchrotron emission rate is computed following theprescription given in Aharonian et al. (2010). This assumes thatthe orientation of the magnetic field is randomized. Observableprofile and spectra are then obtained as for the γ -rays.In Figure 1, we show how the di ff erent profiles and spec-tral slopes for the CRp translate into the γ -ray surface bright-ness profile and spectrum. The top left panel shows the fourdi ff erent radial distributions used for the CRp, in terms of n CRp ( r ) = (cid:82) J CRp ( r , E ) dE . As we can see, the di ff erence be-tween models based on the thermal density or thermal pressureare very close in the case of Coma as expected given the fact thatthe cluster is close to be isothermal. On the top right panel, wecan see the correspondence in terms of the profile of the CRp en-ergy to thermal energy ratio. We note that the normalization wasfixed to X CRp ( R ) = − . The bottom left panel shows the γ -ray surface brightness for the di ff erent models. Because it arisesfrom the product of the CRp and thermal gas density, the pro-files are more peaked than the original CRp distributions. On thebottom right panel, we can see the energy spectrum integratedwithin R for di ff erent slopes α CRp . The CRp slope mainlydrives the γ -ray slope between 1 GeV and 10 GeV. At higher en- ergies, the extragalactic background light induces a cuto ff in thespectrum, but this is not in the range accessible by Fermi -LAT.At lower energies, the energy threshold of the proton-proton col-lisions implies that the spectrum smoothly vanishes. Since thenormalization is fixed, steeper is the spectrum and higher is theflux near the peak. As seen with the dash-dotted lines, the am-plitude decreases when the CRp profiles become flatter becauseit leads to a decrease in the particle collision rate within R . In addition to the spherically symmetric physical models, webuild spatial templates based on the multiwavelength data dis-cussed in Section 2. We consider templates based on the galaxydensity, the tSZ Compton parameter, the X-ray surface bright-ness and the radio surface brightness (halo and relic). The mapsare projected on 5 × maps with 1 arcmin pixel size (wellbelow the Fermi -LAT resolution). In order to reduce the noise,the maps are smoothed so that their angular resolution (FWHM)is 10, 11, 5, and 5 arcmin for the galaxy density, tSZ, X-ray,and radio maps, respectively. This smoothing remains well be-low the
Fermi -LAT angular resolution so that it will be negligiblewhen convolving the maps to the instrument response function(Section 4). In order to minimize bias from noise on large scales,we mask the pixels that are more than 90 arcmin away from thereference center for the galaxy density, tSZ and X-ray maps. Forthe radio map, two cases are considered: the halo for which pix-els more than 60 arcmin from the center are masked (to avoidincluding the relic in the template) and the relic for which we
6. Adam et al.: γ -ray emission toward Coma mask pixels more than 25 arcmin from the coordinates (RA,Dec) = (193.8, 27.2) deg. The mask used for the relic allows us to ex-clude any bright radio galaxy focusing on the emission from therelic only.The resulting templates are shown in Figure 2 on 3 × grids (except for the relic, for which it is 1 × ). In all cases,the cluster is elongated from the northeast to the southwest, withan excess associated with the NGC 4839 group on the southwest.The most compact template is the one based on the X-ray imageas it is proportional to the gas density squared. The tSZ template,proportional to the gas pressure (temperature times density) ismuch more extended. The galaxy density template extension isin between the tSZ and X-ray ones, but it is more elongated.The radio halo template, on the other end, is very spherical andmatches well the tSZ template in terms of the extension. Theradio relic template is very elongated from the southeast to thenorthwest, much smaller than the other templates, and not spa-tially coincident (as will be further discussed in Section 5). Fermi -LAT analysis
In this Section, we describe the
Fermi -LAT γ -ray analysis, per-formed using Fermipy (Wood et al. 2017). After the data selec-tion, we model and fit the ROI in order to extract the signal inthe Coma cluster region under di ff erent modeling assumptions(signal associated with the ICM, a point source, the combinationof both) and compare these scenarios. The cluster model built inSection 3 is used to extract the SED of the source. We also per-form several tests to validate the global ROI model as a functionof energy, radius and check the time stability. Finally, we discussthe systematic e ff ects that might a ff ect our findings. Following the
Fermi -LAT collaboration recommendations foro ff galactic plane compact source analysis, we apply the P8R3 SOURCE V2 selection (event class 128). The energy disper-sion is also accounted for. We select
FRONT+BACK (event type 3)converting photons and apply a cut on zenith angles less than 90degrees to e ff ectively remove photons originating from the Earthlimb. We used the recommended time selection, DATA QUAL>0&& LAT CONFIG==1 , and also apply a cut on rocking angle, (ABS(ROCK ANGLE)<52) .The ROI was defined as a square of 12 ×
12 deg aroundthe Coma center. As a baseline, the data were binned using apixel resolution of 0 . ff ects increase at low energy (see alsoSection 4.7). We start by modeling the ROI accounting for the di ff use back-ground components, corresponding to the isotropic and galac-tic interstellar emission, as well as the 4FGL compact sources.Given the size of our ROI, the most distant pixels are 8.5 de-grees away from the center. Nevertheless, we include all sourceswithin a square with a width of 20 degrees from the referencecenter, because the point spread function (PSF) may extend upto about 10 degrees at low energies.As discussed in Section 1, the 4FGL-DR2 catalog includes asource named 4FGL J1256.9 + = θ . This source is detected with a significance of 4.2,and is modeled by a point source with power law spectrum withindex 2.73 in the 4FGL-DR2 catalog. Its origin is uncertain, butit is given as possibly associated with NGC 4839 (with a prob-ability of 0.24). While this source could be a contaminant forthe di ff use ICM signal, as, e.g., the AGN of NGC 4839, it couldalso corresponds to the peak of the di ff use emission associatedwith the ICM. This second scenario is motivated by the uncer-tain origin of the source due to the limited signal to noise ra-tio and angular resolution of the Fermi-LAT. We thus aim attesting and comparing the two hypothesis. To do so, we con-sider the three following scenarios for modeling the Coma clus-ter region. 1) We assume that 4FGL J1256.9 + ff use emission. In this case,we include it in the ROI model as given by the 4FGL-DR2 cat-alog and do not include any extra di ff use emission associatedwith the ICM. 2) We assume that 4FGL J1256.9 + ff use cluster emission. In this case,the 4FGL-DR2 catalog model is not an appropriate descriptionof the signal and we replace 4FGL J1256.9 + ff useemission model associated with the ICM. 3) We assume that4FGL J1256.9 + ff use emission, but we also consider a di ff use ICM componentas γ -ray emission is expected from the cluster. In this case, weinclude 4FGL J1256.9 + ff use emission associated with the ICM.The di ff use Coma cluster ICM emission is modeled eitherusing the MINOT spatial and spectral physical templates, or usingthe
MINOT spectral template together with the spatial templatesbuilt from other wavelengths (see Section 3). We test di ff erentshapes for the profile of the CRp, and fix the CRp spectral indexslope to α CRp = .
8. This value will be later constrained usingthe measured SED, in Section 6. The nonoptimal choice of thisnumber will slightly reduce the significance associated with theComa di ff use emission, but does not a ff ect the SED fit performedafterward. As a baseline, and unless otherwise stated, we usethe model with η CRp = / Once the overall background model (large-scale di ff use pluspoint source emission) and cluster di ff use emission is defined,we use the following iterative procedure to fit the ROI. 1) Weconstruct a first sky model using the di ff use background and thepoint sources spectral parametrization, with parameters from the4FGL catalog, but exclude the cluster di ff use emission model atthis stage. 2) We run the optimize function of Fermipy , in-cluding all sources in the model (except for the cluster di ff useemission). We thus obtain a first renormalization of the sourcesin the model and a first value of their test statistics (TS, Mattoxet al. 1996), defined asTS = − L − ln L ) , (12)with L the maximum likelihood value for the null hypothesis,and L the maximum likelihood when including the additionalsource. 3) To allow for variability, we free all the model param-eters associated with sources with √ TS >
20, and free only thenormalization of sources with 5 < √ TS <
20. 4) We use the fit function of
Fermipy to perform the maximum likelihood
7. Adam et al.: γ -ray emission toward Coma F G L J . + ( n o r m ) F G L J . + ( n o r m ) F G L J . + ( n o r m ) F G L J . + ( P L i n d e x ) F G L J . + ( n o r m ) F G L J . + ( n o r m ) F G L J . + ( n o r m ) F G L J . + ( n o r m ) F G L J . + ( n o r m ) F G L J . + ( n o r m ) F G L J . + ( L P ) F G L J . + ( L P ) F G L J . + ( n o r m ) F G L J . + ( n o r m ) F G L J . + ( L P ) F G L J . + ( L P ) F G L J . + ( n o r m ) F G L J . + ( L P ) F G L J . + ( L P ) F G L J . + ( n o r m ) F G L J . + ( n o r m ) F G L J . + ( n o r m ) F G L J . + ( L P ) F G L J . + ( L P ) F G L J . + ( n o r m ) F G L J . + ( n o r m ) F G L J . + ( L P ) F G L J . + ( L P ) F G L J . + ( n o r m ) F G L J . + ( n o r m ) F G L J . + ( P L i n d e x ) F G L J . + ( n o r m ) C o m a ( n o r m ) g a l d i ff ( n o r m ) g a l d i ff ( P L i n d e x ) i s o d i ff ( n o r m ) C o rr e l a t i o n m a t r i x Figure 3:
Correlation matrix associated with the likelihood fit. Each entry corresponds to one free parameter of the model. The names isodi ff andgaldi ff correspond to the isotropic and galactic di ff use backgrounds, Coma corresponds to the ICM cluster model, and all the other names to the4FGL sources. The labels of the parameters are as follows: Norm stands for the spectrum normalization; PL index stands for the spectral index ofsources described by a power law spectrum; LP α and LP β stand for the spectral parameters of the sources described by a log-parabola spectrum.This correlation matrix corresponds to the case that includes 4FGL J1256.9 + MINOT cluster model with n CRp ∝ n / e . model fit of all the free parameters in the sky model. 5) We addthe cluster di ff use emission (except in the scenario 1, see sub-section 4.2) in our model, with its normalization being its onlyfree parameter, and refit the sky model as in step 4. In additionto the likelihood fit, we also use the tsmap Fermipy functionto compute TS maps, in the case of radially symmetric clustermodels.In Figure 3, we present the correlation matrix recoveredfrom the likelihood fit for all the fitted parameters includedin our model (i.e., the free parameters). It corresponds to thecase of the scenario 3 (see subsection 4.2), i.e., including boththe Coma cluster di ff use emission and the contribution from4FGL J1256.9 + ff use clusteremission is strongly degenerate (anti-correlated) with the nor-malization of 4FGL J1256.9 + + ff use emission is, on the other hand, notsignificantly correlated with any other compact source in themodel, but is marginally anti-correlated with the isotropic back-ground. We conclude that apart from the presence or not of4FGL J1256.9 + ff ects associatedwith the di ff use background).We present in Table 1 the TS values obtained for the pointsource 4FGL J1256.9 + ff use cluster model com-ponent individually, for all the cases considered (we also in-clude a point source model for the cluster as a spatial tem-plate for comparison). Because of the degeneracy between4FGL J1256.9 + ff use cluster model, we alsoinclude the TS value for the total (i.e., both components in-cluded or excluded when computing the likelihood L and L ,see Equation 12). Indeed, in the case of two degenerate sources,the maximum likelihood will not change much when removingone of the two source since the other source will absorb themissing component, in the null hypothesis. Nonetheless, remov-ing the two sources simultaneously may lead to a large changein the maximum likelihood because at least one component isnecessary to explain the data. Hence, the TS of each individ-ual degenerate source is expected to be low, but that of the twosources taken together might be large. This is expected to bethe case in scenario 3. All the considered MINOT radial mod-els give similar TS values, around TS = −
27, when re-placing 4FGL J1256.9 + + .
61 (sce-nario 1). On the other hand, the cluster emission modeled as apoint source gives a significantly lower value, of 17 .
54. Except
8. Adam et al.: γ -ray emission toward Coma for the radio relic, all the templates based on other data pro-vide a better match to the signal, with the best match reachedfor the galaxy density and tSZ maps (TS ∼ − γ -ray excess isobserved. In the case both 4FGL J1256.9 + ff usecluster emission are included in the ROI model (scenario 3), theTS of each component drastically reduces, highlighting the factthat 4FGL J1256.9 + ff use emission are degener-ate. For instance, the TS value reaches only 12 .
17 and 12 .
44 forthe cluster and 4FGL J1256.9 + + ff use cluster model (scenario 2)generally provides a better description of the data compared to4FGL J1256.9 + n CRp ∝ n / e ), in scenario 2, we obtain a best-fit photonindex of 2 . ± .
19 for a TS value of 27.37. As expected, theoverall spectrum is slightly flatter than in the case of our physicalmodel. Indeed, the photon spectrum vanishes at low energy fora given CRp energy slope with our physical model, thus leadingto a harder spectrum in this regime. The power law model thusaverages the two regimes. Nevertheless, we note that the highenergy photon spectrum does not strictly match the CRp energyspectrum (see Adam et al. 2020, for more details). The TS valueis nearly the same as in the case of the physical spectrum andgiven the available signal-to-noise ratio, it is not possible to dis-criminate the two.The modeling and fitting procedure described in Section 4.2and in this Section is also validated using a null test describedin Appendix C. This shows that no cluster detection is obtainedwhen including a cluster in the sky model close to other sourceswith similar background as around Coma.
In order to check how the data compare to the di ff erent mod-els (in the di ff erent scenarios), we computed maps of the data,model and residual (with and without Coma; with and without4FGL J1256.9 + ff use cluster emission and4FGL J1256.9 + >
16 and TS >
9, respectively). The peak of the excessis slightly o ff set in the southwest direction with respect to ourreference center. The position angle of this elongation agreesin all energy bins, although the elongation itself is more pro-nounced at low energy. When including 4FGL J1256.9 + + Fermi -LAT excess toother data will be discussed further in Section 5. We also show inAppendix B that the o ff set between 4FGL J1256.9 + + ff use ICM emission (i.e., scenario2). In addition to the central excess, we also see two other ex-cesses near RA,Dec = = > = / orspectral shape of this signal well. All these other excesses re-main relatively small and given their significance, it is not clearwhether they correspond to noise fluctuations, the mis-modelingof existing sources, or new sources that are not included in the4FGL catalog. Given their distances to the central excess, we donot expect that they would lead to any significance bias in ouranalysis.We also compute the counts as a function of energy observedwithin 3 × θ from the cluster center and compare it to the dif-ferent components of the model in Figure 5. As we can see, thedominant components are from the isotropic and di ff use galac-tic interstellar emission at almost all energies. The signal fromthe Coma cluster (in the case of scenario 2) is about an order ofmagnitude below depending on the energy, but is the dominantcompact source except at very low energies where the larger PSFleads to leakage from strong sources within the Coma central re-gion (we note that the source 4FGL J1253.8 + × θ ). Similarly in the case of scenario 1,4FGL J1256.9 + + ff erence is observed between the two sce-narios accounting for the error bars. Given this spectrum, it isvery unlikely that other compact sources from the ROI inducethe observed excess within 3 × θ of the Coma cluster.In order to address to which extent the Fermi -LAT data aresensitive to the shape of our radial model, we compute the profileof the data and model in Figure 6. Although each of the 5 firstbins (up to 1 degrees extent) are marginally inconsistent with themodel when the Coma cluster emission or 4FGL J1256.9 + σ , and thus correspond to a clear overall ex-cess. The baseline cluster model provides a good fit to the dataand significantly improves the residual (scenario 2). The pointsource cluster model, also shown for comparison, is too peakedwith respect to the observations. However, it is not in clear dis-agreement with the data given the error bars. This is in agree-ment with the results of Table 1, where the likelihood ratio be-tween the two model is ∆ TS = .
45, providing only a hint thatthe extended model is more appropriate than the point sourcemodel. In the case of scenario 1 (4FGL J1256.9 + Having a model for the sky in hand, we used the sed functionof
Fermipy to extract the SED of the Coma di ff use emission inthe di ff erent cases tested in this paper. The sed function fits forthe flux normalization of a source in independent bins of energy.To do so, we allowed the local photon slope to vary accordingto the MINOT global spectral model and we fixed the background
9. Adam et al.: γ -ray emission toward Coma
200 MeV - 300 GeV 200 MeV - 1 GeV 1 GeV - 300 GeV
Point source 4FGL J1256.9 + D e c . ( d e g r ee s ) C o un t s ( d a t a ) D a t a D e c . ( d e g r ee s ) C o un t s ( d a t a ) D e c . ( d e g r ee s ) C o un t s ( d a t a ) D e c . ( d e g r ee s ) C o un t s ( m o d e l ) M o d e l D e c . ( d e g r ee s ) C o un t s ( m o d e l ) D e c . ( d e g r ee s ) C o un t s ( m o d e l ) D e c . ( d e g r ee s ) C o un t s ( e x c e ss ) R e s i du a l ( w i t h o u t C o m a ) D e c . ( d e g r ee s ) C o un t s ( e x c e ss ) D e c . ( d e g r ee s ) C o un t s ( e x c e ss ) D e c . ( d e g r ee s ) C o un t s ( e x c e ss ) R e s i du a l (t o t a l ) D e c . ( d e g r ee s ) C o un t s ( e x c e ss ) D e c . ( d e g r ee s ) C o un t s ( e x c e ss ) Point source 4FGL J1256.9 + D e c . ( d e g r ee s ) C o un t s ( e x c e ss ) R e s i du a l ( w i t h o u t C o m a ) D e c . ( d e g r ee s ) C o un t s ( e x c e ss ) D e c . ( d e g r ee s ) C o un t s ( e x c e ss ) Figure 4:
Fermi -LAT imaging centered on Coma and comparison to the model (baseline model, n CRp ∝ n / e ). We show the Fermi -LAT data (firstrow), the best-fit model (second row), the residual excluding the cluster model (third row), and the residual with respect to the total model (fourthrow). We also show the residual excluding the cluster when accounting for 4FGL J1256.9 + = [4 , , , γ -ray emission toward Coma Table 1:
List of TS values and flux in the 200 MeV - 300 GeV band for all the models considered.
Scenario ( (cid:63) ) Sky model TS Cluster flux4FGL J1256.9 + ( † ) + − s − cm − )1 included None 25.61 25.61 – 02 replaced Point-source 17.54 – 17.54 11 . ± .
872 replaced Compact model ( n CRp ∝ n gas ) 24.84 – 24.84 13 . ± .
052 replaced Extended model ( n CRp ∝ n / ) 27.00 – 27.00 15 . ± .
252 replaced Flat model ( n CRp = constant) 25.11 – 25.11 18 . ± .
782 replaced Isobar ( n CRp = P gas ) 24.56 – 24.56 13 . ± .
012 replaced tSZ 32.18 – 32.18 17 . ± .
442 replaced X-ray 28.18 – 28.18 14 . ± .
202 replaced Galaxies 34.23 – 34.23 17 . ± .
452 replaced Radio halo 29.65 – 29.65 15 . ± .
252 replaced Radio relic 10.70 – 10.70 8 . ± .
993 included Point-source 32.24 18.75 7.60 4 . ± .
933 included Compact model ( n CRp ∝ n gas ) 34.67 14.29 11.50 6 . ± .
113 included Extended model ( n CRp ∝ n / ) 34.35 12.44 12.17 8 . ± .
233 included Flat model ( n CRp = constant) 34.67 13.02 9.32 10 . ± .
793 included Isobar ( n CRp = P gas ) 34.72 14.53 11.43 6 . ± .
073 included tSZ 34.35 8.06 15.79 10 . ± .
503 included X-ray 35.23 11.79 13.11 7 . ± .
263 included Galaxies 36.10 6.60 17.99 11 . ± .
513 included Radio halo 35.28 10.53 13.95 9 . ± .
313 included Radio relic 26.18 23.25 0.78 2 . ± . Notes. ( (cid:63) ) See subsection 4.2 for the definition. ( † ) The total model TS corresponds to the TS of both the cluster and 4FGL J1256.9 + + + C o un t s Model (all components)Model (without Coma or 4FGL 1256.9+2736)4FGL J1256.9+2736 only (scenario 1)Coma only (scenario 2)isodiffgaldiffData Energy (MeV)505 R e s i d u a l ( N / N ) Figure 5:
Fermi -LAT counts as a function of energy, computed within3 θ from the Coma center. The blue line shows the total model inthe case of scenario 2 (but it is not distinguishable from scenario 1 inthis figure), the red line shows the model when excluding the Comacluster di ff use emission or 4FGL J1256.9 + + ff use contribution in the case ofscenario 2. The contribution from the di ff erent sources in the ROI is alsoindicated as gray lines, except for the isotropic and di ff use backgrounds,given in magenta and purple, respectively. The bottom panel show theresidual between the data and the model, in red when both the Comacluster di ff use emission or 4FGL J1256.9 + component. However, we note that only minor di ff erences areobserved when fixing the slope or leaving the background free.In addition to the flux and error bars in each bin, the sed method Radius (deg)050100150200 P r o f il e ( c o un t s / d e g ) Cluster point source modelCluster baseline model ( n CRp n e )4FGL 1256.9+2736 modelData-Model (without Coma or 4FGL 1256.9+2736) Radius (deg)505
Figure 6:
Radial profile of the
Fermi -LAT excess count in the total en-ergy band (200 MeV - 300 GeV) for di ff erent models. The black datapoints give the excess count when both 4FGL J1256.9 + + ff use model only; baseline in greenand point source in blue). The shaded areas correspond to the expectedPoisson uncertainties given the respective model. The bottom panel pro-vides the residual normalized by the error bar with a similar color code. provides the full likelihood scan for the normalization value ineach bin, and we will use this information later in Section 6.At this stage, we use the SED results to compute the flux,integrated between 200 MeV and 300 GeV, of the Coma di ff useemission. The measured values are listed in Table 1 for the di ff er-ent models and scenario which we test. We note that in the caseof scenario 1 (i.e., no cluster in the model), the flux is by defini-tion equal to zero. The flux ranges from about 13 to 19 × −
11. Adam et al.: γ -ray emission toward Coma F l u x ( p h s c m ) All data fluxLightcurve T S / Figure 7:
Light curve associated with the Coma cluster model fit, forthe full dataset, in bins of about 4.5 months. for clarity, error bars areshown only for points that are mode than 1 σ away from zero, and up-per limits (95% confidence interval) are also given. The bottom panelpresents the square root of the TS associated with the source. ph s − cm − for the radially symmetric models in scenario 2. Themultiwavelength templates lead to similar fluxes, except for theradio relic template for which the flux is only about 9 × − phs − cm − . In case of scenario 3, part of the flux is absorbed in the4FGL J1256.9 + In order to check that the observed signal is consistent with dif-fuse ICM emission, we compute the light curve of the signal.Indeed, while many γ -ray sources are variable, the γ -ray emis-sion from galaxy clusters is expected to be steady, at least overthe timescale of any human observation, and observing a burstwould allow us to rule out the cluster ICM origin of the signal.We thus compute the light curve of the Coma cluster di ff use sig-nal using the lightcurve function from Fermipy , in the case ofscenario 2. We extract the flux within 30 bins, which correspondsto about 4.5 month per bin and report the result in Figure 7.Given the significance of the signal for the full dataset, we onlyexpect upper limits in each bin. As we can see, no significantburst is observed. The upper limit excludes the expected flux at95% confidence limit in one out of 30 bins (i.e., fewer than 5%)in agreement with expectation. Therefore, the light curve is con-sistent with the observed signal being associated with the di ff useComa cluster ICM emission. Our
Fermi -LAT analysis relies on choices that are somewhat ar-bitrary. To further validate the
Fermi -LAT results, we thus checkthe impact of various systematic e ff ects associated with thesechoices. They are listed below and summarized in Table 2 interms of the changes on the flux of the cluster emission. We alsorefer to Xi et al. (2018) and Ackermann et al. (2016) who testedthe impact of similar systematic e ff ects in the region aroundComa. As discussed below, the uncertainty in the di ff use back-ground is the dominant systematic e ff ect; it is expected to remainbelow 40%. Diffuse background emission
Although Coma is located nearthe galactic north pole (i.e., in a very clean region regarding dif-fuse galactic emission), the di ff use background is the dominantcontaminant (see Figure 5) and is slightly correlated with thecluster signal (see Figure 3). Given the location of Coma, wealso note that the galactic emission is nearly isotropic in the ROI,which explains the high degeneracy between the di ff use isotropiccomponent and the di ff use galactic component. Therefore, thedi ff use background modeling might lead to a significant sys-tematic e ff ect. The di ff use background related systematic e ff ectwas investigated by Ackermann et al. (2016) who showed thatabove 300 MeV, the systematic e ff ect was less than 22% (al-though could reach ∼
50% below 300 MeV). Similarly, Xi et al.(2018) concluded that the uncertainty in the di ff use backgroundwas <
30% within 0.2-300 GeV.To test the impact of the di ff use background, we first repro-duce our results by fixing the background to its expected modelvalue (i.e., normalization set to 1 and no spectral index variationallowed). In this case, we obtain TS = . ff set, which is partially ab-sorbed by the cluster model and explains the boost of the signal.Such values for the Coma cluster can thus only be taken as upperlimits.We then reproduce our results using previous backgroundmodels, namely gll iem v06 and iso P8R3 SOURCE V2 for thegalactic and isotropic emission, respectively. In scenario 2, theTS value slightly increases, to 31.93, and the flux remains stablewithin < , given the parameter definition files provided byAckermann et al. (2012). We obtain a set of 16 versions of thedi ff use galactic background components (and the correspondingisotropic di ff use emission that was model simultaneously) in-cluding bremsstrahlung, inverse Compton and pion decay emis-sion. We first fit the sum of the galactic background componentswith a single free normalization and spectral slope parameter.In this case, the TS value is systematically higher than in ourbaseline fit, ranging from 33.19 to 39.26, with a flux from 20 to40% higher. Then we consider a free normalization and spectralslope for each individual components separately, but do not con-sider the subcomponents separately (e.g., inverse Compton fromoptical, far infrared and CMB scattering individually). The TSvalues range from 30.98 to 39.42 and the flux is stable within[ − . , + . Energy range
We test the impact of our choice of the consid-ered energy range. As the systematic e ff ects are expected to bedominant at low energy, we test changing the nominal minimumenergy of 200 MeV to 500 MeV. We also consider the case of100 MeV as the low energy threshold. We extrapolate the fluxassuming our physical model with a CRp spectral index in therange α CRp = [2 . − . See https://galprop.stanford.edu/ , and https://galprop.stanford.edu/webrun.php for the web interface.12. Adam et al.: γ -ray emission toward Coma Table 2:
Estimation of the systematic e ff ects associated with the Fermi -LAT analysis.Type Range explored Change in √ TS Change in the total fluxDi ff use background alternative models 20% 40%Energy threshold ( (cid:63) ) / pixel 1.9% 0.9%ROI size 8-15 degrees 0.9% 12.8%4FGL source selection size 20-30 degrees 0.8% 1.8%Event selection alternative selections 6.3% 12.8%Rocking angle cut yes / no 3.5% 3.6% Notes. ( (cid:63) ) The flux variation is computed by extrapolating the model in the range 200 MeV - 300 GeV. The flux uncertainty is largely driven byuncertainties in the extrapolation. we note that this large number is mainly driven by the uncer-tainty in the extrapolation when converting the flux to the 200MeV - 300 GeV band, and the significance of the detection re-mains stable within 6.3%.
Energy and spatial binning
Our default binning choice is 0.1degree per pixel, and eight energy bins per decade. We vary thepixel resolution from 0.05 to 0.2 degree and the energy binningfrom 5 to 12 bins per decade. The changes are less than 6% over-all.
ROI size
We use a ROI size of 12 by 12 degrees. To make surethat this choice does not introduce any significant bias, e.g., dueto the presence of poorly constrained sources at the peripheryof the ROI, we reproduce the analysis using a ROI of 8 and 15degrees. This does not change our results by more than 13%.
Similarly, the choice of includingsources from the 4FGL catalog that are within a 20 degrees widthregion from the ROI center is checked by increasing this valueby 50%. The changes in our results are less than 2%.
Event class
We reproduce our results by applying the
P8R3 ULTRACLEANVETO V2 selection (event class 1024, i.e., thecleanest Pass 8 event class). We also use the
FRONT only (eventtype 1) converting photons. The changes are less than 13%.
Rocking angle cut
It is recommended by the
Fermi -LAT col-laboration to check the impact on the rocking angle cut on the re-sults. As a baseline, we apply a cut on rocking angle less than 52degrees. We reproduce the results presented here without consid-ering this cut. While the statistics slightly increases, the changesare less than 4% on the results.
In Section 4, we have shown that a significant γ -ray signal is ob-served within the characteristic radius θ of the Coma cluster.This agrees with the results by Xi et al. (2018) and the 4FGL-DR2 catalog (Abdollahi et al. 2020; Ballet et al. 2020).The source 4FGL J1256.9 + ff use cluster models that we have tested, and could correspondto the peak of the di ff use ICM emission. The comparison be-tween the models based on a single cluster di ff use ICM compo- nent and a single point source (4FGL J1256.9 + ff use models based on multi-wavelength templates provide the best agreement with the data.Models with two components (di ff use ICM plus point source)better match the data, but the improvement is marginal given theadditional component involved.In the next sections, we will explore the consequences forthe cosmic ray population in the Coma cluster of the scenariosin which the signal is (at least partly) associated with the di ff useICM. We will thus consider the scenarios 2 in which the sig-nal is entirely associated with the di ff use ICM, and the scenario3 in which the di ff use ICM and an independent point source(4FGL J1256.9 +
5. Multiwavelength comparison
In this Section, we qualitatively compare the
Fermi -LAT excessmap obtained in Section 4 to data at other wavelengths, as de-scribed in Section 2. The interpretation assumes that all the γ -rayemission arises from the di ff use ICM component.First, we briefly compare the Fermi -LAT TS map (baselinemodel) to the optical image constructed using SDSS data inFigure 8. This provides a visual appreciation of the scales probesby
Fermi -LAT. The galaxy density contours are also shown forvisual purpose. As we can see, the TS peak is located about 10arcmin north of the NGC 4839 group and about 20-30 arcminfrom the Coma center and its two brightest galaxies. However,this o ff set is small compared to the Fermi -LAT angular resolu-tion. The TS value remains larger than 16 (about 4 σ ) in most ofthe region where the bulk of the galaxy is located. However, the Fermi -LAT excess also extends further in the southwest direc-tion, as discussed below.In Figure 9, we compare the TS map to images of thetSZ, X-ray, galaxy density, and radio signal. The northwest-southeast elongated morphology of the TS map (and excesscounts) matches well what is seen at other wavelengths. The bestmatch is observed for the galaxy distribution and the tSZ map.The later being proportional to the product between the thermalgas density and temperature, this indicates that the CRp distribu-tion matches the temperature well. The temperature being fairlyconstant up to ∼ Mpc scales, this would favor relatively flat CRpdistributions. Alternatively, it could favor a scenario in which thesignal comes from the sum of unresolved sources, as traced bythe galaxies within the cluster region, and is not necessarily asso-
13. Adam et al.: γ -ray emission toward Coma D e c . ( d e g r ee s ) Figure 8:
SDSS color image of the central 2 × cluster region.The thin gray contours provide the galaxy density as shown in Figure9. The white contours give the TS map levels of 2, 4, 9, 16 and 25. Theimage was constructed using the g , r and i filters of SDSS. ciated with the di ff use ICM component. The X-ray morphologyis more compact than the Fermi -LAT excess, although the di ff er-ence could be largely due to the relatively poor angular resolu-tion of Fermi -LAT. The γ -ray excess matches the radio halo, butalso extend toward the relic and could thus provide a good matchto the combination of the two. In the case of an association withthe relic, it would suggest that a significant fraction of the signalarises from inverse Compton emission toward the relic becausevery little target thermal gas is available for hadronic interac-tion in the relic region, but we leave this interpretation for futurework (see Brown & Rudnick 2011; Ogrean & Br¨uggen 2013;Akamatsu et al. 2013, for discussions about a possible shock atthe relic location and the accretion shock interpretation for itsorigin).As noted in Figure 8, the TS peak is slightly o ff centeredwith respect to the cluster center, and better coincides with thesouthwest extension associated with the merger with the NGC4839 group (between the halo and the relic). It also coincideswell with the location where a radio front was identified (Brown& Rudnick 2011). This front is coincident with a discontinuityseen in the X-ray surface brightness, temperature and entropy(Simionescu et al. 2013; Uchida et al. 2016; Mirakhor & Walker2020) and SZ signal (Planck Collaboration et al. 2013), possi-bly suggesting that the local CR might have been accelerated bya shock. However, given the significance of the excess and the Fermi -LAT angular resolution, this remains di ffi cult to interpretfurther, as also shown in Section 4 and in particular in Figure 6.In conclusion, we find good morphological agreement withdata at other wavelengths. Nevertheless, the interpretation re-mains di ffi cult due to the low significance of the excess and the Fermi -LAT angular resolution. This comparison does not allowus to exclude that the signal observed is entirely, if any, associ-ated with the di ff use ICM emission. It could also be the sum ofseveral components.
6. Implications for the cosmic ray content of theComa cluster
In this Section, we use the γ -ray SED extracted in Section 4,together with the model described in Section 3, in order to con-strain the CR population in the Coma cluster. Although they pro-vide the best match to the data, the multiwavelength data tem-plates are not used because they do not allow us to have a three-dimensional physical model of the cluster, which is needed toconstrain the CR content. We aim at using the
Fermi -LAT extracted SED to constrain theCRp population of our model. As discussed in Section 3, and inmore details in Adam et al. (2020), the hadronically induced γ -ray emission depends on the thermal gas (modeled via the pres-sure and density), and the CRp population (spatial and spectraldistribution). The thermal gas pressure and density are well con-strained from Planck and
ROSAT data, respectively, and are thuskept fixed in this analysis. Since the CRp spatial distribution waskept fixed when extracting the SED, we keep it fixed when con-straining the parameters based on the SED fit. However, we con-sider the di ff erent spatial shapes as in Section 4 because we haveseen that it was not possible to discriminate between the di ff er-ent cases. We also consider the case of scenario 3, where boththe di ff use ICM and 4FGL J1256.9 + X CRp ( R ), which we define at aradius R ; 2) the slope of the CRp energy spectrum α CRp . Thesetwo parameters should be constrained for all the di ff erent mod-els considered (compact, extended, flat, and isobar), which willprovide an assessment of the systematic e ff ect associated withthe spatial model.We employed a Markov chain Monte Carlo (MCMC) ap-proach in order to constrain the parameters space, using the emcee package (Foreman-Mackey et al. 2013). In brief, thechains move in the parameter space according to a proposalfunction and the likelihood of the model given the data. Weadopted flat priors across the range X CRp ( R ) ∈ [0 , .
2] and α CRp ∈ [2 , ff ectour results. This method allows us to e ffi ciently find both thebest-fit parameters (as the parameters that maximize the likeli-hood) and the estimate of the posterior probability distribution.The likelihood function is defined asln L ( (cid:126)θ | D ) = (cid:88) i ln L i ( M i ( (cid:126)θ )) , (13)where i runs over each energy bin and the model parameters are (cid:126)θ ≡ [ X CRp ( R ) , α CRp ]. The value L i is the likelihood of mea-suring a given flux normalization in the energy bin i , dependingon the model flux M i ( (cid:126)θ ). It is obtained by interpolating the like-lihood scan, as provided by Fermipy , when extracting the SEDin Section 4.5.Once the chains have converged, and after removing theburn-in phase, the two-dimensional histogram in the plane X CRp ( R ) − α CRp can be integrated to provide the constraints upto a given confidence interval. The individual chain histogramsalso provide the marginalized posterior probability distributionof each parameter. The integrated posterior probability distribu-tion up to 68% probability gives the uncertainties. In addition
14. Adam et al.: γ -ray emission toward Coma NGC 4839 group
Halo FrontRelic
Figure 9:
Multiwavelength morphological comparison of the Coma cluster signal to the
Fermi -LAT TS map obtained in our baseline model.
Topleft : Planck tSZ.
Top Right : ROSAT
X-ray.
Bottom left : SDSS galaxy density.
Bottom right : WSRT 352 MHz radio signal. The field of view ofall images is 5 × . The white contours give the Fermi -LAT TS map (contours at 4, 9, 16 and 25) for the reference
MINOT model ( n CRp ∝ n / e ).For all panels, the black contours correspond to the maximum of the image divided by 2 i , with i the index of the contours. The dashed graycircle provides the radius θ and 3 × θ . Several relevant features are also indicated in orange. For display purposes, the WSRT image has beenapodized at large radii to reduce the larger noise fluctuations present on the edge of the field. As a complementary figure, Figure 8 provides anoptical image of the central region. to the model parameters themselves we obtain the constraint ofthe spectrum. To do so, we compute the model SED for eachset of parameters and calculate the envelope of all the models asthe 68% confidence limit measured from the model histogramin each energy bin. The same procedure is applied to measurethe γ -ray flux (and luminosity) between 200 MeV and 300 GeVaccording to the MCMC sampling. The SED measured in the case of our baseline model is shownin Figure 10 for both scenario 2 and 3 (see subsection 4.2 forthe detailed definition). Error bars are the 1 σ error on the SEDas evaluated from the likelihood scan curvature. The maximumlikelihood model is also shown in black, as well as 100 mod-els randomly sampled from the MCMC chains, their median andthe associated 68% confidence interval. We can observe that thebest-fit model is in good agreement with the data in both cases, as also highlighted by the residual. The model is relatively wellconstrained around 300 MeV - 1 GeV, but the uncertainties re-main large at larger energies. The peak SED, around 500 MeV,reaches about 3 × − MeV cm − s − . In the case of scenario3 (right panel, both the ICM and point source in the model), wecan see that the spectrum amplitude is reduced and that errorcontours are significantly larger.The corresponding constraint on the posterior likelihood ofthe parameters X CRp ( R ) and α CRp is shown in Figure 11. Theconstraints on the model SED lead to a relatively tight con-straint on the normalization, but the constraint on the slope re-mains fairly loose. The two parameters are degenerate as in-creasing the normalization and the slope simultaneously doesnot strongly change the flux at high energies (see also Figure 1,bottom right panel). The constraint on the CRp to thermal en-ergy is about 1.8% and the slope about 2.8, as also shown onthe marginalized distributions. When both the ICM componentand 4FGL J1256.9 +
15. Adam et al.: γ -ray emission toward Coma E d N d E dSd t ( M e V / c m / s ) Maximum likelihood modelMedian68.0% CLData 10 Energy (MeV)202 10 E d N d E dSd t ( e r g / c m / s ) E d N d E dSd t ( M e V / c m / s ) Maximum likelihood modelMedian68.0% CLData10 Energy (MeV)202 10 E d N d E dSd t ( e r g / c m / s ) Figure 10:
Total SED recovered from the
Fermi -LAT and MCMC constraint in the case of the reference model ( n CRp ∝ n / e ) when4FGL J1256.9 + ff erence between the data and the best-fit model normalizedby the error bars. Table 3:
Constraints on the CRp population and associated flux and luminosity in the case of the radial models. The quoted values and uncertaintiescorrespond to the maximum likelihood and 68% confidence interval of the distributions.Model X CRp (%) α CRp
Flux (10 − ph s − cm − ) Luminosity (10 erg s − )4FGL J1256.9 + n CRp ∝ n gas ) 1 . + . − . . + . − . . + . − . . + . − . Extended model ( n CRp ∝ n / ) 1 . + . − . . + . − . . + . − . . + . − . Flat model ( n CRp = constant) 1 . + . − . . + . − . . + . − . . + . − . Isobar ( n CRp = P gas ) 1 . + . − . . + . − . . + . − . . + . − . Both 4FGL J1256.9 + n CRp ∝ n gas ) 0 . + . − . . + . − . . + . − . . + . − . Extended model ( n CRp ∝ n / ) 1 . + . − . . + . − . . + . − . . + . − . Flat model ( n CRp = constant) 0 . + . − . . + . − . . + . − . . + . − . Isobar ( n CRp = P gas ) 0 . + . − . . + . − . . + . − . . + . − . straint on the normalization shifts toward zero and the posterioronly excludes X CRp ( R ) = σ . The uncertainty onthe slope increases and the best-fit slightly decreases to about α CRp (cid:39) . + The presence of CRp, as traced by the hadronic γ -ray emis-sion, implies the production of secondary CRe, which shouldcontribute to the observed radio emission. We compute both theprofile and the spectrum associated with this population giventhe modeling discussed in Section 3. We stress that this is done assuming a steady state scenario. This calculation is done for the Fermi -LAT SED best-fit model, as well as for the set of param-eters sampled in the MCMC to obtain the 68% confidence limiton the radio emission.As discussed in Section 3, the primary CRe, i.e., CRe , arealso included in our model and we consider two cases for in-terpreting their origin. 1) Since our model does not include ex-plicitly any reacceleration (see Brunetti & Lazarian 2007, 2011,for reacceleration models), the population of CRe that we con-strain is supposed to be independent of the CRp and the CRe .It would correspond, for instance, to a CRe population arisingfrom star formation activity spread over the cluster volume, orthe direct shock acceleration in the ICM. In this case, the CRe and CRe populations coexist and the radio emission accountsfor the sum of the two. 2) We could also interpret the CRe pop-ulation in our model as the reaccelerated CRe population. Inthis case, the CRe would be the seed for the CRe . The CRe population should thus not contribute directly to the radio emis-sion. Instead the ratio CRe / CRe would measure the amount ofreacceleration needed to explain the emission with respect to thepurely hadronic steady state scenario, as a function of energyand radius. We stress that in this second interpretation, we only
16. Adam et al.: γ -ray emission toward Coma ICM model only ( n CRp n e )ICM model ( n CRp n e ) + 4FGL J1256.9+2736 X CRp C R p CRp
Figure 11:
MCMC constraints on the model parameters, in the case of the reference model ( n CRp ∝ n / e ) when 4FGL J1256.9 + X CRp ( R ) − α CRp , where the contourscorrespond to 68% and 95% confidence interval. The marginalized posterior probability distributions are also shown for the parameters X CRp ( R )(top) and α CRp (right), where the shaded area provides the 68% confidence interval. provide a constraint relative to the hadronic steady state scenariosince this is an assumption made when computing the CRe pop-ulation. In contrast, reacceleration models do not assume steadystate, and both CRp and CRe populations evolve together ac-cording to turbulent reacceleration. Nevertheless, we considerthis second case as it is still instructive regarding reaccelerationphysics.The CRe are parametrized using the ExponentialCutoffPowerLaw model (our baseline), the
InitialInjection model, or the
ContinuousInjection model, as described in Section 3.1.4. All models include anormalization X CRe ( R ), a spectral slope α CRe , a breakenergy ( E cut , CRe ), and a spatial scaling relative to the thermaldensity ( η CRe ). We fit for these parameters using simultaneouslythe radio spectrum and the 352 GHz radio profile data (seeSection 2). Regarding the spectrum, we compute our modelusing a cylindrical integration within R = . × R (629kpc, 0.36 deg) as it corresponds to the extent of the radio halo(Brunetti et al. 2013). We note, however, that the radio spectrumdata are not strictly homogeneous in terms of aperture radiusused for flux measurement and our model value only providesan e ff ective radius. Because the profile is extracted from a single instrument (WSRT), which arguably could be the best one interms of capturing the di ff use emission, but which does notnecessarily perfectly agree with other measurements (as seen inFigure 12), we also allow for a cross-calibration of the profilemeasurement by adding an extra normalization which we fitsimultaneously. The fit is performed with the least square func-tion curve fit from the scipy python package. Dependingon the considered case, the radio model includes both the CRe and CRe contributions (case 1, no reacceleration), or only theCRe as the reaccelerated CRe (case 2, pure reacceleration).We compute the error contour on the CRe fitted population byreproducing the fit in the case of the lower and upper boundsfor the CRe population. We thus assume that the uncertaintiesfrom the CRe population (given by the γ -ray) are dominantover the uncertainties associated with the radio data.In Figure 12, we show the constraint on the CRe populationfrom the radio synchrotron spectrum and profile fit in the firstcase using the ExponentialCutoffPowerLaw model. We notethat the Figure is provided for our baseline CRp radial model( n CRp ∝ n / ) in the scenario in which 4FGL J1256.9 +
17. Adam et al.: γ -ray emission toward Coma (MHz)10 F ( J y ) CRe1CRe2TotaltSZ absolute valueWSRT 352 MHzData 10 Radius (kpc)10 S ( J y a r c m i n ) CRe1CRe2TotalData10 Radius (kpc)10 C R d e n s i t y ( c m G e V ) CRe1, 1 GeVCRe2, 1 GeVCRe1, 10 GeVCRe2, 10 GeV 10 Radius (kpc)10 C R e / C R e Figure 12:
Constraint on the CRe populations in the case of scenario 1 (distinct CRe and CRe populations, no reacceleration). Top left :Radio spectrum of Coma, as compiled from Pizzo (2010) and constraint from the reference CRp spatial model ( n CRp ∝ n / e ) and the referenceCRe spectral model ( ExponentialCutoffPowerLaw ). The measurement from the Brown & Rudnick (2011) data is also shown as the magentadiamond. The contribution from the CRe is shown in blue together with its 68% confidence interval, and the remaining contribution from CRe is shown in green. The sum of the two is given as the black line. The dashed gray line provide the expected amplitude of the tSZ signal. All fluxesare computed using cylindrical integration within R = . × R =
629 kpc ≡ .
36 deg.
Top right : Radio profile measured from the WSRT mapat 352 GHz and comparison to the reference model. The contributions from CRe and CRe are as in the left panel. Bottom left : Absolute numberdensity distributions of CRe and CRe taken at 1 GeV and 10 GeV. Bottom right : Ratio between the CRe and CRe number populations. Wenote that in the case of this figure, the confidence limits where computed using a resampling of only 100 parameters, and are thus not very accurate.We also stress that this figure depends on the magnetic field modeling (see Section 7.2 for discussions). els considered. First, we note that the tSZ signal, included as adashed gray line given our pressure model, is not negligible forthe highest frequency data point, but we correct for it (see alsoBrunetti et al. 2013, for a dedicated analysis). Our model pro-vides a reasonable fit to the data for both the spectrum and theprofile (this is also the case for the other considered models, seeAppendix D, and our results do not depend significantly on theconsidered CRe spectral model). The slope of the synchrotronemission from CRe is similar to the one from the CRe , but itis significantly less curved and no high energy cuto ff is presentin the CRe distribution. We can see on the spectrum that the ra-dio emission within 0 . × R (629 kpc, 0.36 deg) from CRe is overall a factor of about four lower than the total emission(except at high frequency where it reaches similar values). Onthe profile, the CRe emission is significantly more concentratedthan the total radio signal in the case of this CRp spatial model(although it is also true for all CRp model, see Appendix D) andlead to slightly over-fitting the total synchrotron emission in thecore when added to the CRe contribution. This high concentra-tion is expected because the CRe profile arises as the product of the gas density and the CRp density. The synchrotron profile,which arises from the product of the magnetic field profile andthe total CRe density, is so flat that it requires a nearly flat CRedistribution given the fixed magnetic field profile (as also notedin Zandanel et al. 2014). Thus, this could be achieved for theCRe only at the cost of an inverted CRp profile (rising with ra-dius). The number density of CRe is at a level of about 10 − and 10 − cm − GeV − at 1 GeV and 10 GeV, respectively, forboth populations, but with opposite radial dependences. Theseconstraints on the CRe populations translate into a ratio CRe toCRe that increases with radius and that does not depend muchon energy up to 10 GeV (before the cuto ff ; best-fit E cut , CRe = / CRe ratio isslightly below unity in the core, and strongly rises to reach about100 at 2 Mpc.In Figure 13, we show the constraint on the CRe populationin the second case (CRe are the seed to the CRe , pure reaccel-eration). As for Figure 12, the alternative models considered inthis paper are shown in Appendix D. The model also providesa good fit to the data as shown in the top panels, where only
18. Adam et al.: γ -ray emission toward Coma Figure 13:
Same as Figure 12 in the case of scenario 2 (CRe interpreted as the reaccelerated CRe population). The CRe prior reaccelerationrefer to the steady state hadronic model. The reacceleration boost is given relative to the steady state hadronic model. the CRe population (as the reaccelerated CRe population) isincluded. On the bottom panels, the interpretation is now di ff er-ent as the amount of CRe now correspond to the CRe popula-tion after reacceleration. As can be observed, the best-fit CRe profile is nearly flat, while the original seed population is moreconcentrated in the core. The amount of reacceleration relativeto the hadronic steady state case is thus relatively low in the core(in agreement with no reacceleration there, or even favoring aboost lower than one for the most compact CRp profile), butstrongly increases in the outskirt reaching about a factor of 100 at2 Mpc. As in the previous case, this radial dependence dependson the considered model for the CRp distribution and we showin Appendix D that flatter is the CRp distribution, flatter will bethe reacceleration boost profile. Nevertheless, in all the consid-ered cases, the reacceleration boost (relative to the steady statehadronic model) increases with radius. The energy dependenceof the boost, comparing the values for 1 GeV and 10 GeV elec-trons, is not much a ff ected by the choice of the spectral modelfor the CRe population.To summarize how the radio data were used, we used theWSRT data to establish the radio profile, but for the spectral de-pendence, we were forced to use only the data within 0 . × R where measurements at other frequencies were available. Thelarger halo size and flux from the WSRT data, which are due toits increased sensitivity, are discussed in detail by Brunetti et al.(2013). Those authors showed (in their Figure 2, right) that theobserved size of the halo was a strong function of the signal to noise ratio, and that the WSRT data were therefore the most reli-able. The magenta diamond in Figure 12 and 13 shows how theWSRT flux is above the rest of the spectral points. In order toexamine the e ff ects of increasing the entire spectrum by a factorof two to match the WSRT point, we added an extra normaliza-tion parameter, as discussed earlier in this Section. We find thatthe ratio between the total radio emission and that arising fromCRe would increase from 4 to 8. However, the e ff ect of chang-ing such normalization would not impact the conclusions of thispaper.In Appendix D, we also provide the results when includingboth the ICM component and 4FGL J1256.9 + Fermi -LAT sky model (scenario 3). While the CRe synchrotron spec-trum is slightly steeper and the amplitude of the CRe compo-nent is reduced, the shapes of the CRe / CRe profile (case 1), orreacceleration boost profile (case 2) are not significantly changedgiven the large uncertainties.Finally, in order to compare the distribution of CRp that wemeasure to that expected in reacceleration models to explain theradio emission, we compare our CRe induced synchrotron spec-trum to that of the model developed by Brunetti & Lazarian(2011) in Figure 14. To best match the model developed byBrunetti & Lazarian (2011), we use the ICM model in which theCRp radial profile is flat, but we note that the two models are notstrictly equal. For instance, Brunetti & Lazarian (2011) use anisothermal β -model, while we set the thermodynamic profiles toX-ray and tSZ data (see Section 4.3.1 from Brunetti & Lazarian
19. Adam et al.: γ -ray emission toward Coma ff (dashed brown line), only the emissionfrom secondaries directly produced from hadronic collisions re-mains, which compares very well to the prediction from ourmodel. This shows that the distribution of CRp that we measureprovides an excellent match to what is needed in the reaccelera-tion model developed by Brunetti & Lazarian (2011) to explainthe overall radio spectrum, when including reacceleration.
7. Discussions
Constraints on the γ -ray emission of the Coma cluster have beenobtained in earlier work using Fermi -LAT data (see Arlen et al.2012; Prokhorov 2014; Zandanel & Ando 2014; Ackermannet al. 2016; Keshet & Reiss 2018; Xi et al. 2018). In general, thesignal is modeled using spatial templates together with powerlaw distribution for the photon spectrum. Our results push for-ward such analysis by directly using a physical model for theCRp and thermal gas. Thus we are able to directly constrain thephysics of CR by comparing model and data, and push the anal-ysis by investigating the implications of γ -ray emission for theCRe population. Our results are in broad agreement with previ-ous searches and limits. Nevertheless, we stress that the reportedCRp to thermal energy (or pressure) are generally obtained underthe assumption of a given spectral distribution (harder the spec-trum and more stringent the limit). Spectral index values usedin the literature correspond to spectra that are generally muchharder ( ∼ . − .
3) than what we obtain here (see also below forfurther discussions).The
Fermi -LAT collaboration observed γ -ray excess withinthe virial radius of Coma using six years of data (Ackermannet al. 2016). However, the signal was too faint for detailed inves-tigation and they published upper limits on the signal for varioustemplates. Our results are consistent with that of Ackermannet al. (2016), although using twice the amount of data and ananalysis that di ff ers in various aspects. They provide an upperlimit of 5 . × − ph cm − s − ( E >
100 MeV). Extrapolatingour best-fit model in the same energy range, we obtain a com-patible total flux of 2 . × − ph cm − s − (in scenario 2, tomatch Ackermann et al. 2016). Among the main di ff erences, wenote that they use a power law emission model with a spectralindex of 2.3 and spatial distribution based on WSRT, while ourbaseline model is directly connected to the underlying CR popu-lation and the spectrum that we measure is significantly steeperthan their model. Given this framework, they obtain a TS valueof 13, compared to about 27 in our case.The detection of γ -ray emission toward Coma was firstclaimed by Xi et al. (2018), who used an unbinned likelihoodapproach (reaching TS values of about 40-50 depending on themodel). While the morphology of the signal and the spectralslope of the emission that we observe agree with their results, weobtain fluxes that are significantly lower (about a factor of two),depending on the exact model used. For instance, they obtainfluxes of about 2 . − . × − ph cm − s − for disk, core, or radioand X-ray based templates ( E = [200 , . − . × − ph cm − s − . We note that their fluxes would also be excluded bythe Ackermann et al. (2016) limit when extrapolating down to100 MeV using a photon power law index of 2.7-2.8, and wouldnearly reach the limit when extrapolating with our best-fit model. Keshet & Reiss (2018) claimed the detection of a ring-likesignal at a position that correspond to the expected location ofthe accretion shock (see also Hurier et al. 2019, for the detectionof such an accretion shock with Planck in A2319). In contrast,our results do not show any ring-like structure, especially whenlooking at the excess profile around 3 Mpc radius (about 2 de-grees).
The amount of CRp and their spectral and spatial distributionshas been predicted from numerical simulation (e.g., Pfrommeret al. 2007; Pinzke & Pfrommer 2010). The amplitude of theCRp that we obtain, relative to the thermal energy, is in line withpredictions from such simulations. We obtain a ratio of about1.5% (about 0.8% when including 4FGL J1256.9 + X CRp ( R = . = . × − (we account for the factthat the pressure ratio that they use is half the energy ration thatwe use). This value is just a factor of ∼ ff erent CRp indexand may also depend on radius so that the comparison dependson the exact CRp spatial distribution. This CRp to thermal en-ergy ratio is related to the CRp injection e ffi ciency for the di ff u-sive shock acceleration process, which increases with the Machnumber. Thus, it gives access, in principle, to the microphysicsof CRp acceleration. We refer to Ackermann et al. (2014) formore discussions.As also found by Xi et al. (2018), the spectral index ofthe CRp that we constrain is significantly larger than what isgenerally assumed ( α CRp ∼ . − . . − .
3; see e.g., Arlen et al. 2012; Zandanel & Ando2014; Pinzke & Pfrommer 2010). For a given shock, the CRpindex is related to the Mach number: higher the Mach numberand harder the spectrum (see Pfrommer et al. 2006, for discus-sions). Therefore, our results could point to shock accelerationfor which the Mach numbers are overall smaller than usually ex-pected.The shape of the CRp is often assumed to follow the shapeof the thermal density profile, with possible di ff usion sometimesincluded (e.g., Zandanel et al. 2014). While our results favormodels with intermediate scaling ( n CRp ∝ n ∼ / e ), the data arenot su ffi cient to firmly constrain the shape of the CRp distribu-tion. In the context of the reacceleration of CRe , a compact CRpprofile would imply that reacceleration strongly increases withincreasing radius, as already highlighted in Pinzke et al. (2017).However, even flat CRp profiles lead to an increasing reacceler-ation boost with radius.Our results are in agreement with earlier work that haveshown that pure hadronic models were excluded given the mag-netic field strength inferred from Faraday rotation measures(Brunetti et al. 2012; Zandanel & Ando 2014; Zandanel et al.2014; Brunetti et al. 2017). Nevertheless, the hadronically in-duced CRe are only a factor of a few lower than the amount re-quired, and could thus provide seeds for turbulent reaccelerationmodels (Brunetti & Lazarian 2007, 2011).Throughout this work, we have fixed the magnetic fieldstrength model to the best-fit result obtained by Bonafede et al.(2010), despite the fact that there are relatively large uncer-tainties in the constraint. We also refer to the recent work byJohnson et al. (2020) who showed that even under ideal condi-
20. Adam et al.: γ -ray emission toward Coma Figure 14:
Comparison between the CRe induced synchrotron spectrum to the reacceleration model developed by Brunetti & Lazarian (2011),in the case of a flat CRp population. The solid brown line corresponds to the full reacceleration model, while the dashed brown line corresponds tothe case where reacceleration is switched o ff (see Brunetti & Lazarian 2011, for more details). Left : Replacing 4FGL J1256.9 + ff use component (scenario 2). Right : Including both 4FGL J1256.9 + ff use component in the sky model (scenario 3). tions, the central magnetic field cannot be determined to betterthan a range of 3, with corresponding uncertainties in the scal-ing parameter η . These uncertainties a ff ect the results presentedin Section 6.3. In our work, increasing (decreasing) the mag-netic field would imply more (fewer) synchrotron emission for agiven CRe population. Thus it would lead to fewer (more) CRe for case 1, or fewer (more) reacceleration boost in case 2. Thise ff ect also applies to the radial dependence. Given the uncer-tainties in the magnetic field measure, these should contributesignificantly to our constraints. The uncertainties in the CRe population are nonetheless expected to dominate, and we leavethe joint investigation of the CRe population and magnetic fielddistribution aside. In this paper, we have generally assumed that the γ -ray emis-sion observed in the direction of the Coma was originating fromhadronic interactions in the ICM. We have also considered thecase where 4FGL J1256.9 + ff use component as fromhadronic interactions.In fact, cluster member and radio galaxies may also con-tribute significantly to the total signal. In Ackermann et al.(2016), a minimum flux was estimated for the two dominant cen-tral radio galaxies NGC 4869 and NGC 4874 assuming that theelectrons responsible for the radio emission also generate inverseCompton emission in the Fermi -LAT band and assuming simplescaling relations for the calculation. They obtained luminositiesof about 6 × and 2 × erg s − (0 . < E <
10 GeV), whichis a factor of ∼
20 lower than what we measure for the hadronicemission at E > . ff use component. In particular star forminggalaxies could lead to γ -ray emission for which the associatedflux is very uncertain and within the range of our constraints (inthe range 3 × − × erg s − for energies in 0.1-100 GeV;see Storm et al. 2012). The cluster member radio and star forming galaxies are alsoexpected to generate CR that di ff use in the ICM and would con-tribute to the population which we model in this paper. Rephaeli& Sadeh (2016) calculated that they might account for a sig-nificant amount of the radio di ff use emission, as well as in the γ -rays.
8. Summary and conclusions
This paper presented the analysis of nearly 12 years of
Fermi -LAT data toward the Coma cluster, together with mul-tiwavelength data and using the
MINOT software in order tomodel the signal. Di ff erent scenarios were considered to modelthe signal: no cluster di ff use emission, replacing the source4FGL J1256.9 + ff use cluster model, or accountingfor a combination of both. Assuming that the di ff use emissionwas associated with the hadronically induced γ -rays in the ICM,we investigated the implications for the CR physics, for bothCRp and CRe.The signal was modeled assuming that the γ -ray emissionarises from hadronic interactions between CRp and the thermalgas. The thermal gas model was set using ROSAT
X-ray and
Planck tSZ data. In this model secondary CRe are also pro-duced, and we fixed the magnetic field to that obtained fromFaraday rotation measurements in order to compute the asso-ciated radio synchrotron emission. The CRp spatial model wascalibrated assuming a scaling relative to the thermal gas pro-files. The CRp normalization and spectrum were defined relativeto the thermal gas energy and using a power law for the spec-trum, respectively. In addition to spherically symmetric models,we built two-dimentional spatial templates based on X-ray, tSZ,radio and galaxy density images to fit the
Fermi -LAT data.We detect γ -ray emission in the direction of the Coma clus-ter. The detection level depends on the model considered, in therange TS ∼ −
34, corresponding to a significance of about4 . − . σ . While extended models provide a better descriptionof the data, it is not possible to strictly exclude that the signal isassociated with a point source, or a combination of the two, andwe include this possibility in our analysis. The morphology ofthe signal is elongated in the northeast to southwest direction, inagreement with other wavelengths. The peak of the emission is
21. Adam et al.: γ -ray emission toward Coma about 0.5 degree o ff set in the southwest direction with respect tothe X-ray peak, and coincides with the well-know merger exten-sion.Using an MCMC approach, we constrained the amplitudeand the spectral index of the CRp population in the cluster as-suming that at least part of the signal was associated with the dif-fuse ICM hadronic interactions. We find that the energy stored inthe CRp is about 1.5% of the thermal energy of the Coma clus-ter (0.8% when including 4FGL J1256.9 + ff use shock acceleration, this could implies that CRp accelera-tion arises in weaker shocks than what is usually assumed.Secondary CRe are also expected from hadronic interactions.Their population was computed in a steady state scenario, lead-ing to a radio synchrotron emission that is about 4 times lowerthan what is observed. While a pure hadronic origin of the radioemission is ruled out, these secondary CRe could serve as theseeds for turbulent reacceleration. In this model, the reacceler-ation should increase with radius, depending on the exact CRpspatial distribution. Alternatively, an independent CRe popula-tion could be at the origin of the remaining radio emission, butit would require a nearly flat radial distribution.Our results show that after almost 12 years of observations,the di ff use γ -ray emission from galaxy clusters might now be-come accessible with the Fermi -LAT. Since the hadronic emis-sion is expected to be a universal property of galaxy clusters,our results might be reproducible for other clusters, renewing theinterest of such analysis, although Coma might be the best tar-get for such searches. In the scenario in which the signal indeedarises from di ff use ICM hadronic interactions, if the large valueof the CRp spectral slope is confirmed and if the CRp contentof clusters is, as expected, a universal property, the very highenergy γ -ray emission ( ∼ TeV) could therefore be much lowerthan usually assumed. In this context, it would be challengingfor future ground-based γ -ray observatories (e.g., the CherenkovTelescope Array, Cherenkov Telescope Array Consortium et al.2019) to detect the di ff use emission associated with hadronic in-teractions in the ICM of galaxy clusters. Acknowledgements. We are thankful to the anonymous referee for useful com-ments that helped improve the quality of the paper. We thank Gianfranco Brunettifor useful discussions and for providing us with the reacceleration model used inFigure 14. Partial support for LR comes from U.S. National Science Foundationgrant AST 17-14205 to the University of Minnesota. This research made use ofAstropy, a community-developed core Python package for Astronomy (AstropyCollaboration et al. 2013), in addition to NumPy (van der Walt et al. 2011), SciPy(Jones et al. 2001) and Ipython (P´erez & Granger 2007). Figures were generatedusing Matplotlib (Hunter 2007) and seaborn (Waskom & the seaborn develop-ment team 2020).
References
Abdollahi, S. et al. 2020, ApJS, 247, 33, 2005.11208Ackermann, M. et al. 2014, ApJ, 787, 18, 1308.5654——. 2016, ApJ, 819, 149, 1507.08995——. 2010, ApJ, 717, L71, 1006.0748——. 2012, ApJ, 750, 3, 1202.4039Adam, R. et al. 2020, arXiv e-prints, arXiv:2009.05373, 2009.05373Aharonian, F. et al. 2009, A&A, 502, 437, 0907.0727Aharonian, F. A., Kelner, S. R., & Prosekin, A. Y. 2010, Phys. Rev. D, 82,043002, 1006.1045Akamatsu, H., Inoue, S., Sato, T., Matsusita, K., Ishisaki, Y., & Sarazin, C. L.2013, PASJ, 65, 89, 1302.2907Aleksi´c, J. et al. 2010, ApJ, 710, 634, 0909.3267Arlen, T. et al. 2012, ApJ, 757, 123, 1208.0676Astropy Collaboration et al. 2013, A&A, 558, A33, 1307.6212Atwood, W. B. et al. 2009, ApJ, 697, 1071, 0902.1089 Ballet, J., Burnett, T. H., Digel, S. W., & Lott, B. 2020, arXiv e-prints,arXiv:2005.11208, 2005.11208Birkinshaw, M. 1999, Phys. Rep., 310, 97, arXiv:astro-ph / / / / ffi n, R. D., Dai, X., & Kochanek, C. S. 2014, ApJ, 795, L21, 1405.7047Huber, B., Tchernin, C., Eckert, D., Farnier, C., Manalaysay, A., Straumann, U.,& Walter, R. 2013, A&A, 560, A64, 1308.6278Hunter, J. D. 2007, Computing In Science & Engineering, 9, 90Hurier, G., Adam, R., & Keshet, U. 2019, A&A, 622, A136Hurier, G., Mac´ıas-P´erez, J. F., & Hildebrandt, S. 2013, A&A, 558, A118,1007.1149Johnson, A. R., Rudnick, L., Jones, T. W., Mendygral, P. J., & Dolag, K. 2020,ApJ, 888, 101, 2001.00903Jones, E., Oliphant, T., Peterson, P., et al. 2001, SciPy: Open source scientifictools for Python, http: // / Kafexhiu, E., Aharonian, F., Taylor, A. M., & Vila, G. S. 2014, Phys. Rev. D, 90,123014, 1406.7369Kelner, S. R., Aharonian, F. A., & Bugayov, V. V. 2006, Phys. Rev. D, 74,034018, astro-ph / / /
22. Adam et al.: γ -ray emission toward Coma Pfrommer, C., Springel, V., Enßlin, T. A., & Jubelgas, M. 2006, MNRAS, 367,113, astro-ph / / / / / / seabornWillson, M. A. G. 1970, MNRAS, 151, 1Wood, M., Caputo, R., Charles, E., Di Mauro, M., Magill, J., Perkins, J. S., &Fermi-LAT Collaboration. 2017, in International Cosmic Ray Conference,Vol. 301, 35th International Cosmic Ray Conference (ICRC2017), 824,1707.09551Xi, S.-Q., Wang, X.-Y., Liang, Y.-F., Peng, F.-K., Yang, R.-Z., & Liu, R.-Y. 2018,Phys. Rev. D, 98, 063006, 1709.08319York, D. G. et al. 2000, AJ, 120, 1579, astro-ph / Appendix A: Inverse Compton emission
This appendix provides an estimate of the inverse Comptonemission associated with CRe and CRe .Firstly, we assume a CRp population with X CRp ( R ) = − and α CRp = .
8, which is consistent with the values found inthis paper for Coma, and use our baseline scaling n CRp ∝ n / e .We compute the associated CRe population in the steady stateapproximation, as in Section 6. We finally compute the inverseCompton emission as explained in Adam et al. (2020). This pro-vides an estimate of the necessary inverse Compton emissionassociated with the hadronic collisions.Secondly, we assume that the CRe population is describedby a power law, which we match to radio data and extrapolate tovery high energies. As strong losses are expected at high energy,this provides an upper limit on the CRe population for energiesthat are large enough to induce inverse Compton emission inthe Fermi -LAT energy range (energies in the range of about 100GeV - 100 TeV, while the radio emission probe CRe up to afew tens of GeV at GHz frequencies). We match the slope andspectrum of the power law to the radio spectrum and computethe associated inverse Compton emission. These two cases (CRe and CRe ) are reported in Figure A.1.On the left panel, we see the contribution of the CRe and CRe to the radio signal. On the right panel, we see their contributionto the inverse Compton emission. Even in the case of boostingthe CRe population (i.e., the turbulent reacceleration model) bya factor of five to ten to match the radio emission, this wouldremain more than an order of magnitude below the hadroni-cally induced γ -rays. An energy-dependent boost, strongly ris-ing with the energy, could lead to the inverse Compton signalbeing significant, but this is not expected. In the case of CRe ,the upper limit provided by the power law model remains be-low the hadronically induced γ -rays at energies above 100 MeV.As losses are expected to strongly reduce the signal as energy in-creases, we expect the emission to be negligible, especially sincethe Coma radio spectrum already present a curvature at GHz fre-quencies.While these results may slightly vary depending on the exactshape of the spatial and spectral distribution of the CR, we do notexpect any major change in the comparison. We conclude thatthe Inverse Compton emission from both CRe and CRe shouldnot contribute significantly to the observed signal at Fermi -LATenergies, and it is neglected when modeling the γ -ray emission. Appendix B: Comparison to Monte Carlorealizations
In this Appendix, we quantify the likelihood of finding a pointsource 0.68 degrees away from the Coma cluster center (as4FGL J1256.9 + ff use cluster ICMemission (i.e., scenario 2). To do so we use an approach that isbased on Monte Carlo realizations, as follows. First, we selectthe best-fit sky model of scenario 2 (with 4FGL J1256.9 + localize method from Fermipy to search for the location which provides the best matchfor such a test point source within 1.5 degree. We repeat the oper-ation 200 times and compute the probability to find such a sourceat distance θ from the cluster.The results are provided in Figure B.1, where we show thebest-fit coordinates of the test point source for all the MonteCarlo realizations and compare it to the position of the Comareference center and 4FGL J1256.9 + ff sets between the recovered test point sourceand the Coma reference center. The results are shown in the caseof our baseline extended ICM model and in the case the spatialICM template is based on the tSZ map, for comparison. The re-sults for the other cases are expected to slightly vary based onthe extension of the spatial template, but are not shown here. Aswe can see, although it is located 0.68 degrees away from theComa reference center, 4FGL J1256.9 + ff use tSZ template case. This is due to the combination of thelow S / N, the limited Fermi-LAT angular resolution and the ex-tension of the ICM emission model. Indeed, the Fermi-LAT PSFvaries with energy and around 1 GeV, where the signal is peak-ing, the 68% containment angle is about 1 deg (about 3 deg for
23. Adam et al.: γ -ray emission toward Coma Frequency (MHz)10 F ( J y ) CRe2CRe1, power-law modelData Energy (GeV)10 E d N d E dSd t ( M e V c m s ) IC (CRe2)IC (CRe1, power-law model)Hadronic emission
Figure A.1:
Estimation of the inverse Compton emission to the γ -ray signal. Left : Radio synchrotron spectrum data to which the CRe populationcan be matched.
Right : Associated emission in the γ -ray energy range. D e c . ( d e g ) Coma reference center and D e c . ( d e g ) Coma reference center and N u m b e r o f M C c o un t s N u m b e r o f M C c o un t s Figure B.1:
Top : sky distribution of the test source coordinates recovered for each Monte Carlo realization and comparison to the Coma referencecenter and the 4FGL J1256.9 + Bottom : distribution of the o ff set from the Coma reference center of the Monte Carlo realizations.The left panels correspond to the extended ICM model and the right panel to the tSZ template model.24. Adam et al.: γ -ray emission toward Coma
95% containment). This can be compared to the o ff set of 0.68deg between 4FGL J1256.9 + / N and the fact the ICM is an extended source (e.g.,surface brightness drops by a factor of ∼ θ = . ff set beingnot significant.Although a di ff use emission plus point source model pro-vides the best match to the data (see Section 4.4), these resultsconfirm that the scenario 2 agrees with the data. This also agreeswith the results presented in the main paper, in which no sig-nificant residual is observed around Coma once the best-fit ICMmodel is subtracted from the data, and where the comparison ofthe TS given in Table 1 does not point toward the need of includ-ing 4FGL J1256.9 + ff useICM component. See also Section 4.4 for discussions. Appendix C: ROI modeling and fitting: A null test inthe direction of point sources
In this Appendix, we present a null test, which we use to val-idate our ROI modeling and fitting procedure. As discussed inthe main text, the source 4FGL J1256.9 + + ff use ICM emission (scenario 2). Here, we checkthat this is not a general feature of our procedure and that it doesnot apply to other sources.We select the following sources from our field because theirlocal background is very close to the Coma region (no closeby source in the 4FGL catalog). However, we note that theirspectra are generally steeper and their TS may di ff er: 4FGLJ1253.8 + = .
57, 2 degrees northwest from the Comacenter, power law index of 1.9), 4FGL J1250.8 + = .
29, 3.8 degrees northwest from the Coma center, power lawindex of 2.16), 4FGL J1316.5 + = .
90, 4.3 degreesnortheast from Coma, power law index of 2.1). To match what isdone for Coma, we define the cluster center 0.68 degrees awayfrom these sources and test four directions for the o ff set (north,east, south, west) to increase statistics. We also perform the testwith the cluster centered on the source, but we note that thisdoes not match the case of Coma. Because the spectra of thesesources are not the same as 4FGL J1256.9 + ff erent TS values that we obtain for the dif-ferent test cases in Table C.1. We can see that for all tests, theTS is much smaller for the cluster than the point source whenboth are included, even when the two are co-aligned. In the caseswhere the point source is replaced by the cluster, the TS of thecluster increases, but remains much smaller than for the pointsource except when the cluster is centered on the source (al-though it always remains lower than for the point source only).We conclude that while the point source4FGL J1256.9 + InitialInjection ContinuousInjection
Figure D.1:
Same as Figure 13 for the
InitialInjection model(left) and the
ContinuousInjection model (right). the other sources 4FGL J1253.8 + + + Appendix D: Constraints on the cosmic rayelectron populations with alternative models
This appendix provides the constraints on the CRe populationsin the case of the alternative models considered in the case ofthe spectral modeling of the CRe in Figure D.1. Similarly,Figure D.2 gives the results for the alternative spatial modelsfor the CRp population. The results are also provided in the caseof including 4FGL J1256.9 + Fermi -LAT sky model(i.e., scenario 3).
25. Adam et al.: γ -ray emission toward Coma Table C.1:
Recovered TS for the di ff erent null tests and comparison to the initial TS of the considered sources.Test case Initial single point source Point source + cluster Single clusterTS point source (TS point source , TS cluster ) TS cluster + + + + + + + + + + + + + + + γ -ray emission toward Coma n CRp ∝ n gas n CRp = constant n CRp ∝ P gas n CRp ∝ n / with 4FGL J1256.9 + + CRe interpretation (MHz)10 F ( J y ) CRe1CRe2TotaltSZ absolute valueWSRT 352 MHzData 10 (MHz)10 F ( J y ) CRe1CRe2TotaltSZ absolute valueWSRT 352 MHzData 10 (MHz)10 F ( J y ) CRe1CRe2TotaltSZ absolute valueWSRT 352 MHzData 10 (MHz)10 F ( J y ) CRe1CRe2TotaltSZ absolute valueWSRT 352 MHzData10 Radius (kpc)10 S ( J y a r c m i n ) CRe1CRe2TotalData 10 Radius (kpc)10 S ( J y a r c m i n ) CRe1CRe2TotalData 10 Radius (kpc)10 S ( J y a r c m i n ) CRe1CRe2TotalData 10 Radius (kpc)10 S ( J y a r c m i n ) CRe1CRe2TotalData10 Radius (kpc)10 C R d e n s i t y ( c m G e V ) CRe1, 1 GeVCRe2, 1 GeVCRe1, 10 GeVCRe2, 10 GeV 10 Radius (kpc)10 C R d e n s i t y ( c m G e V ) CRe1, 1 GeVCRe2, 1 GeVCRe1, 10 GeVCRe2, 10 GeV 10 Radius (kpc)10 C R d e n s i t y ( c m G e V ) CRe1, 1 GeVCRe2, 1 GeVCRe1, 10 GeVCRe2, 10 GeV 10 Radius (kpc)10 C R d e n s i t y ( c m G e V ) CRe1, 1 GeVCRe2, 1 GeVCRe1, 10 GeVCRe2, 10 GeV10 Radius (kpc)10 C R e / C R e Radius (kpc)10 C R e / C R e Radius (kpc)10 C R e / C R e Radius (kpc)10 C R e / C R e CRe + reacceleration interpretation Figure D.2:
Same as Figure 12 (rows 1, 2, 3, and 4) and Figure 13 (rows 6 and 7) for the compact model (first column), the flat model (secondcolumn), isobaric (third column) and the extended model ++