© R. Adam et al. 2020

1. Introduction

Galaxy clusters are the largest gravitationally bound structures that are decoupled from the expansion of the Universe. They form peaks in the matter density field. Their assembly has been driven by the gravitational collapse of dark matter (Kravtsov & Borgani 2012), which is thought to dominate the matter content of clusters (about 80% in mass). Clusters also consist of baryonic matter, essentially in the form of hot ionized thermal plasma, called the intracluster medium (ICM; about 15%), and of galaxies (about 5%). While clusters are used to understand the formation of large-scale structures and to constrain cosmological models, they are also the place of very rich astrophysical processes and excellent targets for testing fundamental physics (see, e.g., Allen et al. 2011, for a review).

Galaxy clusters form through the merging and accretion of other groups and surrounding material (Sarazin 2002). This leads to the propagation of shocks and turbulences in the ICM (Markevitch & Vikhlinin 2007), which can accelerate charged particles to very high energies. These cosmic rays (CRs) interact with the magnetized ICM, generating diffuse radio synchrotron emission (Feretti et al. 2012; van Weeren et al. 2019), and they are expected to produce a γ-ray signal because of the inverse Compton interaction with background light or the decay of pions produced in proton-proton collisions (Brunetti & Jones 2014). In addition, clusters also host active galactic nuclei (AGN), which are known to provide feedback onto the ICM (Fabian 2012). This feedback is only poorly understood, but is expected to have a major effect on the formation and the evolution of galaxy clusters.

Diffuse radio synchrotron emission in galaxy clusters is generally classified as radio halos (including giant and mini-halos), radio relics, and revived AGN fossil plasma source (see van Weeren et al. 2019, for detailed discussions). While relics are thought to be associated with the shock acceleration of electrons in the periphery of clusters, radio halos might originate from turbulent reacceleration of seed electrons and/or secondary electrons produced by hadronic interactions. γ-ray emission is also expected as a result of the inverse Compton emission that arises from the scattering of background photon fields onto relativistic electrons, or the hadronic interaction from CR protons (CRp) and the ICM (see, e.g., Pinzke & Pfrommer 2010, for the signal expected based on numerical simulations).

The annihilation or the decay of dark matter particles might also cause γ-ray emission from galaxy clusters (see, e.g., Combet 2018), and many searches for this signal have been performed (e.g., Ackermann et al. 2010; Aleksić et al. 2010; Arlen et al. 2012; Abramowski et al. 2012; Combet et al. 2012; Cadena 2017; Acciari et al. 2018). In the case of dark matter decay, galaxy clusters are particularly competitive targets because the signal scales linearly to the huge dark matter reservoirs in galaxy clusters. In the case of dark matter annihilation, clusters can be at the same flux level as dwarf galaxies when substructures are accounted for, and they are thus also highly relevant targets (Sánchez-Conde et al. 2011; Moliné et al. 2017). However, the limits that can be set on the properties of dark matter depend on the uncertainties associated with the modeling of the background emission, so that accurate CR modeling is also essential for dark matter searches.

Many attempts to detect the cluster γ-ray emission have been made using ground-based (e.g., H.E.S.S. Collaboration 2009; Aleksić et al. 2012; Arlen et al. 2012; Ahnen et al. 2016, at 50 GeV−10 TeV energies) and space-based observations (e.g., Reimer et al. 2003; Huber et al. 2013; Prokhorov & Churazov 2014; Zandanel & Ando 2014; Ackermann et al. 2014, 2015, 2016, at about 30 MeV−300 GeV). While unsuccessful so far, these searches were very useful to constrain the CR physics and particle acceleration at play in clusters, especially when combined with radio observations (e.g., Vazza et al. 2015; Brunetti et al. 2017). Recently, Xi et al. (2018) claimed the first significant detection of γ-ray signal toward the Coma cluster using data obtained with the Fermi-Large Area Telescope (Fermi-LAT). However, their results might be confounded by a possible point source because the signal-to-noise ratio and angular resolution of the observations are limited. While the Fermi-LAT satellite is to continue to take data for several additional years (compared to about the 12 years of data collected so far), major discoveries concerning galaxy clusters are unlikely given the modest increase in statistics that is expected. From the ground, the Cherenkov Telescope Array (CTA, Cherenkov Telescope Array Consortium 2019) is expected to provide a major improvement in sensitivity in the 100 GeV−100 TeV energy range.

Fig. 1. Overview of the MINOT input modeling, the considered physical processes at play in the ICM, and the observables that are computed. The interdependences are shown by the black arrows.

Fig. 2. Overview of the code structure and interfaces.

In order to address the CR physics in galaxy clusters, multiwavelength observations and analyses are becoming particularly relevant with the construction of such new facilities. While cluster CRs can essentially be accessed the radio and γ-ray bands, their physics is driven by the continuous interaction with the thermal plasma. When data are compared to modeling, or when mock observations are generated, the thermal and the nonthermal components should therefore be modeled together, in a self-consistent way, so that uncertainties and degeneracies between the two can be accounted for. The thermal emission can be probed in particular in the X-ray and at millimeter wavelengths through thermal Bremsstrahlung emission (Sarazin 1986; Böhringer & Werner 2010) and the thermal Sunyaev-Zel’dovich (tSZ) effect (Sunyaev & Zeldovich 1970, 1972). In addition to the primary components, modeling the particle interactions in the ICM relies on particle physics data from accelerators or a sophisticated Monte Carlo code (see, e.g., Kafexhiu et al. 2014, for discussions), and they need to be accounted for carefully. Nevertheless, clusters are commonly modeled focusing on individual (or just a few) components, and no public self-consistent multiwavelength software exists in the literature. For instance, Li et al. (2019) and Brüggen & Vazza (2020) recently modeled the diffuse radio synchrotron emission of radio halos and radio relics, respectively, and employed the Press-Schechter formalism to estimate the statistical properties of the corresponding signal.

Here, we present a software dedicated to the self-consistent modeling of the thermal and nonthermal diffuse components of galaxy clusters, for which the main objective is computing accurate and well-characterized multiwavelength predictions for the radio, millimeter, X-ray, γ-ray, and neutrino emission. This software is called MINOT, modeling the intracluster medium (non-)thermal content and observable prediction tools. It is based on the Python language and is available online¹. MINOT includes various parameterizations for the radial profiles and spectral properties of the different cluster components. The code does not aim at computing the CR production rate from microphysics considerations (e.g., turbulence, shocks, or diffusion), but instead directly models the spatial and spectral distributions of the CRs and the thermal gas. The predictions for associated observables are available in the radio (synchrotron), millimeter (tSZ effect), X-ray (thermal Bremsstrahlung), γ-ray (inverse Compton and hadronic processes), and also for neutrino emissions (hadronic processes). This includes surface brightness profiles or maps, spectra, and integrated flux computed with different options. For γ-rays, CR electrons (CRe), and neutrinos from hadronic origin, MINOT includes the latest description of the hadronic interactions in the ICM, based on the Naima software (Zabalza 2015). The thermal modeling uses the XSPEC software for X-ray predictions (Arnaud 1996), and it includes an accurate description of the tSZ signal up to high plasma temperatures.

This article is organized as follows. In Sect. 2 we provide a general overview of the code and discuss the different interfaces. Section 3 discusses the physical modeling of the cluster components. The physical processes related to particle interactions are detailed in Sect. 4. In Sect. 5 we discuss the prediction of observables in the relevant energy bands. The use of MINOT is illustrated in Sect. 6 for three nearby massive well-known clusters for which multiwavelength data are available in the literature. Finally, Sect. 7 provides a summary and conclusion. Equations are given following the international system of units.

2. General overview and structure of the code

MINOT is a Python-based code². It essentially depends on standard Python libraries, but some functionalities require specific softwares and packages, as discussed below. In this section, we provide a general overview of the working principle of the code, of its structure, and the interactions between the different modules. The list of the code parameters is also discussed, as well as the available functional forms for the radial and spectral models. Figure 1 highlights how the input modeling is used to generate observables via the different plasma processes considered in MINOT, and the general overview of the code is illustrated in Fig. 2.

2.1. Overview of the physical modeling

MINOT was first developed to compute an accurate γ-ray prediction for galaxy clusters. As discussed below and shown in Fig. 1, this requires several key ingredients. Because the same ingredients also provide diagnoses for other observables at different wavelengths through various physical processes, MINOT was further developed to account for them. This allowed us to provide external constraints to a given input modeling that is used to generate γ-ray observables, but also to provide further diagnosis of the physical state of the cluster.

First, the spatial and spectral distributions of primary CRe (CRe₁) and protons are crucial. They generate γ-rays through inverse Compton scattering on the cosmic microwave background (CMB), or through hadronic interactions, respectively. Modeling the thermal gas is also essential because hadronic processes arise from the interaction between CRp and thermal plasma. As we show in Sect. 3, the thermal component is based on the thermal electron pressure and density. Additionally, the normalization of the CR distributions is generally given relative to the thermal energy. The hadronic interactions also generate secondary CRe and positrons (CRe₂). Because they are affected by synchrotron losses, they require that the magnetic field is accounted for as another key ingredient of the input modeling (Sect. 4). These electrons contribute to the inverse Compton emission (Sect. 5). In summary, the necessary input modeling ingredients are CRe₁, magnetic field strength, CRp, and thermal electron pressure and density.

With these ingredients at hand, the radio emission that arises from CRe (primary and secondary) moving in the magnetic field can be modeled. Similarly, the thermal pressure and density allow us to compute the tSZ signal and the thermal Bremsstrahlung X-ray emission. They also provide a complete diagnosis of the thermodynamic properties of the cluster. Neutrinos are also produced during hadronic interactions, and their associated observable is thus available.

2.2. Code structure

In order to model galaxy clusters and predict observables associated with the diffuse thermal and nonthermal components, MINOT is organized in six main parts, each of which is related to specific functions and procedures. We list these parts below.

1. The main class, called Cluster, is written in the file model.py and provides an entry point for the user. It allows defining the model and solve entanglement between parameters.

2. A subclass called Admin allows us to handle administrative tasks, in particular, input and output procedures.

3. The Modpar subclass is dedicated to the model parameterization. It gathers a library of available radial and spectral models.

4. The physical modeling of the cluster is performed in the Physics subclass. It includes many functions for retrieving the desired physical quantities.

5. The Observable subclass allows us to extract the requested cluster observables based on the inner physics encoded in the model.

6. Finally, the subclass Plots is designed for automated plots to provide a cluster diagnostic based on the current modeling.

In addition to these six main parts, MINOT also includes a library called ClusterTools in which numerical tools and astrometric tools are defined. It also includes several classes that are used to compute various physical processes relevant for MINOT. In the following subsections, we describe the working principle of the different functionalities.

2.3. Initialization and parameters

As illustrated in Fig. 2, the user can directly define a cluster object by calling the Cluster (main) class of MINOT,

cluster = minot.Cluster(optional parameters).

Optional parameters, such as cluster name, coordinates, or redshift, can be passed directly to the initialization call. However, any parameter can be modified on the fly, for example,

cluster.redshift = 0.1

The entanglement of parameters is solved in the code. For instance, changing the cluster redshift will automatically change the angular diameter distance of the cluster according to the current cosmological model. Information is provided to the user when the “silent” parameter is set to “False”. The list of parameters that describe the cluster object is available in Table 1. We note that whenever possible, the code uses Astropy units for quantities³.

We can distinguish four types of parameters, as listed in Table 1. The first type corresponds to administration-like parameters (e.g., output directory used to generate products), the second type concerns the global properties of the cluster object (e.g., the redshift), the third type is related to the radial and spectral modeling of the physical quantities of interest of the cluster (e.g., CR number density profile), and the last type allows the user to sample the output observables. In particular, it is possible to set a map header (e.g., obtained for real data) on which the model prediction maps are projected, which facilitates a comparison of data and model.

2.4. Modeling the physical state of the cluster

The parameters describing the cluster can be divided into two types: global properties that apply to the entire cluster (e.g., mass, redshift, and coordinates), and properties that vary as a function of radius or energy. This separation is highlighted in Table 1. Some parameters are assumed to be constant over the entire cluster volume, such as the hydrostatic mass bias or the metal abundances.

In addition to the global properties, the primary quantities that are used to define the physical state of the cluster are (see also Sect. 3 for further details) the gas pressure of thermal electrons, the gas number density of thermal electrons, the CRp number density profile and spectrum, the CRe₁ profile and spectrum, and the profile of the magnetic field strength. The CR distributions are normalized according to the ratio of CR and thermal energy enclosed within a given radius. The physical modeling of the radial and spectral properties of the cluster relies on a library of predefined models in the Modpar subclass of MINOT. The list of models that are currently available in the code is given in Tables 2 and 3 for the spectral and spatial component, respectively. Figure 3 illustrates the shape of the different models; they are further discussed in a more physical context in Sect. 3. We note that the spatial and spectral parts of the modeling are currently decoupled (e.g., the spectrum of CRp does not change with radius), such that a physical quantity f can be expressed as

$$\begin{array}{c}\hfill f(r,E)\propto f_1(r)f_2(E),\end{array}$$(1)

Fig. 3. Left: illustration of the different radial profiles available in the library. The β-model parameters are (n0, rc, β) = (1, 50 kpc, 0.7); the SVM model parameters are (n0, rc, β, rs, γ, ϵ, α) = (1, 50 kpc, 0.7, 1500 kpc, 3, 4, 0.5); the GNFW model parameters are (P0, rp, a, b, c) = (1, 200 kpc, 1.5, 4, 0.3); and the double β-model parameters are (n01, rc1, β1, (n02, rc2, β2) = (1, 50 kpc, 0.7, 5, 5 kpc, 2). Right: illustration of the different spectral models available in the library. The index is set to α = 2.3 for all models, and the cutoff or break energy to Ebreak/cut = 1 GeV. In both cases, the “User” model is not shown, but allows passing any arbitrary function that is interpolated in the code.

where E is the particle energy (only relevant for the CRs), and r is the physical radius in three dimensions. However, functions that couple the radius and the energy dependence are ready to be implemented in the model library because any calculations relying on the modeling of f(r, E) are made on 2D grids (energy versus radius) that are ignorant of the underlying parameterization of the distributions. In addition, it is possible to apply some losses, assuming a given scenario, to the input distribution. In this case, f(r, E) is considered as an injection rate, and the output distribution is affected differently for different energy and radii. The implementation of the losses is discussed in detail in Sect. 4.2.

A new model is set to a given physical property by passing a Python dictionary, such as

cluster.density_gas_model =
----- {``name'':``beta'',
----- ``n_0'':1e-3*u.cm**-3,
----- ``beta'':0.7, ``r_c'':300*u.kpc}

or

cluster.spectrum_crp_model =
----- {``name'':``PowerLaw'',
----- ``Index'': 2.5}

It is also possible to automatically set a parameterization of several quantities to predefined physical states without directly setting the model parameters, for instance, forcing the CRp to follow the radial distribution of the gas density,

cluster.set_density_crp_isodens_scal_param()

or to define the thermal electron number density based on the thermal electron pressure in the case of an isothermal cluster with a given temperature. These functions are written as part of the Modpar subclass.

2.5. Derived physical properties and observables

When the desired physical properties of the cluster are set, functions related to the physical description of the cluster, from the subclass Physics, can be called to extract the thermodynamic and CR properties of the cluster, or the production rate of nonthermal particles (the physical modeling is further detailed in Sect. 3), for example, extracting the hydrostatic mass profile, or the neutrino emission rate via

r, M = cluster.get_hse_mass_profile()
dN_dEdVdt = cluster.get_rate_neutrino().

The user may also generate observables corresponding to the radio synchrotron emission, the tSZ signal in the millimeter, the thermal Bremsstrahlung in the X-ray, the inverse Compton emission, and the hadronic emission in the γ-rays, or the associated neutrino emission (see also Sect. 5 for more details). This is implemented in the subclass Observable of MINOT, for instance,

E, dN_dEdSdt = cluster.get_gamma_spectrum().

We note that many of the functions or subclasses related to the physical processes that are called when a function like this is run are gathered in the ClusterTools subdirectory. This also includes many numerical tools.

2.6. Administrative functions

When it is defined, the cluster object also includes various administrative functions, gathered in the subclass Admin. They can be used to display the current values of the parameters, to save the current status of the cluster object (or load a previously saved model), to generate output observable products automatically (maps, profiles, and spectra), or to get the header of the current map. The generation of automatic plots corresponding to the various observables included in MINOT is also available using the subclass Plots.

2.7. Baseline model and test clusters

In the following sections, we use different cluster models to illustrate the behavior of the MINOT code. First, we define a baseline cluster model using parametric functions in order to show the effect of changes in the modeling on the observables. The baseline properties of the cluster were set using a generalized Navarro-Frenk-White profile (GNFW, Nagai et al. 2007) thermal electron pressure profile with (P₀, c₅₀₀, a, b, c) = (2.2 × 10⁻² keV cm⁻³, 3.2, 1, 5, 3.1, 0.0), and an SVM thermal electron number density profile with (n₀, r_c, β, α, r_s, γ, ϵ) = (3 × 10⁻³ cm⁻³, 290 kpc, 0.6, 0.0, 1000 kpc, 3, 1.7). This corresponds to a typical massive merging cluster. The redshift was set to z = 0.02 and the mass to M₅₀₀ = 7 × 10¹⁴ M_⊙, inspired by the Coma cluster. The CRp followed an exponential cutoff power-law spectrum, with spectral index 2.4 and a cutoff energy of 100 PeV. The normalization was set to have a CRp-to-thermal energy ratio within R₅₀₀ of 10⁻², which corresponds to the typical expected values, according to Pinzke & Pfrommer (2010). The CRe₁ followed a continuous injection spectrum, with an injection spectral index 2.3 and a break energy of 5 GeV. The normalization was set to have a CRe₁ -to-thermal energy ratio within R₅₀₀ of 10⁻⁵, that is, about the proton value scaled by the proton-to-electron mass ratio (i.e., a similar Lorentz factor distribution was assumed for the two). The spatial profiles of both CRp and CRe₁ were set to the same shape as the thermal gas density (see Sect. 3.2). The magnetic field profile was set to follow the square root of the thermal gas density and was normalized to have an amplitude of 5 μG, assuming similar properties as those measured for the Coma cluster (Bonafede et al. 2010). This model, referred to as Baseline in the following, was varied whenever the effect of relevant quantities to the physical state of the cluster or observable are illustrated.

In addition to this baseline model, it is also useful to use real clusters, with the aim of comparing our model predictions to measurements available in the literature. To do so, we need clusters whose thermal properties have been measured and are available, over a wide range of spatial scales, in order to calibrate our model as well as possible. The sources targeted by the XMM Cluster Outskirt project⁴ are perfectly suited for this purpose because the project allowed for the precise measurement of the thermal pressure and density profiles of nearby galaxy clusters from about 10 kpc to the cluster outskirts based on XMM-Newton and Planck data (Tchernin et al. 2016; Eckert et al. 2017; Ghirardini et al. 2019). Because we are interested in the nonthermal component of the ICM, we selected clusters from the 12 XCOP for which a diffuse radio halo has been observed, using the GalaxyCluster database⁵, and for which Fermi-LAT constraints have been obtained by Ackermann et al. (2014). We found three objects: (1) Abell 1795, a relaxed cool-core system; (2) Abell 2142, an elongated, dynamically active cluster with a cool-core; and (3) Abell 2255, a merging cluster with a highly perturbed core (see also the recent work by Botteon et al. 2020). We thus note that these objects also present the advantage of sampling different dynamical states that are generally observed in clusters. In order to model these three clusters, we extrapolated the precisely measured density profile of the thermal plasma with a high-order polynomial function. The pressure profiles were fit with a GNFW model, providing a good extrapolation to the data. The nonthermal properties were set following what we did for our baseline cluster model, except that the magnetic field was normalized to 5 μG at 100 kpc. The redshift, mass, and coordinates of the cluster were taken from the XCOP data. Table 4 summarize the properties of these clusters.

3. Physical modeling of the primary components

In this section, we discuss the physical modeling of the cluster. First, the global cluster properties are briefly discussed, as well as several assumptions employed in the modeling. Then the properties of the thermal and nonthermal components are detailed. The cluster modeling relies on primary base physical quantities from which other cluster properties can be derived, in particular in the case of the thermal component. The choice of the base quantities is discussed. Then, the derivation of the secondary quantities that characterize the cluster are developed both for the thermal and nonthermal components.

3.1. Global properties and assumptions

Before we model the inner structure of the clusters, it is useful to characterize the global cluster properties, as listed in the second block of parameters in Table 1. The cluster location is defined in terms of redshift and sky coordinates. From the redshift, and given a cosmological model, the angular diameter, and luminosity distances are computed and used later in the code. The default cosmological model is based on Planck Collaboration XIII (2016), but can be modified if necessary.

Even if it does not play a direct role in the modeling, the characteristic mass of the cluster M₅₀₀ ⁶ is part of the global parameters. It can be used to set several internal properties of the cluster, to their universal expectation, according to the fact that clusters are at first order self-similar objects (in particular for the thermal pressure, see Arnaud et al. 2010). The value of M₅₀₀ also allows us to set the characteristic radius, R₅₀₀ (see also Eq. (14)).

It is also worth emphasizing one of the global parameters, the truncation radius, which is used in MINOT in order to set a physical boundary to the cluster, beyond which the density drops to zero. This is not only useful for numerical issues when the cluster properties are integrated, but might be associated with the accretion shock radius at which the kinetic energy from accreting structures is converted into thermal energy (see, e.g., Hurier et al. 2019, for the observation of such an accretion shock).

In the modeling, the plasma is assumed to be fully ionized and to follow the ideal gas law. The ions and electrons are assumed to be in thermal equilibrium (see Fox & Loeb 1997, for discussions of the electron and ion temperatures). While it might in principle depend on radius, the hydrostatic mass bias, the helium mass fraction, and the metallicity of the cluster are assumed to be constant (see, e.g., Leccardi & Molendi 2008; Nelson et al. 2014, for measured cluster metallicity profiles and simulations of the nonthermal radial profile). Some of these parameters, related to the global properties of the cluster, are further discussed in the following subsections.

3.2. Thermal component

The base thermal properties are the electron number density and the electron pressure profiles. This choice is motivated by the fact the X-ray emission is directly sensitive to the electron number density, while the tSZ effect probes the electron pressure, but other choices could have been made (e.g., density and temperature). Generic literature parametric models are available, such as the β-model (Cavaliere & Fusco-Femiano 1978) or one of its extensions, the simplified Vikhlinin model (SVM, Vikhlinin et al. 2006), which is generally used to describe thermal density cluster profiles. The GNFW profile is also available, and is generally used to describe the thermal pressure profile (Arnaud et al. 2010; Planck Collaboration V 2013). See Table 3 for the parameterization of these models. In all cases, the different parameters can be used to control the amplitude, the characteristic or transition radius, and the slopes at different radii. In Fig. 3 we show examples of these profiles for a given set of parameters.

Fig. 4. Thermodynamic properties of our cluster sample defined in Table 4, as computed based on the modeling of the thermal electron pressure and density (i.e., the base quantities used to derive the others). The truncation radius is shown at r = 5 Mpc, where all distributions drop, except for the integrated ones that remain constant afterward. For the entropy profile, we also show the expectation for purely gravitational collapse (Voit 2005). For the gas fraction, the ratio between the mean baryon density and the mean matter density is shown as a dashed line. We also report the radius at which a density contrast of 500 (i.e., R500) is reached in the overdensity profile.

With the electron pressure, P_e(r), and electron number density, n_e(r), profiles at hand, it is possible to compute the total gas pressure

$$\begin{array}{c}\hfill P_{\mathrm{gas}}(r)=\frac{{\mu }_{\mathrm{e}}}{{\mu }_{\mathrm{gas}}}{P}_{\mathrm{e}}(r),\end{array}$$(2)

and the thermal proton number density profile as

$$\begin{array}{c}\hfill n_{\mathrm{p}}(r)=\frac{{\mu }_{\mathrm{e}}}{{\mu }_{\mathrm{p}}}{n}_{\mathrm{e}}(r).\end{array}$$(3)

The mean molecular weights, μ_gas, μ_e, μ_p, and μ_He, are computed from the helium primordial abundance and the ICM metallicity, as

$$\begin{array}{c}\hfill \{\begin{array}{cc}& {\mu }_{\mathrm{gas}}=\frac{1}{2(1-Y-Z)+\frac{3}{4}Y+\frac{1}{2}Z}\simeq 0.6\hfill \\ \hfill & {\mu }_{\mathrm{e}}=\frac{1}{(1-\frac{1}{2}Y-\frac{1}{2}Z)}\simeq 1.15\hfill \\ \hfill & {\mu }_{\mathrm{p}}=\frac{1}{1-Y-Z}\simeq 1.35\hfill \\ \hfill & {\mu }_{\mathrm{He}}=\frac{4}{Y}\simeq 14.6,\hfill \end{array}\end{array}$$(4)

where Y ≃ 0.27 is the helium mass fraction and Z ≃ 0.005 is the heavy element mass fraction (defined through the solar reference metallicity multiplied by the metal abundance, see Table 1). Here, we used the approximation that $\frac{N_{\mathrm{charge}}+1}{N_{\mathrm{nucleon}}}\simeq 1/2$ for all metals, with N_charge the number of charge and N_nucleon the number of nucleons.

It is also straightforward to compute the temperature assuming the ideal gas law as

$$\begin{array}{c}\hfill k_{\mathrm{B}}{T}_{\mathrm{e}}(r)=P_{\mathrm{e}}(r)/n_{\mathrm{e}}(r)\equiv k_{\mathrm{B}}{T}_{\mathrm{gas}}(r).\end{array}$$(5)

Similarly, the electron entropy index, which records the thermal history of the cluster (Voit 2005), can be defined as

$$\begin{array}{c}\hfill K_{\mathrm{e}}(r)=\frac{P_{\mathrm{e}}(r)}{n_{\mathrm{e}}{(r)}^{5/3}}·\end{array}$$(6)

Temperature and entropy are useful diagnostics of the ICM. They can show the presence of a cool core (e.g., Cavagnolo et al. 2009), which is itself related to the central AGN activity and possibly to its CR feedback onto the surrounding gas (e.g., Ruszkowski et al. 2017). They provide information on the dynamical state and accretion history of the cluster, which are connected to its CR content (e.g., radio emission that begins during mergers, Rossetti et al. 2011).

The thermal energy density stored in the gas is given by

$$\begin{array}{c}\hfill u_{\mathrm{th}}=\frac{3}{2}{n}_{\mathrm{gas}}{k}_{\mathrm{B}}T=\frac{3}{2}{P}_{\mathrm{gas}}\end{array}$$(7)

and can be integrated over the volume

$$\begin{array}{c}\hfill U_{\mathrm{th}}(R)=4\pi {\displaystyle {\int }_0^R}{u}_{\mathrm{th}}{r}^2\mathrm{d}r\end{array}$$(8)

to obtain the total thermal energy up to radius R. This quantity is very useful when it is compared to the amount of energy stored in the CRs.

The cluster total mass within radius r, under the approximation of hydrostatic equilibrium, is given by

$$\begin{array}{c}\hfill M_{\mathrm{HSE}}(r)=-\frac{r^2}{{\mu }_{\mathrm{gas}}{m}_{\mathrm{p}}{n}_{\mathrm{e}}(r)G}\frac{\mathrm{d}{P}_{\mathrm{e}}(r)}{\mathrm{d}r}.\end{array}$$(9)

The hydrostatic mass is known to be biased with respect to the actual total mass (see, e.g., Pratt et al. 2019, for a review on the cluster mass scale). The two can be related by

$$\begin{array}{c}\hfill M_{\mathrm{tot}}(r)=\frac{M_{\mathrm{HSE}}(r)}{(1-b_{\mathrm{HSE}})},\end{array}$$(10)

where b_HSE is the hydrostatic mass bias (see Table 1), which is assumed to be constant (see Planck Collaboration XX 2014, for detailed discussions of the bias value, which is expected to be b_HSE ∼ 0.2).

Based on the electron number density, the gas mass within radius R can be computed as

$$\begin{array}{c}\hfill M_{\mathrm{gas}}(R)=4\pi {\displaystyle {\int }_0^R}{\mu }_{\mathrm{e}}{m}_{\mathrm{p}}{n}_{\mathrm{e}}(r)r^2\mathrm{d}r,\end{array}$$(11)

and it provides a measurement of the available target mass for the interaction with CRp. The gas fraction can be derived using

$$\begin{array}{c}\hfill f_{\mathrm{gas}}(r)=\frac{M_{\mathrm{gas}}(r)}{M_{\mathrm{tot}}(r)}·\end{array}$$(12)

The overdensity profile is computed using ρ_c(z), the critical density of the Universe, by

$$\begin{array}{c}\hfill \mathrm{\Delta }(R)=\frac{M_{\mathrm{tot}}(R)}{\frac{4}{3}\pi R^3{\rho }_{\mathrm{c}}(z)},\end{array}$$(13)

which allows us to extract the value of the characteristic radius, R_Δ, within which the density of the cluster is Δ times the critical density of the Universe at the cluster redshift. The value of Δ is generally taken to be 500. The enclosed mass within R_Δ is then

$$\begin{array}{c}\hfill M_{\mathrm{\Delta }}=\frac{4\pi }{3}\mathrm{\Delta }{\rho }_{\mathrm{ref}}(z)R_{\mathrm{\Delta }}^3.\end{array}$$(14)

All the quantities defined here can be extracted as a function of radius from MINOT, according to a given cluster model, using the dedicated functions that are located in the Physics subclass. In Fig. 4 we illustrate the main thermodynamic properties of our baseline cluster model and the three Abell cluster models discussed in Sect. 2.7. We show the thermal electron pressure, electron number density, gas temperature, entropy, enclosed thermal energy, enclosed gas mass, enclosed hydrostatic mass, gas fraction, and overdensity contrast.

Depending on the dynamical state and the presence or absence of a cool core, the profiles are different. For instance, A1795 clearly presents a high-density cool core according to its temperature and entropy profiles. Its large-scale electron pressure and density fall quickly, consistent with a compact, relaxed morphology. In contrast, A2255 and our Baseline model show disturbed cores with a high entropy floor, and their pressure and density profiles are much flatter on large scales, consistent with a redistribution of the thermal energy in a merging event. A2142 is an intermediate case. It presents a peaked density profile (showing a compact core), but its pressure profile is relatively flat on large scales, typical of disturbed clusters. Because particle acceleration is expected to depend on the cluster dynamical state, these thermodynamic diagnosis are useful for characterizing individual clusters in the context of understanding cluster CR physics.

The enclosed thermal energy is directly related to the pressure profile. It is particularly relevant here because it provides a normalization for the number of CRs (see Sect. 3.3). Based on its high pressure profile and thus thermal energy, and based on its high density (which implies a high gas mass), A2142 would thus be the best target to search for γ-rays from proton-proton interaction in our sample, assuming the same CR distribution for all clusters.

The hydrostatic mass provides a direct way to measure the cluster total mass profile (given a hydrostatic bias, Eq. (10)). This can be particularly relevant for modeling the γ-ray signal associated with the decay or annihilation of dark matter particles (which is usually done assuming an NFW dark matter density profile with a given concentration and normalization). In addition, the gas fraction gives the ratio between the amount of dark matter and the amount of gas, which is expected to provide a proxy for the dark matter signal to CR background when indirect dark matter searches using clusters are performed, and thus they are an indication for the best targets. The overdensity contrast allows us to measure the radius R₅₀₀ (or R₂₀₀), and thus the corresponding mass.

For further discussions of the thermodynamic properties of galaxy clusters using a similar framework, we refer to the recent work by Tchernin et al. (2016), Ruppin et al. (2018), Ghirardini et al. (2019), Ricci et al. (2020).

Fig. 5. CR properties of the baseline cluster model. Left: CR number density profile (i.e., CR distribution integrated between Emin and Emax, as indicated in the legend) for CRe1 and CRp. Right: CR spectrum integrated over the radius up to a maximum radius, as given in the legend, for CRe1 and CRp.

Fig. 6. Magnetic field and CRp-to-thermal energy properties of the baseline cluster model. Left: magnetic field profile obtained by assuming a fixed magnetic field strength of 5 μG at 100 kpc, and a scaling relative to the thermal density as $B\propto n_{\mathrm{e}}^{0.5}$. Right: CRp-to-thermal energy profile for different models of the CRp distribution. The color lines show each real cluster, using the same colors as in the left panel, and with nCRp ∝ ngas. The black lines correspond to variation in the CRp scaling with respect to the thermal gas for the baseline cluster model, as indicated in the legend.

3.3. Nonthermal component

The base nonthermal properties of the cluster are the magnetic field strength, the spectra and profile of the CRp, and the CRe₁. The CRe₁ differ from the CRe₂ because they correspond to a population that was accelerated from the plasma microphysics (e.g., shocks or turbulences, as for CRp), while the secondaries are the product of hadronic interactions (see Sect. 4 for further details). The radial models available for the nonthermal component are the same as for the thermal component (see Table 3). The spectral distributions for currently available models are listed in Table 2. While all these models might be attributed to either the CRp or the CRe₁, the initial injection (Jaffe & Perola 1973) and continuous injection (Pacholczyk 1970) models are expected to account for electron losses and are thus poorly suited for CRp (see Turner et al. 2018, for discussions about the parameterization). We note that the minimum and maximum energy of the CRs is part of the parameters, and it is possible to use these parameters to truncate the spectra. By default, the minimum energy of CRp and CRe₁ corresponds to the energy threshold of the proton-proton interaction and the rest mass of the electrons, respectively. Finally, it is also possible to use this parametric function to inject CRe₁ and apply losses (see Sect. 4.2 for details) to obtain the actual electron population. In Fig. 3, example spectra are shown for the CRe₁. The effect of these spectra on the observables is shown in Sect. 5.

The radial and spectral distributions of the CRs are currently decoupled, and we can express the CR distribution (i.e., the CR number density per unit energy) as

$$\begin{array}{c}\hfill J_{\mathrm{CR}}(r,E)=A_{\mathrm{CR}}{f}_1(E)f_2(r),\end{array}$$(15)

where A_CR is the normalization, and f₁(E) and f₂(r) are the spectral and radial distributions, respectively (see Tables 2 and 3 for available models). In principle, functions f₁ and f₂ could be merged into f_1, 2(r, E) to include a radial dependence of the spectral component, but this function is not yet implemented. When losses are to be applied to the input distribution, a radial dependence affects the spectrum because the losses themselves depend on the radius (see Sect. 4.2 for details).

In order to normalize the CR distribution, we compute the energy density that is stored between energy E₁ and E₂, which can be expressed by integrating over the energy as

$$\begin{array}{c}\hfill u_{\mathrm{CR}}(r)=3P_{\mathrm{CR}}(r)={\displaystyle {\int }_{E_1}^{E_2}}{E}_{\mathrm{CR}}{J}_{\mathrm{CR}}(E_{\mathrm{CR}})\mathrm{d}{E}_{\mathrm{CR}},\end{array}$$(16)

and it is related to the CR pressure, P_CR. Here we assume that the CRs are ultrarelativistic particles, with adiabatic index Γ = 4/3. The result is only weakly sensitive to the upper bound, E₂ = E_max, CRp/e, because the CR spectrum generally vanishes rapidly for a spectral index higher than 2. The default lower bound is set to the minimum proton energy necessary to trigger the pion production, $E_1\equiv E_{\mathrm{p}}^{\mathrm{th}}$, for protons and to the electron rest mass for the electrons. The total energy stored in CRs enclosed within the radius R can then be computed as

$$\begin{array}{c}\hfill U_{\mathrm{CR}}(R)=4\pi {\displaystyle {\int }_0^R}{u}_{\mathrm{CR}}{r}^2\mathrm{d}r.\end{array}$$(17)

The CR-to-thermal energy density ratio is then given by

$$\begin{array}{c}\hfill x_{\mathrm{CR}}(r)=\frac{u_{\mathrm{CR}}(r)}{u_{\mathrm{th}}(r)},\end{array}$$(18)

or similarly,

$$\begin{array}{c}\hfill X_{\mathrm{CR}}(R)=\frac{U_{\mathrm{CR}}(R)}{U_{\mathrm{th}}(R)},\end{array}$$(19)

when integrated over the volume up to the radius R.

In practice, the CR distribution, A_CR, is normalized by setting the value of X_CR(R) at a given radius (e.g., R₅₀₀). We note that this fraction is defined relative to the enclosed energy here, while it is also common to find this definition in terms of pressure in the literature. The two differ by a factor of 2 because the thermal gas is nonrelativistic, while the CRs are in the relativistic regime.

The CR distributions can be integrated over energy as

$$\begin{array}{c}\hfill n_{\mathrm{CR}}(r)\equiv {\frac{\mathrm{d}N}{\mathrm{d}V}|}_{[E_1,E_2]}={\displaystyle {\int }_{E_1}^{E_2}}{J}_{\mathrm{CR}}(r,E)\mathrm{d}E\end{array}$$(20)

or radius as

$$\begin{array}{c}\hfill \frac{\mathrm{d}N}{\mathrm{d}E}(<R)={\displaystyle {\int }_0^R}4\pi r^2{J}_{\mathrm{CR}}(r,E)\mathrm{d}r\end{array}$$(21)

to compute the number density profile within E₁ and E₂, or the spectrum enclosed within R, respectively.

In Fig. 5 we illustrate the integrated CR number density profiles and spectra for the baseline cluster model. CRe₁ and CRp clearly follow the same profile because they are calibrated to follow the thermal electron number density. The number of electrons and protons is nearly the same given the chosen normalization. The number of CRs drastically decreases when a cut in energy is applied. The spectra show different shapes for the electrons and protons, reflecting our baseline choice (power law for the protons, and continuous injection with a break at 1 GeV for the electrons), and the minimum energy is also different for the two populations. The number of enclosed CRs naturally increases with increasing radius. We note that the figure would be very similar for the real cluster models because the underlying CR modeling is the same.

In Fig. 6 we show the magnetic field profiles of our cluster models in the left panel and the ratio between the CRp energy and the thermal energy in the right panel. The magnetic field was calibrated on the thermal density profile as $B\propto n_{\mathrm{e}}^{0.5}$, and normalized to 5 μG (at the peak for the baseline model, and at 100 kpc for the real clusters). For the CR-to-thermal energy, we also varied the CR number density profiles using different scaling relations with respect to the thermal density and pressure to show the changes in the resulting profiles. However, we note that the exact shape also depends on the shape of the thermal pressure, which is kept fixed here. Because the pressure profile decreases with radius, setting the CR distribution to a flatter profile leads to a deficit in the center and an increase in the outskirt for the energy ratio X_CR. When the CR number density profile follows the thermal density profile, the CR-to-thermal energy of cool-core clusters is boosted in the center because the thermal pressure is low relative to the thermal density in the core (see, e.g., the case of A1795). In all cases, the ratio was set to X_CR(R₅₀₀) = 10⁻².

4. Particle interactions in the ICM

The physical properties of the ICM for its thermal and nonthermal components have been defined in Sect. 3. In this section, we model the hadronic interactions that take place in the plasma and generate secondary particles (see also Fig. 1). We also discuss the loss processes that affect them, in particular, the electrons.

Fig. 7. Production rate of secondary particles from hadronic processes. Left: relative comparison of the production rate of γ-rays in the case of no helium and zero metallicity in the Kafexhiu et al. (2014) and the Kelner et al. (2006) parameterization, using the Pythia8 parameterization from Kafexhiu et al. (2014) as a reference. For Kelner et al. (2006), we also show the low- and high-energy (δ-approximation) limits as dashed lines, which were combined to compute the rate over the total energy range. Middle: effect of including helium and metals in the model. This was done with Pythia8 using the Kafexhiu et al. (2014) parameterization. Right: comparison of the production rate for γ-rays, electrons, positrons, and neutrinos, computed for a radius r = 100 kpc. We also include the power-law function proportional to E−2.4 for comparison because it corresponds to the injected CRp distribution.

4.1. Production rate of secondary particles from hadronic interactions

The collision between high-energy CRp and the thermal ambient gas produces γ-rays, electrons, positrons, and neutrinos, mainly through the production of pions, following these interaction chains:

$$\begin{array}{c}\hfill \{\begin{array}{cc}p+p⟶\hfill & {\pi }^0+{\pi }^{-}+{\pi }^{+}+\mathrm{others}\hfill \\ \hfill & {\pi }^0⟶2\gamma \hfill \\ \hfill & {\pi }^{\pm }⟶{\mu }^{\pm }+{\nu }_{\mu }/{\overline{\nu }}_{\mu }⟶e^{\pm }+{\nu }_e/{\overline{\nu }}_e+{\overline{\nu }}_{\mu }/{\nu }_{\mu .}\hfill \end{array}\end{array}$$(22)

When the thermal plasma is considered to be at rest with respect to the CRs, the CRp-thermal proton collision rate per unit energy of CRp is given by

$$\begin{array}{c}\hfill \frac{\mathrm{d}{N}_{\mathrm{col}}}{\mathrm{d}{E}_{\mathrm{CRp}}\mathrm{d}V\mathrm{d}t}={\sigma }_{\mathrm{pp}}\times v_{\mathrm{CRp}}\times n_{\mathrm{p}}\times J_{\mathrm{CRp}},\end{array}$$(23)

where σ_pp is the proton-proton interaction cross section, v_CRp ≃ c is the speed of CRp, n_p the number density of thermal protons, and J_CRp the number density per unit energy of CRp. The production rate of secondary particles X per unit volume, unit energy, and unit time can then be expressed as

$$\begin{array}{c}\hfill \frac{\mathrm{d}{N}_X}{\mathrm{d}{E}_X\mathrm{d}V\mathrm{d}t}={\displaystyle {\int }_{E_X}^{+\infty }}\frac{\mathrm{d}{N}_{\mathrm{col}}}{\mathrm{d}{E}_{\mathrm{CRp}}\mathrm{d}V\mathrm{d}t}{F}_X(E_X,E_{\mathrm{CRp}})\mathrm{d}{E}_{\mathrm{CRp}}.\end{array}$$(24)

Two ingredients are thus necessary: (1) the total inelastic cross section of the proton-proton interaction and its evolution as a function of energy, σ_pp(E_CRp); and (2) the number of secondary particles produced in a collision per unit energy of the produced particle as a function of the initial energy of the CRp, namely F_X(E_X, E_CRp). These ingredients are usually obtained by fitting parametric functions to accelerator data, together with Monte Carlo simulations performed with sophisticated codes (e.g., Kelner et al. 2006; Kamae et al. 2006; Kafexhiu et al. 2014).

The Kafexhiu et al. (2014) parameterization was implemented in Naima (Zabalza 2015), a publicly available Python package dedicated to the computation of nonthermal radiation from relativistic particle populations. The work presented here is based on Naima, to which the radial dimension was added, and thus, we also use the work by Kafexhiu et al. (2014) as our baseline.

As we show below, heavy elements also contribute significantly to the particle production rate. This contribution is only available in the work by Kafexhiu et al. (2014), which is also expected to have the most current data, in particular at the highest energies. However, Kafexhiu et al. (2014) only focused on the γ-ray production rate, while Kelner et al. (2006) also provided a parameterization for the leptons (electrons, positrons, and neutrinos), but did not include heavy elements. Therefore we employed a hybrid approach. We used the parameterization by Kafexhiu et al. (2014) for the γ-ray production and that from Kelner et al. (2006) for the leptons. To account for heavy elements in the case of leptons, we applied a rescaling of the production rate given by Kelner et al. (2006). To do so, we assumed that the ratio between the production rate of leptons and that of γ-rays does not depend on whether heavy elements are included. This is motivated by the fact that Kafexhiu et al. (2014) accounted for heavy elements using a multiplicative correction to the cross section. The production rate of leptons is finally given by

$$\begin{array}{c}\hfill \frac{\mathrm{d}{N}_{e^{\pm },{\nu }_{\mu ,e}}}{\mathrm{d}E\mathrm{d}V\mathrm{d}t}={\frac{\mathrm{d}{N}_{e^{\pm },{\nu }_{\mu ,e}}}{\mathrm{d}E\mathrm{d}V\mathrm{d}t}|}_{\mathrm{Kelner}2006}\times \frac{{\frac{\mathrm{d}{N}_{\gamma }}{\mathrm{d}E\mathrm{d}V\mathrm{d}t}|}_{\mathrm{Kafexhiu}2014}}{{\frac{\mathrm{d}{N}_{\gamma }}{\mathrm{d}E\mathrm{d}V\mathrm{d}t}|}_{\mathrm{Kelner}2006}}·\end{array}$$(25)

While this approach allowed us to compute the electron, positron, and neutrino production rate in the presence of helium and nonzero metallicity of the ICM, it uses the so-called δ-approximation for the ratio of leptons to γ-rays in the high-energy regime (see Kelner et al. 2006), which is expected to be a relatively crude approximation (see Kafexhiu et al. 2014). Nevertheless, the accuracy of the lepton-to- γ-ray ratio is expected to be a much better fit because biases in the spectra are expected to cancel each other out.

In Fig. 7 we illustrate the computation of the secondary particle production rate in the case of our baseline cluster model (see Sect. 3), with a power-law model with index 2.4 for the CRp. The left panel shows the relative difference between the Kafexhiu et al. (2014) and Kelner et al. (2006) parameterizations for the γ-ray production rate when the helium and metal abundances were set to zero. We used the Pythia8 parametrization from Kafexhiu et al. (2014) as a reference. In practice, the Kelner et al. (2006) parameterization is the combination of a calculation at low energy and the use of the δ-approximation at high energy, which are both shown as dashed lines. The agreement between Pythia8 and Kelner et al. (2006) is relatively good over most of the energy range (lower than 25% for most of it, but a peak reaches more than 100% around 100 MeV). The different high-energy parameterizations available in the work by Kafexhiu et al. (2014), namely using the Monte Carlo codes Pythia8, SIBYLL, QGSJET, or Geant4, are also shown. As expected, the difference with respect to Pythia8 is only large at high energy. It remains below 25% for energies below 1 TeV, and increases to more than 50% above 100 TeV (see Kafexhiu et al. 2014, for further discussions). Based on the comparison of the top panels of Fig. 7, systematic uncertainties in the modeling are expected to be about 30% over most of the energy range probed here.

The middle panel quantifies the effect of accounting for helium and metals in the ICM. The helium mass fraction was chosen to be 0.27, and the metallicity was set to the solar value. The helium can clearly boost the signal by more than 50%, especially at low energies, but has a strong effect over the full energy range where its contribution remains higher than 40%. The metals, on the other hand, only account for percent-level changes in the spectrum. We note that the ratio μ_e/μ_p depends on the ICM composition, which affects the value of n_p for fixed n_e, and explains why the γ-ray production rate can become lower when metals are included compared to the helium-only case; this is visible around 100 MeV. Based on these results, we expect that ignoring the helium contribution will lead to a strong systematic effect in the model amplitude that underestimates the signal by about 40−50%.

Finally, the right panel of Fig. 7 provides the particle injection rate for γ-rays, electrons, and positrons, and both muonic and electronic neutrinos. The high-energy cutoff is due to the maximum energy of the CRp, which was set to 10 PeV, while the decrease below 1 GeV is due to the kinematic production threshold of the proton-proton interaction of about 1.2 GeV. Between these energies, the slope of the secondary particle production rate nearly follows that of the injected CRp, as shown by the dotted black line.

4.2. Energy losses

While the neutrinos and the γ-rays can escape the cluster and be detected by ground- and space-based instruments, the electrons evolve in the ICM and are affected by several sources of energy loss. We considered the main sources of energy loss: synchrotron radiation, inverse Compton interaction, Coulomb losses, and Bremsstrahlung radiation.

Fig. 8. Contributions to the energy loss rate as a function of energy for a radius r = 10 kpc and r = 1000 kpc.

We first define the Lorentz factor of the electrons $\gamma =\frac{E}{m_{\mathrm{e}}{c}^2}$, and the reduced speed of the electrons, $\beta =\sqrt{1-1/{\gamma }^2}$. The synchrotron radiation loss is given (in S.I. units) by Longair (2011),

$$\begin{array}{c}\hfill {\frac{\mathrm{d}E}{\mathrm{d}t}|}_{\mathrm{sync}}=-\frac{4}{3}{\sigma }_{\mathrm{T}}c{\beta }^2{\gamma }^2\frac{B^2}{2{\mu }_0},\end{array}$$(26)

with σ_T the Thomson cross section and μ₀ the vacuum permeability. It is proportional to the amplitude of the magnetic field squared, B². We therefore expect it to be most efficient in the central regions of the cluster. Because of the γ² dependency, the synchrotron loss will be higher at high energy.

Inverse Compton losses can be expressed in a very similar way (Longair 2011) as

$$\begin{array}{c}\hfill {\frac{\mathrm{d}E}{\mathrm{d}t}|}_{\mathrm{IC}}=-\frac{4}{3}{\sigma }_{\mathrm{T}}c{\beta }^2{\gamma }^2{u}_{\mathrm{CMB}},\end{array}$$(27)

where the dependence on magnetic field is replaced by the ambient photon field, assumed to be dominated by the CMB, whose energy density is given by

$$\begin{array}{c}\hfill u_{\mathrm{CMB}}=\frac{8{\pi }^5{\left(k_{\mathrm{B}}{T}_{\mathrm{CMB}}(1+z)\right)}^4}{15{(hc)}^3}·\end{array}$$(28)

The inverse Compton energy losses do not depend on the cluster location, but increase with redshift because of the CMB dependence.

The Coulomb losses are computed as (Gould 1972)

$$\begin{array}{cc}\hfill {\frac{\mathrm{d}E}{\mathrm{d}t}|}_{\mathrm{Coulomb}}=& -\frac{3}{2}{\sigma }_{\mathrm{T}}{n}_{\mathrm{e}}\frac{m_{\mathrm{e}}{c}^3}{\beta }(\ln \left(\frac{m_{\mathrm{e}}{c}^2\beta \sqrt{\gamma -1}}{ħ{\omega }_{\mathrm{p}}}\right)\hfill \\ \hfill & -\ln \left(2\right)(\frac{{\beta }^2}{2}+\frac{1}{\gamma })+{\left(\frac{\gamma -1}{4\gamma }\right)}^2+\frac{1}{2}),\hfill \end{array}$$(29)

where ${\omega }_{\mathrm{p}}=\sqrt{\frac{e^2{n}_{\mathrm{e}}}{m_{\mathrm{e}}{ϵ}_0}}$ is the plasma frequency. The Coulomb losses are proportional to the thermal electron number density and are therefore more effective in the cluster core. They are almost independent of energy.

Finally, the Bremsstrahlung losses are computed as (Blumenthal & Gould 1970)

$$\begin{array}{c}\hfill {\frac{\mathrm{d}E}{\mathrm{d}t}|}_{\mathrm{Brem}.}=-8\alpha c{r}_0^{2}E(n_{\mathrm{p}}+3n_{\mathrm{He}})(\ln \left(2\gamma \right)+\frac{1}{3}),\end{array}$$(30)

where $\alpha =\frac{e^2}{4\pi {ϵ}_0ħc}$ is the fine-structure constant and $r_0=\frac{e^2}{4\pi {ϵ}_0{m}_{\mathrm{e}}{c}^2}$ is the classical electron radius. We neglected elements heavier than helium and used the completely unscreened limit that is appropriate for low-density plasma. As for the Coulomb losses, the Bremsstrahlung losses depend on the number density of thermal nuclei, but they increase with energy.

In Fig. 8 we provide the loss function for the synchrotron, inverse Compton, Bremsstrahlung, and Coulomb contributions for two different radii from the center, 10 and 1000 kpc, for our baseline model. At low energy, the Coulomb losses are expected to dominate, while at high energy, the synchrotron and inverse Compton losses dominate, depending on the relative value of the magnetic field and the CMB photon field. The contribution by Bremsstrahlung is always subdominant.

Fig. 9. Left: profile of the CRe2, taken at 10 MeV and 10 TeV (after applying a multiplicative factor of 1015). Right: spectrum of the secondary electron at 100 and 1000 kpc.

4.3. Secondary electrons in the steady-state approximation

When CRs are injected into the ICM, their evolution is expected to follow the diffusion-loss equation (Berezinskii et al. 1990),

$$\begin{array}{cc}\hfill \frac{\partial n(E,r,t)}{\partial t}=& -\mathbf{\nabla }(n(E,r,t)\mathit{v})\hfill \\ \hfill & +\mathrm{\nabla }[D(E)\mathrm{\nabla }n(E,r,t)]\hfill \\ \hfill & +\frac{\partial }{\partial E}[\ell (E,r)n(E,r)]+q(E,r,t),\hfill \end{array}$$(31)

where $n(E,r,t)\equiv \frac{\mathrm{d}{N}_{\mathrm{CRe}}}{\mathrm{d}E\mathrm{d}V}$ is the number of CRe per unit volume and energy, q(E, r, t) is the injection rate, v is the ICM velocity, and D(E) is the diffusion coefficient, and where the loss function, ℓ(E, r), is given by

$$\begin{array}{c}\hfill \ell (E,r)=-({\frac{\mathrm{d}E}{\mathrm{d}t}|}_{\mathrm{sync}}+{\frac{\mathrm{d}E}{\mathrm{d}t}|}_{\mathrm{IC}}+{\frac{\mathrm{d}E}{\mathrm{d}t}|}_{\mathrm{Coulomb}}+{\frac{\mathrm{d}E}{\mathrm{d}t}|}_{\mathrm{Brem}.}).\end{array}$$(32)

Assuming that the CRe do not significantly diffuse, and assuming steady-state condition, we write the number density of CRe at equilibrium as

$$\begin{array}{c}\hfill \frac{\mathrm{d}{N}_{\mathrm{CRe}}}{\mathrm{d}E\mathrm{d}V}(E,r)=\frac{1}{\ell (E,r)}{\displaystyle {\int }_{\mathrm{E}}^{\infty }}q(ϵ,r)\mathrm{d}ϵ,\end{array}$$(33)

where q(ϵ, r) is computed as the output of Eq. (25).

We note that we did not account explicitly for the possible reacceleration of seed electrons by ICM turbulences (e.g., Brunetti et al. 2017, and references therein). However, because we did not model the details of the microphysics here, such a population might be included in the CRe₁, which are modeled independently, as discussed in Sect. 3.3. By doing so, we would implicitly assume that the reacceleration process, which would effectively contribute to a loss for our secondary electron population (i.e., a population transfer) is subdominant with respect to the losses. In this case the physical consistency between the different CR populations will not necessarily be verified. While reacceleration models are beyond the scope of the current work, we leave room for implementing reacceleration options in the MINOT code in the future.

In Fig. 9 we present the radial profile and the spectra of the CRe₂ in the steady-state approximation, with no diffusion. The profile becomes steeper at higher energy because inverse Compton and synchrotron losses become more important in this regime relative to the Coulomb loss, which is more efficient in the core. Moreover, electrons accumulate around 100 MeV because lower energy electrons quickly disappear because of Coulomb losses and higher energy electrons are more affected by inverse Compton and synchrotron losses. We also provide the same profiles and spectra when the magnetic field and thermal density are higher by a factor of 10. In the latter case, we also decreased the number of CRp by a factor of 10 so that the rate of proton-proton collisions was conserved. An increase in thermal plasma density leads to a much flatter profile because Coulomb losses are far higher in the core where the density is high. An increased magnetic field also flattens the profile, but less drastically because the magnetic field profile itself is flatter. In the spectrum, the magnetic field has a stronger effect at high energy, while the increase in thermal plasma density leads to more losses at low energy. The peak of the secondary electron spectrum therefore depends on the competition between the energy losses in the magnetic field and in thermal plasma.

5. Multiwavelength observables

In this section, the physical properties of the cluster (Sect. 3) and the production of secondary particles in the ICM (Sect. 4) are used to compute the observables of galaxy clusters related to the diffuse gas component. This includes the tSZ effect, the thermal X-ray emission, the radio synchrotron emission, the inverse Compton emission, and the γ and neutrino emission from hadronic processes. We focus here on showing the effect of model changes on the observables using our baseline cluster model.

5.1. General considerations

In general, the cluster observables are associated with physical processes at play in the ICM, which can be described in terms of production rate (this does not strictly apply to the tSZ signal because it is a spectral distortion, as we discuss in Sect. 5.3). We define $Q(r,E)\equiv \frac{\mathrm{d}N}{\mathrm{d}E\mathrm{d}V\mathrm{d}t}$, the emission rate associated with the physical process considered. For instance, in the case of X-ray emission, Q would be the number of X-ray photons emitted per unit volume, per unit of time, and per unit energy in the ICM.

The surface brightness (or flux per solid angle) at a projected distance R from the center is therefore given by integrating Q(r, E) over the line of sight as

$$\begin{array}{cc}\hfill \frac{\mathrm{d}N}{\mathrm{d}E\mathrm{d}S\mathrm{d}t\mathrm{d}\mathrm{\Omega }}(R,E)& =\frac{D_{\mathrm{A}}^2}{4\pi D_{\mathrm{L}}^2}{\displaystyle {\int }_{-\infty }^{+\infty }}Q(r)\mathrm{d}\ell \hfill \\ \hfill & =\frac{D_{\mathrm{A}}^2}{4\pi D_{\mathrm{L}}^2}{\displaystyle {\int }_R^{R_{\max }}}\frac{2rQ(r)}{\sqrt{r^2-R^2}}\mathrm{d}r,\hfill \end{array}$$(34)

where the factor $D_{\mathrm{A}}^2$ accounts for the conversion from physical area into solid angle, and the normalization by $4\pi D_{\mathrm{L}}^2$ assumes that the emission is isotropic. We note that Eq. (34) is valid in the small-angle approximation (i.e., assumes the cluster size to be small against the distance to the observer), which is expected to be accurate for our purpose because the extent of clusters never exceeds a few degrees. It also neglects the redshift extent of the cluster.

Based on Eq. (34), we now wish to compute several quantities accessible from observations: (1) the surface brightness profile; (2) the spectrum, by integrating the signal over the cluster volume or the solid angle; (3) the total flux, by integrating over both the volume and the energy (or at fixed energy); and (4) the map of the signal, which in our case is equivalent to the surface brightness profile because of azimuthal symmetry, but allows us to generate spatial templates for dedicated analysis.

The surface brightness profiles (and maps), are computed by log-log integration of the quantity $\frac{\mathrm{d}N}{\mathrm{d}E\mathrm{d}S\mathrm{d}t\mathrm{d}\mathrm{\Omega }}(R,E)$ over the requested energy range, or simply by fixing the energy to the one required by the observation (as is done, e.g., in the radio and millimeter domain). In the case of the map, the signal is projected on a grid corresponding to the header (or sampling properties) set by the user. An option allows the user to normalize the map to the total flux so that the map only accounts for the spatial dependence of the signal in units proportional to the inverse solid angle. This proves useful for a γ-ray analysis, for example, in which image templates are needed.

The spectra are computed in a similar way, by log-log integration over the volume. In this case, two possibilities are available. (1) Integration over the solid angle within a circle of radius R_max, so that the total integration volume is a cylinder. (2) The emission rate Q(r, E) can be integrated spherically up to R_max before normalization by $4\pi D_{\mathrm{L}}^2$. The two quantities only differ by the definition of the integration volume, and should converge when all the cluster emission is accounted for with increasing R_max. The cylindrical integration resembles more what would be accessible directly from observations, while the spherical integration is a more natural from a physical point of view because it returns a quantity that is computed in a single physical (3D) radius.

In order to compute the flux, we integrated over the cluster volume and the energy. Alternatively, the energy might be fixed, as discussed above. By definition, the luminosity of a given source within the energy band ΔE ≡ [E₁, E₂] is given by

$$\begin{array}{c}\hfill L_{\mathrm{\Delta }\mathrm{E}}=4\pi D_{\mathrm{L}}^2{F}_{\mathrm{\Delta }\mathrm{E}}.\end{array}$$(35)

We also applied the redshift stretching to the energy of the photons. This function can be switched off by the user.

Fig. 10. Observables associated with the X-ray emission. Left: X-ray spectrum within R500 shown for both MEKAL and APEC models, as well as with and without galactic photoelectric absorption. The MEKAL model is difficult to distinguish because it coincides very well with the APEC model. Right: surface brightness profile in the band 0.1−2.4 keV. The dynamical range of the profile amplitude has been set to the same value for all observables.

5.2. Thermal X-ray emission

Because of the gas temperature (a few keV) and the density (typically 10⁻³ − 10⁻⁵ cm⁻³), the leading emission process at X-ray energies is thermal Bremsstrahlung (see Sarazin 1986; Böhringer & Werner 2010, for reviews). The X-ray emission is thus a direct probe for the thermal gas density. It presents a characteristic exponential cutoff at high energies, determined by the gas temperature. Heavy elements also induce a large number of spectral lines. The X-ray surface brightness is generally expressed as

$$\begin{array}{c}\hfill S_{\mathrm{X}}=\frac{1}{4\pi {(1+z)}^4}\int n_{\mathrm{e}}^2\mathrm{\Lambda }(T_{\mathrm{e}},Z)\mathrm{d}l,\end{array}$$(36)

where Λ(T_e, Z) is the cooling function, which varies with temperature.

In practice, MINOT uses the XSPEC software (Arnaud 1996) to directly compute the counts using either the MEKAL or APEC X-ray plasma spectral models. These models require the ICM abundance, the redshift, the temperature, the energy range, and a normalization defined as

$$\begin{array}{c}\hfill \mathrm{norm}=\frac{{10}^{-14}}{4\pi {\left(\left(\frac{D_{\mathrm{A}}}{1\mathrm{cm}}\right)(1+z)\right)}^2}\int \left(\frac{n_{\mathrm{e}}}{1{\mathrm{cm}}^3}\right)\left(\frac{n_{\mathrm{H}}}{1{\mathrm{cm}}^3}\right)\mathrm{d}V.\end{array}$$(37)

MINOT also accounts for the foreground photoelectric absorption using the value of the hydrogen column density at the cluster location. The XSPEC outputs are then normalized to compute the emission rate (counts or energy) per unit volume and time. It is also possible to account for the response function of X-ray satellites, so that the outputs are normalized by the effective area of the observation. The spectrum, surface brightness profile and maps, and flux are then extracted as discussed in Sect. 5.1.

In Fig. 10 we illustrate these X-ray observables in the case of our baseline cluster. The two models MEKAL and APEC agree well. The dashed blue spectrum shows the effect of the photoelectric absorption from the foreground, leading to low energy cuts. At high energy, the exponential cutoff is clearly visible, and the spectral lines are also visible below 10 keV. The raw signal associated with the inverse Compton emission is also shown; it is discussed in Sect. 5.6. We note that it is well below the thermal X-ray emission, but might become significant with increasing energy. The surface brightness profile, computed between 0.1 keV and 2.4 keV, drops very rapidly in the outskirt because the signal is proportional to the density squared. The integrated flux reaches about 1 ph cm⁻² s⁻¹ at large radii.

Fig. 11. Observables associated with the tSZ effect. Left: tSZ spectrum within R500. We also provide the spectrum in the case when relativistic corrections are neglected for illustration. Right: compton parameter profile. The dynamical range of the profile amplitude has been set to the same value for all observables.

5.3. Thermal Sunyaev-Zel’dovich signal

The tSZ effect distorts the CMB blackbody spectrum as a result of inverse Compton scattering onto energetic thermal electrons (see Birkinshaw 1999; Mroczkowski et al. 2019, for reviews). Because it is a spectral distortion, it does not suffer from redshift dimming, and the general considerations of Sect. 5.1 do not strictly apply here. The change in surface brightness is expressed as

$$\begin{array}{c}\hfill \frac{\mathrm{\Delta }{I}_{\mathrm{tSZ}}}{I_0}=yf(x,T_{\mathrm{e}})\end{array}$$(38)

with respect to the CMB, $I_0=\frac{2{(k_{\mathrm{B}}{T}_{\mathrm{CMB}})}^3}{{(hc)}^2}\simeq 270.1$ MJy sr⁻¹. The parameter y is the so-called Compton parameter, which gives the normalization of the tSZ effect. It provides a measurement of the thermal electron pressure integrated along the line of sight as

$$\begin{array}{c}\hfill y=\frac{{\sigma }_{\mathrm{T}}}{m_{\mathrm{e}}{c}^2}\int P_{\mathrm{e}}\mathrm{d}\ell .\end{array}$$(39)

The frequency dependence of the tSZ effect is given by

$$\begin{array}{c}\hfill f(x,T_{\mathrm{e}})=\frac{x^4{e}^x}{{(e^x-1)}^2}(x\coth \left(\frac{x}{2}\right)-4)(1+{\delta }_{\mathrm{tSZ}}(x,T_{\mathrm{e}})),\end{array}$$(40)

where $x=\frac{h\nu }{k_{\mathrm{B}}{T}_{\mathrm{CMB}}}$. The term δ_tSZ(x, T_e) is a relativistic correction that introduces a small temperature dependence to the tSZ effect, and becomes important when the temperature becomes higher than about 10 keV, depending on the frequency. The relativistic correction was implemented following Itoh & Nozawa (2003), who are expected to be accurate at the percent level up to 50 keV. When relativistic corrections are neglected, the tSZ spectrum is null at 217 GHz, negative below (with a minimum around 150 GHz), and positive above (peaking at about 350 GHz).

The tSZ integrated flux, often used to track the cluster mass, can be expressed as

$$\begin{array}{c}\hfill Y_{\mathrm{cyl}}(R)={\displaystyle {\int }_0^R}2\pi ry\mathrm{d}r\end{array}$$(41)

in the case of cylindrical integration, or as

$$\begin{array}{c}\hfill Y_{\mathrm{sph}}(R)=\frac{{\sigma }_{\mathrm{T}}}{m_{\mathrm{e}}{c}^2}{\displaystyle {\int }_0^R}4\pi r^2{P}_{\mathrm{e}}\mathrm{d}r\end{array}$$(42)

for the spherically integrated flux. The integrated flux, Y_cyl, sph, can be expressed in units of surface or normalized by $D_{\mathrm{L}}^2$ to be homogeneous to solid angle, as is commonly done in the literature.

In Fig. 11 we illustrate the spectrum and profile (shown in terms of the Compton parameter) of our reference cluster. The tSZ flux can be significant down to a few GHz, and could thus affect radio synchrotron observations. On the other hand, the synchrotron emission is not shown here, but is far lower than the tSZ signal in the considered frequency range. Given the linear sensitivity of the tSZ signal to the pressure, the profile is relatively flat (e.g., compared to the X-ray surface brightness).

Fig. 12. Effect of the EBL on the normalized spectra of a baseline model cluster at z = 0.02. All the available models are shown as indicated in the legend. Because of the redshift stretching, the shape of the γ-ray spectrum slightly changes even in the case without EBL absorption (“none”).

5.4. γ-ray hadronic emission

In the case of high-energy photons, the absorption by the extragalactic background light (EBL, Dwek & Krennrich 2013) needs to be accounted for, which is thought to have been produced by the sum of all light contributions (e.g., starlight or dust reemission) at all epochs in the Universe. While traveling from the cluster to the Earth, γ-rays may interact with the EBL through electron-positron pair production, and thus be effectively absorbed along the way as

$$\begin{array}{c}\hfill \frac{\mathrm{d}N}{\mathrm{d}E\mathrm{d}S\mathrm{d}t}⟶\frac{\mathrm{d}N}{\mathrm{d}E\mathrm{d}S\mathrm{d}t}\times \exp (-\tau (E)),\end{array}$$(43)

with τ(E) the optical depth. EBL absorption depends on redshift and on the energy of the γ-rays. To account for EBL absorption, we used the ebltable Python package⁷, which reads in and interpolates tables for the photon density of the EBL and the resulting opacity for high-energy γ-rays. This package provides different models for the EBL based on Franceschini et al. (2008), Kneiske & Dole (2010), Finke et al. (2010), Domínguez et al. (2011), and Gilmore et al. (2012). We illustrate the effect of the EBL for various models available in Fig. 12 for a cluster redshift z = 0.02. While it is crucial to account for the EBL, especially at very high energies, the uncertainties associated with the EBL models are expected to be small.

The intergalactic magnetic fields probably also affect the predictions for the γ-ray observables (e.g., Neronov & Semikoz 2009). However, because of the current uncertainties on the properties of the magnetic field, its effect remains uncertain. While this effect is not yet implemented, it is being considered for the future improvement of the MINOT code.

The γ-ray production rate resulting from hadronic interaction is computed following Kafexhiu et al. (2014), as described in Sect. 4.1. We compute the spectrum (within R₅₀₀, using spherical integration) and profile as detailed in Sect. 5.1. These quantities are displayed in Fig. 13 for our baseline cluster model. The spectrum peaks at GeV energies, quickly vanishes at lower energies, and is affected by a cutoff at high energy due to the EBL. We also show the effect of the change in CRp slope on the γ-ray spectrum for a fixed-normalization X_CRp. In the case of a flatter CRp profile, the amplitude of the spectrum is reduced because the number of proton-proton collisions is reduced by the lower spatial coincidence of thermal and CRp. The inverse Compton signal computed for our baseline model is also shown for comparison. It is below the hadronic emission except at low energy. The profile presents a compact signal because it arises from the product of the thermal electron number density and the CR density. When the CR number density profile is flattened, the γ-ray signal itself becomes flatter. Nevertheless, it still decreases with radius even for a completely flat CR number density profile because the thermal density profile remains peaked. The inverse Compton signal is expected to be flatter than the signal of the hadronic emission, but it is also lower in amplitude in our baseline model. The integrated spectrum between 1 GeV and 1 TeV almost reaches 10⁻⁹ photon cm⁻² s⁻¹ in the case of this baseline model.

Fig. 13. Observables associated with the γ-ray hadronic emission. The signal coming from inverse Compton interactions is also shown for comparison. Left: γ-ray spectrum within R500. Right: γ-ray profile integrated between 1 GeV and 1 TeV. The dynamical range of the profile amplitude has been set to the same value for all observables.

5.5. Neutrino (hadronic) emission

In contrast to the γ-rays, the neutrinos are not affected by the EBL absorption. Except for the EBL, their observables are computed in the same way as the γ-rays associated with hadronic interactions (Sect. 5.4); the production rate is computed following a combination of Kelner et al. (2006) and Kafexhiu et al. (2014), as described in Sect. 4.1.

Figure 14 illustrates the neutrino observable both for the muonic and electronic neutrinos in the case of the spectrum and the profile, and for the sum of the two for the flux. Because of the neutrino oscillation, we would in practice expect a mix between the ratio of neutrinos of different flavors. Because the processes associated with the neutrino is the same as that of the γ-rays, their observables are very similar to one another, except for a small difference in the normalization.

Fig. 14. Observables associated with neutrino hadronic emission. Left: neutrino spectrum within R500. Right: neutrino profile integrated between 1 GeV and 1 TeV. The dynamical range of the profile amplitude has been set to the same value for all observables.

Fig. 15. Observables associated with inverse Compton emission. Left: inverse Compton spectrum within R500, with the contribution from both primary and secondary electrons. Right: inverse Compton profile integrated between 1 GeV and 1 TeV. The dynamical range of the profile amplitude has been set to the same value for all observables.

5.6. Inverse Compton emission

The inverse Compton emission is also affected by the EBL absorption in the high-energy limit. We refer to Sect. 5.4 and Fig. 13 for this effect.

We use the analytical approximation for the treatment of inverse Compton scattering of relativistic electrons in the CMB blackbody radiation field given by Khangulyan et al. (2014), which is expected to be accurate within 1% uncertainty throughout the application range. In particular, their Eq. (14) gives us the number of inverse Compton photons produced per unit energy and time per CRe as a function of the CRe energy, $\frac{\mathrm{d}{N}_{\mathrm{IC}}}{\mathrm{d}E\mathrm{d}t}$, which we express as

$$\begin{array}{c}\hfill \frac{\mathrm{d}{N}_{\mathrm{IC}}}{\mathrm{d}{E}_{\mathrm{IC}}\mathrm{d}t}=\frac{2r_0^2{m}_{\mathrm{e}}^3{c}^4}{\pi {ħ}^3}{\left(\frac{k_{\mathrm{B}}{T}_{\mathrm{CMB}}(1+z)}{E_{\mathrm{CRe}}}\right)}^2{G}_{\mathrm{IC}}(E_{\mathrm{IC}},E_{\mathrm{CRe}}),\end{array}$$(44)

with G_IC(E_IC,E_CRe) an analytical function computed following Khangulyan et al. (2014), using the approximation given by their Eq. (24). We thus integrate this quantity over the electron energy, accounting for the amount of CRe in the ICM. The emissivity is expressed as

$$\begin{array}{c}\hfill \frac{\mathrm{d}{N}_{\mathrm{IC}}}{\mathrm{d}{E}_{\mathrm{IC}}\mathrm{d}V\mathrm{d}t}=\int J_{\mathrm{e}}(E_{\mathrm{CRe}})\frac{\mathrm{d}{N}_{\mathrm{IC}}}{\mathrm{d}{E}_{\mathrm{IC}}\mathrm{d}t}\mathrm{d}{E}_{\mathrm{CRe}},\end{array}$$(45)

where $J_{\mathrm{e}}(E_{\mathrm{CRe}})\equiv \frac{\mathrm{d}{N}_{\mathrm{CRe}}}{\mathrm{d}{E}_{\mathrm{CRe}}\mathrm{d}V}$ is the CRe number density, summing the contributions from primary and secondary electrons.

In Fig. 15 we illustrate the observables associated with inverse Compton emission for our baseline cluster model. The spectrum shape reflects the complex processing of secondary electrons by their production rate in hadronic interaction, their losses in the ICM, and the production of inverse Compton after having also summed the CRe₁. In particular, changing the distribution of CRe₁ to an initial injection model removes relativistic electrons with energy higher than E_break, and thus does not remove the inverse Compton emission at energies above ∼E_break/100. In contrast, when a power-law model is used, more CRe are present at low energy, which increases the inverse Compton emission in this regime. The secondary electrons change the slope of the inverse Compton emission according to the CRp slope. In addition, as in the case of hadronic emission, flattening the CRp profile decreases the production of secondary particles, and thus that of inverse Compton associated with secondary electrons. The inverse Compton profile is relatively flat for primary electrons because the signal depends linearly on the CRe spatial distribution. For the secondary electrons, the profile is more compact because the secondary electrons are produced proportionally to the product of the thermal gas and the CRp number density. For both populations, we also show the effect of flattening the profiles of either the CRe₁ or the CRp. The integrated flux reaches more than 10⁻¹¹ ph cm⁻² s⁻¹ between 1 GeV and 100 TeV in our baseline model.

Fig. 16. Observables associated with the synchrotron emission. Left: synchrotron spectrum within R500, with the contribution from both primary and secondary electrons. Right: synchrotron emission profile at 100 MHz. The dynamical range of the profile amplitude has been set to the same value for all observables.

5.7. Radio synchrotron emission

The MINOT code focuses on the diffuse galaxy cluster emission associated with the bulk of the X-ray emitting ICM. In the case of the diffuse synchrotron emission, we therefore focus on the emission associated with radio halos and leave radio shocks (or relics, see van Weeren et al. 2019) aside. Because the orientation of the magnetic field is expected to be chaotic in the bulk ICM regions of a galaxy cluster, we need to average the standard energy distribution of the synchrotron emission over the directions of the field orientation. To do so, we follow the results of Aharonian et al. (2010), Appendix D, in which the orientation of the magnetic field is assumed to be randomized. This provides a convenient approximation with an accuracy better than 0.2% over the entire energy range (see Aharonian et al. 2010, for more details on the approximation and its accuracy).

As in the case of inverse Compton emission, we express

$$\begin{array}{c}\hfill \frac{\mathrm{d}{N}_{\mathrm{sync}}}{\mathrm{d}{E}_{\mathrm{sync}}\mathrm{d}V\mathrm{d}t}=\int J_{\mathrm{e}}(E_{\mathrm{CRe}})\frac{\mathrm{d}{N}_{\mathrm{sync}}}{\mathrm{d}{E}_{\mathrm{sync}}\mathrm{d}t}\mathrm{d}{E}_{\mathrm{CRe}},\end{array}$$(46)

where

$$\begin{array}{c}\hfill \frac{\mathrm{d}{N}_{\mathrm{sync}}}{\mathrm{d}{E}_{\mathrm{sync}}\mathrm{d}t}=\frac{3e^{3}B}{8{\pi }^2{ϵ}_0{m}_{\mathrm{e}}cħE_{\mathrm{sync}}}\stackrel{\sim }{G}(E_{\mathrm{sync}}/E_{\mathrm{c}}).\end{array}$$(47)

The quantity $E_{\mathrm{c}}=\frac{3eBħ{\gamma }^2}{2m_{\mathrm{e}}}$ is the synchrotron characteristic energy and G̃(x) the emissivity function of synchrotron radiation, which quickly increases from x = 0 to x ≃ 0.23 and smoothly vanishes for increasing x (see Aharonian et al. 2010).

In Fig. 16 we illustrate the observables associated with synchrotron emission, including the contribution from primary and secondary electrons. As for the inverse Compton case, the emission reflects the complex processing of secondary electrons, while it is more direct for the primary electrons. In particular, the curvature in the synchrotron emission is due to the losses of CRe at high energies. When the CRe₁ population model is changed to an initial injection scenario, the high-frequency curvature is significantly enhanced by the lack of very high energy electrons. In the case of a power-law CRe population, which extends to a lower energy, the spectrum is enhanced at low frequency. The slope of the CRp for secondary electrons is directly reflected in the synchrotron spectrum slope. Similarly as in the case of hadronic γ-ray emission, flattening the CRp profile decreases the amount of synchrotron emission because fewer secondary electrons are produced. As highlighted in the figure, the tSZ contribution might be high at high frequencies, and it might mimic a curved synchrotron spectrum if not accounted for. The synchrotron profile is more compact than that of the inverse Compton emission because it also depends on the magnetic field, which decreases with radius. For the inverse Compton emission, the profile is more compact for secondary electrons because it depends on the product of the thermal density and the CRp number density. We also show that the flattening of the CR population is reflected in the synchrotron profile. The flux at 100 MHz reaches typical values of 10 Jy in our baseline model.

6. Comparison to the literature

In order to validate the modeling and to further illustrate the use of the MINOT code, we compared the model predictions to results obtained in the literature and existing data. To do so, we used the models of the XCOP clusters defined in Sect. 2.7 and Table 4. We focused on millimeter and X-ray data for the thermal part and on γ-rays and radio data for the nonthermal component. Because the sensitivity of current neutrino telescopes is limited and because of the typical predicted fluxes, we leave the neutrino emission model predictions aside.

6.1. Thermal gas

First we compared the measured X-ray luminosity given by Eckert et al. (2017) to the luminosity we recover, integrated within R₅₀₀. We used the same values of R₅₀₀ and set the MINOT cosmological model to the one used in Eckert et al. (2017) to mitigate differences. The comparison of the obtained rest frame luminosity are given in Table 5. We obtain comparable luminosities, with differences of a few percent for A1795 and A2142, but differences of up to 19% for A2255. Because our model is based on the interpolation of the results by the XCOP project, we expect consistency between the two. However, our model has been defined by extrapolating the thermal plasma density profiles down to small radii. While these profile were measured down to about 10 kpc or even less for A1795 and A2142, it was only measured down to about 30 kpc for A2255. Thus, uncertainties coming from the extrapolation are likely to be larger for this cluster. The differences that we observe in Table 5 are thus likely due to extrapolation uncertainties.

Fig. 17. Comparison between the thermal gas observables and the MINOT model prediction. Left: comparison between ROSAT X-ray images and the MINOT model prediction. The scale is linear from −10 to 10 s−1 deg−2 and logarithmic above. Contours are as ±5 × 2N s−1 deg−2, with N the contour index starting at 0. Right: comparison between the Planck Compton parameter MILCA image and the MINOT model prediction. The contours are multiples of 3σ, where σ is the rms of the residual map. The Planck 10 arcmin FWHM beam is shown as a filled white circle. All maps where smoothed (in addition the intrinsic resolution) with a 2 arcmin FWHM gaussian beam for visual purpose. All the sky patches shown here are 1° × 1°.

We also compared the tSZ flux computed using MINOT to the flux given in the Planck PSZ2 catalog (Planck Collaboration XXVII 2016). Again, our model is based on the XCOP outputs, which are themselves obtained using Planck data, so that we expect consistency. Table 6 shows that the total integrated Compton parameter (i.e., computed within 5 R₅₀₀) agrees well for A1795 and A2142, and it differs by 2.3σ for A2255. Because A2255 is a strongly merging cluster, the flux arising from the cluster outskirt is likely to be higher than for A1795 and A2142, and the truncation radius involved in MINOT (set at 5 Mpc) may not be enough to account for all the tSZ flux. When we increase the truncation radius to 10 R₅₀₀, the flux rises to nearly 18 × 10⁻³ arcmin², which better agrees with the PSZ2 value. As in the case of the X-ray luminosity, the differences that we observe are thus likely due to the interpolation and assumption that we have made in defining the clusters.

In addition to the fluxes, we also directly compared the X-ray and tSZ images to existing data. We used the publicly available ROSAT (Truemper 1993) X-ray pointed data⁸ obtained for our three targets to produce maps in the [0.1, 2.4] keV energy band. The maps were subtracted from the background and normalized by the exposure. The ROSAT PSPC response matrices are accounted for in MINOT through the use of XSPEC, as well as the hydrogen column density taken at the location of each cluster, as obtained by the LAB survey (Kalberla et al. 2005). The model is projected on the same header as the original ROSAT data and accounts for the ROSAT effective area in units of counts per unit of time and solid angle. We accounted for the point spread function (PSF) by smoothing our model with a mean effective Gaussian function with a full width at half maximum (FWHM) of 30 arcsec, which is the typical number expected for ROSAT pointed observations. While the detailed analysis of the X-ray data including all instrumental effects is beyond the scope of this work, this comparison already provides a useful qualitative comparison to our modeling. The data, model, and residual images are displayed in the left panel of Fig. 17. The data and the model are shown on a log scale, and the residual is shown on a linear scale. While the overall agreement is good for all three clusters, several features are evident. First, many point sources that are not accounted for here affect the residual. Then, the real clusters are not perfectly azimuthally symmetric, which is shown for the residual with a positive and negative butterfly shape in the core of all targets. Finally, the model of A2255 slightly overpredicts the signal, which agrees with the prediction of the model luminosity, which is too high, as discussed above.

We also used the MILCA (Hurier et al. 2013) Compton parameter map obtained from Planck (Planck Collaboration XXII 2016) to compare our tSZ model to real data. We extracted a 1° × 1° sky patch centered on the individual clusters, and projected the MINOT Compton parameter map on the same grid for comparison. We smoothed the model with a 10 arcmin FWHM Gaussian beam to account for the angular resolution of Planck. In the right panel of Fig. 17 we show the data, the model, and the residual. The data and the model are shown on a log scale, and the residual is shown on a linear scale. The model and the data agree well for all three clusters. Nevertheless, a low excess in the central part of the A2142 model is evident, which can be explained by the fact that this cluster is slightly elongated, as is also seen in the ROSAT image.

In conclusion, we have compared the prediction from MINOT to the X-ray luminosity and tSZ fluxes, finding an overall good consistency between the two. The differences are likely explained by uncertainties in the model extrapolation. Similarly, MINOT is able with dedicated functions to predict the X-ray and tSZ images associated with a cluster model. The comparison to ROSAT and Planck maps has shown an overall good consistency for the targets tested here. In the context of modeling the nonthermal component of galaxy clusters, it is therefore possible to use MINOT together with X-ray or tSZ data to calibrate the thermal part of the model.

6.2. Comparison to Fermi-LAT constraints

γ-ray constraints have been obtained by Ackermann et al. (2014) for A1795, A2142 and A2255 as part of a larger sample using Fermi-LAT data (Fermi/LAT Collaboration 2009). For all clusters, they predicted the expected fraction of CRp pressure over the thermal pressure based on Pinzke & Pfrommer (2010) and Pinzke et al. (2011), given the mass of the clusters. They used these predictions together with a model for the spatial distribution of the CRs to compute flux prediction for these clusters. In order to compare our predictions to theirs, we first set the value of X_CRp(R₂₀₀) to that of Ackermann et al. (2014), taking into account the fact that the pressure ratio is twice the energy ratio (which we used for our parametrization). Our baseline CRp density model matches their baseline well (based on simulations, Pinzke & Pfrommer 2010) because the CRs are tied to the gas, as they neglect CR transport.

Fig. 18. Comparison between the cluster model flux predictions and the Fermi-LAT upper limit from Ackermann et al. (2014).

Table 7 shows that our predictions are lower than those of Ackermann et al. (2014) by a factor of about 3. However, we note that the masses of our selected clusters used in Ackermann et al. (2014) that were taken from Chen et al. (2007), are higher by a factor of 1.5−2.1 than ours. In order to account for this differences, we therefore rescaled our thermal energy according to self-similarity expectations as

$$\begin{array}{c}\hfill U_{\mathrm{th}}\to U_{\mathrm{th}}{\left(\frac{M_{500,\mathrm{this}\mathrm{work}}}{M_{500,\mathrm{Ackermann}2014}}\right)}^{-5/3}·\end{array}$$(48)

After we applied this rescaling, our fluxes agreed far better, within 20%. The differences may arise from the γ-ray production rate modeling, scatter in the thermal pressure, or differences in the thermal gas distribution that are not necessarily the same. To determine the effect of the latter, we also computed our fluxes by changing our density profiles to the best-fit β-model of the true density. We find that this change leads to differences in the γ-ray flux of up to 7% in the case of these clusters, which is significant.

In Fig. 18 we compute the energy-integrated flux as a function of energy in the case of our reference model (Sect. 2.7) and compare it to the upper limit set by Ackermann et al. (2014). While the predictions are relatively close to the upper limit, they remain below it for all three clusters.

In conclusion, we have shown that our model gives comparable predictions for the γ-ray flux compared to what is used in the literature. However, significant differences may arise as a result of the inner structure modeling of the clusters, which is generally ignored when large samples are used.

6.3. Comparison to radio observations

In this section, we quantitatively compare the radio predictions of our model to measurements available in the literature. To do so, we used the online database⁹, in which available radio data for many clusters are listed. All of our target clusters present diffuse radio emission: a radio mini-halo for A1795 (Giacintucci et al. 2014), a giant radio halo for A2142 (Giovannini & Feretti 2000; Venturi et al. 2017), and a giant radio halo plus a relic for A2255 (Kempner & Sarazin 2001; Govoni et al. 2005). In the case of A2142, we note that two components were distinguished in Venturi et al. (2017), and they likely arise from different physical processes. However, in this qualitative comparison, we only considered the global emission and sum the contribution from the two components.

Fig. 19. Comparison of our model prediction and the radio flux observed in our three clusters (left), and the effect of the radio models on the γ-ray constraints (right). Top: A1795. Middle: A2142. Bottom: A2255.

In Fig. 19 we compare our models (defined in Sect. 2.7) to the flux measured in the literature. In order to compute the model flux emission, we used an aperture radius that matched the signal from the respective articles (100, 500, and 930 kpc, for A1795, A2142, and A2255, respectively) and performed cylindrical integration of the synchrotron emission. For each cluster, we illustrate how it is possible to qualitatively change our model parameters to match the radio data, and also show how these changes translate into the γ-ray prediction and its comparison to the upper limits of Ackermann et al. (2014).

In the case of A1795, our default model underpredicts the radio flux by a factor of about 50%. Increasing the magnetic field by a factor of 1.4 would solve the difference without affecting the γ-ray prediction. Alternatively, the number of CRp might be inceased by a factor of about 1.5, but a the cost of increasing the γ-ray prediction by a similar amount and approaching the Fermi-LAT limit. Finally, increasing the number of CRe₁ by a factor of 3 would also solve the difference, with changes in the γ-ray prediction only at energies below 0.1 GeV, which is barely accessible for Fermi-LAT.

In the case of A2142, our default model overpredicts the radio emission by almost an order of magnitude, and our spectrum appears too flat compared to the data. First, we considered a purely hadronic scenario (i.e., no CRe₁), for which a slope of the CRp spectrum needs to be set to about 2.9 and the amplitude reduced by a factor of two in order to match the data. Another option is to consider only CRe₁, in which case the spectrum slope should be set to about 3.3 and the normalization increased by a factor of 500 to match the data. A combination of the two scenarii might also be used, especially because the radio emission presents two distinct components (see Venturi et al. 2017, for discussions). In both scenarii, the predicted γ-ray emission also decreases, and thus remains in agreement with the Fermi-LAT limit.

In the case of A2255, we observe a disagreement between the 1400 GHz flux from Kempner & Sarazin (2001) and Govoni et al. (2005). For our purpose, we chose to use the value of Govoni et al. (2005) as a reference (flux of the halo alone, without the relic). The amplitude of our default model broadly agrees with the observation, but our spectral index is too low. In the purely hadronic scenario, increasing the slope of the CRp spectrum to 2.9 together with increasing its normalization by a factor of 10 would bring agreement between the model and the data. Alternatively, in a model with only CRe₁, we would need to increase the CR spectrum slope to 3.1 and multiply the normalization by 500 in order to reach agreement. In the two cases, the γ-ray flux remains below the Fermi-LAT upper limit.

In this section, we showed that we can change the model parameters to match the radio data in a qualitative way. We focused on the slope and normalization of the CR content of the cluster, but opening the parameter space to spatial distributions or a functional form of the spectra might also play an important role. In all the considered cases, it is not possible to rule out any model using γ-ray data because the limits remain too high. However, purely hadronic model predictions are just a factor of a few below the Fermi-LAT limits in the case of these clusters.

7. Conclusions and summary

We have provided an exhaustive description of MINOT, a new software dedicated to the modeling the nonthermal components of galaxy clusters and predicting associated observables. While the software was originally developed to describe the γ-ray emission from galaxy clusters, MINOT also accounts for most of the emission associated with the diffuse ICM component: X-rays from thermal bremsstrahlung, tSZ signal in the millimeter band, γ-rays and neutrino emission form hadronic processes, γ-rays from inverse Compton emission, and radio synchrotron emission. Because the γ-ray emission is connected to the same underlying cluster physical properties as these other observables, MNOS provides a useful self-consistent modeling of the signal, and these additional observables can be used, for example, to calibrate a γ-ray model. However, MINOT can also be used to independently model observables in the different bands.

The software is made publicly available¹⁰, and this paper aims at providing a reference for any user of the code. To this aim, we have discussed the structure of the code and the interdependencies of the different modules in Sect. 2, while Sects. 3–5 provided details about the physical processes considered, how they are accounted for, and the way observables are computed. The different functions were illustrated with the use of a reference cluster model. It allowed us to show the dependence of each wavelength on the physical properties of the cluster. In Sect. 6 we also compared the predictions from MINOT to data available in the literature in order to show how the code can be used to model real data. We used Planck and ROSAT data for the thermal component, and Fermi-LAT plus various radio data for the nonthermal component. Finally, we note that the MINOT code is well documented, and many examples are provided in the public repository.

The different assumptions made in the code were discussed. In particular, the modeling relies on primary base quantities that are used to derive secondary properties of the cluster and generate observables under the assumption of spherical symmetry. The primary quantities are the thermal electron pressure and density profiles, the CRe₁ and CRp profiles and spectra, and the profile of the magnetic field strength. Regarding the CRe₂, they are processed assuming no diffusion in a steady-state scenario. However, other electron populations can be accounted for using the CRe₁. In order to set the base physical properties of the cluster, several predefined models are available, but it is also possible to provide any user-defined quantity.

The accuracy of the modeling was addressed. Regarding the thermal component, modeling uncertainties associated with X-rays are below the percent level, and the uncertainties of the tSZ signal are at the percent level when a relativistic correction at high temperature is considered and much smaller otherwise. The hadronic processes (γ-rays, neutrinos, and secondary electrons), on the other hand, present uncertainties at a level of typically 25% above the considered energy range when the latest models available in the literature are considered. In addition, we stress that the effect of helium is about 50% of the signal and should be accounted for (as done in MINOT). The computation of inverse Compton and synchrotron emission is based on analytical approximations whose precision is expected to remain within 1% and 0.2%, respectively. Nevertheless, we stress that the main limitations of the modeling is not the accuracy of the computation, but the underlying assumptions discussed above. In particular the use of spherical symmetry and the assumption of stationarity to compute the distribution of secondary electrons are likely to be the dominant sources of mismodeling.

The MINOT software can be used for a wide variety of operations, and we list just a few examples here:

• The joint modeling of the nonthermal emission in galaxy clusters for which detailed multi-wavelength data are available. The parameters of the model can be fit jointly to such data to constrain different scenarios while accounting for uncertainties in the different components.

• The prediction of the expected signal, based on ancillary data, for observation proposals. For instance, it is possible to predict the tSZ emission associated with that of an X-ray observed cluster, assuming that the pressure follows a universal profile.

• The prediction of the background CR induced γ-rays in the context of dark matter searches. MINOT provides an easy way to model the CR background, which needs to be marginalized over to obtain constraints on the nature of dark matter.

• The simulation of sky maps associated with the observables considered here, given a halo catalog.

• Pedagogical purposes. Because it includes most of the ICM associated processes, MINOT can be used to understand the effect of some given parameters on the observable.

Historically, the understanding of the physical properties of galaxy clusters has strongly benefited from multiwavelength observations and analysis. With the current and upcoming facilities aiming at exploring the nonthermal component of galaxy clusters, in particular in the radio and γ-ray bands, it has become very useful to have an easy-to-use self-consistent modeling software that allows us to predict the expected signal based on some assumptions. The MINOT software provides such a tool.

Acknowledgments

We are thankful to the anonymous referee for useful comments that helped improve the quality of the paper. We would like to thank R. Alves Batista for useful comments. The work of JPR and MASC was supported by the Spanish Agencia Estatal de Investigación through the grants PGC2018-095161-B-I00 and IFT Centro de Excelencia Severo Ochoa SEV-2016-0597, and the Atracción de Talento contract no. 2016-T1/TIC-1542 granted by the Comunidad de Madrid in Spain. We acknowledge the use of HEALPix (Górski et al. 2005). 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). Several modules of MINOT are based on the Naima software (Zabalza 2015).

References

All Tables

All Figures

	Fig. 1. Overview of the MINOT input modeling, the considered physical processes at play in the ICM, and the observables that are computed. The interdependences are shown by the black arrows.
In the text

	Fig. 2. Overview of the code structure and interfaces.
In the text

Fig. 3. Left: illustration of the different radial profiles available in the library. The β-model parameters are (n0, rc, β) = (1, 50 kpc, 0.7); the SVM model parameters are (n0, rc, β, rs, γ, ϵ, α) = (1, 50 kpc, 0.7, 1500 kpc, 3, 4, 0.5); the GNFW model parameters are (P0, rp, a, b, c) = (1, 200 kpc, 1.5, 4, 0.3); and the double β-model parameters are (n01, rc1, β1, (n02, rc2, β2) = (1, 50 kpc, 0.7, 5, 5 kpc, 2). Right: illustration of the different spectral models available in the library. The index is set to α = 2.3 for all models, and the cutoff or break energy to Ebreak/cut = 1 GeV. In both cases, the “User” model is not shown, but allows passing any arbitrary function that is interpolated in the code.

In the text

Fig. 4. Thermodynamic properties of our cluster sample defined in Table 4, as computed based on the modeling of the thermal electron pressure and density (i.e., the base quantities used to derive the others). The truncation radius is shown at r = 5 Mpc, where all distributions drop, except for the integrated ones that remain constant afterward. For the entropy profile, we also show the expectation for purely gravitational collapse (Voit 2005). For the gas fraction, the ratio between the mean baryon density and the mean matter density is shown as a dashed line. We also report the radius at which a density contrast of 500 (i.e., R500) is reached in the overdensity profile.

In the text

	Fig. 5. CR properties of the baseline cluster model. Left: CR number density profile (i.e., CR distribution integrated between Emin and Emax, as indicated in the legend) for CRe1 and CRp. Right: CR spectrum integrated over the radius up to a maximum radius, as given in the legend, for CRe1 and CRp.
In the text

Fig. 6. Magnetic field and CRp-to-thermal energy properties of the baseline cluster model. Left: magnetic field profile obtained by assuming a fixed magnetic field strength of 5 μG at 100 kpc, and a scaling relative to the thermal density as $B\propto n_{\mathrm{e}}^{0.5}$. Right: CRp-to-thermal energy profile for different models of the CRp distribution. The color lines show each real cluster, using the same colors as in the left panel, and with nCRp ∝ ngas. The black lines correspond to variation in the CRp scaling with respect to the thermal gas for the baseline cluster model, as indicated in the legend.

In the text

Fig. 7. Production rate of secondary particles from hadronic processes. Left: relative comparison of the production rate of γ-rays in the case of no helium and zero metallicity in the Kafexhiu et al. (2014) and the Kelner et al. (2006) parameterization, using the Pythia8 parameterization from Kafexhiu et al. (2014) as a reference. For Kelner et al. (2006), we also show the low- and high-energy (δ-approximation) limits as dashed lines, which were combined to compute the rate over the total energy range. Middle: effect of including helium and metals in the model. This was done with Pythia8 using the Kafexhiu et al. (2014) parameterization. Right: comparison of the production rate for γ-rays, electrons, positrons, and neutrinos, computed for a radius r = 100 kpc. We also include the power-law function proportional to E−2.4 for comparison because it corresponds to the injected CRp distribution.

In the text

	Fig. 8. Contributions to the energy loss rate as a function of energy for a radius r = 10 kpc and r = 1000 kpc.
In the text

	Fig. 9. Left: profile of the CRe2, taken at 10 MeV and 10 TeV (after applying a multiplicative factor of 1015). Right: spectrum of the secondary electron at 100 and 1000 kpc.
In the text

Fig. 10. Observables associated with the X-ray emission. Left: X-ray spectrum within R500 shown for both MEKAL and APEC models, as well as with and without galactic photoelectric absorption. The MEKAL model is difficult to distinguish because it coincides very well with the APEC model. Right: surface brightness profile in the band 0.1−2.4 keV. The dynamical range of the profile amplitude has been set to the same value for all observables.

In the text

	Fig. 11. Observables associated with the tSZ effect. Left: tSZ spectrum within R500. We also provide the spectrum in the case when relativistic corrections are neglected for illustration. Right: compton parameter profile. The dynamical range of the profile amplitude has been set to the same value for all observables.
In the text

	Fig. 12. Effect of the EBL on the normalized spectra of a baseline model cluster at z = 0.02. All the available models are shown as indicated in the legend. Because of the redshift stretching, the shape of the γ-ray spectrum slightly changes even in the case without EBL absorption (“none”).
In the text

	Fig. 13. Observables associated with the γ-ray hadronic emission. The signal coming from inverse Compton interactions is also shown for comparison. Left: γ-ray spectrum within R500. Right: γ-ray profile integrated between 1 GeV and 1 TeV. The dynamical range of the profile amplitude has been set to the same value for all observables.
In the text

	Fig. 14. Observables associated with neutrino hadronic emission. Left: neutrino spectrum within R500. Right: neutrino profile integrated between 1 GeV and 1 TeV. The dynamical range of the profile amplitude has been set to the same value for all observables.
In the text

	Fig. 15. Observables associated with inverse Compton emission. Left: inverse Compton spectrum within R500, with the contribution from both primary and secondary electrons. Right: inverse Compton profile integrated between 1 GeV and 1 TeV. The dynamical range of the profile amplitude has been set to the same value for all observables.
In the text

	Fig. 16. Observables associated with the synchrotron emission. Left: synchrotron spectrum within R500, with the contribution from both primary and secondary electrons. Right: synchrotron emission profile at 100 MHz. The dynamical range of the profile amplitude has been set to the same value for all observables.
In the text

Fig. 17. Comparison between the thermal gas observables and the MINOT model prediction. Left: comparison between ROSAT X-ray images and the MINOT model prediction. The scale is linear from −10 to 10 s−1 deg−2 and logarithmic above. Contours are as ±5 × 2N s−1 deg−2, with N the contour index starting at 0. Right: comparison between the Planck Compton parameter MILCA image and the MINOT model prediction. The contours are multiples of 3σ, where σ is the rms of the residual map. The Planck 10 arcmin FWHM beam is shown as a filled white circle. All maps where smoothed (in addition the intrinsic resolution) with a 2 arcmin FWHM gaussian beam for visual purpose. All the sky patches shown here are 1° × 1°.

In the text

	Fig. 18. Comparison between the cluster model flux predictions and the Fermi-LAT upper limit from Ackermann et al. (2014).
In the text

	Fig. 19. Comparison of our model prediction and the radio flux observed in our three clusters (left), and the effect of the radio models on the γ-ray constraints (right). Top: A1795. Middle: A2142. Bottom: A2255.
In the text