# Prospects for γ-ray observations of the Perseus galaxy cluster with the Cherenkov Telescope Array

K Abe, S Abe, Fabio Acero, A Acharyya, R Adam, A Aguasca-Cabot, I Agudo, A Aguirre-Santaella, J Alfaro, R Alfaro, et al.

► To cite this version:

K Abe, S Abe, Fabio Acero, A Acharyya, R Adam, et al.. Prospects for γ-ray observations of the Perseus galaxy cluster with the Cherenkov Telescope Array. JCAP, 2024, 10, pp.004. 10.1088/1475-7516/2024/10/004. hal-04212511

HAL Id: hal-04212511

https://hal.science/hal-04212511v1

Submitted on 23 May 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
PAPER • OPEN ACCESS

# Prospects for γ-ray observations of the Perseus galaxy cluster with the Cherenkov Telescope Array

To cite this article: K. Abe et al JCAP10(2024)004

View the [article online](#) for updates and enhancements.

You may also like

- Performance of the CMS high-level trigger during LHC Run 2

A. Hayrapetyan, A. Tumasyan, W. Adam et al.

- Multi-messenger Observations of a Binary Neutron Star Merger

B. P. Abbott, R. Abbott, T. D. Abbott et al.

- Identification of low-momentum muons in the CMS detector using multivariate techniques in proton-proton collisions at sqrt(s) = 13.6 TeV

V. Chekhovsky, A. Hayrapetyan, V. Makarenko et al.
# Prospects for γ-ray observations of the Perseus galaxy cluster with the Cherenkov Telescope Array

<figure><img src="image_2.png" /></figure>

cherenkov
telescope
array

## The CTA consortium

K. Abe et al.

Full author list at the end of the paper

E-mail: remi.adam@oca.eu, skerzot@ciencias.unam.mx, huetten@icrr.u-tokyo.ac.jp, judit.perez@ung.si, miguel.sanchezconde@uam.es

**ABSTRACT:** Galaxy clusters are expected to be both dark matter (DM) reservoirs and storage rooms for the cosmic-ray protons (CRp) that accumulate along the cluster's formation history. Accordingly, they are excellent targets to search for signals of DM annihilation and decay at γ-ray energies and are predicted to be sources of large-scale γ-ray emission due to hadronic interactions in the intracluster medium (ICM). In this paper, we estimate the sensitivity of the Cherenkov Telescope Array (CTA) to detect diffuse γ-ray emission from the Perseus galaxy cluster. We first perform a detailed spatial and spectral modelling of the expected signal for both the DM and the CRp components. For each case, we compute the expected CTA sensitivity accounting for the CTA instrument response functions. The CTA observing strategy of the Perseus cluster is also discussed. In the absence of a diffuse signal (non-detection), CTA should constrain the CRp to thermal energy ratio $X_{500}$ within the characteristic radius $R_{500}$ down to about $X_{500} < 3 \times 10^{-3}$, for a spatial CRp distribution that follows the thermal gas and a CRp spectral index $\alpha_{CRp} = 2.3$. Under the optimistic assumption of a pure hadronic origin of the Perseus radio mini-halo and depending on the assumed magnetic field profile, CTA should measure $\alpha_{CRp}$ down to about $\Delta\alpha_{CRp} \simeq 0.1$ and the CRp spatial distribution with 10% precision, respectively. Regarding DM, CTA should improve the current ground-based γ-ray DM limits from clusters observations on the velocity-averaged annihilation cross-section by a factor of up to ~5, depending on the modelling of DM halo substructure. In the case of decay of DM particles, CTA will explore a new region of the parameter space, reaching models with $\tau_{\chi} > 10^{27}$ s for DM masses above 1 TeV. These constraints will provide unprecedented sensitivity to the physics of both CRp acceleration and transport at cluster scale and to TeV DM particle models, especially in the decay scenario.

**KEYWORDS:** cosmic ray experiments, dark matter experiments, galaxy clusters, gamma ray experiments

ARXIV ePRINT: 2309.03712

<figure><img src="image_3.png" /></figure>
# Contents

<table><thead><tr><th></th><th></th></tr></thead><tbody><tr><td>1</td><td>Introduction</td><td>2</td></tr><tr><td>1.1</td><td>Diffuse &gamma;-ray emission from cosmic-rays and dark matter in galaxy clusters</td><td>2</td></tr><tr><td>1.2</td><td>The Perseus cluster as a promising &gamma;-ray target</td><td>5</td></tr><tr><td>2</td><td>Modelling the &gamma;-ray emission associated with cosmic-rays</td><td>6</td></tr><tr><td>2.1</td><td>Modelling the intracluster medium components</td><td>7</td></tr><tr><td>2.2</td><td>Non-thermal radiative processes</td><td>11</td></tr><tr><td>2.3</td><td>Calibration of the model parameters</td><td>12</td></tr><tr><td>3</td><td>Modelling the &gamma;-ray emission associated with dark matter</td><td>16</td></tr><tr><td>3.1</td><td>Perseus main halo</td><td>18</td></tr><tr><td>3.2</td><td>The role of substructures</td><td>19</td></tr><tr><td>3.3</td><td>Dark matter annihilation and decay fluxes</td><td>22</td></tr><tr><td>4</td><td>Observation setup</td><td>24</td></tr><tr><td>4.1</td><td>CTA pointing configuration</td><td>24</td></tr><tr><td>4.2</td><td>Background sky</td><td>26</td></tr><tr><td>4.3</td><td>Towards an optimal observing strategy</td><td>28</td></tr><tr><td>5</td><td>CTA sensitivity to CR induced &gamma;-ray emission from Perseus</td><td>31</td></tr><tr><td>5.1</td><td>Data preparation</td><td>31</td></tr><tr><td>5.2</td><td>Sky components and degeneracy</td><td>32</td></tr><tr><td>5.3</td><td>Probing the parameter space with CTA</td><td>34</td></tr><tr><td>5.4</td><td>Measuring the cosmic-ray properties in the case of detection</td><td>36</td></tr><tr><td>5.5</td><td>Discussions</td><td>42</td></tr><tr><td>6</td><td>CTA sensitivity to DM induced &gamma;-ray emission from Perseus</td><td>46</td></tr><tr><td>6.1</td><td>Template-fitting analysis and DM sensitivity</td><td>48</td></tr><tr><td>6.2</td><td>Projected sensitivity to annihilating DM</td><td>49</td></tr><tr><td>6.3</td><td>Projected sensitivity to decaying DM</td><td>52</td></tr><tr><td>7</td><td>Summary and conclusions</td><td>53</td></tr><tr><td>A</td><td>Validation of the thermal model using Planck data</td><td>58</td></tr><tr><td>B</td><td>Impact of the magnetic field on the model calibration</td><td>58</td></tr><tr><td>C</td><td>Fermi-LAT sources around the Perseus cluster</td><td>60</td></tr><tr><td>D</td><td>Galactic foreground estimate</td><td>61</td></tr><tr><td>E</td><td>Convergence of the CR parameter constraints versus the number of simulations</td><td>62</td></tr></tbody></table>
<table><tr><td>F</td><td>Validation of the cosmic-ray analysis with gammapy</td><td>63</td></tr><tr><td>G</td><td>Cosmic-ray parameter constraints in the pure hadronic model</td><td>63</td></tr><tr><td>H</td><td>Interplay between the DM components and the astrophysical γ-ray sources</td><td>64</td></tr><tr><td>I</td><td>CTA sensitivity to DM under the ON-OFF observational setup</td><td>66</td></tr><tr><td>J</td><td>Convergence of the DM fits versus the number of simulations</td><td>70</td></tr><tr><td>K</td><td>Recovered astrophysical parameters from the DM template-fitting</td><td>72</td></tr><tr><td>L</td><td>CTA sensitivity to DM in Perseus with ctools</td><td>73</td></tr><tr><td>L.1</td><td>Data preparation</td><td>73</td></tr><tr><td>L.2</td><td>DM analysis pipeline with ctools</td><td>73</td></tr><tr><td>L.3</td><td>CTA sensitivity to DM under the ON-OFF observational setup with ctools</td><td>74</td></tr><tr><td>L.4</td><td>CTA sensitivity to DM based on template fitting with ctools</td><td>74</td></tr><tr><td>L.5</td><td>Comparison with gammapy</td><td>79</td></tr><tr><td>The CTA consortium</td><td></td><td>96</td></tr></table>

# 1 Introduction

## 1.1 Diffuse γ-ray emission from cosmic-rays and dark matter in galaxy clusters

Clusters of galaxies are the largest virialized structures in the Universe, with masses up to about $10^{15}$ M$_{\odot}$. They are dominated by dark matter (DM; ~85% in mass) and permeated by the intracluster medium (ICM; ~10–15% in mass), whose physical properties are shaped by the hierarchical growth of structures through the merging of subclusters and the smooth accretion of surrounding matter [1, 2]. While the ICM is mostly thermal, these energetic merging events do not only dissipate the kinetic energy into heat via shock waves and turbulence, but may also accelerate cosmic-rays (CR) in the ambient magnetic field [3]. Galaxies in galaxy clusters account only up to a few percent of the total mass, yet they can also directly inject CR via active galactic nuclei (AGN) feedback or star formation activity [4–6].

Direct evidence for the presence of CR electrons (CRe) and magnetic fields have been found in a growing number of galaxy clusters, thanks to the radio observations of diffuse synchrotron emission [7–9]. These sources are classified as radio halos (that roughly follow the thermal ICM) and radio relics (elongated and peripheral [10]). Radio halos are further classified as giant radio halos (~Mpc size), associated with cluster mergers [11], suggesting that they are powered by the energy dissipated during these events, and mini-halos in more relaxed clusters, which extend on 100–300 kpc scales and are generally confined in the core of cool-core clusters [12, 13]. Recent observations have made the phenomenology more complex, showing that in a number of cases mini-halos are surrounded by larger scale (usually very steep spectrum) emission, similar to giant radio halos [14–16]. Furthermore, it has been discovered that radio emission may extend on scales larger than those of giant radio halos,
in the form of radio bridges [17, 18] and mega-halos [19]. These emerging evidences from
observations suggest a more complex picture of non-thermal phenomena in galaxy clusters
that may require a revision of current classification schemes. Galaxy clusters are also expected
to act as storehouses for the CR protons (CRp) and heavier nuclei [20] due to their long
lifetimes, once they are injected in the ICM via several mechanisms (cosmic shocks and
galaxies via supernovae explosions, starbursts or AGN) [3, 21, 22]. These CRp should interact
hadronically in the ICM to produce γ-ray emission [23–28]. Such interactions imply the
production of high-energy secondary electrons, which will eventually contribute to the cluster-
scale radio emission. According to the current theoretical picture for radio halos and relics,
the emitting electrons are (re)accelerated by turbulence and shocks, respectively [29–34].
The large volumes that are probed by radio halos in which the ICM is tenuous disfavour an
important contribution from secondary electrons in these sources [3], although secondaries
can contribute to the population of the seed particles to reaccelerate [33, 35, 36]. On the other
hand, pure hadronic models may still explain the smaller mini-halos in dense cores [37, 38],
although a number of evidences suggests that gas motions may play a major role [39, 40].
Cosmological numerical simulations of CR in clusters have obtained quantitative predictions
for the expected γ-ray emission [41–43], which has proved useful when searching for clusters
in γ-ray observations [44]. Nevertheless, many uncertainties related to acceleration and
transport physics affect the expected γ-ray signal [3, 22, 45] and cosmological simulations
including full turbulent reacceleration physics have not been obtained yet. In this context
the study of γ-ray emission from galaxy clusters plays a central role.

In the past, galaxy clusters have provided strong gravitational evidence in favour of the
existence of DM [46–48]. Thus, being DM-dominated, these objects also represent natural
astrophysical targets for current DM search efforts. Since the underlying nature of DM
(and thus its potential signatures) is still unknown [49], galaxy clusters have been used to
probe the properties of the DM particle with a variety of techniques (e.g., [50–53]). One of
the most promising ones is the search for DM-induced γ-ray signals ([54–56], for reviews),
expected from the annihilation or decay of the so-called Weakly Interacting Massive Particles
(WIMPs), one of the most studied DM particle candidate. WIMPs (e.g., [57–60]) would be
produced thermally in the Universe via the “freeze-out” mechanism and would have masses
$O(0.1–100)$ TeV.<fnref n="1" /> They can arise from several theoretical frameworks, ranging from minimal
extensions of the Standard Model [57, 63, 64] to extra-dimensions [65, 66] and others [67–70].
The expected γ-ray flux from their annihilation or decay in astrophysical objects mainly
depends on the target DM density (squared, for annihilation) and its distance to Earth.
Thus, for DM decay, local galaxy clusters can yield the highest expected fluxes compared
to other possible targets, as they are the most massive structures in the Universe. As for
DM annihilation, clusters can provide fluxes comparable to the ones from dwarf spheroidal
galaxies (dSphs, [71]), as long as the DM interactions expected in their substructures are
taken into account [72–74]. These substructures, usually referred to as subhalos, are a natural
result of the ΛCDM hierarchical growth of structures [75–77], and their abundance is expected
to be comparatively significant in clusters, as they are the largest exponents of structure

<footnote n="1">¹There is an on-going discussion within the community with respect to the viability of these models beyond hundreds of TeV (see e.g., [61, 62]).</footnote>
formation at present. Despite the optimal characteristics of galaxy clusters to be used for
γ-ray DM searches, the main drawback with respect to other targets is the predicted γ-ray
emission from more conventional astrophysical processes. Indeed, the expected γ-rays from
hadronic interactions of the CRp in the ICM can act as a complex background to search
for a DM-induced signal using standard analysis techniques.

The search for diffuse γ-ray emission from galaxy clusters has been going on for over two
decades, both using space-based observations in the GeV band [44, 78–86] and ground-based
observations at energies > 100 GeV [87–91]. Yet, such signal remained elusive.<fnref n="2" /> Nonetheless,
these non-detections allowed to constrain the CRp to thermal energy ratio down to few
percent [95], challenging the understanding of diffusive shock acceleration in the ICM when
combined with radio data [96–98], although the large modelling uncertainty affects the
predictions [21]. The stringent γ-ray limits were also used to disfavour hadronic models for
nearby radio halos [99, 100] and to test models based on the reacceleration of secondary
particles [93, 100].

The non-detection of γ-ray emission from clusters is in agreement with the lack of DM-
induced γ-ray signals from other promising astrophysical targets, especially dSphs [101–103].
DM searches in clusters mainly targeted very massive and local objects [71, 104, 105], since
the expected flux is proportional to the mass of the objects and decreases with the distance
squared. Yet, the lack of a DM signal in clusters has allowed to provide also strong constraints
for annihilating WIMPs from a combined analysis of various clusters [106–113] or single
cluster observations [84, 114, 115]. These DM annihilation limits are nevertheless not at
the level of discarding thermal WIMP models.<fnref n="3" /> The situation comparatively improves for
WIMP DM decay. Indeed, lower limits on the WIMP lifetime derived from the observation of
clusters [107, 116, 117] are among the most constraining ones at present. Some authors have
also studied the possibility that some past hints of detection of γ-rays in clusters were due
to WIMP DM [108, 118, 119], however these works were inconclusive due to the faintness
of such signals, indeed never confirmed.

The Cherenkov Telescope Array<fnref n="4" /> (CTA, [120]) will be the next generation ground-based
γ-ray observatory. It will be amongst the most sensitive γ-ray telescope from 20 GeV to
300 TeV. CTA will be based at two sites: La Palma, in the North, and Paranal in the South,
allowing us to observe sources in a large fraction of the sky. CTA will provide a major
improvement, up to one order of magnitude in sensitivity and up to a factor 2 in angular
resolution, with respect to previous instruments.<fnref n="5" /> The study of CR physics and DM are
among the main drivers of CTA science [121]. In particular, CTA will allow us to search for
diffuse γ-ray emission from galaxy clusters. One of the proposed CTA key science projects
consists in the observation of the Perseus galaxy cluster, which is among the most promising

<footnote n="2">²Note that [92] claimed the detection of the Coma cluster with *Fermi*-LAT. [93] confirmed the signal detection and showed that it would imply a CRp to thermal energy ratio of about 1%. However, they also noted that possible confusion with point sources could not be excluded, so that the first non-ambiguous detection of diffuse γ-ray emission from galaxy clusters is still to be achieved. See also [94] who claimed the detection of extended emission.</footnote>

<footnote n="3">³This is true under the assumption of “vanilla” WIMP models.</footnote>

<footnote n="4">⁴https://www.cta-observatory.org/.</footnote>

<footnote n="5">⁵https://www.cta-observatory.org/science/ctao-performance/.</footnote>
targets for such observations. These observations should complement DM searches in the Galactic center [122], dwarf galaxies, the Large Magellanic Cloud [123] or the search of axion-like particles, and should be used to probe fundamental physics [124].

## 1.2 The Perseus cluster as a promising γ-ray target

The Perseus cluster (Abell 426) is the brightest cluster in the X-ray sky [125]. It is one of the most massive nearby clusters and presents the typical properties of a relaxed cool-core cluster with a dense core. Nonetheless, two main cold fronts have been identified and interpreted as the result from the sloshing due to minor mergers [126, 127]. The Perseus cluster hosts a radio mini-halo [128–136], and X-ray cavities associated with the radio lobes of the central AGN, NGC 1275 (3C84), indicate that the feedback is important in the cluster center [137, 138]. AGN activity may also be responsible for weak shocks and turbulence [139, 140], which could reaccelerate particles, in addition to direct CR injection from the AGN. The contribution to the radio mini-halo from hadronic interactions, direct CR injection from AGN and the role of turbulence remains unknown.

Two AGN from the Perseus cluster are known γ-ray emitters: NGC 1275 [141] and IC 310 [142]. Both sources are variable in time. NGC 1275 is the central galaxy of the cluster and IC 310, located about 0.6 deg southwest from the X-ray peak, is consistent with a narrow-angle tail radio galaxy infalling into the cluster [135]. A few other remarkable radio galaxies are present in the cluster volume: NGC 1265, NGC 1272, CR 15 [135], that were not detected at γ-ray energies.<fnref n="6" /> While these sources are expected to contaminate the CTA data when searching for diffuse emission, they are also contributing to inject CR into the ICM, which could eventually contribute to large-scale γ-ray emission.

The Perseus cluster has been recognised as one of the best targets for searches of CR-induced γ-ray emission [43, 88]. This is because the mini-halo traces a dense region where hadronic collisions and the production of secondaries is maximised. While its central γ-ray bright galaxy, NGC 1275, prohibits reliable constraints on the diffuse ICM component from Fermi-LAT, the better angular resolution of CTA and the larger accessible energy range probed is expected to allow us separating the different sources of a possible γ-ray emission. In fact, the expected mild angular extent of the cluster (about 1 deg) due to its proximity (about 75 Mpc) implies that CTA is expected to resolve the diffuse emission if bright enough, but also that the angular extension of the signal is smaller than the field of view diameter by a factor of 5 to 10, so that systematic effects associated with the background modelling are limited. Along this line, also due to its large mass and proximity, Perseus stands out as one of the best clusters to search for DM-induced γ-ray emission. Indeed, the expected annihilation/decay flux is comparable to the one from other promising local galaxy clusters [143], such as Coma, Fornax, Ophiuchus, Hydra or Centaurus — see [71, 105]<fnref n="7" /> as well as other traditional DM targets such as dwarf satellite galaxies or nearby galaxies.

<footnote n="6">⁶But they could be included in the sky model as point sources and their contribution accounted for.</footnote>

<footnote n="7">⁷The exception is Virgo. The Virgo cluster [144] is sometimes considered as the best cluster for DM searches since it is the closest one to Earth. However, its large angular extension, comparable to the field of view of CTA, and its dynamical condition as a not-yet-virialized object, complicate considerably a potential DM analysis.</footnote>
A previous search for diffuse TeV γ-ray emission towards the Perseus cluster was performed using the MAGIC telescopes [88, 89, 91]. Not having detected the signal, they reported upper limits on the CR to thermal pressure ratio assuming different models. For instance, using a spectral photon index slope of 2.2 and a relatively compact profile for the CRp, they obtained an upper limit of ~1–2% on this ratio. This provided the best limit on the CR content of a cluster obtained from ground-based γ-ray observations so far, at a similar level to that obtained with Fermi-LAT [44].<fnref n="8" /> Fermi-LAT data from Perseus have been recently used to obtain constraints on the velocity-averaged DM annihilation cross-section as well [112]. The obtained limits are more than two orders of magnitude above from the reference value of the thermal relic cross-section. Additionally, MAGIC observations were used to set constraints on the DM decay lifetime [117], these being among the strongest constraints for DM masses in the TeV, range up to date reaching the value of $2 \times 10^{26}$ s.

Finally, the location of the Perseus cluster in the sky allows for low zenith-angle observations from the CTA Northern Array, guaranteeing the best sensitivity of the array over its whole energy range. For all these reasons, the Perseus cluster was selected as the prime target for diffuse γ-ray emission searches from galaxy clusters with CTA, as one of the Key Science Projects (KSP) [121]. We refer to this previous work for further details about this choice.

The paper is organized as follows. Section 2 and section 3 present the cluster modelling of the cluster in the context of CR and DM induced γ-ray emission, respectively. The observation setup and the background sky modelling are discussed in section 4. Section 5 and section 6 provide the results on the CTA sensitivity to CR and DM physics, respectively. We conclude in section 7. A few appendices complement the paper.

Throughout this paper, we assume a flat $\Lambda$CDM cosmology with $H_0 = 70 \text{ km s}^{-1} \text{ Mpc}^{-1}$, $\Omega_M = 0.3$, and $\Omega_\Lambda = 0.7$. The coordinates of the Perseus cluster center are taken as the ones of its central galaxy NGC 1275, R.A., Dec = 49.9507, 41.5117 deg<fnref n="9" /> and its redshift is $z = 0.017284$ [146], corresponding to a luminosity distance of $d_L \approx 75 \text{ Mpc}$. Given our reference cosmological model, 1 deg in the sky corresponds to a 1.265 Mpc distance at the redshift of Perseus. We adopt the characteristic angular radius $\theta_{500} = 59.7 \pm 0.4 \text{ arcmin}^{10}$ obtained by [147] using the fit of the Planck universal pressure profile.<fnref n="11" /> This corresponds to a physical radius $R_{500} = (1.26 \pm 0.01) \text{ Mpc}$ and to a mass $M_{500} = (5.77 \pm 0.12) \times 10^{14} \text{ M}_\odot$. In the paper, International System of Units is used unless specified.

## 2 Modelling the γ-ray emission associated with cosmic-rays

The prediction of the diffuse γ-ray emission induced by CR in the ICM requires the detailed modelling of the physical components at play, their interactions, and the underlying radiative

<footnote n="8">⁸Note that a high-significance image of the Perseus diffuse γ-ray emission was reported in [145] using the SHALON telescopes. However, their flux is in strong disagreement with the limits allowed by the MAGIC observations so that these results are very controversial.</footnote>

<footnote n="9">⁹NASA/IPAC extragalactic database, http://ned.ipac.caltech.edu/.</footnote>

<footnote n="10">¹⁰$R_{500}$ is the radius within which the mean cluster density is 500 times the critical density of the Universe at the cluster's redshift, and $M_{500}$ is the mass within this radius. $\theta_{500}$ is the angle on the sky that corresponds to $R_{500}$. Quantities with subscript 500 (or 200) refers to quantities within $R_{500}$ (or $R_{200}$).</footnote>

<footnote n="11">¹¹This value is in excellent agreement with the one derived from our thermal model and the assumption of hydrostatic equilibrium, assuming an hydrostatic mass bias of 0.1, in section 2.</footnote>
processes. In this section, we describe such a model based on the MINOT software [148].<fnref n="12" /> MINOT is dedicated to compute the observable properties of the ICM (radio synchrotron, thermal Sunyaev-Zel'dovich effect, X-ray thermal bremsstrahlung, inverse-Compton, γ-rays from hadronic interactions, and neutrino emission) given the user-defined physical state of the cluster. The code describes galaxy clusters as spherically symmetric objects. Here, we discuss essentially the calibration of the input model using external data assuming different scenarios, which we feed to MINOT. This includes the thermal gas density and pressure (section 2.1.1), the magnetic field strength (section 2.1.2), and the CR spatial and spectral distributions (section 2.1.3). This allows us to perform predictions for the γ-ray signal and estimate model uncertainties using the method implemented in MINOT.

## 2.1 Modelling the intracluster medium components

We model the Perseus cluster assuming spherical symmetry, considering only radial profiles to describe its ICM components. This assumption is expected to be fairly accurate given the fact that the Perseus cluster is overall a relaxed system.

While clusters are diffuse objects with no clear definitions of their extension, we expect density and pressure discontinuities near the virial radius [149, 150]. We consider a maximum radial extent of the cluster and define this truncation radius as $R_{tr} = 3R_{500} = 3.78$ Mpc $\simeq 2R_{200}$ using the measurement from [151] as a reference. The exact choice of this value does not significantly affect our results, but this allows us to perform numerical integrations over a well defined region.

### 2.1.1 Thermal gas

The modelling of the thermal gas is necessary to compute the CR induced γ-ray emission for two reasons: 1) the nuclei (protons, helium and heavier elements) are involved in hadronic interactions that lead to the emission of γ-rays and secondary electrons, which require modelling of the thermal gas density; 2) the thermal pressure, or thermal energy, is used for the relative normalization of the CR energy.

The thermal electron density is modeled as a double β-model [152], following [153]:

$$n_e(r) = n_{0,1} \left(1 + \left(\frac{r}{r_{c,1}}\right)^2\right)^{-3\beta_1/2} + n_{0,2} \left(1 + \left(\frac{r}{r_{c,2}}\right)^2\right)^{-3\beta_2/2} . \quad (2.1)$$

The core parameters are taken from the XMM-Newton measurement [153, 154]. While [153, 154] used the Einstein telescope results from [155] for the peripheral outskirts region parameters, we use instead the more recent Suzaku measurement [147]. We have $(n_{0,1}, r_{c,1}, \beta_1, n_{0,2}, r_{c,2}, \beta_2) = (4.6 \times 10^{-2} \text{ cm}^{-3}, 57 \text{ kpc}, 1.2, 3.6 \times 10^{-3} \text{ cm}^{-3}, 278 \text{ kpc}, 0.71)$, when accounting for the different cosmological models.

We rely on X-ray spectroscopic measurements to describe the gas temperature, as

$$k_B T(r) = 7 \times \left(1 + \left(\frac{r}{r_{t,1}}\right)^3\right) \left(2.3 + \left(\frac{r}{r_{t,1}}\right)^3\right)^{-1} \left(1 + \left(\frac{r}{r_{t,2}}\right)^{1.7}\right)^{-1} \text{keV}, \quad (2.2)$$

<footnote n="12">¹²The software is available at https://github.com/remi-adam/minot.</footnote>
with $k_B$ the Boltzmann constant. The first terms allow us to describe the core temperature as provided by [153], where $r_{t,1} = 73.8$ kpc, given our cosmological model. We introduce the last term in order to account for the temperature drop in the outskirt as measured by [147], where we set $r_{t,2} = 1600$ kpc.

Given the thermal electron density and temperature, we compute the electron thermal pressure as

$$P_e(r) = n_e(r)k_B T(r). \qquad (2.3)$$

Additionally, we assume an ICM helium mass fraction of 0.2735 and we model the metal abundances as constant, with $\frac{Z}{Z_\odot} = 0.3$ [156], using solar abundances from [157]. The electron number density and thermal pressure, together with the abundances are used to compute the full thermodynamic properties of the thermal component following [148].

In figure 1, we present the thermal electron number density, the temperature and the electron pressure radial profiles. We compare our model to other parameterizations available in the literature. In the core, all density profile models agree since they all rely on XMM-Newton data [153]. In the outskirts, the agreement is good up to $R_{500}$. Beyond, ROSAT [158] and Suzaku [147] agree well but the Einstein telescope model [155] leads to a larger electron number density. The temperature profile is typical of that of a cool-core cluster. Our model matches the [153] model in the core and the [147] model in the outskirt according to their domain of validity. All the pressure profile models agree well in the outskirt except for the model derived from the extrapolation of the [153] temperature in the outskirt, which over-predicts the pressure by a factor of a few. Similarly, in the core, the [147] models do not agree with the more direct measurement based on [153] where they extrapolate the profile with an isothermal component. The observed differences are due to the different scales probed by the respective instruments and the extrapolation of the profiles.

Given the definition of our thermal model and the data used to constrain it, we are able to accurately describe the cluster from its core (10 kpc) to the outskirts (2 Mpc). The uncertainties associated with the thermal model are expected to be negligible compared to the uncertainties associated with the non-thermal component. In appendix A, we also present a validation of our thermal model using the Planck Compton parameter map, showing that our model accurately describes the pressure profile of the cluster. We note that because the temperature model based on [153] is valid only up to about 200 kpc, its use leads to an overestimation of the thermal energy by a factor of a few, depending on the details of the line-of-sight integration of the model, when extrapolated beyond its validity region. This is what is done in [88, 89, 91] up to the virial radius ($\sim R_{200}$), which should affect their constraints on the CR by a similar amount.

### 2.1.2 Magnetic field strength

The modelling of the magnetic field strength is necessary when considering jointly the γ-ray emission and the radio synchrotron emission. However, the magnetic field strength and the structure in galaxy clusters remain poorly known to date [159]. We thus consider several approaches in order to model the magnetic field distribution in the Perseus cluster, which will allow us to quantify the associated systematic effect.
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Perseus cluster thermal gas and magnetic field models. Top left panel: thermal electron number density profile. The model by [153] combines XMM-Newton observations in the core and Einstein observations in the outskirt [155]. The red curve is similar, but we have replaced the outskirt model by that obtained in [158] with ROSAT. Top right panel: gas temperature profile. The model by [153] is constrained by the data up to about 200 kpc. The model by [147] is constrained by the data down to about 200 kpc. Bottom left panel: thermal electron pressure profile. In the case of [147], both the combination of gas temperature and density, and the Planck universal pressure profile fit to their data are shown. Bottom right panel: magnetic field strength models according to scaling from literature measurements. The truncation radius is visible at $3R_{500}$.</figcaption></figure>

Our first approach relies on the scaling of the magnetic field to the thermal gas density. In this case, the magnetic field strength is given by

$$ \langle B \rangle (r) = \langle B_{\text{ref}} \rangle \left( \frac{n_e(r)}{n_{e,\text{ref}}} \right)^{\eta_B}, \qquad (2.4) $$

where $\langle B_{\text{ref}} \rangle$ and $n_{e,\text{ref}}$ are magnetic field and density normalization parameters, respectively. We first compute the magnetic field strength using the rotation measure estimate from [160], which gives $\langle B \rangle (10 \text{ kpc}) \sim 25 \mu\text{G}$, and set $\eta_B = 2/3$ assuming magnetic field flux conservation.

The Coma cluster is one of the only clusters for which the magnetic field strength profile was measured [161]. As a second approach, we thus assume that the Coma and Perseus clusters have the same magnetic field strength to gas density ratio. In this case, we
use the values $\eta_B = [0.4, 0.5, 2/3, 0.9]$, corresponding to $\langle B_{\text{ref}} \rangle = [3.9, 4.7, 5.0, 5.4]$ $\mu$G with $n_{\text{e,ref}} = 3.42 \times 10^{-3}$ $\text{cm}^{-3}$, allowing us to account for the uncertainties in the measurement [161].

The structure of the Kelvin-Helmholtz instability visible in Perseus was used to infer the thermal to magnetic pressure ratio, $\beta_{\text{pl}} = P_{\text{gas}}/P_B \sim 200$, for the overall cluster volume prior to sloshing [126]. In a last approach, we use this measurement to infer the magnetic field strength profile, as

$$ \langle B \rangle (r) = \left( \frac{2\mu_0}{\beta_{\text{pl}}} P_{\text{gas}}(r) \right)^{1/2}, \qquad (2.5) $$

where $\mu_0$ is vacuum permeability.

In figure 1 (bottom right panel), we present our models of the magnetic field strength profiles of the Perseus cluster. All models provide a similar order of magnitude, but the differences reflect the difficulty in having an accurate description. The scatter between models nearly reaches an order of magnitude in the core and in the outskirt, with the best agreement between all models around 200–600 kpc. We note that the model with $\eta_B = 0.9$ and $\langle B_0 \rangle = 5.4 \mu$G implies a core magnetic field that reaches about 50 $\mu$G, corresponding to an energy density that is one third of the thermal pressure ($\beta_{\text{pl}} \sim 3$) and thus very high compared to physical expectations. Assuming the pure hadronic model, a higher magnetic field will lead to a lower $\gamma$-ray flux for a fixed radio emission. Thus, the highest magnetic field model gives a conservative estimate of the expected $\gamma$-ray emission. In the following, as a baseline, we use the model based on [160] with $\eta_B = 2/3$ that corresponds to an intermediate estimate.

### 2.1.3 Cosmic-ray model

The $\gamma$-ray emission induced by hadronic interactions is directly related to the spatial and spectral distribution of CRp in the ICM. They are modeled according to a radial profile and a canonical power-law in momentum space,

$$ \frac{dN_{\text{CRP}}}{dpdV}(E, r) = A_{\text{CRP}} p^{-\alpha_{\text{CRP}}} n_e(r)^{\eta_{\text{CRP}}}. \qquad (2.6) $$

The CR radial profile assumes a scaling with respect to the gas density so that only one parameter, $\eta_{\text{CRP}}$, is necessary to describe the spatial distribution. Physically, this parameter allows us to account for the CR dynamics and the competition between advection and streaming, which remains poorly known. We also consider the case where the CRp profile is scaled with respect to the thermal pressure profile, and thus related to the thermal energy. The value of the power-law slope $\alpha_{\text{CRP}}$ is related to the acceleration of CR such as the associated Mach number distribution [162]. The normalization $A_{\text{CRP}}$ is computed given the value of the CRp to thermal energy ratio, $X_{\text{CRP}}(R) = \frac{U_{\text{CRP}}(R)}{U_{\text{th}}(R)}$, by integrating equation (2.6) accordingly. The thermal energy is computed by integrating the pressure profile. The CR energy is computed by integrating the CR distribution over the volume and between $E_{p,\text{min}} = 1.22$ GeV, i.e., the energy threshold of proton-proton interactions, and $E_{p,\text{max}} = 10$ PeV, above which some CRp could escape the cluster [43], although the exact value does not affect our results. In the following, we refer to $X_{\text{CRP}}(R_{500})$ as $X_{500}$ and use this reference for normalization. The value of this parameter has been predicted using numerical simulations, being $X_{500} \sim 1\%$ [43].
The CRe are modeled accounting for two contributions: the primary electrons that are independent from the CRp, and the secondary electrons that are produced from hadronic interactions assuming stationarity. The primary electrons, whenever considered, are modeled similarly to the CRp (equation (2.6)), but with different spectral model to account for energy losses (e.g., power-law or exponential cutoff power-law). We will show in section 2.3 that in practice, primary CRe are irrelevant for our purpose. Secondary CRe are computed as detailed in [148].

## 2.2 Non-thermal radiative processes

Given the ICM model, we compute the different observables associated with the γ-ray emission from hadronic interactions, the radio synchrotron, and inverse-Compton emission. The calculations are done using the MINOT software [148].

### 2.2.1 Hadronic interactions and γ-ray emission

The hadronic production rate of γ-rays is computed by integrating the collision rate of proton-proton interactions multiplied by the energy distribution of γ-rays produced per collision, over the energy of the CR. The computation is based on the parameterization from [163], and its implementation follows the work by [164]. In this paper, we use the Pythia8 proton-proton interaction model and include corrections for proton-nuclei collision [148].

Once the rest frame production rate of γ-rays is computed, the radial profile and energy spectrum of γ-rays, as would be observed from Earth, are obtained by line of sight integration and eventually integrating over the energy or the solid angle, respectively. We also account for the Universe opacity as a function of photon energy using the extra-galactic background light (EBL) model from [165], but this choice does not affect our results given the fact that the corresponding cutoff is at sufficiently high energy, of about 30 TeV [27, 28].

In proton-proton interactions, in addition to γ-rays, secondary electrons are obtained following the same scheme. In this case, however, the prescription from [166] is used instead of that from [163]. The spatial and spectral distribution of secondary CRe is then obtained from the particle injection rate by applying energy losses assuming a steady state scenario.

Modelling uncertainties associated with the particle interaction rates are expected to be accurate to about 30% at CTA energies, for the same parent CRp population (see figure 7 from [148]). This is estimated by comparing the output from the different parametrizations implemented in MINOT (from [163]: Pythia8, Geant4, QGSJET, SIBYLL, and [166]). See also [167] for another recent determination of the γ-ray production cross-section.

### 2.2.2 Inverse-Compton γ-ray emission

Leptonic γ-ray emission arises via the inverse-Compton scattering of CRe (secondaries, but also primary CRe whenever considered) onto cosmic microwave background (CMB) photons. The inverse-Compton MINOT calculations are based on the analytical approximation given by [168]. We will see later that, in practice, inverse-Compton emission should be negligible (from both primary or secondary electrons) for CTA observations and it will be ignored in the following sections.
<table><thead><tr><th rowspan="2">Model</th><th rowspan="2">X<sub>500</sub> (%)</th><th rowspan="2">α<sub>CRp</sub></th><th rowspan="2">η<sub>CRp</sub></th><th colspan="2">F<sup>(had)</sup><sub>500,E<sub>γ</sub>>150 GeV</sub></th></tr><tr><th>F<sup>(IC)</sup><sub>500,E<sub>γ</sub>>150 GeV</sub></th><th>(10<sup>-14</sup> cm<sup>-2</sup> s<sup>-1</sup>)</th></tr></thead><tbody><tr><td>Baseline</td><td>1.0, [0.0, 20.0]</td><td>2.30, [2.00, 3.00]</td><td>1.00, [0.00, 1.50]</td><td>70.2, [0, 11373.8]</td><td>2.1, [0, 625.4]</td></tr><tr><td>Pure hadronic + Taylor [160]</td><td>4.5, [3.2, 14.6]</td><td>2.36, [2.25, 2.65]</td><td>0.78, [0.67, 0.82]</td><td>159.8, [69.1, 234.6]</td><td>4.3, [0.6, 10.7]</td></tr><tr><td>Pure hadronic + Bonafede upper [161]</td><td>3.8, [2.6, 13.4]</td><td>2.37, [2.24, 2.66]</td><td>0.79, [0.65, 0.85]</td><td>132.2, [65.2, 217.4]</td><td>3.3, [0.6, 10.9]</td></tr><tr><td>Pure hadronic + Walker [126]</td><td>7.5, [4.9, 15.8]</td><td>2.35, [2.19, 2.54]</td><td>0.76, [0.69, 0.80]</td><td>273.0, [157.8, 421.0]</td><td>9.2, [2.4, 28.1]</td></tr><tr><td>Pure leptonic</td><td></td><td>-</td><td></td><td>0</td><td>&lt; 28.1</td></tr></tbody></table>

**Table 1.** Summary of the parameter values and their explored range, and the γ-ray flux at CTA energies for the hadronic and inverse-Compton emission (given as: reference value, [min, max]). The flux $F_{500}$ is computed within $\theta_{500}$ by cylindrical integration for energies above 150 GeV and given in units of $10^{-14} \text{ cm}^{-2} \text{ s}^{-1}$. In the case of pure hadronic and pure leptonic models, the central value corresponds to the best-fit model and the interval corresponds to the 68% confidence level. The changes in $\alpha_{CRp}$ for the different magnetic field models is only due to degeneracies in the parameter space and correlations with $\eta_B$. It is not a physical effect and the difference remains well below statistical uncertainties.

### 2.2.3 Synchrotron radio emission

Synchrotron emission should be modeled when considering radio data in order to calibrate our model or check that it does not imply excess radio signal compared to observations. The MINOT software computes synchrotron emission following [169], for which uncertainties are negligible for this work. This assumes that the orientation of the magnetic field is randomized as it is expected for radio halos.

## 2.3 Calibration of the model parameters

In this section, we discuss the different methodologies employed to set our model parameters.

### 2.3.1 Baseline model

As a baseline, we consider the results obtained from numerical simulations and cosmic-ray transport description (e.g., [43, 45]) in order to set the value of the free parameters of the model:<fnref n="13" /> $\eta_{CRp}$, $\alpha_{CRp}$ and $X_{500}$ (see eq. (2.6)). See also [44] for an application with Fermi-LAT. According to these works, the radial distribution of CR is expected to roughly scale with the thermal gas density ($\eta_{CRp} = 1$) when advection by the turbulent gas dominates, but may flatten if diffusion and streaming become significant. We also refer to, e.g., [26, 170] for prediction of the expected parameters.

The spectral slope is expected to be $\alpha_{CRp} \sim 2.3$ and a normalization, $X_{500}$, of a few percent is usually expected. However, large uncertainties are associated with these values, such as the transport of CR for which the streaming velocity is poorly known [45], or the details of acceleration mechanism [21]. Our baseline parameter set is defined as $(X_{500}, \eta_{CRp}, \alpha_{CRp}) = (10^{-2}, 1.0, 2.3)$, but we explore a large range of values, as given in table 1.

In figure 2, we show the hadronic γ-ray observables associated with the models that we consider and how they vary as a function of the parameters. We also provide the associated radial profile of the ratio between integrated CRp energy and thermal energy, $X_{CRp}(< R)$.

<footnote n="13">¹³We note that the main purpose of this work is to compute CTA predictions for constraining the actual CRp population, without entering on the details of how this may population arise in terms of injection.</footnote>
<figure><img src="image_5.png" /></figure>

**Figure 2.** Observable properties of the hadronic γ-ray emission for different values of the model parameters. The inverse-Compton emission is also reported as a thin line with similar color and style on the top panels. *Top left panel:* γ-ray emission spectrum in the case of hadronic processes (within $R_{500}$). *Top right panel:* γ-ray emission profile, as a function of the projected radius, in the case of hadronic processes (computed within [0.15, 50] TeV). *Bottom left panel:* fraction of the enclosed γ-ray flux, as a function of the projected radius (computed within [0.15, 50] TeV). *Bottom right panel:* enclosed CRp to thermal energy ratio, as a function of physical radius. Note that the two bottom panels do not depend on $\alpha_{CRp}$. The magnetic field model used when computing the secondary CR electrons that induce the inverse-Compton signal is the one from [160].

All the models are calibrated to have the same normalization $X_{500} = 10^{-2}$. In the CTA energy range, the γ-ray flux decreases when increasing the slope $\alpha_{CRp}$ and the profile gets more compact when increasing the scaling $\eta_{CRp}$. We note that for a fixed normalization at $R_{500}$, the flux computed within $R_{500}$ also increases with $\eta_{CRp}$ because of the increased proton-proton collision rate due to the spatial overlap of the CR and the target gas. We note that the inverse-Compton emission arising from secondary CRe in the baseline model is always subdominant compared to the hadronic emission.

### 2.3.2 Pure hadronic scenario

It is also possible to use radio data of the Perseus mini-halo to calibrate our model parameters in the case of a pure hadronic scenario. To do so, we extract the Very Large Array (VLA) surface brightness profile at 1380 MHz from [130] and use the spectrum measured at 327, 609, and 1395 MHz from [171], which were taken from [132]. We use an aperture of 15 arcmin,
as it corresponds to the mini-halo extent given in this work. The profile is affected by the
presence of NGC 1275 in the core. In addition the VLA observation may suffer from the
absence of adequately short baselines that are necessary to properly detect the mini-halo
emission at large scales. For this reason, we exclude the data points below 23 kpc, and those
above 80 kpc in the fitting. We also assume a 10% uncertainty on the measured profile [172].

We fit the radio data using a Markov Chain Monte Carlo (MCMC) technique with
the `emcee` package [173]. We use a Gaussian likelihood function and fit simultaneously the
spectrum and the profile. In addition to our three physical parameters, we consider a nuisance
parameter to account for the inter-calibration uncertainty between the spectrum and the
profile, which were extracted from different instruments, and marginalise over this value using
a Gaussian prior centered on unity and with a standard deviation of 10%. Flat priors are used
for the other parameters: $\eta_{CRp} \in [0, 5]$, $\alpha_{CRp} \in [2, 4]$ and $X_{500} \in [0, 0.2]$. These limits are
defined according to realistic physical expectations of the parameter values. The upper limit
of 0.2 on $X_{500}$ corresponds to the expected hydrostatic mass bias of galaxy clusters, itself
dominated by turbulent motions plus magnetic fields, so that $X_{500} > 0.2$ would be unrealistic
(see [174] for a review about the mass of galaxy clusters). For a given model allowed by
the radio data, we compute the $\gamma$-ray observables associated with hadronic interactions and
inverse-Compton emission. We reproduce the constraint for some of the considered magnetic
field models discussed in section 2.1.2. The magnetic field model based on [160] with $\eta_B$ is
used as a reference and we quantify the implications of this choice in appendix B.

In figure 3, we show the pure hadronic best-fit model, its 68% enclosed confidence interval, 100 model realizations randomly sampled from the MCMC chains, and the considered data points. The pure hadronic model provides a good description of the data that we have considered. Given the reference magnetic field model, it favors parameter values of ($X_{500}$, $\eta_{CRp}$, $\alpha_{CRp}$) ~ (5 × 10$^{-2}$, 0.8, 2.5), which is compatible with the MAGIC upper limit [91], but we stress that uncertainties are large and that the parameters are degenerate. Accordingly, the $\gamma$-ray observables are affected by large uncertainties, especially on the CR spectral slope, which leads to large uncertainties for the expected flux in the CTA energy range. We also stress that systematic effects may affect the radio data that we used. Moreover, this approach relied on a magnetic field model, which limits our prediction. Accounting for different parameterization of the magnetic field, we have bracketed the associated systematic uncertainty. Overall, our estimate is expected to be accurate within a factor of about two in this scenario (see appendix B).

### 2.3.3 Pure leptonic contribution

Although it is not physically motivated, as it would require a continuous in situ injection of CRe without protons, we consider the pure leptonic scenario as an exercise to address the inverse-Compton contribution. We use the radio data including CRe only and compute the $\gamma$-ray emission due to inverse-Compton. Since high-energy CRe suffer from major energy losses [175], the only way to observe inverse-Compton emission at CTA energies would be a scenario in which fresh CRe are injected in the ICM, as for the hadronic model. We consider that all the radio emission arises from a continuous particle injection with a power-law energy spectrum, and we apply energy losses given the magnetic field strength, the CMB at the
<figure><img src="image_6.png" /></figure>

**Figure 3.** Constraints on the radio and γ-ray observables in the case of the pure hadronic model, assuming the magnetic field model based on [160] and $\eta_B = 2/3$. The best-fit models are shown as solid blue lines (pink, in the case of inverse-Compton emission). The dashed lines provide the enclosed 68% confidence region and 100 model realizations are shown as light lines. We also report the baseline model for the γ-ray prediction. Top left panel: radio spectrum within a 15 arcmin aperture diameter. The data were taken from [171]. Top right panel: radio profile at 1380 MHz. The data were taken from [130]. Following [172], we use 10% uncertainties and reject the points below 23 kpc (contamination from NGC 1275) and above 80 kpc (possible issues with large scale flux recovery). The grey areas correspond to the discarded data. Bottom left panel: γ-ray spectrum prediction, computed by cylindrical integration within $\theta_{500}$. Bottom right panel: γ-ray projected profile prediction, computed in the range 150 GeV–50 TeV.

cluster’s redshift and the thermal gas. Given these assumptions, the γ-ray emission at CTA energies is similar to the inverse-Compton generated in the pure hadronic model, and thus much smaller than typical expectations from hadronic interactions in numerical simulations or in the pure hadronic model. In practice, continuous injection is not guaranteed and the amount of inverse-Compton emission may thus drop dramatically in the CTA energy range. This is why our computation only provides an upper limit to the inverse-Compton signal.

### 2.3.4 Turbulent reacceleration model

The most favored scenario for the acceleration of relativistic electrons generating large-scale emission in the giant radio halos in galaxy clusters is based on turbulent reaccelera-
tion [30, 31, 33, 36, 176–179]. Turbulent reacceleration has been proposed also for mini-halos as an alternative to pure hadronic models [40, 171]. Turbulence may reaccelerate both primary (CRp and CRe) and secondary particles if they are present in the ICM. In these models, for a given radio luminosity, the level of γ-rays depends on the ratio between primary and secondary electrons that are reaccelerated. Even in the case where only primary CRp and their secondary products are considered, the γ-ray emission that is generated in these models is smaller than the one in pure hadronic models [93, 100]. As a consequence, if turbulence plays an important role in the acceleration process in the mini-halo volume, the γ-ray spectrum calculated in section 2.3.1 provides an optimistic view of the expected signal. Modeling the expected γ-ray signal in the CTA energy band from the Perseus cluster in the case that the radio mini-halo originates from turbulent reacceleration is beyond the aim of this paper. Yet as a reference we just mention that modellings in the case of the Coma cluster, which hosts a giant radio halo, predict a γ-ray signal 4-6 times smaller than that in the pure hadronic scenario [33, 93]. Note that rescaling the CRp normalization of our pure hadronic model (X$_{500}$ = 4.5 × 10$^{-2}$) by the same amount would give a value that is in line with that of our baseline model (X$_{500}$ = 10$^{-2}$), so that our baseline model is in qualitative agreement with what one can expect from turbulence reacceleration.

In the next sections, we will use the models discussed in section 2.3 to address the CTA sensitivity to the CR physics of the Perseus cluster in the context of the proposed KSP observations.

### 3 Modelling the γ-ray emission associated with dark matter

In this section, we build a DM profile for the Perseus cluster following the most up-to-date results both from observations and numerical cosmological simulations.

Assuming DM is completely composed of WIMPs [58, 180], we can compute the expected DM-induced γ-ray emission from any astrophysical source as

$$ \frac{d\Phi_{\gamma}}{dE}(\Delta\Omega, \text{l.o.s}, E) = \frac{d\phi_{\gamma}}{dE}(E) \times J(\Delta\Omega, \text{l.o.s}). \quad (3.1) $$

For a given energy E, line of sight l.o.s, and solid angle ΔΩ subtended by the region of interest, $\frac{d\Phi_{\gamma}}{dE}$ is the DM-induced γ-ray flux, $\frac{d\phi_{\gamma}}{dE}$ contains the spectral information about the expected emission, and J, referred as the “astrophysical factor”, encloses the details of the spatial morphology of the putative signal. In order to perform this DM flux computation, we will assume that WIMPs annihilate into Standard Model (SM) particles in the halo of Perseus (in the following, the annihilation scenario) or that they decay into SM particles (in the following, the decay scenario). According to these two scenarios, the expression in equation (3.1) becomes these two, respectively:

$$ \begin{aligned} \frac{d\Phi_{\gamma}^{\text{ann}}}{dE}(\Delta\Omega, \text{l.o.s}, E) &= \left. \frac{\langle\sigma v\rangle}{8\pi m_{\chi}^2} \frac{dN_{\gamma}}{dE_0} \right|_{E_0=(1+z)E} \times J_{\text{ann}}(\Delta\Omega, \text{l.o.s}); \\ \frac{d\Phi_{\gamma}^{\text{dec}}}{dE}(\Delta\Omega, \text{l.o.s}, E) &= \left. \frac{1}{4\pi m_{\chi}\tau_{\chi}} \frac{dN_{\gamma}}{dE_0} \right|_{E_0=(1+z)E} \times J_{\text{dec}}(\Delta\Omega, \text{l.o.s}), \end{aligned} \quad (3.2) $$
where $\frac{dN_γ}{dE}$ is the number of photons per unit energy we expect from a given annihilation channel and DM particle mass $m_χ$, $<σv>$ is the velocity-averaged annihilation cross-section in the annihilation scenario and $τ_χ$ is the mean lifetime of the DM particle in the decay scenario. In deriving these expressions, we have assumed Majorana WIMPs (for Dirac WIMPs it would be necessary to multiply by a factor 1/2 [180]). To obtain the corresponding fluxes, we use for $\frac{dN_γ}{dE}$ the results from [181], including electro-weak corrections.<fnref n="14" /> For the thermally-averaged cross-section, we expect a value around $3 \times 10^{-26} \text{ cm}^3 \text{ s}^{-1}$ [183, 184] to produce the observed DM relic abundance [185], while the main bound for the decay lifetime to account for the observed DM density is the current age of the Universe $\sim 10^{17}$s, despite the efforts of trying to obtain tighter constraints from just gravitational probes [186]. We recall that for decaying DM only half of the energy budget stored in $m_χ$ is available per SM particle produced in two-body decays. We can also define the spatial morphology of the signal, the so-called J-factor, as:

$$ 
\begin{align} 
J_{\text{ann}}(\Delta\Omega, \text{l.o.s}) &= \int_{\Delta\Omega} d\Omega \int_{l.o.s} \rho_{\text{DM}}^2(r, l)dl; \\ 
J_{\text{dec}}(\Delta\Omega, \text{l.o.s}) &= \int_{\Delta\Omega} d\Omega \int_{l.o.s} \rho_{\text{DM}}(r, l)dl, 
\end{align} 
\tag{3.3} $$

where $\Delta\Omega = 2\pi(1 - \cos\alpha_{\text{int}})$, being $\alpha_{\text{int}}$ the integration angle and $\rho_{\text{DM}}(r)$ the DM density profile.

From equation (3.3) the main dependencies in both scenarios can be deduced [73]:

$$ J_{\text{ann}} \propto \frac{M_{200}c_{200}^3}{d_L^2}; \quad J_{\text{dec}} \propto \frac{M_{200}}{d_L^2}, \tag{3.4} $$

which are the mass $M_{200}$,<fnref n="15" /> the luminosity distance to the Earth $d_L$ and the concentration $c_{200}$,<fnref n="16" /> for the annihilation case. From these dependencies, one can clearly see why local galaxy clusters are the best targets to consider to test decaying DM [187], since they are the most massive objects in the Universe. For the case of annihilation it may not seem so straightforward. Yet, in the last years some studies where performed [71] comparing the suitability of galaxy clusters and dwarf spheroidal galaxies (dSphs) for γ-ray DM searches, concluding that galaxy clusters can also provide competitive results. One key element for this outcome is to consider the distribution of smaller halos that galaxy clusters should host according to the ΛCDM structure formation theory [77], usually called subhalos. The role of these substructures come into play through the $c_{200}$ dependency, since we expect subhalos of a given mass to be much more concentrated than main halos of the same mass. In contrast, in the decay scenario (see equation (3.4)) the subhalos do not provide a sizable contribution since their masses are really low compared to the distance. Thus, we will not consider them for this scenario.

<footnote n="14">¹⁴Ref. [181] only includes prompt emission of γ-rays, not secondary radiation as the one from inverse-Compton processes. However, electromagnetic cascades can be relevant for hard cosmic-ray spectra or some DM annihilation channel [182].</footnote>

<footnote n="15">¹⁵Mass enclosed in a radius where an spherical overdensity is 200 times the critical density of the Universe $\rho_{\text{crit}}$.</footnote>

<footnote n="16">¹⁶Dimensionless quantity that describes how much concentrated halos are and it is dependent on the mass of the object.</footnote>
The morphology of the expected DM γ-ray emission is strictly determined by the DM density profile $\rho_{DM}$ and in case of annihilation, particularly by the DM substructure. In the next sections, we present the models we build for the DM distribution in the main halo of Perseus as well as a model for the population of subhalos.

## 3.1 Perseus main halo

The main DM halo of Perseus can be described through a Navarro-Frenk-White (NFW) density profile [188, 189]:

$$ \rho_{\text{NFW}}(r) = \frac{\rho_0}{\left(\frac{r}{r_s}\right)\left(1+\frac{r}{r_s}\right)^2}, \qquad (3.5) $$

where $r_s$ is the scale radius and $\rho_0$ the normalization of the DM density. Introduced as the result of DM-only N-body simulations, the NFW model falls in the family of “cuspy”-like profiles.

To obtain the parameters of the NFW profile, we need to consider a concentration-mass ($c-M$) relation. Consistently, this relation should comprehend the mass scales involved in cluster physics. Taking this into account, we use the parametrization developed in [73] for main halos. This concentration-mass relation was found to have an associated $1\sigma$ scatter of 0.14 dex. We build the DM density profile starting from the measured mass of Perseus. We consider the halo to fulfill the spherical collapse model for overdensities $\Delta = 200$ times the critical density of the Universe $\rho_{crit}$. This allows us to calculate $M_{200}$ shown in table 2 from the measured X-ray mass, $M_{500}$, provided in [147].

We can also obtain the radius $R_{200}$ that contains an enclosed mass $M_{200}$, as

$$ R_{200} = \left( \frac{M_{200}}{\frac{4}{3}\pi\Delta\rho_{\text{crit}}} \right)^{1/3}, \qquad (3.6) $$

whose corresponding projected angle is $\theta_{200} = \arctan\left(\frac{R_{200}}{d_L}\right)$. Now we can easily compute the scale radius in equation (3.5) as

$$ r_s = \frac{R_{200}}{c_{200}}. \qquad (3.7) $$

Finally, the normalization of the density profile $\rho_0$ is obtained by imposing the recovery of $M_{200}$ after the integration over the cluster volume of the profile ($M_{200} = \int_0^{R_{200}} \int_0^{4\pi} \rho_{\text{NFW}}(r)r^2 dr d\Omega$) and isolating $\rho_0$:

$$ \rho_0 = \frac{2 \Delta_{200} \rho_{\text{crit}} c_{200}}{3 f(c_{200})}, \qquad (3.8) $$

where $f(c_{200}) = \frac{2}{c_{200}^2} (\ln(1+c_{200}) - \frac{c_{200}}{1+c_{200}})$. Table 2 lists the corresponding parameters of the obtained NFW DM density profile following the above description.

It is interesting to compare the value we obtained for $c_{200}$ and reported in table 2 with that obtained from observational data. Indeed, from X-ray observations an observational value of $c_{200} = 5.0 \pm 0.5$ was inferred [190], which well matches our adopted $c_{200}$.
<table><caption>**Table 2.** NFW density profile parameters for the Perseus galaxy cluster; see equations (3.5)–(3.8) for the definition of each parameter and accompanying text for details on their derivation.</caption><thead><tr><th>M<sub>200</sub><br>[10<sup>14</sup> M<sub>☉</sub>]</th><th>c<sub>200</sub></th><th>r<sub>s</sub><br>[kpc]</th><th>log<sub>10</sub> ρ<sub>0</sub><br>[M<sub>☉</sub>/kpc<sup>3</sup>]</th><th>R<sub>200</sub><br>[kpc]</th><th>θ<sub>200</sub><br>[deg]</th></tr></thead><tbody><tr><td>7.5</td><td>5.0</td><td>370.8</td><td>6.1</td><td>1865.0</td><td>1.4</td></tr></tbody></table>

## 3.2 The role of substructures

According to the ΛCDM hierarchical bottom-up structure formation scenario [191, 192], galaxy clusters should host a significant amount of substructure or subhalos. Due to stripping processes, these subhalos are expected to have higher concentrations than main field halos of the same mass (e.g., [74]). Because of these high concentrations, we expect them to have a considerable impact in the annihilation fluxes (equation (3.4)). This enhancement effect in $J_{\text{ann}}$-factor, usually called *boost factor B*, was already estimated in [71] for a selected list of local clusters. These authors found that halo substructure in galaxy clusters could provide boost factors of the order $B \sim 40$, in general agreement with current, more refined boost calculations from N-body simulations [74, 193]. In order to account for the contribution of these substructures to the $J_{\text{ann}}$-factor, we factorize the computation in the following way:

$$ 
\begin{align} 
J_{\text{ann}}(\Delta\Omega, \text{l.o.s}) &= \int_{0}^{\Delta\Omega} d\Omega \int_{\text{l.o.s}} \left( \rho_{\text{main}}(r) + \sum_{i}^{N_{\text{sub}}} \rho_{\text{sub}}(r) \right)^2 dl = \tag{3.9} \\ 
&= J_{\text{main}}^{\text{ann}} + \langle J_{\text{sub}}^{\text{ann}} \rangle + \langle J_{\text{cross-prod}}^{\text{ann}} \rangle, \nonumber 
\end{align} 
$$

where $N_{\text{sub}}$ is the total number of subhalos, $\rho_{\text{main}}$ is the DM profile of Perseus' main halo obtained in section 3.1 and $\rho_{\text{sub}}$ the DM profile for each subhalo. Due to the complex dynamics that substructures undergo (e.g., tidal stripping, dynamical friction, interaction with baryons, etc. — see [192, 194–196]), the precise survival probability of the smallest subhalos is not yet known [197–200]. In any case, when it comes to the subhalo density profile $\rho_{\text{sub}}$, deviations from NFW are expected to be pronounced mainly in the outermost subhalo regions, where mass losses are severe. Therefore we keep adopting the NFW profile also to model the subhalo inner structure. We remark that due to their small extension, we do not expect individual subhalos to be spatially resolved by CTA. Therefore, we decided to average their contribution in $\langle J_{\text{sub}}^{\text{ann}} \rangle$. We then model the population of subhalos as:

$$ \frac{d^3 N}{dV dM dc} = N_{\text{sub}} \frac{dP_V}{dV} (R) \frac{dP_M}{dM} (M) \frac{dP_c}{dc} (M, c), \qquad (3.10) $$

where $P_i$ with $i = V, M, c$ represents the probability distribution in the volume of the main halo $V$, of the subhalo masses $M$ and the subhalo concentrations $c$. The probability distributions in equation (3.10) have been modeled based on results of DM-only cosmological simulations, namely from Via Lactea-II (VL-II, [201, 202]), Aquarius [203, 204] and the work of [74]. We describe each of these distributions in the following:

* Subhalo distribution within the main halo $\frac{dP_V}{dV}(r)$: As we assume that Perseus main halo is spherically symmetric, the distribution within the volume collapses to the distribution
in radius, so in the following we refer to it simply as the Subhalo Radial Distribution (SRD). Based on N-body simulations, we choose to follow the *antibiasing relation* [202],<fnref n="17" /> defined as:

$$ \frac{dP_V}{dV}(r) = \rho_{\text{main}}(r) \frac{r/r_a}{1 + \frac{r}{r_a}}, \qquad (3.11) $$

where $r_a$ works as a scale radius and it is defined by the fraction of the total mass that will be in form of subhalos $f_{\text{sub}}$. For the distribution of subhalos, we also consider the number of substructure levels to $N_{\text{lvl}} = 2$ (subhalos inside subhalos).<fnref n="18" />

* Subhalo mass distribution $\frac{dP_M}{dM}(M)$: Also known as the subhalo mass function (SHMF), is usually modelled using a power-law (e.g., [207]):

$$ \frac{dP_M}{dM} \propto M^{-\alpha}. \qquad (3.12) $$

Values range between $\alpha = 1.9$ [203, 204] and $\alpha = 2.0$ [201]. This slope is key to evaluate the contribution of substructure to the $J_{\text{ann}}$-factor [74]. For higher values of $\alpha$, i.e., a more numerous population of small subhalos, we obtain proportionally higher boosts. As a way to account for this uncertainty, we will define different benchmark models covering different physical scenarios for the abundance of subhalos. A value of $\alpha = 1.9$, more in line with latest simulations results [192, 208], will lead to conservative values for the $J_{\text{ann}}$-factor, while $\alpha = 2.0$ will yield an upper bound to subhalo boost values. The mass budget in both cases is a fraction of the total mass, $f_{\text{sub}}$. Assuming a value for the minimum subhalo mass of $M_{\text{min}} = 10^{-6} M_{\odot}$ and that the maximum subhalo mass in terms of the host is $M_{\text{max}} = 0.01$,<fnref n="19" /> we need different values of $f_{\text{sub}}$ to conserve the total mass. Integrating the SHMF with the different values of $\alpha$ for the selected $M_{\text{min}}$ and $M_{\text{max}}$, we obtain a value of $f_{\text{sub}} = 0.182$ for $\alpha = 1.9$, and $f_{\text{sub}} = 0.319$ for $\alpha = 2.0$.

* Subhalo concentrations $\frac{dP_c}{dc}(c, M)$: In the same way we assumed a concentration-mass relation for the main halo, we select a $(c - M)$ relation to describe subhalo DM profiles. We adopt the state-of-the-art $(c - M)$ subhalo model by [74], which includes a radial dependence of the concentration within the main halo:

$$ c_{200}(M_{200}, x_{\text{sub}}) = c_0 \left[ 1 + \sum_{i=0}^{3} \left[ a_i \log_{10} \left( \frac{M_{200}}{10^8 h^{-1} M_{\odot}} \right) \right]^i \right] \times [1 + b \log_{10}(x_{\text{sub}})], \quad (3.13) $$

where $c_0 = 19.9$, $a_i = [-0.195, 0.089, 0.089]$, $b = -0.54$ and $x_{\text{sub}}$ refers to subhalo distance with respect to the center of the host halo. The importance of including a

<footnote n="17">¹⁷There are other options available for the SRD modelling, e.g., the one based on results from the Aquarius simulations, where they fit the SRD to an Einasto profile [203, 204]. Yet, this provides a subhalo distribution that is compatible compared to the VL-II results, which is the one used in this work.</footnote>

<footnote n="18">¹⁸In [73, 205], they concluded that adding more levels of substructure only contribute up to a 5% to the total flux. See, however, the recent [206].</footnote>

<footnote n="19">¹⁹We note that these values are generally accepted and standard in the community as today, representatives and consistent for the defined models. We could have selected a higher (lower) value for $M_{\text{min}}$, which would translate into slightly lower (higher) values of the boost. For instance, for a value of $M_{\text{min}} = 10^{-9}$, the boost only increases less than 5%. The reason for this small variation is the flattening of the concentration-mass relation at lower halo and subhalo masses, see [73]. All in all, this exercise proves that the variation of these parameters ($M_{\text{min}}, M_{\text{max}}$) yield $J_{\text{ann}}$-factors already encapsulated by our benchmark models.</footnote>
<table><caption>**Table 3.** Summary of the defined benchmark DM models for the annihilation interaction scenario. From left to right, $\rho_{\text{sub}}$ is the main DM profile for both the main halo and the subhalos, $(c - M)_{\text{main}}$ is the concentration-mass relation used for Perseus main halo, SRD is the Subhalo Radial Distribution, $(c - M)_{\text{sub}}$ is the concentration-mass relation for the subhalos, $\alpha$ is the slope of the SHMF (equation (3.12)) and $f_{\text{sub}}$ is the fraction of the total mass bound in substructures, see text for more details. Regarding the models, “SC+14” refers to the concentration model of [73], *antibias ed* to the SRD from [202] from equation (3.10), “M+17” to the $(c - M)$ subhalo model of [74] from equation (3.13).</caption><thead><tr><th>Model</th><th>ρ<sub>DM</sub></th><th>(c − M)<sub>main</sub></th><th>SRD</th><th>(c − M)<sub>sub</sub></th><th>α</th><th>f<sub>sub</sub></th></tr></thead><tbody><tr><td>MIN</td><td>NFW</td><td>SC+14</td><td>-</td><td>-</td><td>-</td><td>0</td></tr><tr><td>MED</td><td>NFW</td><td>SC+14</td><td><i>Antibias ed</i></td><td>M+17</td><td>1.9</td><td>0.182</td></tr><tr><td>MAX</td><td>NFW</td><td>SC+14</td><td><i>Antibias ed</i></td><td>M+17</td><td>2.0</td><td>0.319</td></tr></tbody></table>

$(c - M)$ relation specifically derived for subhalos resides, as already stated, in the fact that subhalo concentrations are known to be higher than that of field halos of the same mass [74]. Thus, their contribution to the $J_{\text{ann}}$-factor is expected to be critical. Note that, in deriving this relation, [74] assume NFW profiles for subhalos, keeping our modelling consistent. Every concentration-mass relation is known to exhibit an intrinsic scatter [74, 209], yet it is highly computationally expensive to take it into account for each subhalo in the field halo. Then, for the sake of this study we decide to neglect it, as its impact on the $J_{\text{ann}}$-factor will lie within the spread we have by considering different values for $\alpha$ in the SHMF.

After the detailed discussion on the parametrization followed for the subhalo population modelling, we establish three benchmark models for the computation of the $J$-factors. Each of them represents an expected different level of contribution of the subhalo population to the annihilation flux:

* MIN: considers the Perseus main halo and neglects the existence of substructures.

* MED: the SRD follows the *antibias ed* relation. We adopt $\alpha = 1.9$ for the slope of the SHMF (equation (3.12)), using a coherent subhalo mass fraction of $f_{\text{sub}} = 0.182$.

* MAX: similar to the MED model but we choose $\alpha = 2.0$ for the slope of the SHMF (equation (3.12)), with a coherent adjustment of the subhalo mass fraction to $f_{\text{sub}} = 0.319$.

These models and their values are summarized in table 3. With the definition of the above benchmark models, we aim to bracket a wide range of possible substructure scenarios. This will translate into a bracketing for the possible values of the $J_{\text{ann}}$-factor, being the MED model our best guess and for which we will produce the main results, and the MIN and MAX models as realistic lower and upper limits to the contribution levels of substructures. For decaying DM, we will use a realistic value for $J_{\text{dec}}$ assuming the MIN model, since subhalos do not provide a sizeable contribution to it (see equation (3.4)).
<table><caption>**Table 4.** Perseus integrated J-factors in the case of the three benchmark annihilation scenarios (MIN, MED, MAX — see table 3) and the $J_{\text{dec}}$-factor for the decay scenario. In all cases, we integrate the DM signal up to $R_{200}$ ($\alpha_{\text{int}} = \theta_{200}$). See text for details.</caption><thead><tr><th></th><th>Annihilation</th><th>Decay</th></tr><tr><th></th><th>log<sub>10</sub> J<sub>ann</sub> (GeV<sup>2</sup> cm<sup>-5</sup>)</th><th>log<sub>10</sub> J<sub>dec</sub> (GeV cm<sup>-2</sup>)</th></tr></thead><tbody><tr><th>MIN</th><td>17.42</td><td>19.20</td></tr><tr><th>MED</th><td>18.43</td><td>-</td></tr><tr><th>MAX</th><td>19.20</td><td>-</td></tr></tbody></table>

### 3.3 Dark matter annihilation and decay fluxes

Once defined the DM models we use for Perseus cluster and for its expected subhalo population, we are interested in computing their γ-ray induced fluxes, the J-factors. We use the publicly available CLUMPY code [205, 210, 211] to compute them. The obtained integrated J-factors for the cluster and for each benchmark model are summarized in table 4 and shown in the left panels of figure 4. The right panels in this same figure show the differential J-factors as a function of the angle from the center of Perseus.

From the $J_{\text{ann}}$-factor values given in table 4 for each benchmark model and from figure 4 upper panels, we can better understand the impact of including subhalos in our calculations. Indeed, substructures do not only boost the expected annihilation signal but also modify its spatial morphology, subhalos being particularly important in the outskirts of the cluster. The latter is the main reason to also take into account the expected extension of the cluster DM emission in the analysis itself, as it will be described later on. Although the precise properties of the subhalo population within the cluster is the main origin of uncertainty in the computation of the $J_{\text{ann}}$-factors, we also wanted to check the impact of the uncertainty in the $M_{200}$ value and the one from the scatter in the concentration-mass relation. The mass uncertainty from the study of [147] is too small to have an impact itself, thus we compute two different values of $M_{200}$ assuming extreme values for the hydrostatic mass bias ($b_{\text{HSE}} = 0, 0.3$) and propagate this to the $J_{\text{ann}}$-factor computation. This results into $\sigma_{J_{\text{ann}}} = 0.002-0.005$ dex correspondingly, which is two orders of magnitude smaller than the variation introduced from the substructure benchmark models ($\sigma_{J_{\text{ann}}} \sim 1$ dex). If we also include the 0.14 dex from the scatter of the concentration-mass relation for each of the above extreme mass values, we obtain $\sigma_{J_{\text{ann}}} = 0.2$ dex, a $J_{\text{ann}}$-factor uncertainty similar to what is typically obtained for dSphs [101]. Note though that this $\sigma_{J_{\text{ann}}}$ is still well within the uncertainty related to the different considered models for the substructure population. Therefore, from now on we only consider as a theoretical uncertainty the spread among the three substructure benchmark models. As for decay, we compute the uncertainty of $J_{\text{dec}}$ from only the use of different mass values and from the scatter in the concentration-mass relation as well, since the substructures do not play a role for this case. We found $\sigma_{J_{\text{dec}}} \lesssim 10^{-3}$ dex for all cases, thus we decided to disregard this in the following.

To quantify the enhancement in the total annihilation flux due to the presence of subhalos, we define a subhalo *boost factor* as:

$$B = J_{\text{ann}}^{X} / J_{\text{ann}}^{\text{MIN}} - 1, \qquad (3.14)$$
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Integrated J-factors (equation (3.3)) versus the integration angle $\alpha_{\text{int}}$ (left panels) and differential J-factors versus the radial angle $\theta$ (right panels) for the annihilation scenario (top panels) and decay scenario (bottom panels) for the three substructure benchmark models correspondingly. In all the panels $\theta_{200}$ and the projected angle of the scale radius, $\theta_s$, are shown as a measure of the containment of the emission and the extension of the cluster.</figcaption></figure>

where $X$ can be either MED or MAX ($B = 0$ means no substructure contribution in our definition).<fnref n="20" /> For Perseus and our MED and MAX models, we obtain boost factors of $B = 9.2$ and $B = 59.3$, respectively.

We can compare our obtained J-factors for the annihilation scenario with the results in [71].<fnref n="21" /> For the MIN scenario (no substructure inclusion), these authors obtain $\log_{10}(J_{\text{ann}}) = 17.35$ (units of GeV² cm⁻⁵), very similar to our results (table 4). For the case where they include substructure (which, according to their modelling, would be roughly similar to our MED case), the $J_{\text{ann}}$-factors increase up to $\log_{10}(J_{\text{ann}}) = 18.23$, also compatible with our MED results (table 4). The *boost factor* that they obtain between their two models is $B = 6.6$ (following our definition in equation (3.14)), a value just around one point lower than our results for the MED model (table 4). This slight mismatch is expected due to some differences in the subhalo modelling in each case. In order to perform a comparison with more recent

<footnote n="20">²⁰We recall that other works may use a different definition for the boost, where $B = 1$ means no substructure contribution.</footnote>

<footnote n="21">²¹We caution that there exists a factor $\frac{1}{4\pi}$ difference with their J-factor definition.</footnote>
results, we can examine the boost values we should expect from using the $(c-M)$ relation of [74]. In this case, for a Perseus-like galaxy cluster halo, they obtain $B \sim 9$ for $\alpha = 1.9$ and $B \sim 72$ for $\alpha = 2.0$. The agreement with our MED scenario is absolute, while for MAX we get a relative discrepancy of $\sim 20\%$.<fnref n="22" />

A similar comparison can be performed for the decay scenario as well. Given their large masses and relatively close distances, galaxy clusters are the perfect laboratories for probing DM decay and, accordingly, there have been several previous studies in this matter (see refs. in section 1). One of the latest works on DM decay in the TeV range focused on Perseus and consisted of an observational campaign performed by the MAGIC collaboration [117]. In this work, authors used the main halo model introduced in [71] to compute the expected decay flux, obtaining a value of $\log_{10}(J_{\text{dec}}) = 19.18$ (units of GeV cm$^{-2}$), compatible within 5% with our value from table 4.

As a last step, we create 2D spatial templates of the expected DM annihilation/decay emissions using CLUMPY. Four maps are created, three for annihilation (each one corresponding to one of our benchmark models in table 3), and one more for decay. These maps are all shown in figure 5. As seen already in figure 4, the presence of substructures mainly impacts the morphology of the DM signal in the outer regions of the cluster.

# 4 Observation setup

Before considering the CTA data analysis, we discuss the configuration of the observations. In addition to the instrumental background due to CR induced air showers, we present the modelling of the known sources located in the Perseus cluster region which may affect the analysis. We also investigate how the observing strategy will impact the CTA sensitivity to a putative cluster signal.

## 4.1 CTA pointing configuration

The CTA observing setup will consist in a set of pointings, with a radial offset with respect to the Perseus cluster center, $\theta_{\text{pointing}}$. Such offset is necessary to analyse the data using classical ON-OFF techniques used in Imaging Air Cherenkov Telescopes (IACTs), in whose case the background is estimated using a mirror region equivalent to the region of interest assuming the azimuthal symmetry of the background with respect to the pointing center [213]. For template-based analysis, in which the full region of interest is modeled using dedicated templates that describe the different sky components and are fit to the data, such offset could be set to zero. Various position angles of the pointing offset can also be considered. This would allow us to cover a larger field around the cluster and thus to increase the chance of serendipitous discovery of new sources. The CTA field of view depends on the energy, given the field of view of the respective telescope class dominating each energy regime: about 4.3 deg diameter for the large-sized telescope (LST) and 7.5 deg diameter for the medium-sized telescope (MST). See [120] and [121].

<footnote n="22">²²In their calculation the mass in form of subhalos in the cluster is not subtracted from the smooth DM distribution, hence this fact leading to an overestimation of the boost, which becomes more important as more mass is modeled as substructures.</footnote>
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Two-dimensional spatial templates of the expected DM emission in Perseus, quantified in terms of the differential J-factors on the z axis, for the three considered MIN, MED, MAX annihilation scenarios (top panels and left bottom panel) and for decay (bottom right panel). See tables 3 and 4 for details on each DM model.</figcaption></figure>

A total of 300 hours are proposed for observations of the Perseus cluster in the course of the CTA key science program and this is what we assume as our ON-source time for the remaining. Nonetheless, we consider a dead-time fraction of 5% so that the effective observing time is 285 hours. The observing time will be split into the different pointings. The pointing directions are split between 8 positions, equally spaced from 0 to 315 deg around the cluster reference center, with steps of 45 deg. Nevertheless, we stress that our results do not significantly depend on the number of position directions that we consider.

Because the Perseus cluster field is only observable from the CTA north site (La Palma), we consider the instrument response function (IRF, which connects the sky γ-ray signal to the actual measurement) associated with long observations [214] from the north site.<fnref n="23" /> We neglect any variation of the zenith angle of the target source during the observations and assume that the zenith angle is fixed to 20 deg for all pointings. For comparison, the zenith angle ranges from 12 to 36 deg for the MAGIC observations presented in [89], also performed from La Palma and within which range there is no significant difference. The IRFs are also used to model the instrumental background. Note that IRFs describing the initial CTA array configuration are now available: 9 MST and 4 LST instead of the 15 MST and 4 LST used

<footnote n="23">²³The IRFs are available at https://www.cta-observatory.org/science/ctao-performance/ for “Alpha configuration” or prod5 and https://zenodo.org/record/5163273#Y5nebbKZNQ0 for prod3b.</footnote>
<table><tr><td>Observing time</td><td>300 hours</td></tr><tr><td>IRF</td><td>prod3b-v2 <sup>†</sup>, North_z20_50h</td></tr><tr><td>Pointing offset from cluster center</td><td>1.0 deg</td></tr><tr><td>Azimuthal position offset</td><td>[0,45,90,135,180,225,270,315] deg</td></tr></table>

Notes. †This IRF corresponds to the full configuration with 15 MST and 4 LST. More recent IRFs are now also available, corresponding to the CTA initial configuration including 9 MST and 4 LST (prod5 [212]). The corresponding on axis sensitivity is reduced by a factor of about 30% at high energies ($E \gtrsim 500$ GeV), but the off-axis sensitivity is slightly improved.

<caption>**Table 5.** CTA observation configuration summary.</caption>

for the present work. They would imply a reduced sensitivity by about 30% at high energies, compared to the IRFs used here, although the off-axis sensitivity improved and could make up for some of the deficit in telescope number.

We provide a summary of the CTA observation setup in table 5. A schematic view of the CTA observations configuration is provided in figure 6. We include NGC 1275 and IC 310 (see discussions in section 4.2). The X-ray peak (reference center) is aligned with NGC 1275 and IC 310 is offset toward the southwest by about 0.6 deg but remains within $\theta_{500}$. For a given pointing, the OFF regions are spread over a circle around the pointing for which the radius corresponds to the pointing offset from the source. This is illustrated by the yellow (OFF) and orange (ON) filled regions for which we have used a radius of 0.5 deg. Also shown are the locations of those sources in the Fermi-LAT 4FGL catalog [215, 216] within five degrees from the cluster center.

## 4.2 Background sky

In addition to the diffuse emission from the Perseus cluster, it is necessary to model other sources that are located around the target cluster. At very high energies, two $\gamma$-ray point sources have been detected by the MAGIC telescopes, namely NGC 1275 and IC 310.

NGC 1275, the brightest galaxy in the Perseus cluster, is a variable source presenting flare activities up to fifty times its mean flux at $E > 100$ GeV, with observed day-by-day variability [217]. It is not possible to predict what will be the exact state of the source at the time of CTA observations. Yet, we expect that a large majority of the data will be obtained when the source is in a quiescent state and periods with high intensity can be removed from analysis. We therefore assume that the spectrum of NGC 1275 follows the quiescent energy spectrum obtained from MAGIC observations as presented in [91]. This spectrum is well described, as a function of energy $E$, by a simple power-law given by

$$ \frac{dN}{dEdSdt} = 2.1 \times 10^{-11} \left( \frac{E}{200 \text{ GeV}} \right)^{-3.6} \text{cm}^{-2}\text{s}^{-1}\text{TeV}^{-1}. \quad (4.1) $$

IC 310, a member galaxy of the Perseus cluster, is also a variable source presenting high-amplitude and short duration flares [218]. While its spectral shape was observed to not significantly change, its amplitude does by a factor of up to $\sim 7$. Similarly to NGC 1275,
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Schematic view of the Perseus cluster region and the assumed observation setup. An illustration of the ON-source and OFF-source regions is shown for one pointing, assuming a 0.5 deg aperture radius. IC 310 and NGC 1275 are observed at very high energy and are shown as cyan stars. 4FGL sources within a 5 deg radius around the cluster center are also depicted as blue crosses. For comparison, the point spread function (PSF) of CTA at 1 TeV is shown in the bottom left corner.</figcaption></figure>

its spectrum is modeled as a power-law,

$$ \frac{dN}{dEdSdt} = 0.741 \times 10^{-12} \left( \frac{E}{1 \text{ TeV}} \right)^{-1.81} \text{cm}^{-2}\text{s}^{-1}\text{TeV}^{-1}, \quad (4.2) $$

corresponding to the quiescent state of the source [218]. This source state corresponds to the large majority of the reported observations. We note that, because we apply EBL absorption a posteriori when performing simulations, we use the intrinsic MAGIC best-fit model. The equatorial coordinates of IC 310 are set to (R.A., Dec.) = (49.179, +41.325) deg.

In addition to NGC 1275 and IC 310, several sources have been detected at lower energies with Fermi-LAT around the Perseus cluster. In appendix C we list these sources and estimate their fluxes in the CTA energy range. Given their location, none of them should affect the analysis.
Finally, we consider the contribution from the Galactic diffuse emission in the region around the Perseus cluster using the work by [219]. Within $\theta_{500}$, the emission is expected to be smooth and present a soft gradient. Compared to other contaminant sources, the instrumental background, and the cluster diffuse emission (e.g., baseline CR model), the Galactic $\gamma$-ray foreground should be largely subdominant over all the considered energy range. We refer to appendix D for more details. Because of the uncertainties in the foreground model, its spatial structure and the fact that it is subdominant, it is not considered in the following.

## 4.3 Towards an optimal observing strategy

In this section, we explore how the observation setup parameters described in section 4.1, together with the background sky discussed in section 4.2, will affect the expected detection significance. To do so, and as it will be the case in section 5, we use the ctools software$^{24}$ [220]. We focus on the CR component since the observation setup selection will be driven by the CR case. We have also checked the influence of different setups for the DM case (appendix I) and the chosen one is appropriate for both science cases. Our background model includes both the instrumental background and the two AGNs NGC 1275 and IC 310. We consider the case of ON-OFF analysis, and the case of template fitting.

**ON-OFF analysis.** We first consider the case of a classical ON-OFF analysis. A pointing offset, $\theta_{pointing}$, is necessary to define the OFF regions that will be used to monitor the background (see figure 6). Indeed, given the azimuthal symmetry of the instrumental background, the definition of an OFF region, symmetric to the ON region in terms of background, is not possible if the pointing center coincides with the region of interest (ON region). In addition to the pointing offset, the aperture radius of the ON region will affect the detection significance depending on the shape of the signal. For instance, a small aperture will be favored for compact sources because it will lead to a lower instrumental background while keeping most of the signal photon counts. By contrast, a larger aperture will be preferred for more extended sources. The size of the ON region is also limited by the pointing offset, for which we aim at testing the impact. Here, we request at least 3 independent OFF regions so that the radius of the ON region, $\theta_{ON}$, is constrained to $\theta_{ON} \le \sqrt{2}\theta_{pointing}$. Such condition will also affect the maximum significance for a given pointing offset. We also consider the possibility of masking the known point sources. We use an energy-dependent aperture radius proportional to the point spread function (PSF) 68% containment angle $\theta_{PS}(E) = N_{PSF} \times \text{PSF}(E)$, where $N_{PSF}$ can also be varied and for which the optimal value may depend on the diffuse cluster model.

In order to optimise the observation setup, we compute the expected significance according to [221]:

$$ \sigma = \sqrt{2 \times \left( d \times \ln\left(\frac{d}{m}\right) + m - d \right)}. \qquad (4.3) $$

The quantity $d$ corresponds to the expected photon counts (between 150 GeV and 50 TeV) of the diffuse cluster signal plus background in the ON region, after accounting for the

<footnote>$^{24}$http://cta.irap.omp.eu/ctools/.</footnote>
point-source mask. The quantity $m$ corresponds to the expected photon counts associated with the background only for the same region. The background includes both the instrumental background and the point sources, assuming that their spectral energy distribution (SED) is known (see section 4.2) and kept fix. The cluster and point-source models are convolved with the IRFs. We compute $\sigma$ for a set of pointing offset values, for each of which we vary the values of the ON region aperture $\theta_{ON}$ and the size of the point-source mask $\theta_{PS}$. The counts are summed over all the energy bins.

**Template fitting.** Beside the classical ON-OFF analysis technique, we consider the case of template fitting, which might be more appropriate for diffuse sources assuming that the background can be properly modeled [122], as it will be further explored in section 5. To do so we perform simulations of the observations by varying the pointing offset, and use the maximum likelihood technique to fit the data and recover the test statistic value. We fit for the normalization of the cluster model, for the instrumental background normalization and spectral tilt, and for the normalization and spectral index of the two point sources. We use the square root of the recovered test statistic ($\sqrt{TS}$) as an estimate of the significance, with

$$TS = 2 \times \ln \frac{\mathcal{L}(H_1)}{\mathcal{L}(H_0)}, \qquad (4.4)$$

where $\mathcal{L}(H_1)$ is the Poisson likelihood associated with the assumed emission model (i.e., including the cluster) and $\mathcal{L}(H_0)$ corresponds to the Poisson likelihood of the null-hypothesis (i.e., without including the cluster) [222]. We perform this test for the different diffuse cluster models. We perform several simulations for each case to obtain a measurement of the uncertainty on the recovered significance.

**Results and discussions.** In figure 7, we show the expected significance of the ON-OFF analysis as a function of $\theta_{PS}$ and $N_{PSF}$ in the case of $\theta_{pointing} = 1.5$ deg. Note that the point-source mask aperture cannot be larger than the size of the ON region because NGC 1275 sits at the center of the cluster. Although the best significance always corresponds to no point-source mask ($\theta_{PS} = 0$), because the point source was accounted for in the background, we can observe that the recovered significance as a function of the mask and the ON region aperture highly depends on the cluster signal. This implies that the precision in the modeling of NGC 1275 will be more critical as the cluster diffuse emission is more compact (higher value of $\eta_{CRp}$). For instance, when fixing the ON region aperture to its optimal value, masking the point source with $N_{PSF} = 1$ will reduce the significance by about 1% for $\eta_{CRp} = 0$ and about 50% for $\eta_{CRp} = 1.5$. In practice, the optimization of the mask aperture will also depend on how well the central point source, NGC 1275, can be modeled. For instance, understanding the tails of the PSF could be a major issue when searching for diffuse emission, but we leave these investigations for future work, when the real CTA data will be available. We also observe that the best ON region aperture increases with the extent of the cluster (parameter $\eta_{CRp}$), going from about 0.4 deg to 1 deg for the considered cases and using 1.5 deg as the pointing offset. We also compute the maximal value of the significance as a function of the pointing offset, as shown in figure 8, left panel. As expected, the shape of the curve depends on the diffuse cluster model and we can observe maximal values between about 0.2 deg and 1.5 deg for all the cases.
<figure><img src="image_10.png" /></figure>

**Figure 7.** Normalized signal-to-noise ratio as a function of the radius of the ON region and the radius used to define the point-source mask (in units of $N_{\text{PSF}}$), for different cluster models defined via the $\eta_{\text{CRp}}$ parameter. Black lines correspond to 90%, 99% and 99.9% of the peak significance. The pointing offset was set to 1.5 deg in the present case.

<figure><img src="image_11.png" /></figure>

**Figure 8.** Normalized significance as a function of the pointing offset in the case of the ON-OFF analysis (left) and the template analysis (right).

The results of the template fitting analysis are presented in figure 8, right panel, where we can observe the normalized significance as a function of the pointing offset for all four cluster models considered (fixing $X_{500} = 10^{-2}$, $\alpha_{\text{CRp}} = 2.3$, and varying the spatial parameter $\eta_{\text{CRp}} = (0, 0.5, 1, 1.5)$). The significance smoothly decreases with the pointing offset, depending on the combination of the effective area and the instrumental background. Unlike for the ON-OFF analysis, no significant significance decrease is observed for small pointing offsets. We note that the background associated with NGC 1275, relative to the cluster signal, does
not depend on the offset and acts as an extra background. For this reason, the sensitivity as a function of the radius is flattened by the presence of the central galaxy. We do not observe any significant impact of the cluster diffuse model in this case. The significance remains relatively constant up to a radius of about 1 deg and vanishes beyond.

According to the results presented above and in order to allow for ON-OFF analyses while maximizing the significance of the template analyses, the pointing offset will be set to 1.0 deg in the following. This would allow us to perform both an ON-OFF and a template-based analysis after the actual observations have been made without a significant lose of sensitivity in each case.

# 5 CTA sensitivity to CR induced γ-ray emission from Perseus

In this section, we use the cluster models defined in section 2 together with the observation setup and background sky defined in section 4 in order to predict the CTA sensitivity to CR induced γ-ray emission from the Perseus cluster. This is done via the use of the publicly available package KESACCO (Keen Event Simulation and Analysis for CTA Cluster Observations).<fnref n="25" /> KESACCO is a python package that is based on the ctools software [220], and dedicated to perform the simulation and analysis of galaxy clusters observed with CTA. In appendix F, we validate our results using the gammapy software [223]. We first discuss the data preparation. We then quantify the degeneracy between the different signal and background components. Finally, we perform the analysis of simulated data using different methods and under different assumption for the signal to address the sensitivity of CTA to the ICM induced γ-ray emission.

## 5.1 Data preparation

The analysis and results presented in this section are based on event files simulated according to the setup defined above. The events are binned in energy and according to their sky coordinates. We consider 30 bins in energy and 0.02 deg sky pixels to make sure that the PSF and spectrum sampling is sufficient, unless otherwise specified. The region of interest is centered on the Perseus cluster and is 3 deg wide in diameter. Although we consider eight different pointing positions, we use a stacked analysis in which all data from multiple observations are stacked into a single counts cube for each sky and energy bin. Unless otherwise stated, we focus on 3D template analysis as it is expected to be more appropriate in the case of such extended sources and given the presence of the two AGN. Indeed, this method will allow us to constrain the shape of the diffuse Perseus cluster γ-ray emission in a straightforward way.

The CTA sensitivity is expected to be best around a few TeV, but covers energies from 20 GeV to 300 TeV. Because of the EBL attenuation, the differential flux of all sources which we account for in this paper (all located at the same redshift) is expected to drop drastically around 30 TeV (see figure 2, top left panel). Therefore, we consider a maximum energy of 50 TeV for the analysis. We also conservatively increase the low energy threshold to 150 GeV. This choice is driven by the fact that we aim at constraining the spectrum and the shape of the

<footnote n="25">²⁵https://github.com/remi-adam/kesacco.</footnote>
<figure><img src="image_12.png" /><figcaption>**Figure 9.** Left panel: combined counts per pixel map, over the full energy range (150 GeV–50 TeV). The bottom left grey circle gives the PSF at 1 TeV, and also accounts for an extra 0.1 deg smoothing used for visual purpose. The grey crosses give the pointing coordinates, the green circles show the position of the point sources, the grey dashed lines show the radius $\theta_{500}$ and the location of the cluster center is indicated by the cyan cross. Right panel: counts per bin measured as a function of energy within one degree from the cluster center. Contributions from the instrumental background, the point sources and the cluster diffuse emission are shown. This is in the case of the baseline model. The normalized residual is also shown at the bottom of the figure.</figcaption></figure>

ICM induced emission in the presence of the steep spectrum source, NGC 1275, which strongly dominates at low energy where the PSF is also larger. Such a threshold will allow us to avoid any significant bias that could arise due to mismodeling, but its value could be adapted when the true data will be available. This choice will not apply to the DM search in section 6.

## 5.2 Sky components and degeneracy

We first investigate how the components of the sky model compare to each other and how degenerate they are. We focus on the region enclosed within 1 deg from the cluster center($\sim \theta_{500}$). This region contains the two point sources and a large fraction of the cluster flux. In figure 9 (left panel), we present the stacked sky map including both the astrophysical and instrumental background, computed over the full energy range. The location of the different pointing centers is shown, falling very close to the radius $\theta_{500}$. The location of the two point sources are indicated as green circles and they are both clearly visible in the map. The cluster diffuse emission (in the case of our baseline model here) is blended with NGC 1275 and cannot be clearly distinguished at this stage. Figure 9 (right panel) shows how the different model components compare to each other. The instrumental background is the dominant component, by about an order of magnitude, depending on the energy. NGC 1275 dominates at low energy and IC 310 dominates at high energy. The cluster emission is subdominant in this case (baseline cluster model) given the area that is considered. Note that the signal from IC 310 can be easily masked given its location.

Secondly, we perform the joint likelihood fit of all the model free parameters: normalization and spectral index of the instrumental background, normalization and spectral index of
<figure><img src="image_13.png" /><figcaption>**Figure 10.** Correlation matrix between the different free parameters included in our model. The cluster model corresponds to the baseline model.</figcaption></figure>

both NGC 1275 and IC 310, and amplitude of the cluster model. We extract the correlation matrix between these components as reported in figure 10 for the baseline cluster model. The parameters associated with IC 310 are not significantly degenerate with any other model components given its sky location. Additionally, we note that IC 310 is sufficiently far from the cluster center to be masked without losing significant information on the cluster, so that any mismodeling systematic effect (e.g., different spectral shape) should have a negligible impact on the results. On the other hand, small degeneracies are observed between the cluster amplitude and the instrumental background (correlation ~ 0.10), as well as the parameters associated with NGC 1275. We note that the degeneracy depends on the cluster model parameters. For instance, a more extended cluster presents less degeneracy with NGC 1275 but is more degenerate with the instrumental background. We also note that increasing the size of the region of interest leads to a slightly lower instrumental background degeneracy because more leverage is available. The degeneracy also depends on the energy range used to select the data.

According to the degeneracy between the model components and the amplitude of these components within the central region of the field of view, we expect that the instrumental background will lead to most of the uncertainties associated with the cluster model constraints that can be extracted. The presence of NGC 1275 is also significantly affecting the properties that can be extracted for the Perseus cluster.
<figure><img src="image_14.png" /><figcaption>**Figure 11.** Exclusion limit at 95% confidence level for the parameter $X_{CRp}$, as a function of $\alpha_{CRp}$ and for different values of $\eta_{CRp}$, assuming non-detection of the signal. The shaded area represent the 68% confidence interval. The limits obtained by MAGIC [91] are reported for their isobaric model (roughly corresponding to our $\eta_{CRp} = 1$ model) and their extended model (roughly corresponding to our $\eta_{CRp} = 0.5$ model). We also report the baseline model and the best-fit parameters in the case of the pure hadronic scenario (where $\eta_{CRp} = 0.78$, and given the reference magnetic field model based on [160]) in black. To highlight the statistical uncertainties in the pure hadronic model, a set of 500 samples randomly extracted from the MCMC chains is shown, color-coded by the value of $\eta_{CRp}$, together with the 68 and 95% confidence intervals reported as grey lines. The change of the pure hadronic model constraint in the parameter space for different magnetic field models is given in figure 26. It would imply a shift along the $X_{CRp}$ axis by a factor of about 0.8 and 1.7, for the model given by [161] (upper $\eta$ scaling) and [126], respectively.</figcaption></figure>

## 5.3 Probing the parameter space with CTA

We first estimate the parameter space constraint that can be obtained with CTA. We perform event simulations without including the cluster signal, but still account for it afterwards.

We use our `ctools` framework to fit for the cluster together with the other sky model components and extract the upper limit on the cluster flux normalization. We also extract the corresponding upper limit on the overall cluster $\gamma$-ray flux. This is reproduced for 50 simulations in order to compute the mean upper limit and its standard deviation (implying an uncertainty on the mean and standard deviation of the upper limit of less than 5%, see appendix E), and for the different models that we test: spectral indices $\alpha_{CRp} = [2.0, 2.2, 2.4, 2.6, 2.8, 3.0]$ and spatial scalings $\eta_{CRp} = [1.5, 1.0, 0.5, 0.0]$. A higher spatial scaling indicates a more compact cluster, as can be seen by referring to section 2 and figure 2.

Figures 11 and 12 show the expected parameter space constraints in case of non-detection, and the flux upper limit in the range [0.15–50] TeV for 285 hours of observations, respectively.
<figure><img src="image_15.png" /><figcaption>**Figure 12.** Same as figure 11 for the flux upper limit for energies in the range [0.15, 50] TeV, with the same color code.</figcaption></figure>

In both figures, the exclusion limit at 95% confidence interval is shown, with the shaded areas representing the standard deviation. Figure 11 shows the upper limit in γ-ray flux as a function of the spectral index $\alpha_{CRp}$, for different spatial scalings. A flux upper limit down to $10^{-13} \text{ cm}^{-2} \text{ s}^{-1}$ is expected with CTA. Figure 12 shows $X_{CRp}(R_{500})$, the CR to thermal energy density ratio, as a function of spectral index. In the case of non-detection we are able to put a constraint down to $X_{CRp}(R_{500}) \sim 10^{-4}$ for a very compact CRp distribution. The most recent constraints put forth on diffuse γ-ray emission from Perseus come from the MAGIC telescope observations and analysis by [91]. As mentioned in section 2.3.3, the γ-ray and thermal emission modelling used by MAGIC have overestimated the thermal energy, meaning that the MAGIC constraints in the left panel are actually too optimistic. Either way, it is evident that CTA will be able to put a constraint that is about an order of magnitude deeper than MAGIC.

Both figures show that a more compact cluster would give better constraints. From the previous sections we know that in the CTA energy range, a lower spectral index gives more diffuse γ-ray emission (see figure 2), hence the non-detection of diffuse γ-ray emission puts a better constraint for lower spectral index numbers.

We also show that CTA will be able to exclude the pure hadronic best-fit model, assuming [160] magnetic field model with $\eta_B = 2/3$. Accounting for statistical uncertainties and assuming the most pessimistic magnetic field model (about 20% change, see section B), CTA should still be able to exclude significantly most of the parameters allowed in the pure hadronic scenario, although a small region of the parameter space with high $\alpha_{CRp}$ and low $\eta_{CRp}$ may still be viable in the case of non-detection, despite the higher value of $X_{500}$ that would be implied.
## 5.4 Measuring the cosmic-ray properties in the case of detection

In this section, we consider the baseline and the pure hadronic models, and produce the corresponding event simulations. We then investigate the sensitivity of CTA to constrain the model parameters. We consider two different scenarios. In the first one (section 5.4.1), we focus only on the recovered spectral properties by constraining the background model independently from the ICM and assuming that the cluster spatial shape is known. Then, in section 5.4.2, we consider simultaneously both the spectral and spatial properties together with all background components. The two methods are completely independent from each other. In the following, we analyse and discuss representative simulations of the data corresponding to our two models.

### 5.4.1 Spectral constraints only

We first focus on constraining the spectral energy distribution only. We start by using a likelihood fit to constrain the Perseus cluster normalization and all the other free sky parameters (see figure 10). To obtain the detection significance, we apply the Test Statistics defined in equation (4.4). As a reference, we obtain a test statistic value of $\sqrt{TS} \sim 6.5$ in the case of the baseline model, and $\sqrt{TS} \sim 11$ for the pure hadronic model. We then extract the spectrum of the diffuse cluster component thanks to the dedicated `csspec` function (from `ctools`) using the true input cluster spatial template to do so, which we keep fixed. The instrumental background is left free but the point sources are kept fixed to their maximum likelihood values. This procedure allows us to extract the Perseus cluster flux normalization in each independent energy bin. Additionally, the `csspec` function is used to recover the full likelihood scan for the normalization in each bin.

In order to measure the posterior likelihood in the parameter space, we employ an MCMC approach using the `emcee` package [173], following the method used by [93]. The sampling is performed according to the log likelihood function defined as

$$ \ln \mathcal{L}(\vec{\theta}|D) = \sum_i \text{interp} \left( \ln \mathcal{L}_i(F_{\text{scan}}, F_i(\vec{\theta})) \right). \quad (5.1) $$

The parameters $\vec{\theta} = (X_{500}, \alpha_{CRp})$ correspond to the normalization and the CRp spectral index, respectively, and $D$ to the data points. In each energy bin $i$, $\mathcal{L}_i$ is the likelihood scan extracted with `csspec` as a function of the flux $F_{\text{scan}}$. This quantity is interpolated at the location of any model to be tested against the data, $F_i(\vec{\theta})$, and summed over all bins. We use flat priors ($X_{500} > 0, 2 < \alpha_{CRp} < 5$) on the parameters and check a posteriori that the limits do not affect our results.

Once the MCMC has converged, we remove the burn-in phase and the chains provide us also with the posterior probability distribution function. We also recompute 100 SED models using parameters randomly sampled from the chains and measure the median and 68% confidence limit envelop on the recovered spectrum. The best-fit model is computed using the parameters that correspond to the maximum likelihood point in the chains.

In figure 13, we present the recovered spectrum and the constraints in the parameter space obtained in the case of the baseline and pure hadronic best-fit models. In both cases the best-fit describes the data well, as can be observed in the residual. The constraints in
<figure><img src="image_16.png" /><figcaption>**Figure 13.** Spectral constraints. The left panels correspond to the baseline model while the right panels correspond to the pure hadronic best-fit model (see table 1). **Top panels:** recovered spectrum (red points) and their error bars (according to the likelihood scan curvature), best-fit model (black), 68% confidence level (blue region) and a set of 100 realizations (light blue). The residual normalized by the error bars, $\chi$, is shown at the bottom. The dashed lines correspond to $\chi = \pm 2$. **Bottom panels:** constraints on the parameter space. Contours provide the 68% and 95% confidence region and the black star represents the injected model. The blue shaded region on the one dimensional histograms provides the 68% confidence region after marginalization.</figcaption></figure>

the parameter space show that both input models are recovered with the 68% confidence interval. The pure hadronic model allows us to better constrain the model given the higher signal-to-noise ratio. We note that the normalization and the spectral index of the CRp distributions are highly degenerate. This is particularly true for CTA because no leverage is available near the spectral bump around 1 GeV. Because of this degeneracy, an amplitude up to four times the input value is still allowed in the case of the baseline model. This reduces to about twice the input value in the case of the pure hadronic model, for which the signal is stronger. The error on the spectral slope is about 0.1 and 0.07 for the baseline and pure hadronic models, respectively. While the details are beyond the scope of this paper, we note that combining lower energy data (e.g., Fermi-LAT) to CTA, even in the case of non-detection, could help break the parameter degeneracy.
### 5.4.2 Joint spectral-imaging constraints

The analysis described in section 5.4.1 is standard, but intrinsically limited by the fact that the cluster profile is a priori unknown. This might lead to significant systematic effects when recovering the spectral energy distribution if it is improperly modeled. In addition, the uncertainties in the background (diffuse and point sources) are expected to affect the constraints on the cluster diffuse emission, as discussed in section 5.2. Although the error bars on the extracted cluster spectrum are expected to account for the correlation between the different components, the respective parameters are not co-varied when constraining the model to reduce the computing time. To mitigate these possible biases, we consider the joint-fit of the full sky model simultaneously, as described below. This allows us to sample the full parameter space. While the analysis described in section 5.4.1 is relatively quick, the one described hereafter is significantly more demanding in term of computing time, and the two are therefore complementary.

We compute the IRF convolved model counts cube, in spatial and energy bin $i$, as

$$M_i(\vec{\theta}) = C_i(X_{500}, \alpha_{CRp}, \eta_{CRp}) + \sum_{j \in [1,2]} \mathrm{PS}_i^{(j)}(A_{\mathrm{PS}}^{(j)}, \alpha_{\mathrm{PS}}^{(j)}) + B_i(A_{\mathrm{bkg}}, \alpha_{\mathrm{bkg}}). \quad (5.2)$$

The component $C_i(X_{500}, \alpha_{CRp}, \eta_{CRp})$ corresponds to the cluster model, which depends on the parameters $X_{500}$, $\alpha_{CRp}$, and $\eta_{CRp}$. The components $\mathrm{PS}_i^{(1,2)}(A_{\mathrm{PS}}^{(1,2)}, \alpha_{\mathrm{PS}}^{(1,2)})$ represent the point-source contribution, NGC 1275 and IC 310, respectively. The component $B_i(A_{\mathrm{bkg}}, \alpha_{\mathrm{bkg}})$ is the diffuse instrumental background, for which we assume a fixed shape with a floating spectral index and normalization. The parameter space includes nine parameters,

$$\vec{\theta} \equiv (X_{500}, \alpha_{CRp}, \eta_{CRp}, A_{\mathrm{PS}}^{(1,2)}, \alpha_{\mathrm{PS}}^{(1,2)}, A_{\mathrm{bkg}}, \alpha_{\mathrm{bkg}}) \quad (5.3)$$

As done in section 5.4.1, we use an MCMC method to sample the parameter space. The likelihood function is defined according to Poisson statistics as

$$\ln\mathcal{L}(\vec{\theta}|D) = \sum_{i} \tilde{M}_{i}(\vec{\theta}) - d_{i}\ln(\tilde{M}_{i}(\vec{\theta})), \quad (5.4)$$

with $d_i$ the data counts cube. Since the convolution of a model with the IRF is computationally expensive, it is not possible to compute $M_i(\vec{\theta})$ for all the MCMC test parameters. Therefore, we compute the models associated to all the components on predefined grids and $\tilde{M}_i(\vec{\theta})$ corresponds to the model $M_i(\vec{\theta})$ after interpolation at the exact location of the requested parameters. Each model component amplitude only accounts for a linear scaling and does not need to be included into the grid. The cluster model grid is thus bi-dimensional, defined as a function of the non-trivial parameters $\alpha_{CRp}$ and $\eta_{CRp}$ (typically 10 × 10 points allows us to recover any model with percent level precision). For the point sources and the instrumental background, the slopes $\alpha_{PS,bkg}$ are the only relevant parameters and the grids are therefore one-dimensional. The range over which the model grids are computed correspond to flat prior and we check a posteriori that it does not affect our result. We note that an extra prior is used for the cluster normalization, as $X_{500} \in [0, 20 \times X_{500}^{(\mathrm{input model})}]$ to avoid sampling nonphysical models in the case the signal-to-noise ratio is too low to provide a reliable constraint on this parameter. In order to highlight the impact of the point sources on the
<table><caption>**Table 6.** Recovered parameters of the joint spectral-imaging analysis for the baseline model, baseline model plus prior on NGC 1275, and the pure hadronic best-fit model cases.</caption><thead><tr><th>X<sub>CRp</sub><br>X<sub>CRp,0</sub></th><th>&eta;<sub>CRp</sub></th><th>&alpha;<sub>CRp</sub></th><th>A<sub>bkg</sub></th><th>&Delta;&alpha;<sub>bkg</sub></th><th>A<sub>ps,1</sub></th><th>&Delta;&alpha;<sub>ps,1</sub></th><th>A<sub>ps,2</sub></th><th>&Delta;&alpha;<sub>ps,2</sub></th></tr></thead><tbody><tr><td colspan="9" style="text-align: center; font-weight: bold;">Baseline model</td></tr><tr><td>1.69<sup>+3.59</sup><sub>-1.21</sub></td><td>1.17<sup>+0.25</sup><sub>-0.24</sub></td><td>2.45<sup>+0.18</sup><sub>-0.20</sub></td><td>0.9992<sup>+0.0016</sup><sub>-0.0016</sub></td><td>-0.0003<sup>+0.0011</sup><sub>-0.0011</sub></td><td>0.932<sup>+0.042</sup><sub>-0.060</sub></td><td>0.027<sup>+0.079</sup><sub>-0.102</sub></td><td>1.010<sup>+0.011</sup><sub>-0.011</sub></td><td>-0.010<sup>+0.010</sup><sub>-0.010</sub></td></tr><tr><td colspan="9" style="text-align: center; font-weight: bold;">Baseline model plus NGC 1275 prior</td></tr><tr><td>0.89<sup>+1.48</sup><sub>-0.57</sub></td><td>1.10<sup>+0.19</sup><sub>-0.17</sub></td><td>2.35<sup>+0.14</sup><sub>-0.16</sub></td><td>0.9993<sup>+0.0015</sup><sub>-0.0016</sub></td><td>-0.0002<sup>+0.0010</sup><sub>-0.0011</sub></td><td>0.963<sup>+0.030</sup><sub>-0.033</sub></td><td>0.019<sup>+0.043</sup><sub>-0.043</sub></td><td>1.007<sup>+0.010</sup><sub>-0.010</sub></td><td>-0.010<sup>+0.010</sup><sub>-0.010</sub></td></tr><tr><td colspan="9" style="text-align: center; font-weight: bold;">Pure hadronic best-fit model</td></tr><tr><td>0.85<sup>+0.77</sup><sub>-0.38</sub></td><td>0.87<sup>+0.11</sup><sub>-0.11</sub></td><td>2.38<sup>+0.10</sup><sub>-0.09</sub></td><td>1.000<sup>+0.0015</sup><sub>-0.0015</sub></td><td>-0.0014<sup>+0.0010</sup><sub>-0.0010</sub></td><td>0.969<sup>+0.037</sup><sub>-0.043</sub></td><td>0.0078<sup>+0.071</sup><sub>-0.081</sub></td><td>1.011<sup>+0.011</sup><sub>-0.011</sub></td><td>-0.015<sup>+0.010</sup><sub>-0.010</sub></td></tr></tbody></table>

cluster parameters constraint, we also perform the fit using a prior on the point-source parameters. The latter is Gaussian, centered on the input parameter values and with a standard deviation corresponding to 5% of the parameter values.

Figure 14 presents the posterior distribution in the full parameter space for the baseline cluster model (the case of the pure hadronic model is available in appendix G). Table 6 also reports the recovered parameters and their uncertainties in the different cases. The cluster normalisation, the point sources and the instrumental background parameters are shown relative to their input values. We can see that for this simulated data, the input model is recovered within 95% confidence interval for all the parameters (and 68% for most of them). While the observed posterior distribution is non-trivial for most of the parameters, we have tested that in the limit of high S/N, the parameter space is well described by a multivariate Gaussian function. As discussed in the case of the spectral only analysis (figure 13), strong degeneracy is observed between the cluster normalization and the CRp slope. On the other hand, the cluster profile parameter is not strongly degenerate with the other cluster parameters, but is anti-correlated with the parameters associated with NGC 1275. Indeed, as the cluster profile gets more compact, it resembles more to the point source so that increasing $\eta_{CRp}$ leads to a smaller point-source normalization. As the amplitude and slope of NGC 1275 are positively correlated, the same trend is observed between $\eta_{CRp}$ and $\alpha_{PS}^{(1)}$. The cluster normalization is also anti-correlated with NGC 1275 parameters for similar reasons. As already pointed out in section 5.2, no significant correlation is observed between the cluster model and IC 310, as expected given its sky coordinates. IC 310 could also be masked without losing significant information on the cluster properties and any mismodeling of its spectral shape is not expected to affect our results. The instrumental background normalization and slope is slightly correlated with the cluster parameters, in particular with the profile parameter, since lowering the value of $\eta_{CRp}$ leads to a flatter signal that slightly mimics the instrumental background. The correlation is thus positive between $\eta_{CRp}$ and $A_{bkg}$. As expected, the amplitude and spectral slope of the instrumental background and point sources are correlated. Compared to the spectral constraints only, the error bars on the parameter $\alpha_{CRp}$ increase from 0.10 to 0.17, highlighting the importance of constraining simultaneously all the components. The error on the cluster profile parameter, $\eta_{CRp}$, is about 0.26. The normalization is constrained to $X_{500} \lesssim 0.1$ (at 68% limit), which is significantly larger than in the case of spectral constraints only.
<figure><img src="image_17.png" /><figcaption>**Figure 14.** Posterior constraint on the full parameter space in the case of the baseline model. Contours provide the 68% and 95% confidence interval. The black star shows the input model parameters. The constraints are reported in the standard case (blue) and in the case where we assume a 5% prior on the flux and spectral indices of point sources (orange). The CRp normalisation is given relative to the input one, $X_{CRp,0} = 10^{-2}$.</figcaption></figure>

When using the prior on point source, we observe a strong reduction of the volume allowed by the cluster parameters (especially avoiding too high values of the normalization), in agreement with the degeneracy between components. On the other hand, the instrumental background parameters are nearly unchanged. The posterior constraints on NGC 1275 parameters is strongly improved, as can be expected due to the prior, but the constraints on IC 310 is nearly unchanged because the CTA data already provide a constraint better than 5% on the corresponding parameters.

We note that the steep spectrum of NGC 1275 makes it very sensitive to any uncertainty in IRF convolution of the model, which in turns may affect the cluster reconstruction when
<figure><img src="image_18.png" /></figure>

**Figure 15.** Comparison between the data and the best-fit model, for the full energy range. Units are in counts per pixel. *Top left panel:* data (baseline model only). *Top right panel:* best-fit model (baseline model only). *Bottom and middle left panel:* residual between the data and the background components (including point sources) for the baseline and pure hadronic models, respectively. *Bottom and middle right panel:* total residual between the data and the model for the baseline and the pure hadronic models, respectively. Contours provide the signal-to-noise ratio with $2\sigma$ spacing. The maps where smoothed with a Gaussian kernel with 0.15 deg for visualization purpose.

all parameters are allowed to vary, due to degeneracy. This may cause notable systematic effects if improperly accounted for in the case of a low significance, but we leave its detailed investigation for future work when the true data will be taken.

In figure 15, we show the input data, model and best-fit residual (including the cluster or not) after summing over all the energy bins. We can see that the best-fit model describes
<figure><img src="image_19.png" /></figure>

**Figure 16.** Left panel: cluster diffuse emission counts recovered within an aperture radius of 1 deg, after instrumental background and point-source subtraction. Right panel: cluster diffuse emission profile computed over all the energy bins, after instrumental background and point-source subtraction. The top panel is in the case of the baseline model and the bottom panel for the best-fit pure hadronic model.

the data well. The signal is dominated by the instrumental background on large scales, and by the point sources on small scales. Once these components are subtracted, we can observe the faint signal associated with the cluster emission.

In figure 16, we show the recovered cluster diffuse emission spectrum (left) and profile (right) after subtracting the point sources and the instrumental background. The spectrum is computed within 1 deg radius, which nearly corresponds to $\theta_{500}$ and the profile is computed by summing over the energy bins. The spectrum is overall well constrained over all the CTA energy range for this specific model, allowing us to recover the spectral distribution of the CRp. The cluster emission is detected on the profile up to about 0.2 deg radius (given the binning of 0.05 deg per radial bin), which allows us to constrain the shape of the signal well.

## 5.5 Discussions

The CTA observations are expected to provide unprecedented constraints on the diffuse γ-ray emission from the Perseus cluster. In this section we review its sensitivity to cluster scale CR physics and discuss the assumptions of the present analysis and caveats.
<figure><img src="image_20.png" /></figure>

**Figure 17.** Comparison between the CTA simulated map (background subtracted, pure hadronic best-fit model case, with [160] plus $\eta_B = 2/3$ magnetic field) and multi-wavelength data. **Left panel:** CTA surface brightness together with signal-to-noise ratio (with $2\sigma$ steps). The map was smoothed with a 0.15 deg Gaussian for visual purposes and the bottom left grey circle accounts for the effective map resolution. The top left subplot shows a zoom comparison between CTA and the radio mini-halo as imaged by the VLA at 1.4 GHz (extracted from the NASA/IPAC extragalactic database). VLA contours are logarithmically spaced. **Top right panel:** ROSAT [224] X-ray image extracted from the ROSAT All Sky Survey. Contours are logarithmically spaced. **Bottom right panel:** Compton parameter from [225]. The central region appears negative due to the strong contamination from NGC 1275, and another radio source is visible on the West. Contours are linearly spaced.

**Toward a multi-wavelength view of the Perseus cluster from radio to $\gamma$-rays.**
We first highlight how CTA observations would compare to data at other wavelengths, in particular in terms of the scales that are probed, in figure 17. We consider the most optimistic case that is the pure hadronic model (see figure 15 for the baseline model). We compare a realization of the CTA (background subtracted) image to radio data (VLA, 1.4 GHz), X-ray (ROSAT All Sky Survey), and tSZ (Planck). The CTA signal, while slightly resolved, is expected to be relatively compact and is only detected up to about $\theta_{500}/4$. The mini-halo is also very compact (at least at 1.4 GHz), slightly extended when compared to the CTA image angular resolution. In contrast, the X-ray and tSZ images that trace the thermal gas density squared and the pressure, respectively. We note that the tSZ signal is contaminated by two radio galaxies that appear negative on the map, including NGC 1275 at the center.
The γ-ray signal morphology is expected to compare well with X-ray data assuming that the CR follow the thermal gas well. In contrast, if the CR follow a distribution comparable to the temperature it would resemble more the tSZ image ($\propto P_e = n_e k_B T$).

**The energy budget of the Perseus cluster.** The cluster ICM is dominated by its thermal component, but non-thermal contributions (turbulence, magnetic fields, cosmic-rays) should be non-negligible and are only poorly known to date. The amplitude of the CTA detection (or upper limit) will provide a direct measurement (or a further constraint) of the energy stored in the CR component, and thus the CR pressure. Noteworthy, this will also provide a measurement of one of the contributions to the hydrostatic mass bias, which is one of the key ingredients for cluster cosmology (see [174] for a review). Nonetheless, even in the case of high significance detection, the CTA constraints on this parameter would have large uncertainties due to the degeneracy with other parameters (see figure 13 and 14), unless strong priors are available for the shape and the spectrum of CRp via other wavelengths. This could be the case thanks to high quality radio observations, assuming specific scenario (e.g., pure hadronic model). The use of lower energy data, such as Fermi-LAT, could also help breaking the degeneracy between the parameters of the cluster components, but also provide a useful prior on the spectrum of NGC 1275. The measurement of the CR to thermal energy will be easier in the case of a harder CRp spectrum and a more compact profile.

**Cosmic-ray acceleration mechanism and transport.** If the cluster is detected, the CTA data will allow us constraining the spectral index of the CRp that induce the γ-ray emission. The CRp spectrum (and amplitude) reflects the acceleration mechanism at play and gives estimates of the CR acceleration efficiency. Such measurement will therefore have important implications to discriminate between the various models that remain poorly constrained to date [3, 21]. The CTA sensitivity to diffuse emission will, however, loosen as $\alpha_{CRp}$ increases. For instance, assuming that the CR spatial distribution follows that of the thermal gas, it will not be possible to detect the cluster for a CR slope larger than 2.6 if the CR to thermal energy ratio is lower than about 5% (figure 11 and 12). Assuming that the CR pressure support within $R_{500}$ is less than 20%, as can be reasonably expected (e.g., [174]), CTA will be limited to probe $\alpha_{CRp} < [2.6, 2.9]$ depending on the spatial distribution ($\eta_{CRp} \in [0, 1.5]$).

Thanks to the resolved observations provided by CTA, the data will also be sensitive to the spatial profile of the γ-ray emission, which is related to the spatial distribution of the CRp in the cluster. In turn, the cluster CR profile will give precious information about the production sites and transport processes (advection and turbulent motion imply a more compact profile, while diffusion and streaming flatten it). In this paper, we have chosen to model the CR profile using a scaling with respect to the thermal gas density, which allows us to test the effect of transport via a single parameter, $\eta_{CRp}$. We stress that this choice is not unique and the present work could be extended to test any model. However, given the current limited knowledge about the CR production and transport in the ICM, simple models might already be sufficient for such investigations. The precision on the CRp spatial distribution will depend on CR to thermal energy ratio, the spectral distribution of CR, and the shape of the CR profile itself since a more compact distribution will lead to a larger signal-to-noise ratio. For hard spectra, CTA should be able to detect the diffuse Perseus
emission even in the case of flat CR distributions, provided that the CR to thermal energy ratio is sufficient ($X_{500} \gtrsim 10^{-2}$).

**Probing the time integrated AGN and starburst activity in Perseus.** In this paper, we have considered that the cluster ICM was filled with CRp following a given spectral and spatial distribution, regardless of their detailed origin. While they may arise from shocks and turbulence, it is also possible that accumulated injection from compact sources play an important role [226]. Therefore, the CTA observations will be sensitive to the time integrated AGN cosmic-ray outburst and starburst activity. Since the number of AGNs and starbursts in the Perseus cluster can be estimated, this will provide a unique way to estimate the average energy injection via CRp from such sources in galaxy clusters. The level of γ-ray emission will depend on how CR from AGN and starbursts are processed, so that acceleration, transport and injection are closely linked to each other.

**Constraints on the magnetic field strength in the pure hadronic model.** The cluster non-thermal diffuse emission that we aim at probing with CTA is physically connected with the diffuse synchrotron emission observed at radio wavelengths, the latter being produced by the entangled combination of CRe and magnetic field (see section 2.2). Therefore, it is in principle possible to constrain the magnetic field strength using the joint analysis of radio and γ-ray data, assuming a given physical scenario such as the pure hadronic model (as done for the Coma cluster [99, 100]). In the present analysis, we have assumed a fixed magnetic field model ([160] plus $\eta_B = 2/3$ as a reference, with $<B(10\,\text{kpc})> \sim 25\,\mu\text{G}$) when considering the pure hadronic scenario. For instance, neglecting uncertainties arising from the radio data, our reference magnetic field model implies a test statistic $TS = 125$ ($\equiv 11\sigma$ detection) for CTA in the pure hadronic scenario. Lowering the γ-ray emission, so that the cluster is not detected, would imply lowering the amount of secondary CRe by the same amount. The magnetic field can then be adjusted so that the radio synchrotron emission remains the same.

In appendix B, we show how the different models of magnetic field considered in this paper affect our model parameters and translate into γ-ray emission (table 1). For example, the model giving the lowest magnetic field profile, by a factor of about 2.5 within 100 kpc from the center compared to our reference, implies a normalization that is 60% larger with nearly unchanged CR slope and spatial distribution, and therefore an increased γ-ray emission by about 60% as well. This shows that in turn, constraints can be obtained on the magnetic field strength when fixing the radio and γ-ray observables. Given the current uncertainties in the radio spatial and spectral distribution of the Perseus cluster, we only provide a qualitative discussion. However, we note that current and future radio observations<fnref n="26" /> should provide an unprecedented measurement of the diffuse radio emission from the Perseus cluster, which could then be combined with CTA observations to measure the magnetic field. This could be further combined with independent measurements of the magnetic field (e.g., Faraday rotation measure) to test CR physics and magnetic field in the ICM and particle acceleration mechanisms.

<footnote n="26">²⁶E.g., with the Square Kilometer Array (https://www.skatelescope.org/) and its precursors.</footnote>
**Analysis choices and CTA configuration.** The analysis presented in this paper rely on several choices that are expected to affect our sensitivity results. We discuss the most relevant ones hereafter.

Most of the analysis described in this paper has been performed using a template fitting analysis. Indeed, it was necessary in order to not only test for the detection of the cluster, but also model and constrain the spectral and spatial distribution of the emission. Because it uses more information than classical ON-OFF technique, the template fitting analysis is expected to provide more accurate constraints. On the other hand, it might be more sensitive to systematic effects due to the mismodelling of the data.

We decided to select events in the energy range between 150 GeV and 50 TeV. While the exact value of the upper bound does not affect the results, because no signal is expected above ~30 TeV, reducing the minimal energy is expected to increase statistics and improve the CTA constraints. On one hand, we use a conservative value that is relatively large compared to what can be expected with CTA. On the other hand, low energy data may be more difficult to model and in addition to the instrumental background, NGC 1275 strongly dominates at low energies. Therefore, it is likely that the improvement in sensitivity will not be major when pushing to lower energy thresholds.

We conservatively selected events within a 3 × 3 deg² region around the cluster. While no significant cluster diffusion emission is expected beyond this region, it is likely that the corresponding data will help constraining the instrumental background, provided that it can be properly modeled. This might, in turn, help breaking the cluster — instrumental background degeneracy and slightly improve the CTA constraints.

A new version of the CTA IRFs has been released since the analysis presented here was completed. They correspond to the initial array configuration, made of 4 LST plus 9 MST (the so-called CTAO Alpha configuration, prod5-IRFs) instead of the one with 4 LST plus 15 MST (prod3-IRFS) used here. The reduction in sensitivity is mostly effective at high energy ($E \gtrsim 500$ GeV), by a factor of up to 30%, but improved off-axis sensitivity. Therefore we would expect the results presented in this paper to be affected by such sensitivity reduction. Assuming the initial array configuration would imply a reduction of the upper limits (figure 11 and 12) and a widening of the uncertainties on the cluster recovered parameters (figures 13 and 14) by up to 30%, although the improved off-axis sensitivity should mitigate this effect. Cluster models with harder spectra (small $\alpha_{CRp}$) should be the most affected because a larger part of their flux arises from the high energy part of their spectrum.

# 6 CTA sensitivity to DM induced γ-ray emission from Perseus

In this section we investigate the sensitivity of CTA to the DM emission model presented in section 3 for the cluster. The results shown in the following have been obtained using the gammapy (v.0.18.2) open-source code [223].<fnref n="27" />

Perseus represents a complex environment where multiple γ-ray sources coexist. Indeed, with the current generations of IACTs, these astrophysical sources have posed large difficulties to the search for a putative DM signal in the area, as it is hard to avoid ‘contamination’ from

<footnote n="27">²⁷All the notebooks and scripts developed for the DM analysis can be obtained in https://github.com/peroju/dmtools_gammapy.</footnote>
the astrophysical sources and to disentangle between these and DM, both spectrally and
spatially. In contrast, due to its expected improvement in sensitivity, angular resolution and
larger FoV, CTA will allow for a much better discrimination by means of analysis techniques
that were traditionally discarded for this kind of telescopes and used only for all-sky γ-ray
facilities such as Fermi-LAT. In particular, in our DM search we perform a template-based
analysis, which allows to take into account the rest of expected γ-ray emissions as backgrounds,
and is in line with the analysis previously performed in section 5 for CRs and also with
recent studies in the Galactic center region by the CTA Consortium [122]. This type of
analysis, although computationally expensive, is expected to provide not only the highest DM
sensitivity but also to yield the most realistic one.<fnref n="28" /> Yet, mainly for the sake of comparison
with previous IACT works and in order to provide a first-order sensitivity study in our work
as well, we also carry out the so-called ON-OFF analysis, i.e. the most standard analysis
strategy adopted by some current IACTs. Full details of such analysis are left for appendix I
with the corresponding comparison with template-based results.

To perform the simulation of the DM-induced γ-ray emission from Perseus, we adopt the
observation setup described in section 4, i.e. a pointing with one degree offset from the cluster
center. Similarly to the CR analysis, we stack the data from the different pointings in a single
3D cube, containing the proposed 300h observation time. We adopt ten energy bins, starting
from 50 GeV up to 100 TeV. This allows us to exploit the whole range of CTA's sensitivity as
well as to overlap with the DM parameter space previously explored by the Fermi-LAT (at the
lowest considered energies) and by the H.E.S.S and MAGIC telescopes (in the TeV regime).
We use a spatial binning of 0.02 deg per pixel and our field of view diameter is set to 5 deg.

To model the DM distribution in the cluster, we use the 2D templates described in
section 3.3 for each of the three benchmark models in the case of annihilation and for the one
of decay. For annihilation, we simulate a putative γ-ray signal assuming a branching ratio
of 100% to the bbar annihilation channel, a DM mass m<sub>χ</sub> = 10 TeV<sup>29</sup> and a velocity-averaged
cross-section matching the thermal value of &lt; σv &gt; = 3 × 10<sup>-26</sup> cm<sup>3</sup> s<sup>-1</sup>. For decay, we
adopt τ<sub>χ</sub> = 10<sup>27</sup> s as the DM particle lifetime. The expected γ-ray contribution from the
existing AGNs in the area is accounted for via the models described in section 4.2, while the
CR-induced γ-ray emission is described using the spectral and spatial template corresponding
to the 'Baseline' benchmark model (section 2.3.1). We remark that those γ-rays originated
from both the AGNs and CRs are considered as backgrounds in our DM search, together
with the CTA instrumental background modelled by the IRFs. In total, we create four sets of
simulations, one per each DM template (MIN, MED, MAX and DEC — see table 3), and
always adopting the same combination of backgrounds. An example of a simulation can
be seen in figure 18 for the MED annihilation model. We also compare in the same figure
the spectra of the different γ-ray sources in Perseus to the ones for annihilation and decay.
This comparison shows that in the annihilation scenario the γ-ray flux is expected to be
orders of magnitude below the astrophysical backgrounds for realistic cross-section values.

<footnote n="28"><sup>28</sup>There is an on-going debate on the possible introduction of bias due to the use of theoretical models in the template-based methods.</footnote>

<footnote n="29"><sup>29</sup>The effect of the EBL on the DM simulated spectrum will be much more significant if we consider m<sub>χ</sub> > 30 TeV, in which case the characteristic DM cutoff would be substituted by the EBL cutoff due to γ-ray attenuation.</footnote>
<figure><img src="image_21.png" /><figcaption>**Figure 18.** Left panel: Simulation of an observation of the Perseus region with CTA. The color informs on the expected counts. This example is for the MED annihilation model, with a DM particle of $m_\chi = 10$ TeV, $b\bar{b}$ channel and $\langle \sigma v \rangle = 3 \times 10^{-26} \text{ cm}^2 \text{ s}^{-1}$. We adopted the ‘Baseline’ model for the CR-induced emission (section 2.3.1) and included the two AGNs in the region via the descriptions in equations (4.1) and (4.2). We only show the first four energy bins of our analysis, each of them corresponding to one panel. The white dashed circle goes over the different proposed pointings (all of them at 1 deg from the cluster’s center). Right panel: Comparison of the spectra of the different simulated $\gamma$-ray components integrated up to $R_{200}$. We also include the spectrum of a $m_\chi = 10$ TeV DM particle decaying via the $b\bar{b}$ channel. We recall that all these emissions do not happen through the entire CTA FoV.</figcaption></figure>

Yet, in the case of decay, we can have comparable fluxes (adopting $\tau_\chi = 10^{27}$ s, which is indeed very close to current constraints, [227]).

Finally, in order to obtain statistically meaningful results, we create 100 different simulations with the same observation setup, and use the corresponding mean for the data analysis itself. A study of the stability of our simulation results with respect to the number of simulations is included in appendix J, which shows that a reasonable compromise in terms of convergence and computation time is indeed achieved after a few tens of simulations in almost all cases. Thus, in total, we produce 400 simulations, 100 per each benchmark DM model (MIN, MED, MAX and DEC), which we now proceed to analyze.

## 6.1 Template-fitting analysis and DM sensitivity

To search for DM emission, we first define for convenience the normalization of the DM-induced $\gamma$-ray flux $A_\chi$, as $\langle \sigma v \rangle = A_\chi \times \langle \sigma v \rangle_{\text{thermal}}$, for the three annihilation models, and $\tau_\chi = \tau_{\text{ref}}/A_\chi$ (equation (3.2)), where $\tau_{\text{ref}} = 10^{27}$ s, for the decay scenario. We fit this normalization for each simulation using the `iminuit` [228] backend of `gammapy` to the theoretical emission of different DM particle candidates. To cover the available mass range for WIMPs, we perform the fit to fourteen DM masses across the defined energy range and for two representative annihilation/decay channels, i.e. $b\bar{b}$ and $\tau^+\tau^-$.<fnref n="30" />

<footnote n="30">³⁰Initially we also included the $W^+W^-$ channel and noticed that it produced similar results to $b\bar{b}$, as expected, so we decided not to show it in the following for the sake of clarity.</footnote>
Together with the DM normalization, we apply a joint-likelihood fitting including the relevant parameters of the considered astrophysical backgrounds. For example, since the CR-related flux is considered as a background for this analysis, we only allow for an overall normalization on the CR-induced γ-ray flux for the ‘Baseline’ model. Then, the complete flux model to fit the simulations is defined, in spatial and energy bin i, as:

$$M_i(\vec{\theta}) = DM_i(A_{\chi}) + CR_i(A_{CR}) + \sum_{j \in [1,2]} PS_i^{(j)}(A_{PS}^{(j)}, \alpha_{PS}^{(j)}) + B_i(A_{bkg}, \alpha_{bkg}). \quad (6.1)$$

In total, we fit eight different parameters:

$$\vec{\theta} = (A_{\chi}, A_{CR}, A_{PS}^{(1,2)}, \alpha_{PS}^{(1,2)}, A_{bkg}, \alpha_{bkg}), \quad (6.2)$$

where all the parameters have been previously defined (see section 5.4.2 and equation (5.2)). The likelihood function is the one corresponding to Poissonian data (equation (5.4)). We also adopt a flat prior on $A_{\chi}$, to allow only for zero or positive values.<fnref n="31" /> This joint-likelihood fit will enable to check for intrinsic correlations and to account for them when estimating the size of the uncertainties. To obtain the detection significance of the DM-induced γ-ray emission, we apply the Test Statistics as defined in equation (4.4). We define a detection when $TS = 25$, which roughly corresponds to $5\sigma$.

In neither of the MIN, MED, MAX annihilation or decay scenarios, and for none of the channels considered, a statistically significant detection is predicted. We proceed to compute the 95% confidence level (C.L.) upper limits (lower limits for decay) of the DM normalization parameter $A_{\chi}$ using the likelihood profile method [222, 229]. For the DM normalization, we assume a one-sided distribution (we recall the flat prior $A_{\chi} \ge 0$), thus the 95% C.L. corresponds to a change of $\Delta TS = 2.71$, with respect to the best fit.

## 6.2 Projected sensitivity to annihilating DM

Assuming the templates and spectra for DM annihilation (section 3), we compute the 95% C.L. upper limits for $A_{\chi}$ for all considered cases, convert this value to $< \sigma v >$ (the standard parameter to be constrained for annihilating DM) and average the results for 100 simulations.

In figure 19, we show the CTA projected sensitivity as 95% C.L. upper limits on the velocity-averaged cross-section versus the DM mass, for the MED annihilation model and the two considered representative annihilation channels. As it can be seen, the most constraining limits are obtained for the $\tau^+\tau^-$ channel, although also the $b\bar{b}$ channel yields similar results for DM masses above $\sim 10$ TeV. In all cases, our exclusion limits are more than $\sim O(10^2)$ above the thermal relic cross-section value. Yet, we note that they will be the most constraining ones (in the TeV energy range) considering galaxy clusters as DM targets. This can also be seen in both panels of figure 19, where we show a comparison to recent works. The latest cluster DM limits obtained by currently operating IACTs come from the observation of the Fornax cluster (14.5 hours) by the H.E.S.S Collaboration [115], and are around one order of magnitude weaker than our CTA predictions at a few TeV (i.e. at the peak of H.E.S.S. sensitivity). We also compare our CTA predictions with limits from Fermi-LAT in the sub-TeV WIMP mass range. In [112], authors analyze 12 years of Fermi-LAT data for a

<footnote n="31">³¹This method is usually known as the bounded likelihood method [229].</footnote>
<figure><img src="image_22.png" /><figcaption>**Figure 19.** Sensitivity of CTA to a DM annihilation signal from the Perseus cluster. Curves represent the 95% C.L. upper limits on the velocity-averaged cross-section versus the DM mass for the MED annihilation model. The green (yellow) band shows the $1\sigma$ (2$\sigma$) scatter of the projected limits. The black dashed line is the thermal relic cross-section ($<\sigma v>_{thermal} = 3 \times 10^{-26} \text{ cm}^3 \text{ s}^{-1}$). Left panel: Cross-section upper limits for the $b\bar{b}$ channel (right panel for $\tau^+\tau^-$ channel) in comparison with the most recent results on DM-annihilation searches in galaxy clusters using Fermi-LAT ([112] with orange dotted lines; [107] with blue dot-dashed lines) and H.E.S.S. ([115], with purple dashed lines). All of our results are shown for the template fitting analysis, unless stated otherwise.</figcaption></figure>

sample of five clusters and, in the absence of a signal, set constraints only for the $b\bar{b}$ channel. Perseus is the cluster yielding their weakest limits while Fornax gives the most constraining ones. The CTA sensitivity predictions for Perseus are always better than the obtained in [112] in the whole range but below 100 GeV, where both results are still comparable. But even for Fornax, CTA’s improvement in sensitivity will allow us to set the tightest DM constraints for DM masses above 1 TeV. In the case of the $\tau^+\tau^-$ annihilation channel, we can compare our results with those in [107], where authors combine 3 years of Fermi-LAT data from 8 galaxy clusters to derive the corresponding DM limits. As a consequence of using a considerable reduced amount of data compared to [112], their limits are considerably weaker in comparison and cannot compete to those from CTA. In conclusion, CTA DM limits from Perseus will be the most constraining ones from these class of targets in the TeV energy range.

By far, the largest uncertainty on the obtained cross-section upper limits is the modelling and contribution of cluster DM subhalos to the total $J_{\text{ann}}$-factor. This was the reason to build three different annihilation benchmark models in the first place (MIN, MED and MAX — section 3.2), each of them representative of very diverse levels of substructure and their contribution to the annihilation flux. To properly quantify the impact of this uncertainty in our limits, in addition to the 95% C.L mean upper limit for the MED case presented above, we also compute limits for the MIN and MAX benchmark models. These results are shown in figure 20 for the $\tau^+\tau^-$ annihilation channel, and reveal that our limits can be modified substantially depending on the considered subhalo scenario. In particular, the *boost factor* associated to our MED model ($B_{\text{MED}} = 9.2$) translates into an improvement of a factor $\sim 7$ with respect to the MIN benchmark model. In the case of the MAX model ($B_{\text{MAX}} = 59.3$), the limits improve up to $\mathcal{O}(10)$ times the MED model constraints.<fnref n="32" /> Previous studies have

<footnote n="32">³²Note that the improvement in the limits takes into account the precise impact of the subhalo boost on the</footnote>
<figure><img src="image_23.png" /><figcaption>**Figure 20.** 95% C.L. upper limits of the velocity-averaged cross-section versus the DM mass for the MIN (dot-dashed line), MED (solid) and MAX (dashed) (see table 4) subhalo models and the $\tau^+\tau^-$ annihilation channel. The green (yellow) band represents the $1\sigma$ (2$\sigma$) scatter of the projected limits. The black dashed line represents the thermal relic cross-section ($<\sigma v>_{thermal} = 3 \times 10^{-26} \text{ cm}^3 \text{ s}^{-1}$). All of our results are shown for the template fitting analysis, unless stated otherwise.</figcaption></figure>

also accounted for the subhalo contribution in their limits. In [112], for instance, authors show that their boosted DM model improves their limits only by a factor $\sim 1.5$, indeed a natural consequence of the low contribution of their substructure model to the annihilation flux. In earlier work by [107], authors gained up to $\sim \mathcal{O}(10^2)$ of improvement in their constraints due to subhalos, however they assumed boosts $B = 500-1200$, which from recent studies can be considered as overly optimistic (e.g., [193] and references therein). As for the constraints from the H.E.S.S Collaboration [115], their boosted MED ($B=10$) and MAX ($B=100$) models improved their limits by a factor $\sim 7$ and $\sim 75$, respectively, in reasonable agreement with our results (we note though that, even with these improvements in their limits, our CTA projected constraints for the MED model are more restrictive). All in all, our results in figure 20 illustrate the key role that halo substructure can have to discover/rule-out WIMP DM.

Other uncertainties affect the results in this section as well, namely the one coming from the scatter in the concentration-mass relation and the uncertainty in the estimation of Perseus mass. Yet, the level of these uncertainties ($\sigma_{J_{\text{ann}}} = 0.2$ dex; see section 3.3) is well below the one induced by the modelling of the subhalo population ($\sigma_{B_{\text{MAX}}} = 2$ dex). As a result, we note that the variation of upper limits shown in figure 20 among our three different subhalo models encompasses by far these other, second-order uncertainties. Finally, we remark that these results are for the use of the “Baseline” CR model, since different models for the CR densities may also impact our results (see appendix H).

<footnote>spatial morphology of the DM signal as well; thus the relation of this improvement with the boost value is not linear, as naively expected.</footnote>
<figure><img src="image_24.png" /></figure>

**Figure 21.** Sensitivity of CTA to a DM decay signal from the Perseus cluster, at 95% C.L., in terms of the mean lower limits of the lifetime of the DM particle versus the DM mass. The green (yellow) band shows the $1\sigma$ (2$\sigma$) scatter of the projected limits. Left panel: Mean lifetime lower limits for the $b\bar{b}$ channel (right panel for $\tau^+\tau^-$ channel) in comparison with the most recent results on DM decay in galaxy clusters using MAGIC data (olive long-dashed lines; [117]), Fermi-LAT data (blue dot-dashed lines; [107]) and H.E.S.S data (purple dashed lines; [116]). All of our results are shown for the template fitting analysis, unless stated otherwise.

## 6.3 Projected sensitivity to decaying DM

Similarly to the annihilating DM case, we compute 95% C.L. upper limits for $A_\chi$, convert this value to $\tau_\chi$ (the standard parameter to be constrained for decaying DM) and average the results of 100 simulations. We note that as $A_\chi$ and $\tau_\chi$ are inversely proportional, the upper limits on $A_\chi$ translating into lower limits for $\tau_\chi$.

We show, in figure 21, the projected sensitivity of CTA as the 95% lower limits on the DM particle lifetime versus the DM mass for our DEC model. As for the annihilation constraints, the most constraining limits are for the $\tau^+\tau^-$ channel, reaching similar values than the $b\bar{b}$ channel for masses above $\sim 10$ TeV. We also compare our CTA projections to the most up-to-date constraints on DM decay in clusters. The most recent work comes from the MAGIC Collaboration [117], where authors analyze 202h of data targeting the Perseus cluster. Their limits yield the tightest constraints for masses above $\sim 1$ TeV for the $\tau^+\tau^-$ channel, yet CTA will improve these limits by more than two orders of magnitude at lower masses and up to a factor 7 for masses higher than 5 TeV. This large improvement in the low mass range is probably twofold: i) the improvement of up to one order of magnitude in the sensitivity of CTA with respect to MAGIC$^{33}$ in the lower energy range; ii) the adoption of a mask of 0.1 deg in the MAGIC analysis to avoid contamination from NGC 1275, thus losing a considerable amount of data. Another comparable work from existing IACTs comes from [116], where authors analyze 14.5h of H.E.S.S. observations of the Fornax cluster. Though covering a smaller range of DM masses, their DM decay limits, both for the $b\bar{b}$ and $\tau^+\tau^-$ channels, are of the same order of magnitude than the MAGIC limits for Perseus [117]. Lastly, we also compare in figure 21 our CTA predictions to the constraints obtained using 3 years of Fermi-LAT data from 8 galaxy clusters [107]. These limits prevail up to DM masses of a few

<footnote>$^{33}$https://www.cta-observatory.org/science/ctao-performance/.</footnote>
<figure><img src="image_25.png" /><figcaption>**Figure 22.** 95% C.L. mean lower limits of the lifetime of the DM particle versus the DM mass (black lines) for $b\bar{b}$ (left panel) and $\tau^+\tau^-$ (right panel) decay channels, in comparison with the limits obtained by [227] (corresponding to NFW $5\sigma$ free source fits, in dashed red lines). The green (yellow) band shows the $1\sigma$ ($2\sigma$) scatter of the projected limits. All of our results are shown for the template fitting analysis, unless stated otherwise.</figcaption></figure>

TeV, where CTA starts to yield first comparable and then, quickly, stronger constraints for both decay channels. From all these comparisons, we conclude that CTA will not only be able to test an unexplored region in the DM decay parameter space, but should yield the most stringent constraints from $\gamma$-ray DM decay searches above 1 TeV.<fnref n="34" />

Finally, we compare in figure 22 our CTA projections to the most constraining DM decay results from $\gamma$-ray observations so far [227]. In the latter work, the *Fermi*-LAT Collaboration analyzes 1 year of LAT data from the Milky Way halo. The CTA exclusion curves, as already seen in previous comparisons with *Fermi*-LAT limits, become more constraining above DM masses around a few TeV (few hundreds of GeV for the case of the $\tau^+\tau^-$ channel).<fnref n="35" /> From these comparisons, we conclude that the null-detection of a DM decay signal in Perseus with CTA would still provide unprecedented constraints on DM decay models at the TeV mass scale.

# 7 Summary and conclusions

In this work, we have analyzed the CTA sensitivity to detect diffuse $\gamma$-ray emission from the Perseus galaxy cluster. We have considered the possible contribution from CR and from DM annihilation and decay. We have assumed 300 hours of observations as proposed by the CTA consortium as a key science project [121].

We built a CR-induced $\gamma$-ray emission model using the MINOT software [148] (section 2). Our model relies on a description of the thermal gas pressure and density, the magnetic field strength, and a parameterization of the spatial and spectral distribution of the CRp. We calibrated the model components according to available data in the literature (figure 1). The thermal part is expected to be accurate from the core to the outskirt and was kept fixed. The magnetic field was modeled by considering available measurements and different assumptions

<footnote n="34">³⁴The IceCube Collaboration provides more constraining limits above 100 TeV [230].</footnote>

<footnote n="35">³⁵*Fermi*-LAT observations of the Extragalactic $\gamma$-Ray Background [107] also yield DM decay limits that are comparable to the ones from the Milky Way Halo.</footnote>
and scaling to estimate the associated systematic uncertainty. The CRp parameters were calibrated according to several scenarios including a baseline model and the pure hadronic model (figure 1 and table 1). The description of the background sky was made according to the literature, noting that at least two point sources will affect the CTA observations (section 4.2). The observation setup was investigated focusing on the offset between the cluster center and the pointing (figure 8). Finally, the CTA sensitivity to diffuse γ-ray emission was investigated (section 5). We considered both the case of non-detection in which we computed the expected exclusion limit that CTA will obtain (figure 11), and the case of specific models for which a detection is expected. In the case of detection, we investigated the constraints that CTA will be able to provide regarding the spatial and spectral properties of the γ-ray diffuse emission (figures 13 and 14).

Our main findings with regard to γ-ray emission are as follows.

* The pure hadronic model, as calibrated using existing radio data, implies about 5% of CRp energy relative to the thermal energy within $R_{500}$, a CRp spectral index of ~2.3 and a CRp profile slightly shallower than the thermal gas density, assuming a magnetic field strength model based on existing data, although we note that parameters are strongly degenerate. Even when fixing the magnetic field strength, the uncertainty on the γ-ray flux corresponds to a factor of about two. This large uncertainty is mainly due to the limited spectral coverage of the available radio data that we used. The systematic uncertainty associated with the magnetic field implies an additional factor of a few in the uncertainty, depending on the radius.

* According to our modelling and the CTA IRF, we found that a pointing offset of about 1 deg was an optimal choice. This is obtained by considering both ON-OFF analysis and template-fitting techniques and for both CR and DM related analysis.

* In the case of non detection, we find that CTA should improve the current limits on the CRp content of the Perseus cluster by about an order of magnitude. Assuming a standard scenario ($\alpha_{CRp} = 2.3$ and $\eta_{CRp} = 1$), CTA should be able to constrain the CRp to thermal energy ratio within $R_{500}$, $X_{500}$, down to about a $3 \times 10^{-3}$.

* Assuming the pure hadronic model, CTA should allow us to detect the ICM induced diffuse emission with a high significance. The spectral index of the CRp should be constrained to an uncertainty of about ±0.1 and the spatial distribution down to about 10% precision.

CTA observations of Perseus will allow us to address fundamental questions related to the underlying mechanism that accelerates particles in the ICM (such as shocks, turbulent reacceleration, the direct injection of CR from AGN), the direct injection of CR from AGN. CTA will also allow us to test the physics associated with CR transport in the ICM. According to the results presented in this paper, CTA will provide unprecedented constraints on the physics associated with particle acceleration in galaxy clusters.

In our work, we also investigated the potential of CTA to search for DM in Perseus. To model its DM content, we assume the DM to be entirely composed of WIMPs. We consider two different DM scenarios, where the WIMP annihilates or decays into SM particles. We
build the DM density profile of the cluster accounting for the smooth DM distribution in the main halo plus the abundant DM in the form of substructures, or subhalos, that these massive objects are expected to host. The main halo is modelled starting from X-ray data, that allows for an estimate of the cluster’s mass. For that mass, we then build a tailored NFW profile for Perseus following results from numerical cosmological simulations at these scales. The contribution of the subhalo population to the expected DM flux is the main uncertainty in our DM annihilation model (subhalos do not play a role for decay), and simply reflects the existing ongoing debates in the literature on the precise properties of these objects. To bracket this uncertainty, we define three benchmark models (table 3), each of them representative of the different contribution that subhalos may have in the computation of the DM annihilation flux: the MIN model, where substructures are completely neglected; the MED model, our best guess according to the most recent results from numerical simulations; and the MAX model, defined to provide an upper limit to the role that substructures could play for the DM flux. From these, we compute the corresponding $J_{\text{ann}}$ factors for the MIN, MED and MAX model and $J_{\text{dec}}$ for the MIN case (summarized in table 4) and produce 2D spatial templates of the expected emissions (figure 5). The effect of substructures in the $J_{\text{ann}}$-factor is quantified in terms of the so-called subhalo *boost factor*. More precisely, we find $B_{\text{MED}} \sim 9$ and $B_{\text{MAX}} \sim 59$ for the MED and MAX models, respectively. As expected, substructure affects mainly the morphology and strength of the DM signal in the outskirts of the cluster (figure 4).

Adopting these DM templates and a particular DM spectrum (figure 18), we then create simulations of CTA observations of Perseus, that take also into account the γ-ray emission from the CR (baseline model), the two AGNs in the area (NGC 1275 and IC 310) and the instrumental background. All of these components act as background in our DM search analysis. In a next step, we fit these observations to $b\bar{b}$ and $\tau^+\tau^-$ annihilation/decay channels and several DM masses using a template-fitting approach. In the absence of a DM signal in all the considered DM scenarios, we find the following results:

* DM annihilation: We obtain 95% C.L. upper limits on the annihilation cross-section for different DM masses (figure 19). The best results are for the $\tau^+\tau^-$ annihilation channel, reaching values of $<\sigma v> \sim 5 \times 10^{-24} \text{ cm}^3 \text{ s}^{-1}$ for DM masses up to $\sim 1$ TeV. For masses above few TeV, our limits become less constraining as the sensitivity of the CTA weakens, reaching values up to $10^{-21} \text{ cm}^3 \text{ s}^{-1}$. These limits are between 2–4 orders of magnitude above the value of the thermal-relic cross-section. Different prescriptions for the subhalo population in Perseus weaken (MIN) or strengthen (MAX) our results by a factor $\sim \mathcal{O}(10)$ (figure 20).

* DM decay: We obtain 95% C.L. lower limits on the DM particle lifetime for different DM masses (figure 21). The best results are again for the $\tau^+\tau^-$ decay channel, reaching values of $\tau_\chi \sim 10^{27} \text{ s}$ for DM masses in the range of 10-30 TeV.

To put these DM results into the more general context: for WIMP annihilation, our prospects show that CTA will be able to provide the best constraints from γ-ray DM searches in galaxy clusters above 1 TeV (figure 19). Note, though, that CTA is not expected to reach
the thermal relic cross-section value from these observations.<fnref n="36" /> On the other hand, the sensitivity of CTA to DM decay in Perseus will be particularly significant. Indeed, comparing our predictions with other existing constraints, CTA will test an unexplored region of the DM decay parameter space for TeV WIMPs. This will allow CTA to set unprecedented constraints on the DM particle lifetime at these masses (figures 21 and 22). All these DM results not only demonstrate the superb capabilities of CTA to search for DM and to test the preferred DM models, especially for heavy WIMPs, but also the excellent potential of galaxy clusters as excellent targets for γ-ray DM searches.

## Acknowledgments

This work was conducted in the context of the CTA Consortium (mainly the CTA DMEP and CR Working Groups and the CTA Galaxy Clusters Task Force). We gratefully acknowledge financial support from the following agencies and organizations:

State Committee of Science of Armenia, Armenia; The Australian Research Council, Astronomy Australia Ltd, The University of Adelaide, Australian National University, Monash University, The University of New South Wales, The University of Sydney, Western Sydney University, Australia; Federal Ministry of Education, Science and Research, and Innsbruck University, Austria; Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq), Fundação de Amparo à Pesquisa do Estado do Rio de Janeiro (FAPERJ), Fundação de Amparo à Pesquisa do Estado de São Paulo (FAPESP), Fundação de Apoio à Ciência, Tecnologia e Inovação do Paraná — Fundação Araucária, Ministry of Science, Technology, Innovations and Communications (MCTIC), Brasil; Ministry of Education and Science, National RI Roadmap Project DO1-153/28.08.2018, Bulgaria; The Natural Sciences and Engineering Research Council of Canada and the Canadian Space Agency, Canada; CONICYT-Chile grants CATA AFB 170002, ANID PIA/APOYO AFB 180002, ACT 1406, FONDECYT-Chile grants, 1161463, 1170171, 1190886, 1171421, 1170345, 1201582, Gemini-ANID 32180007, Chile, W.M. gratefully acknowledges support by the ANID BASAL projects ACE210002 and FB210003, and FONDECYT 11190853; Croatian Science Foundation, Rudjer Boskovic Institute, University of Osijek, University of Rijeka, University of Split, Faculty of Electrical Engineering, Mechanical Engineering and Naval Architecture, University of Zagreb, Faculty of Electrical Engineering and Computing, Croatia; Ministry of Education, Youth and Sports, MEYS LM2015046, LM2018105, LTT17006, EU/MEYS CZ.02.1.01/0.0/0.0/16_013/0001403, CZ.02.1.01/0.0/0.0/18_046/0016007 and CZ.02.1.01/0.0/0.0/16_019/0000754, Czech Republic; Academy of Finland (grant nr.317636 and 320045), Finland; Ministry of Higher Education and Research, CNRS-INSU and CNRS-IN2P3, CEA-Irfu, ANR, Regional Council Ile de France, Labex ENIGMASS, OCEVU, OSUG2020 and P2IO, France; Max Planck Society, BMBF, DESY, Helmholtz Association, TU Dortmund University grant DFG SFB 1491, Germany; Department of Atomic Energy, Department of Science and Technology, India; Istituto Nazionale di Astrofisica (INAF), Istituto Nazionale di Fisica Nucleare (INFN), MIUR, Istituto Nazionale di Astrofisica (INAF-OABRERA) Grant Fondazione Cariplo/Regione Lombardia ID 2014-1980/RST_ERC, Italy; ICRR, University of Tokyo, JSPS, MEXT, Japan; Netherlands

<footnote n="36">³⁶Yet, we did not include possible flux enhancements due to the Sommerfeld effect [231, 232], that could drastically impact the constraints especially in the TeV range. This will be done elsewhere.</footnote>
Research School for Astronomy (NOVA), Netherlands Organization for Scientific Research (NWO), Netherlands; University of Oslo, Norway; Ministry of Science and Higher Education, DIR/WK/2017/12, the National Centre for Research and Development and the National Science Centre, UMO-2016/22/M/ST9/00583, Poland; Slovenian Research Agency, grants P1-0031, P1-0385, I0-0033, J1-9146, J1-1700, N1-0111, and the Young Researcher program, Slovenia; South African Department of Science and Technology and National Research Foundation through the South African Gamma-Ray Astronomy Programme, South Africa; The Spanish groups acknowledge the Spanish Ministry of Science and Innovation and the Spanish Research State Agency (AEI) through the government budget lines PGE2021/28.06.000X.411.01, PGE2022/28.06.000X.411.01 and PGE2022/28.06.000X.711.04, and grants PGC2018-095161-B-I00, PGC2018-095512-B-I00, PID2019-104114RB-C31, PID2019-107847RB-C44, PID2019-104114RB-C32, PID2019-105510GB-C31, PID2019-104114RB-C33, PID2019-107847RB-C41, PID2019-107847RB-C43, PID2019-107847RB-C42, PID2019-107988GB-C22; the "Centro de Excelencia Severo Ochoa" program through grants no. SEV-2017-0709, CEX2019-000920-S, CEX2020-001007-S; the "Unidad de Excelencia María de Maeztu" program through grants no. CEX2019-000918-M, CEX2020-001058-M; the "Ramón y Cajal" program through grant RYC-2017-22665; the "Juan de la Cierva-Incorporación" program through grants no. IJC2018-037195-I, IJC2019-040315-I. They also acknowledge the La Caixa Banking Foundation, grant no. LCF/BQ/PI21/11830030; the "Programa Operativo" FEDER 2014-2020, Consejería de Economía y Conocimiento de la Junta de Andalucía (Ref. 1257737), PAIDI 2020 (Ref. P18-FR-1580) and Universidad de Jaén; "Programa Operativo de Crecimiento Inteligente" FEDER 2014-2020 (Ref. ESFRI-2017-IAC-12), Ministerio de Ciencia e Innovación, 15% co-financed by Consejería de Economía, Industria, Comercio y Conocimiento del Gobierno de Canarias; the "CERCA" program of the Generalitat de Catalunya; and the European Union's "Horizon 2020" GA:824064 and NextGenerationEU. Swedish Research Council, Royal Physiographic Society of Lund, Royal Swedish Academy of Sciences, The Swedish National Infrastructure for Computing (SNIC) at Lunarc (Lund), Sweden; State Secretariat for Education, Research and Innovation (SERI) and Swiss National Science Foundation (SNSF), Switzerland; Durham University, Leverhulme Trust, Liverpool University, University of Leicester, University of Oxford, Royal Society, Science and Technology Facilities Council, UK; U.S. National Science Foundation, U.S. Department of Energy, Argonne National Laboratory, Barnard College, University of California, University of Chicago, Columbia University, Georgia Institute of Technology, Institute for Nuclear and Particle Astrophysics (INPAC-MRPI program), Iowa State University, the Smithsonian Institution, V.V.D. is funded by NSF grant AST-1911061, Washington University McDonnell Center for the Space Sciences, The University of Wisconsin and the Wisconsin Alumni Research Foundation, USA.

The research leading to these results has received funding from the European Union's Seventh Framework Programme (FP7/2007-2013) under grant agreements No 262053 and No 317446. This project is receiving funding from the European Union's Horizon 2020 research and innovation programs under agreement No 676134. This research has made use of the CTA instrument response functions provided by the CTA Consortium and Observatory, see http://www.cta-observatory.org/science/ctao-performance/ for more details.

JPR work was supported by grant SEV-2016-0597-17-2 funded by MCIN/AEI/10.13039/501100011033 and "ESF Investing in your future". JPR research was co-funded by the
Slovenian Research Agency, grants P1-0031 and J1-3011. MASC was also supported by the *Atracción de Talento* contracts no. 2016-T1/TIC-1542 and 2020-5A/TIC-19725 granted by the Comunidad de Madrid in Spain. The work of JPR and MASC was additionally supported by the grants PGC2018-095161-B-I00 and CEX2020-001007-S, both funded by MCIN/AEI/10.13039/501100011033 and by "ERDF A way of making Europe". MH acknowledges funding from the Max Planck Society, the University of Tokyo, and the ICRR Inter-University Research Program in the Fiscal Years 2021 and 2022. PTL is supported by the Swedish Research Council under contract 2019-05135. SHC work is supported by UNAM-PAPIIT IG101323.

This research made use of *gammapy*,<fnref n="37" /> a community-developed core Python package for TeV γ-ray astronomy [223]. This research made use of Astropy, a community-developed core Python package for Astronomy [233], in addition to NumPy [234], SciPy [235], Heapy [236] and Ipython [237]. Figures were generated using Matplotlib [238].

We acknowledge support from the CNRS/IN2P3 Computing Center (Lyon — France) for providing computing and data-processing resources needed for this work. We thank the support of the Hydra HPC cluster in Instituto de Física Teórica (IFT UAM-CSIC) for the computing time and resources.

## A Validation of the thermal model using Planck data

The thermal Sunyaev-Zel'dovich (tSZ) effect provides a direct measurement of the electron thermal pressure in galaxy clusters. In order to validate our thermal gas pressure model calibration, which is based on the indirect inference of the pressure from X-ray observations, we compare the tSZ prediction of our model to Planck data [225] obtained with MILCA [239].

We extract a MILCA sky patch centered on the Perseus cluster and compute its profile. Uncertainties are computed according to the standard deviation of the map in regions free from emission. We project our tSZ model on the same sky patch, convolve the map to the 10 arcmin Planck beam and compare it to the data in figures 23 (map) and 24 (profile). Radio sources (including NGC 1275) are masked within 30 arcmin.

Our model, as calibrated using X-ray information, provides an excellent match to the Planck Compton parameter map and profile, and thus to the cluster thermal pressure and thermal energy distribution. The model is validated up to $\gtrsim 2R_{500}$. In the core the bright radio source associated with NGC 1275 prohibits the validation below 30 arcmin, but this is where high quality X-ray data are available. We verify that extrapolation of the model by [153] leads to an overestimation of the tSZ signal by a factor of a few.

## B Impact of the magnetic field on the model calibration

In figure 25, we show how the choice of the magnetic field strength model affects the observable in terms of radio synchrotron emission and γ-ray emission in the case of the pure hadronic scenario. The model parameters were fixed to {$X_{CRp,500}$, $\alpha_{CRp}$, $\eta_{CRp}$} = {0.075, 2.5, 0.75}. Since the radio data are limited to the cluster core, the observable profile is essentially sensitive to the core magnetic field. While the hadronic emission does not depend on the

<footnote n="37">³⁷https://gammapy.org/.</footnote>
<figure><img src="image_26.png" /></figure>

**Figure 23.** Comparison between the Compton parameter map measured by *Planck* and that inferred from our thermal model. *Left panel:* raw data. *Middle:* model. *Right panel:* residual. The dashed circle indicates $\theta_{500}$. Radio sources appear as negative point sources on the Compton parameter map. The two white circles are masks for NGC 1275 and another radio source.

<figure><img src="image_27.png" /></figure>

**Figure 24.** Comparison between the Compton parameter profile measured by *Planck* and that inferred from our thermal model. We also show the result of the extrapolation of the model by [153].

magnetic field, the inverse-Compton signal is affected by the magnetic field via the synchrotron losses that modify the spectral and spatial distribution of the secondary CRe. However, matching the CRp model normalization so that a given magnetic field model would fit the radio data would also imply changing the normalisation of the hadronic and inverse-Compton model, accordingly.

We also note that in the CTA energy range, the radio data implies that the inverse-Compton signal is too small to be detected, and much smaller than the expected hadronic signal. The predicted level of emission would be similar in a pure leptonic model assuming continuous injection (as it is the case in the hadronic model, via secondary particle production). However, in the case of discontinuous particle injection, energy losses are expected to lead to a cutoff in the CR spectrum, which will drastically reduce the amount of inverse-Compton emission at CTA energies. The level of emission seen in figure 25 can thus be considered as an upper limit for the inverse-Compton emission, and it is therefore neglected in the present paper.

In figure 26, instead we fit for the CRp parameters {$X_{CRp,500}$, $\alpha_{CRp}$, and $\eta_{CRp}$} by matching the radio synchrotron emission to the data in the pure hadronic scenario (as in section 2.3.2). We show the constraints on the recovered pure hadronic parameter space for three extreme magnetic field models. As we can see, changing the magnetic field model
<figure><img src="image_28.png" /></figure>

**Figure 25.** Same as figure 3, but showing the dependence of the pure hadronic model observables on the magnetic field model (see section 2.1.2 and figure 1). The color code is the same for all panels. *Top left panel:* radio synchrotron spectrum. *Top right panel:* radio synchrotron profile. *Bottom left panel:* γ-ray spectrum. The hadronic component, which does not depend on the magnetic field, is shown in black and the inverse-Compton in color, for the different magnetic field models. The shaded pink area corresponds to the CTA energy range. *Bottom right panel:* γ-ray profile. The dashed lines correspond to the CTA PSF and θ$_{500}$, respectively.

essentially affects the normalization, and only minor variations are observed on the other
parameters. This also mainly leads to a change in the normalization of the γ-ray observables.

C Fermi-LAT sources around the Perseus cluster

No sources other than NGC 1275 and IC 310 has been detected within the CTA field of view in the direction of Perseus at very high γ-ray energies yet. Thus, we do not consider any other object in the modelling of the sky. Yet, the improved CTA sensitivity and increased field of view may allow us to detect other γ-ray emitters. To highlight such opportunities, we use the Fermi-LAT catalog, 4FGL-DR2 [215, 216], and select all sources detected within 5 deg from the Perseus reference center. These sources are listed in table 7 and their locations are reported in figure 6. NGC 1275 and IC 310 are included in the 4FGL source list. The spectra of all of these sources, at Fermi-LAT energies, are best described by power-laws, except for the pulsar PSR J0340+4130 and NGC 1275. For reference, and comparison to detected VHE sources NGC 1275 and IC 310, we extrapolate their fluxes up to the CTA energy range at 100 GeV. Most of the obtained values correspond to $\lesssim 10^{-11}$ cm$^{-2}$ TeV$^{-1}$ s$^{-1}$, which is smaller but comparable to the extrapolation of the fluxes of NGC 1275 and IC 310 at the same energy
<figure><img src="image_29.png" /><figcaption>**Figure 26.** Constraints on the parameter space of the pure hadronic model, for three extreme magnetic field models. See section 2.1.2 and figure 1.</figcaption></figure>

(about $25 \text{ and } 5 \times 10^{-11} \text{ cm}^{-2} \text{ TeV}^{-1} \text{ s}^{-1}$) according to equations (4.1) and (4.2). However, we stress that a cutoff in the spectra is possible, especially for distant sources affected by EBL absorption, and these numbers provide optimistic fluxes. Additionally, some of the sources may be variable so that these values only provide rough estimates. Even in the case of a detection, none of these sources is expected to affect the results presented in this paper because of the large angular separation from the cluster center.

## D Galactic foreground estimate

The Galactic foreground is estimated using the models developed by [219, 240]. We consider the *min* and *max* models as two different estimates. In figure 27, we present the spatial distribution of the diffuse foreground signal at 1 TeV, and the spectrum integrated in different apertures around the cluster center. We can observe that the signal presents a soft gradient over the considered region, plus a few clumps. The latter are not correlated with the cluster and none of them is located within $\theta_{500}$. Even in the case of a large aperture of about 1 deg ($\simeq \theta_{500}$), the integrated diffuse foreground is expected to be subdominant compared to the cluster CR induced signal within the same region (baseline CR model). When reducing the
<table><thead><tr><th>4FGL name</th><th>Alt. name</th><th>z<sup>†</sup></th><th>Offset (deg)</th><th>Type</th><th>Spectrum model</th><th>PL index</th><th>PL extrapolation at 100 GeV (cm<sup>-2</sup> TeV<sup>-1</sup> s<sup>-1</sup>)</th></tr></thead><tbody><tr><td>J0319.8+4130</td><td>NGC 1275*</td><td>0.01756</td><td>0.0</td><td>RDG</td><td>LP</td><td>—</td><td>—</td></tr><tr><td>J0316.8+4120</td><td>IC 310*</td><td>0.01894</td><td>0.6</td><td>RDG</td><td>PL</td><td>1.85</td><td>4.5 × 10<sup>-11</sup></td></tr><tr><td>J0315.5+4231</td><td>NVSS J031527+423249</td><td>N/A</td><td>1.3</td><td>BCU</td><td>PL</td><td>1.95</td><td>1.3 × 10<sup>-11</sup></td></tr><tr><td>J0312.9+4119</td><td>B3 0309+411B</td><td>0.134</td><td>1.3</td><td>RDG</td><td>PL</td><td>2.69</td><td>1.5 × 10<sup>-12</sup></td></tr><tr><td>J0311.6+4134</td><td>B3 0308+413</td><td>N/A</td><td>1.5</td><td>BCU</td><td>PL</td><td>1.95</td><td>1.2 × 10<sup>-11</sup></td></tr><tr><td>J0333.8+4007</td><td>B3 0330+399</td><td>N/A</td><td>3.0</td><td>BCU</td><td>PL</td><td>2.17</td><td>1.3 × 10<sup>-11</sup></td></tr><tr><td>J0334.3+3920</td><td>4C +39.12</td><td>0.02059</td><td>3.5</td><td>RDG</td><td>PL</td><td>1.81</td><td>4.0 × 10<sup>-11</sup></td></tr><tr><td>J0310.9+3815</td><td>B3 0307+380</td><td>0.816</td><td>3.7</td><td>FSRQ</td><td>PL</td><td>2.40</td><td>8.3 × 10<sup>-12</sup></td></tr><tr><td>J0340.3+4130</td><td>PSR J0340+4130</td><td>0</td><td>3.9</td><td>PSR</td><td>PLSEC</td><td>—</td><td>—</td></tr><tr><td>J0342.2+3858</td><td>GB6 J0342+3858</td><td>0.945</td><td>5.0</td><td>FSRQ</td><td>PL</td><td>2.26</td><td>1.2 × 10<sup>-11</sup></td></tr></tbody></table>

Notes. *Detected at very high γ-ray energies; † Taken from the NASA/IPAC extragalactic database; PL: power-law; LP: LogParabola; PLSEC: PLSuperExpCutoff; RDG: radio galaxies; BCU: blazar candidates of uncertain type; FSRQ: flat-spectrum radio quasars type of blazar; PSR: pulsar.

<figure><img src="image_30.png" /><figcaption>**Table 7.** 4FGL-DR2 sources located within 5 deg from the Perseus cluster center.</figcaption></figure>

**Figure 27.** Left panel: Surface brightness image of the Galactic foreground contribution at 1 TeV (min model). Right panel: spectral energy distribution of the Galactic foreground estimates and comparison to the baseline cluster model and the contribution from NGC 1275. Three different aperture radius are used, as indicated in the legend. The min model is given in red and the max model is given in pink. See [219, 240] for more details about the modelling.

aperture, the Galactic foreground becomes completely negligible compared to the cluster emission since the latter is much more compact.

## E Convergence of the CR parameter constraints versus the number of simulations

In section 5.3, we derived exclusion limits on the CR normalization as a function of the CR spectral and spatial distribution parameters. This was done by averaging the upper limits obtained over a given number of simulations in order to reduce the Poissonian noise.

Figures 28 shows the evolution of the constraint on the CR normalisation upper limit (mean and standard deviation), as a function of the number of simulations performed. We can see that averaging the results over 50 simulations is enough to obtain uncertainties
<figure><img src="image_31.png" /></figure>

**Figure 28.** Evolution of the upper limit on the parameters $X_{500} \equiv X_{\text{CRp}}(R_{500})$ (95% confidence level), as a function of the number of simulations performed and for the different model parameters ($\alpha_{\text{CRp}}$, $\eta_{\text{CRp}}$) presented in figure 11. Left: normalized mean upper limit. Right: normalized mean standard deviation. Note that the results obtained for the different models are correlated because the same underlying Poissonian realization is used for the different models.

on the mean and the standard deviation of the exclusion limit lower than 10%, which is negligible compared to other sources of uncertainties discussed in the paper. Therefore, we use a set of 50 simulations in the analysis.

## F Validation of the cosmic-ray analysis with gammapy

In this appendix, we validate our baseline cosmic-ray analysis, which relies on `ctools` [220], using the alternative software `gammapy` [223]. First, we compute the number count prediction given the sky model discussed in section 4 using the `gammapy` framework and a similar setup as the one used for `ctools`, but we use only one fixed pointing direction here. We verify that for this setup the two count cubes agree within a few percent, so that the residual difference is negligible compared to other sources of uncertainties. We then reproduce the results presented in section 5.3 and figure 11 with `gammapy`. In this case, the 95% upper limits are computed using the `iminuit` software [228]. To reduce the computing time, we free only the normalization of each sky component in the fit and compute how the upper limit changes when considering the full parameter space (normalization plus spectral parameters) from test cases. The corresponding correction, of about 50% depending on the model, are applied to the results. As shown in figure 29, `gammapy` and `ctools` give consistent results on the constraints obtained in the cosmic-ray parameter space within error bars. The residual differences are smaller than other sources of uncertainties presented in the paper (e.g., the cluster modeling) and are likely due to the fact that only the normalization parameters are fitted, or the slight differences in the setup. We conclude that using either `ctools` or `gammapy` for the cosmic-ray analysis will not affect the results presented in the paper.

## G Cosmic-ray parameter constraints in the pure hadronic model

In this appendix, we report the constraints in the parameter space obtained as discussed in section 5.4.2, in figure 30. In the case of this data realization, the input parameters are recovered within 68% confidence, or right at the limit, for all parameters.
<figure><img src="image_32.png" /><figcaption>**Figure 29.** Comparison of the exclusion limit on the parameter $X_{\text{CRp}}$ obtained with ctools and gammapy. The four panels give the results for the different values of $\eta_{\text{CRp}}$. The error bars represent the standard deviation of the upper limits and the data points give the means. A set of 50 simulations is used both for ctools and gammapy.</figcaption></figure>

## H Interplay between the DM components and the astrophysical $\gamma$-ray sources

In sections 6.2 and 6.3, we have discussed the results of the template-fitting analysis in terms of the DM parameters. Now we will investigate the results obtained for the rest of free parameters (equation (6.2)) and scrutinize the correlations between them.

Appendix K (table 9) shows the recovered mean values for the parameters regarding the astrophysical backgrounds (CRs, NGC 1275, IC 310, instrumental). In all of the cases we recover, within $1\sigma$, the input values that were used to produce the simulations in the first place, independently of the DM channel or mass that we fit. These results reinforce the idea that the DM component seems to be, at most, mildly correlated with the parameters describing the astrophysical sources in the area.

To properly quantify the correlation between the different parameters we compute the correlation matrix. For this, we use an MCMC method based on the `emcee` python package [173]. We follow the same methodology explained in section 5.4.2 for the astrophysical parameters ($A_{\text{CR}}$, $A_{\text{PS1}}$, $\alpha_{\text{PS1}}$, $A_{\text{PS2}}$, $\alpha_{\text{PS2}}$, $A_{\text{bkg}}$, $\alpha_{\text{bkg}}$). We simplify our analysis only including the CR normalization due large amount of computation time that it is needed. For the DM normalization ($A_{\chi}$), we use a flat prior ($A_{\chi} \ge 0$) to avoid sampling non-physical values. Since this approach is computationally very expensive, we only perform the fit for a representative DM mass of $m_{\chi} = 10$ TeV, the $\tau^{+}\tau^{-}$ channel and the MED annihilation model, and use a total of 30 simulated observations. The obtained correlation matrix is shown in figure 31. Indeed, the DM normalization is only mildly anticorrelated with the CR normalization, showing little
<figure><img src="image_33.png" /><figcaption>**Figure 30.** Example of posterior constraint on the full parameter space, as in figure 14, but in the case of the pure hadronic model. Contours provide the 68% and 95% confidence interval. The black star shows the input model parameters.</figcaption></figure>

correlation with the parameters describing NGC 1275. This implies that the results obtained for the limits on DM cross-section and DM lifetime should be robust. We caution that this correlation matrix may be dependent on the considered DM scenario (annihilation/decay), channel or $m_χ$, that were not explored in their totality given the constraints on the computation time. In figure 31, we can also appreciate that, with respect to the correlation matrix shown in figure 10, which did not include DM fluxes, the correlation of the CR normalization with the normalization of NGC 1275 is strengthened. Also, NGC 1275 parameters change their correlation sign. In the case of a detection of γ-rays from Perseus, these correlations should be investigated in depth and properly taken into account in the analysis.
<figure><img src="image_34.png" /><figcaption>**Figure 31.** Correlation matrix obtained for a fit to the MED annihilation model, a DM mass of $m_\chi = 10$ TeV and the $\tau^+\tau^-$ channel. The color scale ranges from -1 for parameters completely degenerated and anticorrelated (in red), to 1 for those degenerated and correlated (in blue). See text for discussion on the found correlations.</figcaption></figure>

# I CTA sensitivity to DM under the ON-OFF observational setup

The ON-OFF method is one of the most standard observational and analysis approaches adopted by the existing IACTs to obtain DM constraints (e.g., [115, 117, 121]). This choice is mainly motivated by their limited FoV and angular resolution, which makes it very difficult a (optimal) use of spatial information in large regions of the sky, as needed for a template-fitting analysis. In fact, a large number of IACT DM analyses for clusters even consider the DM-induced $\gamma$-ray fluxes as a point-like source in the target center containing the integrated DM flux of the whole object, this way neglecting its extension as a first approximation to the DM search.

For the sake of comparison, in this appendix we first compute the CTA sensitivity to DM in Perseus assuming the DM-induced emission as a point-like source. We will only include the instrumental background and will neglect the rest of $\gamma$-ray sources in the area. This methodology, although over-simplistic, will provide a first-order evaluation of the CTA sensitivity to DM in the object and is expected to yield the most stringent (unrealizable) constraints. Yet, as described in section 4.2, the Perseus cluster hosts a very bright AGN, NGC 1275, in its center, which in a real observation may be difficult to neglect. Thus, in a second step, and in order to perform a more refined, realistic ON-OFF analysis, we use the DM templates to account for the spatial extension of the DM emission and, as adopted in [117], we place a circular mask of 0.1 deg. radius in the center of Perseus to avoid AGN