Gamma-ray detection toward the Coma cluster with Fermi-LAT: Implications for the cosmic ray content in the hadronic scenario

The presence of relativistic electrons within the diffuse gas phase of galaxy clusters is now well established, but their detailed origin remains unclear. Cosmic ray protons are also expected to accumulate during the formation of clusters and would lead to gamma-ray emission through hadronic interactions within the thermal gas. Recently, the detection of gamma-ray emission has been reported toward the Coma cluster with Fermi-LAT. Assuming that this gamma-ray emission arises from hadronic interactions in the ICM, we aim at exploring the implication of this signal on the cosmic ray content of the Coma cluster. We use the MINOT software to build a physical model of the cluster and apply it to the Fermi-LAT data. We also consider contamination from compact sources and the impact of various systematic effects. We confirm that a significant gamma-ray signal is observed within the characteristic radius $\theta_{500}$ of the Coma cluster, with a test statistic TS~27 for our baseline model. The presence of a possible point source may account for most of the observed signal. However, this source could also correspond to the peak of the diffuse emission of the cluster itself and extended models match the data better. We constrain the cosmic ray to thermal energy ratio within $R_{500}$ to $X_{\rm CRp}=1.79^{+1.11}_{-0.30}$\% and the slope of the energy spectrum of cosmic rays to $\alpha=2.80^{+0.67}_{-0.13}$. Finally, we compute the synchrotron emission associated with the secondary electrons produced in hadronic interactions assuming steady state. This emission is about four times lower than the overall observed radio signal, so that primary cosmic ray electrons or reacceleration of secondary electrons is necessary to explain the total emission. Assuming an hadronic origin of the signal, our results provide the first quantitative measurement of the cosmic ray proton content in a cluster.[Abridged]


Introduction
Galaxy clusters form hierarchically via the smooth accretion of surrounding material and the merging of subclusters and groups (Kravtsov & Borgani 2012). During these processes, most of the binding gravitational energy is dissipated into the hot, thermal, ionized gas phase. However, shock waves propagating in the intracluster medium (ICM) and turbulence are also expected to accelerate both electrons and protons at relativistic energies, leading to a nonthermal population of cosmic rays (CR) that are confined within the cluster magnetic fields (see Brunetti & Jones 2014, for a review). In addition, CR may also arise from direct injection in the ICM, trough active galactic nucleus (AGN) outbursts (e.g., Bonafede et al. 2014a), or galactic winds asso-Corresponding author: Rémi Adam,remi.adam@llr.in2p3.fr ciated with star formation activity in cluster member galaxies (e.g., Rephaeli & Sadeh 2016). While high energy cosmic ray electrons (CRe) should quickly lose their energy (∼ 10 8 years at 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 remain largely unknown to date.
In fact, the presence of CRe and magnetic fields in galaxy clusters is now well established, thanks to the deep observations of diffuse radio synchrotron emission obtained over the last decade (see van Weeren et al. 2019, for a review). We can distinguish between radio relics, believed to trace merger shock acceleration at the periphery of clusters (e.g., van Weeren et al. 2010), and radio halos, which are spatially coincident with the thermal gas in massive clusters. Radio halos are generally fur-ther divided into mini-halos, associated with relaxed cool-core clusters (Giacintucci et al. 2017), and giant radio halos, which are megaparsec-size sources found in merging clusters (Cassano et 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 to explain the origin of CRe that generate radio halos: the hadronic model (Dennison 1980;Blasi & Colafrancesco 1999;Dolag & Enßlin 2000), in which secondary CRe are the products of pion decay generated in collisions between the CRp and the thermal gas, and the turbulent reacceleration of a seed population of electrons (possibly secondary CRe; Brunetti & Lazarian 2007. In both scenarios, CR might show up as a γ-ray signal due to the inverse Compton interaction of CRe with background light. In the hadronic model, the decays of pions should also lead to additional γ-ray (and neutrino) emission, which is expected to be the dominant component of the γ-ray flux at energies above ∼ 100 MeV (Pinzke & Pfrommer 2010).
Attempts to detect the γ-ray emission from galaxy clusters have been carried out over the last two decades using individual targets, stacking analysis, and cross correlations to external datasets (see, e.g., Reimer et al. 2003;Ackermann et al. 2010;Aleksić et al. 2010;Dutson et al. 2013;Huber et al. 2013;Griffin et al. 2014;Branchini et al. 2017;Colavincenzo et al. 2020). Despite the lack of detection, these searches proved extremely useful in constraining the CRp content of clusters to below a few percent relative to the thermal pressure (Ackermann et al. 2014).
Among the relevant individual targets, the Coma cluster has been one of the most promising sources to search for γ-ray emission. Indeed, it is a massive (M 500 7 × 10 14 M ) and nearby system (redshift z = 0.023, about 100 Mpc); as such its γ-ray flux is expected to be large (Pinzke & Pfrommer 2010). Additionally, it is located near the galactic north pole and thus benefits from a low galactic background. The signal is expected to be extended even for the Fermi-LAT (characteristic radius θ 500 ∼ 0.75 deg), but not so extended that the analysis suffer strong systematic effects in the background modeling. The region around Coma is not affected by the presence of bright γray compact sources. Finally, the Coma cluster is a well-known ongoing merger, which has been extensively observed at various 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 relic have been observed, proving the presence of CRe (e.g., Willson 1970; Giovannini et al. 1993;Thierbach et al. 2003;Brown & Rudnick 2011;Bonafede et al. 2020).
Several analyses have focussed on the search for γ-rays toward the Coma cluster, using both ground-based and spacebased instruments (e.g., Perkins 2008;Aharonian et al. 2009;Arlen et al. 2012;Prokhorov 2014;Zandanel & Ando 2014). It is noteworthy that the Coma cluster was analyzed by the Fermi-LAT collaboration (Ackermann et al. 2016), with six years of data, who found some residual emission in the direction of the cluster, though not enough to claim a detection. Despite these unsuccessful searches, the combination of the obtained γ-ray upper limits with radio synchrotron data was used to show that pure hadronic models cannot explain all the observed radio emission in the case of Coma (Brunetti et al. 2012), and this allowed constraints to be set on turbulent reacceleration models . Recently, Xi et al. (2018) claimed the first significant 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 for the CR physics were carried out. In 2020, the Fermi-LAT collaboration released an update of the 4FGL catalog (4FGL-DR2, Abdollahi et al. 2020;Ballet et al. 2020), also indicating that a source was detected in the direction of the Coma cluster, named 4FGL J1256.9+2736. While this source could be a contaminant, it could also be associated with the cluster's diffuse emission itself.
In this paper, we explore the consequences on the CR physics of the γ-ray emission observed in the direction of the Coma cluster under the assumption that this signal is associated with the diffuse ICM. To do so, we present an analysis of the Fermi-LAT data (Atwood et al. 2009), to search for γ-ray emission within the Coma cluster region, and use this measurement to constrain the properties of the CR content of the cluster. We construct a model for the expected signal, in which we set the thermal gas properties using thermal Sunyaev-Zel'dovich (tSZ) effect and X-ray observations. The Fermi-LAT data are analyzed, and we extract the spectral energy distribution (SED) of the emission assuming different spatial templates. The Fermi-LAT maps are compared to data at other wavelengths for better interpretation of the results. In particular, we compare the γ-ray signal to the radio synchrotron emission, to the thermal pressure traced by the tSZ signal, to the galaxy distribution, and to the thermal gas distribution from X-ray data. We use our model to constrain the CRp population in the Coma cluster assuming that all the signal observed is due to the hadronic interactions between CRp and the ambient thermal gas, but also consider possible point source contamination. Finally, we compute the amount of secondary CRe associated with the γ-ray emission assuming steady state. We use this model to constrain the fraction of primary to secondary CRe in the absence of reacceleration, and the boost required from reacceleration with respect to the steady state hadronic model, to explain the radio emission assuming that reacceleration of secondary electrons is dominant, respectively. This paper is organized as follows. Section 2 describes the multiwavelength dataset that is used through the paper. In Section 3, we present the physical modeling of the cluster and the computation of the observables. The Fermi-LAT analysis is presented in Section 4 and the outputs are compared to multiwavelength data in Section 5. Section 6 discusses the implications of the observed signal for the CR content of the Coma cluster. In Section 7, we discuss the results presented in this paper and compare them to the literature. Finally, a summary and conclusions are given in Section 8. Throughout this paper, we assume a flat ΛCDM cosmology according to Planck results (Planck Collaboration et al. 2016b) with H 0 = 67.8 km s −1 Mpc −1 , Ω M = 0.308, and Ω Λ = 0.692. The reference coordinates of the Coma cluster are RA, Dec = 194.9118,27.9537 degrees and the redshift z = 0.0231, based on the Planck catalog (PSZ2, Planck Collaboration et al. 2016a). Given the reference cosmological model, 1 degree corresponds to 1.73 Mpc at the redshift of Coma. We use the value of R 500 = 1310 kpc based on Planck Collaboration et al. (2013), corresponding to M 500 = 6.13 × 10 14 M given our cosmological model.

Data
While this paper focuses mainly on the γ-ray analysis of the Coma cluster, we also use complementary data both for multiwavelength comparison of the signal and in order to build templates used to model the γ-ray emission. In addition to the Fermi-LAT γ-ray data, we thus use thermal Sunyaev-Zel'dovich data from Planck, X-ray data from ROSAT, radio data from the Westerbork Synthesis Radio Telescope and optical data from the Sloan Digital Sky Survey. This section presents a brief description of these data.

Fermi Large Area Telescope γ-ray data
The γ-ray analysis is based on the publicly available 1 Fermi Large Area Telescope data (Fermi-LAT, Atwood et al. 2009), which has been collecting all sky γ-ray data at energies from about 20 MeV to more than 300 GeV since June 2009. In this paper, we use almost 12 years of Pass 8 data (P8R3), collected from August 4, 2008, to April 2, 2020. Events within a radius of 10 degrees from the cluster center were collected for the analysis (see also Section 4 for the data selection). As a description of the compact sources around the cluster, we use the second release of the 4FGL catalog (see Abdollahi et al. 2020 for the 4FGL catalog and Ballet et al. 2020 for the 4FGL-DR2 update). The 4FGL-DR2 catalog is an incremental version of the 4FGL catalog. It is based on 10 years of survey data in the 50 MeV -1 TeV range (see Abdollahi et al. 2020, for the analysis improvement relative to previous catalogs). The diffuse isotropic background and the galactic interstellar emission are modeled using the latest available templates: iso P8R3 SOURCE V2 v1 and gll iem v07, respectively.

Planck thermal Sunyav-Zel'dovich data
The tSZ effect (Sunyaev & Zeldovich 1972) is a relevant probe for comparison to the γ-ray signal as it provides a direct measurement of the integrated thermal gas pressure along the lineof-sight (see Birkinshaw 1999;Mroczkowski et al. 2019, for reviews). We used the MILCA (Hurier et al. 2013) Compton parameter map obtained from Planck (Planck Collaboration et al. 2016c) to extract a high signal-to-noise image of the signal at an angular resolution of 10 arcmin (full width at half maximum, FWHM). In the region of Coma, the signal is not significantly affected by any artifacts, e.g., from compact sources, or remaining diffuse galactic emission (see also Planck Collaboration et al. 2013, for a Planck analysis of the Coma cluster). The data are publicly available at the Planck Legacy Archive 2 .

ROSAT All Sky Survey X-ray data
The X-ray diffuse cluster emission traces the thermal gas density (see Böhringer & Werner 2010, for a review) and is expected to be correlated with the γ-ray signal. Given the large size of the region of interest (ROI) and the fact that the angular resolution is not critical for our purpose, we used ROSAT (Truemper 1993) data as obtained from the ROSAT All Sky Survey (RASS) 3 to image the cluster at a resolution of 1.8 arcmin (FWHM). To do so, we subtracted the background and normalized the data using the exposure maps. In addition, we subtracted the compact sources present in the field using both the RASS bright and faint source catalogs (Voges et al. 1999(Voges et al. , 2000.

Westerbork Synthesis Radio Telescope data
The radio emission traces relativistic electrons via synchrotron emission and could be associated with the γ-ray signal. We used the Westerbork Synthesis Radio Telescope (WSRT) radio data at 352 MHz (Brown & Rudnick 2011) in order to dispose of an image of the cluster halo and relic. Compact sources were subtracted using the method described in Rudnick (2002). The residual contamination was estimated by injecting fake point sources into the image, applying the filter to the image, and measuring the residual flux. For sources that were detected by more than 2σ, the residual flux was less than 3% and as such the contamination is expected to be less than 3%. Nevertheless, some of the emission from the bright central radio galaxy remains blended with the signal, and we masked it in the following analysis using a threshold of 100 mJy/beam. Additionally, the bright source Coma A north of the radio relic causes some artifacts that should be accounted for in the morphological comparison (Section 5). The image angular resolution is 4 × 3 arcmin 2 (FWHM).
The WSRT map is also used to extract the average radial profile of the radio halo. To do so, we first subtract the zero level (i.e., offset), and define bins of 3 arcmin width. The zero level arises from the standard clean deconvolution that creates a negative offset around bright sources such as those found in the Coma cluster. After filtering out the point sources and convolving the image to a lower resolution, this resulted in a local background of -20 mJy/beam. Error bars on the profile are computed using inverse variance Gaussian weighting with a constant noise root-mean-square of 5 mJy/beam. The pixels from the masked regions are excluded from the profile. We also convert from Jy/beam to Jy/sr assuming that the beam is Gaussian. We stress that we are considering the total (i.e., average) profile, but the variations in profile from one direction to another (Brown & Rudnick 2011) should not have a significant effect on our analysis given the purpose of this paper. The profile will be used in Section 6 as a comparison to our model.
In addition, we also use the integrated fluxes of the Coma radio halo as compiled in Pizzo (2010). Given the diffuse nature of the signal, it is important to make sure that the flux we used was measured in a consistent aperture for the different observations. In this respect, we used a radius of 0.48 × R 500 (i.e., 629 kpc or 0.36 deg, see Brunetti et al. 2013, for details).

Sloan Digital Sky Survey optical data
The spatial distribution of galaxies is informative of the different groups in the cluster, and the γ-ray emission might be associated with CR escaping from cluster member galaxies. We thus use Sloan Digital Sky Survey (SDSS, York et al. 2000) optical data 4 to construct a color image of the Coma region. In addition, we select galaxies with spectroscopic information from the SDSS database 5 and use this catalog to construct a galaxy density map of the Coma region. To to so, we bin the galaxy catalog on a grid, and we weight them with w = exp − (z gal c−z Coma c) 2 ) 2σ 2 v to minimize the background. The quantity c is the speed of light, z gal and z Coma the redshift of each galaxy and that of the cluster, respectively, and σ v = 2128 km s −1 (i.e., FWHM of 5000 km s −1 ). The map is convolved with a Gaussian kernel with a FWHM of 10 arcmin to enhance the signal-to-noise. The background is then estimated as the mean of the map at radii above 90 arcmin from the cluster, and subtracted from the map.

Modeling
The physical modeling of the signal, as expected in the γ-ray band at Fermi-LAT energies, is necessary for two reasons. First, we aim at constructing a template (spatial and spectral) of the diffuse emission in order to include the cluster when fitting the data (as done in Section 4). Then, the physical model is needed to connect the observations to the CR physics and will be used to constrain the CR content of the Coma cluster (Section 6). Additionally, the synchrotron emission should be included in the modeling when comparing the implication of the Fermi-LAT constraint on the radio signal.
In this paper, we model the Coma cluster using the MINOT software (Adam et al. 2020) 6 . This allows us to have a complete description 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 using multiwavelength data. While this approach may provide a better match to the spatial shape of the signal due to deviations from spherical symmetry, it does not allow us to constrain the CR content of the cluster because of sky projection and degeneracies between the different physical components.

Thermal component
The MINOT software models the thermal ICM component using the gas electron pressure and density profiles. The gas density profile is key to the γ-ray emission as it provides the target for proton-proton interactions. The thermal pressure is also particularly important, because it will be used to provide a normalization 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 with parameters taken from Planck Collaboration et al. (2013) and rescaled to our cosmological model (P 0 , r p , a, b, c = 0.022 keV cm −3 , 466.8 kpc, 1.8, 3.1, 0.0). The electron density is described by a β-model (Cavaliere & Fusco-Femiano 1978) according to with parameters taken from Briel et al. (1992) and rescaled to our cosmological model (n 0 , r c , β = 3.36 × 10 −3 cm −3 , 310 kpc, 0.75). From the electron pressure and electron density, we also derive the total gas pressure and thermal proton (and heavier elements) density assuming a constant helium mass fraction of Y He = 0.2735 and metal abundances of 0.3, relative to the solar value (see Adam et al. 2020, for more details). Finally, we note that the radial models are truncated at a radius r > 3 × R 500 . This also apply for nonthermal quantities discussed below. 6 https://github.com/remi-adam/minot

Magnetic field
The synchrotron emission, from a given CRe population, depends on the magnetic field. The later is thus an important component of our model when comparing the implications of our constraints with radio data (Section 6.3). The magnetic field strength profile of the Coma cluster was inferred by Bonafede et al. (2010) using Faraday rotation measures and we used their result to model it as where B 0 = 4.7 µG and η B = 0.5. We note, however, that uncertainties on these parameters are relatively large and the parameters are degenerate (see also the recent work by Johnson et al. 2020). They are constrained in the range ( B 0 = 3.9 µG; η B = 0.4) to ( B 0 = 5.4 µG; η B = 0.7) (see Bonafede et al. 2010, for more details). The impact of the magnetic field modeling will be further discussed in Section 7.2.

Cosmic ray protons
The CRp number density, per unit volume and energy, is modeled as where f 1 (E) is the energy spectrum, f 2 (r) the radial dependence, and A CRp is the normalization. The energy spectrum is modeled as a power law (expected from Fermi acceleration), as with α CRp the slope of the CRp energy spectrum. The radial dependence is modeled assuming that the CRp number density scales with either the thermal density or pressure, as Since the radial dependence of the CRp is not well-known, we test different values for the scaling, η CRp = 0, 1 2 , 1 , corresponding to flat, extended, and compact distributions. Finally, the normalization A CRp is computed relative to the thermal energy enclosed within the radius R 500 . In the following, we use the parameter where U CRp (R) is the total CRp energy enclosed within R computed by integrating Equation 4 and U th (R) the total thermal energy, obtained by integrating the gas pressure profile (see Adam et al. 2020, for more details).
In the end, three parameters are used to describe the CRp: the spatial scaling relative to the thermal gas η CRp , the CRp spectrum slope α CRp , and the normalization encoded in X CRp (R 500 ). The value of these parameters will be further discussed in Section 4 and 6.

Primary cosmic ray electrons
At the Fermi-LAT energies, the γ-ray emission is expected to be dominated by hadronic processes, and we neglect the inverse Compton emission due to secondary or primary CRe interacting with background light (see Appendix A for more discussions). However, CRe need to be included in our model when computing the radio synchrotron emission in Section 6. The primary cosmic ray electrons, CRe 1 , are modeled similarly as the CRp. However, CRe 1 are affected by losses and we therefore account for this using three different models: the ExponentialCutoffPowerLaw, the InitialInjection, and respectively. Using different parametrizations will allow us to estimate the systematic effects associated with the model.
There are therefore four parameters used to describe the primary CRe 1 : the spatial scaling of the CRe 1 population relative to the thermal gas, η CRe 1 (as in Equation 6), the CRe 1 spectrum slope α CRe 1 and break energy E cut,CRe 1 , and the normalization encoded in X CRe 1 (R 500 ) (as in Equation 7). The value of these parameters will be further discussed when comparing our model to radio data in Section 6.

γ-ray and synchrotron signal
Given the thermal gas, magnetic field and CRp modeling, we compute the γ-ray emission due to hadronic collisions, the secondary cosmic ray electrons assuming a steady state scenario (CRe 2 ), and the radio synchrotron emission (from both CRe 1 and CRe 2 ) according to Adam et al. (2020). We do not account for inverse Compton emission at Fermi-LAT energies because it is expected to be negligible, as detailed in Appendix A.
In brief, the production rate of γ-rays and CRe 2 are computed as where dN col dE CRp dVdt ∝ n e (r) × J CRp (r, E) is the CRp-ICM collision rate and the function F E, E CRp gives the number of secondary particles (electrons or γ-rays) produced in a collision as a function of the initial energy of the CRp. This computation is based on the parametrization by Kelner et al. (2006) and Kafexhiu et al. (2014) following the work by Zabalza (2015), as detailed in Adam et al. (2020). The distribution of CRe 2 is then obtained 5 R. Adam et al.: γ-ray emission toward Coma based on the distribution of galaxies, the tSZ signal, the X-ray emission, the radio halo emission, and the radio relic emission. Images are 3 × 3 deg 2 , except for the radio relic (bottom right), which is only 1 × 1 deg 2 . Units are arbitrary, but the integral of all maps over the solid angle is set to 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 integrating Equation 11 along the line-of-sight, accounting for the distance to the Coma cluster, and integrating over the energy or the solid angle, respectively. The γ-ray attenuation due to interaction with the extragalactic background light (EBL) is also included and induces a cutoff in the spectrum above E 10 4 GeV. These spatial (profile) and spectral templates are used in Section 4 to account for the cluster in the modeling of the ROI. The synchrotron emission rate is computed following the prescription given in Aharonian et al. (2010). This assumes that the orientation of the magnetic field is randomized. Observable profile and spectra are then obtained as for the γ-rays.
In Figure 1, we show how the different profiles and spectral slopes for the CRp translate into the γ-ray surface brightness profile and spectrum. The top left panel shows the four different radial distributions used for the CRp, in terms of n CRp (r) = J CRp (r, E) dE. As we can see, the difference between models based on the thermal density or thermal pressure are very close in the case of Coma as expected given the fact that the cluster is close to be isothermal. On the top right panel, we can see the correspondence in terms of the profile of the CRp energy to thermal energy ratio. We note that the normalization was fixed to X CRp (R 500 ) = 10 −2 . The bottom left panel shows the γray surface brightness for the different models. Because it arises from the product of the CRp and thermal gas density, the profiles are more peaked than the original CRp distributions. On the bottom right panel, we can see the energy spectrum integrated within R 500 for different slopes α CRp . The CRp slope mainly drives the γ-ray slope between 1 GeV and 10 4 GeV. At higher en-ergies, the extragalactic background light induces a cutoff in the spectrum, but this is not in the range accessible by Fermi-LAT. At lower energies, the energy threshold of the proton-proton collisions implies that the spectrum smoothly vanishes. Since the normalization is fixed, steeper is the spectrum and higher is the flux near the peak. As seen with the dash-dotted lines, the amplitude decreases when the CRp profiles become flatter because it leads to a decrease in the particle collision rate within R 500 .

Construction of spatial templates from multiwavelength data
In addition to the spherically symmetric physical models, we build spatial templates based on the multiwavelength data discussed in Section 2. We consider templates based on the galaxy density, the tSZ Compton parameter, the X-ray surface brightness and the radio surface brightness (halo and relic). The maps are projected on 5 × 5 deg 2 maps with 1 arcmin pixel size (well below 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 below the Fermi-LAT angular resolution so that it will be negligible when 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 the reference center for the galaxy density, tSZ and X-ray maps. For the radio map, two cases are considered: the halo for which pixels more than 60 arcmin from the center are masked (to avoid including the relic in the template) and the relic for which we 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 exclude any bright radio galaxy focusing on the emission from the relic only.
The resulting templates are shown in Figure 2 on 3 × 3 deg 2 grids (except for the relic, for which it is 1 × 1 deg 2 ). In all cases, the cluster is elongated from the northeast to the southwest, with an excess associated with the NGC 4839 group on the southwest. The most compact template is the one based on the X-ray image as it is proportional to the gas density squared. The tSZ template, proportional to the gas pressure (temperature times density) is much more extended. The galaxy density template extension is in between the tSZ and X-ray ones, but it is more elongated. The radio halo template, on the other end, is very spherical and matches well the tSZ template in terms of the extension. The radio relic template is very elongated from the southeast to the northwest, much smaller than the other templates, and not spatially coincident (as will be further discussed in Section 5).

Fermi-LAT analysis
In this Section, we describe the Fermi-LAT γ-ray analysis, performed using Fermipy (Wood et al. 2017). After the data selection, we model and fit the ROI in order to extract the signal in the Coma cluster region under different modeling assumptions (signal associated with the ICM, a point source, the combination of both) and compare these scenarios. The cluster model built in Section 3 is used to extract the SED of the source. We also perform several tests to validate the global ROI model as a function of energy, radius and check the time stability. Finally, we discuss the systematic effects that might affect our findings.

Data selection and binning
Following the Fermi-LAT collaboration recommendations for off galactic plane compact source analysis, we apply the P8R3 SOURCE V2 selection (event class 128). The energy dispersion is also accounted for. We select FRONT+BACK (event type 3) converting photons and apply a cut on zenith angles less than 90 degrees to effectively remove photons originating from the Earth limb. 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 2 around the Coma center. As a baseline, the data were binned using a pixel resolution of 0.1 deg and 8 energy bins per decade between 200 MeV and 300 GeV. The choice of the low energy threshold is a compromise between count statistics and robustness of the results since systematic effects increase at low energy (see also Section 4.7).

ROI modeling
We start by modeling the ROI accounting for the diffuse background components, corresponding to the isotropic and galactic interstellar emission, as well as the 4FGL compact sources. Given the size of our ROI, the most distant pixels are 8.5 degrees away from the center. Nevertheless, we include all sources within a square with a width of 20 degrees from the reference center, because the point spread function (PSF) may extend up to about 10 degrees at low energies.
As discussed in Section 1, the 4FGL-DR2 catalog includes a source named 4FGL J1256.9+2736 that lies close to the Coma cluster peak. It is located at (RA,Dec) = 194.2417, 27.6076, corresponding to about 0.68 degrees from the cluster center and within θ 500 . This source is detected with a significance of 4.2, and is modeled by a point source with power law spectrum with index 2.73 in the 4FGL-DR2 catalog. Its origin is uncertain, but it is given as possibly associated with NGC 4839 (with a probability of 0.24). While this source could be a contaminant for the diffuse ICM signal, as, e.g., the AGN of NGC 4839, it could also corresponds to the peak of the diffuse emission associated with the ICM. This second scenario is motivated by the uncertain origin of the source due to the limited signal to noise ratio and angular resolution of the Fermi-LAT. We thus aim at testing and comparing the two hypothesis. To do so, we consider the three following scenarios for modeling the Coma cluster region. 1) We assume that 4FGL J1256.9+2736 is a point source independent of the cluster diffuse emission. In this case, we include it in the ROI model as given by the 4FGL-DR2 catalog and do not include any extra diffuse emission associated with the ICM. 2) We assume that 4FGL J1256.9+2736 is in fact associated with the diffuse cluster emission. In this case, the 4FGL-DR2 catalog model is not an appropriate description of the signal and we replace 4FGL J1256.9+2736 by a diffuse emission model associated with the ICM. 3) We assume that 4FGL J1256.9+2736 is a point source independent of the cluster diffuse emission, but we also consider a diffuse ICM component as γ-ray emission is expected from the cluster. In this case, we include 4FGL J1256.9+2736 in the ROI model as given by the 4FGL-DR2 catalog and add an extra component to account for the diffuse emission associated with the ICM.
The diffuse Coma cluster ICM emission is modeled either using the MINOT spatial and spectral physical templates, or using the MINOT spectral template together with the spatial templates built from other wavelengths (see Section 3). We test different shapes for the profile of the CRp, and fix the CRp spectral index slope to α CRp = 2.8. This value will be later constrained using the measured SED, in Section 6. The nonoptimal choice of this number will slightly reduce the significance associated with the Coma diffuse emission, but does not affect the SED fit performed afterward. As a baseline, and unless otherwise stated, we use the model with η CRp = 1/2 scaled with respect the gas thermal density (extended model).

ROI fitting
Once the overall background model (large-scale diffuse plus point source emission) and cluster diffuse emission is defined, we use the following iterative procedure to fit the ROI. 1) We construct a first sky model using the diffuse background and the point sources spectral parametrization, with parameters from the 4FGL catalog, but exclude the cluster diffuse emission model at this stage. 2) We run the optimize function of Fermipy, including all sources in the model (except for the cluster diffuse emission). We thus obtain a first renormalization of the sources in the model and a first value of their test statistics (TS, Mattox et al. 1996), defined as with L 0 the maximum likelihood value for the null hypothesis, and L the maximum likelihood when including the additional source. 3) To allow for variability, we free all the model parameters associated with sources with √ TS > 20, and free only the normalization of sources with 5 < √ TS < 20. 4) We use the fit function of Fermipy to perform the maximum likelihood 7 R. Adam et al.: γ-ray emission toward Coma 4 F G L J 1 2 2 4 . 9 + 2 1 2 2 ( n o r m ) 4 F G L J 1 2 3 0 . 2 + 2 5 1 7 ( n o r m ) 4 F G L J 1 2 3 1 . 7 + 2 8 4 7 ( n o r m ) 4 F G L J 1 2 3 1 . 7 + 2 8 4 7 ( P L i n d e x ) 4 F G L J 1 2 3 7 . 0 + 3 0 1 9 ( n o r m ) 4 F G L J 1 2 5 0 . 8 + 3 1 1 7 ( n o r m ) 4 F G L J 1 2 5 3 . 8 + 2 9 2 9 ( n o r m ) 4 F G L J 1 2 5 4 . 5 + 2 2 1 0 ( n o r m ) 4 F G L J 1 2 5 6 . 9 + 2 7 3 6 ( n o r m ) 4 F G L J 1 2 5 7 . 8 + 3 2 2 8 ( n o r m ) 4 F G L J 1 2 5 7 . 8 + 3 2 2 8 ( L P 4 F G L J 1 3 2 0 . 7 + 3 3 1 4 ( n o r m ) 4 F G L J 1 3 2 1 . 1 + 2 2 1 6 ( n o r m ) 4 F G L J 1 3 2 1 . 1 + 2 2 1 6 ( L P ) 4 F G L J 1 3 2 1 . 1 + 2 2 1 6 ( L P ) 4 F G L J 1 3 2 1 . 9 + 3 2 1 9 ( n o r m ) 4 F G L J 1 3 2 3 . 0 + 2 9 4 1 ( n o r m ) 4 F G L J 1 3 2 3 . 0 + 2 9 4 1 ( P L i n d e x ) 4 F G L J 1 3 2 6 . 9 + 2 2 1 0 ( n o r m ) galdiff correspond to the isotropic and galactic diffuse backgrounds, Coma corresponds to the ICM cluster model, and all the other names to the 4FGL sources. The labels of the parameters are as follows: Norm stands for the spectrum normalization; PL index stands for the spectral index of sources 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+2736, and for the MINOT cluster model with n CRp ∝ n 1/2 e . model fit of all the free parameters in the sky model. 5) We add the cluster diffuse emission (except in the scenario 1, see subsection 4.2) in our model, with its normalization being its only free parameter, and refit the sky model as in step 4. In addition to the likelihood fit, we also use the tsmap Fermipy function to compute TS maps, in the case of radially symmetric cluster models.
In Figure 3, we present the correlation matrix recovered from the likelihood fit for all the fitted parameters included in our model (i.e., the free parameters). It corresponds to the case of the scenario 3 (see subsection 4.2), i.e., including both the Coma cluster diffuse emission and the contribution from 4FGL J1256.9+2736. As we can observe the diffuse cluster emission is strongly degenerate (anti-correlated) with the normalization of 4FGL J1256.9+2736, which is expected given its coordinate, as the signal from one component can be absorbed in the other one. This shows that a scenario in which 4FGL J1256.9+2736 is in fact associated with the cluster diffuse emission should be considered, as done in the following. The Coma cluster diffuse emission is, on the other hand, not significantly correlated with any other compact source in the model, but is marginally anti-correlated with the isotropic background. We conclude that apart from the presence or not of 4FGL J1256.9+2736, the Coma cluster model fit is expected to be robust with respect to mis-modeling of other components of the ROI (see also Section 4.7 for systematic effects associated with the diffuse background).
We present in Table 1 the TS values obtained for the point source 4FGL J1256.9+2736 and the diffuse cluster model component individually, for all the cases considered (we also include a point source model for the cluster as a spatial template for comparison). Because of the degeneracy between 4FGL J1256.9+2736 and the diffuse cluster model, we also include the TS value for the total (i.e., both components included or excluded when computing the likelihood L and L 0 , see Equation 12). Indeed, in the case of two degenerate sources, the maximum likelihood will not change much when removing one of the two source since the other source will absorb the missing component, in the null hypothesis. Nonetheless, removing the two sources simultaneously may lead to a large change in the maximum likelihood because at least one component is necessary to explain the data. Hence, the TS of each individual degenerate source is expected to be low, but that of the two sources taken together might be large. This is expected to be the case in scenario 3. All the considered MINOT radial models give similar TS values, around TS = 24 − 27, when replacing 4FGL J1256.9+2736 (scenario 2). This is similar to the TS value for 4FGL J1256.9+2736 alone, which is 25.61 (scenario 1). On the other hand, the cluster emission modeled as a point source gives a significantly lower value, of 17.54. Except for the radio relic, all the templates based on other data provide a better match to the signal, with the best match reached for the galaxy density and tSZ maps (TS ∼ 32 − 34). This is likely because of the elongation of the Coma cluster, especially toward the southwest, where most of the γ-ray excess is observed. In the case both 4FGL J1256.9+2736 and the diffuse cluster emission are included in the ROI model (scenario 3), the TS of each component drastically reduces, highlighting the fact that 4FGL J1256.9+2736 and the diffuse emission are degenerate. For instance, the TS value reaches only 12.17 and 12.44 for the cluster and 4FGL J1256.9+2736, respectively, in the case of the baseline model. However, the total TS value is higher compared to scenario 1 and 2, reaching about 32-36, but the improvement is only marginal considering the fact that two components, instead of one, are included. Again, this highlight the degeneracy between 4FGL J1256.9+2736 and the expected cluster signal. We conclude that while a diffuse cluster model (scenario 2) generally provides a better description of the data compared to 4FGL J1256.9+2736 alone (scenario 1), the two scenarios cannot be significantly discriminate based on their likelihood fit. In Appendix B, we also further discuss the agreement between scenario 2 and the data using a Monte Carlo realization.
Finally, we also consider the case of using a power law for the photon spectrum of Coma. In the case of our baseline spatial model (n CRp ∝ n 1/2 e ), in scenario 2, we obtain a best-fit photon index of 2.45 ± 0.19 for a TS value of 27.37. As expected, the overall spectrum is slightly flatter than in the case of our physical model. Indeed, the photon spectrum vanishes at low energy for a given CRp energy slope with our physical model, thus leading to a harder spectrum in this regime. The power law model thus averages the two regimes. Nevertheless, we note that the high energy photon spectrum does not strictly match the CRp energy spectrum (see Adam et al. 2020, for more details). The TS value is nearly the same as in the case of the physical spectrum and given the available signal-to-noise ratio, it is not possible to discriminate the two.
The modeling and fitting procedure described in Section 4.2 and in this Section is also validated using a null test described in Appendix C. This shows that no cluster detection is obtained when including a cluster in the sky model close to other sources with similar background as around Coma.

Comparison between data and model
In order to check how the data compare to the different models (in the different scenarios), we computed maps of the data, model and residual (with and without Coma; with and without 4FGL J1256.9+2736) in three energy bins: from 200 MeV to 300 GeV (the total considered range), from 200 MeV to 1 GeV, and from 1 GeV to 300 GeV. This is shown in Figure 4 in the case of the baseline model.
When excluding both the diffuse cluster emission and 4FGL J1256.9+2736, we observe a significant excess near the center of the map (where the green cross indicates the Coma reference center). This excess is also visible independently in the high and low energy bands, although with lower significance (TS > 16 and TS > 9, respectively). The peak of the excess is slightly offset in the southwest direction with respect to our reference center. The position angle of this elongation agrees in all energy bins, although the elongation itself is more pronounced at low energy. When including 4FGL J1256.9+2736 in the sky model, the central excess disappears almost entirely (especially since 4FGL J1256.9+2736 can absorb most of the sig-nal from Coma, if any, as indicated by the correlation matrix and discussed above). The comparison of the Fermi-LAT excess to other data will be discussed further in Section 5. We also show in Appendix B that the offset between 4FGL J1256.9+2736 and the Coma reference center agrees with 4FGL J1256.9+2736 corresponding to the peak of the diffuse ICM emission (i.e., scenario 2).
In addition to the central excess, we also see two other excesses near RA,Dec = 200, 25.5 degrees and RA,Dec = 202, 29.5 degrees (TS > 9). Another excess count is seen near RA,Dec = 197.5, 32 degrees (in the high energy bin), but its TS remains low because our cluster model does not match the spatial and/or spectral shape of this signal well. All these other excesses remain relatively small and given their significance, it is not clear whether they correspond to noise fluctuations, the mis-modeling of existing sources, or new sources that are not included in the 4FGL catalog. Given their distances to the central excess, we do not expect that they would lead to any significance bias in our analysis.
We also compute the counts as a function of energy observed within 3 × θ 500 from the cluster center and compare it to the different components of the model in Figure 5. As we can see, the dominant components are from the isotropic and diffuse galactic interstellar emission at almost all energies. The signal from the Coma cluster (in the case of scenario 2) is about an order of magnitude below depending on the energy, but is the dominant compact source except at very low energies where the larger PSF leads to leakage from strong sources within the Coma central region (we note that the source 4FGL J1253.8+2929, northwest of Coma, falls within 3 × θ 500 ). Similarly in the case of scenario 1, 4FGL J1256.9+2736 is the dominant compact source at low energy but not at high energy as its spectrum is steeper. Including the Coma cluster or 4FGL J1256.9+2736 in the model improves the residual, as expected from the result of the likelihood fit, but no significant difference is observed between the two scenarios accounting for the error bars. Given this spectrum, it is very unlikely that other compact sources from the ROI induce the observed excess within 3 × θ 500 of the Coma cluster.
In order to address to which extent the Fermi-LAT data are sensitive to the shape of our radial model, we compute the profile of the data and model in Figure 6. Although each of the 5 first bins (up to 1 degrees extent) are marginally inconsistent with the model when the Coma cluster emission or 4FGL J1256.9+2736 are not accounted for, they all point to an excess emission at a level of 0.5-3.5 σ, and thus correspond to a clear overall excess. The baseline cluster model provides a good fit to the data and significantly improves the residual (scenario 2). The point source cluster model, also shown for comparison, is too peaked with respect to the observations. However, it is not in clear disagreement with the data given the error bars. This is in agreement with the results of Table 1, where the likelihood ratio between the two model is ∆TS = 9.45, providing only a hint that the extended model is more appropriate than the point source model. In the case of scenario 1 (4FGL J1256.9+2736 only), the agreement with the data is also good.

Extraction of the cluster spectral energy distribution
Having a model for the sky in hand, we used the sed function of Fermipy to extract the SED of the Coma diffuse emission in the different cases tested in this paper. The sed function fits for the flux normalization of a source in independent bins of energy. To do so, we allowed the local photon slope to vary according to the MINOT global spectral model and we fixed the background   provides the full likelihood scan for the normalization value in each 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 diffuse emission. The measured values are listed in Table 1 for the different models and scenario which we test. We note that in the case of scenario 1 (i.e., no cluster in the model), the flux is by definition equal to zero. The flux ranges from about 13 to 19 × 10 −10 ph s −1 cm −2 for the radially symmetric models in scenario 2. The multiwavelength templates lead to similar fluxes, except for the radio relic template for which the flux is only about 9 × 10 −10 ph s −1 cm −2 . In case of scenario 3, part of the flux is absorbed in the 4FGL J1256.9+2736 component. This leads to a flux lower by a factor of about two, depending on the cluster model considered.

Light curve
In order to check that the observed signal is consistent with diffuse ICM emission, we compute the light curve of the signal. Indeed, while many γ-ray sources are variable, the γ-ray emission from galaxy clusters is expected to be steady, at least over the timescale of any human observation, and observing a burst would allow us to rule out the cluster ICM origin of the signal. We thus compute the light curve of the Coma cluster diffuse signal using the lightcurve function from Fermipy, in the case of scenario 2. We extract the flux within 30 bins, which corresponds to about 4.5 month per bin and report the result in Figure 7. Given the significance of the signal for the full dataset, we only expect upper limits in each bin. As we can see, no significant burst is observed. The upper limit excludes the expected flux at 95% confidence limit in one out of 30 bins (i.e., fewer than 5%) in agreement with expectation. Therefore, the light curve is consistent with the observed signal being associated with the diffuse Coma cluster ICM emission.

Systematic effects
Our Fermi-LAT analysis relies on choices that are somewhat arbitrary. To further validate the Fermi-LAT results, we thus check the impact of various systematic effects associated with these choices. They are listed below and summarized in Table 2 in terms of the changes on the flux of the cluster emission. We also refer to Xi et al. (2018) and Ackermann et al. (2016) who tested the impact of similar systematic effects in the region around Coma. As discussed below, the uncertainty in the diffuse background is the dominant systematic effect; it is expected to remain below 40%.
Diffuse background emission Although Coma is located near the galactic north pole (i.e., in a very clean region regarding diffuse galactic emission), the diffuse background is the dominant contaminant (see Figure 5) and is slightly correlated with the cluster signal (see Figure 3). Given the location of Coma, we also note that the galactic emission is nearly isotropic in the ROI, which explains the high degeneracy between the diffuse isotropic component and the diffuse galactic component. Therefore, the diffuse background modeling might lead to a significant systematic effect. The diffuse background related systematic effect was investigated by Ackermann et al. (2016) who showed that above 300 MeV, the systematic effect was less than 22% (although could reach ∼ 50% below 300 MeV). Similarly, Xi et al. (2018) concluded that the uncertainty in the diffuse background was < 30% within 0.2-300 GeV.
To test the impact of the diffuse background, we first reproduce our results by fixing the background to its expected model value (i.e., normalization set to 1 and no spectral index variation allowed). In this case, we obtain TS = 60.0 and the flux of Coma increases by 90% (scenario 2). Such an increase is expected as the background mis-modeling can be partially absorbed by the cluster template. In fact, we note that in this case the residual image presents a significant positive offset, which is partially absorbed by the cluster model and explains the boost of the signal. Such values for the Coma cluster can thus only be taken as upper limits.
We then reproduce our results using previous background models, namely gll iem v06 and iso P8R3 SOURCE V2 for the galactic and isotropic emission, respectively. In scenario 2, the TS value slightly increases, to 31.93, and the flux remains stable within < 4%.
We also consider the 16 alternative background models discussed in Xi et al. (2018). They are obtained using the GALPROP web interface for computation (Vladimirov et al. 2011) 7 , given the parameter definition files provided by Ackermann et al. (2012). We obtain a set of 16 versions of the diffuse galactic background components (and the corresponding isotropic diffuse emission that was model simultaneously) including bremsstrahlung, inverse Compton and pion decay emission. We first fit the sum of the galactic background components with a single free normalization and spectral slope parameter. In this case, the TS value is systematically higher than in our baseline fit, ranging from 33.19 to 39.26, with a flux from 20 to 40% higher. Then we consider a free normalization and spectral slope for each individual components separately, but do not consider the subcomponents separately (e.g., inverse Compton from optical, far infrared and CMB scattering individually). The TS values range from 30.98 to 39.42 and the flux is stable within [−8.5, +21.8]%. We also note that the morphology of the residual signal does not change significantly depending on the considered background model.

Energy range
We test the impact of our choice of the considered energy range. As the systematic effects are expected to be dominant at low energy, we test changing the nominal minimum energy of 200 MeV to 500 MeV. We also consider the case of 100 MeV as the low energy threshold. We extrapolate the flux assuming our physical model with a CRp spectral index in the range α CRp = [2.6 − 3.2], and also considering power law extrapolation with photon index from 2.25 to 2.65. The statistics is slightly reduced and our results change within 63%. However, Notes. ( ) The flux variation is computed by extrapolating the model in the range 200 MeV -300 GeV. The flux uncertainty is largely driven by uncertainties in the extrapolation.
we note that this large number is mainly driven by the uncertainty in the extrapolation when converting the flux to the 200 MeV -300 GeV band, and the significance of the detection remains stable within 6.3%.
Energy and spatial binning Our default binning choice is 0.1 degree per pixel, and eight energy bins per decade. We vary the pixel resolution from 0.05 to 0.2 degree and the energy binning from 5 to 12 bins per decade. The changes are less than 6% overall.
ROI size We use a ROI size of 12 by 12 degrees. To make sure that this choice does not introduce any significant bias, e.g., due to the presence of poorly constrained sources at the periphery of the ROI, we reproduce the analysis using a ROI of 8 and 15 degrees. This does not change our results by more than 13%.
4FGL source selection size Similarly, the choice of including sources from the 4FGL catalog that are within a 20 degrees width region from the ROI center is checked by increasing this value by 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., the cleanest Pass 8 event class). We also use the FRONT only (event type 1) converting photons. The changes are less than 13%.
Rocking angle cut It is recommended by the Fermi-LAT collaboration to check the impact on the rocking angle cut on the results. As a baseline, we apply a cut on rocking angle less than 52 degrees. We reproduce the results presented here without considering this cut. While the statistics slightly increases, the changes are less than 4% on the results.

On the nature of the signal toward the Coma cluster
In Section 4, we have shown that a significant γ-ray signal is observed within the characteristic radius θ 500 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+2736, as modeled in the 4FGL-DR2 catalog (a point source), is strongly degenerate with the diffuse cluster models that we have tested, and could correspond to the peak of the diffuse ICM emission. The comparison between the models based on a single cluster diffuse ICM compo-nent and a single point source (4FGL J1256.9+2736) cannot be strongly discriminated, although diffuse models based on multiwavelength templates provide the best agreement with the data. Models with two components (diffuse ICM plus point source) better match the data, but the improvement is marginal given the additional component involved.
In the next sections, we will explore the consequences for the cosmic ray population in the Coma cluster of the scenarios in which the signal is (at least partly) associated with the diffuse ICM. We will thus consider the scenarios 2 in which the signal is entirely associated with the diffuse ICM, and the scenario 3 in which the diffuse ICM and an independent point source (4FGL J1256.9+2736) both account for the signal. However, given the data, it is not possible to exclude that the signal arises essentially from a point source independent from the cluster diffuse ICM (scenario 1), but this will not be considered in the following since no CRp are needed in this scenario.

Multiwavelength comparison
In this Section, we qualitatively compare the Fermi-LAT excess map obtained in Section 4 to data at other wavelengths, as described in Section 2. The interpretation assumes that all the γ-ray emission arises from the diffuse ICM component.
First, we briefly compare the Fermi-LAT TS map (baseline model) to the optical image constructed using SDSS data in Figure 8. This provides a visual appreciation of the scales probes by Fermi-LAT. The galaxy density contours are also shown for visual purpose. As we can see, the TS peak is located about 10 arcmin north of the NGC 4839 group and about 20-30 arcmin from the Coma center and its two brightest galaxies. However, this offset is small compared to the Fermi-LAT angular resolution. The TS value remains larger than 16 (about 4σ) in most of the region where the bulk of the galaxy is located. However, the Fermi-LAT excess also extends further in the southwest direction, as discussed below.
In Figure 9, we compare the TS map to images of the tSZ, X-ray, galaxy density, and radio signal. The northwestsoutheast elongated morphology of the TS map (and excess counts) matches well what is seen at other wavelengths. The best match is observed for the galaxy distribution and the tSZ map. The later being proportional to the product between the thermal gas density and temperature, this indicates that the CRp distribution matches the temperature well. The temperature being fairly constant up to ∼Mpc scales, this would favor relatively flat CRp distributions. Alternatively, it could favor a scenario in which the signal comes from the sum of unresolved sources, as traced by the galaxies within the cluster region, and is not necessarily asso- The thin gray contours provide the galaxy density as shown in Figure  9. The white contours give the TS map levels of 2, 4, 9, 16 and 25. The image was constructed using the g, r and i filters of SDSS. ciated with the diffuse ICM component. The X-ray morphology is more compact than the Fermi-LAT excess, although the difference could be largely due to the relatively poor angular resolution of Fermi-LAT. The γ-ray excess matches the radio halo, but also extend toward the relic and could thus provide a good match to the combination of the two. In the case of an association with the relic, it would suggest that a significant fraction of the signal arises from inverse Compton emission toward the relic because very little target thermal gas is available for hadronic interaction in the relic region, but we leave this interpretation for future work (see Brown & Rudnick 2011;Ogrean & Brüggen 2013;Akamatsu et al. 2013, for discussions about a possible shock at the relic location and the accretion shock interpretation for its origin).
As noted in Figure 8, the TS peak is slightly off centered with respect to the cluster center, and better coincides with the southwest extension associated with the merger with the NGC 4839 group (between the halo and the relic). It also coincides well with the location where a radio front was identified (Brown & Rudnick 2011). This front is coincident with a discontinuity seen in the X-ray surface brightness, temperature and entropy (Simionescu et al. 2013;Uchida et al. 2016;Mirakhor & Walker 2020) and SZ signal (Planck Collaboration et al. 2013), possibly suggesting that the local CR might have been accelerated by a shock. However, given the significance of the excess and the Fermi-LAT angular resolution, this remains difficult to interpret further, as also shown in Section 4 and in particular in Figure 6.
In conclusion, we find good morphological agreement with data at other wavelengths. Nevertheless, the interpretation remains difficult due to the low significance of the excess and the Fermi-LAT angular resolution. This comparison does not allow us to exclude that the signal observed is entirely, if any, associated with the diffuse ICM emission. It could also be the sum of several components.

Implications for the cosmic ray content of the Coma cluster
In this Section, we use the γ-ray SED extracted in Section 4, together with the model described in Section 3, in order to constrain the CR population in the Coma cluster. Although they provide the best match to the data, the multiwavelength data templates are not used because they do not allow us to have a threedimensional physical model of the cluster, which is needed to constrain the CR content.

Methodology
We aim at using the Fermi-LAT extracted SED to constrain the CRp population of our model. As discussed in Section 3, and in more details in Adam et al. (2020), the hadronically induced γray emission depends on the thermal gas (modeled via the pressure and density), and the CRp population (spatial and spectral distribution). The thermal gas pressure and density are well constrained from Planck and ROSAT data, respectively, and are thus kept fixed in this analysis. Since the CRp spatial distribution was kept fixed when extracting the SED, we keep it fixed when constraining the parameters based on the SED fit. However, we consider the different spatial shapes as in Section 4 because we have seen that it was not possible to discriminate between the different cases. We also consider the case of scenario 3, where both the diffuse ICM and 4FGL J1256.9+2736 are included in the sky model. We are left with two parameters to be constrained: 1) the CRp normalization defined as the energy stored in the CRp relative to the thermal energy X CRp (R 500 ), which we define at a radius R 500 ; 2) the slope of the CRp energy spectrum α CRp . These two parameters should be constrained for all the different models considered (compact, extended, flat, and isobar), which will provide an assessment of the systematic effect associated with the spatial model. We employed a Markov chain Monte Carlo (MCMC) approach in order to constrain the parameters space, using the emcee package (Foreman-Mackey et al. 2013). In brief, the chains move in the parameter space according to a proposal function and the likelihood of the model given the data. We adopted flat priors across the range X CRp (R 500 ) ∈ [0, 0.2] and α CRp ∈ [2, 5], which corresponds to the physically acceptable parameter range, but we checked that this limit does not affect our results. This method allows us to efficiently find both the best-fit parameters (as the parameters that maximize the likelihood) and the estimate of the posterior probability distribution.
The likelihood function is defined as where i runs over each energy bin and the model parameters are θ ≡ [X CRp (R 500 ), α CRp ]. The value L i is the likelihood of measuring a given flux normalization in the energy bin i, depending on the model flux M i ( θ). It is obtained by interpolating the likelihood scan, as provided by Fermipy, when extracting the SED in Section 4.5. Once the chains have converged, and after removing the burn-in phase, the two-dimensional histogram in the plane X CRp (R 500 ) − α CRp can be integrated to provide the constraints up to a given confidence interval. The individual chain histograms also provide the marginalized posterior probability distribution of each parameter. The integrated posterior probability distribution up to 68% probability gives the uncertainties. In addition 14 R. Adam et al.: γ-ray emission toward Coma . 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 gray circle provides the radius θ 500 and 3 × θ 500 . Several relevant features are also indicated in orange. For display purposes, the WSRT image has been apodized at large radii to reduce the larger noise fluctuations present on the edge of the field. As a complementary figure, Figure 8 provides an optical image of the central region.
to the model parameters themselves we obtain the constraint of the spectrum. To do so, we compute the model SED for each set of parameters and calculate the envelope of all the models as the 68% confidence limit measured from the model histogram in each energy bin. The same procedure is applied to measure the γ-ray flux (and luminosity) between 200 MeV and 300 GeV according to the MCMC sampling.

Constraints on the cosmic ray proton population
The SED measured in the case of our baseline model is shown in Figure 10 for both scenario 2 and 3 (see subsection 4.2 for the detailed definition). Error bars are the 1σ error on the SED as evaluated from the likelihood scan curvature. The maximum likelihood model is also shown in black, as well as 100 models randomly sampled from the MCMC chains, their median and the associated 68% confidence interval. We can observe that the best-fit model is in good agreement with the data in both cases, as also highlighted by the residual. The model is relatively well constrained around 300 MeV -1 GeV, but the uncertainties remain large at larger energies. The peak SED, around 500 MeV, reaches about 3 × 10 −7 MeV cm −2 s −1 . In the case of scenario 3 (right panel, both the ICM and point source in the model), we can see that the spectrum amplitude is reduced and that error contours are significantly larger.
The corresponding constraint on the posterior likelihood of the parameters X CRp (R 500 ) and α CRp is shown in Figure 11. The constraints on the model SED lead to a relatively tight constraint on the normalization, but the constraint on the slope remains fairly loose. The two parameters are degenerate as increasing the normalization and the slope simultaneously does not strongly change the flux at high energies (see also Figure 1, bottom right panel). The constraint on the CRp to thermal energy is about 1.8% and the slope about 2.8, as also shown on the marginalized distributions. When both the ICM component and 4FGL J1256.9+2736 are included in the sky model, the con- The best-fit model is shown in black and the 68% confidence interval is show in blue, together with 100 models randomly sampled from the MCMC parameter chains. The median of these models is also shown as a dashed blue line. The residual provides the difference between the data and the best-fit model normalized by the error bars. straint on the normalization shifts toward zero and the posterior only excludes X CRp (R 500 ) = 0 at about 2σ. The uncertainty on the slope increases and the best-fit slightly decreases to about α CRp 2.6. In Table 3 we provide the MCMC constraints on these parameters in the case of all tested spatial models. The fluxes and corresponding luminosities are also constrained according to the MCMC fit of the model. Given the uncertainties, all models are in agreement and the associated systematic shift on the parameter is about 20% on the CRp to thermal energy ratio, and about 7% on the slope. When including the source 4FGL J1256.9+2736 and the cluster simultaneously in the sky model, the slope is slightly reduced (albeit with increased uncertainty) and the normalization is reduced by about 70%.

Implications for the cosmic ray electrons
The presence of CRp, as traced by the hadronic γ-ray emission, implies the production of secondary CRe, which should contribute to the observed radio emission. We compute both the profile and the spectrum associated with this population given the 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 parameters sampled in the MCMC to obtain the 68% confidence limit on the radio emission.
As discussed in Section 3, the primary CRe, i.e., CRe 1 , are also included in our model and we consider two cases for interpreting their origin. 1) Since our model does not include explicitly any reacceleration (see Brunetti & Lazarian 2007, for reacceleration models), the population of CRe 1 that we constrain is supposed to be independent of the CRp and the CRe 2 . It would correspond, for instance, to a CRe 1 population arising from star formation activity spread over the cluster volume, or the direct shock acceleration in the ICM. In this case, the CRe 1 and CRe 2 populations coexist and the radio emission accounts for the sum of the two. 2) We could also interpret the CRe 1 population in our model as the reaccelerated CRe 2 population. In this case, the CRe 2 would be the seed for the CRe 1 . The CRe 2 population should thus not contribute directly to the radio emission. Instead the ratio CRe 1 /CRe 2 would measure the amount of reacceleration needed to explain the emission with respect to the purely hadronic steady state scenario, as a function of energy and radius. We stress that in this second interpretation, we only  Figure 11: MCMC constraints on the model parameters, in the case of the reference model (n CRp ∝ n 1/2 e ) when 4FGL J1256.9+2736 is excluded from the sky model (green) or included (orange). The bottom left panel provides the constraint in the plane X CRp (R 500 ) − α CRp , where the contours correspond to 68% and 95% confidence interval. The marginalized posterior probability distributions are also shown for the parameters X CRp (R 500 ) (top) and α CRp (right), where the shaded area provides the 68% confidence interval. provide a constraint relative to the hadronic steady state scenario since this is an assumption made when computing the CRe 2 population. In contrast, reacceleration models do not assume steady state, and both CRp and CRe populations evolve together according to turbulent reacceleration. Nevertheless, we consider this second case as it is still instructive regarding reacceleration physics.
The CRe 1 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 a normalization X CRe 1 (R 500 ), a spectral slope α CRe 1 , a break energy (E cut,CRe 1 ), and a spatial scaling relative to the thermal density (η CRe 1 ). We fit for these parameters using simultaneously the radio spectrum and the 352 GHz radio profile data (see Section 2). Regarding the spectrum, we compute our model using a cylindrical integration within R = 0.48 × R 500 (629 kpc, 0.36 deg) as it corresponds to the extent of the radio halo (Brunetti et al. 2013). We note, however, that the radio spectrum data are not strictly homogeneous in terms of aperture radius used for flux measurement and our model value only provides an effective radius. Because the profile is extracted from a single instrument (WSRT), which arguably could be the best one in terms of capturing the diffuse emission, but which does not necessarily perfectly agree with other measurements (as seen in Figure 12), we also allow for a cross-calibration of the profile measurement by adding an extra normalization which we fit simultaneously. The fit is performed with the least square function curve fit from the scipy python package. Depending on the considered case, the radio model includes both the CRe 1 and CRe 2 contributions (case 1, no reacceleration), or only the CRe 1 as the reaccelerated CRe 2 (case 2, pure reacceleration). We compute the error contour on the CRe 1 fitted population by reproducing the fit in the case of the lower and upper bounds for the CRe 2 population. We thus assume that the uncertainties from the CRe 2 population (given by the γ-ray) are dominant over the uncertainties associated with the radio data.
In Figure 12, we show the constraint on the CRe population from the radio synchrotron spectrum and profile fit in the first case using the ExponentialCutoffPowerLaw model. We note that the Figure is provided for our baseline CRp radial model (n CRp ∝ n 1/2 gas ) in the scenario in which 4FGL J1256.9+2736 was replace by the ICM component (scenario 2). In Appendix D, we also provide these constraints in the case of the alternative mod- Radio spectrum of Coma, as compiled from Pizzo (2010) and constraint from the reference CRp spatial model (n CRp ∝ n 1/2 e ) and the reference CRe 1 spectral model (ExponentialCutoffPowerLaw). The measurement from the Brown & Rudnick (2011) data is also shown as the magenta diamond. The contribution from the CRe 2 is shown in blue together with its 68% confidence interval, and the remaining contribution from CRe 1 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 fluxes are computed using cylindrical integration within R = 0.48 × R 500 = 629 kpc ≡ 0.36 deg. Top right: Radio profile measured from the WSRT map at 352 GHz and comparison to the reference model. The contributions from CRe 2 and CRe 1 are as in the left panel. Bottom left: Absolute number density distributions of CRe 1 and CRe 2 taken at 1 GeV and 10 GeV. Bottom right: Ratio between the CRe 1 and CRe 2 number populations. We note 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 a dashed gray line given our pressure model, is not negligible for the highest frequency data point, but we correct for it (see also Brunetti et al. 2013, for a dedicated analysis). Our model provides a reasonable fit to the data for both the spectrum and the profile (this is also the case for the other considered models, see Appendix D, and our results do not depend significantly on the considered CRe 1 spectral model). The slope of the synchrotron emission from CRe 2 is similar to the one from the CRe 1 , but it is significantly less curved and no high energy cutoff is present in the CRe 2 distribution. We can see on the spectrum that the radio emission within 0.48 × R 500 (629 kpc, 0.36 deg) from CRe 2 is overall a factor of about four lower than the total emission (except at high frequency where it reaches similar values). On the profile, the CRe 2 emission is significantly more concentrated than the total radio signal in the case of this CRp spatial model (although it is also true for all CRp model, see Appendix D) and lead to slightly over-fitting the total synchrotron emission in the core when added to the CRe 1 contribution. This high concentration is expected because the CRe 2 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 and the total CRe density, is so flat that it requires a nearly flat CRe distribution given the fixed magnetic field profile (as also noted in . Thus, this could be achieved for the CRe 2 only at the cost of an inverted CRp profile (rising with radius). The number density of CRe is at a level of about 10 −14 and 10 −17 cm −3 GeV −1 at 1 GeV and 10 GeV, respectively, for both populations, but with opposite radial dependences. These constraints on the CRe populations translate into a ratio CRe 1 to CRe 2 that increases with radius and that does not depend much on energy up to 10 GeV (before the cutoff; best-fit E cut,CRe 1 = 17 GeV). Given this CRp spatial model, the CRe 1 /CRe 2 ratio is slightly below unity in the core, and strongly rises to reach about 100 at 2 Mpc.
In Figure 13, we show the constraint on the CRe population in the second case (CRe 2 are the seed to the CRe 1 , pure reacceleration). As for Figure 12, the alternative models considered in this paper are shown in Appendix D. The model also provides a good fit to the data as shown in the top panels, where only 18 R. Adam et al.: γ-ray emission toward Coma Figure 13: Same as Figure 12 in the case of scenario 2 (CRe 1 interpreted as the reaccelerated CRe 2 population). The CRe 2 prior reacceleration refer to the steady state hadronic model. The reacceleration boost is given relative to the steady state hadronic model.
the CRe 1 population (as the reaccelerated CRe 2 population) is included. On the bottom panels, the interpretation is now different as the amount of CRe 1 now correspond to the CRe 2 population after reacceleration. As can be observed, the best-fit CRe 1 profile is nearly flat, while the original seed population is more concentrated in the core. The amount of reacceleration relative to the hadronic steady state case is thus relatively low in the core (in agreement with no reacceleration there, or even favoring a boost lower than one for the most compact CRp profile), but strongly increases in the outskirt reaching about a factor of 100 at 2 Mpc. As in the previous case, this radial dependence depends on the considered model for the CRp distribution and we show in Appendix D that flatter is the CRp distribution, flatter will be the reacceleration boost profile. Nevertheless, in all the considered cases, the reacceleration boost (relative to the steady state hadronic model) increases with radius. The energy dependence of the boost, comparing the values for 1 GeV and 10 GeV electrons, is not much affected by the choice of the spectral model for the CRe 1 population.
To summarize how the radio data were used, we used the WSRT data to establish the radio profile, but for the spectral dependence, we were forced to use only the data within 0.48 × R 500 where measurements at other frequencies were available. The larger halo size and flux from the WSRT data, which are due to its increased sensitivity, are discussed in detail by Brunetti et al. (2013). Those authors showed (in their Figure 2, right) that the observed size of the halo was a strong function of the signal to noise ratio, and that the WSRT data were therefore the most reliable. The magenta diamond in Figure 12 and 13 shows how the WSRT flux is above the rest of the spectral points. In order to examine the effects of increasing the entire spectrum by a factor of two to match the WSRT point, we added an extra normalization parameter, as discussed earlier in this Section. We find that the ratio between the total radio emission and that arising from CRe 2 would increase from 4 to 8. However, the effect of changing such normalization would not impact the conclusions of this paper.
In Appendix D, we also provide the results when including both the ICM component and 4FGL J1256.9+2736 in the Fermi-LAT sky model (scenario 3). While the CRe 2 synchrotron spectrum is slightly steeper and the amplitude of the CRe 2 component is reduced, the shapes of the CRe 2 /CRe 1 profile (case 1), or reacceleration boost profile (case 2) are not significantly changed given the large uncertainties.
Finally, in order to compare the distribution of CRp that we measure to that expected in reacceleration models to explain the radio emission, we compare our CRe 2 induced synchrotron spectrum to that of the model developed by Brunetti & Lazarian (2011) in Figure 14. To best match the model developed by Brunetti & Lazarian (2011), we use the ICM model in which the CRp radial profile is flat, but we note that the two models are not strictly equal. For instance, Brunetti & Lazarian (2011) use an isothermal β-model, while we set the thermodynamic profiles to X-ray and tSZ data (see Section 4.3.1 from Brunetti & Lazarian 2011, for more details). As the reacceleration model was calibrated on the Coma cluster, the total radio synchrotron (solid brown line) matches well the data, as expected. When the reacceleration is switched off (dashed brown line), only the emission from secondaries directly produced from hadronic collisions remains, which compares very well to the prediction from our model. This shows that the distribution of CRp that we measure provides an excellent match to what is needed in the reacceleration model developed by Brunetti & Lazarian (2011) to explain the overall radio spectrum, when including reacceleration.

Comparison to previous analysis
Constraints on the γ-ray emission of the Coma cluster have been obtained in earlier work using Fermi-LAT data (see Arlen et al. 2012;Prokhorov 2014;Zandanel & Ando 2014;Ackermann et al. 2016;Keshet & Reiss 2018;Xi et al. 2018). In general, the signal is modeled using spatial templates together with power law distribution for the photon spectrum. Our results push forward such analysis by directly using a physical model for the CRp and thermal gas. Thus we are able to directly constrain the physics of CR by comparing model and data, and push the analysis by investigating the implications of γ-ray emission for the CRe population. Our results are in broad agreement with previous searches and limits. Nevertheless, we stress that the reported CRp to thermal energy (or pressure) are generally obtained under the assumption of a given spectral distribution (harder the spectrum and more stringent the limit). Spectral index values used in the literature correspond to spectra that are generally much harder (∼ 2.1 − 2. 3) than what we obtain here (see also below for further discussions).
The Fermi-LAT collaboration observed γ-ray excess within the virial radius of Coma using six years of data (Ackermann et al. 2016). However, the signal was too faint for detailed investigation and they published upper limits on the signal for various templates. Our results are consistent with that of Ackermann et al. (2016), although using twice the amount of data and an analysis that differs in various aspects. They provide an upper limit of 5.2 × 10 −9 ph cm −2 s −1 (E > 100 MeV). Extrapolating our best-fit model in the same energy range, we obtain a compatible total flux of 2.1 × 10 −9 ph cm −2 s −1 (in scenario 2, to match Ackermann et al. 2016). Among the main differences, we note that they use a power law emission model with a spectral index of 2.3 and spatial distribution based on WSRT, while our baseline model is directly connected to the underlying CR population and the spectrum that we measure is significantly steeper than their model. Given this framework, they obtain a TS value of 13, compared to about 27 in our case.
The detection of γ-ray emission toward Coma was first claimed by Xi et al. (2018), who used an unbinned likelihood approach (reaching TS values of about 40-50 depending on the model). While the morphology of the signal and the spectral slope of the emission that we observe agree with their results, we obtain fluxes that are significantly lower (about a factor of two), depending on the exact model used. For instance, they obtain fluxes of about 2.3−3.1×10 −9 ph cm −2 s −1 for disk, core, or radio and X-ray based templates (E = [200, 300] GeV). For the same energy range and similar models, our fluxes are 1.3 − 1.6 × 10 −9 ph cm −2 s −1 . We note that their fluxes would also be excluded by the Ackermann et al. (2016) limit when extrapolating down to 100 MeV using a photon power law index of 2.7-2.8, and would nearly reach the limit when extrapolating with our best-fit model. Keshet & Reiss (2018) claimed the detection of a ring-like signal at a position that correspond to the expected location of the accretion shock (see also Hurier et al. 2019, for the detection of such an accretion shock with Planck in A2319). In contrast, our results do not show any ring-like structure, especially when looking at the excess profile around 3 Mpc radius (about 2 degrees).

Cosmic ray physics
The amount of CRp and their spectral and spatial distributions has been predicted from numerical simulation (e.g., Pfrommer et al. 2007;Pinzke & Pfrommer 2010). The amplitude of the CRp that we obtain, relative to the thermal energy, is in line with predictions from such simulations. We obtain a ratio of about 1.5% (about 0.8% when including 4FGL J1256.9+2736), while simulations suggest a few percent. For instance, in Ackermann et al. (2014), based on the work by Pinzke & Pfrommer (2010), the CRp to thermal energy ratio of Coma was expected to be X CRp (R 200 = 2.2 Mpc) = 2.4 × 10 −2 (we account for the fact that the pressure ratio that they use is half the energy ration that we use). This value is just a factor of ∼ 2 above the one that we constrain, although it was obtained for a different CRp index and may also depend on radius so that the comparison depends on the exact CRp spatial distribution. This CRp to thermal energy ratio is related to the CRp injection efficiency for the diffusive shock acceleration process, which increases with the Mach number. Thus, it gives access, in principle, to the microphysics of CRp acceleration. We refer to Ackermann et al. (2014) for more discussions.
As also found by Xi et al. (2018), the spectral index of the CRp that we constrain is significantly larger than what is generally assumed (α CRp ∼ 2.5 − 3.2 for our constraints versus 2.1 − 2.3; see e.g., Arlen et al. 2012;Zandanel & Ando 2014;Pinzke & Pfrommer 2010). For a given shock, the CRp index is related to the Mach number: higher the Mach number and harder the spectrum (see Pfrommer et al. 2006, for discussions). Therefore, our results could point to shock acceleration for which the Mach numbers are overall smaller than usually expected.
The shape of the CRp is often assumed to follow the shape of the thermal density profile, with possible diffusion sometimes included (e.g., . While our results favor models with intermediate scaling (n CRp ∝ n ∼1/2 e ), the data are not sufficient to firmly constrain the shape of the CRp distribution. In the context of the reacceleration of CRe 2 , a compact CRp profile would imply that reacceleration strongly increases with increasing radius, as already highlighted in Pinzke et al. (2017). However, even flat CRp profiles lead to an increasing reacceleration boost with radius.
Our results are in agreement with earlier work that have shown that pure hadronic models were excluded given the magnetic field strength inferred from Faraday rotation measures (Brunetti et al. 2012;Zandanel & Ando 2014;Brunetti et al. 2017). Nevertheless, the hadronically induced CRe are only a factor of a few lower than the amount required, and could thus provide seeds for turbulent reacceleration models (Brunetti & Lazarian 2007. Throughout this work, we have fixed the magnetic field strength model to the best-fit result obtained by Bonafede et al. (2010), despite the fact that there are relatively large uncertainties in the constraint. We also refer to the recent work by Johnson et al. (2020) who showed that even under ideal condi- tions, the central magnetic field cannot be determined to better than a range of 3, with corresponding uncertainties in the scaling parameter η. These uncertainties affect the results presented in Section 6.3. In our work, increasing (decreasing) the magnetic field would imply more (fewer) synchrotron emission for a given CRe population. Thus it would lead to fewer (more) CRe 1 for case 1, or fewer (more) reacceleration boost in case 2. This effect also applies to the radial dependence. Given the uncertainties in the magnetic field measure, these should contribute significantly to our constraints. The uncertainties in the CRe 2 population are nonetheless expected to dominate, and we leave the joint investigation of the CRe population and magnetic field distribution aside.

Contamination from discrete sources
In this paper, we have generally assumed that the γ-ray emission observed in the direction of the Coma was originating from hadronic interactions in the ICM. We have also considered the case where 4FGL J1256.9+2736 was included in the overall model as a point source contaminant in addition to the diffuse emission, but still modeling the diffuse component as from hadronic interactions.
In fact, cluster member and radio galaxies may also contribute significantly to the total signal. In Ackermann et al. (2016), a minimum flux was estimated for the two dominant central radio galaxies NGC 4869 and NGC 4874 assuming that the electrons responsible for the radio emission also generate inverse Compton emission in the Fermi-LAT band and assuming simple scaling relations for the calculation. They obtained luminosities of about 6×10 40 and 2×10 40 erg s −1 (0.1 < E < 10 GeV), which is a factor of ∼ 20 lower than what we measure for the hadronic emission at E > 0.2 GeV (thus even lower in our energy range for usual spectra).
However, other sources may also contribute and might lead to an unresolved diffuse component. In particular star forming galaxies could lead to γ-ray emission for which the associated flux is very uncertain and within the range of our constraints (in the range 3 × 10 40 − 3 × 10 42 erg s −1 for energies in 0.1-100 GeV; see Storm et al. 2012).
The cluster member radio and star forming galaxies are also expected to generate CR that diffuse in the ICM and would contribute to the population which we model in this paper. Rephaeli & Sadeh (2016) calculated that they might account for a significant amount of the radio diffuse emission, as well as in the γ-rays.

Summary and conclusions
This paper presented the analysis of nearly 12 years of Fermi-LAT data toward the Coma cluster, together with multiwavelength data and using the MINOT software in order to model the signal. Different scenarios were considered to model the signal: no cluster diffuse emission, replacing the source 4FGL J1256.9+2736 by a diffuse cluster model, or accounting for a combination of both. Assuming that the diffuse emission was associated with the hadronically induced γ-rays in the ICM, we investigated the implications for the CR physics, for both CRp and CRe.
The signal was modeled assuming that the γ-ray emission arises from hadronic interactions between CRp and the thermal gas. The thermal gas model was set using ROSAT X-ray and Planck tSZ data. In this model secondary CRe are also produced, and we fixed the magnetic field to that obtained from Faraday rotation measurements in order to compute the associated radio synchrotron emission. The CRp spatial model was calibrated assuming a scaling relative to the thermal gas profiles. The CRp normalization and spectrum were defined relative to the thermal gas energy and using a power law for the spectrum, 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 cluster. The detection level depends on the model considered, in the range TS ∼ 24 − 34, corresponding to a significance of about 4.9 − 5.8σ. While extended models provide a better description of the data, it is not possible to strictly exclude that the signal is associated with a point source, or a combination of the two, and we include this possibility in our analysis. The morphology of the signal is elongated in the northeast to southwest direction, in agreement with other wavelengths. The peak of the emission is about 0.5 degree offset in the southwest direction with respect to the X-ray peak, and coincides with the well-know merger extension.
Using an MCMC approach, we constrained the amplitude and the spectral index of the CRp population in the cluster assuming that at least part of the signal was associated with the diffuse ICM hadronic interactions. We find that the energy stored in the CRp is about 1.5% of the thermal energy of the Coma cluster (0.8% when including 4FGL J1256.9+2736 and the cluster in the model). The slope is larger than what is usually assumed, around 2.8, although with large error bars. In the framework of diffuse shock acceleration, this could implies that CRp acceleration 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, leading to a radio synchrotron emission that is about 4 times lower than what is observed. While a pure hadronic origin of the radio emission is ruled out, these secondary CRe could serve as the seeds for turbulent reacceleration. In this model, the reacceleration should increase with radius, depending on the exact CRp spatial distribution. Alternatively, an independent CRe population could be at the origin of the remaining radio emission, but it would require a nearly flat radial distribution.
Our results show that after almost 12 years of observations, the diffuse γ-ray emission from galaxy clusters might now become accessible with the Fermi-LAT. Since the hadronic emission is expected to be a universal property of galaxy clusters, our results might be reproducible for other clusters, renewing the interest of such analysis, although Coma might be the best target for such searches. In the scenario in which the signal indeed arises from diffuse ICM hadronic interactions, if the large value of the CRp spectral slope is confirmed and if the CRp content of clusters is, as expected, a universal property, the very high energy γ-ray emission (∼TeV) could therefore be much lower than usually assumed. In this context, it would be challenging for future ground-based γ-ray observatories (e.g., the Cherenkov Telescope Array, Cherenkov Telescope Array Consortium et al. 2019) to detect the diffuse emission associated with hadronic interactions in the ICM of galaxy clusters. 10 3 10 4 Frequency (MHz) 95% containment). This can be compared to the offset of 0.68 deg between 4FGL J1256.9+2736 and the cluster center. Given the low S/N and the fact the ICM is an extended source (e.g., surface brightness drops by a factor of ∼ 2 at θ = 0.2 degrees, in the extended case, see Figure 1), this agrees with the offset being not significant.
Although a diffuse emission plus point source model provides the best match to the data (see Section 4.4), these results confirm that the scenario 2 agrees with the data. This also agrees with the results presented in the main paper, in which no significant residual is observed around Coma once the best-fit ICM model is subtracted from the data, and where the comparison of the TS given in Table 1 does not point toward the need of including 4FGL J1256.9+2736 in the sky model if we include a diffuse ICM component. See also Section 4.4 for discussions.

Appendix C: ROI modeling and fitting: A null test in the direction of point sources
In this Appendix, we present a null test, which we use to validate our ROI modeling and fitting procedure. As discussed in the main text, the source 4FGL J1256.9+2736 can be better described by an extended ICM model ( Table 1) that is centered 0.68 degrees away from the source, leading us to consider the scenario in which 4FGL J1256.9+2736 corresponds in fact to the peak of the diffuse ICM emission (scenario 2). Here, we check that this is not a general feature of our procedure and that it does not apply to other sources. We select the following sources from our field because their local background is very close to the Coma region (no close by source in the 4FGL catalog). However, we note that their spectra are generally steeper and their TS may differ: 4FGL J1253.8+2929 (TS = 27.57, 2 degrees northwest from the Coma center, power law index of 1.9), 4FGL J1250.8+3117 (TS = 137.29, 3.8 degrees northwest from the Coma center, power law index of 2.16), 4FGL J1316.5+3013 (TS = 20.90, 4.3 degrees northeast from Coma, power law index of 2.1). To match what is done for Coma, we define the cluster center 0.68 degrees away from these sources and test four directions for the offset (north, east, south, west) to increase statistics. We also perform the test with the cluster centered on the source, but we note that this does not match the case of Coma. Because the spectra of these sources are not the same as 4FGL J1256.9+2736, we model the cluster with a power law of free index instead of a fixed physical spectrum. However, we note that a flatter spectrum will lead to an improved angular resolution since the PSF decreases with increasing energy. We then apply the same fitting procedure as the one done for Coma, i.e., we add the cluster template in the sky model either by replacing the point source with the cluster, or including both the point source and the cluster. We use the baseline extended cluster model for the test performed here.
We report the different TS values that we obtain for the different test cases in Table C.1. We can see that for all tests, the TS is much smaller for the cluster than the point source when both are included, even when the two are co-aligned. In the cases where the point source is replaced by the cluster, the TS of the cluster increases, but remains much smaller than for the point source except when the cluster is centered on the source (although it always remains lower than for the point source only).
We conclude that while the point source 4FGL J1256.9+2736 can be better explained by an extended ICM model (even when both cluster and point source are included, see Section 4.3 and Table 1), this is not the case for the other sources 4FGL J1253.8+2929, 4FGL J1250.8+3117 and 4FGL J1316.5+3013. This provides us with a null test that allows us to strengthen the modeling and fitting procedure described in the main text.

Appendix D: Constraints on the cosmic ray electron populations with alternative models
This appendix provides the constraints on the CRe populations in the case of the alternative models considered in the case of the spectral modeling of the CRe 1 in Figure