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

The radio galaxy M87 is a variable very-high energy (VHE) gamma-ray source, exhibiting three major flares reported in 2005, 2008, and 2010. Despite extensive studies, the origin of the VHE gamma-ray emission is yet to be understood. In this study, we investigate the VHE gamma-ray spectrum of M87 during states of high gamma-ray activity, utilizing 20.2 hours the 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.3TeV at the 4 σ 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 mainly sensitive in the wavelength range 12.4 µ m – 40 µ m.


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 (we refer to 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 H.E.S.S. Collaboration (2006a), Acciari et al. (2009), and 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 helps shed light on the radiation mechanisms in blazars.Various mechanisms have been proposed to explain the production of veryhigh energy (VHE, E>100 GeV) gamma-ray flares from M 87.These include, e.g., leptonic models (for a review, Rieger & Levinson 2018), magnetospheric models (Ansoldi et al. 2018), multi-blob synchrotron self-Compton (SSC; Lenain et al. 2008), photo-hadronic model (Sahu et al. 2019;Alfaro et al. 2022), 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 lead 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 might 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 identify three photon fields (and regions) that could significantly alter the VHE gamma-ray spectrum of M 87.Firstly, internal absorption, i.e., the absorption of gamma rays at the core of M 87 by the photons emitted in the accretion region (Brodatzki et al. 2011).Secondly, 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, therefore, shedding light on the origin of the emission and on the mechanism of particle acceleration.
One approach to investigate spectral imprints is to examine the γ-γ absorption of TeV photons by photon fields from the region.Brodatzki et al. (2011) suggest that the theoretical optical depth τ(E, z) from internal γ-γ absorption reaches τ≳1 at ≳20 TeV, assuming for that the VHE gamma rays are produced within 25r S from the central black hole, where r S is the Schwarzschild radius of the super-massive black hole (SMBH).This estimate assumes a Shakura-Sunyaev disk and a mass of the SMBH of (6.4 ± 0.5) × 10 9 M ⊙ , where M ⊙ is the solar mass, which 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), i.e., 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, one 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.
The VHE flux from gamma-ray sources quickly falls off at higher energies, decreasing typically as a power law (PL) with photon index Γ∼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 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 to study the VHE gamma-ray spectrum of M 87 up to the highest energies.Therefore, we select 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.No significant evidence for such oscillations was found.
This publication is organized as follows: in Sect. 2 we introduce the H.E.S.S. experiment, 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, before summarizing them in Sect. 5.

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 minutes.This study utilizes H.E.S.S. observations of M 87 from 2004 to 2020.Although some of the data were already used in previous publications (H.E.S.S. Collaboration 2006a; Acciari et al. 2009;Abramowski et al. 2012; EHT MWL Science Working Group 2021; H.E.S.S. Collaboration 2023), we focus here, for the first time, specifically on the high state, combining data from individual observational runs with elevated flux.To have a consistent data set, we only consider data from the 4 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) is in operation only since 2012, it has not observed the aforementioned flaring states of M 87.Hence, we choose to exclude its data from this study, ensuring a consistent representation of M 87 in VHE gamma rays over time.We have applied the quality criteria for a spectral analysis (H.E.S.S. Collaboration 2006b), selecting only the best quality observations with at least three telescopes participating in the observation.
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 cm −2 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 select 10% of the runs with the highest flux, i.e., 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 hours 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 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 use the likelihood L= i P(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: where the normalization ϕ 0 and the spectral index Γ are free parameters, while the reference energy E 0 is fixed to the decorrelation energy, i.e., to the energy at which the correlation between the spectral model parameters is minimized.
To account for potential curvature in the spectrum, we consider a Log-Parabola (LP) spectral model: where the coefficient β, the normalization ϕ 0 , and the spectral index Γ are 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 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: 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 utilize 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 lnL).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 T S = −2 ln(L 0 /L 1 ).For large event statistics, T S follows a χ 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

√
T S .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, though 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.

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 hours.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 consider the effects of EBL absorption for τ< 0.2 as small and focus 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.
In a second step, we consider 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 −12 cm −2 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 in the last bin of the low-state spectrum is defined at a 95% confidence level (c.l.).
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 loglikelihood values for the parameter space of α norm and β within the LPxEBL model.The result for the finke2022 EBL model is shown in terms of

√
T S 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 β=0, the bestfit 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 1 , 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.
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 PLx-ECxEBL model, taking into consideration the effects of the EBL.1: Best fit results for the spectral VHE gamma-ray distribution of the high state of M 87.We compare 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 ϕ 0 is given in units of 10 −13 cm −2 s −1 TeV −1 .The decorrelation energy of the PL spectral model was used as the reference energy, given in TeV in the table.

Discussion and conclusions
A curved gamma-ray spectrum has been already observed up to TeV energies in the radio galaxy NGC1275 (Ansoldi et al. 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 consider four possibilities and discuss them in the subsequent paragraphs based on the details of our analysis: -External absorption: EBL photons; -Multi-component high state: a composition of different particle populations; -Internal absorption: by local photon fields, e.g., the accretion flow and the star light; -Intrinsic curved electron spectrum: due to the nature of the mechanism of acceleration, e.g., 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 Fig. 3 might suggest.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, examine but do not find a strong correlation between the location of voids and hard gamma-ray sources.However, they 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.09and 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.3TeV).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 employ the surface brightness of M 87 (UGC 7654) in the 1.4 µm -1.7 µm range, denoted as µ λ =15.51 mag/arcsec 2 from Capetti & Balmaverde (2006).Using Eq.A5 of Zacharias et al. (2017), we convert the surface brightness to emissivity j λ : where we set the galactic extinction A to zero, as it has already been corrected for the µ λ estimate and r 0 is the characteristic radius of the source.We use r 0 ≈ 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 0 centered at λ is estimated through: 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, i.e. 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 gammaray 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).

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 hours of observations.
We first consider the energy range between 0.3 and 10 TeV, in which the EBL absorption is small, i.e., τ<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 confirm 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)x10−12 cm −2 s −1 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, i.e. 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. .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

√
T S 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.  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 ϕ 0 is given in units of 10 −13 cm −2 s −1 TeV −1 .The decorrelation energy of the PL spectral model was used as the reference energy, given in TeV in the   .2displays 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.

Model
To conclude the systematic checks, we conducted an analysis of the high state as defined in H.E.S.S. Collaboration (2023), utilizing 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 gammaray 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): 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 gammaray high state of M 87.We refrain from fitting more complex models, i.e. a PLxECxEBL model with free β c , as such models inevitably lead to degeneracy and inconclusive results for the current dataset with limited event statistics.
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.).

Fig. 3 :
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).
Fig. A.2:The EBL photon distribution based on the three models employed in this studykneiske, 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.

Fig
Fig. B.2:The best spectral fits for the LPxEBL models are illustrated for the various datasets, revealing the emergence of a curved spectrum in the high state.

Fig. B
Fig. B.2 displays the best-fitLPxEBL 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.To conclude the systematic checks, we conducted an analysis of the high state as defined in H.E.S.S. Collaboration (2023), utilizing 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.

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