© The Authors 2024

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 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 M 87, 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 M 87. Firstly, internal absorption, namely: the absorption of gamma rays at the core of M 87 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 M 87 (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 τ(E, z) from internal γ-γ absorption reaches τ ≳ 1 at ≳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 ± 0.5)×10⁹ M_⊙, where M_⊙ 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 M 87 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 Γ ∼ 2–3. Given that the optical depth τ increases monotonically with the gamma-ray energy, efforts to investigate its influence on the gamma-ray spectrum at the highest energies are restricted to bright and nearby VHE gamma-ray sources.

Due to its close proximity to Earth, M 87 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 M 87 at ≳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 (Γ ≈ 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 M 87 up to the highest energies. Therefore, we selected high-state data sets from H.E.S.S. observations of M 87 and search for absorption imprints in it. The VHE gamma-ray observations of M 87 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 M 87 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 M 87 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 (CT 1–4), which have been in operation since the earliest detected flare. Given that the 28 m diameter mirror telescope (CT 5) has only been in operation since 2012, it has not observed the aforementioned flaring states of M 87. Hence, we chose to exclude its data from this study, ensuring a consistent representation of M 87 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 ≤0.11° around the position of the core of M87 at right ascension (RA) of 187.7059° and declination (Dec) 12.3911° for J2000 (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⁻¹³ cm⁻² s⁻¹ 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).

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.

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 Multiwavelength 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 ℒ = ∏_iP(n_i|μ_i), and P(n_i|μ_i) as the Poisson probability of observing n_i events in energy bin, i, given the predicted number of events, μ_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, ϕ(E), as a function of energy, E, each with increasing complexity. The first model is a power law (PL) spectrum, defined as:

$$\begin{array}{c}\hfill \varphi (E)={\varphi }_0{\left(\frac{E}{E_0}\right)}^{-\mathrm{\Gamma }},\end{array}$$(1)

where the normalization, ϕ₀, and the spectral index, Γ, are free parameters, while the reference energy, E₀, 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:

$$\begin{array}{c}\hfill \varphi (E)={\varphi }_0{\left(\frac{E}{E_0}\right)}^{-\mathrm{\Gamma }-\beta \ln (E/E_0)},\end{array}$$(2)

where the coefficient is β, the normalization is ϕ₀, and the spectral index, Γ, 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 Domínguez et al. (2011). Collectively, these three models capture the diverse models from literature regarding the optical depth, τ.

The models characterizing the EBL spectrum, presented as the optical depth, τ, 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 ≈ 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:

$$\begin{array}{c}\hfill \varphi (E)={\varphi }_0{\left(\frac{E}{E_0}\right)}^{-\mathrm{\Gamma }}exp(-{\alpha }_{\text{norm}}\tau (E,z)),\end{array}$$(3)

$$\begin{array}{c}\hfill \varphi (E)={\varphi }_0{\left(\frac{E}{E_0}\right)}^{-\mathrm{\Gamma }-\beta \ln (E/E_0)}exp(-{\alpha }_{\text{norm}}\tau (E,z)),\end{array}$$(4)

where τ(E, z) is the model-dependent EBL optical depth, parameterized by the photon energy and redshift, and α_norm is the EBL normalization. While τ(E, z) is defined by the EBL model, α_norm can either be set to vary freely, in case one is interested in constraining the EBL, or be fixed to the expected value α_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 (−2lnℒ). 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 = −2ln(ℒ₀/ℒ₁). For large event statistics, TS follows a χ² 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{\mathit{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σ 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  ±  28 excess gamma-ray events and 376 ± 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σ. For energies above 10 TeV, we find 15 ± 4 gamma-ray excess events for 2.3 ± 0.4 expected background events, which leads to a significant detection of 6.0σ.

Given the distance to M 87, we find that the optical depth τ < 0.2 holds for energies < 10 TeV for all three considered EBL models. We considered the effects of EBL absorption for τ < 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. 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σ. We further discuss the implications of this curvature in the following section.

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σ uncertainty contours. The upper limit (UL) in the last bin of the low-state spectrum is defined at a 95% confidence level (c.l.).

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 α_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σ, 4.2σ, and 3.6σ 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  ±  0.5)×10⁻¹² cm⁻² s⁻¹ for the finke2022 LPxEBL model.

Lastly, we evaluated the LPxEBL spectral model, as previously defined, in conjunction with the same three EBL models and with the EBL normalization α_norm allowed to vary freely. Upon fitting the LPxEBL with a variable α_norm to the full energy range dataset, the fit converged to a pure LP model. This suggests a degeneracy between the curvature parameter β and the EBL normalization α_norm.

To elucidate the observed degeneracy, we present the log-likelihood values for the parameter space of α_norm and β within the LPxEBL model. The result for the finke2022 EBL model is shown in terms of $\sqrt{\mathit{TS}}$ from the best fit position in Fig. 3. Based on the intersection of the α_norm with the 2σ contour in Fig. 3, we derive the maximum allowed EBL normalization α_norm within 95% c.l. of α_norm < 5.5, for the finke2022 EBL model. Similarly, for the kneiske, and dominguez-upper EBL models, we derive α_norm < 8.7 and α_norm < 2.0, respectively. In the case of β = 0, the best-fit position in the parameter space is always more than 4σ away from the fit for an EBL intensity for α_norm ≲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 ≈0.4 μm–40 μm EBL photons¹, although our analysis is more sensitive for gamma rays ≳10 TeV, i.e., for EBL photons ≳12.4 μm, where their optical depth reaches τ ≳ 0.2.

Fig. 3. Parameter space around the best-fit LPxEBL spectral model for the αnorm and β 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σ, 2σ, 3σ, and 4σ uncertainty contours. The vertical green lines show the αnorm = 1 and αnorm = 5.5, which provides the 95% c.l. UL on the α parameter. For comparison, the dotted red line shows the αnorm that corresponds to the UL above which the finke2022 EBL model would start to contradict the ULs from Biteau & Williams (2015).

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 PLxECxEBL 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 (α_norm) and the intrinsic curvature coefficient from the LP model (β) 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 α_norm is expected to be at least close to one and, therefore, it cannot be set to zero in expense of β as may be suggested by Fig. 3. It is clear through the 4σ statistical uncertainty contours from Fig. 3 that the curvature in the spectrum is strongly detected at a 4σ level even by considering α_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 (Γ ≈ 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 the selected high state data (Fig. 1), separately. The results show a hint of 1.5σ and 2.9σ 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 ± 0.09 and 1.9 ± 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 ≈ 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, Γ = 2.63 ± 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 ≈ 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 MeV < E < 300 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 ≈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 ≲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 μm–1.7 μm range, denoted as μ_λ = 15.51 mag arcsec⁻² from Capetti & Balmaverde (2006). Using Eq. (A5) of Zacharias et al. (2017), we convert the surface brightness to emissivity j_λ:

$$\begin{array}{c}\hfill j_{\lambda }=2.32\times {10}^{-5}{\left(\frac{r_0}{\mathrm{cm}}\right)}^{-1}{10}^{0.4\times (25.992-({\mu }_{\lambda }-A))}\left(\frac{\mathrm{erg}}{{\mathrm{cm}}^3{\mathrm{s}}^1}\right),\end{array}$$(5)

where we set the galactic extinction A to zero, as it has already been corrected for the μ_λ estimate and r₀ is the characteristic radius of the source. We used r₀ ≈ 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 τ_M87, λ for the galactic emission up to r₀ centered at λ is estimated through:

$$\begin{array}{c}\hfill {\tau }_{\mathrm{M}87,\lambda }=\frac{{\sigma }_{\gamma \gamma }{j}_{\lambda }\lambda r_0^2}{h{c}^2},\end{array}$$(6)

where σ_γγ is the cross section for γ − γ 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 cosinus 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 λ = 1.6 μm peaks for a gamma-ray energy of ≈1 TeV, although τ_M87, λ < 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, τ < 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σ) 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 ≈ 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 confirme that the curvature in the spectrum is also strongly detected considering the full energy range as the LPxEBL models are preferred at ≈4σ over the PLxEBL models. The total flux resulted from the analysis above 0.3 TeV is estimated to be (5.1 ± 0.5) × 10⁻¹² cm⁻² s⁻¹ for the finke2022 EBL model. The curvature is characterized by a spectral index Γ = 1.80 ± 0.08 and a curvature coefficient of β = 0.27 ± 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 (α_norm = 1). Furthermore, we derive the maximum allowed EBL normalization, α_norm, within 95% confidence level of α_norm < 8.7, α_norm < 5.5, α_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.

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

Appendix A: EBL models and derived ULs

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.

Fig. A.1. Optical depth (top) and the relative intensity of the gamma-ray flux after absorption (bottom) at z≈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.

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 α_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 showcased 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).

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.

Appendix B: Systematic analysis of the high-state selection

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{\mathit{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 σ, 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 σ 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{\mathit{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.

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{\mathit{TS}}$ denoted by black crosses, while the solid line illustrates the dataset’s livetime distribution.

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 curvature in the high-state VHE gamma-ray emissions of M 87 in the main text.

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.

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 σ for an LPxEBL finke2022 model compared to the PLxEBL finke2022 model. Importantly, this result is within 1 σ 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):

$$\begin{array}{c}\hfill \varphi (E)={\varphi }_0{\left(\frac{E}{E_0}\right)}^{-\mathrm{\Gamma }}exp[-{\left(\frac{E}{E_c}\right)}^{{\beta }_c}-{\alpha }_{\text{norm}}\tau (E,z)],\end{array}$$(C.1)

where E_c is the critical energy and β_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 α_norm = 1 and fixed β_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 Γ = 1.0 ± 0.3 appears exceptionally hard for an intrinsic PL distribution and likely nonphysical. 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 β_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 β_c, as such models inevitably lead to degeneracy and inconclusive results for the current dataset with limited event statistics.

All Tables

All Figures

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

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σ uncertainty contours. The upper limit (UL) in the last bin of the low-state spectrum is defined at a 95% confidence level (c.l.).

In the text

Fig. 3. Parameter space around the best-fit LPxEBL spectral model for the αnorm and β 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σ, 2σ, 3σ, and 4σ uncertainty contours. The vertical green lines show the αnorm = 1 and αnorm = 5.5, which provides the 95% c.l. UL on the α parameter. For comparison, the dotted red line shows the αnorm that corresponds to the UL above which the finke2022 EBL model would start to contradict the ULs from Biteau & Williams (2015).

In the text

Fig. A.1. Optical depth (top) and the relative intensity of the gamma-ray flux after absorption (bottom) at z≈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.

In the text

	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 the text

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{\mathit{TS}}$ denoted by black crosses, while the solid line illustrates the dataset’s livetime distribution.

In the text

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