© R. Adam et al. 2021

1. 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 associated 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⁸ 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 further 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, 2011). 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; MAGIC Collaboration 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₅₀₀ ≃ 7 × 10¹⁴ 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 θ₅₀₀ ∼ 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 X 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. 2021).

Several analyses have focussed on the search for γ-rays toward the Coma cluster, using both ground-based and space-based instruments (e.g., Perkins 2008; H.E.S.S. Collaboration 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 (Brunetti et al. 2017). 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 context 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 (Fermi/LAT Collaboration 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 Sect. 3, we present the physical modeling of the cluster and the computation of the observables. The Fermi-LAT analysis is presented in Sect. 4 and the outputs are compared to multiwavelength data in Sect. 5. Section 6 discusses the implications of the observed signal for the CR content of the Coma cluster. In Sect. 7, we discuss the results presented in this paper and compare them to the literature. Finally, a summary and conclusions are given in Sect. 8. Throughout this paper, we assume a flat ΛCDM cosmology according to Planck results (Planck Collaboration XIII 2016) with H₀ = 67.8 km s⁻¹ Mpc⁻¹, Ω_M = 0.308, and Ω_Λ = 0.692. The reference coordinates of the Coma cluster are RA, Dec = 194.9118, 27.9537° and the redshift z = 0.0231, based on the Planck catalog (PSZ2, Planck Collaboration XXVII 2016). Given the reference cosmological model, 1° corresponds to 1.73 Mpc at the redshift of Coma. We use the value of R₅₀₀ = 1310 kpc based on Planck Collaboration X (2013), corresponding to M₅₀₀ = 6.13 × 10¹⁴ M_⊙ given our cosmological model.

2. 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.

2.1. Fermi Large Area Telescope γ-ray data

The γ-ray analysis is based on the publicly available¹Fermi Large Area Telescope data (Fermi-LAT, Fermi/LAT Collaboration 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° from the cluster center were collected for the analysis (see also Sect. 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.

2.2. 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 line-of-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 XXII 2016) 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 X 2013, for a Planck analysis of the Coma cluster). The data are publicly available at the Planck Legacy Archive².

2.3. 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)³ 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, 2000).

2.4. 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 (Sect. 5). The image angular resolution is 4 × 3 arcmin² (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 Sect. 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₅₀₀ (i.e., 629 kpc or 0.36 deg, see Brunetti et al. 2013, for details).

2.5. 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⁴ to construct a color image of the Coma region. In addition, we select galaxies with spectroscopic information from the SDSS database⁵ 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 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⁻¹ (i.e., FWHM of 5000 km s⁻¹). 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.

3. 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 Sect. 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 (Sect. 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)⁶. 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.

3.1. Physically motivated modeling using MINOT

3.1.1. 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

(1)

with parameters taken from Planck Collaboration X (2013) and rescaled to our cosmological model (P₀, r_p, a, b, c = 0.022 keV cm⁻³, 466.8 kpc, 1.8, 3.1, 0.0).

The electron density is described by a β-model (Cavaliere & Fusco-Femiano 1978) according to

(2)

with parameters taken from Briel et al. (1992) and rescaled to our cosmological model (n₀, r_c, β = 3.36 × 10⁻³ cm⁻³, 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₅₀₀. This also apply for nonthermal quantities discussed below.

3.1.2. 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 (Sect. 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

(3)

where ⟨B₀⟩ = 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₀⟩ = 3.9 μG; η_B = 0.4) to (⟨B₀⟩ = 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 Sect. 7.2.

3.1.3. Cosmic ray protons

The CRp number density, per unit volume and energy, is modeled as

(4)

where f₁(E) is the energy spectrum, f₂(r) the radial dependence, and A_CRp is the normalization. The energy spectrum is modeled as a power law (expected from Fermi acceleration), as

(5)

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

(6)

Since the radial dependence of the CRp is not well-known, we test different values for the scaling, , corresponding to flat, extended, and compact distributions. Finally, the normalization A_CRp is computed relative to the thermal energy enclosed within the radius R₅₀₀. In the following, we use the parameter

(7)

where U_CRp(R) is the total CRp energy enclosed within R computed by integrating Eq. (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₅₀₀). The value of these parameters will be further discussed in Sects. 4 and 6.

3.1.4. 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 Sect. 6. The primary cosmic ray electrons, CRe₁, are modeled similarly as the CRp. However, CRe₁ are affected by losses and we therefore account for this using three different models: the ExponentialCutoffPowerLaw, the InitialInjection, and the ContinuouslInjection models as implemented in the MINOT software (Adam et al. 2020). They are expressed as

(8)

(9)

and

(10)

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₁: the spatial scaling of the CRe₁ population relative to the thermal gas, η_CRe1 (as in Eq. (6)), the CRe₁ spectrum slope α_CRe1 and break energy E_cut, CRe1, and the normalization encoded in X_CRe1(R₅₀₀) (as in Eq. (7)). The value of these parameters will be further discussed when comparing our model to radio data in Sect. 6.

3.1.5. γ-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₂), and the radio synchrotron emission (from both CRe₁ and CRe₂) 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₂ are computed as

(11)

where 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₂ is then obtained accounting for losses under the steady state assumption. The γ-ray emission profile and spectrum are computed by integrating Eq. (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⁴ GeV. These spatial (profile) and spectral templates are used in Sect. 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 Fig. 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₅₀₀) = 10⁻². 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₅₀₀ for different slopes α_CRp. The CRp slope mainly drives the γ-ray slope between 1 GeV and 10⁴ GeV. At higher energies, 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₅₀₀.

Fig. 1. MINOT templates used to describe the cluster, for different assumptions regarding the CRp distributions. All models are normalized so that XCRp(R500) = 10−2. Top left: radial profile of the CRp distribution. Top right: enclosed CRp to thermal energy profile. Bottom left: γ-ray surface brightness profile. Bottom right: γ-ray energy spectrum integrated within R500.

3.2. 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 Sect. 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² 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 (Sect. 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 Fig. 2 on 3 × 3 deg² grids (except for the relic, for which it is 1 × 1 deg²). 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 Sect. 5).

Fig. 2. Template images used to describe the morphology of the γ-ray emission. From left to right and top to bottom: templates 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 deg2, except for the radio relic (bottom right), which is only 1 × 1 deg2. 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.

4. 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 Sect. 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.

4.1. 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° 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² 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 Sect. 4.7).

4.2. 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° away from the center. Nevertheless, we include all sources within a square with a width of 20° from the reference center, because the point spread function (PSF) may extend up to about 10° at low energies.

As discussed in Sect. 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° from the cluster center and within θ₅₀₀. 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 Sect. 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 Sect. 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).

4.3. 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

(12)

with ℒ₀ the maximum likelihood value for the null hypothesis, and ℒ the maximum likelihood when including the additional source. (3) To allow for variability, we free all the model parameters associated with sources with , and free only the normalization of sources with . (4) We use the fit function of Fermipy to perform the maximum likelihood model fit of all the free parameters in the sky model. (5) We add the cluster diffuse emission (except in the scenario 1, see Sect. 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 Fig. 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 Sect. 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 Sect. 4.7 for systematic effects associated with the diffuse background).

Fig. 3. Correlation matrix associated with the likelihood fit. Each entry corresponds to one free parameter of the model. The names isodiff and 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 .

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 ℒ and ℒ₀, see Eq. (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 (), 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 Sect. 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.

4.4. 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 Fig. 4 in the case of the baseline model.

Fig. 4. Fermi-LAT imaging centered on Coma and comparison to the model (baseline model, ). We show the Fermi-LAT data (first row), the best-fit model (second row), the residual excluding the cluster model (third row), and the residual with respect to the total model (fourth row). We also show the residual excluding the cluster when accounting for 4FGL J1256.9+2736 (fifth row). Left, middle and right columns: total (200 MeV−300 GeV), low (200 MeV−1 GeV), and high (1 GeV−300 GeV) energy range, respectively. The white contours correspond to TS = [4, 9, 16, 25]. The gray crosses provide the location of the 4FGL sources and the green cross the Coma center.

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 signal 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 Sect. 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° and RA, Dec = 202, 29.5° (TS > 9). Another excess count is seen near RA, Dec = 197.5, 32° (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 × θ₅₀₀ from the cluster center and compare it to the different components of the model in Fig. 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 × θ₅₀₀). 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 × θ₅₀₀ of the Coma cluster.

Fig. 5. Fermi-LAT counts as a function of energy, computed within 3θ500 from the Coma center. The blue line shows the total model in the case of scenario 2 (but it is not distinguishable from scenario 1 in this figure), the red line shows the model when excluding the Coma cluster diffuse emission or 4FGL J1256.9+2736, the orange line show the contribution from 4FGL J1256.9+2736 in the case of scenario 1, and the green line show the Coma cluster diffuse contribution in the case of scenario 2. The contribution from the different sources in the ROI is also indicated as gray lines, except for the isotropic and diffuse backgrounds, given in magenta and purple, respectively. Bottom panel: residual between the data and the model, in red when both the Coma cluster diffuse emission or 4FGL J1256.9+2736 are excluded from the model, in orange in the case of scenario 1 and in green in the case of scenario 2.

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 Fig. 6. Although each of the 5 first bins (up to 1° 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.

Fig. 6. Radial profile of the Fermi-LAT excess count in the total energy band (200 MeV−300 GeV) for different models. The black data points give the excess count when both 4FGL J1256.9+2736 or a cluster model is excluded from the model. We show the best-fit model in the case of scenario 1 (4FGL J1256.9+2736 only, orange line) and in the case of scenario 2 (cluster diffuse model only; baseline in green and point source in blue). The shaded areas correspond to the expected Poisson uncertainties given the respective model. Bottom panel: residual normalized by the error bar with a similar color code.

4.5. 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 component. However, we note that only minor differences are observed when fixing the slope or leaving the background free. In addition to the flux and error bars in each bin, the sed method provides the full likelihood scan for the normalization value in each bin, and we will use this information later in Sect. 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⁻¹⁰ ph s⁻¹ cm⁻² 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⁻¹⁰ ph s⁻¹ cm⁻². 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.

4.6. 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 Fig. 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.

Fig. 7. Light curve associated with the Coma cluster model fit, for the full dataset, in bins of about 4.5 months. For clarity, error bars are shown only for points that are mode than 1σ away from zero, and upper limits (95% confidence interval) are also given. Bottom panel: square root of the TS associated with the source.

4.7. 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 Fig. 5) and is slightly correlated with the cluster signal (see Fig. 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)⁷, 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, 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° per pixel, and eight energy bins per decade. We vary the pixel resolution from 0.05 to 0.2° 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°. 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°. 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° 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°. We reproduce the results presented here without considering this cut. While the statistics slightly increases, the changes are less than 4% on the results.

4.8. On the nature of the signal toward the Coma cluster

In Sect. 4, we have shown that a significant γ-ray signal is observed within the characteristic radius θ₅₀₀ of the Coma cluster. This agrees with the results by Xi et al. (2018) and the 4FGL-DR2 catalog (Abdollahi et al. 2020; Ballet et al. 2020).

The source 4FGL J1256.9+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 component 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.

5. Multiwavelength comparison

In this section, we qualitatively compare the Fermi-LAT excess map obtained in Sect. 4 to data at other wavelengths, as described in Sect. 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 Fig. 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.

Fig. 8. SDSS color image of the central 2 × 2 deg2 cluster region. The thin gray contours provide the galaxy density as shown in Fig. 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.

In Fig. 9, we compare the TS map to images of the tSZ, X-ray, galaxy density, and radio signal. The northwest-southeast 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 associated 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).

Fig. 9. Multiwavelength morphological comparison of the Coma cluster signal to the Fermi-LAT TS map obtained in our baseline model. Top left: Planck tSZ. Top right: ROSAT X-ray. Bottom left: SDSS galaxy density. Bottom right: WSRT 352 MHz radio signal. The field of view of all images is 5 × 5 deg2. The white contours give the Fermi-LAT TS map (contours at 4, 9, 16, and 25) for the reference MINOT model (). For all panels, the black contours correspond to the maximum of the image divided by 2i, 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, Fig. 8 provides an optical image of the central region.

As noted in Fig. 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 X 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 Sect. 4 and in particular in Fig. 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.

6. Implications for the cosmic ray content of the Coma cluster

In this section, we use the γ-ray SED extracted in Sect. 4, together with the model described in Sect. 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 three-dimensional physical model of the cluster, which is needed to constrain the CR content.

6.1. Methodology

We aim at using the Fermi-LAT extracted SED to constrain the CRp population of our model. As discussed in Sect. 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 Sect. 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₅₀₀), which we define at a radius R₅₀₀; (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₅₀₀) ∈ [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

(13)

where i runs over each energy bin and the model parameters are θ ≡ [X_CRp(R₅₀₀),α_CRp]. The value ℒ_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 Sect. 4.5.

Once the chains have converged, and after removing the burn-in phase, the two-dimensional histogram in the plane X_CRp(R₅₀₀)−α_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 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.

6.2. Constraints on the cosmic ray proton population

The SED measured in the case of our baseline model is shown in Fig. 10 for both scenario 2 and 3 (see Sect. 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⁻⁷ MeV cm⁻² s⁻¹. 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.

Fig. 10. Total SED recovered from the Fermi-LAT and MCMC constraint in the case of the reference model () when 4FGL J1256.9+2736 is replaced by the ICM component in the sky model (left) or both are included (right). 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.

The corresponding constraint on the posterior likelihood of the parameters X_CRp(R₅₀₀) and α_CRp is shown in Fig. 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 Fig. 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 constraint on the normalization shifts toward zero and the posterior only excludes X_CRp(R₅₀₀) = 0 at about 2σ. The uncertainty on the slope increases and the best-fit slightly decreases to about α_CRp ≃ 2.6.

Fig. 11. MCMC constraints on the model parameters, in the case of the reference model () when 4FGL J1256.9+2736 is excluded from the sky model (green) or included (orange). Bottom left panel: constraint in the plane XCRp(R500)−αCRp, where the contours correspond to 68% and 95% confidence interval. The marginalized posterior probability distributions are also shown for the parameters XCRp(R500) (top) and αCRp (right), where the shaded area provides the 68% confidence interval.

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%.

6.3. 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 Sect. 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 Sect. 3, the primary CRe, i.e., CRe₁, 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, 2011, for reacceleration models), the population of CRe₁ that we constrain is supposed to be independent of the CRp and the CRe₂. It would correspond, for instance, to a CRe₁ population arising from star formation activity spread over the cluster volume, or the direct shock acceleration in the ICM. In this case, the CRe₁ and CRe₂ populations coexist and the radio emission accounts for the sum of the two. (2) We could also interpret the CRe₁ population in our model as the reaccelerated CRe₂ population. In this case, the CRe₂ would be the seed for the CRe₁. The CRe₂ population should thus not contribute directly to the radio emission. Instead the ratio CRe₁/CRe₂ 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 provide a constraint relative to the hadronic steady state scenario since this is an assumption made when computing the CRe₂ 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₁ are parametrized using the ExponentialCutoffPowerLaw model (our baseline), the InitialInjection model, or the ContinuousInjection model, as described in Sect. 3.1.4. All models include a normalization X_CRe1(R₅₀₀), a spectral slope α_CRe1, a break energy (E_cut, CRe1), and a spatial scaling relative to the thermal density (η_CRe1). We fit for these parameters using simultaneously the radio spectrum and the 352 GHz radio profile data (see Sect. 2). Regarding the spectrum, we compute our model using a cylindrical integration within R = 0.48 × R₅₀₀ (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 Fig. 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₁ and CRe₂ contributions (case 1, no reacceleration), or only the CRe₁ as the reaccelerated CRe₂ (case 2, pure reacceleration). We compute the error contour on the CRe₁ fitted population by reproducing the fit in the case of the lower and upper bounds for the CRe₂ population. We thus assume that the uncertainties from the CRe₂ population (given by the γ-ray) are dominant over the uncertainties associated with the radio data.

Fig. 12. Constraint on the CRe populations in the case of scenario 1 (distinct CRe1 and CRe2 populations, no reacceleration). Top left: radio spectrum of Coma, as compiled from Pizzo (2010) and constraint from the reference CRp spatial model () and the reference CRe1 spectral model (ExponentialCutoffPowerLaw). The measurement from the Brown & Rudnick (2011) data is also shown as the magenta diamond. The contribution from the CRe2 is shown in blue together with its 68% confidence interval, and the remaining contribution from CRe1 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 × R500 = 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 CRe2 and CRe1 are as in the left panel. Bottom left: absolute number density distributions of CRe1 and CRe2 taken at 1 GeV and 10 GeV. Bottom right: ratio between the CRe1 and CRe2 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 Sect. 7.2 for discussions).

In Fig. 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 () 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 models 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₁ spectral model). The slope of the synchrotron emission from CRe₂ is similar to the one from the CRe₁, but it is significantly less curved and no high energy cutoff is present in the CRe₂ distribution. We can see on the spectrum that the radio emission within 0.48 × R₅₀₀ (629 kpc, 0.36 deg) from CRe₂ is overall a factor of about four lower than the total emission (except at high frequency where it reaches similar values). On the profile, the CRe₂ 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₁ contribution. This high concentration is expected because the CRe₂ profile arises as the product of the gas density and the CRp density. The synchrotron profile, which arises from the product of the magnetic field profile 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 Zandanel et al. 2014). Thus, this could be achieved for the CRe₂ only at the cost of an inverted CRp profile (rising with radius). The number density of CRe is at a level of about 10⁻¹⁴ and 10⁻¹⁷ cm⁻³ GeV⁻¹ 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₁ to CRe₂ that increases with radius and that does not depend much on energy up to 10 GeV (before the cutoff; best-fit E_cut, CRe1 = 17 GeV). Given this CRp spatial model, the CRe₁/CRe₂ ratio is slightly below unity in the core, and strongly rises to reach about 100 at 2 Mpc.

In Fig. 13, we show the constraint on the CRe population in the second case (CRe₂ are the seed to the CRe₁, pure reacceleration). As for Fig. 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 the CRe₁ population (as the reaccelerated CRe₂ population) is included. On the bottom panels, the interpretation is now different as the amount of CRe₁ now correspond to the CRe₂ population after reacceleration. As can be observed, the best-fit CRe₁ 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₁ population.

Fig. 13. Same as Fig. 12 in the case of scenario 2 (CRe1 interpreted as the reaccelerated CRe2 population). The CRe2 prior reacceleration refer to the steady state hadronic model. The reacceleration boost is given relative to the steady state hadronic model.

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₅₀₀ 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 Fig. 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 Figs. 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₂ 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₂ synchrotron spectrum is slightly steeper and the amplitude of the CRe₂ component is reduced, the shapes of the CRe₂/CRe₁ 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₂ induced synchrotron spectrum to that of the model developed by Brunetti & Lazarian (2011) in Fig. 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 Sect. 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.

Fig. 14. Comparison between the CRe2 induced synchrotron spectrum to the reacceleration model developed by Brunetti & Lazarian (2011), in the case of a flat CRp population. The solid brown line corresponds to the full reacceleration model, while the dashed brown line corresponds to the case where reacceleration is switched off (see Brunetti & Lazarian 2011, for more details). Left: replacing 4FGL J1256.9+2736 by the cluster diffuse component (scenario 2). Right: including both 4FGL J1256.9+2736 and the cluster diffuse component in the sky model (scenario 3).

7. Discussions

7.1. 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⁻⁹ ph cm⁻² s⁻¹ (E > 100 MeV). Extrapolating our best-fit model in the same energy range, we obtain a compatible total flux of 2.1 × 10⁻⁹ ph cm⁻² s⁻¹ (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⁻⁹ ph cm⁻² s⁻¹ 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⁻⁹ ph cm⁻² s⁻¹. 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°).

7.2. 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₂₀₀ = 2.2 Mpc) = 2.4 × 10⁻² (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., Zandanel et al. 2014). While our results favor models with intermediate scaling (), the data are not sufficient to firmly constrain the shape of the CRp distribution. In the context of the reacceleration of CRe₂, 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, 2017; Zandanel & Ando 2014; Zandanel et al. 2014). 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, 2011).

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 conditions, 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 Sect. 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₁ 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₂ population are nonetheless expected to dominate, and we leave the joint investigation of the CRe population and magnetic field distribution aside.

7.3. 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⁴⁰ and 2 × 10⁴⁰ erg s⁻¹ (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⁴⁰ − 3 × 10⁴² erg s⁻¹ 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.

8. 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° 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 2019) to detect the diffuse emission associated with hadronic interactions in the ICM of galaxy clusters.

Acknowledgments

We are thankful to the anonymous referee for useful comments that helped improve the quality of the paper. We thank Gianfranco Brunetti for useful discussions and for providing us with the reacceleration model used in Fig. 14. Partial support for LR comes from U.S. National Science Foundation grant AST 17-14205 to the University of Minnesota. This research made use of Astropy, a community-developed core Python package for Astronomy (Astropy Collaboration 2013), in addition to NumPy (van der Walt et al. 2011), SciPy (Jones et al. 2001) and Ipython (Pérez & Granger 2007). Figures were generated using Matplotlib (Hunter 2007) and seaborn (Waskom et al. 2020).

References

Appendix A: Inverse Compton emission

This appendix provides an estimate of the inverse Compton emission associated with CRe₂ and CRe₁.

Firstly, we assume a CRp population with X_CRp(R₅₀₀) = 10⁻² and α_CRp = 2.8, which is consistent with the values found in this paper for Coma, and use our baseline scaling . We compute the associated CRe₂ population in the steady state approximation, as in Sect. 6. We finally compute the inverse Compton emission as explained in Adam et al. (2020). This provides an estimate of the necessary inverse Compton emission associated with the hadronic collisions.

Secondly, we assume that the CRe₁ population is described by a power law, which we match to radio data and extrapolate to very high energies. As strong losses are expected at high energy, this provides an upper limit on the CRe₁ population for energies that are large enough to induce inverse Compton emission in the Fermi-LAT energy range (energies in the range of about 100 GeV−100 TeV, while the radio emission probe CRe₁ up to a few tens of GeV at GHz frequencies). We match the slope and spectrum of the power law to the radio spectrum and compute the associated inverse Compton emission.

These two cases (CRe₁ and CRe₂) are reported in Fig. A.1. On the left panel, we see the contribution of the CRe₁ and CRe₂ to the radio signal. On the right panel, we see their contribution to the inverse Compton emission. Even in the case of boosting the CRe₂ population (i.e., the turbulent reacceleration model) by a factor of five to ten to match the radio emission, this would remain more than an order of magnitude below the hadronically induced γ-rays. An energy-dependent boost, strongly rising with the energy, could lead to the inverse Compton signal being significant, but this is not expected. In the case of CRe₁, the upper limit provided by the power law model remains below the hadronically induced γ-rays at energies above 100 MeV. As losses are expected to strongly reduce the signal as energy increases, we expect the emission to be negligible, especially since the Coma radio spectrum already present a curvature at GHz frequencies.

Fig. A.1. Estimation of the inverse Compton emission to the γ-ray signal. Left: radio synchrotron spectrum data to which the CRe population can be matched. Right: associated emission in the γ-ray energy range.

While these results may slightly vary depending on the exact shape of the spatial and spectral distribution of the CR, we do not expect any major change in the comparison. We conclude that the Inverse Compton emission from both CRe₁ and CRe₂ should not contribute significantly to the observed signal at Fermi-LAT energies, and it is neglected when modeling the γ-ray emission.

Appendix B: Comparison to Monte Carlo realizations

In this appendix, we quantify the likelihood of finding a point source 0.68° away from the Coma cluster center (as 4FGL J1256.9+2736), in the hypothesis in which the point source corresponds in fact to the peak of the diffuse cluster ICM emission (i.e., scenario 2). To do so we use an approach that is based on Monte Carlo realizations, as follows. First, we select the best-fit sky model of scenario 2 (with 4FGL J1256.9+2736 replaced by the cluster), compute the number count prediction model, and use it to perform Monte Carlo realizations of the data using poisson sampling. For each realization, we fit the model as done for the real data when excluding the cluster ICM from the sky model (see Sects. 4.2 and 4.3). We then include a point source in the sky model, starting at the coordinates of the Coma reference center, and use the localize method from Fermipy to search for the location which provides the best match for such a test point source within 1.5°. We repeat the operation 200 times and compute the probability to find such a source at distance θ from the cluster.

The results are provided in Fig. B.1, where we show the best-fit coordinates of the test point source for all the Monte Carlo realizations and compare it to the position of the Coma reference center and 4FGL J1256.9+2736. We also display the histogram of the offsets between the recovered test point source and the Coma reference center. The results are shown in the case of our baseline extended ICM model and in the case the spatial ICM template is based on the tSZ map, for comparison. The results for the other cases are expected to slightly vary based on the extension of the spatial template, but are not shown here. As we can see, although it is located 0.68° away from the Coma reference center, 4FGL J1256.9+2736 does not appear as an outlier compared to our Monte Carlo. We find that 4.8% of the recovered sources are located further than 0.68° from the Coma center in the extended model, and 9.5% in the more diffuse tSZ template case. This is due to the combination of the low S/N, the limited Fermi-LAT angular resolution and the extension of the ICM emission model. Indeed, the Fermi-LAT PSF varies with energy and around 1 GeV, where the signal is peaking, the 68% containment angle is about 1 deg (about 3 deg for 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°, in the extended case, see Fig. 1), this agrees with the offset being not significant.

Fig. B.1. Top: sky distribution of the test source coordinates recovered for each Monte Carlo realization and comparison to the Coma reference center and the 4FGL J1256.9+2736 location. Bottom: distribution of the offset from the Coma reference center of the Monte Carlo realizations. Left panels: extended ICM model and right panel: tSZ template model.

Although a diffuse emission plus point source model provides the best match to the data (see Sect. 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 Sect. 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° 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° northwest from the Coma center, power law index of 1.9), 4FGL J1250.8+3117 (TS = 137.29, 3.8° northwest from the Coma center, power law index of 2.16), 4FGL J1316.5+3013 (TS = 20.90, 4.3° northeast from Coma, power law index of 2.1). To match what is done for Coma, we define the cluster center 0.68° 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 Sect. 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₁ in Fig. D.1. Similarly, Fig. D.2 gives the results for the alternative spatial models for the CRp population. The results are also provided in the case of including 4FGL J1256.9+2736 in the Fermi-LAT sky model (i.e., scenario 3).

Fig. D.1. Same as Fig. 13 for the InitialInjection model (left) and the ContinuousInjection model (right).

Fig. D.2. Same as Fig. 12 (rows 1, 2, 3, and 4) and Fig. 13 (rows 6 and 7) for the compact model (first column), the flat model (second column), isobaric (third column) and the extended model + point source (fourth column).

All Tables

All Figures

	Fig. 1. MINOT templates used to describe the cluster, for different assumptions regarding the CRp distributions. All models are normalized so that XCRp(R500) = 10−2. Top left: radial profile of the CRp distribution. Top right: enclosed CRp to thermal energy profile. Bottom left: γ-ray surface brightness profile. Bottom right: γ-ray energy spectrum integrated within R500.
In the text

Fig. 2. Template images used to describe the morphology of the γ-ray emission. From left to right and top to bottom: templates 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 deg2, except for the radio relic (bottom right), which is only 1 × 1 deg2. 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.

In the text

Fig. 3. Correlation matrix associated with the likelihood fit. Each entry corresponds to one free parameter of the model. The names isodiff and 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 .

In the text

Fig. 4. Fermi-LAT imaging centered on Coma and comparison to the model (baseline model, ). We show the Fermi-LAT data (first row), the best-fit model (second row), the residual excluding the cluster model (third row), and the residual with respect to the total model (fourth row). We also show the residual excluding the cluster when accounting for 4FGL J1256.9+2736 (fifth row). Left, middle and right columns: total (200 MeV−300 GeV), low (200 MeV−1 GeV), and high (1 GeV−300 GeV) energy range, respectively. The white contours correspond to TS = [4, 9, 16, 25]. The gray crosses provide the location of the 4FGL sources and the green cross the Coma center.

In the text

Fig. 5. Fermi-LAT counts as a function of energy, computed within 3θ500 from the Coma center. The blue line shows the total model in the case of scenario 2 (but it is not distinguishable from scenario 1 in this figure), the red line shows the model when excluding the Coma cluster diffuse emission or 4FGL J1256.9+2736, the orange line show the contribution from 4FGL J1256.9+2736 in the case of scenario 1, and the green line show the Coma cluster diffuse contribution in the case of scenario 2. The contribution from the different sources in the ROI is also indicated as gray lines, except for the isotropic and diffuse backgrounds, given in magenta and purple, respectively. Bottom panel: residual between the data and the model, in red when both the Coma cluster diffuse emission or 4FGL J1256.9+2736 are excluded from the model, in orange in the case of scenario 1 and in green in the case of scenario 2.

In the text

Fig. 6. Radial profile of the Fermi-LAT excess count in the total energy band (200 MeV−300 GeV) for different models. The black data points give the excess count when both 4FGL J1256.9+2736 or a cluster model is excluded from the model. We show the best-fit model in the case of scenario 1 (4FGL J1256.9+2736 only, orange line) and in the case of scenario 2 (cluster diffuse model only; baseline in green and point source in blue). The shaded areas correspond to the expected Poisson uncertainties given the respective model. Bottom panel: residual normalized by the error bar with a similar color code.

In the text

	Fig. 7. Light curve associated with the Coma cluster model fit, for the full dataset, in bins of about 4.5 months. For clarity, error bars are shown only for points that are mode than 1σ away from zero, and upper limits (95% confidence interval) are also given. Bottom panel: square root of the TS associated with the source.
In the text

	Fig. 8. SDSS color image of the central 2 × 2 deg2 cluster region. The thin gray contours provide the galaxy density as shown in Fig. 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.
In the text

Fig. 9. Multiwavelength morphological comparison of the Coma cluster signal to the Fermi-LAT TS map obtained in our baseline model. Top left: Planck tSZ. Top right: ROSAT X-ray. Bottom left: SDSS galaxy density. Bottom right: WSRT 352 MHz radio signal. The field of view of all images is 5 × 5 deg2. The white contours give the Fermi-LAT TS map (contours at 4, 9, 16, and 25) for the reference MINOT model (). For all panels, the black contours correspond to the maximum of the image divided by 2i, 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, Fig. 8 provides an optical image of the central region.

In the text

Fig. 10. Total SED recovered from the Fermi-LAT and MCMC constraint in the case of the reference model () when 4FGL J1256.9+2736 is replaced by the ICM component in the sky model (left) or both are included (right). 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.

In the text

Fig. 11. MCMC constraints on the model parameters, in the case of the reference model () when 4FGL J1256.9+2736 is excluded from the sky model (green) or included (orange). Bottom left panel: constraint in the plane XCRp(R500)−αCRp, where the contours correspond to 68% and 95% confidence interval. The marginalized posterior probability distributions are also shown for the parameters XCRp(R500) (top) and αCRp (right), where the shaded area provides the 68% confidence interval.

In the text

Fig. 12. Constraint on the CRe populations in the case of scenario 1 (distinct CRe1 and CRe2 populations, no reacceleration). Top left: radio spectrum of Coma, as compiled from Pizzo (2010) and constraint from the reference CRp spatial model () and the reference CRe1 spectral model (ExponentialCutoffPowerLaw). The measurement from the Brown & Rudnick (2011) data is also shown as the magenta diamond. The contribution from the CRe2 is shown in blue together with its 68% confidence interval, and the remaining contribution from CRe1 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 × R500 = 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 CRe2 and CRe1 are as in the left panel. Bottom left: absolute number density distributions of CRe1 and CRe2 taken at 1 GeV and 10 GeV. Bottom right: ratio between the CRe1 and CRe2 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 Sect. 7.2 for discussions).

In the text

	Fig. 13. Same as Fig. 12 in the case of scenario 2 (CRe1 interpreted as the reaccelerated CRe2 population). The CRe2 prior reacceleration refer to the steady state hadronic model. The reacceleration boost is given relative to the steady state hadronic model.
In the text

Fig. 14. Comparison between the CRe2 induced synchrotron spectrum to the reacceleration model developed by Brunetti & Lazarian (2011), in the case of a flat CRp population. The solid brown line corresponds to the full reacceleration model, while the dashed brown line corresponds to the case where reacceleration is switched off (see Brunetti & Lazarian 2011, for more details). Left: replacing 4FGL J1256.9+2736 by the cluster diffuse component (scenario 2). Right: including both 4FGL J1256.9+2736 and the cluster diffuse component in the sky model (scenario 3).

In the text

	Fig. A.1. Estimation of the inverse Compton emission to the γ-ray signal. Left: radio synchrotron spectrum data to which the CRe population can be matched. Right: associated emission in the γ-ray energy range.
In the text

	Fig. B.1. Top: sky distribution of the test source coordinates recovered for each Monte Carlo realization and comparison to the Coma reference center and the 4FGL J1256.9+2736 location. Bottom: distribution of the offset from the Coma reference center of the Monte Carlo realizations. Left panels: extended ICM model and right panel: tSZ template model.
In the text

	Fig. D.1. Same as Fig. 13 for the InitialInjection model (left) and the ContinuousInjection model (right).
In the text

	Fig. D.2. Same as Fig. 12 (rows 1, 2, 3, and 4) and Fig. 13 (rows 6 and 7) for the compact model (first column), the flat model (second column), isobaric (third column) and the extended model + point source (fourth column).
In the text