# Constraining the cosmic-ray pressure in the inner Virgo Cluster using H.E.S.S. observations of M 87

F Aharonian, F. Ait Benkhali, C Arcaro, J Aschersleben, M Backes, V

Barbosa Martins, R Batzofin, Y Becherini, D Berge, K Bernlöhr, et al.

► To cite this version:

F Aharonian, F. Ait Benkhali, C Arcaro, J Aschersleben, M Backes, et al. Constraining the cosmic-ray pressure in the inner Virgo Cluster using H.E.S.S. observations of M 87. Astronomy and Astrophysics - A&A, 2023, 675, pp.A138. 10.1051/0004-6361/202346056. hal-04111745

HAL Id: hal-04111745

https://hal.science/hal-04111745

Submitted on 5 May 2024

**HAL** is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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
# Constraining the cosmic-ray pressure in the inner Virgo Cluster using H.E.S.S. observations of M87

H.E.S.S. Collaboration: F. Aharonian<sup>1,2,3</sup>, F. Ait Benkhali<sup>4</sup>, C. Arcaro<sup>7,* Ⓐ</sup>, J. Aschersleben<sup>5</sup>, M. Backes<sup>6,7</sup>, V. Barbosa Martins<sup>8,* Ⓑ</sup>, R. Batzofin<sup>9</sup>, Y. Becherini<sup>10,11</sup>, D. Berge<sup>8,12</sup>, K. Bernlöhr<sup>2</sup>, B. Bi<sup>13</sup>, M. Böttcher<sup>7</sup>, C. Boisson<sup>14</sup>, J. Bolmont<sup>15</sup>, J. Borowska<sup>12</sup>, F. Bradascio<sup>16</sup>, M. Breuhaus<sup>2</sup>, R. Brose<sup>1</sup>, F. Brun<sup>16</sup>, B. Bruno<sup>17</sup>, T. Bulik<sup>18</sup>, C. Burger-Scheidlin<sup>1</sup>, T. Bylund<sup>11</sup>, S. Caroff<sup>19</sup>, S. Casanova<sup>20</sup>, R. Cecil<sup>21</sup>, J. Celic<sup>17</sup>, M. Cerruti<sup>10</sup>, T. Chand<sup>7</sup>, S. Chandra<sup>7</sup>, A. Chen<sup>22</sup>, J. Chibueze<sup>7</sup>, O. Chibueze<sup>7</sup>, G. Cotter<sup>23</sup>, J. Damascene Mbarubucyeye<sup>8</sup>, A. Djannati-Ataï<sup>10</sup>, K. Egberts<sup>9</sup>, J.-P. Ernenwein<sup>24</sup>, G. Fichet de Clairfontaine<sup>14</sup>, M. Filipovic<sup>25</sup>, G. Fontaine<sup>26</sup>, M. Füßling<sup>8</sup>, S. Funk<sup>17</sup>, S. Gabici<sup>10</sup>, S. Ghafourizadeh<sup>4</sup>, G. Giavitto<sup>8</sup>, D. Glawion<sup>17</sup>, J. F. Glicenstein<sup>16</sup>, P. Goswami<sup>7</sup>, G. Grolleron<sup>15</sup>, M.-H. Grondin<sup>27</sup>, L. Haerer<sup>2</sup>, M. Haupt<sup>8</sup>, G. Hermann<sup>2</sup>, J. A. Hinton<sup>2</sup>, T. L. Holch<sup>8</sup>, D. Horns<sup>21</sup>, M. Jamrozy<sup>28</sup>, F. Jankowsky<sup>4</sup>, V. Joshi<sup>17</sup>, I. Jung-Richardt<sup>18</sup>, E. Kasai<sup>6</sup>, K. Katarzyński<sup>29</sup>, R. Khatoon<sup>7</sup>, B. Khélifi<sup>10</sup>, W. Kluźniak<sup>30</sup>, Nu. Komin<sup>22</sup>, K. Kosack<sup>16</sup>, D. Kostunin<sup>8</sup>, R. G. Lang<sup>17</sup>, S. Le Stum<sup>24</sup>, F. Leitl<sup>17</sup>, A. Lemière<sup>10</sup>, M. Lemoine-Goumard<sup>27</sup>, J.-P. Lenain<sup>15</sup>, F. Leuschner<sup>13</sup>, T. Lohse<sup>12</sup>, A. Luashvili<sup>14</sup>, I. Lypova<sup>4</sup>, J. Mackey<sup>1</sup>, D. Malyshev<sup>13</sup>, D. Malyshev<sup>17</sup>, V. Marandon<sup>2</sup>, P. Marchegiani<sup>22</sup>, A. Marcowith<sup>31</sup>, P. Marinos<sup>32</sup>, G. Martí-Devesa<sup>33</sup>, R. Marx<sup>4</sup>, M. Meyer<sup>21</sup>, A. Mitchell<sup>17</sup>, R. Moderski<sup>30</sup>, L. Mohrmann<sup>2</sup>, A. Montanari<sup>16</sup>, E. Moulin<sup>16</sup>, J. Muller<sup>26</sup>, K. Nakashima<sup>17</sup>, M. de Naurois<sup>26,* Ⓐ</sup>, J. Niemiec<sup>20</sup>, A. Priyana Noel<sup>28</sup>, P. O'Brien<sup>34</sup>, S. Ohm<sup>8,* Ⓑ</sup>, L. Olivera-Nieto<sup>2</sup>, E. de Ona Wilhelmi<sup>8</sup>, S. Panny<sup>33</sup>, M. Panter<sup>2</sup>, R. D. Parsons<sup>12</sup>, G. Peron<sup>10</sup>, S. Pita<sup>10</sup>, D. A. Prokhorov<sup>35</sup>, H. Prokoph<sup>8</sup>, G. Pühlhofer<sup>13</sup>, M. Punch<sup>10</sup>, A. Quirrenbach<sup>4</sup>, P. Reichherzer<sup>16</sup>, A. Reimer<sup>33</sup>, O. Reimer<sup>33</sup>, M. Renaud<sup>31</sup>, F. Rieger<sup>2</sup>, G. Rowell<sup>32</sup>, B. Rudak<sup>30</sup>, E. Ruiz-Velasco<sup>2</sup>, V. Sahakian<sup>36</sup>, H. Salzmann<sup>13</sup>, D. A. Sanchez<sup>19</sup>, A. Santangelo<sup>13</sup>, M. Sasaki<sup>17</sup>, J. Schäfer<sup>17</sup>, F. Schüssler<sup>16</sup>, U. Schwanke<sup>12</sup>, J. N. S. Shapopi<sup>6</sup>, H. Sol<sup>14</sup>, A. Specovius<sup>17</sup>, S. Spencer<sup>17</sup>, Ł. Stawarz<sup>28</sup>, R. Steenkamp<sup>6</sup>, S. Steinmassl<sup>2</sup>, C. Steppa<sup>9</sup>, I. Sushch<sup>7</sup>, H. Suzuki<sup>37</sup>, T. Takahashi<sup>38</sup>, T. Tanaka<sup>37</sup>, A. M. Taylor<sup>8</sup>, R. Terrier<sup>10</sup>, M. Tsrhou<sup>8</sup>, N. Tsuji<sup>39</sup>, Y. Uchiyama<sup>40</sup>, C. van Eldik<sup>17</sup>, B. van Soelen<sup>41</sup>, M. Vecchi<sup>5</sup>, J. Veh<sup>17</sup>, C. Venter<sup>7</sup>, J. Vink<sup>35</sup>, T. Wach<sup>17</sup>, S. J. Wagner<sup>4</sup>, R. White<sup>2</sup>, A. Wierzcholska<sup>20</sup>, Y. W. Wong<sup>17</sup>, M. Zacharias<sup>4,7</sup>, D. Zargaryan<sup>1</sup>, A. A. Zdziarski<sup>30</sup>, A. Zech<sup>14</sup>, S. Zouari<sup>10</sup>, and N. Żywucka<sup>7,* Ⓑ</sup>

*(Affiliations can be found after the references)*

Received 1 February 2023 / Accepted 15 May 2023

## ABSTRACT

The origin of the gamma-ray emission from M87 is currently a matter of debate. This work aims to localize the very high-energy (VHE; 100 GeV – 100 TeV) gamma-ray emission from M87 and probe a potential extended hadronic emission component in the inner Virgo Cluster. The search for a steady and extended gamma-ray signal around M87 can constrain the cosmic-ray energy density and the pressure exerted by the cosmic rays onto the intracluster medium and allow us to investigate the role of cosmic rays in the active galactic nucleus feedback as a heating mechanism in the Virgo Cluster. The High Energy Stereoscopic System (H.E.S.S.) telescopes are sensitive to VHE gamma rays and have been used to observe M87 since 2004. We utilized a Bayesian block analysis to identify M87 emission states with H.E.S.S. observations from 2004 to 2021, dividing them into low, intermediate, and high states. Because of the causality argument, an extended ($\ge$1 kpc) signal is allowed only in steady emission states. Hence, we fitted the morphology of the 120 h low-state data and find no significant gamma-ray extension. Therefore, we derive for the low state an upper limit of 58" (corresponding to $\approx$4.6 kpc) in the extension of a single-component morphological model described by a rotationally symmetric 2D Gaussian model at the 99.7% confidence level. Our results exclude the radio lobes ($\approx$30 kpc) as the principal component of the VHE gamma-ray emission from the low state of M87. The gamma-ray emission is compatible with a single emission region at the radio core of M87. These results, with the help of two multiple-component models, constrain the maximum cosmic-ray to thermal pressure ratio to $X_{CR,max} \le 0.32$ and the total energy in cosmic-ray protons to $U_{CR} \le 5 \times 10^{58}$ erg in the inner 20 kpc of the Virgo Cluster for an assumed cosmic-ray proton power-law distribution in momentum with spectral index $\alpha_p = 2.1$.

**Key words.** astroparticle physics – gamma rays: galaxies: clusters – galaxies: clusters: intracluster medium – radio continuum: galaxies

## 1. Introduction

Galaxy clusters are the largest gravitationally bound structures in the Universe. The Virgo Cluster is a massive cluster of galaxies centered around the radio galaxy M87, which is 16.5 Mpc (Blakeslee et al. 2009; Bird et al. 2010) from Earth, that extends up to $\approx 2.3^\circ$ from its center ($r_{500} \approx 662.6$ kpc; Simionescu et al. 2017). It is known as a cool core (CC) cluster, meaning its central region ($\le 0.01 r_{180}$; Peterson & Fabian 2006; Urban et al. 2011) is filled with a plasma that is colder and denser than the surrounding

<footnote>* Corresponding authors;
e-mail: contact.hess@hess-experiment.eu</footnote>
# Constraining the cosmic-ray pressure in the inner Virgo Cluster using H.E.S.S. observations of M 87

H.E.S.S. Collaboration: F. Aharonian<sup>1,2,3</sup>, F. Ait Benkhali<sup>4</sup>, C. Arcaro<sup>7,* Ⓐ</sup>, J. Aschersleben<sup>5</sup>, M. Backes<sup>6,7</sup>, V. Barbosa Martins<sup>8,* Ⓑ</sup>, R. Batzofin<sup>9</sup>, Y. Becherini<sup>10,11</sup>, D. Berge<sup>8,12</sup>, K. Bernlöhr<sup>2</sup>, B. Bi<sup>13</sup>, M. Böttcher<sup>7</sup>, C. Boisson<sup>14</sup>, J. Bolmont<sup>15</sup>, J. Borowska<sup>12</sup>, F. Bradascio<sup>16</sup>, M. Breuhaus<sup>2</sup>, R. Brose<sup>1</sup>, F. Brun<sup>16</sup>, B. Bruno<sup>17</sup>, T. Bulik<sup>18</sup>, C. Burger-Scheidlin<sup>1</sup>, T. Bylund<sup>11</sup>, S. Caroff<sup>19</sup>, S. Casanova<sup>20</sup>, R. Cecil<sup>21</sup>, J. Celic<sup>17</sup>, M. Cerruti<sup>10</sup>, T. Chand<sup>7</sup>, S. Chandra<sup>7</sup>, A. Chen<sup>22</sup>, J. Chibueze<sup>7</sup>, O. Chibueze<sup>7</sup>, G. Cotter<sup>23</sup>, J. Damascene Mbarubucyeye<sup>8</sup>, A. Djannati-Ataï<sup>10</sup>, K. Egberts<sup>9</sup>, J.-P. Ernenwein<sup>24</sup>, G. Fichet de Clairfontaine<sup>14</sup>, M. Filipovic<sup>25</sup>, G. Fontaine<sup>26</sup>, M. Füßling<sup>8</sup>, S. Funk<sup>17</sup>, S. Gabici<sup>10</sup>, S. Ghafourizadeh<sup>4</sup>, G. Giavitto<sup>8</sup>, D. Glawion<sup>17</sup>, J. F. Glicenstein<sup>16</sup>, P. Goswami<sup>7</sup>, G. Grolleron<sup>15</sup>, M.-H. Grondin<sup>27</sup>, L. Haerer<sup>2</sup>, M. Haupt<sup>8</sup>, G. Hermann<sup>2</sup>, J. A. Hinton<sup>2</sup>, T. L. Holch<sup>8</sup>, D. Horns<sup>21</sup>, M. Jamrozy<sup>28</sup>, F. Jankowsky<sup>4</sup>, V. Joshi<sup>17</sup>, I. Jung-Richardt<sup>18</sup>, E. Kasai<sup>6</sup>, K. Katarzyński<sup>29</sup>, R. Khatoon<sup>7</sup>, B. Khélifi<sup>10</sup>, W. Kluźniak<sup>30</sup>, Nu. Komin<sup>22</sup>, K. Kosack<sup>16</sup>, D. Kostunin<sup>8</sup>, R. G. Lang<sup>17</sup>, S. Le Stum<sup>24</sup>, F. Leitl<sup>17</sup>, A. Lemière<sup>10</sup>, M. Lemoine-Goumard<sup>27</sup>, J.-P. Lenain<sup>15</sup>, F. Leuschner<sup>13</sup>, T. Lohse<sup>12</sup>, A. Luashvili<sup>14</sup>, I. Lypova<sup>4</sup>, J. Mackey<sup>1</sup>, D. Malyshev<sup>13</sup>, D. Malyshev<sup>17</sup>, V. Marandon<sup>2</sup>, P. Marchegiani<sup>22</sup>, A. Marcowith<sup>31</sup>, P. Marinos<sup>32</sup>, G. Martí-Devesa<sup>33</sup>, R. Marx<sup>4</sup>, M. Meyer<sup>21</sup>, A. Mitchell<sup>17</sup>, R. Moderski<sup>30</sup>, L. Mohrmann<sup>2</sup>, A. Montanari<sup>16</sup>, E. Moulin<sup>16</sup>, J. Muller<sup>26</sup>, K. Nakashima<sup>17</sup>, M. de Naurois<sup>26,* Ⓐ</sup>, J. Niemiec<sup>20</sup>, A. Priyana Noel<sup>28</sup>, P. O'Brien<sup>34</sup>, S. Ohm<sup>8,* Ⓑ</sup>, L. Olivera-Nieto<sup>2</sup>, E. de Ona Wilhelmi<sup>8</sup>, S. Panny<sup>33</sup>, M. Panter<sup>2</sup>, R. D. Parsons<sup>12</sup>, G. Peron<sup>10</sup>, S. Pita<sup>10</sup>, D. A. Prokhorov<sup>35</sup>, H. Prokoph<sup>8</sup>, G. Pühlhofer<sup>13</sup>, M. Punch<sup>10</sup>, A. Quirrenbach<sup>4</sup>, P. Reichherzer<sup>16</sup>, A. Reimer<sup>33</sup>, O. Reimer<sup>33</sup>, M. Renaud<sup>31</sup>, F. Rieger<sup>2</sup>, G. Rowell<sup>32</sup>, B. Rudak<sup>30</sup>, E. Ruiz-Velasco<sup>2</sup>, V. Sahakian<sup>36</sup>, H. Salzmann<sup>13</sup>, D. A. Sanchez<sup>19</sup>, A. Santangelo<sup>13</sup>, M. Sasaki<sup>17</sup>, J. Schäfer<sup>17</sup>, F. Schüssler<sup>16</sup>, U. Schwanke<sup>12</sup>, J. N. S. Shapopi<sup>6</sup>, H. Sol<sup>14</sup>, A. Specovius<sup>17</sup>, S. Spencer<sup>17</sup>, Ł. Stawarz<sup>28</sup>, R. Steenkamp<sup>6</sup>, S. Steinmassl<sup>2</sup>, C. Steppa<sup>9</sup>, I. Sushch<sup>7</sup>, H. Suzuki<sup>37</sup>, T. Takahashi<sup>38</sup>, T. Tanaka<sup>37</sup>, A. M. Taylor<sup>8</sup>, R. Terrier<sup>10</sup>, M. Tsrhou<sup>8</sup>, N. Tsuji<sup>39</sup>, Y. Uchiyama<sup>40</sup>, C. van Eldik<sup>17</sup>, B. van Soelen<sup>41</sup>, M. Vecchi<sup>5</sup>, J. Veh<sup>17</sup>, C. Venter<sup>7</sup>, J. Vink<sup>35</sup>, T. Wach<sup>17</sup>, S. J. Wagner<sup>4</sup>, R. White<sup>2</sup>, A. Wierzcholska<sup>20</sup>, Y. W. Wong<sup>17</sup>, M. Zacharias<sup>4,7</sup>, D. Zargaryan<sup>1</sup>, A. A. Zdziarski<sup>30</sup>, A. Zech<sup>14</sup>, S. Zouari<sup>10</sup>, and N. Żywucka<sup>7,* Ⓑ</sup>

*(Affiliations can be found after the references)*

Received 1 February 2023 / Accepted 15 May 2023

## ABSTRACT

The origin of the gamma-ray emission from M 87 is currently a matter of debate. This work aims to localize the very high-energy (VHE; 100 GeV – 100 TeV) gamma-ray emission from M 87 and probe a potential extended hadronic emission component in the inner Virgo Cluster. The search for a steady and extended gamma-ray signal around M 87 can constrain the cosmic-ray energy density and the pressure exerted by the cosmic rays onto the intracluster medium and allow us to investigate the role of cosmic rays in the active galactic nucleus feedback as a heating mechanism in the Virgo Cluster. The High Energy Stereoscopic System (H.E.S.S.) telescopes are sensitive to VHE gamma rays and have been used to observe M 87 since 2004. We utilized a Bayesian block analysis to identify M 87 emission states with H.E.S.S. observations from 2004 to 2021, dividing them into low, intermediate, and high states. Because of the causality argument, an extended ($\ge$1 kpc) signal is allowed only in steady emission states. Hence, we fitted the morphology of the 120 h low-state data and find no significant gamma-ray extension. Therefore, we derive for the low state an upper limit of 58" (corresponding to $\approx$4.6 kpc) in the extension of a single-component morphological model described by a rotationally symmetric 2D Gaussian model at the 99.7% confidence level. Our results exclude the radio lobes ($\approx$30 kpc) as the principal component of the VHE gamma-ray emission from the low state of M 87. The gamma-ray emission is compatible with a single emission region at the radio core of M 87. These results, with the help of two multiple-component models, constrain the maximum cosmic-ray to thermal pressure ratio to $X_{CR,max} \le 0.32$ and the total energy in cosmic-ray protons to $U_{CR} \le 5 \times 10^{58}$ erg in the inner 20 kpc of the Virgo Cluster for an assumed cosmic-ray proton power-law distribution in momentum with spectral index $\alpha_p = 2.1$.

**Key words.** astroparticle physics – gamma rays: galaxies: clusters – galaxies: clusters: intracluster medium – radio continuum: galaxies

## 1. Introduction

Galaxy clusters are the largest gravitationally bound structures in the Universe. The Virgo Cluster is a massive cluster of galaxies centered around the radio galaxy M 87, which is 16.5 Mpc (Blakeslee et al. 2009; Bird et al. 2010) from Earth, that extends up to $\approx 2.3^\circ$ from its center ($r_{500} \approx 662.6$ kpc; Simionescu et al. 2017). It is known as a cool core (CC) cluster, meaning its central region ($\le 0.01 r_{180}$; Peterson & Fabian 2006; Urban et al. 2011) is filled with a plasma that is colder and denser than the surrounding

<footnote>* Corresponding authors;
e-mail: contact.hess@hess-experiment.eu</footnote>
gas. The cooling flow model has been proposed to explain the formation of CC clusters based on the inward flow of radiatively cooled material, which, in the absence of heating mechanisms, results in a mass deposition of up to ~1000 $M_{\odot}$ yr$^{-1}$. In contrast to this prediction, one to two orders of magnitude lower mass deposition rates are observed in CC clusters, which results in the rise of the cooling flow problem. The discrepancy points to the need for an additional heating mechanism to counterbalance the radiative cooling of the intracluster medium (ICM; Churazov et al. 2002; Jacob & Pfrommer 2017a,b).

Thermal conduction is not capable of solely heating the ICM, since the necessary ICM conductivity would exceed the theoretical maximum, the Spitzer conductivity (Zakamska & Narayan 2003). The dissipation of sound waves and turbulent motions have also been proposed as heating sources, though their contribution is likely not sufficient to account for the missing heating source (Ruszkowski et al. 2004; Zhuravleva et al. 2014). Furthermore, the mechanical heating by hot bubbles (Brüggen et al. 2002; Mathews et al. 2006) can also contribute to the ICM heating, though the energy available is limited due to the disruption of the bubbles (Pfrommer 2013).

In addition to the aforementioned feedback channels, the active galactic nucleus (AGN) feedback via cosmic rays (CRs; Churazov et al. 2002; Peterson & Fabian 2006; Guo & Oh 2008) accelerated in the central AGN region could provide the necessary heat to prevent the ICM from cooling down further. Jacob & Pfrommer (2017a,b) propose that the heating by the AGN feedback and the thermal conduction counterbalances the radiative cooling at any distance from the cluster center (steady-state system). The AGN accretes cooled gas and launches relativistic jets of particles. This process transfers energy to the surrounding gas and delays its radiative cooling. The cosmic ray electrons (CRE) injected by the AGN suffer severe energy losses via synchrotron and inverse Compton emission at GeV-TeV energies. In fact, 10 GeV electrons in ~10 µG magnetic fields, as in the M 87 radio lobes (De Gasperin et al. 2012), have a radiative lifetime of ~10 Myr (Enßlin et al. 2011). Cosmic ray protons and nuclei (CRp), also accelerated in the central AGN, accumulate and fill the cluster over cosmic timescales due to the persistent infall of cooled gas and the long radiative lifetimes of CRp. In reality, the radiative lifetime of CRp above 10 GeV in the ICM is at least 60 times longer than the lifetime of CRE at any energy in the same ICM (Fig. 2 from Enßlin et al. 2011). Hadronic interactions of the CRp with the ICM lead to the production of charged and neutral pions. While charged pions decay into electrons and positrons (hereafter secondary electrons), neutral pions decay into gamma rays. Secondary electrons are responsible for at least part of the extended radio halo emission seen around several galaxy clusters (Jacob & Pfrommer 2017b). In fact, cluster radio halos are the primary evidence for the existence of CRs in galaxy clusters (Enßlin et al. 2011). In dense environments (~$10^{-3}$–$10^{-1}$ cm$^{-3}$) filled with CRp, such as the central region of galaxy clusters, neutral pion decay is expected to produce a spatially extended and steady gamma-ray signal.

There have been many attempts to predict and observe diffuse GeV gamma-ray emission from galaxy clusters, for instance, Pinzke et al. (2011), Huber et al. (2013), Prokhorov & Churazov (2014), Ackermann et al. (2014, 2015). Xi et al. (2018) are the first to claim the detection of GeV gamma rays from the Coma Cluster, and Baghmanyan et al. (2022) are the first to claim that the gamma-ray signal from the Coma Cluster is extended at GeV energies. The extended GeV emission from the Coma Cluster is better fitted in a hadronic ICM scenario, although other models are also able to explain the data (Adam et al. 2021). In the

TeV regime, searches have been conducted for GeV gamma-ray emission from the Coma, Abell 496, Abell 85, and Perseus galaxy clusters (Aharonian et al. 2009a,b; Domainko et al. 2009; Galante 2009; Arlen et al. 2012; Aleksić et al. 2012) without success.

Detection of extended very high-energy (VHE; 100 GeV–100 TeV) gamma-ray emission from clusters would not only establish a new class of VHE gamma-ray emitter but also support the hypothesis of AGN feedback by streaming CRs in CC clusters by revealing the total energy accelerated in CRp. The presence of buoyant rising bubbles in M 87 indicates that AGN feedback plays a significant role in the Virgo Cluster (Churazov et al. 2001). Due to its active nucleus and its proximity to Earth, M 87 is the best candidate for searching for steady and extended VHE gamma-ray emission from a galaxy cluster. However, the highly variable gamma-ray emission from its AGN (Aharonian et al. 2006a; Acciari et al. 2009; Abramowski et al. 2012) dominates over a potential steady emission component. This poses a challenge when it comes to accessing an underlying steady component (Barbiellini et al. 2014). Nevertheless, only the cluster diffuse emission, that is, the emission that is of hadronic origin, should extend beyond the jet and fill the inner region of the cluster (Arlen et al. 2012). Therefore, we aim to reduce the contribution of the variable gamma-ray emission to probe the cluster diffuse emission through an extended gamma-ray signal in the low state of M 87.

Imaging atmospheric Cherenkov telescopes (IACTs) are sensitive to VHE gamma rays from astrophysical sources (Hinton & Hofmann 2009). These telescopes can detect extragalactic extended gamma-ray emission (Abdalla et al. 2020a) with the High Energy Stereoscopic System (H.E.S.S.). This study uses observations of M 87 with H.E.S.S. from 2004 to 2021 to probe the origin and size of the gamma-ray emission in the low state of M 87. We first tested for deviations from a single point source by fitting a point-like model and a rotationally symmetric 2D Gaussian model to the emission. Afterward, we compared the morphology results to known features from the radio and X-ray emissions. Finally, we considered two scenarios for a multiple-component emission to interpret the results in terms of the CR pressure in the inner Virgo Cluster.

This paper is organized as follows. In Sect. 2 we introduce the H.E.S.S. observations of M 87 and the analysis methods. In Sect. 3 we present the results regarding the origin and morphology of M 87 gamma-ray emission. In Sect. 4 we interpret the results in terms of the CRp pressure in the inner Virgo Cluster and the total energy in CRp. Finally, in Sect. 5 we summarize the results and consider future prospects.

## 2. Methods

### 2.1. H.E.S.S. experiment and data analysis

H.E.S.S. (Aharonian et al. 2006b) is an array of five IACTs located in the Khomas Highland in Namibia, 1800 m above sea level. The experiment has operated since 2003 with four 12 m telescopes (CT 1-4) and since 2012 with an additional 28 m telescope (CT 5), not utilized in this analysis. In 2017 a major upgrade of the cameras of the first four telescopes (CT 1-4) improved the read-out and the stability of the system (Ashton et al. 2020).

M 87 has been observed with H.E.S.S. since 2004 during several monitoring campaigns and in reaction to alerts of flaring activities. The data utilized in this work span from 2004 up to 2021 and were selected considering the following selection cuts:
gas. The cooling flow model has been proposed to explain the formation of CC clusters based on the inward flow of radiatively cooled material, which, in the absence of heating mechanisms, results in a mass deposition of up to ~1000 $M_{\odot}$ yr$^{-1}$. In contrast to this prediction, one to two orders of magnitude lower mass deposition rates are observed in CC clusters, which results in the rise of the cooling flow problem. The discrepancy points to the need for an additional heating mechanism to counterbalance the radiative cooling of the intracluster medium (ICM; Churazov et al. 2002; Jacob & Pfrommer 2017a,b).

Thermal conduction is not capable of solely heating the ICM, since the necessary ICM conductivity would exceed the theoretical maximum, the Spitzer conductivity (Zakamska & Narayan 2003). The dissipation of sound waves and turbulent motions have also been proposed as heating sources, though their contribution is likely not sufficient to account for the missing heating source (Ruszkowski et al. 2004; Zhuravleva et al. 2014). Furthermore, the mechanical heating by hot bubbles (Brüggen et al. 2002; Mathews et al. 2006) can also contribute to the ICM heating, though the energy available is limited due to the disruption of the bubbles (Pfrommer 2013).

In addition to the aforementioned feedback channels, the active galactic nucleus (AGN) feedback via cosmic rays (CRs; Churazov et al. 2002; Peterson & Fabian 2006; Guo & Oh 2008) accelerated in the central AGN region could provide the necessary heat to prevent the ICM from cooling down further. Jacob & Pfrommer (2017a,b) propose that the heating by the AGN feedback and the thermal conduction counterbalances the radiative cooling at any distance from the cluster center (steady-state system). The AGN accretes cooled gas and launches relativistic jets of particles. This process transfers energy to the surrounding gas and delays its radiative cooling. The cosmic ray electrons (CRE) injected by the AGN suffer severe energy losses via synchrotron and inverse Compton emission at GeV-TeV energies. In fact, 10 GeV electrons in ~10 µG magnetic fields, as in the M 87 radio lobes (De Gasperin et al. 2012), have a radiative lifetime of ~10 Myr (Enßlin et al. 2011). Cosmic ray protons and nuclei (CRp), also accelerated in the central AGN, accumulate and fill the cluster over cosmic timescales due to the persistent infall of cooled gas and the long radiative lifetimes of CRp. In reality, the radiative lifetime of CRp above 10 GeV in the ICM is at least 60 times longer than the lifetime of CRE at any energy in the same ICM (Fig. 2 from Enßlin et al. 2011). Hadronic interactions of the CRp with the ICM lead to the production of charged and neutral pions. While charged pions decay into electrons and positrons (hereafter secondary electrons), neutral pions decay into gamma rays. Secondary electrons are responsible for at least part of the extended radio halo emission seen around several galaxy clusters (Jacob & Pfrommer 2017b). In fact, cluster radio halos are the primary evidence for the existence of CRs in galaxy clusters (Enßlin et al. 2011). In dense environments (~$10^{-3}$–$10^{-1}$ cm$^{-3}$) filled with CRp, such as the central region of galaxy clusters, neutral pion decay is expected to produce a spatially extended and steady gamma-ray signal.

There have been many attempts to predict and observe dif-
fuse GeV gamma-ray emission from galaxy clusters, for instance,
Pinzke et al. (2011), Huber et al. (2013), Prokhorov & Churazov
(2014), Ackermann et al. (2014, 2015). Xi et al. (2018) are the
first to claim the detection of GeV gamma rays from the Coma
Cluster, and Baghmanyan et al. (2022) are the first to claim that
the gamma-ray signal from the Coma Cluster is extended at GeV
energies. The extended GeV emission from the Coma Cluster is
better fitted in a hadronic ICM scenario, although other mod-
els are also able to explain the data (Adam et al. 2021). In the

There have been many attempts to predict and observe dif-
fuse GeV gamma-ray emission from galaxy clusters, for instance,
Pinzke et al. (2011), Huber et al. (2013), Prokhorov & Churazov
(2014), Ackermann et al. (2014, 2015). Xi et al. (2018) are the
first to claim the detection of GeV gamma rays from the Coma
Cluster, and Baghmanyan et al. (2022) are the first to claim that
the gamma-ray signal from the Coma Cluster is extended at GeV
energies. The extended GeV emission from the Coma Cluster is
better fitted in a hadronic ICM scenario, although other mod-
els are also able to explain the data (Adam et al. 2021). In the

TeV regime, searches have been conducted for GeV gamma-
ray emission from the Coma, Abell 496, Abell 85, and Perseus
galaxy clusters (Aharonian et al. 2009a,b; Domainko et al. 2009;
Galante 2009; Arlen et al. 2012; Aleksić et al. 2012) without suc-
cess.

Detection of extended very high-energy (VHE; 100 GeV–100 TeV) gamma-ray emission from clusters would not only establish a new class of VHE gamma-ray emitter but also support the hypothesis of AGN feedback by streaming CRs in CC clusters by revealing the total energy accelerated in CRp. The presence of buoyant rising bubbles in M 87 indicates that AGN feedback plays a significant role in the Virgo Cluster (Churazov et al. 2001). Due to its active nucleus and its proximity to Earth, M 87 is the best candidate for searching for steady and extended VHE gamma-ray emission from a galaxy cluster. However, the highly variable gamma-ray emission from its AGN (Aharonian et al. 2006a; Acciari et al. 2009; Abramowski et al. 2012) dominates over a potential steady emission component. This poses a challenge when it comes to accessing an underlying steady component (Barbiellini et al. 2014). Nevertheless, only the cluster diffuse emission, that is, the emission that is of hadronic origin, should extend beyond the jet and fill the inner region of the cluster (Arlen et al. 2012). Therefore, we aim to reduce the contribution of the variable gamma-ray emission to probe the cluster diffuse emission through an extended gamma-ray signal in the low state of M 87.

Imaging atmospheric Cherenkov telescopes (IACTs) are sensitive to VHE gamma rays from astrophysical sources (Hinton & Hofmann 2009). These telescopes can detect extragalactic extended gamma-ray emission (Abdalla et al. 2020a) with the High Energy Stereoscopic System (H.E.S.S.). This study uses observations of M 87 with H.E.S.S. from 2004 to 2021 to probe the origin and size of the gamma-ray emission in the low state of M 87. We first tested for deviations from a single point source by fitting a point-like model and a rotationally symmetric 2D Gaussian model to the emission. Afterward, we compared the morphology results to known features from the radio and X-ray emissions. Finally, we considered two scenarios for a multiple-component emission to interpret the results in terms of the CR pressure in the inner Virgo Cluster.

This paper is organized as follows. In Sect. 2 we introduce the H.E.S.S. observations of M 87 and the analysis methods. In Sect. 3 we present the results regarding the origin and morphology of M 87 gamma-ray emission. In Sect. 4 we interpret the results in terms of the CRp pressure in the inner Virgo Cluster and the total energy in CRp. Finally, in Sect. 5 we summarize the results and consider future prospects.

## 2. Methods

2.1. H.E.S.S. experiment and data analysis

H.E.S.S. (Aharonian et al. 2006b) is an array of five IACTs
located in the Khomas Highland in Namibia, 1800 m above
sea level. The experiment has operated since 2003 with four
12 m telescopes (CT 1-4) and since 2012 with an additional
28 m telescope (CT 5), not utilized in this analysis. In 2017
a major upgrade of the cameras of the first four telescopes
(CT 1-4) improved the read-out and the stability of the system
(Ashton et al. 2020).

M 87 has been observed with H.E.S.S. since 2004 during
several monitoring campaigns and in reaction to alerts of flaring
activities. The data utilized in this work span from 2004 up to
2021 and were selected considering the following selection cuts:
at least three telescopes participating in the observation, a maximum zenith angle of observation of 50°, and an energy threshold of 300 GeV for all the reconstructed events. These selection cuts ensure that the later reconstructed gamma-ray events have improved angular resolution. Since the majority of the observations are from the first phase of the experiment (CT 1-4), CT5 data were not utilized in the analysis to avoid including different sources of systematic uncertainties. Finally, a total of 194 h of observations was obtained after the application of the quality criteria suitable for spectral and morphological analyses (Aharonian et al. 2006b).

We analyzed the data with the H.E.S.S. analysis software applying the template-based reconstruction technique ImPACT (Parsons & Hinton 2014). The gamma-hadron separation was performed through a multivariate analysis method (Ohm et al. 2009). The ring background and the reflected-region background (Berge et al. 2007) techniques provided an estimation of the remaining background. While the first method is ideal for studying the morphology of the emission, the second method is optimized for spectral studies. The point spread function (PSF), that is, the response of the instrument to a point-like source, was estimated based on Monte Carlo simulations that accounted for the observation conditions. With the template reconstruction, the PSF reaches a 68% containment radius of approximately 0.05° for energies above 1 TeV. We estimated the flux and spectrum of M 87 from a region within $\theta \le 0.071°$ around the radio core of M 87, that is, where the supermassive black hole (SMBH) is located. The exact position of the M 87 radio core is at right ascension (RA) 12h30m49.423s and declination (Dec) 12°23′28.04″ (Massaro et al. 2013). The region is defined to optimize the telescope's sensitivity toward a point-like source (Parsons & Hinton 2014). A forward-folding method (Piron et al. 2001) yielded the best power-law (PL) spectrum that fits the data above a safe energy threshold (Aharonian et al. 2006b).

## 2.2. The light curve and Bayesian blocks

Based on the gamma-ray flux derived for the individual observations, we binned the flux data points in time, weighting them according to their statistical uncertainties. The result is the long-term light curve of the source.

To identify the low flux state, we utilized a Bayesian blocks algorithm (Scargle 1998; Scargle et al. 2013) on the 30-day-binned light curve. We also investigated other bin sizes for the light curve; these are shown in Appendix A. The Bayesian block analysis is a method for detecting statistically significant changes in data from counting detectors and can be utilized to estimate flux levels assuming a piecewise-constant representation of time series data. It is effective also for non-constant sampling rates and hence is a very useful method for VHE gamma-ray astronomy (Ahnen et al. 2016). The method requires setting a prior function that directly influences the expected total number of blocks. We chose the prior function such that the false-positive rate, that is, the chance of wrongly detecting a flux change, is at the level of 5% (Scargle et al. 2013). Studies that directly rely on the data points to define the low state, that is, without the use of the Bayesian blocks, cannot assure a low false-positive rate, leading to an unreliable low-state data set.

To assign the Bayesian blocks to distinct source states, we inspected Fig. 1 and identified that Block 3 contains the 2008 VHE gamma-ray flare (Acciari et al. 2009). Furthermore, block 8 contains one data point more than $1\sigma$ above the average flux. Therefore, blocks 3 and 8 are not allowed in the low state, which was then defined as the blocks below the average flux. Blocks

<figure><img src="image_1.png" /><figcaption>**Fig. 1.** Thirty-day-binned H.E.S.S. light curve above 300 GeV with its Bayesian blocks and source states. The H.E.S.S. data points are represented by the black data points and a constant fit by the solid orange line. Dashed orange lines indicate a change in the flux level. The blue blocks represent the low state, the light gray the intermediate state, and the dark gray blocks the high state blocks. The orange numbers at the bottom indicate the labels of the blocks.</figcaption></figure>

3 and 8 were assigned to the intermediate state since they are slightly (<30%) higher than the average flux. The remaining blocks were assigned to the high state. Dividing the Bayesian blocks into three states helps us achieve a clear distinction between the low and high-state data sets, which minimizes the variable emission present in the low state. The definition of source states utilized in this study is focused on retrieving a low-state data set from the 30-day-binned light curve with a minimum contribution from a variable component and might not be representative for flaring studies.

The low state is our prime target for probing extended gamma-ray emission in M 87 data. In the low state, the contribution from the tails of the PSF distribution from a point-like component at the (variable) core of M 87 is minimized. Therefore, H.E.S.S. is most sensitive to detecting a diffuse emission from the inner Virgo Cluster in the low state. In addition to probing for extended emission in the low state of M 87, we also tested if the source states have compatible emission regions through the morphology fit (Sect. 2.3).

Since the choice of the light curve bin size also influences the defined source states, we also analyzed the data and fitted the morphology of a low state derived from a daily binned light curve. The results were rather poor due to the lack of event statistics in the low state (≈24 h lifetime of observations). This test and the results shown in Appendix A favor the choice of a bin size larger than 15 days. We decided on a light curve with a bin size of 30 days since it provides good event statistics in the low state for a morphology fit, that is, the statistical uncertainties are at the same level as the systematic uncertainties in the position determination of 20" (Acero et al. 2010). Therefore, by binning the flux points in time we intentionally allowed a certain contribution from variable emission into the low state in favor of larger event statistics. We later accounted for and estimated the contribution of a varying source component to the low state through the introduction of hybrid models (Sect. 4).

## 2.3. The morphology fit procedure

After dividing the data set into three different source states, we analyzed the individual states separately with the H.E.S.S. analysis software (Sect. 2.1). We obtained 2D distributions (sky maps)
at least three telescopes participating in the observation, a maximum zenith angle of observation of 50°, and an energy threshold of 300 GeV for all the reconstructed events. These selection cuts ensure that the later reconstructed gamma-ray events have improved angular resolution. Since the majority of the observations are from the first phase of the experiment (CT 1-4), CT5 data were not utilized in the analysis to avoid including different sources of systematic uncertainties. Finally, a total of 194 h of observations was obtained after the application of the quality criteria suitable for spectral and morphological analyses (Aharonian et al. 2006b).

We analyzed the data with the H.E.S.S. analysis software applying the template-based reconstruction technique ImPACT (Parsons & Hinton 2014). The gamma-hadron separation was performed through a multivariate analysis method (Ohm et al. 2009). The ring background and the reflected-region background (Berge et al. 2007) techniques provided an estimation of the remaining background. While the first method is ideal for studying the morphology of the emission, the second method is optimized for spectral studies. The point spread function (PSF), that is, the response of the instrument to a point-like source, was estimated based on Monte Carlo simulations that accounted for the observation conditions. With the template reconstruction, the PSF reaches a 68% containment radius of approximately 0.05° for energies above 1 TeV. We estimated the flux and spectrum of M 87 from a region within $\theta \le 0.071°$ around the radio core of M 87, that is, where the supermassive black hole (SMBH) is located. The exact position of the M 87 radio core is at right ascension (RA) 12h30m49.423s and declination (Dec) 12°23′28.04″ (Massaro et al. 2013). The region is defined to optimize the telescope's sensitivity toward a point-like source (Parsons & Hinton 2014). A forward-folding method (Piron et al. 2001) yielded the best power-law (PL) spectrum that fits the data above a safe energy threshold (Aharonian et al. 2006b).

## 2.2. The light curve and Bayesian blocks

Based on the gamma-ray flux derived for the individual observations, we binned the flux data points in time, weighting them according to their statistical uncertainties. The result is the long-term light curve of the source.

To identify the low flux state, we utilized a Bayesian blocks algorithm (Scargle 1998; Scargle et al. 2013) on the 30-day-binned light curve. We also investigated other bin sizes for the light curve; these are shown in Appendix A. The Bayesian block analysis is a method for detecting statistically significant changes in data from counting detectors and can be utilized to estimate flux levels assuming a piecewise-constant representation of time series data. It is effective also for non-constant sampling rates and hence is a very useful method for VHE gamma-ray astronomy (Ahnen et al. 2016). The method requires setting a prior function that directly influences the expected total number of blocks. We chose the prior function such that the false-positive rate, that is, the chance of wrongly detecting a flux change, is at the level of 5% (Scargle et al. 2013). Studies that directly rely on the data points to define the low state, that is, without the use of the Bayesian blocks, cannot assure a low false-positive rate, leading to an unreliable low-state data set.

To assign the Bayesian blocks to distinct source states, we inspected Fig. 1 and identified that Block 3 contains the 2008 VHE gamma-ray flare (Acciari et al. 2009). Furthermore, block 8 contains one data point more than $1\sigma$ above the average flux. Therefore, blocks 3 and 8 are not allowed in the low state, which was then defined as the blocks below the average flux. Blocks

<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Thirty-day-binned H.E.S.S. light curve above 300 GeV with its Bayesian blocks and source states. The H.E.S.S. data points are represented by the black data points and a constant fit by the solid orange line. Dashed orange lines indicate a change in the flux level. The blue blocks represent the low state, the light gray the intermediate state, and the dark gray blocks the high state blocks. The orange numbers at the bottom indicate the labels of the blocks.</figcaption></figure>

3 and 8 were assigned to the intermediate state since they are slightly (<30%) higher than the average flux. The remaining blocks were assigned to the high state. Dividing the Bayesian blocks into three states helps us achieve a clear distinction between the low and high-state data sets, which minimizes the variable emission present in the low state. The definition of source states utilized in this study is focused on retrieving a low-state data set from the 30-day-binned light curve with a minimum contribution from a variable component and might not be representative for flaring studies.

The low state is our prime target for probing extended gamma-ray emission in M 87 data. In the low state, the contribution from the tails of the PSF distribution from a point-like component at the (variable) core of M 87 is minimized. Therefore, H.E.S.S. is most sensitive to detecting a diffuse emission from the inner Virgo Cluster in the low state. In addition to probing for extended emission in the low state of M 87, we also tested if the source states have compatible emission regions through the morphology fit (Sect. 2.3).

Since the choice of the light curve bin size also influences the defined source states, we also analyzed the data and fitted the morphology of a low state derived from a daily binned light curve. The results were rather poor due to the lack of event statistics in the low state (≈24 h lifetime of observations). This test and the results shown in Appendix A favor the choice of a bin size larger than 15 days. We decided on a light curve with a bin size of 30 days since it provides good event statistics in the low state for a morphology fit, that is, the statistical uncertainties are at the same level as the systematic uncertainties in the position determination of 20" (Acero et al. 2010). Therefore, by binning the flux points in time we intentionally allowed a certain contribution from variable emission into the low state in favor of larger event statistics. We later accounted for and estimated the contribution of a varying source component to the low state through the introduction of hybrid models (Sect. 4).

## 2.3. The morphology fit procedure

After dividing the data set into three different source states, we analyzed the individual states separately with the H.E.S.S. analysis software (Sect. 2.1). We obtained 2D distributions (sky maps)
for the estimated background, the exposure time of the observations, the H.E.S.S. PSF, and the detected events. The PSF was convolved with the H.E.S.S. systematic uncertainty of 20" (Acero et al. 2010), although we tested the effect of the convolution on the final results and determined that it is minor. The sky maps, derived using a pixel size of 36" centered at the M87 radio core, were included in the spatial model of the emission.

We utilized predefined spatial models for the distribution of the intrinsic gamma-ray emission: a point-like model and a rotationally symmetric 2D Gaussian model given by

$$f(x, y) = A \cdot \exp\left(-\frac{(x - x_0)^2 + (y - y_0)^2}{2\sigma_G^2}\right), \quad (1)$$

where x and y are variables, $x_0$ and $y_0$ are the center of the function, A is the amplitude, and $\sigma_G$ is the extension (width). While the point-like model represents the emission limited to the core region, the centrally peaked Gaussian model allows us to probe the radiative cooling of the CRs during their transport from the AGN core toward the ICM. Alternatively, extended models such as a top-hat function could also be utilized, although they are less physically motivated than a Gaussian distribution for the case of CRp accelerated in the central AGN.

The region inside a circle of 0.5° radius from the radio core defines the region of interest in the M87 low-state data, that is, the region considered in the fit. Since counts are sampled from Poisson distributions, the fit algorithm is based on the Cash statistic (Cash 1979). Afterward, the model that results in the smallest test statistic (TS) is determined to be the one that best describes the data<fnref n="1" />.

A set of systematic checks confirms the stability of the results toward different configurations<fnref n="2" />. We analyzed the M87 low-state data several times under different circumstances: (1) with the maximum zenith angle of observation set at 45°; (2) with the energy threshold set at 0.7 TeV; (3) with a sky map bin size of 18" and (4) with a shift of 18" in the center of the sky map. The checks show that the results are stable against these changes. Furthermore, we reanalyzed the gamma-ray extension of the Crab Nebula (Abdalla et al. 2020b) with the same procedure described here and obtained compatible (within $1\sigma$) results.

A cross-check analysis with an independent analysis chain based on a semi-analytical shower model (de Naurois & Rolland 2009) also confirms the robustness of the results of this work from the H.E.S.S. data analysis (Sect. 2.1) up to the morphology fit (Sect. 2.3).

## 3. Results

In this section we present the results of the analysis of the H.E.S.S. observations of M87 following the procedure described in Sect. 2. First, we derive the light curve and source states in Sect. 3.1. Afterward, we present the results of the morphology fit of the low state in Sect. 3.2 and derive conclusions about its origin in Sect. 3.3.

<table><caption>**Table 1.** Bayesian blocks of the 30-day-binned H.E.S.S. light curve with its classification into low, intermediate, and high flux states.</caption><thead><tr><th>Block</th><th>State</th><th>Start date</th><th>End date</th></tr></thead><tbody><tr><td>0</td><td>Low</td><td>2004-02-16</td><td>2005-02-25</td></tr><tr><td>1</td><td>High</td><td>2005-02-25</td><td>2005-03-27</td></tr><tr><td>2</td><td>High</td><td>2005-03-27</td><td>2005-11-07</td></tr><tr><td>3</td><td>Intermediate</td><td>2005-11-07</td><td>2008-03-11</td></tr><tr><td>4</td><td>Low</td><td>2008-03-11</td><td>2009-03-21</td></tr><tr><td>5</td><td>High</td><td>2009-03-21</td><td>2012-10-31</td></tr><tr><td>6</td><td>Low</td><td>2012-10-31</td><td>2017-10-20</td></tr><tr><td>7</td><td>High</td><td>2017-10-20</td><td>2018-10-15</td></tr><tr><td>8</td><td>Intermediate</td><td>2018-10-15</td><td>2021-04-17</td></tr></tbody></table>

### 3.1. Bayesian blocks and source states

We applied the Bayesian blocks technique to the 30-day-binned light curve as discussed in Sect. 2.2. Figure 1 shows the M87 light curve with the derived Bayesian blocks and source states. The start and end dates of the blocks are given in Table 1 and the average gamma-ray flux above 300 GeV throughout the entire period is $\approx 1.8 \times 10^{-12} \text{ cm}^{-2} \text{ s}^{-1}$.

We stacked data from the individual blocks according to the defined source states and analyzed them separately. The results are summarized in Table 2. The low state has the largest data set ($\approx 120$ h), while the high state has the highest significance of detection ($\approx 19\sigma$). The spectral index of a PL energy distribution<fnref n="3" /> hardens from the low to the high state, consistent with previous results (Aharonian et al. 2006a; MAGIC Collaboration 2020).

### 3.2. The morphology of the low state of M87

We performed a morphology fit to the M87 low-state gamma-ray emission to probe its extension. No indication of a spatially extended emission component is found when comparing a point-like and an extended model as described in Sect. 2.3. An upper limit (UL) on the Gaussian $\sigma_G$ extension of 58" is derived at the 99.7% confidence level (c.l.). This translates to a physical extension limit of $\approx 4.6$ kpc, assuming a distance to M87 of $\approx 16.5$ Mpc (Blakeslee et al. 2009; Bird et al. 2010). This result improves by a factor of approximately two over the latest results (MAGIC Collaboration 2020). Previous measurements for the best-fit position (Aharonian et al. 2006a; Acciari et al. 2008) are in agreement with our results as shown in Fig. 2 (left).

### 3.3. The origin of the gamma-ray emission of M87

The results of the best-fit position of a point-like and Gaussian model are given in Table 3 for the low state. A slight shift of $\approx 25''$ in the best-fit position of the point-like model from the radio core is present. To investigate this apparent shift, we derived the $3\sigma$ uncertainty contours. Figure 2 (right) shows the best-fit position of the point-like model with its $3\sigma$ statistical uncertainty contours and the 99.7% c.l. extension UL of the low state. We also plot the Very Large Array (VLA) radio 21 cm emission in color scale, since it traces the energetic electrons in the inner radio cocoon ($\lesssim 2$ kpc). The shift between the best-fit position of the low state and the radio core is less than $3\sigma$ even without including the systematic uncertainty of 20". Therefore,

<footnote n="1">¹ We utilized the following Python packages throughout the analysis: NumPy 1.17.2 (Harris et al. 2020), SciPy 1.3.1 (Virtanen et al. 2020), Matplotlib 3.5.2 (Hunter 2007), Astropy 3.2.2 (Astropy Collaboration 2013, 2018), Gammapy 0.17 (Deil et al. 2017; Nigro et al. 2019), and Sherpa 4.12.0 (Freeman et al. 2001; Refsdal et al. 2009).</footnote>

<footnote n="2">² More details in the Ph.D. Thesis Barbosa Martins (2022).</footnote>

<footnote n="3">³ $d\phi/dE = \phi_0(E/\text{TeV})^{-\Gamma}$, where $\Gamma$ is the spectral index, E the gamma-ray energy and $\phi_0$ the normalization at 1 TeV.</footnote>
for the estimated background, the exposure time of the obser-
vations, the H.E.S.S. PSF, and the detected events. The PSF
was convolved with the H.E.S.S. systematic uncertainty of 20"
(Acero et al. 2010), although we tested the effect of the convo-
lution on the final results and determined that it is minor. The
sky maps, derived using a pixel size of 36" centered at the M 87
radio core, were included in the spatial model of the emission.

We utilized predefined spatial models for the distribution of
the intrinsic gamma-ray emission: a point-like model and a rota-
tionally symmetric 2D Gaussian model given by

$$
f(x, y) = A \cdot \exp \left( - \frac{(x - x_0)^2 + (y - y_0)^2}{2\sigma_G^2} \right), \quad (1)
$$

where x and y are variables, x₀ and y₀ are the center of the func-
tion, A is the amplitude, and σG is the extension (width). While
the point-like model represents the emission limited to the core
region, the centrally peaked Gaussian model allows us to probe
the radiative cooling of the CRs during their transport from the
AGN core toward the ICM. Alternatively, extended models such
as a top-hat function could also be utilized, although they are less
physically motivated than a Gaussian distribution for the case of
CRp accelerated in the central AGN.

The region inside a circle of 0.5° radius from the radio core
defines the region of interest in the M87 low-state data, that
is, the region considered in the fit. Since counts are sampled
from Poisson distributions, the fit algorithm is based on the Cash
statistic (Cash 1979). Afterward, the model that results in the
smallest test statistic (TS) is determined to be the one that best
describes the data<fnref n="1" />.

A set of systematic checks confirms the stability of the results
toward different configurations<fnref n="2" />. We analyzed the M87 low-state
data several times under different circumstances: (1) with the
maximum zenith angle of observation set at 45°; (2) with the
energy threshold set at 0.7 TeV; (3) with a sky map bin size of
18" and (4) with a shift of 18" in the center of the sky map. The
checks show that the results are stable against these changes. Fur-
thermore, we reanalyzed the gamma-ray extension of the Crab
Nebula (Abdalla et al. 2020b) with the same procedure described
here and obtained compatible (within 1 σ) results.

A cross-check analysis with an independent analysis chain
based on a semi-analytical shower model (de Naurois & Rolland
2009) also confirms the robustness of the results of this work
from the H.E.S.S. data analysis (Sect. 2.1) up to the morphology
fit (Sect. 2.3).

3. Results

In this section we present the results of the analysis of
the H.E.S.S. observations of M87 following the procedure
described in Sect. 2. First, we derive the light curve and source
states in Sect. 3.1. Afterward, we present the results of the mor-
phology fit of the low state in Sect. 3.2 and derive conclusions
about its origin in Sect. 3.3.

<table><caption>**Table 1.** Bayesian blocks of the 30-day-binned H.E.S.S. light curve with its classification into low, intermediate, and high flux states.</caption>
  <thead>
    <tr>
      <th>Block</th>
      <th>State</th>
      <th>Start date</th>
      <th>End date</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <td>0</td>
      <td>Low</td>
      <td>2004-02-16</td>
      <td>2005-02-25</td>
    </tr>
    <tr>
      <td>1</td>
      <td>High</td>
      <td>2005-02-25</td>
      <td>2005-03-27</td>
    </tr>
    <tr>
      <td>2</td>
      <td>High</td>
      <td>2005-03-27</td>
      <td>2005-11-07</td>
    </tr>
    <tr>
      <td>3</td>
      <td>Intermediate</td>
      <td>2005-11-07</td>
      <td>2008-03-11</td>
    </tr>
    <tr>
      <td>4</td>
      <td>Low</td>
      <td>2008-03-11</td>
      <td>2009-03-21</td>
    </tr>
    <tr>
      <td>5</td>
      <td>High</td>
      <td>2009-03-21</td>
      <td>2012-10-31</td>
    </tr>
    <tr>
      <td>6</td>
      <td>Low</td>
      <td>2012-10-31</td>
      <td>2017-10-20</td>
    </tr>
    <tr>
      <td>7</td>
      <td>High</td>
      <td>2017-10-20</td>
      <td>2018-10-15</td>
    </tr>
    <tr>
      <td>8</td>
      <td>Intermediate</td>
      <td>2018-10-15</td>
      <td>2021-04-17</td>
    </tr>
  </tbody>
</table>

3.1. Bayesian blocks and source states

We applied the Bayesian blocks technique to the 30-day-binned
light curve as discussed in Sect. 2.2. Figure 1 shows the M87
light curve with the derived Bayesian blocks and source states.
The start and end dates of the blocks are given in Table 1 and the
average gamma-ray flux above 300 GeV throughout the entire
period is ≈1.8 × 10⁻¹² cm⁻<fnref n="2" /> s⁻<fnref n="1" />.

We stacked data from the individual blocks according to the
defined source states and analyzed them separately. The results
are summarized in Table 2. The low state has the largest data
set (≈120 h), while the high state has the highest significance
of detection (≈19σ). The spectral index of a PL energy distri-
bution<fnref n="3" /> hardens from the low to the high state, consistent with
previous results (Aharonian et al. 2006a; MAGIC Collaboration
2020).

3.2. *The morphology of the low state of M87*

We performed a morphology fit to the M87 low-state gamma-
ray emission to probe its extension. No indication of a spa-
tially extended emission component is found when comparing
a point-like and an extended model as described in Sect. 2.3. An
upper limit (UL) on the Gaussian σG extension of 58" is derived
at the 99.7% confidence level (c.l.). This translates to a physi-
cal extension limit of ≈4.6 kpc, assuming a distance to M87 of
≈16.5 Mpc (Blakeslee et al. 2009; Bird et al. 2010). This result
improves by a factor of approximately two over the latest results
(MAGIC Collaboration 2020). Previous measurements for the
best-fit position (Aharonian et al. 2006a; Acciari et al. 2008) are
in agreement with our results as shown in Fig. 2 (left).

3.3. The origin of the gamma-ray emission of M87

The results of the best-fit position of a point-like and Gaus-
sian model are given in Table 3 for the low state. A slight shift
of ≈25" in the best-fit position of the point-like model from
the radio core is present. To investigate this apparent shift, we
derived the 3σ uncertainty contours. Figure 2 (right) shows the
best-fit position of the point-like model with its 3σ statistical
uncertainty contours and the 99.7% c.l. extension UL of the low
state. We also plot the Very Large Array (VLA) radio 21 cm
emission in color scale, since it traces the energetic electrons
in the inner radio cocoon (≤2 kpc). The shift between the best-fit
position of the low state and the radio core is less than 3σ even
without including the systematic uncertainty of 20". Therefore,

<footnote n="1">¹ We utilized the following Python packages throughout the analysis: NumPy 1.17.2 (Harris et al. 2020), SciPy 1.3.1 (Virtanen et al. 2020), Matplotlib 3.5.2 (Hunter 2007), Astropy 3.2.2 (Astropy Collaboration 2013, 2018), Gammapy 0.17 (Deil et al. 2017; Nigro et al. 2019), and Sherpa 4.12.0 (Freeman et al. 2001; Refsdal et al. 2009).</footnote>

<footnote n="2">² More details in the Ph.D. Thesis Barbosa Martins (2022).</footnote>

<footnote n="3">³ dφ/dE = φ₀(E/TeV)⁻⁴, where Γ is the spectral index, E the gamma-ray energy and φ₀ the normalization at 1 TeV.</footnote>
<table><caption>**Table 2.** Results of the H.E.S.S. analysis for the low, intermediate, and high flux states.</caption><thead><tr><th>State</th><th>Excess (counts)</th><th>Excess-to-background ratio</th><th>Significance (σ)</th><th>Lifetime (h)</th><th>Flux at 1 TeV (10<sup>-13</sup> TeV<sup>-1</sup> cm<sup>-2</sup> s<sup>-1</sup>)</th><th>Spectral index</th><th>Integrated flux (10<sup>-12</sup> cm<sup>-2</sup> s<sup>-1</sup>)</th></tr></thead><tbody><tr><td>Low</td><td>593.6</td><td>0.49</td><td>15.6</td><td>120.4</td><td>3.4±0.2</td><td>2.63±0.09</td><td>1.50±0.13</td></tr><tr><td>Intermediate</td><td>198.8</td><td>0.68</td><td>10.5</td><td>28.5</td><td>4.1±0.9</td><td>2.36±0.10</td><td>1.57±0.20</td></tr><tr><td>High</td><td>397.1</td><td>1.34</td><td>19.4</td><td>29.0</td><td>10.4±0.5</td><td>2.25±0.05</td><td>3.76±0.22</td></tr></tbody></table>

**Notes.** The integrated flux is given above 300 GeV.

<figure><img src="image_2.png" /><figcaption>**Fig. 2.** Results of the morphology fit in comparison to the shape of the radio emission and previous results. Left: H.E.S.S. low-state excess counts map with the derived extension UL plotted alongside previous ULs and the radio contour from the VLA 90 cm map (Thopson et al. 1980; Owen et al. 2000). The H.E.S.S. map is smoothed with a rotationally symmetric 2D Gaussian function with $\sigma_G = 18''$ for better visualization. The best-fit position is shown by the blue marker for the Gaussian model, with its $\sigma_G$ extension UL at the 99.7% c.l given by the dashed dark blue circle. The blue error bars in the position include the $1\sigma$ statistical uncertainty from the fit and 20'' systematic uncertainty (Acero et al. 2010). Previous results by H.E.S.S. (Aharonian et al. 2006a), VERITAS (Acciari et al. 2008), and MAGIC (MAGIC Collaboration 2020) are shown with the orange, light blue, and purple circles, respectively. The radio core of M 87 is marked by the black cross (Massaro et al. 2013). The dotted black square indicates the region depicted in the right figure. Right: The VLA 21 cm (Becker et al. 1995; The FIRST Project Team 1994) radio map of M 87 with the best-fit position for the point-like model for low, intermediate, and high flux states, shown with the blue, green, and pink triangles, respectively. The solid line contours represent $3\sigma$ statistical uncertainties for the respective source states. The dashed blue circle shows the resulting 99.7% c.l. of the Gaussian $\sigma_G$ extension UL in the low state. The best-fit positions of the source states are consistent with one another and with the M 87 radio core, and they agree with previous results.</figcaption></figure>

<table><caption>**Table 3.** Best-fit parameters of the morphology fit.</caption><thead><tr><th>Model</th><th>RA(°)</th><th>Dec (°)</th><th>σ<sub>G</sub>(")</th></tr></thead><tbody><tr><td>Point-like</td><td>187.707±0.005</td><td>12.398±0.003</td><td>—</td></tr><tr><td>Gaussian</td><td>187.711±0.002</td><td>12.397<sup>+0.003</sup><sub>-0.002</sub></td><td>12<sup>+15</sup><sub>-12</sub></td></tr></tbody></table>

**Notes.** The best-fit parameters of the low state of M 87 (see Table 1) considering the point-like and the rotationally symmetric 2D Gaussian models. $1\sigma$ statistical errors are given. The lower limit in the $\sigma_G$ does not reach the $1\sigma$ level within the valid interval ($\sigma_G > 0$).

the best-fit position of the point-like model is consistent with the radio core.

Our extension UL on the low state of M 87 excludes the radio lobes ($\approx 30$ kpc; black contours in Fig. 2 left) as the principal component of the low state of M 87 gamma-ray emission. Since the radio emission of the inner cocoon is still contained within the extension UL (Fig. 2 right), we conclude that the inner radio cocoon cannot be excluded as the principal component. Further observations of the low state of M 87 will improve the extension UL and probe the region inside the inner radio cocoon.

Since the sizes of the kiloparsec-scale jet ($\le 1$ kpc) and the X-ray knots are smaller than the H.E.S.S. extension UL, they could still contribute to (part of) M 87 gamma-ray low-state emission. Furthermore, the results from The EHT MWL Science Working Group (2021) show that the VHE emission during the low state of M 87 cannot originate from a single zone leptonic scenario in the very close vicinity of the SMBH ($\le 10r_g \approx 0.003$ pc, where $r_g$ is the SMBH gravitational radius). On the other hand, the addition of a hadronic emission component in the close vicinity of the SMBH can indeed explain the broadband spectral emission (Alfaro et al. 2022; Boughelilba et al. 2022; Xue et al. 2022). Figure 3 summarizes the scales of some known structures in M 87, highlighting the region excluded by this work.

Despite the lack of an extension detected in the M 87 low-state gamma-ray emission, we tested the hypothesis that the same emission region is responsible for the gamma rays in the different source states. Hence, we also fitted the morphology of the intermediate and high-flux states. None of the source states show significant extended emission and the best-fit positions of the point-like model of the M 87 source states are consistent with each other.
<table><caption>**Table 2.** Results of the H.E.S.S. analysis for the low, intermediate, and high flux states.</caption><thead><tr><th>State</th><th>Excess (counts)</th><th>Excess-to-background ratio</th><th>Significance (σ)</th><th>Lifetime (h)</th><th>Flux at 1 TeV (10<sup>-13</sup> TeV<sup>-1</sup> cm<sup>-2</sup> s<sup>-1</sup>)</th><th>Spectral index</th><th>Integrated flux (10<sup>-12</sup> cm<sup>-2</sup> s<sup>-1</sup>)</th></tr></thead><tbody><tr><td>Low</td><td>593.6</td><td>0.49</td><td>15.6</td><td>120.4</td><td>3.4±0.2</td><td>2.63±0.09</td><td>1.50±0.13</td></tr><tr><td>Intermediate</td><td>198.8</td><td>0.68</td><td>10.5</td><td>28.5</td><td>4.1±0.9</td><td>2.36±0.10</td><td>1.57±0.20</td></tr><tr><td>High</td><td>397.1</td><td>1.34</td><td>19.4</td><td>29.0</td><td>10.4±0.5</td><td>2.25±0.05</td><td>3.76±0.22</td></tr></tbody></table>

**Notes.** The integrated flux is given above 300 GeV.

<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Results of the morphology fit in comparison to the shape of the radio emission and previous results. Left: H.E.S.S. low-state excess counts map with the derived extension UL plotted alongside previous ULs and the radio contour from the VLA 90 cm map (Thopson et al. 1980; Owen et al. 2000). The H.E.S.S. map is smoothed with a rotationally symmetric 2D Gaussian function with $\sigma_G = 18''$ for better visualization. The best-fit position is shown by the blue marker for the Gaussian model, with its $\sigma_G$ extension UL at the 99.7% c.l given by the dashed dark blue circle. The blue error bars in the position include the $1\sigma$ statistical uncertainty from the fit and 20'' systematic uncertainty (Acero et al. 2010). Previous results by H.E.S.S. (Aharonian et al. 2006a), VERITAS (Acciari et al. 2008), and MAGIC (MAGIC Collaboration 2020) are shown with the orange, light blue, and purple circles, respectively. The radio core of M 87 is marked by the black cross (Massaro et al. 2013). The dotted black square indicates the region depicted in the right figure. Right: The VLA 21 cm (Becker et al. 1995; The FIRST Project Team 1994) radio map of M 87 with the best-fit position for the point-like model for low, intermediate, and high flux states, shown with the blue, green, and pink triangles, respectively. The solid line contours represent $3\sigma$ statistical uncertainties for the respective source states. The dashed blue circle shows the resulting 99.7% c.l. of the Gaussian $\sigma_G$ extension UL in the low state. The best-fit positions of the source states are consistent with one another and with the M 87 radio core, and they agree with previous results.</figcaption></figure>

<table><caption>**Table 3.** Best-fit parameters of the morphology fit.</caption><thead><tr><th>Model</th><th>RA(°)</th><th>Dec (°)</th><th>σ<sub>G</sub>(")</th></tr></thead><tbody><tr><td>Point-like</td><td>187.707±0.005</td><td>12.398±0.003</td><td>—</td></tr><tr><td>Gaussian</td><td>187.711±0.002</td><td>12.397<sup>+0.003</sup><sub>-0.002</sub></td><td>12<sup>+15</sup><sub>-12</sub></td></tr></tbody></table>

**Notes.** The best-fit parameters of the low state of M 87 (see Table 1) considering the point-like and the rotationally symmetric 2D Gaussian models. $1\sigma$ statistical errors are given. The lower limit in the $\sigma_G$ does not reach the $1\sigma$ level within the valid interval ($\sigma_G > 0$).

the best-fit position of the point-like model is consistent with the radio core.

Our extension UL on the low state of M 87 excludes the radio lobes ($\approx 30$ kpc; black contours in Fig. 2 left) as the principal component of the low state of M 87 gamma-ray emission. Since the radio emission of the inner cocoon is still contained within the extension UL (Fig. 2 right), we conclude that the inner radio cocoon cannot be excluded as the principal component. Further observations of the low state of M 87 will improve the extension UL and probe the region inside the inner radio cocoon.

Since the sizes of the kiloparsec-scale jet ($\le 1$ kpc) and the X-ray knots are smaller than the H.E.S.S. extension UL, they could still contribute to (part of) M 87 gamma-ray low-state emission. Furthermore, the results from The EHT MWL Science Working Group (2021) show that the VHE emission during the low state of M 87 cannot originate from a single zone leptonic scenario in the very close vicinity of the SMBH ($\le 10r_g \approx 0.003$ pc, where $r_g$ is the SMBH gravitational radius). On the other hand, the addition of a hadronic emission component in the close vicinity of the SMBH can indeed explain the broadband spectral emission (Alfaro et al. 2022; Boughelilba et al. 2022; Xue et al. 2022). Figure 3 summarizes the scales of some known structures in M 87, highlighting the region excluded by this work.

Despite the lack of an extension detected in the M 87 low-state gamma-ray emission, we tested the hypothesis that the same emission region is responsible for the gamma rays in the different source states. Hence, we also fitted the morphology of the intermediate and high-flux states. None of the source states show significant extended emission and the best-fit positions of the point-like model of the M 87 source states are consistent with each other.
<figure><img src="image_3.png" /><figcaption>**Fig. 3.** Physical scales and prominent (jet) components in M 87 on a logarithmic scale with the region excluded by this work shaded in gray. This work (Sect. 3.2) excludes the region $\ge$4.6 kpc as the origin of the VHE gamma-ray emission of the low state of M 87. The half-light radius of the stars, $R_{1/2} \approx 7.2$ kpc (Weil et al. 1997), and the dimensions of the radio lobes are also shown. The order of magnitude of the size of the EHT ring-like structure ($\sim 10^{-3}$ pc) is shown (The Event Horizon Telescope Collaboration 2019; The EHT MWL Science Working Group 2021). The measured extension from Centaurus A is shown for comparison (Abdalla et al. 2020a).</figcaption></figure>

# 4. Discussion

Based on the extension UL derived in the previous section, we explore two plausible multiple-component scenarios that could lead to extended gamma-ray emission in the inner Virgo Cluster. Physical properties such as the energy content in CRp and the pressure they exert on the ICM are constrained and the role CRp play in feedback processes is discussed.

## 4.1. Physical scenarios for an extended gamma-ray emission

The M 87 radio lobes extend up to $\approx 30$ kpc from the core and trace mildly relativistic electrons. The contribution from secondary electrons, that is, from the decay of charged pions, to M 87 radio-lobe emission is subdominant (Pfrommer 2013). Reacceleration of electrons by plasma waves (Enßlin et al. 2011) has been proposed to boost electrons with energies between 0.1 and 10 GeV to tens of GeV, which would complement the radio-lobe emission. In fact, CRp would also be boosted by this reacceleration mechanism. Given their longer radiative lifetime in the ICM, the CRp dominate the pressure in the inner Virgo Cluster over CRe. In general, reacceleration through plasma waves is a rather inefficient process and it is unlikely to boost electrons up to TeV energies beyond the kiloparsec-jet region. Therefore, we can ignore the VHE gamma-ray contribution from reaccelerated mildly relativistic electrons in the cluster region.

On the other hand, the streaming CRp might hadronically interact with the local target material and produce pions. The neutral pions could generate a detectable and extended VHE gamma-ray signal within the cluster region (Sect. 1). However, for this to occur, a dense target material and a strong CRp component must be present. The depletion of target material in the X-ray cavities (Abdulla et al. 2019) could locally impede pion production, and subsequently reduce the total gamma-ray signal. While gamma-ray production through neutral pion decay is a promising scenario, it may face challenges due to the possible presence of material-depleted regions.

Particle acceleration could also take place in the lobes of radio galaxies as evidence from Fermi Large Area Telescope (LAT) observations of Fornax A indicates (Ackermann et al. 2016). The upscattering of the synchrotron and cosmic microwave background photons in the radio lobes by the locally accelerated CRe could also contribute to the VHE gamma-ray emission of M 87 low state. As the H.E.S.S. extension UL derived here is significantly smaller than the size of the radio lobes, this process is very likely not the principal contributor to the low-state gamma-ray emission. On the other hand, electrons accelerated in the central AGN could significantly contribute to the gamma-ray signal through inverse Compton scattering in the photon fields of the core region. The cooling time of 10 GeV CRe is typically $\sim 10$ Myr in a $\sim 10$ µG magnetic field and hence orders of magnitude shorter than the $p-p$ cooling time (see Fig. 2 in Enßlin et al. 2011). Therefore, primary CRe are expected to mainly contribute to the emission in the close vicinity of the SMBH, which H.E.S.S. cannot resolve. CRe accelerated in the central AGN lose most of their energy before reaching kiloparsec distances, but could be reaccelerated via stochastic and/or shear particle acceleration to TeV energies (Rieger et al. 2007). These reaccelerated electrons would scatter via the inverse Compton process in the photon fields of the jet and produce an extended gamma-ray jet emission (Abdalla et al. 2020a). The jet extension at radio to X-ray wavelengths is $\approx 1$ kpc ($\approx 13''$), and also cannot be resolved by H.E.S.S. Therefore, a potential leptonic gamma-ray emission from the M 87 kiloparsec jet will appear point-like and indistinguishable from the emission from the core in this work.

Gamma rays can also interact with extragalactic background light photons on their way to Earth and produce electron-positrons pairs. These will initiate electromagnetic cascades, potentially producing an extended gamma-ray halo signal (Anderhub et al. 2010). Given the close distance of M 87 to Earth ($\approx 16.5$ Mpc; Blakeslee et al. 2009; Bird et al. 2010), effects from the background light would become relevant to the
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Physical scales and prominent (jet) components in M 87 on a logarithmic scale with the region excluded by this work shaded in gray. This work (Sect. 3.2) excludes the region $\ge$4.6 kpc as the origin of the VHE gamma-ray emission of the low state of M 87. The half-light radius of the stars, $R_{1/2} \approx 7.2$ kpc (Weil et al. 1997), and the dimensions of the radio lobes are also shown. The order of magnitude of the size of the EHT ring-like structure ($\sim 10^{-3}$ pc) is shown (The Event Horizon Telescope Collaboration 2019; The EHT MWL Science Working Group 2021). The measured extension from Centaurus A is shown for comparison (Abdalla et al. 2020a).</figcaption></figure>

# 4. Discussion

Based on the extension UL derived in the previous section, we explore two plausible multiple-component scenarios that could lead to extended gamma-ray emission in the inner Virgo Cluster. Physical properties such as the energy content in CRp and the pressure they exert on the ICM are constrained and the role CRp play in feedback processes is discussed.

## 4.1. Physical scenarios for an extended gamma-ray emission

The M 87 radio lobes extend up to $\approx 30$ kpc from the core and trace mildly relativistic electrons. The contribution from secondary electrons, that is, from the decay of charged pions, to M 87 radio-lobe emission is subdominant (Pfrommer 2013). Reacceleration of electrons by plasma waves (Enßlin et al. 2011) has been proposed to boost electrons with energies between 0.1 and 10 GeV to tens of GeV, which would complement the radio-lobe emission. In fact, CRp would also be boosted by this reacceleration mechanism. Given their longer radiative lifetime in the ICM, the CRp dominate the pressure in the inner Virgo Cluster over CRe. In general, reacceleration through plasma waves is a rather inefficient process and it is unlikely to boost electrons up to TeV energies beyond the kiloparsec-jet region. Therefore, we can ignore the VHE gamma-ray contribution from reaccelerated mildly relativistic electrons in the cluster region.

On the other hand, the streaming CRp might hadronically interact with the local target material and produce pions. The neutral pions could generate a detectable and extended VHE gamma-ray signal within the cluster region (Sect. 1). However, for this to occur, a dense target material and a strong CRp component must be present. The depletion of target material in the X-ray cavities (Abdulla et al. 2019) could locally impede pion production, and subsequently reduce the total gamma-ray signal. While gamma-ray production through neutral pion decay is a promising scenario, it may face challenges due to the possible presence of material-depleted regions.

Particle acceleration could also take place in the lobes of radio galaxies as evidence from Fermi Large Area Telescope (LAT) observations of Fornax A indicates (Ackermann et al. 2016). The upscattering of the synchrotron and cosmic microwave background photons in the radio lobes by the locally accelerated CRe could also contribute to the VHE gamma-ray emission of M 87 low state. As the H.E.S.S. extension UL derived here is significantly smaller than the size of the radio lobes, this process is very likely not the principal contributor to the low-state gamma-ray emission. On the other hand, electrons accelerated in the central AGN could significantly contribute to the gamma-ray signal through inverse Compton scattering in the photon fields of the core region. The cooling time of 10 GeV CRe is typically $\sim 10$ Myr in a $\sim 10$ µG magnetic field and hence orders of magnitude shorter than the $p-p$ cooling time (see Fig. 2 in Enßlin et al. 2011). Therefore, primary CRe are expected to mainly contribute to the emission in the close vicinity of the SMBH, which H.E.S.S. cannot resolve. CRe accelerated in the central AGN lose most of their energy before reaching kiloparsec distances, but could be reaccelerated via stochastic and/or shear particle acceleration to TeV energies (Rieger et al. 2007). These reaccelerated electrons would scatter via the inverse Compton process in the photon fields of the jet and produce an extended gamma-ray jet emission (Abdalla et al. 2020a). The jet extension at radio to X-ray wavelengths is $\approx 1$ kpc ($\approx 13''$), and also cannot be resolved by H.E.S.S. Therefore, a potential leptonic gamma-ray emission from the M 87 kiloparsec jet will appear point-like and indistinguishable from the emission from the core in this work.

Gamma rays can also interact with extragalactic background light photons on their way to Earth and produce electron-positrons pairs. These will initiate electromagnetic cascades, potentially producing an extended gamma-ray halo signal (Anderhub et al. 2010). Given the close distance of M 87 to Earth ($\approx 16.5$ Mpc; Blakeslee et al. 2009; Bird et al. 2010), effects from the background light would become relevant to the
gamma-ray spectrum at energies $\ge$10 TeV (Franceschini et al. 2019). In fact, only 9 events above the background level are reconstructed with energies above 10 TeV in the low state, 4 in the intermediate state, and 13 in the high state. Therefore, the contribution of such an extended gamma-ray component in the low state of M87 can be safely neglected for morphological studies.

The CRp from the jet could also interact with the local photon fields and produce neutral pions, which could lead to an extended gamma-ray signal. However, Boughelilba et al. (2022) have shown that the accretion flow and the disk around the SMBH can be neglected as targets for particle-photon interactions in the jet.

Another strong candidate for accelerating CRs and producing an extended gamma-ray signal is cluster mergers. Giant radio halos ($\ge$200 kpc) and relics (Enßlin et al. 2011) are typically found in clusters with recent merging activities. Nevertheless, the Virgo Cluster does not have any of the aforementioned signals, although the cluster is also not yet completely dynamically relaxed (Kashibadze et al. 2020). During merger events, CRp are expected to be accelerated at shocks near the cluster's virial radius ($\approx$1.7 Mpc, that is, $\approx$6° for the Virgo Cluster; Ackermann et al. 2015), and they may be transported toward the cluster center, depending on the cluster's turbulent history. In the absence of turbulent advective transport that counterbalances diffusion, the CRp profile in galaxy clusters tends to flatten (Enßlin et al. 2011). Hadronically interacting with the local ICM, these CRp can generate a diffuse gamma-ray signal up to the cluster's virial radius. However, it would be difficult for H.E.S.S. to detect such an extended gamma-ray signal, as its sensitivity degrades significantly to almost 80% of the Crab flux for a source diameter of $\approx$2° (Casanova & Dingus 2008; Mitchell et al. 2023).

In contrast to merger events, accretion events, that is, merging smaller virialized objects, might produce stronger shocks and accelerate CRs more efficiently. Inoue et al. (2005) have shown that for a ~µG magnetic field in the shock region, CRp can be accelerated up to ~$10^{18}$ eV. These CRp would produce pairs of electron-positron, which would rapidly cool through inverse Compton scattering and synchrotron losses, ultimately generating a signal in the form of gamma rays with an energy flux of $10^{-12}-10^{-11}$ erg cm$^{-2}$ s$^{-1}$, marginally detectable by H.E.S.S given the extension of the signal. For the Virgo Cluster the shock diameter can reach up to $\approx$1.9° around the virial radius ($\approx$6°), that is, outside the region of interest of our analysis.

Finally, the annihilation of hypothetical weakly interacting massive particles in the dark matter halo around M87 could also contribute to an extended gamma-ray signal, as investigated in Ackermann et al. (2015) at GeV energies. The size of the gamma-ray emission would strongly depend on the dark matter particle model. Therefore, for the sake of simplicity, the dark matter scenario is not covered in this study.

Among the discussed scenarios for an extended VHE gamma-ray signal in the low state of M87, we consider the neutral pion decay gamma-ray emission from p-p interactions as the most likely one. To model the CR pressure and the gamma-ray emission from pion decay in the inner Virgo Cluster, the CRp energy and spatial distributions as well as the ICM density distribution have to be considered. Very little is known about CRs in galaxy clusters and assumptions have to be made about their energy and spatial distribution. In contrast, the ICM density distribution can be estimated using X-ray measurements assuming an element composition of the plasma (Appendix B). However, X-ray measurements provide no information about the

spatial distribution of the emission along the line of sight, which makes it challenging to construct a 3D model of the ICM without making further assumptions. Additionally, the X-ray surface brightness of the inner Virgo Cluster exhibits a complex morphology, featuring X-ray cavities that coincide with the inner radio lobes, as well as X-ray arcs with the matter at lower temperatures ($\approx$1 keV) than the surrounding gas (Young et al. 2002). The composition of X-ray cavities and the nature of the sustaining pressure in radio galaxies are still not fully understood (Abdulla et al. 2019). Therefore, we adopt a simplified assumption of a radially symmetric ICM distribution, which does not account for the complex X-ray morphology of the cavities in the inner radio lobes. This approach may result in an overestimation of the gamma rays produced by neutral pion decay in the cavities if they are depleted of target material, which could ultimately affect the final UL on the CRp pressure. Furthermore, averaging the 2D X-ray brightness in azimuth leads to a $1\sigma$ uncertainty on the radial X-ray brightness of approximately 85% in the region close to the M87 core (within 4 kpc) and approximately 30-40% in the region up to 0.5° from the M87 core. Since these uncertainties propagate to the ICM density distribution, the limits derived in Sect. 4.2 and Sect. 4.3 are approximations for the case of a radially symmetric ICM distribution.

We assume that CRp are distributed in momentum space according to a PL with spectral index $\alpha_p$, which we vary from 2.1 to 2.6. PL distributions with spectral indices $\alpha_p \ge$2 are characteristic of first-order Fermi acceleration (Fermi 1949; Rieger et al. 2007). Due to the long cooling time of CRp in the ICM (Enßlin et al. 2011), no exponential cutoff is expected in the highest-energy part of the particle spectrum. The formalism utilized in this work to interpret the morphology results in terms of the CR to thermal pressure ratio is presented in Appendix B and is largely based on Appendix B2 of Jacob & Pfrommer (2017b). The parametrizations for the ICM density and temperature distributions in the Virgo Cluster shown in Appendix B are extracted from Jacob & Pfrommer (2017a,b).

Finally, we consider two multiple-component models that lead to an extended gamma-ray component in the M87 low state: the magnetic confinement and the steady-state models. Both models are discussed and derived in Appendices C.1 and C.2, respectively. We define hybrid models composed of one diffuse model and a point-like component to account for the leptonic emission from the AGN. The hybrid models have two free parameters: the amplitudes of the template and the point-like component. The center of the point-like component is fixed to the best-fit position of the pure point-like model, given in Table 3. In the following, we show the results of the fit of the hybrid models to the M87 low-state data and discuss the implications of the results.

## 4.2. The magnetic confinement model

We fitted the hybrid model composed of the magnetic confinement model and a point-like source at the core of M87 to the H.E.S.S. data. The fit converges to a point-like model, showing that an additional diffuse component in the morphology model does not improve the fit. To place an UL on the gamma-ray flux from the diffuse component, we first defined a set of hybrid models with increasing fixed amplitude of the diffuse component. The new hybrid models have now only one free parameter accounting for the intensity of the point-like component. Finally, we fitted the H.E.S.S. data to this defined set of hybrid models and compared their TS with the TS of the fit of a pure point-like model. The maximum intensity of the diffuse component
gamma-ray spectrum at energies $\ge$10 TeV (Franceschini et al. 2019). In fact, only 9 events above the background level are reconstructed with energies above 10 TeV in the low state, 4 in the intermediate state, and 13 in the high state. Therefore, the contribution of such an extended gamma-ray component in the low state of M87 can be safely neglected for morphological studies.

The CRp from the jet could also interact with the local photon fields and produce neutral pions, which could lead to an extended gamma-ray signal. However, Boughelilba et al. (2022) have shown that the accretion flow and the disk around the SMBH can be neglected as targets for particle-photon interactions in the jet.

Another strong candidate for accelerating CRs and producing an extended gamma-ray signal is cluster mergers. Giant radio halos ($\ge$200 kpc) and relics (Enßlin et al. 2011) are typically found in clusters with recent merging activities. Nevertheless, the Virgo Cluster does not have any of the aforementioned signals, although the cluster is also not yet completely dynamically relaxed (Kashibadze et al. 2020). During merger events, CRp are expected to be accelerated at shocks near the cluster's virial radius ($\approx$1.7 Mpc, that is, $\approx$6° for the Virgo Cluster; Ackermann et al. 2015), and they may be transported toward the cluster center, depending on the cluster's turbulent history. In the absence of turbulent advective transport that counterbalances diffusion, the CRp profile in galaxy clusters tends to flatten (Enßlin et al. 2011). Hadronically interacting with the local ICM, these CRp can generate a diffuse gamma-ray signal up to the cluster's virial radius. However, it would be difficult for H.E.S.S. to detect such an extended gamma-ray signal, as its sensitivity degrades significantly to almost 80% of the Crab flux for a source diameter of $\approx$2° (Casanova & Dingus 2008; Mitchell et al. 2023).

In contrast to merger events, accretion events, that is, merging smaller virialized objects, might produce stronger shocks and accelerate CRs more efficiently. Inoue et al. (2005) have shown that for a ~µG magnetic field in the shock region, CRp can be accelerated up to ~$10^{18}$ eV. These CRp would produce pairs of electron-positron, which would rapidly cool through inverse Compton scattering and synchrotron losses, ultimately generating a signal in the form of gamma rays with an energy flux of $10^{-12}-10^{-11}$ erg cm$^{-2}$ s$^{-1}$, marginally detectable by H.E.S.S given the extension of the signal. For the Virgo Cluster the shock diameter can reach up to $\approx$1.9° around the virial radius ($\approx$6°), that is, outside the region of interest of our analysis.

Finally, the annihilation of hypothetical weakly interacting massive particles in the dark matter halo around M87 could also contribute to an extended gamma-ray signal, as investigated in Ackermann et al. (2015) at GeV energies. The size of the gamma-ray emission would strongly depend on the dark matter particle model. Therefore, for the sake of simplicity, the dark matter scenario is not covered in this study.

Among the discussed scenarios for an extended VHE gamma-ray signal in the low state of M87, we consider the neutral pion decay gamma-ray emission from p-p interactions as the most likely one. To model the CR pressure and the gamma-ray emission from pion decay in the inner Virgo Cluster, the CRp energy and spatial distributions as well as the ICM density distribution have to be considered. Very little is known about CRs in galaxy clusters and assumptions have to be made about their energy and spatial distribution. In contrast, the ICM density distribution can be estimated using X-ray measurements assuming an element composition of the plasma (Appendix B). However, X-ray measurements provide no information about the

spatial distribution of the emission along the line of sight, which makes it challenging to construct a 3D model of the ICM without making further assumptions. Additionally, the X-ray surface brightness of the inner Virgo Cluster exhibits a complex morphology, featuring X-ray cavities that coincide with the inner radio lobes, as well as X-ray arcs with the matter at lower temperatures ($\approx$1 keV) than the surrounding gas (Young et al. 2002). The composition of X-ray cavities and the nature of the sustaining pressure in radio galaxies are still not fully understood (Abdulla et al. 2019). Therefore, we adopt a simplified assumption of a radially symmetric ICM distribution, which does not account for the complex X-ray morphology of the cavities in the inner radio lobes. This approach may result in an overestimation of the gamma rays produced by neutral pion decay in the cavities if they are depleted of target material, which could ultimately affect the final UL on the CRp pressure. Furthermore, averaging the 2D X-ray brightness in azimuth leads to a $1\sigma$ uncertainty on the radial X-ray brightness of approximately 85% in the region close to the M87 core (within 4 kpc) and approximately 30-40% in the region up to 0.5° from the M87 core. Since these uncertainties propagate to the ICM density distribution, the limits derived in Sect. 4.2 and Sect. 4.3 are approximations for the case of a radially symmetric ICM distribution.

We assume that CRp are distributed in momentum space according to a PL with spectral index $\alpha_p$, which we vary from 2.1 to 2.6. PL distributions with spectral indices $\alpha_p \ge$2 are characteristic of first-order Fermi acceleration (Fermi 1949; Rieger et al. 2007). Due to the long cooling time of CRp in the ICM (Enßlin et al. 2011), no exponential cutoff is expected in the highest-energy part of the particle spectrum. The formalism utilized in this work to interpret the morphology results in terms of the CR to thermal pressure ratio is presented in Appendix B and is largely based on Appendix B2 of Jacob & Pfrommer (2017b). The parametrizations for the ICM density and temperature distributions in the Virgo Cluster shown in Appendix B are extracted from Jacob & Pfrommer (2017a,b).

Finally, we consider two multiple-component models that lead to an extended gamma-ray component in the M87 low state: the magnetic confinement and the steady-state models. Both models are discussed and derived in Appendices C.1 and C.2, respectively. We define hybrid models composed of one diffuse model and a point-like component to account for the leptonic emission from the AGN. The hybrid models have two free parameters: the amplitudes of the template and the point-like component. The center of the point-like component is fixed to the best-fit position of the pure point-like model, given in Table 3. In the following, we show the results of the fit of the hybrid models to the M87 low-state data and discuss the implications of the results.

## 4.2. The magnetic confinement model

We fitted the hybrid model composed of the magnetic confinement model and a point-like source at the core of M87 to the H.E.S.S. data. The fit converges to a point-like model, showing that an additional diffuse component in the morphology model does not improve the fit. To place an UL on the gamma-ray flux from the diffuse component, we first defined a set of hybrid models with increasing fixed amplitude of the diffuse component. The new hybrid models have now only one free parameter accounting for the intensity of the point-like component. Finally, we fitted the H.E.S.S. data to this defined set of hybrid models and compared their TS with the TS of the fit of a pure point-like model. The maximum intensity of the diffuse component
allowed by the morphology fit (at a 99.7% c.l.) is found when the ΔTS of the morphology fit reaches $3\sigma$ preference for the point-like model in comparison to the hybrid model. The results of our analysis show that the $3\sigma$ preference is reached when the template component accounts for 45% of the total low-state flux of M87. This translates into a maximum gamma-ray flux above 300 GeV of $\le 6.7 \times 10^{-13} \text{ cm}^{-2}\text{s}^{-1}$ for an extended gamma-ray component in the magnetic confinement scenario. Likewise, the contribution of the point-like component is estimated to account for $\ge 55\%$ of the M87 low-state emission.

To derive the maximum CR pressure allowed by the morphology fit of the hybrid model, the CRp distribution has to be taken into account. A centrally peaked CRp distribution is expected since we probe the source of CRp at the center of the cluster, that is, the AGN. Based on the equations from Appendix B (Jacob & Pfrommer 2017a,b), which describes the assumed ICM composition, density, and temperature distributions, we estimated the CRp distribution such that the CR to thermal pressure ratio ($X_{CR}(r) = P_{CR}(r)/P_{th}(r)$) is constant in the inner region of the cluster. A constant $X_{CR}$ in the inner parts of the cluster is characteristic of an equilibrium state between the heating of the ICM by the streaming CRp and the ICM thermal cooling (Pfrommer 2013). The results show that the CR to thermal pressure ratio is $X_{CR} \le 0.17$ for a proton distribution with $\alpha_p = 2.1$ at 99.7% c.l. This estimate is directly influenced by the assumptions on the CRp and ICM spatial distributions. If large regions are material-depleted, as is possibly the case in the X-ray cavities (Young et al. 2002), the UL on the $X_{CR}$ obtained in this study is underestimated, at least in the region of the cavities.

The CRp energy density $e_{CR}$ was estimated and integrated within the volume to yield the total energy in CRp in the inner 20 kpc of the Virgo Cluster of $U_{CR} \le 5 \times 10^{58}$ erg. This is twice as much energy as the total energy estimated from a theoretical model consisting of shocks produced by outbursts that explains the M87 radio and X-ray emissions (Forman et al. 2017). Brüggen et al. (2002) has shown that the buoyant gas in a galaxy cluster can reach a distance of $\approx 20$ kpc after $\approx 15$ Myr. Therefore, for the streaming CRp to reach $U_{CR} \le 5 \times 10^{58}$ erg in 15 Myr, an average jet power of $P_j \le 6 \times 10^{43}$ erg s$^{-1}$ is necessary, considering the extreme case of 100% efficiency in the CRp acceleration. Even though the estimates on $U_{CR}$ and $P_j$ are rather uncertain and model dependent, the previous arguments show that these ULs are larger than (but on the order of magnitude of) the values expected from the literature (Forman et al. 2017). Furthermore, the UL on the $U_{CR}$ derived in this work is twice as constraining as the result from Nizamov & Pshirkov (2022), which placed an UL on $U_{CR}$ up to 35 kpc from the M87 core based on Fermi-LAT data (Abdollahi et al. 2020). Our estimates above rely on the assumption that the spectral index of the CRp is $\alpha_p = 2.1$. Since the spectral index of the gamma rays from neutral pions is expected to follow the index of the CRp population, an additional brighter and steeper component in the M87 low-state emission would be necessary to account for the overall gamma-ray flux estimated with H.E.S.S. ($\alpha_\gamma = 2.63 \pm 0.09$; Table 2), possibly, the kiloparsec jet as in Centaurus A (Abdalla et al. 2020a). For steeper CRp spectral indices our ULs on $X_{CR}$ and $U_{CR}$ become less constraining as seen in Fig. 4 by the blue ULs.

### 4.3. The steady-state model

The fit of the M87 low-state data to the steady-state hybrid model converged into a pure point-like model. Nevertheless, by fitting a set of hybrid models with increasing fixed amplitude of

<figure><img src="image_4.png" /><figcaption>**Fig. 4.** UL of the maximum CR to thermal pressure ratio, $X_{CR,max.} = \max(P_{CR}(r)/P_{th}(r))$, within the inner 20 kpc of the cluster for different proton spectral indices. Two different CRp spatial distributions were considered for a set of CRp spectral indices. The blue ULs are the result of the magnetic confinement approach (Sect. 4.2), and the orange ULs are the result of the steady-state approach (Sect. 4.3). The prediction by Jacob & Pfrommer (2017a) is also shown, in dashed black.</figcaption></figure>

the diffuse component, we derive a 99.7% c.l. UL on the total flux of the diffuse component of $\le 8.1 \times 10^{-13} \text{ cm}^{-2}\text{s}^{-1}$, that is, $\le 55\%$ of the low state of M87 VHE gamma-ray flux (Table 2). This UL is a factor of $\approx 4$ larger than predicted by the model. Hence, the morphology fit was not able to probe the model proposed in Jacob & Pfrommer (2017a) for the Virgo Cluster. However, we could probe a scenario in which the gamma-ray distribution from neutral pion decay resembles the one from the steady-state model but is normalized to $\approx 55\%$ of the flux of M87 low state. Afterward, we derived an UL on the maximum CR to thermal pressure ratio ($X_{CR,max.}$) in the inner Virgo Cluster for the set of CRp spectral indices. The results are similar to the results of the first approach as shown by the orange ULs in Fig. 4. The same approach as in Sect. 4.2 is used here to account for the ICM and CRp spatial distributions. For a proton distribution with $\alpha_p = 2.1$, $X_{CR,max.} \le 0.32$ and the total energy in CRp in the inner 20 kpc of the Virgo Cluster is also constrained to $U_{CR} \le 5 \times 10^{58}$ erg. While the steady-state model predicts a level of $X_{CR,max.} \approx 0.10$ (Fig. A1 in Jacob & Pfrommer 2017a), the H.E.S.S. UL is $X_{CR,max.} \le 0.32$ for $\alpha_p = 2.1$ and becomes less constraining for steeper proton distributions, where the expected flux in VHE gamma rays is lower. Therefore, the H.E.S.S. $X_{CR,max.}$ UL derived from the morphology fit of the M87 gamma-ray low state does not rule out the steady-state model regardless of the spectral index of the proton distribution.

## 5. Summary and conclusions

In this work we aim to localize the VHE gamma-ray emission from M87, probing an extended emission in its low state and testing its hadronic origin. Neutral pions are produced in p-p interactions between the relativistic protons (CRp) from the jet and the ICM (Churazov et al. 2002; Peterson & Fabian 2006). The neutral pions decay almost immediately to gamma rays, which can be detected with H.E.S.S. as an extended and steady gamma-ray signal (Jacob & Pfrommer 2017a,b).

First, we investigated the VHE gamma-ray flux of M87 with H.E.S.S. between 2004 and 2021 and classified the source emission into low, intermediate, and high flux states based on a Bayesian block analysis (Scargle 1998; Scargle et al. 2013;
allowed by the morphology fit (at a 99.7% c.l.) is found when the ΔTS of the morphology fit reaches $3\sigma$ preference for the point-like model in comparison to the hybrid model. The results of our analysis show that the $3\sigma$ preference is reached when the template component accounts for 45% of the total low-state flux of M87. This translates into a maximum gamma-ray flux above 300 GeV of $\le 6.7 \times 10^{-13} \text{ cm}^{-2}\text{s}^{-1}$ for an extended gamma-ray component in the magnetic confinement scenario. Likewise, the contribution of the point-like component is estimated to account for $\ge 55\%$ of the M87 low-state emission.

To derive the maximum CR pressure allowed by the morphology fit of the hybrid model, the CRp distribution has to be taken into account. A centrally peaked CRp distribution is expected since we probe the source of CRp at the center of the cluster, that is, the AGN. Based on the equations from Appendix B (Jacob & Pfrommer 2017a,b), which describes the assumed ICM composition, density, and temperature distributions, we estimated the CRp distribution such that the CR to thermal pressure ratio ($X_{CR}(r) = P_{CR}(r)/P_{th}(r)$) is constant in the inner region of the cluster. A constant $X_{CR}$ in the inner parts of the cluster is characteristic of an equilibrium state between the heating of the ICM by the streaming CRp and the ICM thermal cooling (Pfrommer 2013). The results show that the CR to thermal pressure ratio is $X_{CR} \le 0.17$ for a proton distribution with $\alpha_p = 2.1$ at 99.7% c.l. This estimate is directly influenced by the assumptions on the CRp and ICM spatial distributions. If large regions are material-depleted, as is possibly the case in the X-ray cavities (Young et al. 2002), the UL on the $X_{CR}$ obtained in this study is underestimated, at least in the region of the cavities.

The CRp energy density $e_{CR}$ was estimated and integrated within the volume to yield the total energy in CRp in the inner 20 kpc of the Virgo Cluster of $U_{CR} \le 5 \times 10^{58}$ erg. This is twice as much energy as the total energy estimated from a theoretical model consisting of shocks produced by outbursts that explains the M87 radio and X-ray emissions (Forman et al. 2017). Brüggen et al. (2002) has shown that the buoyant gas in a galaxy cluster can reach a distance of $\approx 20$ kpc after $\approx 15$ Myr. Therefore, for the streaming CRp to reach $U_{CR} \le 5 \times 10^{58}$ erg in 15 Myr, an average jet power of $P_j \le 6 \times 10^{43}$ erg s$^{-1}$ is necessary, considering the extreme case of 100% efficiency in the CRp acceleration. Even though the estimates on $U_{CR}$ and $P_j$ are rather uncertain and model dependent, the previous arguments show that these ULs are larger than (but on the order of magnitude of) the values expected from the literature (Forman et al. 2017). Furthermore, the UL on the $U_{CR}$ derived in this work is twice as constraining as the result from Nizamov & Pshirkov (2022), which placed an UL on $U_{CR}$ up to 35 kpc from the M87 core based on Fermi-LAT data (Abdollahi et al. 2020). Our estimates above rely on the assumption that the spectral index of the CRp is $\alpha_p = 2.1$. Since the spectral index of the gamma rays from neutral pions is expected to follow the index of the CRp population, an additional brighter and steeper component in the M87 low-state emission would be necessary to account for the overall gamma-ray flux estimated with H.E.S.S. ($\alpha_\gamma = 2.63 \pm 0.09$; Table 2), possibly, the kiloparsec jet as in Centaurus A (Abdalla et al. 2020a). For steeper CRp spectral indices our ULs on $X_{CR}$ and $U_{CR}$ become less constraining as seen in Fig. 4 by the blue ULs.

### 4.3. The steady-state model

The fit of the M87 low-state data to the steady-state hybrid model converged into a pure point-like model. Nevertheless, by fitting a set of hybrid models with increasing fixed amplitude of

<figure><img src="image_5.png" /><figcaption>**Fig. 4.** UL of the maximum CR to thermal pressure ratio, $X_{CR,max.} = \max(P_{CR}(r)/P_{th}(r))$, within the inner 20 kpc of the cluster for different proton spectral indices. Two different CRp spatial distributions were considered for a set of CRp spectral indices. The blue ULs are the result of the magnetic confinement approach (Sect. 4.2), and the orange ULs are the result of the steady-state approach (Sect. 4.3). The prediction by Jacob & Pfrommer (2017a) is also shown, in dashed black.</figcaption></figure>

the diffuse component, we derive a 99.7% c.l. UL on the total flux of the diffuse component of $\le 8.1 \times 10^{-13} \text{ cm}^{-2}\text{s}^{-1}$, that is, $\le 55\%$ of the low state of M87 VHE gamma-ray flux (Table 2). This UL is a factor of $\approx 4$ larger than predicted by the model. Hence, the morphology fit was not able to probe the model proposed in Jacob & Pfrommer (2017a) for the Virgo Cluster. However, we could probe a scenario in which the gamma-ray distribution from neutral pion decay resembles the one from the steady-state model but is normalized to $\approx 55\%$ of the flux of M87 low state. Afterward, we derived an UL on the maximum CR to thermal pressure ratio ($X_{CR,max.}$) in the inner Virgo Cluster for the set of CRp spectral indices. The results are similar to the results of the first approach as shown by the orange ULs in Fig. 4. The same approach as in Sect. 4.2 is used here to account for the ICM and CRp spatial distributions. For a proton distribution with $\alpha_p = 2.1$, $X_{CR,max.} \le 0.32$ and the total energy in CRp in the inner 20 kpc of the Virgo Cluster is also constrained to $U_{CR} \le 5 \times 10^{58}$ erg. While the steady-state model predicts a level of $X_{CR,max.} \approx 0.10$ (Fig. A1 in Jacob & Pfrommer 2017a), the H.E.S.S. UL is $X_{CR,max.} \le 0.32$ for $\alpha_p = 2.1$ and becomes less constraining for steeper proton distributions, where the expected flux in VHE gamma rays is lower. Therefore, the H.E.S.S. $X_{CR,max.}$ UL derived from the morphology fit of the M87 gamma-ray low state does not rule out the steady-state model regardless of the spectral index of the proton distribution.

## 5. Summary and conclusions

In this work we aim to localize the VHE gamma-ray emission from M87, probing an extended emission in its low state and testing its hadronic origin. Neutral pions are produced in p-p interactions between the relativistic protons (CRp) from the jet and the ICM (Churazov et al. 2002; Peterson & Fabian 2006). The neutral pions decay almost immediately to gamma rays, which can be detected with H.E.S.S. as an extended and steady gamma-ray signal (Jacob & Pfrommer 2017a,b).

First, we investigated the VHE gamma-ray flux of M87 with H.E.S.S. between 2004 and 2021 and classified the source emission into low, intermediate, and high flux states based on a Bayesian block analysis (Scargle 1998; Scargle et al. 2013;
Ahnen et al. 2016). We focus our studies on the low state since the detection of steady and extended emission could point to a hadronic origin and provide an estimate of the CRp pressure in the inner Virgo Cluster.

We do not detect extended emission via the morphology fit of the low state. Nevertheless, we derive an UL on the $\sigma_G$ of a rotationally symmetric 2D Gaussian model of $58''$ ($\approx 4.6$ kpc) at the 99.7% c.l. The best-fit position of the point-like source model is compatible with the radio core (Fig. 2) within $3\sigma$ statistical uncertainty. Furthermore, our extension UL is twice as constraining as the latest result (MAGIC Collaboration 2020), and, considering the uncertainties in the best-fit position, it excludes for the first time the M 87 radio lobes ($\approx 30$ kpc; Owen et al. 2000) as the main contributor to the low state of M 87 gamma-ray emission (Fig. 2, left). On the other hand, the inner radio cocoon (VLA 21 cm), as shown in Fig. 2 (right), cannot be ruled out as the principal component. Our UL on the extension already lies within the optical extent of M 87 ($R_{1/2} \approx 7.2$ kpc; Weil et al. 1997). The origin of the VHE gamma-ray emission from the M 87 source states is, within the uncertainties, consistent with a single origin at the M 87 core.

We consider two plausible multiple-component scenarios to explain the gamma-ray emission in the low state of M 87. We derived templates for the diffuse gamma-ray components of these two approaches and defined two hybrid models composed of a point-like component centered at M 87 to account for the AGN emission and each of the respective diffuse templates. The morphology fit of the hybrid models allow us to constrain the CRp pressure in the inner Virgo Cluster. The first approach (Appendix C.1) is based on the magnetic confinement model, and the second approach (Appendix C.2) is based on the steady-state model from Jacob & Pfrommer (2017a). The results show that the contribution of the diffuse component is constrained at the 99.7% c.l. to $\le 45\%$ of the VHE gamma-ray flux detected from the low state of M 87 for the magnetic confinement hybrid model and $\le 55\%$ for the steady-state hybrid model. To interpret these limits in terms of CRp pressure in the inner Virgo Cluster, we first considered the CRp distributed as a PL in momentum with a spectral index, $\alpha_p$, that we varied from 2.1 to 2.6. We utilized the formalism and the ICM parametrizations from Jacob & Pfrommer (2017a,b) for the ICM distribution to derive ULs on the maximum CR to thermal pressure ratio ($X_{CR,max.}$) and the maximum energy in CRp ($U_{CR}$) in the region up to 20 kpc from the cluster center. For a CRp distribution with spectral index $\alpha_p = 2.1$, $X_{CR,max.} \le 0.17$ for the hybrid models with the magnetic confinement template, while $X_{CR,max.} \le 0.32$ for the hybrid model with the steady-state template. For steeper CRp distributions, the H.E.S.S. UL is less constraining. The energy in CRp, assuming $\alpha_p = 2.1$, is constrained at the 99.7% c.l. to $U_{CR} \le 5 \times 10^{58}$ erg up to 20 kpc from the M 87 core in both approaches. This limit is larger than, but of the same order of magnitude as, the total energy expected from a theoretical model consisting of shocks produced by outbursts from M 87 (Forman et al. 2017). Our UL on the $U_{CR}$ is also twice as constraining as the UL based on Fermi-LAT data (Nizamov & Pshirkov 2022).

**Acknowledgements.** The support of the Namibian authorities and of the University of Namibia in facilitating the construction and operation of H.E.S.S. is gratefully acknowledged, as is the support by the German Ministry for Education and Research (BMBF), the Max Planck Society, the German Research Foundation (DFG), the Helmholtz Association, the Alexander von Humboldt Foundation, the French Ministry of Higher Education, Research and Innovation, the Centre National de la Recherche Scientifique (CNRS/IN2P3 and CNRS/INSU), the Commissariat à l'énergie atomique et aux énergies alternatives (CEA), the UK Science and Technology Facilities Council (STFC), the Irish Research

Council (IRC) and the Science Foundation Ireland (SFI), the Knut and Alice Wallenberg Foundation, the Polish Ministry of Education and Science, agreement no. 2021/WK/06, the South African Department of Science and Technology and National Research Foundation, the University of Namibia, the National Commission on Research, Science & Technology of Namibia (NCRST), the Austrian Federal Ministry of Education, Science and Research and the Austrian Science Fund (FWF), the Australian Research Council (ARC), the Japan Society for the Promotion of Science, the University of Amsterdam and the Science Committee of Armenia grant 21AG-1C085. We appreciate the excellent work of the technical support staff in Berlin, Zeuthen, Heidelberg, Palaiseau, Paris, Saclay, Tübingen and in Namibia in the construction and operation of the equipment. This work benefited from services provided by the H.E.S.S. Virtual Organisation, supported by the national resource providers of the EGI Federation. Finally, we thank De Gasperin who kindly provided the LOFAR radio emission of M 87 at 140 MHz (Fig. 5 from De Gasperin et al. 2012) and Pfrommer for the discussions on his model.

## References

* Abdalla, H., Adam, R., Aharonian, F., et al. 2020a, Nature, 582, 356

* Abdalla, H., Aharonian, F., Ait Benkhali, F., et al. 2020b, Nat. Astron., 4, 167

* Abdollahi, S., Acero, F., Ackermann, M., et al. 2020, ApJS, 247, 33

* Abdulla, Z., Carlstrom, J. E., Mantz, A. B., et al. 2019, ApJ, 871, 195

* Abramowski, A., Acero, F., Aharonian, F., et al. 2012, ApJ, 746, 151

* Acciari, V. A., Beilicke, M., Blaylock, G., et al. 2008, ApJ, 679, 397

* Acciari, V. A., Aliu, E., Arlen, T., et al. 2009, Science, 325, 444

* Acero, F., Aharonian, F., Akhperjanian, A. G., et al. 2010, MNRAS, 402, 1877

* Ackermann, M., Ajello, M., Albert, A., et al. 2014, ApJ, 787, 18

* Ackermann, M., Ajello, M., Albert, A., et al. 2015, ApJ, 812, 159

* Ackermann, M., Ajello, M., Baldini, L., et al. 2016, ApJ, 826, 1

* Adam, R., Goksu, H., Brown, S., Rudnick, L., & Ferrari, C. 2021, A&A, 648, A60

* Aharonian, F., Akhperjanian, A. G., Bazer-Bachi, A. R., et al. 2006a, Science, 314, 1424

* Aharonian, F., Akhperjanian, A. G., Bazer-Bachi, A. R., et al. 2006b, A&A, 457, 899

* Aharonian, F., Akhperjanian, A. G., Anton, G., et al. 2009a, A&A, 495, 27

* Aharonian, F., Akhperjanian, A. G., Anton, G., et al. 2009b, A&A, 502, 437

* Ahnen, M. L., Ansoldi, S., Antonelli, L. A., et al. 2016, A&A, 593, A91

* Aleksić, J., Alvarez, E. A., Antonelli, L. A., et al. 2012, A&A, 541, A99

* Alfaro, R., Alvarez, C., Arteaga-Velázquez, J. C., et al. 2022, ApJ, 934, 158

* Anderhub, H., Antonelli, L. A., Antoranz, P., et al. 2010, ApJ, 710, 828

* Arlen, T., Aune, T., Beilicke, M., et al. 2012, ApJ, 757, 123

* Ashton, T., Backes, M., Balzer, A., et al. 2020, Astropart. Phys., 118, 102425

* Astropy Collaboration (Robitaille, T. P., et al.) 2013, A&A, 558, A33

* Astropy Collaboration (Price-Whelan, A. M., et al.) 2018, AJ, 156, 123

* Baghmanyan, V., Zargaryan, D., Aharonian, F., et al. 2022, MNRAS, 516, 562

* Barbiellini, G., Bastieri, D., Bechtol, K., et al. 2014, ApJ, 784, 118

* Barbosa Martins, V. 2022, Ph.D. Thesis, Humboldt-Universität zu Berlin, Germany, https://edoc.hu-berlin.de/handle/18452/25715

* Becker, H. B., White, R. L., & Helfand, D. J. 1995, ApJ, 450, 559

* Berge, D., Funk, S., & Hinton, J. 2007, A&A, 466, 1219

* Bird, S., Harris, W. E., Blakeslee, J. P., & Flynn, C. 2010, A&A, 524, A71

* Blakeslee, J. P., Jordn, A., Mei, S., et al. 2009, ApJ, 694, 556

* Boughelilba, M., Reimer, A., & Merten, L. 2022, ApJ, 938, 79

* Brüggen, M., Kaiser, C. R., Churazov, E., & Enßlin, T. A. 2002, MNRAS, 331, 545

* Casanova, S., & Dingus, B. L. 2008, Astropart. Phys., 29, 63

* Cash, W. 1979, ApJ, 228, 939

* Churazov, E., Brüggen, M., Kaiser, C. R., Bohringer, H., & Forman, W. 2001, ApJ, 554, 261

* Churazov, E., Sunyaev, R., Forman, W., & Böhringer, H. 2002, MNRAS, 332, 729

* De Gasperin, F., Orrú, E., Murgia, M., et al. 2012, A&A, 547, A56

* Deil, C., Lefaucheur, J., Zanin, R., et al. 2017, in Proc. 35th Int. Cosmic RayConf. – PoS(ICRC2017), 35, 766

* de Naurois, M., & Rolland, L. 2009, Astropart. Phys., 32, 231

* Domainko, W., Nedbal, D., Hinton, J. A., & Martineau-Huynh, O. 2009, Int. J. Mod. Phys. D, 18, 1627

* Enßlin, T., Pfrommer, C., Miniati, F., & Subramanian, K. 2011, A&A, 527, 1

* Fermi, E. 1949, Phys. Rev., 75, 1169

* Forman, W., Churazov, E., Jones, C., et al. 2017, ApJ, 844, 122

* Franceschini, A., Foffano, L., Prandini, E., & Tavecchio, F. 2019, A&A, 629, A2

* Freeman, P., Doe, S., & Siemiginowska, A. 2001, SPIE, 4477, 76

* Galante, N. 2009, in Proc. 31th Int. Cosmic Ray Conf. – PoS(ICRC2009), 37, 4

* Guo, F., & Oh, S. P. 2008, MNRAS, 384, 251
Ahnen et al. 2016). We focus our studies on the low state since the detection of steady and extended emission could point to a hadronic origin and provide an estimate of the CRp pressure in the inner Virgo Cluster.

We do not detect extended emission via the morphology fit of the low state. Nevertheless, we derive an UL on the $\sigma_G$ of a rotationally symmetric 2D Gaussian model of $58''$ ($\approx 4.6$ kpc) at the 99.7% c.l. The best-fit position of the point-like source model is compatible with the radio core (Fig. 2) within $3\sigma$ statistical uncertainty. Furthermore, our extension UL is twice as constraining as the latest result (MAGIC Collaboration 2020), and, considering the uncertainties in the best-fit position, it excludes for the first time the M 87 radio lobes ($\approx 30$ kpc; Owen et al. 2000) as the main contributor to the low state of M 87 gamma-ray emission (Fig. 2, left). On the other hand, the inner radio cocoon (VLA 21 cm), as shown in Fig. 2 (right), cannot be ruled out as the principal component. Our UL on the extension already lies within the optical extent of M 87 ($R_{1/2} \approx 7.2$ kpc; Weil et al. 1997). The origin of the VHE gamma-ray emission from the M 87 source states is, within the uncertainties, consistent with a single origin at the M 87 core.

We consider two plausible multiple-component scenarios to explain the gamma-ray emission in the low state of M 87. We derived templates for the diffuse gamma-ray components of these two approaches and defined two hybrid models composed of a point-like component centered at M 87 to account for the AGN emission and each of the respective diffuse templates. The morphology fit of the hybrid models allow us to constrain the CRp pressure in the inner Virgo Cluster. The first approach (Appendix C.1) is based on the magnetic confinement model, and the second approach (Appendix C.2) is based on the steady-state model from Jacob & Pfrommer (2017a). The results show that the contribution of the diffuse component is constrained at the 99.7% c.l. to $\le 45\%$ of the VHE gamma-ray flux detected from the low state of M 87 for the magnetic confinement hybrid model and $\le 55\%$ for the steady-state hybrid model. To interpret these limits in terms of CRp pressure in the inner Virgo Cluster, we first considered the CRp distributed as a PL in momentum with a spectral index, $\alpha_p$, that we varied from 2.1 to 2.6. We utilized the formalism and the ICM parametrizations from Jacob & Pfrommer (2017a,b) for the ICM distribution to derive ULs on the maximum CR to thermal pressure ratio ($X_{CR,max}$) and the maximum energy in CRp ($U_{CR}$) in the region up to 20 kpc from the cluster center. For a CRp distribution with spectral index $\alpha_p = 2.1$, $X_{CR,max} \le 0.17$ for the hybrid models with the magnetic confinement template, while $X_{CR,max} \le 0.32$ for the hybrid model with the steady-state template. For steeper CRp distributions, the H.E.S.S. UL is less constraining. The energy in CRp, assuming $\alpha_p = 2.1$, is constrained at the 99.7% c.l. to $U_{CR} \le 5 \times 10^{58}$ erg up to 20 kpc from the M 87 core in both approaches. This limit is larger than, but of the same order of magnitude as, the total energy expected from a theoretical model consisting of shocks produced by outbursts from M 87 (Forman et al. 2017). Our UL on the $U_{CR}$ is also twice as constraining as the UL based on Fermi-LAT data (Nizamov & Pshirkov 2022).

**Acknowledgements.** The support of the Namibian authorities and of the University of Namibia in facilitating the construction and operation of H.E.S.S. is gratefully acknowledged, as is the support by the German Ministry for Education and Research (BMBF), the Max Planck Society, the German Research Foundation (DFG), the Helmholtz Association, the Alexander von Humboldt Foundation, the French Ministry of Higher Education, Research and Innovation, the Centre National de la Recherche Scientifique (CNRS/IN2P3 and CNRS/INSU), the Commissariat à l'énergie atomique et aux énergies alternatives (CEA), the UK Science and Technology Facilities Council (STFC), the Irish Research

Council (IRC) and the Science Foundation Ireland (SFI), the Knut and Alice Wallenberg Foundation, the Polish Ministry of Education and Science, agreement no. 2021/WK/06, the South African Department of Science and Technology and National Research Foundation, the University of Namibia, the National Commission on Research, Science & Technology of Namibia (NCRST), the Austrian Federal Ministry of Education, Science and Research and the Austrian Science Fund (FWF), the Australian Research Council (ARC), the Japan Society for the Promotion of Science, the University of Amsterdam and the Science Committee of Armenia grant 21AG-1C085. We appreciate the excellent work of the technical support staff in Berlin, Zeuthen, Heidelberg, Palaiseau, Paris, Saclay, Tübingen and in Namibia in the construction and operation of the equipment. This work benefited from services provided by the H.E.S.S. Virtual Organisation, supported by the national resource providers of the EGI Federation. Finally, we thank De Gasperin who kindly provided the LOFAR radio emission of M 87 at 140 MHz (Fig. 5 from De Gasperin et al. 2012) and Pfrommer for the discussions on his model.

## References

* Abdalla, H., Adam, R., Aharonian, F., et al. 2020a, Nature, 582, 356

* Abdalla, H., Aharonian, F., Ait Benkhali, F., et al. 2020b, Nat. Astron., 4, 167

* Abdollahi, S., Acero, F., Ackermann, M., et al. 2020, ApJS, 247, 33

* Abdulla, Z., Carlstrom, J. E., Mantz, A. B., et al. 2019, ApJ, 871, 195

* Abramowski, A., Acero, F., Aharonian, F., et al. 2012, ApJ, 746, 151

* Acciari, V. A., Beilicke, M., Blaylock, G., et al. 2008, ApJ, 679, 397

* Acciari, V. A., Aliu, E., Arlen, T., et al. 2009, Science, 325, 444

* Acero, F., Aharonian, F., Akhperjanian, A. G., et al. 2010, MNRAS, 402, 1877

* Ackermann, M., Ajello, M., Albert, A., et al. 2014, ApJ, 787, 18

* Ackermann, M., Ajello, M., Albert, A., et al. 2015, ApJ, 812, 159

* Ackermann, M., Ajello, M., Baldini, L., et al. 2016, ApJ, 826, 1

* Adam, R., Goksu, H., Brown, S., Rudnick, L., & Ferrari, C. 2021, A&A, 648, A60

* Aharonian, F., Akhperjanian, A. G., Bazer-Bachi, A. R., et al. 2006a, Science, 314, 1424

* Aharonian, F., Akhperjanian, A. G., Bazer-Bachi, A. R., et al. 2006b, A&A, 457, 899

* Aharonian, F., Akhperjanian, A. G., Anton, G., et al. 2009a, A&A, 495, 27

* Aharonian, F., Akhperjanian, A. G., Anton, G., et al. 2009b, A&A, 502, 437

* Ahnen, M. L., Ansoldi, S., Antonelli, L. A., et al. 2016, A&A, 593, A91

* Aleksić, J., Alvarez, E. A., Antonelli, L. A., et al. 2012, A&A, 541, A99

* Alfaro, R., Alvarez, C., Arteaga-Velázquez, J. C., et al. 2022, ApJ, 934, 158

* Anderhub, H., Antonelli, L. A., Antoranz, P., et al. 2010, ApJ, 710, 828

* Arlen, T., Aune, T., Beilicke, M., et al. 2012, ApJ, 757, 123

* Ashton, T., Backes, M., Balzer, A., et al. 2020, Astropart. Phys., 118, 102425

* Astropy Collaboration (Robitaille, T. P., et al.) 2013, A&A, 558, A33

* Astropy Collaboration (Price-Whelan, A. M., et al.) 2018, AJ, 156, 123

* Baghmanyan, V., Zargaryan, D., Aharonian, F., et al. 2022, MNRAS, 516, 562

* Barbiellini, G., Bastieri, D., Bechtol, K., et al. 2014, ApJ, 784, 118

* Barbosa Martins, V. 2022, Ph.D. Thesis, Humboldt-Universität zu Berlin, Germany, https://edoc.hu-berlin.de/handle/18452/25715

* Becker, H. B., White, R. L., & Helfand, D. J. 1995, ApJ, 450, 559

* Berge, D., Funk, S., & Hinton, J. 2007, A&A, 466, 1219

* Bird, S., Harris, W. E., Blakeslee, J. P., & Flynn, C. 2010, A&A, 524, A71

* Blakeslee, J. P., Jordn, A., Mei, S., et al. 2009, ApJ, 694, 556

* Boughelilba, M., Reimer, A., & Merten, L. 2022, ApJ, 938, 79

* Brüggen, M., Kaiser, C. R., Churazov, E., & Enßlin, T. A. 2002, MNRAS, 331, 545

* Casanova, S., & Dingus, B. L. 2008, Astropart. Phys., 29, 63

* Cash, W. 1979, ApJ, 228, 939

* Churazov, E., Brüggen, M., Kaiser, C. R., Bohringer, H., & Forman, W. 2001, ApJ, 554, 261

* Churazov, E., Sunyaev, R., Forman, W., & Böhringer, H. 2002, MNRAS, 332, 729

* De Gasperin, F., Orrú, E., Murgia, M., et al. 2012, A&A, 547, A56

* Deil, C., Lefaucheur, J., Zanin, R., et al. 2017, in Proc. 35th Int. Cosmic RayConf. – PoS(ICRC2017), 35, 766

* de Naurois, M., & Rolland, L. 2009, Astropart. Phys., 32, 231

* Domainko, W., Nedbal, D., Hinton, J. A., & Martineau-Huynh, O. 2009, Int. J. Mod. Phys. D, 18, 1627

* Enßlin, T., Pfrommer, C., Miniati, F., & Subramanian, K. 2011, A&A, 527, 1

* Fermi, E. 1949, Phys. Rev., 75, 1169

* Forman, W., Churazov, E., Jones, C., et al. 2017, ApJ, 844, 122

* Franceschini, A., Foffano, L., Prandini, E., & Tavecchio, F. 2019, A&A, 629, A2

* Freeman, P., Doe, S., & Siemiginowska, A. 2001, SPIE, 4477, 76

* Galante, N. 2009, in Proc. 31th Int. Cosmic Ray Conf. – PoS(ICRC2009), 37, 4

* Guo, F., & Oh, S. P. 2008, MNRAS, 384, 251
Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357

Hinton, J. A., & Hofmann, W. 2009, ARA&A, 47, 523

Huber, B., Tchernin, C., Eckert, D., et al. 2013, A&A, 560, A64

Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90

Inoue, S., Aharonian, F. A., & Sugiyama, N. 2005, ApJ, 628, L9

Jacob, S., & Pfrommer, C. 2017a, MNRAS, 467, 1449

Jacob, S., & Pfrommer, C. 2017b, MNRAS, 467, 1478

Kashibadze, O. G., Karachentsev, I. D., & Karachentseva, V. E. 2020, A&A, 635, 1

MAGIC Collaboration (Acciari, V. A., et al.) 2020, MNRAS, 492, 5354

Massaro, F., Harris, D. E., Tremblay, G. R., et al. 2013, ApJS, 206, 7

Mathews, W. G., Faltenbacher, A., & Brighenti, F. 2006, ApJ, 638, 659

Mitchell, A. M., Caroff, S., Hinton, J., et al. 2023, A&A, 395, 18248039

Nigro, C., Deil, C., Zanin, R., et al. 2019, A&A, 625, A10

Nizamov, B. A., & Pshirkov, M. S. 2022, JETP, 115, 245

Ohm, S., van Eldik, C., & Egberts, K. 2009, Astropart. Phys., 31, 383

Owen, F. N., Eilek, J. A., & Kassim, N. E. 2000, ApJ, 543, 611

Parsons, R. D., & Hinton, J. A. 2014, Astropart. Phys., 56, 26

Peterson, J. R., & Fabian, A. C. 2006, Phys. Rep., 427, 1

Pfrommer, C. 2013, ApJ, 779, 10

Pfrommer, C., & Enßlin, T. A. 2004, A&A, 426, 777

Pinzke, A., Pfrommer, C., & Bergström, L. 2011, Phys. Rev. D, 84, 123509

Piron, F., Djannati-Atai, A., Punch, M., et al. 2001, A&A, 374, 895

Prokhorov, D. A., & Churazov, E. M. 2014, A&A, 567, A93

Refsdal, B. L., Doe, S. M., Nguyen, D. T., et al. 2009, in Proc. 8th Python Sci. Conf. (SciPy2009), Pasadena, CA, USA, eds. G. Varoquaux, S. van der Walt, & J. Millman, 51

Rieger, F. M., Bosch-Ramon, V., & Duffy, P. 2007, Ap&SS, 309, 119

Ruszkowski, M., Bruggen, M., & Begelman, M. C. 2004, ApJ, 611, 158

Scargle, J. D. 1998, ApJ, 504, 405

Scargle, J. D., Norris, J. P., Jackson, B., & Chiang, J. 2013, ApJ, 764, 167

Simionescu, A., Werner, N., Mantz, A., Allen, S. W., & Urban, O. 2017, MNRAS, 469, 1476

The EHT MWL Science Working Group (Algaba, J. C., et al.) 2021, ApJ, 911, L11

The Event Horizon Telescope Collaboration (Akiyama, K., et al.) 2019, ApJ, 875, L1

The FIRST Project Team 1994, VLA FIRST (1.4 GHz): FIRST, https://skyview.gsfc.nasa.gov/

Thopson, A. R., Clark, B. G., Wade, C. M., & Napier, P. J. 1980, ApJS, 44, 151

Urban, O., Werner, N., Simionescu, A., Allen, S. W., & Böhringer, H. 2011, MNRAS, 414, 2101

Van Haarlem, M. P., Wise, M. W., Gunst, A. W., et al. 2013, A&A, 556, A2

Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nat. Methods, 17, 261

Weil, M. L., Bland-Hawthorn, J., & Malin, D. F. 1997, ApJ, 490, 664

Weisskopf, M. C., Tananbaum, H. D., Speybroeck, L. P. V., & O'Dell, S. L. 2000, SPIE, 4012, 2

Xi, S. Q., Wang, X. Y., Liang, Y. F., et al. 2018, PRD, 98, 063006

Xue, R., Wang, Z. R., & Li, W. J. 2022, Phys. Rev. D, 106, 103021

Young, A. J., Wilson, A. S., & Mundell, C. G. 2002, ApJ, 579, 560

Zakamska, N. L., & Narayan, R. 2003, ApJ, 582, 162

Zhuravleva, I., Churazov, E., Schekochihin, A. A., et al. 2014, Nature, 515, 85

11 Department of Physics and Electrical Engineering, Linnaeus University, 351 95 Växjö, Sweden

12 Institut für Physik, Humboldt-Universität zu Berlin, Newtonstr. 15, 12489 Berlin, Germany

13 Institut für Astronomie und Astrophysik, Universität Tübingen, Sand 1, 72076 Tübingen, Germany

14 Laboratoire Univers et Théories, Observatoire de Paris, Université PSL, CNRS, Université de Paris, 5 Pl. Jules Janssen, 92190 Meudon, France

15 Sorbonne Université, Université Paris Diderot, Sorbonne Paris Cité, CNRS/IN2P3, Laboratoire de Physique Nucléaire et de Hautes Energies, LPNHE, 4 Place Jussieu, 75252 Paris, France

16 IRFU, CEA, Université Paris-Saclay, Bât. 141, 91191 Gif-sur-Yvette, France

17 Friedrich-Alexander-Universität Erlangen-Nürnberg, Erlangen Centre for Astroparticle Physics, Nikolaus-Fiebiger-Str. 2, 91058 Erlangen, Germany

18 Astronomical Observatory, The University of Warsaw, Al. Ujazdowskie 4, 00-478 Warsaw, Poland

19 Université Savoie Mont Blanc, CNRS, Laboratoire d'Annecy de Physique des Particules – IN2P3, 9 Chem. de Bellevue, 74000 Annecy, France

20 Instytut Fizyki Jądrowej PAN, ul. Radzikowskiego 152, 31-342 Kraków, Poland

21 Universität Hamburg, Institut für Experimentalphysik, Luruper Chaussee 149, 22761 Hamburg, Germany

22 School of Physics, University of the Witwatersrand, 1 Jan Smuts Avenue, Braamfontein, Johannesburg 2050, South Africa

23 University of Oxford, Department of Physics, Denys Wilkinson Building, Keble Road, Oxford OX1 3RH, UK

24 Aix Marseille Université, CNRS/IN2P3, CPPM, Marseille, France

25 School of Science, Western Sydney University, Locked Bag 1797, Penrith South DC, NSW 2751, Australia

26 Laboratoire Leprince-Ringuet, École Polytechnique, CNRS, Institut Polytechnique de Paris, 91128 Palaiseau, France

27 Université Bordeaux, CNRS, LP2I Bordeaux, UMR 5797, 19 Chem. du Solarium, 33170 Gradignan, France

28 Obserwatorium Astronomiczne, Uniwersytet Jagielloński, ul. Orla 171, 30-244 Kraków, Poland

29 Institute of Astronomy, Faculty of Physics, Astronomy and Informatics, Nicolaus Copernicus University, Grudziadzka 5, 87-100 Torun, Poland

30 Nicolaus Copernicus Astronomical Center, Polish Academy of Sciences, ul. Bartycka 18, 00-716 Warsaw, Poland

31 Laboratoire Univers et Particules de Montpellier, Université Montpellier, CNRS/IN2P3, CC 72, Place Eugène Bataillon, 34095 Montpellier Cedex 5, France

32 School of Physical Sciences, University of Adelaide, N1.36, Engineering North, North Terrace Campus, Adelaide 5005, Australia

33 Leopold-Franzens-Universität Innsbruck, Institut für Astro- und Teilchenphysik, Technikerstraße 25, 6020 Innsbruck, Austria

34 Department of Physics and Astronomy, The University of Leicester, University Road, Leicester LE1 7RH, UK

35 GRAPPA, Anton Pannekoek Institute for Astronomy, University of Amsterdam, Science Park 904, 1098 XH Amsterdam, The Netherlands

36 Yerevan Physics Institute, 2 Alikhanian Brothers St., 375036 Yerevan, Armenia

37 Department of Physics, Konan University, 8-9-1 Okamoto, Higashinada, Kobe, Hyogo 658-8501, Japan

38 Kavli Institute for the Physics and Mathematics of the Universe (WPI), The University of Tokyo Institutes for Advanced Study (UTIAS), The University of Tokyo, 5-1-5 Kashiwa-no-Ha, Kashiwa, Chiba 277-8583, Japan

39 RIKEN, 2-1 Hirosawa, Wako, Saitama 351-0198, Japan

40 Department of Physics, Rikkyo University, 3-34-1 Nishi-Ikebukuro, Toshima-ku, Tokyo 171-8501, Japan

41 Department of Physics, University of the Free State, PO Box 339, Bloemfontein 9300, South Africa

<footnote n="1">1 Dublin Institute for Advanced Studies, 31 Fitzwilliam Place, Dublin 2, Ireland</footnote>

<footnote n="2">2 Max-Planck-Institut für Kernphysik, PO Box 103980, 69029 Heidelberg, Germany</footnote>

<footnote n="3">3 Yerevan State University, 1 Alek Manukyan St, Yerevan 0025, Armenia</footnote>

<footnote n="4">4 Landessternwarte, Universität Heidelberg, Königstuhl 12, 69117 Heidelberg, Germany</footnote>

<footnote n="5">5 Kapteyn Astronomical Institute, University of Groningen, Landleven 12, 9747 AD Groningen, The Netherlands</footnote>

<footnote n="6">6 University of Namibia, Department of Physics, Private Bag 13301, Windhoek 10005, Namibia</footnote>

<footnote n="7">7 Centre for Space Research, North-West University, Potchefstroom 2520, South Africa</footnote>

<footnote n="8">8 Deutsches Elektronen-Synchrotron DESY, Platanenallee 6, 15738 Zeuthen, Germany</footnote>

<footnote n="9">9 Institut für Physik und Astronomie, Universität Potsdam, Karl-Liebknecht-Strasse 24/25, 14476 Potsdam, Germany</footnote>

<footnote n="10">10 Université de Paris, CNRS, Astroparticule et Cosmologie, 10 Rue Alice Domon et Léonie Duquet, 75013 Paris, France</footnote>
Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357

Hinton, J. A., & Hofmann, W. 2009, ARA&A, 47, 523

Huber, B., Tchernin, C., Eckert, D., et al. 2013, A&A, 560, A64

Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90

Inoue, S., Aharonian, F. A., & Sugiyama, N. 2005, ApJ, 628, L9

Jacob, S., & Pfrommer, C. 2017a, MNRAS, 467, 1449

Jacob, S., & Pfrommer, C. 2017b, MNRAS, 467, 1478

Kashibadze, O. G., Karachentsev, I. D., & Karachentseva, V. E. 2020, A&A, 635, 1

MAGIC Collaboration (Acciari, V. A., et al.) 2020, MNRAS, 492, 5354

Massaro, F., Harris, D. E., Tremblay, G. R., et al. 2013, ApJS, 206, 7

Mathews, W. G., Faltenbacher, A., & Brighenti, F. 2006, ApJ, 638, 659

Mitchell, A. M., Caroff, S., Hinton, J., et al. 2023, A&A, 395, 18248039

Nigro, C., Deil, C., Zanin, R., et al. 2019, A&A, 625, A10

Nizamov, B. A., & Pshirkov, M. S. 2022, JETP, 115, 245

Ohm, S., van Eldik, C., & Egberts, K. 2009, Astropart. Phys., 31, 383

Owen, F. N., Eilek, J. A., & Kassim, N. E. 2000, ApJ, 543, 611

Parsons, R. D., & Hinton, J. A. 2014, Astropart. Phys., 56, 26

Peterson, J. R., & Fabian, A. C. 2006, Phys. Rep., 427, 1

Pfrommer, C. 2013, ApJ, 779, 10

Pfrommer, C., & Enßlin, T. A. 2004, A&A, 426, 777

Pinzke, A., Pfrommer, C., & Bergström, L. 2011, Phys. Rev. D, 84, 123509

Piron, F., Djannati-Atai, A., Punch, M., et al. 2001, A&A, 374, 895

Prokhorov, D. A., & Churazov, E. M. 2014, A&A, 567, A93

Refsdal, B. L., Doe, S. M., Nguyen, D. T., et al. 2009, in Proc. 8th Python Sci. Conf. (SciPy2009), Pasadena, CA, USA, eds. G. Varoquaux, S. van der Walt, & J. Millman, 51

Rieger, F. M., Bosch-Ramon, V., & Duffy, P. 2007, Ap&SS, 309, 119

Ruszkowski, M., Bruggen, M., & Begelman, M. C. 2004, ApJ, 611, 158

Scargle, J. D. 1998, ApJ, 504, 405

Scargle, J. D., Norris, J. P., Jackson, B., & Chiang, J. 2013, ApJ, 764, 167

Simionescu, A., Werner, N., Mantz, A., Allen, S. W., & Urban, O. 2017, MNRAS, 469, 1476

The EHT MWL Science Working Group (Algaba, J. C., et al.) 2021, ApJ, 911, L11

The Event Horizon Telescope Collaboration (Akiyama, K., et al.) 2019, ApJ, 875, L1

The FIRST Project Team 1994, VLA FIRST (1.4 GHz): FIRST, https://skyview.gsfc.nasa.gov/

Thopson, A. R., Clark, B. G., Wade, C. M., & Napier, P. J. 1980, ApJS, 44, 151

Urban, O., Werner, N., Simionescu, A., Allen, S. W., & Böhringer, H. 2011, MNRAS, 414, 2101

Van Haarlem, M. P., Wise, M. W., Gunst, A. W., et al. 2013, A&A, 556, A2

Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nat. Methods, 17, 261

Weil, M. L., Bland-Hawthorn, J., & Malin, D. F. 1997, ApJ, 490, 664

Weisskopf, M. C., Tananbaum, H. D., Speybroeck, L. P. V., & O'Dell, S. L. 2000, SPIE, 4012, 2

Xi, S. Q., Wang, X. Y., Liang, Y. F., et al. 2018, PRD, 98, 063006

Xue, R., Wang, Z. R., & Li, W. J. 2022, Phys. Rev. D, 106, 103021

Young, A. J., Wilson, A. S., & Mundell, C. G. 2002, ApJ, 579, 560

Zakamska, N. L., & Narayan, R. 2003, ApJ, 582, 162

Zhuravleva, I., Churazov, E., Schekochihin, A. A., et al. 2014, Nature, 515, 85

11 Department of Physics and Electrical Engineering, Linnaeus University, 351 95 Växjö, Sweden

12 Institut für Physik, Humboldt-Universität zu Berlin, Newtonstr. 15, 12489 Berlin, Germany

13 Institut für Astronomie und Astrophysik, Universität Tübingen, Sand 1, 72076 Tübingen, Germany

14 Laboratoire Univers et Théories, Observatoire de Paris, Université PSL, CNRS, Université de Paris, 5 Pl. Jules Janssen, 92190 Meudon, France

15 Sorbonne Université, Université Paris Diderot, Sorbonne Paris Cité, CNRS/IN2P3, Laboratoire de Physique Nucléaire et de Hautes Energies, LPNHE, 4 Place Jussieu, 75252 Paris, France

16 IRFU, CEA, Université Paris-Saclay, Bât. 141, 91191 Gif-sur-Yvette, France

17 Friedrich-Alexander-Universität Erlangen-Nürnberg, Erlangen Centre for Astroparticle Physics, Nikolaus-Fiebiger-Str. 2, 91058 Erlangen, Germany

18 Astronomical Observatory, The University of Warsaw, Al. Ujazdowskie 4, 00-478 Warsaw, Poland

19 Université Savoie Mont Blanc, CNRS, Laboratoire d'Annecy de Physique des Particules – IN2P3, 9 Chem. de Bellevue, 74000 Annecy, France

20 Instytut Fizyki Jądrowej PAN, ul. Radzikowskiego 152, 31-342 Kraków, Poland

21 Universität Hamburg, Institut für Experimentalphysik, Luruper Chaussee 149, 22761 Hamburg, Germany

22 School of Physics, University of the Witwatersrand, 1 Jan Smuts Avenue, Braamfontein, Johannesburg 2050, South Africa

23 University of Oxford, Department of Physics, Denys Wilkinson Building, Keble Road, Oxford OX1 3RH, UK

24 Aix Marseille Université, CNRS/IN2P3, CPPM, Marseille, France

25 School of Science, Western Sydney University, Locked Bag 1797, Penrith South DC, NSW 2751, Australia

26 Laboratoire Leprince-Ringuet, École Polytechnique, CNRS, Institut Polytechnique de Paris, 91128 Palaiseau, France

27 Université Bordeaux, CNRS, LP2I Bordeaux, UMR 5797, 19 Chem. du Solarium, 33170 Gradignan, France

28 Obserwatorium Astronomiczne, Uniwersytet Jagielloński, ul. Orla 171, 30-244 Kraków, Poland

29 Institute of Astronomy, Faculty of Physics, Astronomy and Informatics, Nicolaus Copernicus University, Grudziadzka 5, 87-100 Torun, Poland

30 Nicolaus Copernicus Astronomical Center, Polish Academy of Sciences, ul. Bartycka 18, 00-716 Warsaw, Poland

31 Laboratoire Univers et Particules de Montpellier, Université Montpellier, CNRS/IN2P3, CC 72, Place Eugène Bataillon, 34095 Montpellier Cedex 5, France

32 School of Physical Sciences, University of Adelaide, N1.36, Engineering North, North Terrace Campus, Adelaide 5005, Australia

33 Leopold-Franzens-Universität Innsbruck, Institut für Astro- und Teilchenphysik, Technikerstraße 25, 6020 Innsbruck, Austria

34 Department of Physics and Astronomy, The University of Leicester, University Road, Leicester LE1 7RH, UK

35 GRAPPA, Anton Pannekoek Institute for Astronomy, University of Amsterdam, Science Park 904, 1098 XH Amsterdam, The Netherlands

36 Yerevan Physics Institute, 2 Alikhanian Brothers St., 375036 Yerevan, Armenia

37 Department of Physics, Konan University, 8-9-1 Okamoto, Higashinada, Kobe, Hyogo 658-8501, Japan

38 Kavli Institute for the Physics and Mathematics of the Universe (WPI), The University of Tokyo Institutes for Advanced Study (UTIAS), The University of Tokyo, 5-1-5 Kashiwa-no-Ha, Kashiwa, Chiba 277-8583, Japan

39 RIKEN, 2-1 Hirosawa, Wako, Saitama 351-0198, Japan

40 Department of Physics, Rikkyo University, 3-34-1 Nishi-Ikebukuro, Toshima-ku, Tokyo 171-8501, Japan

41 Department of Physics, University of the Free State, PO Box 339, Bloemfontein 9300, South Africa

<footnote n="1">1 Dublin Institute for Advanced Studies, 31 Fitzwilliam Place, Dublin 2, Ireland</footnote>

<footnote n="2">2 Max-Planck-Institut für Kernphysik, PO Box 103980, 69029 Heidelberg, Germany</footnote>

<footnote n="3">3 Yerevan State University, 1 Alek Manukyan St, Yerevan 0025, Armenia</footnote>

<footnote n="4">4 Landessternwarte, Universität Heidelberg, Königstuhl 12, 69117 Heidelberg, Germany</footnote>

<footnote n="5">5 Kapteyn Astronomical Institute, University of Groningen, Landleven 12, 9747 AD Groningen, The Netherlands</footnote>

<footnote n="6">6 University of Namibia, Department of Physics, Private Bag 13301, Windhoek 10005, Namibia</footnote>

<footnote n="7">7 Centre for Space Research, North-West University, Potchefstroom 2520, South Africa</footnote>

<footnote n="8">8 Deutsches Elektronen-Synchrotron DESY, Platanenallee 6, 15738 Zeuthen, Germany</footnote>

<footnote n="9">9 Institut für Physik und Astronomie, Universität Potsdam, Karl-Liebknecht-Strasse 24/25, 14476 Potsdam, Germany</footnote>

<footnote n="10">10 Université de Paris, CNRS, Astroparticule et Cosmologie, 10 Rue Alice Domon et Léonie Duquet, 75013 Paris, France</footnote>
## Appendix A: The long-term light curve of M 87

<figure><img src="image_5.png" /><figcaption>**Fig. A.1.** M 87 long-term light curves with their Bayesian blocks derived according to Sect. 2.2 for 1-day- (top), 7-day- (middle), and 15-day-binned flux points (bottom). The H.E.S.S. data points are shown in black, and a constant fit is shown with the solid orange line. The blocks with their flux levels are given by the blue rectangles. Dashed orange lines indicate a change in the flux level. The orange numbers at the bottom indicate the labels of the blocks.</figcaption></figure>

M 87 is known to be highly variable in VHE gamma-rays, with flares on the order of one day (Abramowski et al. 2012). Therefore, a daily binned light curve could well identify the flares and isolate the low-state periods. Nevertheless, M 87 is a rather weak source in the TeV regime at low states. In fact, at least 5 hours of

observations are needed for a significant source detection (see The EHT MWL Science Working Group 2021), hardly achievable in a single night. The choice of a daily binned light curve would lead to large statistical uncertainties in the data points, and hence, an uncertain definition of the source states. On the other hand, the choice of very large bin sizes, for instance 60 days, would mix high and low states and result in a flattened light curve. Given the trade-off between reasonable statistical uncertainties in the flux points and a low state defined with the least variable emission in it, we chose the compromise of having a bin size equal to 30 days. The published flares from 2005 (Aharonian et al. 2006a), 2008 (Acciari et al. 2009), and 2010 (Abramowski et al. 2012) are well visible in the 30-day-binned light curve as shown in Fig. 1. Furthermore, we also produced the 1-, 7-, and 15-day-binned light curves and their Bayesian blocks as shown in Fig. A.1. In fact, the blocks with elevated flux in these light curves mostly coincided with the elevated blocks from the 30-day-binned light curve (Fig 1). Nevertheless, very long blocks, for instance blocks 5 and 9 in the one-day light curve and blocks 4 and 8 in the 7-day light curve, show that the statistical uncertainties of the flux points are still too high for new flux levels to emerge from the Bayesian block analysis. More importantly, the VHE gamma-ray flare from 2008 falls within blocks below the average in the light curves of Fig. A.1. Therefore, the light curves with bin sizes up to 15-day long have shown to be insufficient for a reliable Bayesian block analysis with M 87 data.

## Appendix B: The interplay between CR pressure and gamma-ray emission

Following the formalism from Jacob & Pfrommer (2017b), a single beta profile function is utilized to parametrize the electron distribution as seen by the Chandra X-Ray Observatory (Weisskopf et al. 2000) in the inner Virgo Cluster:

$$n_e(r) = n_0 \left[ 1 + (r/r_c)^2 \right]^{-3\beta/2}, \qquad (B.1)$$

where $n_0=0.230 \text{ cm}^{-3}$, $\beta=0.29$ and $r_c=0.6 \text{ kpc}$. The equation is valid for $r<44 \text{ kpc}$. For a fully ionized ICM composed of a hydrogen mass fraction of X=0.7 and helium mass fraction of Y=0.28, the ICM density distribution is given by

$$n_{\text{ICM}}(r) = \frac{\mu_e}{\mu} n_e(r), \qquad (B.2)$$

where $\mu_e=1.18$ and $\mu=0.62$ are the mean molecular weight per electron and per particle in the ICM, respectively.

The CRp energy and spatial distributions can be described as a function of the distance from the cluster center, r, and the CRp dimensionless momentum, $p_p$:

$$n_p(r, p_p) = C_p(r) p_p^{-\alpha_p} H(p_p - q_p), \qquad (B.3)$$

where $p_p = P_p/(m_p c)$, $P_p$ is the CRp momentum, $m_p$ the proton mass and c the speed of light, $C_p(r)$ is the spatial distribution of CRp, H is the Heaviside step function, $q_p = 0.5$ is the dimensionless lower momentum cut-off and $\alpha_p$ is the spectral index of the CRp distribution.

The gamma-ray source density distribution $s_\gamma(E_\gamma, r)$ produced locally by the neutral pion decay can be derived based on the $n_{\text{ICM}}(r)$ and $n_p(r)$ distributions. The integration of $s_\gamma(E_\gamma, r)$ from $E_1=300 \text{ GeV}$ to $E_2=\infty$, considering the p-p cross-section
## Appendix A: The long-term light curve of M 87

<figure><img src="image_6.png" /><figcaption>**Fig. A.1.** M 87 long-term light curves with their Bayesian blocks derived according to Sect. 2.2 for 1-day- (top), 7-day- (middle), and 15-day-binned flux points (bottom). The H.E.S.S. data points are shown in black, and a constant fit is shown with the solid orange line. The blocks with their flux levels are given by the blue rectangles. Dashed orange lines indicate a change in the flux level. The orange numbers at the bottom indicate the labels of the blocks.</figcaption></figure>

M 87 is known to be highly variable in VHE gamma-rays, with flares on the order of one day (Abramowski et al. 2012). Therefore, a daily binned light curve could well identify the flares and isolate the low-state periods. Nevertheless, M 87 is a rather weak source in the TeV regime at low states. In fact, at least 5 hours of

observations are needed for a significant source detection (see The EHT MWL Science Working Group 2021), hardly achievable in a single night. The choice of a daily binned light curve would lead to large statistical uncertainties in the data points, and hence, an uncertain definition of the source states. On the other hand, the choice of very large bin sizes, for instance 60 days, would mix high and low states and result in a flattened light curve. Given the trade-off between reasonable statistical uncertainties in the flux points and a low state defined with the least variable emission in it, we chose the compromise of having a bin size equal to 30 days. The published flares from 2005 (Aharonian et al. 2006a), 2008 (Acciari et al. 2009), and 2010 (Abramowski et al. 2012) are well visible in the 30-day-binned light curve as shown in Fig. 1. Furthermore, we also produced the 1-, 7-, and 15-day-binned light curves and their Bayesian blocks as shown in Fig. A.1. In fact, the blocks with elevated flux in these light curves mostly coincided with the elevated blocks from the 30-day-binned light curve (Fig 1). Nevertheless, very long blocks, for instance blocks 5 and 9 in the one-day light curve and blocks 4 and 8 in the 7-day light curve, show that the statistical uncertainties of the flux points are still too high for new flux levels to emerge from the Bayesian block analysis. More importantly, the VHE gamma-ray flare from 2008 falls within blocks below the average in the light curves of Fig. A.1. Therefore, the light curves with bin sizes up to 15-day long have shown to be insufficient for a reliable Bayesian block analysis with M 87 data.

## Appendix B: The interplay between CR pressure and gamma-ray emission

Following the formalism from Jacob & Pfrommer (2017b), a single beta profile function is utilized to parametrize the electron distribution as seen by the Chandra X-Ray Observatory (Weisskopf et al. 2000) in the inner Virgo Cluster:

$$n_e(r) = n_0 \left[ 1 + (r/r_c)^2 \right]^{-3\beta/2}, \qquad (B.1)$$

where $n_0=0.230 \text{ cm}^{-3}$, $\beta=0.29$ and $r_c=0.6 \text{ kpc}$. The equation is valid for $r<44 \text{ kpc}$. For a fully ionized ICM composed of a hydrogen mass fraction of X=0.7 and helium mass fraction of Y=0.28, the ICM density distribution is given by

$$n_{\text{ICM}}(r) = \frac{\mu_e}{\mu} n_e(r), \qquad (B.2)$$

where $\mu_e=1.18$ and $\mu=0.62$ are the mean molecular weight per electron and per particle in the ICM, respectively.

The CRp energy and spatial distributions can be described as a function of the distance from the cluster center, r, and the CRp dimensionless momentum, $p_p$:

$$n_p(r, p_p) = C_p(r) p_p^{-\alpha_p} H(p_p - q_p), \qquad (B.3)$$

where $p_p = P_p/(m_p c)$, $P_p$ is the CRp momentum, $m_p$ the proton mass and c the speed of light, $C_p(r)$ is the spatial distribution of CRp, H is the Heaviside step function, $q_p = 0.5$ is the dimensionless lower momentum cut-off and $\alpha_p$ is the spectral index of the CRp distribution.

The gamma-ray source density distribution $s_\gamma(E_\gamma, r)$ produced locally by the neutral pion decay can be derived based on the $n_{\text{ICM}}(r)$ and $n_p(r)$ distributions. The integration of $s_\gamma(E_\gamma, r)$ from $E_1=300 \text{ GeV}$ to $E_2=\infty$, considering the p-p cross-section
$\sigma_{pp}$ is given by

$$ \lambda_{\gamma}(r) = \int_{300\text{GeV}}^{\infty} dE_{\gamma}s_{\gamma}(E_{\gamma}, r) = \frac{4C_{\text{p}}(r) \frac{m_{\pi^0} c \sigma_{\text{pp}} n_{\text{ICM}}(r)}{3\alpha_{\text{p}}\delta_{\gamma}} \left(\frac{m_{\text{p}}}{2m_{\pi^0}}\right)^{\alpha_{\text{p}}} \left[B_{\chi}\left(\frac{\alpha_{\text{p}}+1}{2\delta_{\gamma}}, \frac{\alpha_{\text{p}}-1}{2\delta_{\gamma}}\right)\right]_{\chi_1}^{\chi_2}}{\frac{m_{\pi^0} c \sigma_{\text{pp}} n_{\text{ICM}}(r)}{3\alpha_{\text{p}}\delta_{\gamma}} \left(\frac{m_{\text{p}}}{2m_{\pi^0}}\right)^{\alpha_{\text{p}}} \left[B_{\chi}\left(\frac{\alpha_{\text{p}}+1}{2\delta_{\gamma}}, \frac{\alpha_{\text{p}}-1}{2\delta_{\gamma}}\right)\right]_{\chi_1}^{\chi_2}}, \quad (B.4) $$

in units of cm⁻³s⁻¹, where $\sigma_{pp} = 3.2 \cdot 10^{-26} (0.96 + e^{4.4-2.4\alpha_{\text{p}}})$ cm⁻² (Pfrommer & Enßlin 2004), $m_{\pi^0}$ is the neutral pion mass, $B_{\chi}(a,b)$ is the incomplete beta function, $\chi_i$ is given by

$$ \chi_i = \left[ 1 + \left( \frac{m_{\pi^0} c^2}{2E_{\gamma,i}} \right)^{2\delta_{\gamma}} \right]^{-1}, \quad (B.5) $$

where $E_{\gamma}$ is the gamma-ray energy calculated at $i=E_1$ and $i=E_2$, $\delta_{\gamma} \approx 0.14\alpha_{\text{p}}^{-1.6} + 0.44$ is the shape factor, and $B_{\chi}[(a, b)]_{\chi_1}^{\chi_2} = B_{\chi_2}(a, b) - B_{\chi_1}(a, b)$.

The uncertainty in the particle interaction model, specifically in the parametrization of the cross-section $\sigma_{pp}$, results in an uncertainty of approximately 50% in the energy-integrated gamma-ray source function for energies ranging from 300 GeV to tens of TeV. Additionally, the systematic uncertainty associated with the detected gamma-ray flux with H.E.S.S. is approximately 20% (Aharonian et al. 2006b). As a result, estimates on the CRp pressure and total energy in CRs should be considered order-of-magnitude estimates. As for the gamma-ray flux estimates, the CRp pressure, $P_{\text{CR}}(r)$, depends on the CRp spatial distribution, $C_{\text{p}}(r)$, and on the proton spectral index, $\alpha_{\text{p}}$:

$$ P_{\text{CR}}(r) = \frac{1}{6} m_{\text{p}} c^2 C_{\text{p}}(r) \left[ \mathcal{B} \left( \frac{1}{1+q_{\text{p}}^2} \right) \left( \frac{\alpha_{\text{p}}-2}{2}, \frac{3-\alpha_{\text{p}}}{2} \right) \right]. \quad (B.6) $$

The CR pressure distribution can be also represented as an energy density distribution, considering an effective adiabatic index for the CRs (fully relativistic value) of $\gamma_{\text{CR}}=4/3$ (Jacob & Pfrommer 2017a):

$$ P_{\text{CR}}(r) = (\gamma_{\text{CR}} - 1)e_{\text{CR}}(r). \quad (B.7) $$

The integral of $e_{\text{CR}}(r)$ in the volume around the source gives the total energy in CRp in units of erg:

$$ U_{\text{CR}} = \int_V e_{\text{CR}}(r) dV. \quad (B.8) $$

The representation of the CR pressure $P_{\text{CR}}(r)$ in terms of the thermal pressure $P_{\text{th.}}(r)$ is useful to characterize the steady state model, in which the heating of the ICM counterbalances the cooling:

$$ X_{\text{CR}}(r) = \frac{P_{\text{CR}}(r)}{P_{\text{th.}}(r)}, \quad (B.9) $$

The thermal pressure $P_{\text{th.}}(r)$ is given by

$$ P_{\text{th.}}(r) = \frac{\mu_e}{\mu} n_e(r) k_B T(r), \quad (B.10) $$

where $k_B$ is the Boltzmann constant and T(r) is the temperature profile of the X-ray emitting electrons. Like the density distribution of the electrons in the plasma, $n_e(r)$, the T(r) is

also parametrized from X-ray Chandra data (Jacob & Pfrommer 2017a):

$$ T(r) = T_0 + (T_1 - T_0) \left[ 1 + \left( \frac{r}{r_T} \right)^{-\eta} \right]^{-1} \left[ 1 + \left( \frac{r}{a r_{200}} \right)^2 \right]^{-0.32}, \quad (B.11) $$

where $T_0=1.9$ keV, $T_1=3.1$ keV, $r_i=28$ kpc, $\eta=1.4$, $a=0.2$ and $r_{200}=1.08$ Mpc.

# Appendix C: The diffuse emission templates

## C.1. The magnetic confinement template

We consider the LOw Frequency Array (LOFAR; Van Haarlem et al. 2013) study of M 87. The 140 MHz radio emission (De Gasperin et al. 2012) traces relativistic electrons, which emit synchrotron photons giving rise to the ≈30 kpc micro-halos (Jacob & Pfrommer 2017a). The radio micro-halos are very well confined within boundaries that have the same dimensions at all radio frequencies down to 25 MHz (De Gasperin et al. 2012). This indicates that the distribution of nonthermal electrons is energy independent and that they are magnetically confined in the lobes. CRp accelerated in the central AGN populate the cluster and could also be present up to the same boundaries as seen in the radio band. Apart from X-ray-depleted regions (cavities) in radio galaxies, where the composition of the material is yet unknown (Abdulla et al. 2019), the ICM at the inner Virgo cluster has a rather high density of 0.1-0.01 cm⁻³ (Jacob & Pfrommer 2017a). Therefore, assuming that the CRp mix with the ICM, the CRp could hadronically interact with it and produce pions. While the charged pions decay to electrons and positrons that likely contribute to part of the radio synchrotron emission from the micro halo, the neutral pions decay to gamma rays. The morphology of the hadronic gamma-ray emission depends not only on the ICM density distribution but also on the CRp energy and spatial distributions (Appendix B). A complete model of this emission would demand a large number of assumptions, for instance, on the poorly known content of the X-ray cavities and the distribution of CRp in the cluster. Therefore, for the sake of simplicity, we assume in this approach that the gamma rays follow the same spatial distribution as the radio emission detected by LOFAR at 140 MHz. This allows us to probe the contribution of an extended hadronic component to M 87 low-state emission despite the different radiation mechanisms. Based on this assumption, we generated a 2D template for the gamma-ray emission, the magnetic confinement template, shown in Fig. C.1 as it would be seen with H.E.S.S. Naturally, we do not expect that the magnetic confinement template explains the entirety of M 87 low-state VHE emission because the template does not account for the emission from the AGN and this study has not measured an extension in the Gaussian model (Sect. 3.2).

Although we utilized the template shown in Fig. C.1 for the main analysis in Sect. 4, we also tested an alternative template. This alternative template attempted to reduce the contribution of the kpc-jet to the overall shape of the emission by masking the central bin in Fig. C.1 (top) before convolving it with the H.E.S.S. PSF. We found that the morphology fit of the alternative template showed a maximum allowed contribution of 30% to the total low-state gamma-ray emission of M87 (≤4.5×10⁻¹³ cm⁻²s⁻¹), compared to 45% when using the original template. However, we opted to proceed with the original approach since it yields a more conservative UL on the CR pressure.
$\sigma_{pp}$ is given by

$$ \lambda_{\gamma}(r) = \int_{300\text{GeV}}^{\infty} dE_{\gamma}s_{\gamma}(E_{\gamma}, r) = \frac{4C_{\text{p}}(r) \frac{m_{\pi^0} c \sigma_{\text{pp}} n_{\text{ICM}}(r)}{3\alpha_{\text{p}}\delta_{\gamma}} \left(\frac{m_{\text{p}}}{2m_{\pi^0}}\right)^{\alpha_{\text{p}}} \left[B_{\chi}\left(\frac{\alpha_{\text{p}}+1}{2\delta_{\gamma}}, \frac{\alpha_{\text{p}}-1}{2\delta_{\gamma}}\right)\right]_{\chi_1}^{\chi_2}}{\frac{m_{\pi^0} c \sigma_{\text{pp}} n_{\text{ICM}}(r)}{3\alpha_{\text{p}}\delta_{\gamma}} \left(\frac{m_{\text{p}}}{2m_{\pi^0}}\right)^{\alpha_{\text{p}}} \left[B_{\chi}\left(\frac{\alpha_{\text{p}}+1}{2\delta_{\gamma}}, \frac{\alpha_{\text{p}}-1}{2\delta_{\gamma}}\right)\right]_{\chi_1}^{\chi_2}}, \quad (B.4) $$

in units of cm⁻³s⁻¹, where $\sigma_{pp} = 3.2 \cdot 10^{-26}(0.96+e^{4.4-2.4\alpha_{\text{p}}})$ cm⁻² (Pfrommer & Enßlin 2004), $m_{\pi^0}$ is the neutral pion mass, $B_{\chi}(a,b)$ is the incomplete beta function, $\chi_i$ is given by

$$ \chi_i = \left[ 1 + \left( \frac{m_{\pi^0} c^2}{2E_{\gamma,i}} \right)^{2\delta_{\gamma}} \right]^{-1}, \quad (B.5) $$

where $E_{\gamma}$ is the gamma-ray energy calculated at $i=E_1$ and $i=E_2$, $\delta_{\gamma} \approx 0.14\alpha_p^{-1.6} + 0.44$ is the shape factor, and $B_{\chi}[(a, b)]_{\chi_1}^{\chi_2} = B_{\chi_2}(a, b) - B_{\chi_1}(a, b)$.

The uncertainty in the particle interaction model, specifically in the parametrization of the cross-section $\sigma_{pp}$, results in an uncertainty of approximately 50% in the energy-integrated gamma-ray source function for energies ranging from 300 GeV to tens of TeV. Additionally, the systematic uncertainty associated with the detected gamma-ray flux with H.E.S.S. is approximately 20% (Aharonian et al. 2006b). As a result, estimates on the CRp pressure and total energy in CRs should be considered order-of-magnitude estimates. As for the gamma-ray flux estimates, the CRp pressure, $P_{\text{CR}}(r)$, depends on the CRp spatial distribution, $C_{\text{p}}(r)$, and on the proton spectral index, $\alpha_{\text{p}}$:

$$ P_{\text{CR}}(r) = \frac{1}{6} m_{\text{p}} c^2 C_{\text{p}}(r) \left[ \mathcal{B} \left( \frac{1}{1+q_{\text{p}}^2} \right) \left( \frac{\alpha_{\text{p}}-2}{2}, \frac{3-\alpha_{\text{p}}}{2} \right) \right]. \quad (B.6) $$

The CR pressure distribution can be also represented as an energy density distribution, considering an effective adiabatic index for the CRs (fully relativistic value) of $\gamma_{\text{CR}}=4/3$ (Jacob & Pfrommer 2017a):

$$ P_{\text{CR}}(r) = (\gamma_{\text{CR}} - 1)e_{\text{CR}}(r). \quad (B.7) $$

The integral of $e_{\text{CR}}(r)$ in the volume around the source gives the total energy in CRp in units of erg:

$$ U_{\text{CR}} = \int_V e_{\text{CR}}(r) dV. \quad (B.8) $$

The representation of the CR pressure $P_{\text{CR}}(r)$ in terms of the thermal pressure $P_{\text{th.}}(r)$ is useful to characterize the steady state model, in which the heating of the ICM counterbalances the cooling:

$$ X_{\text{CR}}(r) = \frac{P_{\text{CR}}(r)}{P_{\text{th.}}(r)}, \quad (B.9) $$

The thermal pressure $P_{\text{th.}}(r)$ is given by

$$ P_{\text{th.}}(r) = \frac{\mu_e}{\mu} n_e(r) k_B T(r), \quad (B.10) $$

where $k_B$ is the Boltzmann constant and $T(r)$ is the temperature profile of the X-ray emitting electrons. Like the density distribution of the electrons in the plasma, $n_e(r)$, the $T(r)$ is

also parametrized from X-ray Chandra data (Jacob & Pfrommer 2017a):

$$ T(r) = T_0 + (T_1 - T_0) \left[ 1 + \left( \frac{r}{r_T} \right)^{-\eta} \right]^{-1} \left[ 1 + \left( \frac{r}{a r_{200}} \right)^2 \right]^{-0.32}, \quad (B.11) $$

where $T_0=1.9$ keV, $T_1=3.1$ keV, $r_i=28$ kpc, $\eta=1.4$, $a=0.2$ and $r_{200}=1.08$ Mpc.

# Appendix C: The diffuse emission templates

## C.1. The magnetic confinement template

We consider the LOw Frequency Array (LOFAR; Van Haarlem et al. 2013) study of M 87. The 140 MHz radio emission (De Gasperin et al. 2012) traces relativistic electrons, which emit synchrotron photons giving rise to the ≈30 kpc micro-halos (Jacob & Pfrommer 2017a). The radio micro-halos are very well confined within boundaries that have the same dimensions at all radio frequencies down to 25 MHz (De Gasperin et al. 2012). This indicates that the distribution of nonthermal electrons is energy independent and that they are magnetically confined in the lobes. CRp accelerated in the central AGN populate the cluster and could also be present up to the same boundaries as seen in the radio band. Apart from X-ray-depleted regions (cavities) in radio galaxies, where the composition of the material is yet unknown (Abdulla et al. 2019), the ICM at the inner Virgo cluster has a rather high density of 0.1-0.01 cm⁻³ (Jacob & Pfrommer 2017a). Therefore, assuming that the CRp mix with the ICM, the CRp could hadronically interact with it and produce pions. While the charged pions decay to electrons and positrons that likely contribute to part of the radio synchrotron emission from the micro halo, the neutral pions decay to gamma rays. The morphology of the hadronic gamma-ray emission depends not only on the ICM density distribution but also on the CRp energy and spatial distributions (Appendix B). A complete model of this emission would demand a large number of assumptions, for instance, on the poorly known content of the X-ray cavities and the distribution of CRp in the cluster. Therefore, for the sake of simplicity, we assume in this approach that the gamma rays follow the same spatial distribution as the radio emission detected by LOFAR at 140 MHz. This allows us to probe the contribution of an extended hadronic component to M 87 low-state emission despite the different radiation mechanisms. Based on this assumption, we generated a 2D template for the gamma-ray emission, the magnetic confinement template, shown in Fig. C.1 as it would be seen with H.E.S.S. Naturally, we do not expect that the magnetic confinement template explains the entirety of M 87 low-state VHE emission because the template does not account for the emission from the AGN and this study has not measured an extension in the Gaussian model (Sect. 3.2).

Although we utilized the template shown in Fig. C.1 for the main analysis in Sect. 4, we also tested an alternative template. This alternative template attempted to reduce the contribution of the kpc-jet to the overall shape of the emission by masking the central bin in Fig. C.1 (top) before convolving it with the H.E.S.S. PSF. We found that the morphology fit of the alternative template showed a maximum allowed contribution of 30% to the total low-state gamma-ray emission of M87 (≤4.5×10⁻¹³ cm⁻²s⁻¹), compared to 45% when using the original template. However, we opted to proceed with the original approach since it yields a more conservative UL on the CR pressure.
<figure><img src="image_6.png" /><figcaption>**Fig. C.1.** 2D templates for the gamma-ray emission in the inner Virgo Cluster, following the square root of the radio emission intensity normalized at its maximum as seen by LOFAR (Fig. 5 from De Gasperin et al. 2012) for better visualization: the model displayed with a bin size of 36" (top) and the intrinsic emission convolved with the H.E.S.S. PSF (bottom).</figcaption></figure>

<caption>**Fig. C.2.** 2D templates for the gamma-ray emission in the inner Virgo Cluster, following the square root of the gamma-ray intensity due to $\pi^0$ decay considering the angular resolution of the H.E.S.S. sky maps: the model displayed with a bin size of 36" (top) and the intrinsic emission convolved with the H.E.S.S. PSF (bottom).</caption>

## C.2. The steady-state template

In the second approach, we adopted the model from Jacob & Pfrommer (2017a,b) for the AGN feedback. The authors propose a steady-state model, in which the CRp pressure is sufficient to counterbalance, at every distance from the cluster center, the radiative cooling of the ICM. The CRp excite Alfvén waves through the streaming instability and the nonlinear Landau damping of these waves provides an efficient mechanism for heating the ICM. In addition to the CR pressure, the thermal pressure helps halt the cooling flow toward the cluster center. The model solves the cooling flow problem and predicts a steady and extended gamma-ray signal due to neutral pion decay. The assumption of a steady state leads to an analytical estimation of

the CRp pressure distribution in the inner cluster. Using the CRp pressure distribution in Fig. A.1 from Jacob & Pfrommer (2017a) and equations from Appendix B we derived the distribution of gamma rays above 300 GeV produced by the decay of neutral pions for the same set of CRp distributions utilized in the first approach of our study. Afterward, we integrated the gamma-ray distribution along the line of sight to generate a second 2D gamma-ray emission template as it would be seen with H.E.S.S. The template is shown in Fig. C.2. We also integrated the gamma-ray distribution in the volume around the source ($\le 20$ kpc) and accounted for the distance to M 87 to obtain a predicted VHE gamma-ray flux of $2.2 \times 10^{-13} \text{cm}^{-2}\text{s}^{-1}$ by the steady-state model with $\alpha_p=2.1$, that is, $\approx 15\%$ of the low state of M 87 VHE emission (Table 2). The predicted flux becomes smaller for steeper indices and is, therefore, less relevant for this study.
<figure><img src="image_7.png" /><figcaption>**Fig. C.1.** 2D templates for the gamma-ray emission in the inner Virgo Cluster, following the square root of the radio emission intensity normalized at its maximum as seen by LOFAR (Fig. 5 from De Gasperin et al. 2012) for better visualization: the model displayed with a bin size of 36" (top) and the intrinsic emission convolved with the H.E.S.S. PSF (bottom).</figcaption></figure>

<caption>**Fig. C.2.** 2D templates for the gamma-ray emission in the inner Virgo Cluster, following the square root of the gamma-ray intensity due to $\pi^0$ decay considering the angular resolution of the H.E.S.S. sky maps: the model displayed with a bin size of 36" (top) and the intrinsic emission convolved with the H.E.S.S. PSF (bottom).</caption>

## C.2. The steady-state template

In the second approach, we adopted the model from Jacob & Pfrommer (2017a,b) for the AGN feedback. The authors propose a steady-state model, in which the CRp pressure is sufficient to counterbalance, at every distance from the cluster center, the radiative cooling of the ICM. The CRp excite Alfvén waves through the streaming instability and the nonlinear Landau damping of these waves provides an efficient mechanism for heating the ICM. In addition to the CR pressure, the thermal pressure helps halt the cooling flow toward the cluster center. The model solves the cooling flow problem and predicts a steady and extended gamma-ray signal due to neutral pion decay. The assumption of a steady state leads to an analytical estimation of

the CRp pressure distribution in the inner cluster. Using the CRp pressure distribution in Fig. A.1 from Jacob & Pfrommer (2017a) and equations from Appendix B we derived the distribution of gamma rays above 300 GeV produced by the decay of neutral pions for the same set of CRp distributions utilized in the first approach of our study. Afterward, we integrated the gamma-ray distribution along the line of sight to generate a second 2D gamma-ray emission template as it would be seen with H.E.S.S. The template is shown in Fig. C.2. We also integrated the gamma-ray distribution in the volume around the source ($\le 20$ kpc) and accounted for the distance to M 87 to obtain a predicted VHE gamma-ray flux of $2.2 \times 10^{-13} \text{cm}^{-2}\text{s}^{-1}$ by the steady-state model with $\alpha_p=2.1$, that is, $\approx 15\%$ of the low state of M 87 VHE emission (Table 2). The predicted flux becomes smaller for steeper indices and is, therefore, less relevant for this study.