# Curvature in the very-high energy gamma-ray spectrum of M 87*

H.E.S.S. Collaboration: F. Aharonian<sup>1,2,3</sup>, F. Ait Benkhali<sup>5</sup>, J. Aschersleben<sup>6</sup>, H. Ashkar<sup>7</sup>, M. Backes<sup>8,9</sup>, V. Barbosa Martins<sup>10,**</sup>, R. Batzofin<sup>11</sup>, Y. Becherini<sup>12,13</sup>, D. Berge<sup>10,14</sup>, K. Bernlöhr<sup>2</sup>, M. Böttcher<sup>9</sup>, C. Boisson<sup>15</sup>, J. Bolmont<sup>16</sup>, M. de Bony de Lavergne<sup>17</sup>, F. Bradascio<sup>17</sup>, R. Brose<sup>1</sup>, F. Brun<sup>17</sup>, B. Bruno<sup>18</sup>, T. Bulik<sup>19</sup>, C. Burger-Scheidlin<sup>1</sup>, T. Bylund<sup>17</sup>, S. Casanova<sup>20</sup>, R. Cecil<sup>21,**</sup>, J. Celic<sup>18</sup>, M. Cerruti<sup>12</sup>, T. Chand<sup>9</sup>, S. Chandra<sup>9</sup>, A. Chen<sup>22</sup>, J. Chibueze<sup>9</sup>, O. Chibueze<sup>9</sup>, T. Collins<sup>11</sup>, G. Cotter<sup>23</sup>, J. Damascene Mbarubucyeye<sup>10</sup>, A. Djannati-Ataï<sup>12</sup>, J. Djuvsland<sup>2</sup>, A. Dmytriiev<sup>9</sup>, K. Egberts<sup>11</sup>, S. Einecke<sup>24</sup>, S. Fegan<sup>7</sup>, G. Fontaine<sup>7</sup>, S. Funk<sup>18</sup>, S. Gabici<sup>12</sup>, J. F. Glicenstein<sup>17</sup>, J. Glombitza<sup>18</sup>, G. Grolleron<sup>16</sup>, L. Haerer<sup>2</sup>, W. Hofmann<sup>2</sup>, T. L. Holch<sup>10</sup>, M. Holler<sup>25</sup>, D. Horns<sup>21</sup>, M. Jamrozy<sup>26</sup>, F. Jankowsky<sup>5</sup>, V. Joshi<sup>18</sup>, I. Jung-Richardt<sup>18</sup>, E. Kasai<sup>8</sup>, K. Katarzyński<sup>27</sup>, R. Khatoon<sup>9</sup>, B. Khélifi<sup>12</sup>, W. Kluźniak<sup>28</sup>, Nu. Komin<sup>22</sup>, K. Kosack<sup>17</sup>, D. Kostunin<sup>10</sup>, A. Kundu<sup>9</sup>, R. G. Lang<sup>18</sup>, S. Le Stum<sup>29</sup>, F. Leitl<sup>18</sup>, A. Lemière<sup>12</sup>, M. Lemoine-Goumard<sup>30</sup>, J.-P. Lenain<sup>16</sup>, I. Lypova<sup>5,**</sup>, A. Luashvili<sup>15</sup>, J. Mackey<sup>1</sup>, D. Malyshev<sup>31</sup>, G. Martí-Devesa<sup>25</sup>, R. Marx<sup>5</sup>, A. Mehta<sup>10</sup>, M. Meyer<sup>21,**</sup>, A. Mitchell<sup>18</sup>, R. Moderski<sup>28</sup>, M. O. Moghadam<sup>11</sup>, L. Mohrmann<sup>2</sup>, A. Montanari<sup>5</sup>, E. Moulin<sup>17</sup>, T. Murach<sup>10</sup>, M. de Naurois<sup>7</sup>, J. Niemiec<sup>20</sup>, S. Ohm<sup>10</sup>, L. Olivera-Nieto<sup>2</sup>, E. de Ona Wilhelmi<sup>10</sup>, S. Panny<sup>25</sup>, M. Panter<sup>2</sup>, R. D. Parsons<sup>14</sup>, U. Pensec<sup>16</sup>, S. Pita<sup>12</sup>, G. Pühlhofer<sup>31</sup>, M. Punch<sup>12</sup>, A. Quirrenbach<sup>5</sup>, M. Regeard<sup>12</sup>, A. Reimer<sup>25</sup>, O. Reimer<sup>25</sup>, H. Ren<sup>2</sup>, B. Reville<sup>2</sup>, F. Rieger<sup>2</sup>, B. Rudak<sup>28</sup>, E. Ruiz-Velasco<sup>2</sup>, V. Sahakian<sup>4</sup>, H. Salzmann<sup>31</sup>, A. Santangelo<sup>31</sup>, M. Sasaki<sup>18</sup>, F. Schüssler<sup>17</sup>, H. M. Schutte<sup>9</sup>, J. N. S. Shapopi<sup>8</sup>, H. Sol<sup>15</sup>, S. Spencer<sup>18</sup>, Ł. Stawarz<sup>26</sup>, R. Steenkamp<sup>8</sup>, S. Steinmassl<sup>2</sup>, C. Steppa<sup>11</sup>, K. Streil<sup>18</sup>, T. Takahashi<sup>32</sup>, T. Tanaka<sup>33</sup>, A. M. Taylor<sup>10</sup>, R. Terrier<sup>12</sup>, M. Tsirou<sup>10</sup>, C. van Eldik<sup>18</sup>, C. Venter<sup>9</sup>, J. Vink<sup>34</sup>, T. Wach<sup>18</sup>, S. J. Wagner<sup>5</sup>, A. Wierzcholska<sup>20</sup>, M. Zacharias<sup>5,9</sup>, A. A. Zdziarski<sup>28</sup>, A. Zech<sup>15</sup>, P. Zilberman<sup>10,35,**</sup>, and N. Żywucka<sup>9</sup>

*(Affiliations can be found after the references)*

Received 11 December 2023 / Accepted 15 February 2024

## ABSTRACT

The radio galaxy M 87 is a variable very-high energy (VHE) gamma-ray source, exhibiting three major flares, which were reported in 2005, 2008, and 2010. Despite extensive studies, the origin of the VHE gamma-ray emission is not yet fully understood. In this study, we investigate the VHE gamma-ray spectrum of M 87 during states of high gamma-ray activity, utilizing 20.2 h of H.E.S.S. observations. Our findings indicate a preference for a curved spectrum, characterized by a log-parabola model with extra-galactic background light (EBL) model above 0.3 TeV at the $4\sigma$ level, compared to a power-law spectrum with EBL. We investigate the degeneracy between the absorption feature and the EBL normalization and derive upper limits on EBL models that are mainly sensitive in the wavelength range of 12.4 µm–40 µm.

**Key words.** astroparticle physics – galaxies: active – galaxies: jets – gamma rays: galaxies – infrared: general

## 1. Introduction

Messier 87 (M 87) is a bright nearby radio galaxy at 16.8 Mpc from Earth (Akiyama et al. 2019), one of the rare radio galaxies from which TeV photons have been detected (Aharonian et al. 2003). This galaxy has been the subject of extensive study across the entire electromagnetic spectrum (see EHT Mwl Science Working Group 2021 for a review of recent observations). M 87 is a variable VHE gamma-ray source with three major TeV flares reported to date (see Aharonian et al. 2006a; Acciari et al. 2009; Abramowski et al. 2012). The potential origin of such rapid flares (~1 day) is linked to either the core of M 87 or to very compact regions in the jet, for instance, the HST-1 knot (Harris et al. 2006).

Understanding the TeV radiation from radio galaxies can help shed light on the radiation mechanisms in blazars. Various mechanisms have been proposed to explain the production of very-high-energy (VHE, $E > 100$ GeV) gamma-ray flares from M 87. These include leptonic models (for a review, see Rieger & Levinson 2018), magnetospheric models (MAGIC Collaboration 2018), multi-blob synchrotron self-Compton (SSC; Lenain et al. 2008), a photo-hadronic model (Sahu et al. 2019; Alfaro et al. 2022), or a cloud-jet interaction model (Barkov et al. 2012). Stochastic acceleration (Massaro et al. 2004; Tramacere et al. 2011) may lead, for instance, to the production of curved electron spectra, which ultimately leads to a curved gamma-ray spectrum. Close to the maximum particle energy in the accelerator, the resultant gamma-ray spectrum might follow a stretched cut-off shape, mimicking a curvature (Romoli et al. 2017). Furthermore, Klein-Nishina effects might also induce a cut-off in the gamma-ray spectrum at the highest energies (Lefa et al. 2012). Such features in the

<footnote>* Figures data points, FITS, and code example are available at https://zenodo.org/records/11066183</footnote>

<footnote>** Corresponding authors; e-mail: contact.hess@hess-experiment.eu</footnote>
spectrum of the detected gamma rays could point to the mechanism and location of acceleration of the particles, as well as the type of relativistic particle.

Regardless of the origin of the VHE flares from M87, the emitted gamma rays propagate to Earth, potentially interacting with photon fields en route, producing electron-positron pairs and leaving an imprint in the spectrum. We identified three photon fields (and regions) that could significantly alter the VHE gamma-ray spectrum of M87. Firstly, internal absorption, namely: the absorption of gamma rays at the core of M87 by the photons emitted in the accretion region (Brodatzki et al. 2011). Secondly, there is the gamma-ray absorption by the host galaxy's starlight (Stawarz et al. 2006; Zacharias et al. 2017) and, thirdly, the extra-galactic background light (EBL; e.g., Franceschini et al. 2019), all of which may leave detectable imprints in the observed spectrum of M87 (e.g., Fig. 11.5 from Aharonian 2004). These imprints help us identify the photon fields in the neighborhood of the acceleration site, thereby shedding light on the origin of the emission and on the mechanism of particle acceleration.

One approach to investigating spectral imprints is to examine the γ-γ absorption of TeV photons by photon fields from the region. Brodatzki et al. (2011) suggested that the theoretical optical depth $\tau(E, z)$ from internal γ-γ absorption reaches $\tau \gtrsim 1$ at $\gtrsim 20$ TeV, assuming that the VHE gamma rays are produced within $25r_s$ from the central black hole, where $r_s$ is the Schwarzschild radius of the supermassive black hole (SMBH). This estimate assumes a Shakura-Sunyaev disk and a mass of the SMBH of $(6.4 \pm 0.5) \times 10^9 M_\odot$, where $M_\odot$ is the solar mass. This result is in agreement with the latest estimates of the SMBH mass (Akiyama et al. 2019).

As we move further away from the SMBH, starlight becomes the dominant photon field. The galaxy extends up to a half-light radius of $R_{1/2} = 7.2$ kpc (Weil et al. 1997), namely, the radius within which half of the total starlight emitted by the galaxy is contained. By considering the presence of starlight from the host galaxy, we can estimate the effect of photon-photon absorption on the gamma-ray spectrum from the central AGN. (Zacharias et al. 2017). However, the lack of information about the spatial distribution of photons along the line of sight may complicate the interpretation of the results.

Beyond the M87 optical galaxy, the EBL becomes the dominant photon field in the µm wavelength range. The EBL encompasses all light emitted from stars and dust in galaxies throughout the Universe's history. It contains valuable information about stellar and galactic evolution. Accurate estimation of the EBL is crucial for studying extragalactic gamma-ray sources. Conversely, if an extragalactic source's intrinsic spectrum is well known, it can also provide a valuable avenue for estimating the EBL.

Attempts to estimate the EBL via attenuation of VHE photons primarily probe the 1 µm–50 µm range (see e.g., Aharonian et al. 2007; H.E.S.S. Collaboration 2013, 2017; Ahnen et al. 2016; Acciari et al. 2019; Abeysekara et al. 2019 or Biteau & Meyer 2022 for a review). This range overlaps with the range probed by direct measurements. Since direct measurements suffer from foreground contamination (Hauser et al. 1998), VHE gamma-ray absorption by the EBL provides complementary, indirect estimates.

The VHE flux from gamma-ray sources quickly falls off at higher energies, typically decreasing as a power law (PL) with a photon index of $\Gamma \sim 2-3$. Given that the optical depth $\tau$ increases monotonically with the gamma-ray energy, efforts to investigate its influence on the gamma-ray spectrum at the high-

est energies are restricted to bright and nearby VHE gamma-ray sources.

Due to its close proximity to Earth, M87 is a unique object that can be used to probe the 10 µm to 50 µm region of the EBL spectrum. According to several models, the theoretical optical depth should become measurable in the spectrum of M87 at $\gtrsim 10$ TeV (e.g., Fig. 11.5 from Aharonian 2004; Franceschini et al. 2019). There is an observable hardening of the spectrum of the source during high states ($\Gamma \approx 2.2$, Abramowski et al. 2012; H.E.S.S. Collaboration 2023), which presents us with an opportunity to probe the EBL even at the source's relative proximity.

The High Energy Stereoscopic System (H.E.S.S.) has good coverage of all the aforementioned TeV gamma-ray flares, which makes it the ideal instrument for studying the VHE gamma-ray spectrum of M87 up to the highest energies. Therefore, we selected high-state data sets from H.E.S.S. observations of M87 and search for absorption imprints in it. The VHE gamma-ray observations of M87 can also be used to search for physics beyond the standard model. Using the same data set, Cecil et al. (2023) searched for spectral signatures of oscillations between photons and axion-like particles in the magnetic field of the Virgo Cluster. They found no significant evidence for such oscillations.

This publication is organized as follows. In Sect. 2, we introduce the H.E.S.S. experiment, the data set, and the methods used in the data analysis. In Sect. 3, we present the results of the spectral fits and in Sect. 4, we discuss the implications of our findings. We give a summary in Sect. 5.

## 2. Observations and data analysis

H.E.S.S. is an array of five Imaging Atmospheric Cherenkov Telescopes (IACTs) situated in the Khomas Highland of Namibia. The H.E.S.S. telescopes have been regularly used, initially with four telescopes, to observe M87 since 2004, both in monitoring campaigns and in response to triggers from external instruments reporting flaring activities. Observations are conducted in individual runs, each typically lasting up to 28 min. This study utilizes H.E.S.S. observations of M87 from 2004 to 2020. Although some of the data have already been used in previous publications (Aharonian et al. 2006a; Acciari et al. 2009; Abramowski et al. 2012; EHT Mwl Science Working Group 2021; H.E.S.S. Collaboration 2023), here we focus, for the first time, specifically on the high state by combining data from individual observational runs with elevated flux. To have a consistent data set, we only considered data from the four telescopes with 12 m diameter mirrors (CT1–4), which have been in operation since the earliest detected flare. Given that the 28 m diameter mirror telescope (CT5) has only been in operation since 2012, it has not observed the aforementioned flaring states of M87. Hence, we chose to exclude its data from this study, ensuring a consistent representation of M87 in VHE gamma rays over time. We applied the quality criteria for a spectral analysis (Aharonian et al. 2006b), selecting only the best quality observations with at least three telescopes participating in the observation.

The data were processed using *gammapy v.1.1* (Donath et al. 2023), based on the generation of FITS files (Nigro et al. 2019) according to a multivariate gamma-selection and geometrical event reconstruction (Aharonian et al. 2006b; Ohm et al. 2009). We selected gamma-ray-like events within a cone of radius $\le 0.11^\circ$ around the position of the core of M87 at right ascension (RA) of 187.7059° and declination (Dec) 12.3911° for J2000
<figure><img src="image_1.png" /><figcaption>Fig. 1. Gamma-ray flux above 1 TeV per observation run. The high state is defined by the 10% highest flux above 1 TeV, as indicated by the red data points (above the 90% percentile of the distribution). In the inset on the right side, a histogram with number of observational runs below and above the 90% percentile is shown.</figcaption></figure>

(Lambert & Gontier 2009). The background was estimated utilizing the reflected-region background method (Berge et al. 2007), with a minimum angular separation to M 87 of 0.15°.

In Fig. 1, we present the run-wise flux of the selected data above 1 TeV, emphasizing observational runs with elevated flux. All data points exhibiting a flux exceeding $10^{-13} \text{ cm}^{-2} \text{ s}^{-1}$ are included. For clarity, uncertainties associated with the flux points are omitted in order to maintain a clean visualization. It is important to note that each flux point does not necessarily denote a detection but rather represents the estimated flux for each observation. Based on an empirical cumulative distribution function, we selected 10% of the runs with the highest flux, namely, the runs with the flux above the 90% percentile of the distribution to compose the high state. The choice of the 90% percentile is based on an expected effective duty cycle of the source, which is also influenced by the H.E.S.S. observation strategies in the past. The high state runs are marked in red in Fig. 1. This amounts to 20.2 h of collected data. The definition of high state ensures that we selected only the observations from the tail of the high-state flux distribution (see Fig. 1 inset on the right panel).

We regard the selection of the high state based on the 10% highest flux above 1 TeV as stringent, evidenced by the exclusion of the 2008 flare from our selection, which H.E.S.S. did not observe at its peak intensity. Some observation runs conducted during the 2010 flare were excluded from the dataset. This exclusion was due to these observations being performed on the target without the usual wobbling motion. The typical wobble is essential for obtaining reliable background and spectral estimates (Berge et al. 2007). In 2018, M 87 underwent a VHE gamma-ray flare, which is to be reported in depth in a forthcoming publication within the Event Horizon Telescope Multi-wavelength Group. Hence, our high state dataset predominantly comprises data from the 2005 and the 2018 flares. As a systematic check, we study the influence of the high state selection into the final spectral results and the rise of the curvature in the high state data set as shown in Appendix B.

We used the likelihood $\mathcal{L} = \prod_i P(n_i|\mu_i)$, and $P(n_i|\mu_i)$ as the Poisson probability of observing $n_i$ events in energy bin, $i$, given the predicted number of events, $\mu_i$, from the background and the source models, to fit the data to spectral models. We explore several spectral models to describe the differential gamma-ray flux, $\phi(E)$, as a function of energy, $E$, each with increasing complexity. The first model is a power law (PL) spectrum, defined as:

$$ \phi(E) = \phi_0 \left(\frac{E}{E_0}\right)^{-\Gamma}, \qquad (1) $$

where the normalization, $\phi_0$, and the spectral index, $\Gamma$, are free parameters, while the reference energy, $E_0$, is fixed to the decorrelation energy, that is, to the energy at which the correlation between the spectral model parameters is minimized.

To account for potential curvature in the spectrum, we considered a log-parabola (LP) spectral model:

$$ \phi(E) = \phi_0 \left(\frac{E}{E_0}\right)^{-\Gamma - \beta \ln(E/E_0)}, \qquad (2) $$

where the coefficient is $\beta$, the normalization is $\phi_0$, and the spectral index, $\Gamma$, represents the free parameters.

To incorporate the EBL absorption into our spectral models, we initially reviewed a range of models from the literature, eventually selecting three that comprehensively represent the range of options. We started with the kneiske model (Kneiske & Dole 2010), tailored to emulate the lower boundary of the EBL flux. For a contemporary estimate of the EBL founded on data, we turned to the finke2022 model, specifically referencing model A from Finke et al. (2022a). To round off our selection with a model reflecting high-intensity EBL, we incorporated the upper bounds of the Dominguez model, termed dominguez-upper, which mirrors the upper uncertainties outlined in Dominguez et al. (2011). Collectively, these three models capture the diverse models from literature regarding the optical depth, $\tau$.

The models characterizing the EBL spectrum, presented as the optical depth, $\tau$, in relation to the gamma-ray energy, are commonly expressed as functions of redshift. We used the tables from Meyer (2022) considering the redshift of M 87 as $z \approx 0.0042$. However, considering the potential influence of the local Hubble flow on the redshift measurements of nearby sources, we also took into account a direct distance measurement to M87 of 16.8 Mpc (Bird et al. 2010; Akiyama et al. 2019). This distance was converted to a redshift $z = 0.00389$ using the Planck18 cosmology (Planck Collaboration VI 2020) with the assistance of Astropy v.5.3.3 (The Astropy Collaboration 2022). The resulting variations between the two redshift values in terms of EBL optical depth are below 10%. We base our estimates on the originally estimated redshift.

Finally, we incorporate an EBL component into the PL and LP spectral models and define the notation of the new models as PLxEBL and LPxEBL spectral models, respectively:

$$ \phi(E) = \phi_0 \left(\frac{E}{E_0}\right)^{-\Gamma} \exp(-\alpha_{\text{norm}} \tau(E, z)), \qquad (3) $$

$$ \phi(E) = \phi_0 \left(\frac{E}{E_0}\right)^{-\Gamma - \beta \ln(E/E_0)} \exp(-\alpha_{\text{norm}} \tau(E, z)), \qquad (4) $$

where $\tau(E, z)$ is the model-dependent EBL optical depth, parameterized by the photon energy and redshift, and $\alpha_{\text{norm}}$ is the EBL
<table><caption>**Table 1.** Best-fit results for the spectral VHE gamma-ray distribution of the high state of M 87.</caption><thead><tr><th>Model</th><th>-2lnL</th><th>φ<sub>0</sub></th><th>Γ</th><th>β</th><th>Ref. energy</th><th>Energy range</th></tr></thead><tbody><tr><td>PLxEBL finke2022</td><td>30.68</td><td>4.7 ± 0.3</td><td>1.81 ± 0.08</td><td>–</td><td>1.96</td><td>0.3–10</td></tr><tr><td>LPxEBL finke2022</td><td>18.69</td><td>6.0 ± 0.5</td><td>1.76 ± 0.10</td><td>0.38 ± 0.13</td><td>1.96</td><td>0.3–10</td></tr><tr><td>PLxEBL kneiske</td><td>36.23</td><td>3.7 ± 0.2</td><td>1.94 ± 0.06</td><td>–</td><td>2.15</td><td>0.3–32</td></tr><tr><td>LPxEBL kneiske</td><td>16.86</td><td>4.8 ± 0.4</td><td>1.81 ± 0.08</td><td>0.29 ± 0.07</td><td>2.15</td><td>0.3–32</td></tr><tr><td>PLxEBL finke2022</td><td>34.38</td><td>3.8 ± 0.2</td><td>1.92 ± 0.06</td><td>–</td><td>2.15</td><td>0.3–32</td></tr><tr><td>LPxEBL finke2022</td><td>17.04</td><td>4.9 ± 0.4</td><td>1.80 ± 0.08</td><td>0.27 ± 0.08</td><td>2.15</td><td>0.3–32</td></tr><tr><td>PLxEBL dominguez-upper</td><td>31.28</td><td>3.9 ± 0.2</td><td>1.88 ± 0.06</td><td>–</td><td>2.15</td><td>0.3–32</td></tr><tr><td>LPxEBL dominguez-upper</td><td>18.20</td><td>4.9 ± 0.4</td><td>1.79 ± 0.08</td><td>0.25 ± 0.08</td><td>2.15</td><td>0.3–32</td></tr></tbody></table>

**Notes.** We compared the PL and the respective LP best-fit spectral model for the reduced energy range and for the full energy range for each EBL model considered. The energy range is given in TeV and $\phi_0$ is given in units of $10^{-13} \text{ cm}^{-2} \text{ s}^{-1} \text{ TeV}^{-1}$. The decorrelation energy of the PL spectral model was used as the reference energy, given in TeV in the table.

normalization. While $\tau(E, z)$ is defined by the EBL model, $\alpha_{norm}$ can either be set to vary freely, in case one is interested in constraining the EBL, or be fixed to the expected value $\alpha_{norm} = 1$ to probe the intrinsic spectral model. As we explain in Sect. 3, we utilized both approaches.

We used the implementation from gammapy to fit the spectral models to the data following the forward-folded approach (Piron et al. 2001) and obtain the best-fit parameters with their loglikelihood values ($-2 \ln L$). For nested models, the preference of a model 1 in comparison to model 0, where model 1 has more free parameters, is given by the test statistic $TS = -2 \ln(L_0/L_1)$. For large event statistics, $TS$ follows a $\chi^2$ distribution with $k$ degrees of freedom, where $k$ is the difference of free parameters between the two nested models. This representation allows us to express the results in terms of the significance of the fit. Since $k = 1$ for the PL and the LP spectral model comparison, the significance of the LP in comparison to the PL spectral model is approximated by $\sqrt{TS}$.

In addition to the main analysis described in this section, we performed two cross-checks to assess the reliability of the results. The first cross-check utilized an alternative high-state definition, based on H.E.S.S. Collaboration (2023). This approach led to a slightly different selection of observations, although the resulting best-fit spectral models are within $1\sigma$ statistical uncertainties consistent with each other. The second cross-check utilized the same observations as the first cross-check with an independent analysis chain (de Naurois & Rolland 2009). The results were also found to be consistent with the main analysis.

### 3. H.E.S.S. analysis results

For the high-state dataset defined in Sect. 2, we obtained $390 \pm 28$ excess gamma-ray events and $376 \pm 5$ expected background events for a livetime of 20.2 h. The significance of detection is calculated based on Eq. (17) of Li & Ma (1983) and yields $16.9\sigma$. For energies above 10 TeV, we find $15 \pm 4$ gamma-ray excess events for $2.3 \pm 0.4$ expected background events, which leads to a significant detection of $6.0\sigma$.

Given the distance to M 87, we find that the optical depth $\tau < 0.2$ holds for energies <10 TeV for all three considered EBL models. We considered the effects of EBL absorption for $\tau < 0.2$ as small and focused first on the analysis of the gamma rays with energy between 0.3 TeV and 10 TeV. Based on the methodology described in Sect. 2, we fit the PLxEBL finke2022 and the LPxEBL spectral models to the data in the reduced energy range.

<figure><img src="image_2.png" /><figcaption>**Fig. 2.** Detected spectral energy distribution of M 87 data with H.E.S.S. The best-fit LPxEBL finke2022 spectral models are shown in yellow for the reduced energy range (0.3–10 TeV) and in purple for the full energy range (0.3–32 TeV). The best-fit PL spectrum for the low-state data (based on H.E.S.S. Collaboration 2023) is shown for comparison in blue. The color bands show the $1\sigma$ uncertainty contours. The upper limit (UL) in the last bin of the low-state spectrum is defined at a 95% confidence level (c.l.).</figcaption></figure>

The results, shown in the first two rows of Table 1 and in yellow in Fig. 2, indicate a preference for an LPxEBL model over the PLxEBL model with a significance of $3.5\sigma$. We further discuss the implications of this curvature in the following section.

In a second step, we considered the full energy range (0.3 TeV–32 TeV). We fit PLxEBL and LPxEBL models now for the three EBL models kneiske, finke2022, and dominguez-upper, with fixed $\alpha_{norm} = 1$ to the data. The results are given in Table 1 for all the models and in Fig. 2 by the red curve and data points. The curvature in the spectrum for the full energy range (0.3–32 TeV) is detected, given that the LPxEBL models are preferred over the PLxEBL models with $4.4\sigma$, $4.2\sigma$, and $3.6\sigma$ for the kneiske, finke2022, and dominguez-upper EBL models, respectively. For the sake of readability, in the following discussions, we primarily refer to the finke2022 LPxEBL model, as it stands as the most recent and representative EBL model among the three considered. The energy-integrated flux above 0.3 TeV amounts to $(5.1 \pm 0.5) \times 10^{-12} \text{ cm}^{-2} \text{ s}^{-1}$ for the finke2022 LPxEBL model.

Lastly, we evaluated the LPxEBL spectral model, as previously defined, in conjunction with the same three EBL models
<figure><img src="image_3.png" /><figcaption>**Fig. 3.** Parameter space around the best-fit LPxEBL spectral model for the $\alpha_{\text{norm}}$ and $\beta$ parameters of the finke2022 EBL model. The color map gives the distance in standard deviations from the best-fit parameters (shown by the black cross), i.e., the lowest value provides the best fit. The black contours give the $1\sigma$, $2\sigma$, $3\sigma$, and $4\sigma$ uncertainty contours. The vertical green lines show the $\alpha_{\text{norm}} = 1$ and $\alpha_{\text{norm}} = 5.5$, which provides the 95% c.l. UL on the $\alpha$ parameter. For comparison, the dotted red line shows the $\alpha_{\text{norm}}$ that corresponds to the UL above which the finke2022 EBL model would start to contradict the ULs from Biteau & Williams (2015).</figcaption></figure>

and with the EBL normalization $\alpha_{\text{norm}}$ allowed to vary freely. Upon fitting the LPxEBL with a variable $\alpha_{\text{norm}}$ to the full energy range dataset, the fit converged to a pure LP model. This suggests a degeneracy between the curvature parameter $\beta$ and the EBL normalization $\alpha_{\text{norm}}$.

To elucidate the observed degeneracy, we present the log-likelihood values for the parameter space of $\alpha_{\text{norm}}$ and $\beta$ within the LPxEBL model. The result for the finke2022 EBL model is shown in terms of $\sqrt{\text{TS}}$ from the best fit position in Fig. 3. Based on the intersection of the $\alpha_{\text{norm}}$ with the $2\sigma$ contour in Fig. 3, we derive the maximum allowed EBL normalization $\alpha_{\text{norm}}$ within 95% c.l. of $\alpha_{\text{norm}} < 5.5$, for the finke2022 EBL model. Similarly, for the kneiske, and dominguez-upper EBL models, we derive $\alpha_{\text{norm}} < 8.7$ and $\alpha_{\text{norm}} < 2.0$, respectively. In the case of $\beta = 0$, the best-fit position in the parameter space is always more than $4\sigma$ away from the fit for an EBL intensity for $\alpha_{\text{norm}} \le 1.5$. Since the gamma-ray energy range used is between 0.3 TeV and 32 TeV, our limits on the EBL intensity probe the $\approx 0.4$ µm–40 µm EBL photons<fnref n="1" />, although our analysis is more sensitive for gamma rays $\ge 10$ TeV, i.e., for EBL photons $\ge 12.4$ µm, where their optical depth reaches $\tau \gtrsim 0.2$.

The curvature observed in the spectrum of the high state of M 87 was analyzed in this study using a LP distribution. Nevertheless, it is worth noting that an exponential cut-off model, despite having an equivalent number of degrees of freedom, may suggest different underlying physics. Hence, in Appendix C, we explore the PL with an exponential cut-off, defined as the PLxECEBL model, taking into consideration the effects of the EBL.

## 4. Discussion and conclusions

A curved gamma-ray spectrum has been already observed up to TeV energies in the radio galaxy NGC1275 (MAGIC Collaboration 2018), as well as in other gamma-ray sources (Zdziarski et al. 2020). Various mechanisms could underlie this spectral feature, suggesting its potential common occurrence in nature. In this section, we analyze possible physics scenarios that could lead to a curvature in the VHE gamma-ray spectrum of the high state of M 87. We considered four possibilities and discuss them in the subsequent paragraphs based on the details of our analysis: 1) external absorption: EBL photons; 2) multi-component high state: a composition of different particle populations; 3) internal absorption: by local photon fields, for instance, the accretion flow and the star light; 4) intrinsic curved electron spectrum: due to the nature of the mechanism of acceleration, for instance, stochastic acceleration, or due to a cut-off in the electron spectrum.

The presence of the curvature in the high-state spectrum of M 87 at VHE energies in the reduced energy range (0.3–10 TeV; first two rows in Table 1) indicates that the curvature is not related to the EBL absorption, which is negligible in this energy range. This hypothesis gained further support from the minimal discrepancies observed when fitting the spectral data within the reduced energy range, taking into account the EBL absorption within the model. Although showing a compatible curvature, the spectral curvature is more strongly detected in the full energy range (0.3–32 TeV; Table 1). Hypotheses with more complex spectral behavior cannot be tested due to the limited event statistics in the dataset (Appendix C).

For the full energy range (0.3–32 TeV), we find that the EBL normalization ($\alpha_{\text{norm}}$) and the intrinsic curvature coefficient from the LP model ($\beta$) are statistically degenerate. However, given that the EBL is in fact directly obtained from observations at optical and infrared wavelengths, (e.g., Domínguez et al. 2011; Saldana-Lopez et al. 2021), the EBL normalization $\alpha_{\text{norm}}$ is expected to be at least close to one and, therefore, it cannot be set to zero in expense of $\beta$ as may be suggested by Fig. 3. It is clear through the $4\sigma$ statistical uncertainty contours from Fig. 3 that the curvature in the spectrum is strongly detected at a $4\sigma$ level even by considering $\alpha_{\text{norm}} = 1$. Therefore, based on current available EBL models, we rule out the external absorption by EBL photons as the exclusive cause for the observed curvature in this study. The existing EBL models utilize data collected from all directions in the sky and do not consider potential anisotropies in the EBL at the same redshift. This lack of consideration may influence the predicted absorption feature in gamma-ray spectra. Furniss et al. (2015), for instance, did not find a strong correlation between the location of voids and hard gamma-ray sources in their investigation. However, they did estimate a 2% variation in the density of EBL photons from voids compared to the nominal EBL density. Notably, EBL anisotropy remains a relatively unexplored subject (Hervet et al. 2023).

Since our data set is a collection of high emission states, the overall VHE gamma-ray emission from the high state of M 87 (as defined here) results from the contribution of the individual flares. Although the flaring states of M 87 are usually accompanied by a hard spectral index ($\Gamma \approx 2$), it is not trivial to conclude from it that the individual flares originate from the same region (Abramowski et al. 2012). Therefore, different emission regions with distinguished electron populations might be behind the signal. We tested this hypothesis by fitting the data around the 2005 flare and the 2018 flare, as they are the largest contributors to

<footnote n="1">¹ This wavelength range is obtained from the peak of the pair production cross section at the wavelength $\lambda_{\text{max.}} \approx 1.24(E_{\gamma}/\text{TeV})$ µm, where $E_{\gamma}$ is the gamma-ray energy (Franceschini et al. 2019).</footnote>
the selected high state data (Fig. 1), separately. The results show a hint of $1.5\sigma$ and $2.9\sigma$ preference for an LPxEBL finke2022 model in comparison to the PLxEBL finke2022 model, respectively, in the 2005 and 2018 flares. Furthermore, the spectral indices of the PL models from the 2005 and the 2018 flares are $2.01 \pm 0.09$ and $1.9 \pm 0.1$, respectively, which is consistent with each other. Therefore, the hint of curvature, especially from the 2018 dataset, and the similar spectral shapes suggest that the curvature detected in the combined dataset is not due exclusively to the superposition of different spectral components.

The presence of a persistent steady emission component even during flaring episodes could also contribute to the spectral shape, especially in its lower energy part ($E \approx 0.3$ TeV). In fact, the low state of M 87 is well-described by a PL distribution with a steep spectral index (H.E.S.S. Collaboration 2023, $\Gamma = 2.63 \pm 0.09$), where no significant spectral curvature was found in the H.E.S.S. data. If the low and high states of the VHE gamma-ray emission of M 87 are indeed due to different components, the high-state component should be curved and dominate over the steady PL component during flares to yield the observed SED. However, the low-state component would still contribute, at the same level as the high-state component to the flux of gamma rays with energies around $E \approx 0.3$ TeV during flaring activities (see Fig. 2). This hypothesis is also supported by the fact that Fermi-LAT observations did not detect a significant increase in the high-energy ($100 \text{ MeV} < E < 300 \text{ GeV}$) flux during the 2010 VHE gamma-ray flare (Abramowski et al. 2012). Furthermore, the Fermi-LAT flux points (regular state in Benkhali et al. 2019) smoothly connect to the H.E.S.S. flux points of the low state PL distribution at $\approx 0.3$ TeV. Therefore, even in high emission states, we expect an up-rise in the curvature (Fig. B.2) for energies below $0.3$ TeV, where the low-state would dominate the high-state component. The previous arguments support the hypothesis that the VHE gamma-ray emission of M 87 is composed separately of a low-energy and a high-energy component, whereas the VHE high-energy component appears during flaring activities. Nonetheless, a single emission region responsible for the emission in both source emission states cannot yet be discarded.

As for the internal absorption, Brodatzki et al. (2011) estimates that the inner region of M 87 should be transparent to gamma rays with energies $\le 20$ TeV. The absorption by the starlight from the host galaxy can also be estimated based on its emissivity. We employed the surface brightness of M 87 (UGC 7654) in the $1.4 \mu\text{m}$-$1.7 \mu\text{m}$ range, denoted as $\mu_\lambda = 15.51 \text{ mag arcsec}^{-2}$ from Capetti & Balmaverde (2006). Using Eq. (A5) of Zacharias et al. (2017), we convert the surface brightness to emissivity $j_\lambda$:

$$j_\lambda = 2.32 \times 10^{-5} \left(\frac{r_0}{\text{cm}}\right)^{-1} 10^{0.4 \times (25.992 - (\mu_\lambda - A))} \left(\frac{\text{erg}}{\text{cm}^3 \text{s}^{-1}}\right), \quad (5)$$

where we set the galactic extinction A to zero, as it has already been corrected for the $\mu_\lambda$ estimate and $r_0$ is the characteristic radius of the source. We used $r_0 \approx 432$ pc, half the size of the image in Capetti & Balmaverde (2006) used for the brightness estimate, converted to parsec by considering the distance to M 87 of 16.8 Mpc. For a radial distribution of light, the optical depth $\tau_{M87,\lambda}$ for the galactic emission up to $r_0$ centered at $\lambda$ is estimated through:

$$\tau_{M87,\lambda} = \frac{\sigma_{\gamma\gamma} j_{\lambda} \lambda r_0^2}{hc^2}, \quad (6)$$

where $\sigma_{\gamma\gamma}$ is the cross section for $\gamma-\gamma$ absorption, c is the speed of light, and h is the Planck constant. We utilized the full cross

section (Gould & Schréder 1967) for an average cosine angle of interaction between the incoming gamma rays and the photons from the star light of 0.5. Our estimates indicate that the galactic optical depth at $\lambda = 1.6 \mu\text{m}$ peaks for a gamma-ray energy of $\approx 1$ TeV, although $\tau_{M87,\lambda} < 0.001$ in the entire gamma-ray energy range, that is, from the threshold energy of pair production up to 32 TeV. Therefore, gamma-ray absorption by starlight is weak and unlikely to explain the reported curvature here.

The curvature in the gamma-ray spectrum can also be associated with an intrinsically curved particle spectrum, typically explained by stochastic particle accelerations (Massaro et al. 2004; Tramacere et al. 2011). The event statistics of our dataset do not allow us to distinguish between a cutoff on the electron spectrum due to the maximum electron energy in the source or a cutoff due to the onset of the Klein-Nishina regime (Aghanim et al. 2020) from an intrinsic parabolic curvature. The rapid VHE gamma-ray flares undergone by M 87 with hard spectral indices point to an efficient acceleration mechanism. However, we are not able to distinguish the acceleration mechanism based solely on the spectral shape of the VHE gamma rays. In addition to that, the exact acceleration site also remains an open question (see Rieger & Aharonian 2012).

## 5. Summary

In this work, we study the spectrum of the high-state emission of M 87 with H.E.S.S. data. We define the high state as the 10% of the observation runs with the highest flux above 1 TeV (Fig. 1). This accounts for a total of 20.2 h of observations.

We first consider the energy range between 0.3 and 10 TeV, in which the EBL absorption is small, namely, $\tau < 0.2$. We then fit a power-law (PL) distribution with the finke2022 (Finke et al. 2022a) EBL model, as the null hypothesis, to describe the spectrum and test alternative spectral models with curvature. The result shows a strong ($4.2\sigma$) preference for a log parabola (LP) spectral model with the EBL (LPxEBL) in comparison to the PL model with EBL (PLxEBL), which establishes a curved spectrum in the VHE gamma-ray emission from the high state of M 87.

In a second analysis step, we consider the full energy range (0.3–32 TeV) and fit an LPxEBL to the data. We utilized three EBL models at $z \approx 0.0042$ that cover a multitude of other models. Therefore, our selection is representative (Meyer 2022): the kneiske model (Kneiske & Dole 2010), the finke2022 model (model A in Finke et al. 2022b), and the dominguez-upper model (the upper uncertainties in Domínguez et al. 2011). We confirm that the curvature in the spectrum is also strongly detected considering the full energy range as the LPxEBL models are preferred at $\approx 4\sigma$ over the PLxEBL models. The total flux resulted from the analysis above 0.3 TeV is estimated to be $(5.1 \pm 0.5) \times 10^{-12} \text{ cm}^{-2} \text{s}^{-1}$ for the finke2022 EBL model. The curvature is characterized by a spectral index $\Gamma = 1.80 \pm 0.08$ and a curvature coefficient of $\beta = 0.27 \pm 0.08$ for the LPxEBL finke2022 model, which is also compatible with the results of the LPxEBL with dominguez-upper and kneiske EBL models. Although the spectral curvature is somewhat more significant in analysis of the full energy range in comparison to the reduced energy range (Table 1), they are compatible with each other.

In this study, we also investigate the degeneracy between the EBL absorption feature and an intrinsic curved spectrum (Fig. 3). We conclude that, despite not being able to separate both dependencies completely, the spectral curvature is detected even when considering the EBL with normalization as predicted by the models ($\alpha_{\text{norm}}=1$). Furthermore, we derive the maximum
allowed EBL normalization, $\alpha_{\text{norm}}$, within 95% confidence level of $\alpha_{\text{norm}} < 8.7$, $\alpha_{\text{norm}} < 5.5$, $\alpha_{\text{norm}} < 2.0$, for the kneiske, finke2022, and dominguez-upper EBL models, respectively. In addition to the LPxEBL model, we also fitted the high-state data with a spectrum described by an PL with exponential cut-off (Appendix C), which leads to a slightly worse description of the data.

In summary, we measure for the first time a curvature in the VHE gamma-ray spectrum of the high state of M 87. We rule out external absorption by EBL photons as its sole sources and we deem multiple particle populations in the high state as an unlikely explanation for the detected curvature. Furthermore, our estimates show that absorption by star light at 1.6 µm is minimal and unlikely to explain the reported curvature. Finally, intrinsic curved radiation spectrum, namely, due to a cut off in the particle spectrum at the particle accelerator or due to the onset of the Klein-Nishina effect, also stand as plausible explanations for the observed curvature.

**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 U.K. 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.

## References

Abeysekara, A. U., Archer, A., Benbow, W., et al. 2019, ApJ, 885, 150

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

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

Acciari, V. A., Ansoldi, S., Antonelli, L. A., et al. 2019, MNRAS, 486, 4233

Aghanim, N., Akrami, Y., Ashdown, M., et al. 2020, A&A, 641, A6

Aharonian, F. A. 2004, Very High Energy Cosmic Gamma Radiation - A Crucial Window on the Extreme Universe (World Scientific Publishing Co., Pte. Ltd.)

Aharonian, F., Akhperjanian, A., Beilicke, M., et al. 2003, A&A, 403, L1

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., Barres de Almeida, U., et al. 2007, A&A, 475, 9

Ahnen, M. L., Ansoldi, S., Antonelli, L. A., et al. 2016, A&A, 590, A24

Akiyama, K., Alberdi, A., Alef, W., et al. 2019, ApJ, 875, 44

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

Barkov, M. V., Bosch-Ramon, V., & Aharonian, F. A. 2012, ApJ, 755, 170

Benkhali, F. A., Chakraborty, N., & Rieger, F. M. 2019, A&A, 623, A2

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

Biteau, J., & Meyer, M. 2022, Galaxies, 10, 1

Biteau, J., & Williams, D. A. 2015, ApJ, 812, 26

Brodatzki, K. A., Pardy, D. J. S., Becker, J. K., & Schlickeiser, R. 2011, ApJ, 736, 98

Capetti, A., & Balmaverde, B. 2006, A&A, 453, 27

Cecil, R., Martins Barbos, V., Lypova, I., et al. 2023, PoS(ICRC2023), 444, 908 de Neurois, M., & Rolland, L. 2009, Astropart. Phys., 32, 231

Domínguez, A., Primack, J. R., Rosario, D. J., et al. 2011, MNRAS, 410, 2556

Donath, A., Sinha, Q., Remy, F., et al. 2023, A&A, 678, A157

EHT Mwl Science Working Group 2021, ApJ, 911, L11

Finke, J. D., Ajello, M., Domínguez, A., et al. 2022a, ApJ, 941, 33

Finke, J. D., Ajello, M., Domínguez, A., et al. 2022b, https://zenodo.org/records/7023073

Franceschini, L., Foffano, E., Prandini, F., et al. 2019, A&A, 629, A2

Furniss, A., Sutter, P. M., Primack, J. R., & Domínguez, A. 2015, MNRAS, 446, 2277

Gould, R. J., & Schréder, G. P. 1967, Phys. Rev., 155, 1404

Harris, D. E., Cheung, C. C., Biretta, J. A., et al. 2006, ApJ, 640, 211

Hauser, M., Arendt, R., Kelsall, T., et al. 1998, ApJ, 508, 25

Hervet, O., Williams, D., & Furniss, A. 2023, PoS, ICRC2023, 753

H.E.S.S. Collaboration (Abramowski, A., et al.) 2013, A&A, 550, A4

H.E.S.S. Collaboration (Abdalla, H., et al.) 2017, A&A, 606, A59

H.E.S.S. Collaboration 2023, A&A, 675, A138

Kneiske, T. M., & Dole, H. 2010, A&A, 515, A19

Lambert, S. B., & Gontier, A. M. 2009, A&A, 493, 317

Lefa, E., Kelner, S. R., & Aharonian, F. A. 2012, ApJ, 753, 176

Lenain, J.-P., Boisson, C., Sol, H., & Katarzyński, K. 2008, A&A, 478, 111

Li, T.-P., & Ma, Y.-Q. 1983, ApJ, 272, 317

MAGIC Collaboration (Ansoldi, S., et al.) 2018, A&A, 617, A91

Massaro, E., Perri, M., Giommi, P., et al. 2004, A&A, 413, 489

Meyer, M. 2022, https://doi.org/10.5281/zenodo.7312062

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

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

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

Planck Collaboration VI. 2020, A&A, 641, A6

Rieger, F. M., & Aharonian, F. 2012, Mod. Phys. Lett. A, 27, 1

Rieger, F. M., & Levinson, A. 2018, Galaxies, 6, 1

Romoli, C., Taylor, A. M., & Aharonian, F. 2017, Astropart. Phys., 88, 38

Sahu, S., Fortín, C. E. L., & Nagataki, S. 2019, ApJ, 884, L17

Saldana-Lopez, A., Domínguez, A., Pérez-González, P. G., et al. 2021, MNRAS, 507, 5144

Stawarz, L., Aharonian, F., Wagner, S., & Ostrowski, M. 2006, MNRAS, 371, 1705

The Astropy Collaboration 2022, ApJ, 935, 167

Tramacere, A., Massaro, E., & Taylor, A. M. 2011, ApJ, 739, 2

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

Zacharias, M., Chen, X., & Wagner, S. J. 2017, MNRAS, 465, 3767

Zdziarski, A. A., Malyshev, D., de Oña Wilhelmi, E., et al. 2020, MNRAS, 455, 1451

<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 Yerevan Physics Institute, 2 Alikhanian Brothers St., 0036 Yerevan, Armenia</footnote>

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

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

<footnote n="7">7 Laboratoire Leprince-Ringuet, École Polytechnique, CNRS, Institut Polytechnique de Paris, 91128 Palaiseau, France</footnote>

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

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

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

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

<footnote n="12">12 Université de Paris, CNRS, Astroparticule et Cosmologie, 75013 Paris, France</footnote>

<footnote n="13">13 Department of Physics and Electrical Engineering, Linnaeus University, 351 95 Växjö, Sweden</footnote>

<footnote n="14">14 Institut für Physik, Humboldt-Universität zu Berlin, Newtonstr. 15, 12489 Berlin, Germany</footnote>

<footnote n="15">15 Laboratoire Univers et Théories, Observatoire de Paris, Université PSL, CNRS, Université Paris Cité, 5 Pl. Jules Janssen, 92190 Meudon, France</footnote>
16 Sorbonne Université, Université Paris Diderot, Sorbonne Paris Cité, CNRS/IN2P3, Laboratoire de Physique Nucléaire et de Hautes Énergies, LPNHE, 4 Place Jussieu, 75252 Paris, France

17 IRFU, CEA, Université Paris-Saclay, 91191 Gif-sur-Yvette, France

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

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

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 School of Physical Sciences, University of Adelaide, Adelaide 5005, Australia

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

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

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

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

29 Aix-Marseille Université, CNRS/IN2P3, CPPM, Marseille, France

30 Université Bordeaux, CNRS, LP2I Bordeaux, UMR 5797, 33170 Gradignan, France

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

32 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

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

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

35 Department of Physics, University of Wisconsin-Madison, Madison, WI 53706, USA
# Appendix A: EBL models and derived ULs

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

**Fig. A.1.** Optical depth (top) and the relative intensity of the gamma-ray flux after absorption (bottom) at $z \approx 0.0042$ of the three EBL models employed in this study are presented based on the information provided in the table from Meyer (2022). The two dashed lines in the upper panel indicate the energy and attenuation levels at which we deem the influence of the EBL spectrum on the overall spectrum as non-negligible. Additionally, for comparison, the remaining available EBL models are depicted in faded black. In the bottom panel, a reference dashed line is provided to guide the eyes in the case of no absorption.

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

**Fig. A.2.** EBL photon distribution based on the three models employed in this study - kneiske, finke2022, and dominguez-upper models - is depicted in orange, blue, and green, respectively. The 95% confidence level upper limits (ULs) derived from fitting the high state, as defined in the main text, to LPxEBL models are represented by dashed lines.

In this section, we delve into the optical depths of the three EBL models employed in this study and present the derived upper limits (ULs). The chosen models were selected to encompass a wide range of potential attenuation scenarios and are representative, as illustrated in Fig. A.1, within the gamma-ray energy range considered in our investigation.

Based on the results obtained from fitting the LPxEBL models to the high-state data of M,87 and exploring the interplay between the curvature parameter and the EBL normalization, we have derived, in the main text, 95% confidence level upper limits (ULs) on $\alpha_{\text{norm}}$. Concluding this study, in Fig. A.2, we present the spectral distribution of the EBL photons along with the corresponding ULs. The ULs on the EBL distribution are show-

cased within the gamma-ray energy range of 10 TeV to 32 TeV (EBL photon energy range between 12.4 and 40 µm), where our work exhibits the highest sensitivity. It is noteworthy that our strongest UL pertains to the brightest EBL model (dominguez-upper). Despite not very constraining, ULs in these regions of the EBL spectrum are relatively uncommon (Biteau & Williams 2015).

# Appendix B: Systematic analysis of the high-state selection

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

**Fig. B.1.** TS distribution for the PLxEBL (blue) and LPxEBL (purple) models, incorporating the finke2022 EBL model, is analyzed in relation to datasets defined by the flux distribution above 1 TeV. Each dataset encompasses observations with flux surpassing a specified percentile. The lower panel presents the distribution of $\sqrt{TS}$ denoted by black crosses, while the solid line illustrates the dataset's livetime distribution.

The choice of the high-state dataset significantly influences the outcomes of the final spectral fitting results. To explore this influence, we established nine distinct datasets. Each dataset comprises observations with a flux above 1 TeV specifically exceeding a predefined percentile ranging from 0% to 90%, with intervals of 10%. Subsequently, we applied fits to these datasets using PLxEBL and LPxEBL models, employing the finke2022 EBL model. The results, depicting the likelihood distribution and $\sqrt{TS}$ of the curved model in comparison to the PL model, are presented in Fig. B.1. Notably, we observe the detection of curvature at a significance level exceeding $4\sigma$, commencing with datasets featuring flux above the 30% percentile and extending to those above the 90% percentile. Moreover, the curvature is evident at a $5\sigma$ level in datasets ranging from flux above the 30% percentile to those above the 60% percentile. Consequently, we consider the measurement of curvature in the highest flux observation runs to be robust. Variations in $\sqrt{TS}$ stem from both event statistics, which systematically decrease with higher cuts, and the presence of low-state runs in the dataset. The latter is particularly relevant for datasets up to the 40% percentile cut.

Figure B.2 displays the best-fit LPxEBL model for diverse datasets categorized by their percentile cuts. The spectrum derived from the entire dataset (0% percentile) suggests a tendency for curvature, which becomes more pronounced in datasets focusing solely on observations with the highest fluxes. These outcomes reinforce the previously reported detection of
<table><caption>**Table B.1.** Best fit results for the spectral VHE gamma-ray distribution of the high state of M 87. In addition to the models presented in Table 1, we include the best-fit PLxECxEBL *finke2022* for the reduced energy range and for full energy range. The energy range and the critical energy are given in TeV and $\phi_0$ is given in units of $10^{-13} \text{cm}^{-2} \text{s}^{-1} \text{TeV}^{-1}$. The decorrelation energy of the PL spectral model was used as the reference energy, given in TeV in the table.</caption><thead><tr><th>Model</th><th>-2lnL</th><th>φ<sub>0</sub></th><th>Γ</th><th>E<sub>c</sub></th><th>Ref. energy</th><th>Energy range</th></tr></thead><tbody><tr><td>PLxECxEBL finke2022</td><td>18.0</td><td>11±3</td><td>1.02±0.26</td><td>3±1</td><td>1.96</td><td>0.3 – 10</td></tr><tr><td>PLxECxEBL finke2022</td><td>20.8</td><td>5.9±0.9</td><td>1.45±0.16</td><td>8±3</td><td>2.15</td><td>0.3 – 32</td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>**Fig. B.2.** Best spectral fits for the LPxEBL models are illustrated for the various datasets, revealing the emergence of a curved spectrum in the high state.</figcaption></figure>

curvature in the high-state VHE gamma-ray emissions of M 87 in the main text.

To conclude the systematic checks, we conducted an analysis of the high state as defined in H.E.S.S. Collaboration (2023), utilizing a Bayesian block analysis. It is worth noting that this definition may not be ideal for studying high states, given that flares are smoothed out in the presence of nearby low emission states. However, despite this limitation, the spectral fit based on the high state definition exhibits a preference of $3.6 \sigma$ for an LPxEBL *finke2022* model compared to the PLxEBL *finke2022* model. Importantly, this result is within $1 \sigma$ consistency with the findings presented in the main text.

## Appendix C: The spectrum according to an exponential cut-off function

The emergence of both the Klein-Nishina regime and the spectral region close to the maximum particle energy achievable by the accelerator might give rise to an exponential cut-off gamma-ray distribution in a synchrotron-self Compton scenario, potentially resembling curvature. To explore this possibility, we define a PL with an exponential cut-off model including the EBL absorption (PLxECxEBL *finke2022*):

$$ \phi(E) = \phi_0 \left(\frac{E}{E_0}\right)^{-\Gamma} \exp\left[-\left(\frac{E}{E_c}\right)^{\beta_c} - \alpha_{\text{norm}}\tau(E, z)\right], \quad (C.1) $$

where $E_c$ is the critical energy and $\beta_c$ determines the steepness of the cut-off.

We fitted the high state data as defined in the main text to Eq. C.1 with fixed $\alpha_{\text{norm}} = 1$ and fixed $\beta_c = 1$ both for the reduced (0.3 - 10 TeV) and for the complete (0.3 - 32 TeV) energy range. The results are presented in Table B.1. In the reduced energy range, the fit statistic of the best-fit PLxECxEBL model is indistinguishable from the best-fit LPxEBL model (Table 1). However, the spectral index of $\Gamma=1.0 \pm 0.3$ appears exceptionally hard for an intrinsic PL distribution and likely non-physical. In the full energy range, despite the more realistic best-fit spectral index, the fit statistic of the PLxECxEBL model is marginally poorer than that of the LPxEBL *finke2022* model in the same energy range, implying that a cut-off with $\beta_c = 1$ is not the preferred explanation for the detected curvature in the VHE gamma-ray high state of M 87. We refrain from fitting more complex models, that is, a PLxECxEBL model with free $\beta_c$, as such models inevitably lead to degeneracy and inconclusive results for the current dataset with limited event statistics.