Open Access
Issue
A&A
Volume 675, July 2023
Article Number A138
Number of page(s) 13
Section Extragalactic astronomy
DOI https://doi.org/10.1051/0004-6361/202346056
Published online 11 July 2023

© The Authors 2023

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.

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 ≈2.3° from its center (r500 ≈ 662.6 kpc; Simionescu et al. 2017). It is known as a cool core (CC) cluster, meaning its central region (≲0.01r180; Peterson & Fabian 2006; Urban et al. 2011) is filled with a plasma that is colder and denser than the surrounding 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 ∼1000M 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: 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), CT 5 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 θ ≤ 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σ 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 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.

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

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 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 rotationally symmetric 2D Gaussian model given by

(1)

where x and y are variables, x0 and y0 are the center of the function, 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 M 87 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 data1.

A set of systematic checks confirms the stability of the results toward different configurations2. We analyzed the M 87 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 σ) 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 M 87 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.

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 M 87 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−12 cm−2 s−1.

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.

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 distribution3 hardens from the low to the high state, consistent with previous results (Aharonian et al. 2006a; MAGIC Collaboration 2020).

Table 2.

Results of the H.E.S.S. analysis for the low, intermediate, and high flux states.

3.2. The morphology of the low state of M 87

We performed a morphology fit to the M 87 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 σG extension of 58″ is derived at the 99.7% confidence level (c.l.). This translates to a physical extension limit of ≈4.6 kpc, assuming a distance to M 87 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).

thumbnail 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 σG = 18″ for better visualization. The best-fit position is shown by the blue marker for the Gaussian model, with its σ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σ 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σ statistical uncertainties for the respective source states. The dashed blue circle shows the resulting 99.7% c.l. of the Gaussian σ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.

3.3. The origin of the gamma-ray emission of M 87

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 ≈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, the best-fit position of the point-like model is consistent with the radio core.

Table 3.

Best-fit parameters of the morphology fit.

Our extension UL on the low state of M 87 excludes the radio lobes (≈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 (≲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 (≲10rg ≈ 0.003 pc, where rg 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.

thumbnail 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 ≳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, R1/2 ≈ 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 (∼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).

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.

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 ≈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 ∼10 Myr in a ∼10 μG magnetic field and hence orders of magnitude shorter than the pp 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 ≈1 kpc (≈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 (≈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 ≳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 M 87 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 (≳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 (≈1.7 Mpc, that is, ≈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 ≈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 ∼1018 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 ≈1.9° around the virial radius (≈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 M 87 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 M 87, we consider the neutral pion decay gamma-ray emission from pp 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 (≈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σ uncertainty on the radial X-ray brightness of approximately 85% in the region close to the M 87 core (within 4 kpc) and approximately 30–40% in the region up to 0.5° from the M 87 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 αp, which we vary from 2.1 to 2.6. PL distributions with spectral indices αp ≳ 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 M 87 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 M 87 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 M 87 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σ preference for the point-like model in comparison to the hybrid model. The results of our analysis show that the 3σ 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 ≲6.7 × 10−13 cm−2s−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 ≥55% of the M 87 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 (XCR(r) = PCR(r)/Pth.(r)) is constant in the inner region of the cluster. A constant XCR 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 XCR ≲ 0.17 for a proton distribution with α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 XCR obtained in this study is underestimated, at least in the region of the cavities.

The CRp energy density eCR was estimated and integrated within the volume to yield the total energy in CRp in the inner 20 kpc of the Virgo Cluster of UCR ≲ 5 × 1058 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 M 87 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 ≈20 kpc after ≈15 Myr. Therefore, for the streaming CRp to reach UCR  ≲  5  ×  1058 erg in 15 Myr, an average jet power of Pj ≲ 6 × 1043 erg s−1 is necessary, considering the extreme case of 100% efficiency in the CRp acceleration. Even though the estimates on UCR and Pj 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 UCR derived in this work is twice as constraining as the result from Nizamov & Pshirkov (2022), which placed an UL on UCR up to 35 kpc from the M 87 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 α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 M 87 low-state emission would be necessary to account for the overall gamma-ray flux estimated with H.E.S.S. (αγ = 2.63 ± 0.09; Table 2), possibly, the kiloparsec jet as in Centaurus A (Abdalla et al. 2020a). For steeper CRp spectral indices our ULs on XCR and UCR become less constraining as seen in Fig. 4 by the blue ULs.

thumbnail Fig. 4.

UL of the maximum CR to thermal pressure ratio, XCR, max. = max(PCR(r)/Pth.(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.

4.3. The steady-state model

The fit of the M 87 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 the diffuse component, we derive a 99.7% c.l. UL on the total flux of the diffuse component of ≲8.1 × 10−13 cm−2 s−1, that is, ≲55% of the low state of M 87 VHE gamma-ray flux (Table 2). This UL is a factor of ≈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 ≈55% of the flux of M 87 low state. Afterward, we derived an UL on the maximum CR to thermal pressure ratio (XCR, 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 αp = 2.1, XCR, max. ≲ 0.32 and the total energy in CRp in the inner 20 kpc of the Virgo Cluster is also constrained to UCR ≲ 5 × 1058 erg. While the steady-state model predicts a level of XCR, max. ≈ 0.10 (Fig. A1 in Jacob & Pfrommer 2017a), the H.E.S.S. UL is XCR, max. ≲ 0.32 for α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. XCR, max. UL derived from the morphology fit of the M 87 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 M 87, probing an extended emission in its low state and testing its hadronic origin. Neutral pions are produced in pp 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 M 87 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 σG of a rotationally symmetric 2D Gaussian model of 58″ (≈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σ 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 (≈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 (R1/2 ≈ 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 ≲45% of the VHE gamma-ray flux detected from the low state of M 87 for the magnetic confinement hybrid model and ≲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, α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 (XCR, max.) and the maximum energy in CRp (UCR) in the region up to 20 kpc from the cluster center. For a CRp distribution with spectral index αp = 2.1, XCR, max. ≲ 0.17 for the hybrid models with the magnetic confinement template, while XCR, max. ≲ 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 αp = 2.1, is constrained at the 99.7% c.l. to UCR  ≲  5  ×  1058 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 UCR is also twice as constraining as the UL based on Fermi-LAT data (Nizamov & Pshirkov 2022).


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

2

More details in the Ph.D. Thesis Barbosa Martins (2022).

3

dϕ/dE = ϕ0(E/TeV)−Γ, where Γ is the spectral index, E the gamma-ray energy and ϕ0 the normalization at 1 TeV.

Acknowledgments

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

  1. Abdalla, H., Adam, R., Aharonian, F., et al. 2020a, Nature, 582, 356 [NASA ADS] [CrossRef] [Google Scholar]
  2. Abdalla, H., Aharonian, F., Ait Benkhali, F., et al. 2020b, Nat. Astron., 4, 167 [NASA ADS] [CrossRef] [Google Scholar]
  3. Abdollahi, S., Acero, F., Ackermann, M., et al. 2020, ApJS, 247, 33 [Google Scholar]
  4. Abdulla, Z., Carlstrom, J. E., Mantz, A. B., et al. 2019, ApJ, 871, 195 [NASA ADS] [CrossRef] [Google Scholar]
  5. Abramowski, A., Acero, F., Aharonian, F., et al. 2012, ApJ, 746, 151 [NASA ADS] [CrossRef] [Google Scholar]
  6. Acciari, V. A., Beilicke, M., Blaylock, G., et al. 2008, ApJ, 679, 397 [NASA ADS] [CrossRef] [Google Scholar]
  7. Acciari, V. A., Aliu, E., Arlen, T., et al. 2009, Science, 325, 444 [Google Scholar]
  8. Acero, F., Aharonian, F., Akhperjanian, A. G., et al. 2010, MNRAS, 402, 1877 [NASA ADS] [CrossRef] [Google Scholar]
  9. Ackermann, M., Ajello, M., Albert, A., et al. 2014, ApJ, 787, 18 [Google Scholar]
  10. Ackermann, M., Ajello, M., Albert, A., et al. 2015, ApJ, 812, 159 [NASA ADS] [CrossRef] [Google Scholar]
  11. Ackermann, M., Ajello, M., Baldini, L., et al. 2016, ApJ, 826, 1 [NASA ADS] [CrossRef] [Google Scholar]
  12. Adam, R., Goksu, H., Brown, S., Rudnick, L., & Ferrari, C. 2021, A&A, 648, A60 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  13. Aharonian, F., Akhperjanian, A. G., Bazer-Bachi, A. R., et al. 2006a, Science, 314, 1424 [Google Scholar]
  14. Aharonian, F., Akhperjanian, A. G., Bazer-Bachi, A. R., et al. 2006b, A&A, 457, 899 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  15. Aharonian, F., Akhperjanian, A. G., Anton, G., et al. 2009a, A&A, 495, 27 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  16. Aharonian, F., Akhperjanian, A. G., Anton, G., et al. 2009b, A&A, 502, 437 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  17. Ahnen, M. L., Ansoldi, S., Antonelli, L. A., et al. 2016, A&A, 593, A91 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  18. Aleksić, J., Alvarez, E. A., Antonelli, L. A., et al. 2012, A&A, 541, A99 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  19. Alfaro, R., Alvarez, C., Arteaga-Velázquez, J. C., et al. 2022, ApJ, 934, 158 [NASA ADS] [CrossRef] [Google Scholar]
  20. Anderhub, H., Antonelli, L. A., Antoranz, P., et al. 2010, ApJ, 710, 828 [NASA ADS] [CrossRef] [Google Scholar]
  21. Arlen, T., Aune, T., Beilicke, M., et al. 2012, ApJ, 757, 123 [Google Scholar]
  22. Ashton, T., Backes, M., Balzer, A., et al. 2020, Astropart. Phys., 118, 102425 [NASA ADS] [CrossRef] [Google Scholar]
  23. Astropy Collaboration (Robitaille, T. P., et al.) 2013, A&A, 558, A33 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  24. Astropy Collaboration (Price-Whelan, A. M., et al.) 2018, AJ, 156, 123 [Google Scholar]
  25. Baghmanyan, V., Zargaryan, D., Aharonian, F., et al. 2022, MNRAS, 516, 562 [NASA ADS] [CrossRef] [Google Scholar]
  26. Barbiellini, G., Bastieri, D., Bechtol, K., et al. 2014, ApJ, 784, 118 [NASA ADS] [CrossRef] [Google Scholar]
  27. Barbosa Martins, V. 2022, Ph.D. Thesis, Humboldt-Universität zu Berlin, Germany, https://edoc.hu-berlin.de/handle/18452/25715 [Google Scholar]
  28. Becker, H. B., White, R. L., & Helfand, D. J. 1995, ApJ, 450, 559 [NASA ADS] [CrossRef] [Google Scholar]
  29. Berge, D., Funk, S., & Hinton, J. 2007, A&A, 466, 1219 [CrossRef] [EDP Sciences] [Google Scholar]
  30. Bird, S., Harris, W. E., Blakeslee, J. P., & Flynn, C. 2010, A&A, 524, A71 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  31. Blakeslee, J. P., Jordn, A., Mei, S., et al. 2009, ApJ, 694, 556 [NASA ADS] [CrossRef] [Google Scholar]
  32. Boughelilba, M., Reimer, A., & Merten, L. 2022, ApJ, 938, 79 [NASA ADS] [CrossRef] [Google Scholar]
  33. Brüggen, M., Kaiser, C. R., Churazov, E., & Enßlin, T. A. 2002, MNRAS, 331, 545 [CrossRef] [Google Scholar]
  34. Casanova, S., & Dingus, B. L. 2008, Astropart. Phys., 29, 63 [NASA ADS] [CrossRef] [Google Scholar]
  35. Cash, W. 1979, ApJ, 228, 939 [Google Scholar]
  36. Churazov, E., Bruggen, M., Kaiser, C. R., Bohringer, H., & Forman, W. 2001, ApJ, 554, 261 [NASA ADS] [CrossRef] [Google Scholar]
  37. Churazov, E., Sunyaev, R., Forman, W., & Böhringer, H. 2002, MNRAS, 332, 729 [NASA ADS] [CrossRef] [Google Scholar]
  38. De Gasperin, F., Orrú, E., Murgia, M., et al. 2012, A&A, 547, A56 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  39. Deil, C., Lefaucheur, J., Zanin, R., et al. 2017, in Proc. 35th Int. Cosmic RayConf.– PoS(ICRC2017), 35, 766 [Google Scholar]
  40. de Naurois, M., & Rolland, L. 2009, Astropart. Phys., 32, 231 [Google Scholar]
  41. Domainko, W., Nedbal, D., Hinton, J. A., & Martineau-Huynh, O. 2009, Int. J. Mod. Phys. D, 18, 1627 [NASA ADS] [CrossRef] [Google Scholar]
  42. Enßlin, T., Pfrommer, C., Miniati, F., & Subramanian, K. 2011, A&A, 527, 1 [Google Scholar]
  43. Fermi, E. 1949, Phys. Rev., 75, 1169 [NASA ADS] [CrossRef] [Google Scholar]
  44. Forman, W., Churazov, E., Jones, C., et al. 2017, ApJ, 844, 122 [NASA ADS] [CrossRef] [Google Scholar]
  45. Franceschini, A., Foffano, L., Prandini, E., & Tavecchio, F. 2019, A&A, 629, A2 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  46. Freeman, P., Doe, S., & Siemiginowska, A. 2001, SPIE, 4477, 76 [Google Scholar]
  47. Galante, N. 2009, in Proc. 31th Int. Cosmic Ray Conf.– PoS(ICRC2009), 37, 4 [Google Scholar]
  48. Guo, F., & Oh, S. P. 2008, MNRAS, 384, 251 [NASA ADS] [CrossRef] [Google Scholar]
  49. Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357 [Google Scholar]
  50. Hinton, J. A., & Hofmann, W. 2009, ARA&A, 47, 523 [NASA ADS] [CrossRef] [Google Scholar]
  51. Huber, B., Tchernin, C., Eckert, D., et al. 2013, A&A, 560, A64 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  52. Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90 [Google Scholar]
  53. Inoue, S., Aharonian, F. A., & Sugiyama, N. 2005, ApJ, 628, L9 [NASA ADS] [CrossRef] [Google Scholar]
  54. Jacob, S., & Pfrommer, C. 2017a, MNRAS, 467, 1449 [Google Scholar]
  55. Jacob, S., & Pfrommer, C. 2017b, MNRAS, 467, 1478 [NASA ADS] [Google Scholar]
  56. Kashibadze, O. G., Karachentsev, I. D., & Karachentseva, V. E. 2020, A&A, 635, 1 [Google Scholar]
  57. MAGIC Collaboration (Acciari, V. A., et al.) 2020, MNRAS, 492, 5354 [NASA ADS] [CrossRef] [Google Scholar]
  58. Massaro, F., Harris, D. E., Tremblay, G. R., et al. 2013, ApJS, 206, 7 [Google Scholar]
  59. Mathews, W. G., Faltenbacher, A., & Brighenti, F. 2006, ApJ, 638, 659 [CrossRef] [Google Scholar]
  60. Mitchell, A. M., Caroff, S., Hinton, J., et al. 2023, A&A, 395, 18248039 [Google Scholar]
  61. Nigro, C., Deil, C., Zanin, R., et al. 2019, A&A, 625, A10 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  62. Nizamov, B. A., & Pshirkov, M. S. 2022, JETP, 115, 245 [NASA ADS] [Google Scholar]
  63. Ohm, S., van Eldik, C., & Egberts, K. 2009, Astropart. Phys., 31, 383 [NASA ADS] [CrossRef] [Google Scholar]
  64. Owen, F. N., Eilek, J. A., & Kassim, N. E. 2000, ApJ, 543, 611 [NASA ADS] [CrossRef] [Google Scholar]
  65. Parsons, R. D., & Hinton, J. A. 2014, Astropart. Phys., 56, 26 [Google Scholar]
  66. Peterson, J. R., & Fabian, A. C. 2006, Phys. Rep., 427, 1 [Google Scholar]
  67. Pfrommer, C. 2013, ApJ, 779, 10 [NASA ADS] [CrossRef] [Google Scholar]
  68. Pfrommer, C., & Enßlin, T. A. 2004, A&A, 426, 777 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  69. Pinzke, A., Pfrommer, C., & Bergström, L. 2011, Phys. Rev. D, 84, 123509 [NASA ADS] [CrossRef] [Google Scholar]
  70. Piron, F., Djannati-Atai, A., Punch, M., et al. 2001, A&A, 374, 895 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  71. Prokhorov, D. A., & Churazov, E. M. 2014, A&A, 567, A93 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  72. 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 [Google Scholar]
  73. Rieger, F. M., Bosch-Ramon, V., & Duffy, P. 2007, Ap&SS, 309, 119 [NASA ADS] [CrossRef] [Google Scholar]
  74. Ruszkowski, M., Bruggen, M., & Begelman, M. C. 2004, ApJ, 611, 158 [NASA ADS] [CrossRef] [Google Scholar]
  75. Scargle, J. D. 1998, ApJ, 504, 405 [Google Scholar]
  76. Scargle, J. D., Norris, J. P., Jackson, B., & Chiang, J. 2013, ApJ, 764, 167 [Google Scholar]
  77. Simionescu, A., Werner, N., Mantz, A., Allen, S. W., & Urban, O. 2017, MNRAS, 469, 1476 [Google Scholar]
  78. The EHT MWL Science Working Group (Algaba, J. C., et al.) 2021, ApJ, 911, L11 [NASA ADS] [CrossRef] [Google Scholar]
  79. The Event Horizon Telescope Collaboration (Akiyama, K., et al.) 2019, ApJ, 875, L1 [Google Scholar]
  80. The FIRST Project Team 1994, VLA FIRST (1.4 GHz): FIRST, https://skyview.gsfc.nasa.gov/ [Google Scholar]
  81. Thopson, A. R., Clark, B. G., Wade, C. M., & Napier, P. J. 1980, ApJS, 44, 151 [NASA ADS] [CrossRef] [Google Scholar]
  82. Urban, O., Werner, N., Simionescu, A., Allen, S. W., & Böhringer, H. 2011, MNRAS, 414, 2101 [NASA ADS] [CrossRef] [Google Scholar]
  83. Van Haarlem, M. P., Wise, M. W., Gunst, A. W., et al. 2013, A&A, 556, A2 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  84. Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nat. Methods, 17, 261 [Google Scholar]
  85. Weil, M. L., Bland-Hawthorn, J., & Malin, D. F. 1997, ApJ, 490, 664 [Google Scholar]
  86. Weisskopf, M. C., Tananbaum, H. D., Speybroeck, L. P. V., & O’Dell, S. L. 2000, SPIE, 4012, 2 [NASA ADS] [Google Scholar]
  87. Xi, S. Q., Wang, X. Y., Liang, Y. F., et al. 2018, PRD, 98, 063006 [NASA ADS] [CrossRef] [Google Scholar]
  88. Xue, R., Wang, Z. R., & Li, W. J. 2022, Phys. Rev. D, 106, 103021 [NASA ADS] [CrossRef] [Google Scholar]
  89. Young, A. J., Wilson, A. S., & Mundell, C. G. 2002, ApJ, 579, 560 [NASA ADS] [CrossRef] [Google Scholar]
  90. Zakamska, N. L., & Narayan, R. 2003, ApJ, 582, 162 [NASA ADS] [CrossRef] [Google Scholar]
  91. Zhuravleva, I., Churazov, E., Schekochihin, A. A., et al. 2014, Nature, 515, 85 [Google Scholar]

Appendix A: The long-term light curve of M 87

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.

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

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:

(B.1)

where n0 = 0.230 cm−3, β = 0.29 and rc = 0.6 kpc. The equation is valid for r< 44 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

(B.2)

where μe = 1.18 and μ = 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, pp:

(B.3)

where pp = Pp/(mpc), Pp is the CRp momentum, mp the proton mass and c the speed of light, Cp(r) is the spatial distribution of CRp, H is the Heaviside step function, qp = 0.5 is the dimensionless lower momentum cut-off and αp is the spectral index of the CRp distribution.

The gamma-ray source density distribution sγ(Eγ, r) produced locally by the neutral pion decay can be derived based on the nICM(r) and np(r) distributions. The integration of sγ(Eγ, r) from E1 = 300  GeV to E2=∞, considering the p-p cross-section σpp is given by

(B.4)

in units of cm−3s−1, where σpp = 3.2 ⋅ 10−26(0.96 + e4.4 − 2.4αp) cm−2 (Pfrommer & Enßlin 2004), mπ0 is the neutral pion mass, B𝒳(a,b) is the incomplete beta function, 𝒳i is given by

(B.5)

where Eγ is the gamma-ray energy calculated at i=E1 and i=E2, is the shape factor, and .

The uncertainty in the particle interaction model, specifically in the parametrization of the cross-section σ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, PCR(r), depends on the CRp spatial distribution, Cp(r), and on the proton spectral index, αp:

(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 γCR = 4/3 (Jacob & Pfrommer 2017a):

(B.7)

The integral of eCR(r) in the volume around the source gives the total energy in CRp in units of erg:

(B.8)

The representation of the CR pressure PCR(r) in terms of the thermal pressure Pth.(r) is useful to characterize the steady state model, in which the heating of the ICM counterbalances the cooling:

(B.9)

The thermal pressure Pth.(r) is given by

(B.10)

where kB 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, ne(r), the T(r) is also parametrized from X-ray Chandra data (Jacob & Pfrommer 2017a):

(B.11)

where T0 = 1.9 keV, T1 = 3.1 keV, rt = 28 kpc, η = 1.4, a = 0.2 and r200 = 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−3 (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).

thumbnail 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).

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−13 cm−2s−1), 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.

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 (≲20 kpc) and accounted for the distance to M 87 to obtain a predicted VHE gamma-ray flux of 2.2 × 10−13cm−2s−1 by the steady-state model with αp = 2.1, that is, ≈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.

thumbnail 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 π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).

All Tables

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.

Table 2.

Results of the H.E.S.S. analysis for the low, intermediate, and high flux states.

Table 3.

Best-fit parameters of the morphology fit.

All Figures

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

In the text
thumbnail 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 σG = 18″ for better visualization. The best-fit position is shown by the blue marker for the Gaussian model, with its σ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σ 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σ statistical uncertainties for the respective source states. The dashed blue circle shows the resulting 99.7% c.l. of the Gaussian σ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.

In the text
thumbnail 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 ≳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, R1/2 ≈ 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 (∼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).

In the text
thumbnail Fig. 4.

UL of the maximum CR to thermal pressure ratio, XCR, max. = max(PCR(r)/Pth.(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.

In the text
thumbnail 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.

In the text
thumbnail 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).

In the text
thumbnail 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 π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).

In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.