Very high energy gamma-ray emission beyond 10 TeV from GRB 221009A

One sentence summary: LHAASO detected the highest energy gamma-rays, beyond 10 TeV, from the brightest gamma-ray burst GRB 221009A.

Introduction

Gamma-ray bursts (GRBs) are sudden explosions of gamma-rays that occur in random directions from cosmological distances. They are the most luminous fireworks in the Universe and are thought to originate from collapsing massive stars and compact star mergers. Prompt flashes and long-lasting afterglows have been thoroughly studied, with thousands of GRBs observed in a wide range of energies from radio waves to megaelectronvolt (MeV) gamma-rays. A fraction of GRBs are also observed with gigaelectronvolt (GeV) gamma-rays. Recently, a handful of GRBs have been observed in teraelectronvolt (TeV) gamma-rays during the afterglow period [1, 2, 3, 4]. The synchrotron self-Compton (SSC) process of relativistic electrons in the afterglow has been proposed to explain the origin of the TeV emission [2, 5, 6]. However, challenges [4] have been reported for such a scenario. Observation of the highest energy emission from GRBs is important to probe the gamma-ray emission mechanism and particle acceleration processes in GRBs.

The exceptionally bright GRB 221009A ($\mathrm{RA}=288.264^{\circ}$, $\mathrm{Dec}=19.768^{\circ}$) was detected by Fermi-GBM on October 9, 2022, at 13:16:59.99 UT (denoted as $T_0$ hereafter) [7], about one hour earlier than the Swift trigger [8]. The GBM light curve consists of two emission episodes: a single isolated peak from $T_0$ to $T_0+20$s, followed by a longer, extremely bright, multi-pulsed emission episode from about $T_0+220$s to $T_0+550$s [7]. It was also detected by the Large Area Telescope (Fermi-LAT) with the highest-energy gamma-ray of 99.3 GeV observed 240 seconds after the GBM trigger [9]. This GRB was also clearly observed by many other detectors, such as Insight-HXMT, GECAM-C [10], and Konus-Wind [11]. Follow-up optical observations revealed that the redshift of this GRB is $z=0.151$ (corresponding to a distance $\sim753$ Mpc) [12]. The estimated isotropic energy release is about $10^{55}$ ergs [7, 10, 11], which denotes it as an extremely energetic GRB.

Results

LHAASO observed GRB 221009A at the highest energy band. LHAASO [13] consists of three interconnected detectors located at 4,410 m above sea level in Sichuan Province, China. The sub-arrays, a 78,000 m$^2$ Water Cherenkov Detector Array (WCDA) and a 1.3 km$^2$ Kilometer Squared Array (KM2A), are dedicated to gamma-ray observations with a wide field-of-view (FOV). Both arrays can monitor the sky for zenith angles less than 50°. GRB 221009A entered the FOV of LHAASO at around $T_{0}$-20000s, culminating at a zenith angle of $9.5^\circ$ at $T_{0}$-7000s, and left the FOV at $T_{0}$+6000s. For the first time for a TeV-detected GRB, the observations cover both the prompt period and the afterglow period. The zenith angles were $28^\circ$ at $T_{0}$ and $31.5^\circ$ at $T_{0}$+1000s, which are favorable positions for LHAASO observation. Within 2000s after $T_{0}$, more than 5000 gamma-rays at $>0.5$ TeV were detected by WCDA from GRB 221009A with significance above 100 s.d. (14). A previous report has focused on the temporal characteristics of the TeV emission based on WCDA data at lower energies, which revealed the earliest TeV afterglow from GRB 221009A (15). The light curve includes a sharp rise at $T_{0}$+230s reaching the peak at about $T_{0}$+245s and a following smooth decay with the jet break. Here, we focus on the spectral results extending to energies above 10 TeV, as measured by KM2A.

KM2A is optimized for gamma-rays with energies from 10 TeV to 10 PeV (16). It also has sensitivity to detect gamma-rays below 10 TeV with an effective area of 10,000 $\mathrm{m}^{2}$ at 4 TeV and 100,000 $\mathrm{m}^{2}$ at 7 TeV, which overlaps with WCDA. Below 10 TeV, KM2A can still remove 98% of the cosmic ray background using a “muon-less” content criterion. The angular resolution is about $1^\circ$ and the energy resolution is around 40%. The count-rate light curve observed by KM2A is shown in the left panel of Figure 1. The emission started at $T_{0}$+230s with flux peak around $T_{0}$+245s. The emission gradually decreased since $T_{0}$+300s and faded out after $T_{0}$+900s. During the period from $T_{0}$+230s to $T_{0}$+900s, 142 events with energies above 3 TeV from the direction of GRB 221009A were detected, and the cosmic ray background was estimated to be 16.7. Therefore, the gamma-ray emission was detected with a significance of 20.6 $\sigma$. The position of the GRB was estimated to be $\mathrm{RA}=288.26^\circ\pm0.07^\circ$, $\mathrm{Dec}=19.83^\circ\pm0.07^\circ$ using the KM2A data, which is consistent with that of Fermi-LAT. The significance map can be found in the right panel of Figure 1.

The KM2A data are divided into two intervals for spectral analysis, according to the light curve shown in Figure 1, i.e., interval 1 from $T_{0}$+230s to $T_{0}$+300s and interval 2 from $T_{0}$+300s to $T_{0}$+900s. We combined the WCDA and KM2A data using a joint forward-folded fit to determine the spectral energy distribution (SED) of GRB 221009A, which is shown in Figure 2, assuming a log-parabola (denoted as LP) or a power-law with exponential cutoff (denoted as PLEC) function. WCDA covers the energy range from 0.2 TeV to 7 TeV, and KM2A covers the energy range from 3 TeV to 20 TeV. Details about the analysis of KM2A data are presented in the Materials and methods section. More details about the analysis of WCDA data can be found elsewhere (15). The two measurements are consistent with each other in the overlapping energies from 3 TeV to 7 TeV. Fitting an LP function, the yielded $\chi^{2}/ndf$, where $ndf$ is the number of degrees of freedom, are 14.1/9 and 37.3/10, respectively. The probability for the second interval is very low which indicate the LP function is not favored. If fitting an PLEC function, the yielded $\chi^2/\mathit{ndf}$ become much better, with 10.1/9 and 18.3/10 for the two time intervals, respectively. Therefore, the fitting using the PLEC function is better than that using the LP function. The detailed parameters yielded by these two functions for the two intervals are listed in Table 1. The spectrum at interval 2 is slightly harder than that of interval 1. Besides the statistic uncertainty, an addition 7% systematic uncertainty should also exist for the flux according to (16). It is worth noting that the precise measurement by LHAASO is the first detection of a cutoff at the high-energy end of a gamma-ray spectrum from a GRB.

PLEC	$J_0(10^{-8}\ \mathrm{TeV}^{-1}\ \mathrm{cm}^{-2}\ \mathrm{s}^{-1})$	$a$	$E_{\mathrm{cut}}\ (\mathrm{TeV})$	$\chi^2/\mathrm{ndf}$	Probability
230-300s	$51.9\pm3.3$	$2.87\pm0.05$	$2.62\pm0.26$	$10.1/9$	$0.34$
300-900s	$10.5\pm0.6$	$2.65\pm0.05$	$1.99\pm0.15$	$18.3/10$	$0.050$
LP			$b$
230-300s	$35.7\pm1.3$	$3.46\pm0.03$	$0.58\pm0.05$	$14.1/9$	$0.12$
300-900s	$6.37\pm0.16$	$3.38\pm0.02$	$0.73\pm0.05$	$37.3/10$	$5.0\times10^{-5}$
Saldana-Lopez EBL
230-300s	$144\pm4$	$2.35\pm0.03$	$0\ (\text{fixed})$	$11.0/10$	$0.36$
300-900s	$24.5\pm0.5$	$2.26\pm0.02$	$0\ (\text{fixed})$	$6.6/11$	$0.83$
LHAASO constrained EBL
230-300s	$214\pm6$	$2.12\pm0.03$	$0\ (\text{fixed})$	$5.9/10$	$0.82$
300-900s	$37.8\pm0.9$	$2.03\pm0.02$	$0.15\pm0.06$	$5.5/10$	$0.86$

With the KM2A energy resolution of about 40% at 10 TeV, the reconstruction of energy of an event is non-trivial, in particular for very steep spectra, e.g. in the region of an energy cutoff. To correctly estimate the true energy and corresponding errors of each high energy event, the probability function of the true energy is constructed for each event, which strongly depends on the assumption of the spectral function. More details can be found in the material and method section.

Using the LP function spectra shown in panel A of Figure 2, eight events with reconstructed energy above 10 TeV were observed during the period from $T_0$+230s to $T_0$+900s. The event with the maximum energy is at $17.8^{+7.4}_{-5.1}$ TeV. This energy is similar with that reported in the GCN (14), which is preliminary estimated using the KM2A alone data assuming a power-law function spectrum. However, if using the PLEC function spectra shown in panel B of Figure 2, the energy of the maximum energy event is $12.2^{+3.5}_{-2.4}$ TeV. If the spectral shape is governed by the strongly energy-dependent absorption of gamma-rays in collisions with photons of the Extragalactic Background Light (EBL), which will be discussed later, the highest energy photon is reconstructed at $12.5^{+3.2}_{-2.4}$ TeV. The energies of each event shown in Figure 1 were obtained under the PLEC function spectra. The energies of the nine gamma-ray-like events with the highest energies have been estimated under three spectra, i.e., LP, PLEC, and an EBL model. It is worth to note that the LP spectrum is not favored by the data.

There is furthermore a small probability that an event is a misidentified cosmic-ray event; the chance probability has been estimated for each event individually and varies for the nine events between 4.5% and 17%. More details can be found from the Table S2 in the Supplementary Materials section. Gamma-ray events with energy above 100 TeV were searched for during a long period from $T_0$ to $T_0$+6000s; however, no event was detected.

$T_{\mathrm{event}}$(s)	$E_{LP}$ (TeV)	$E_{PLEC}$ (TeV)	$E_{EBL}$ (TeV)	$N_{e}$	$N_{\mu}$	$\theta$ (°)	$\Delta\psi$ (°)	$D_{edge}$ (m)	P (%)
236.6	$12.7^{+6.2}_{-3.8}$	$9.7^{+3.3}_{-2.1}$	$9.8^{+3.1}_{-2.3}$	60.6	0	28.5	0.46	77	7.0
242.5	$10.5^{+5.0}_{-3.2}$	$8.3^{+3.0}_{-2.1}$	$8.4^{+3.2}_{-2.2}$	57.4	0	28.8	0.45	111	10
262.4	$12.6^{+5.5}_{-3.8}$	$9.5^{+3.4}_{-2.3}$	$9.6^{+3.3}_{-2.4}$	57.3	0	28.6	0.53	180	5.7
358.1	$10.0^{+4.8}_{-3.2}$	$7.4^{+3.1}_{-1.8}$	$7.9^{+3.3}_{-2.2}$	46.0	0	28.7	0.54	119	6.0
571.1	$9.4^{+5.1}_{-3.0}$	$7.4^{+2.6}_{-2.5}$	$7.7^{+3.0}_{-2.5}$	45.7	0	29.5	0.52	99	7.8
643.0	$17.8^{+7.4}_{-5.1}$	$12.2^{+3.5}_{-2.4}$	$12.5^{+3.2}_{-2.4}$	81.8	0.3	29.7	0.62	181	4.5
812.4	$11.1^{+5.9}_{-4.3}$	$7.4^{+3.6}_{-2.8}$	$7.6^{+3.9}_{-3.0}$	68.0	0	30.3	0.66	112	11
863.8	$12.9^{+6.1}_{-3.9}$	$9.2^{+3.0}_{-2.8}$	$9.7^{+3.2}_{-3.1}$	100.2	0.8	30.1	1.07	81	17
894.1	$13.6^{+6.1}_{-4.2}$	$9.7^{+3.4}_{-2.5}$	$10.4^{+3.3}_{-3.0}$	60.5	0	31.8	0.83	214	16

EBL absorption of emission from GRB 221009A

It is known that very high-energy (VHE) gamma-rays emitted from distant astronomical sources are affected by the EBL absorption through photon-photon interaction, inducing a cutoff in the energy spectrum in the form of $e^{-\tau(E)}$ with $\tau(E)$ being the energy-dependent optical depth. For GRB 221009A at a redshift of $z=0.151$, the absorption at 1 TeV is modest with a survival fraction ranging from 18% to 21%, while it becomes very heavy at 10 TeV ranging from 0.5% to 0.05%, depending on the EBL models (see Figure S1). The cutoff energies in the two intervals measured by LHAASO are consistent with each other within a 2.1$\sigma$ error range, which is expected due to EBL absorption. Therefore, the spectra of GRB 221009A measured by LHAASO provide an excellent test of EBL models.

To quantitatively test different EBL models, we utilize a log-parabolic function to describe the intrinsic GRB spectrum and then fit the EBL-attenuated gamma-ray spectrum in the form $dN/dE = J_0 E^{a+b\cdot\log(E)}e^{-\tau(E)}$ to the data. The fitting results are given in the Supplementary Materials section. The goodness of fit combining the two intervals is comparable for different EBL models, i.e., Saldana-Lopez et al. 2021 (17), Gilmore et al. 2012 (18), Dominguez et al. 2011 (19), and Finke et al. 2010 (20), as shown in Table S3. The intrinsic GRB spectrum obtained using the recent EBL model of Saldana-Lopez et al. 2021 (17) is shown in panel A of Figure 3. The spectral curvature is insignificant and therefore a simple power-law (PL) function is used. The detailed parameters yielded by the fitting for the two intervals are listed in Table 1.

EBL model	230-300s	300-900s	total
Saldana-Lopez et al. 2021	11.02/9	5.44/10	16.46/19
Gilmore et al. 2012	3.53/9	15.29/10	18.82/19
Dominguez et al. 2011	4.16/9	11.33/10	15.49/19
Finke et al. 2010	6.12/9	13.51/10	19.63/19
LHAASO constrained	5.93/-	5.49/-	11.42/16

The intrinsic SED of the two intervals follows a single power-law model with an index of -(2.35$\pm$0.03) for interval 1 and -(2.26$\pm$0.02) for interval 2. For interval 1, the $\chi^2$/ndf of the fit is 11.0/10, and the maximum deviation is at 5.7 TeV, which deviates from the fit line by 2.3$\sigma$. For interval 2, the $\chi^2$/ndf of the fit is 6.6/11, while a clear deviation appears at the highest energy point. The measured flux is 11 times higher than the fit line, however, the deviation is not statistically significant due to large flux error (panel A of Figure 3). This situation is similar to that obtained for the blazar Mrk501 during its huge flare in 1997 (21). Even though some hypotheses were proposed to produce a non-classical sharp pile-up spectrum (22), a straightforward approach is to reduce the deviation at the highest energy point by decreasing the optical depth at high energy. An empirical method would be to decrease the EBL intensity, the measurement of which still contains large uncertainties due to the heavy contamination caused by foregrounds of different origin, predominantly by the zodiacal (interplanetary dust) light.

The absorption of gamma-rays at a specific energy $E$ is mainly due to the EBL photons within a narrow interval with $\Delta\lambda\sim(1 \pm1/2)\lambda$ centered on $\lambda\approx1.5(E/1\mathrm{TeV})\,\mu\mathrm{m}$ (23). Therefore, fore, the measured gamma-ray spectra from GRB 221009A can be used to constrain the EBL intensities at different wavelengths. For this purpose, the EBL model of Saldana-Lopez et al. 2021 (17) is independently re-scaled at three wavelength ranges, i.e., $\lambda< 8\mu\mathrm{m}$, $8\mu\mathrm{m} < \lambda< 28\mu\mathrm{m}$, and $\lambda> 28\mu\mathrm{m}$, which make the dominant contribution to the optical depth for gamma-rays below 5 TeV, between 5 to 10 TeV, and above 10 TeV, respectively. Assuming a log-parabolic intrinsic spectral function by fitting to the LHAASO data, we obtain the re-scaled factors for the three EBL ranges as $1.30^{+0.33}_{-0.20}$, $1.20^{+0.23}_{-0.20}$, and $0.40^{+0.44}_{-0.16}$, respectively. The constrained EBL intensity is shown in Figure 4.

To check the effect of possible absolute energy uncertainty, we have shifted the energy by $\pm10\%$. Then the scale factor of the last EBL bin changes from 0.4 to 0.35 and 0.60, respectively. It should be emphasized that the EBL intensity at large wavelengths is now close to the lower bound set by galaxy count observations (24).

The intrinsic spectra using the re-scaled best-fit EBL model are shown in panel B of Figure 3. The detailed parameters yielded by the fitting for the two intervals are listed in Table 1. For interval 1, the intrinsic SED follows a single power-law model with an index of $-(2.12\pm0.03)$, with a $\chi^2/\mathrm{ndf}$ of 5.9/10 yielded by the fitting. For interval 2, the intrinsic SED favors a log-parabolic form with an index of $-(2.03\pm0.02)-(0.15\pm0.06)\log(E/\mathrm{TeV})$. The corresponding $\chi^2/\mathrm{ndf}$ of the fit is 5.5/10. We note that the flux at the highest energy point is much lower now and consistent with the fit spectrum since the EBL absorption is reduced after re-scaling EBL for wavelengths $\lambda> 28\mu\mathrm{m}$. Using the LHAASO constrained EBL model and assuming the intrinsic SED shown in panel B of Figure 3, the energy of the maximum energy event is estimated to be $12.5^{+3.2}_{-2.4}$ TeV.

EBL absorption corrected light curve from GRB 221009A

The light curve of GRB 221009A in 0.2-5 TeV, as observed by WCDA, is characterized by a rapid rise to a peak, followed by a decay that persists for at least 3000 s after the peak (15). The smooth temporal profile suggests an external shock origin for the emission, caused by the interaction of the GRB ejecta with the ambient medium. To compare the temporal variation between WCDA and KM2A emissions, the count rate light curve shown in Figure 1 is converted to the energy flux light curve assuming the spectra shown in the right panel of Figure 3. A direct comparison between the energy flux light curves observed by WCDA and that of KM2A is shown in Figure 5. The solid line indicates the best-fit four-segment light curve observed by WCDA. Obviously, the light curve observed by KM2A is similar to that of WCDA, indicating that the bulk of KM2A emission should have the same origin as that of the WCDA.

Limits on new physics using the observation of GRB 221009A

The detection of TeV $\gamma$-rays from a cosmological distance has important implications for possible new physics beyond the Standard Model (SM) of particle physics. New physics scenarios, such as axions and LIV, may affect the process of EBL absorption. To simplify the discussion, we adopt the EBL model of Saldana-Lopez et al. 2021 (17) here.

The axion and axion-like particles (ALPs) are hypothetical particles introduced by theoretical models developed beyond the SM. They are among the most attractive cold dark matter candidates. The most distinctive feature of axion is its coupling vertex with two photons, leading to axion-gamma-ray conversion in external magnetic fields for astrophysical detection. The axion-gamma-ray conversion leads to a suppression of EBL absorption of high-energy gamma-rays and higher transparency of space. We still utilize a log-parabolic form to describe the intrinsic GRB spectrum in this scenario. The $\chi^2$ of the best fitting for the spectra during the two intervals reduces by only 1.7 compared to that without axion-gamma-ray conversion. Therefore, we give constraints on the parameter space of axion coupling with 95% confidence level (C.L.). Details of the calculation are given in the Materials and methods section. In Figure 6, the constraint on the axion-gamma-ray coupling constant is shown. It shows that the constraint on the coupling constant by LHAASO observation of GRB 221009A is improved compared with the constraint set by CAST (25) (CERN Axion Solar Telescope) and is comparable to that derived from the observation of PKS 2155-304 energy spectrum by HESS (26) and a long-term monitoring of Mrk 421 by ARGO-YBJ (27).

Lorentz invariance is a fundamental principle in modern physics. However, some theories that attempt to unify quantum mechanics and general relativity may suggest LIV at the energy scale approaching the Planck scale $M_{\mathrm{Pl}} = 1.22 \times10^{28}$ eV. This may allow many interesting phenomena which are otherwise forbidden processes (28). One consequence is that the threshold of an interaction may be altered. For GRB gamma-rays, if the threshold of interacting with EBL is slightly increased, the process $\gamma\gamma\to e^{-}e^{+}$ is suppressed, and the universe becomes more transparent. In such a scenario, the detection of gamma rays beyond 10 TeV is easier to understand. However, the reduce of the $\chi^2$ of the fitting is only 2.6, which is turned to set a constraint on the LIV effect. We derive that for the first-order LIV, the effective energy scale of LIV should be greater than $1.5M_{\mathrm{Pl}}$. This constraint is comparable to that derived by exploring the energy dependence of the propagation speed of gamma-rays from GRB090510 by Fermi-LAT (29).

Discussion

In this work, we report the first detection of gamma-rays beyond 10 TeV from a gamma-ray burst, GRB 221009A, during the early afterglow period by LHAASO. The measured energy spectrum has a clear cutoff at the high-energy end, which is consistent with the standard picture of gamma-rays propagating and interacting with EBL. The derived intrinsic energy spectrum shows a power-law shape extending beyond 10 TeV, indicating a highly efficient particle acceleration mechanism working within the relativistic shock of the GRB. The LHAASO data prefer a suppression of the EBL flux at mid-IR wavelengths. These observations also place strong constraints on new physics parameters, such as an axion origin of the signal or LIV. In the following, we will discussion the implication for GRB physics and the possible reasons account for the highest energy gamma-ray.

Implication for GRB physics

Figure 3 shows that the intrinsic spectrum of GRB 221009A does not exhibit softening up to at least 10 TeV. For previous TeV afterglows, the SSC process of relativistic electrons has been proposed to produce the TeV emission [2, 5, 6]. However, since the Klein-Nishina (KN) effect becomes increasingly significant at higher energies, one would expect a softening of the SSC spectrum towards higher energies (30). This can be seen from the comparison between the model that explains the WCDA data (15) and the data of KM2A of GRB 221009A (see Figure 7). This may be similar to the challenge in explaining the TeV data of GRB 190829A with the SSC scenario [4]. Such a discrepancy would suggest more complicated processes occurring during the early afterglow phase.

One possible solution to this discrepancy for GRB 221009A is to assume that an additional hard spectral component emerges at the highest energy end. In the ”standard” paradigm of an external shock, two types of shocks are formed when the ultrarelativistic ejecta are decelerated by the swept-up ambient medium: a forward shock propagating into the external medium that produces the afterglow emission, and a reverse shock propagating back into the ejecta shell. In the external shocks, both electrons and protons are accelerated. Proton synchrotron emission from accelerated ultra-high-energy cosmic rays (UHECRs) in the external reverse shock has been suggested to produce a hard spectral component in GRB 221009A (31). The magnetic field in the reverse shock could be larger than that of the forward shock, so protons could be accelerated to ultra-high energies by the reverse shock and thus produce synchrotron emission above 10 TeV (31). For a proton spectrum $dN_p/dE \sim E^{-2}$, the spectral index of the proton synchrotron emission would be $-1.5$, much harder than the SSC emission. In addition, the intergalactic electromagnetic cascade due to the propagation of UHECRs (32) can also produce a hard spectral component at the highest energy, because $>10$ TeV gamma-rays could be produced by UHECRs that have propagated to a nearby distance to us. To explain the observed flux of $>10$ TeV gamma-rays in the hadronic scenario, both models require efficient production of UHECRs from GRBs (31, 32).

An alternative explanation could be to invoke an additional hard leptonic component. This could be realized in the multi-zone models where the magnetic field is inhomogeneous throughout the emitting volume (33). Synchrotron photons from the strong magnetic field zone provide the dominant target for cooling of the electrons in the weak magnetic field zone, which could result in a formation of a hard electron distribution due to the Klein-Nishina scattering effect. A hard electron spectral component could also be formed by the hydrodynamical turbulence that is excited in the GRB forward shock and stochastically accelerates protons and electrons (34). The stochastic acceleration can yield a hard electron spectrum with $p<2$, which has also been discussed as a possible electron acceleration mechanism in AGN jets.

Possible reasons account for the highest energy gamma-ray event

The detection of a gamma-ray event beyond 10 TeV from a source at a cosmological distance of $z=0.151$ seems unlikely, considering the heavy absorption by EBL. As shown above, it may hint a suppression of the EBL intensity at large wavelengths. Several possible reasons may account for the observation.

The first possibility is that the highest gamma-ray event is actually a cosmic ray background event with unusually low muon content. The gamma-ray and cosmic ray events recorded by LHAASO are mainly discriminated by the muon content in the air shower. Therefore, a less muonic cosmic ray event could mimic a gamma-ray event. The chance of misidentification has been estimated to be 4.5%, taking into account the arrival time, direction, and the muon content. The second possibility is that this event is actually a lower energy gamma-ray but is reconstructed to above 10 TeV because of the poor energy resolution. If the intrinsic SED follows the fitting line shown in panel A of Figure 3 and taking the EBL model of Saldana-Lopez et al. 2021 (17), the expected number of events beyond the highest energy during the interval from $T_0+300\mathrm{s}$ to $T_0+900\mathrm{s}$ is 0.009, which will be enlarged to 0.015 if using the interval from $T_0+230\mathrm{s}$ to $T_0+900\mathrm{s}$. The third possibility is that the absorption of EBL is weaker than the current models predict. As shown in panel B of Figure 3, the situation improves after decreasing the EBL intensity at wavelengths above 28 $\mu\mathrm{m}$.

Besides these, other interesting possibilities that avoid EBL absorption may be due to some “exotic” physics, such as Lorentz invariance violation (LIV) or an axion origin of the signal. However, detailed discussions about these scenarios will be plagued with uncertainties of EBL models. Simple constraints on these new physics scenarios have been discussed.

Materials and Methods

KM2A is composed of 5216 electromagnetic particle detectors (EDs) and 1188 muon detectors (MDs), which are distributed in an area of $1.3\ \mathrm{km}^2$. The EDs are designed to detect the gamma-ray/cosmic ray incident showers with determining their directions and energies. The MDs are designed to detect the muon component of showers, which is used to discriminate between gamma-ray and hadron-induced showers. The whole KM2A detector was completed and operational on July 19th, 2021, and the duty cycle is about 99%. A trigger is generated when 20 EDs are fired within a 400 ns window, and the trigger rate is about 2.5 kHz. The performance, including angular resolution, energy resolution, and gamma-ray/cosmic-ray discrimination power, of KM2A for gamma-rays has been thoroughly tested using the observation of the Crab Nebula (16). The following will show some key information and methods for the observation of GRB 221009A with KM2A.

Effective area of KM2A for gamma-rays

The LHAASO observatory arrays can monitor the overhead sky with a zenith angle less than $50^\circ$. The zenith angle of GRB 221009A as a function of time is shown in panel A of Figure 8. The duty cycle of KM2A was 100% during this period. The zenith angles were $28.8^\circ$ at $T_0+230\ \mathrm{s}$, $31.2^\circ$ at $T_0+900\ \mathrm{s}$, and $35.1^\circ$ at $T_0+2000\ \mathrm{s}$. The effective area of KM2A for gamma-rays with different energies at these zenith angles is shown in the panel B of Figure 8. At a zenith angle of $31.2^\circ$, the effective area is $10{,}000\ \mathrm{m}^2$ at 4 TeV and $100{,}000\ \mathrm{m}^2$ at 7 TeV, reaching a roughly constant value of $900{,}000\ \mathrm{m}^2$ around 20 TeV. For comparison, the detector area of WCDA is $78{,}000\ \mathrm{m}^2$.

Energy reconstruction for SED measurement

In the normal KM2A data analysis pipeline, the particle density at 50 m (denoted as $\rho_{50}$) from the shower core location is used to evaluate the gamma-ray energy. For a certain zenith angle $\theta$ and a certain $\rho_{50}$ value, the probability distribution of the true energy can be achieved using the Bayes theorem:

where $f(E)$ is spectral function, $A_{\mathrm{eff}}(E,\theta)$ is the effective area at true energy $E$ and zenith angle $\theta$, $P(\rho_{50}\mid(E,\theta))$ is the probability to measure $\rho_{50}$ for a gamma-ray event with energy $E$ and zenith angle $\theta$. Obviously, the distribution of $P(E\mid(\rho_{50},\theta))$ depends on the assumption of $f(E)$. Usually, the median value of the distribution of $P(E\mid(\rho_{50},\theta))$ is used as the reconstructed energy $E_{\mathrm{rec}}$ for a gamma-ray event with a value of $\rho_{50}$ at zenith angle $\theta$ assuming a power-law spectral function with index of -3.0 (16). Taking into account the derived spectrum for GRB 221009A using the KM2A data, the relation function between $\rho_{50}$ and $E_{\mathrm{rec}}$ was renewed, assuming a power-law energy spectrum with an index of $\alpha=-4.7$. This updating will significantly reduces the value of $E_{\mathrm{rec}}$. In this work, events with reconstructed energy ($E_{\mathrm{rec}}$) above 2.5 TeV are divided into five bins per decade. The distributions of primary true energy ($E_{\mathrm{true}}$) of the simulated gamma-ray events in each reconstructed energy bin are shown in panel C of Figure 8, which roughly presents the energy resolution of KM2A at energies below 25 TeV.

Gamma-ray and cosmic ray background discrimination

Most of the events recorded by KM2A are cosmic ray-induced showers, which constitute the major background for gamma-ray observations. Considering that gamma-ray induced showers are muon-poor and cosmic ray-induced showers are muon-rich, the ratio $R=\log((N_{\mu}+0.0001)/N_{\mathrm{e}})$ between the measured muons and electrons is used to discriminate primary gamma-rays from cosmic nuclei. The panel D of Figure 8 shows the distribution of this ratio for all events from the GRB direction and the background regions. The background regions are 20 off-source regions with the same zenith angle as the GRB. The ratio distribution is about the same for the source region and background region at $R>-1.9$, which is dominated by the cosmic ray background. The source region is clearly higher than the background region at $R<-1.9$, which is mainly due to gamma-ray signals. The events with $R<-1.9$ are used to select gamma-ray-like events to analyze the emission from the GRB in this work. This selection criterion can remove 98% of the cosmic ray background. The survival fraction of gamma-rays is 74% according to MC simulation.

Spectrum and flux determination

The gamma-ray flux from GRB 221009A is estimated using the number of excess events and the corresponding statistical uncertainty in each energy bin. Combined with the WCDA measurement, the gamma-ray emission from the GRB is assumed to follow a log-parabola (LP) spectrum or a power-law with an exponential cutoff (PLEC) spectrum. The response of the KM2A detector was simulated by tracing the trajectory of the GRB within the FOV of KM2A. The expected energy distribution detected by detector $N(E_{\mathrm{rec}})$ can be achieved:

where $A_{\mathrm{eff}}(E,\theta)$ is effective area at true energy $E$ and zenith angle $\theta$, $P(E_{\mathrm{rec}}|(E,\theta))$ is the probability to measure $E_{\mathrm{rec}}$ for an event with energy $E$ and zenith angle $\theta$, $\theta(t)$ is the zenith angle of the GRB at time t. With this energy distribution, the number of signals expected by the MC simulation $N_{\mathrm{MC}i}(f(E))$ using the spectrum $f(E)$ in each energy bin can be achieved. Then, the best-fit values of the spectrum $f(E)$ can be obtained using a forward-folded method to minimize a $\chi^2$ function for all energy bins:

where $N_{S_i}$ is the number of excess events, and $\sigma_{N_{S_i}}$ is the uncertainty of $N_{S_i}$ in the ith energy bin. $n$ denotes the number of energy bins. In this forward-folded method, the biases and energy resolution in the energy assignments shown in the Figure 8 are taken into account. The resulting differential flux has been shown in the Figure 2. The detailed information about these results are listed in Table S1. The median energy $E_{\mathrm{LP}}$ and $E_{\mathrm{PLEC}}$ is the median value of the probability distribution of the true energy $P(E|(E_1<E_{\mathrm{rec}}<E_2))$ for a reconstructed energy bin $[E_1,E_2]$ assuming corresponding LP and PLEC spectrum, respectively.

Time after $T_{0}$	$\log(E_{\mathrm{rec}}/TeV)$	$N_{\mathrm{on}}$	$N_{b}$	$E_{LP}$ (TeV)	Flux (LP/$10^{-10}$) (ergs cm$^{-2}$ s$^{-1}$)	$E_{PLEC}$ (TeV)	Flux (PLEC/$10^{-10}$) (ergs cm$^{-2}$ s$^{-1}$)
230-300s	0.4–0.6	8	0.40	4.84	$549^{+240}_{-195}$	4.32	$548^{+239}_{-194}$
	0.6–0.8	24	1.20	6.10	$179^{+42}_{-37}$	5.31	$207^{+48}_{-42}$
	0.8–1.0	29	1.55	7.85	$111^{+23}_{-20}$	6.53	$156^{+32}_{-28}$
	1.0–1.2	3	0.15	11.6	$12.8^{+10.3}_{-8.2}$	8.61	$28.5^{+23.1}_{-18.4}$
	1.2–1.4	0	0	18.8	$<11$	13.0	$<46$
	1.4–1.6	0	0	29.2	$<15$	18.8	$<152$
300-900s	0.4–0.6	10	1.25	4.62	$103^{+45}_{-38}$	4.12	$100^{+44}_{-37}$
	0.6–0.8	34	6.05	5.82	$35.2^{+8.0}_{-7.2}$	4.95	$41.0^{+9.3}_{-8.3}$
	0.8–1.0	28	5.45	7.50	$12.9^{+3.3}_{-3.0}$	5.96	$20.4^{+5.3}_{-4.7}$
	1.0–1.2	5	1.40	10.8	$2.42^{+1.89}_{-1.48}$	8.04	$5.46^{+4.26}_{-3.35}$
	1.2–1.4	1	0.15	17.6	$0.59^{+1.22}_{-0.42}$	11.9	$3.01^{+6.22}_{-2.16}$
	1.4–1.6	0	0	27.2	$<2.2$	17.2	$<29$

The flux at corresponding median energy $E_m$ is achieved using $\frac{N_{S_i}}{N_{\mathrm{MC}_i}(f(E))}f(E_m)$.

Energy reconstruction for the highest energy events

Due to the poor energy resolution of KM2A at the low-energy band, the distribution of $P(E_{\mathrm{rec}\mid}(E,\theta))$ is wide across different $E_{\mathrm{rec}}$ for a given true energy $E$. Hence, the probability distribution of the true energy $E$ for a reconstructed energy $E_{\mathrm{rec}}$ is also much wide, which can be achieved using the Bayes theorem:

Obviously, the distribution of $P(E\mid(E_{\mathrm{rec}},\theta))$ depends on the assumption of the spectral function $f(E)$. To correctly estimate the true energies and corresponding energy errors of the highest energy events observed by KM2A, the $P(E\mid(E_{\mathrm{rec}},\theta))$ for each event using different $f(E)$ is achieved. The median energy and corresponding errors can be achieved by integrating $P(E\mid(E_{\mathrm{rec}},\theta))$ from 0 TeV to $E_{\xi}$ when the corresponding value equal to $\xi$:

The median energy is estimated using $\xi= 0.5$. The corresponding errors are estimated using $\xi= 0.16$ and $\xi= 0.84$, respectively.

C. Bigongiari, A. Biland, O. Blanch, G. Bonnoli, Ž. Bošnjak, G. Busetto, R. Carosi, G. Ceribella, Y. Chai, A. Chilingaryan, S. Cikota, S. M. Colak, U. Colin, E. Colombo, J. L. Contreras, J. Cortina, S. Covino, V. D’Elia, P. da Vela, F. Dazzi, A. de Angelis, B. de Lotto, M. Delfino, J. Delgado, D. Depaoli, F. di Pierro, L. di Venere, E. Do Souto Espiñeira, D. Dominis Prester, A. Donini, D. Dorner, M. Doro, D. Elsaesser, V. Fallah Ramazani, A. Fattorini, G. Ferrara, D. Fidalgo, L. Foffano, M. V. Fonseca, L. Font, C. Fruck, S. Fukami, R. J. García López, M. Garczarczyk, S. Gasparyan, M. Gaug, N. Giglietto, F. Giordano, N. Godinović, D. Green, D. Guberman, D. Hadasch, A. Hahn, J. Herrera, J. Hoang, D. Hrupec, M. Hütten, T. Inada, S. Inoue, K. Ishio, Y. Iwamura, L. Jouvin, D. Kerszberg, H. Kubo, J. Kushida, A. Lamastra, D. Lelas, F. Leone, E. Lindfors, S. Lombardi, F. Longo, M. López, R. López-Coto, A. López-Oramas, S. Loporchio, B. Machado de Oliveira Fraga, C. Maggio, P. Majumdar, M. Makariev, M. Mallamaci, G. Maneva, M. Manganaro, K. Mannheim, L. Maraschi, M. Mariotti, M. Martínez, D. Mazin, S. Mićanović, D. Miceli, M. Minev, J. M. Miranda, R. Mirzoyan, E. Molina, A. Moralejo, D. Morcuende, V. Moreno, E. Moretti, P. Munar-Adrover, V. Neustroev, C. Nigro, K. Nilsson, D. Ninci, K. Nishijima, K. Noda, L. Nogués, S. Nozaki, S. Paiano, M. Palatiello, D. Paneque, R. Paoletti, J. M. Paredes, P. Peñil, M. Peresano, M. Persic, P. G. Prada Moroni, E. Prandini, I. Puljak, W. Rhode, M. Ribó, J. Rico, C. Righi, A. Rugliancich, L. Saha, N. Sahakyan, T. Saito, S. Sakurai, K. Satalecka, K. Schmidt, T. Schweizer, J. Sitarek, I. Šnidarić, D. Sobczynska, A. Somero, A. Stamerra, D. Strom, M. Strzys, Y. Suda, T. Surić, M. Takahashi, F. Tavecchio, P. Temnikov, T. Terzić, M. Teshima, N. Torres-Albà, L. Tosti, V. Vagelli, J. van Scherpenberg, G. Vanzo, M. Vazquez Acosta, C. F. Vigorito, V. Vitale, I. Vovk, M. Will, D. Zarić, L. Nava, P. Veres, P. N. Bhat, M. S. Briggs, W. H. Cleveland, R. Hamburg, C. M. Hui, B. Mailyan, R. D. Preece, O. J. Roberts, A. von Kienlin, C. A. Wilson-Hodge, D. Kocevski, M. Arimoto, D. Tak, K. Asano, M. Axelsson, G. Barbiellini, E. Bissaldi, F. Fana Dirirsa, R. Gill, J. Granot, J. McEnery, N. Omodei, S. Razzaque, F. Piron, J. L. Racusin, D. J. Thompson, S. Campana, M. G. Bernardini, N. P. M. Kuin, M. H. Siegel, S. B. Cenko, P. O’Brien, M. Capalbi, A. Dai, M. de Pasquale, J. Gropp, N. Klingler, J. P. Osborne, M. Perri, R. L. C. Starling, G. Tagliaferri, A. Tohuvavohu, A. Ursi, M. Tavani, M. Cardillo, C. Casentini, G. Piano, Y. Evangelista, F. Verrecchia, C. Pittori, F. Lucarelli, A. Bulgarelli, N. Parmiggiani, G. E. Anderson, J. P. Anderson, G. Bernardi, J. Bolmer, M. D. Caballero-García, I. M. Carrasco, A. Castellón, N. Castro Segura, A. J. Castro-Tirado, S. V. Cherukuri, A. M. Cockeram, P. D’Avanzo, A. di Dato, R. Diretse, R. P. Fender, E. Fernández-García, J. P. U. Fynbo,

Acknowledgments

We would like to thank all staff members who work at the LHAASO site above 4400 meters above sea level year-round to maintain the detector and keep the water recycling system, electricity power supply and other components of the experiment operating smoothly. We are grateful to Chengdu Management Committee of Tianfu New Area for the constant financial support for research with LHAASO data. We deeply appreciate the computing and data service support provided by the National High Energy Physics Data Center for the data analysis in this paper. Funding: This research work is supported by the following grants: The National Key R&D program of China No.2018YFA0404201, No.2018YFA0404202, No.2018YFA0404203, No.2018YFA0404204, National Natural Science Foundation of China No.12022502, No.U1831208, No.12205314, No.12105301, No.12261160362, No.12105294, No.U1931201, No.12005246, No.12173039, No.12121003, No.12333006, Department of Science and Technology of Sichuan Province, China No.2021YFSY0030, Project for Young Scientists in Basic Research of Chinese Academy of Sciences No.YSBR-061, and in Thailand by the NSRF via the Program Management Unit for Human Resources & Institutional Development, Research and Innovation (No. B37G660015).

Author contributions: S.Z. Chen and X.J. Bi led the writing of the text about the data analysis and interpretation, respectively. S.Z. Chen performed the data analysis of KM2A and S. Wu provided the cross-check. S.C. Hu (supervised by M. Zha) performed the spectrum analysis of WCDA. X.J. Bi led the interpretation on EBL, LIV and axions. X.Y. Wang and J. H. Zheng performed modelling of the SED using a GRB afterglow model. Zhen Cao is the spokesperson of the LHAASO Collaboration and the principal investigator of the LHAASO project and coordinated the specific working group for this paper involving all corresponding authors. F. Aharonian provided crucial comments on drafting of the manuscript. All other authors participated in data analysis, including detector calibration, data processing, event reconstruction, data quality check, and various simulations, and provided comments on the manuscript.

Competing Interests: There is no conflict of interest of the collaboration members. All relevant funding grants are listed in the Acknowledgments section.

Data and materials availability: The data supporting the conclusions of this paper are available in the Supplementary Materials section. The data of figures can be found at https://cstr.cn/CSTR:17081.11.opendata.LHAASO.20231121112426

• Corresponding authors: S.Z. Chen([email protected]), X.J. Bi ([email protected]), S.C. Hu([email protected]), X.Y. Wang ([email protected])

This PDF file includes:

• LHAASO Collaboration Author List

• Supplementary Materials

• Figures S1 to S3

• Tables S1 to S3

LHAASO Collaboration authors and affiliations

Zhen Cao$^{1,2,3}$, F. Aharonian$^{4,5}$, Q. An$^{6,7}$, Axikegu$^{8}$, Y.X. Bai$^{1,3}$, Y.W. Bao$^{9}$, D. Bastieri$^{10}$, X.J. Bi$^{1,2,3}$, Y.J. Bi$^{1,3}$, J.T. Cai$^{10}$, Q. Cao$^{11}$, W.Y. Cao$^{7}$, Zhe Cao$^{6,7}$, J. Chang$^{12}$, J.F. Chang$^{1,3,6}$, A.M. Chen$^{13}$, E.S. Chen$^{1,2,3}$, Liang Chen$^{14}$, Lin Chen$^{8}$, Long Chen$^{8}$, M.J. Chen$^{1,3}$, M.L. Chen$^{1,3,6}$, Q.H. Chen$^{8}$, S.H. Chen$^{1,2,3}$, S.Z. Chen$^{1,3}$, T.L. Chen$^{15}$, Y. Chen$^{9}$, N. Cheng$^{1,3}$, Y.D. Cheng$^{1,3}$, M.Y. Cui$^{12}$, S.W. Cui$^{11}$, X.H. Cui$^{16}$, Y.D. Cui$^{17}$, B.Z. Dai$^{18}$, H.L. Dai$^{1,3,6}$, Z.G. Dai$^{7}$, Danzengluobu$^{15}$, D. della Volpe$^{19}$, X.Q. Dong$^{1,2,3}$, K.K. Duan$^{12}$, J.H. Fan$^{10}$, Y.Z. Fan$^{12}$, J. Fang$^{18}$, K. Fang$^{1,3}$, C.F. Feng$^{20}$, L. Feng$^{12}$, S.H. Feng$^{1,3}$, X.T. Feng$^{20}$, Y.L. Feng$^{15}$, S. Gabici$^{21}$, B. Gao$^{1,3}$, C.D. Gao$^{20}$, L.Q. Gao$^{1,2,3}$, Q. Gao$^{15}$, W. Gao$^{1,3}$, W.K. Gao$^{1,2,3}$, M.M. Ge$^{18}$, L.S. Geng$^{1,3}$, G. Giacinti$^{13}$, G.H. Gong$^{22}$, Q.B. Gou$^{1,3}$, M.H. Gu$^{1,3,6}$, F.L. Guo$^{14}$, X.L. Guo$^{8}$, Y.Q. Guo$^{1,3}$, Y.Y. Guo$^{12}$, Y.A. Han$^{23}$, H.H. He$^{1,2,3}$, H.N. He$^{12}$, J.Y. He$^{12}$, X.B. He$^{17}$, Y. He$^{8}$, M. Heller$^{19}$, Y.K. Hor$^{17}$, B.W. Hou$^{1,2,3}$, C. Hou$^{1,3}$, X. Hou$^{24}$, H.B. Hu$^{1,2,3}$, Q. Hu$^{7,12}$, S.C. Hu$^{1,2,3}$, D.H. Huang$^{8}$, T.Q. Huang$^{1,3}$, W.J. Huang$^{17}$, X.T. Huang$^{20}$, X.Y. Huang$^{12}$, Y. Huang$^{1,2,3}$, Z.C. Huang$^{8}$, X.L. Ji$^{1,3,6}$, H.Y. Jia$^{8}$, K. Jia$^{20}$, K. Jiang$^{6,7}$, X.W. Jiang$^{1,3}$, Z.J. Jiang$^{18}$, M. Jin$^{8}$, M.M. Kang$^{25}$, T. Ke$^{1,3}$, D. Kuleshov$^{26}$, K. Kurinov$^{26}$, B.B. Li$^{11}$, Cheng Li$^{6,7}$, Cong Li$^{1,3}$, D. Li$^{1,2,3}$, F. Li$^{1,3,6}$, H.B. Li$^{1,3}$, H.C. Li$^{1,3}$, H.Y. Li$^{7,12}$, J. Li$^{7,12}$, Jian Li$^{7}$, Jie Li$^{1,3,6}$, K. Li$^{1,3}$, W.L. Li$^{20}$, W.L. Li$^{13}$, X.R. Li$^{1,3}$, Xin Li$^{6,7}$, Y.Z. Li$^{1,2,3}$, Zhe Li$^{1,3}$, Zhuo Li$^{27}$, E.W. Liang$^{28}$, Y.F. Liang$^{28}$, S.J. Lin$^{17}$, B. Liu$^{7}$, C. Liu$^{1,3}$, D. Liu$^{20}$, H. Liu$^{8}$, H.D. Liu$^{23}$, J. Liu$^{1,3}$, J.L. Liu$^{1,3}$, J.Y. Liu$^{1,3}$, M.Y. Liu$^{15}$, R.Y. Liu$^{9}$, S.M. Liu$^{8}$, W. Liu$^{1,3}$, Y. Liu$^{10}$, Y.N. Liu$^{22}$, R. Lu$^{18}$, Q. Luo$^{17}$, H.K. Lv$^{1,3}$, B.Q. Ma$^{27}$, L.L. Ma$^{1,3}$, X.H. Ma$^{1,3}$, J.R. Mao$^{24}$, Z. Min$^{1,3}$, W. Mitthumsiri$^{29}$, H.J. Mu$^{23}$, Y.C. Nan$^{1,3}$, A. Neronov$^{21}$, Z.W. Ou$^{17}$, B.Y. Pang$^{8}$, P. Pattarakijwanich$^{29}$, Z.Y. Pei$^{10}$, M.Y. Qi$^{1,3}$, Y.Q. Qi$^{11}$, B.Q. Qiao$^{1,3}$, J.J. Qin$^{7}$, D. Ruffolo$^{29}$, A. Saíz$^{29}$, D. Semikoz$^{21}$, C.Y. Shao$^{17}$, L. Shao$^{11}$, O. Shchegolev$^{26,30}$, X.D. Sheng$^{1,3}$, F.W. Shu$^{31}$, H.C. Song$^{27}$, Yu.V. Stenkin$^{26,30}$, V. Stepanov$^{26}$, Y. Su$^{12}$, Q.N. Sun$^{8}$, X.N. Sun$^{28}$, Z.B. Sun$^{32}$, P.H.T. Tam$^{17}$, Q.W. Tang$^{31}$, Z.B. Tang$^{6,7}$, W.W. Tian$^{2,16}$, C. Wang$^{32}$, C.B. Wang$^{8}$, G.W. Wang$^{7}$, H.G. Wang$^{10}$, H.H. Wang$^{17}$, J.C. Wang$^{24}$, K. Wang$^{9}$, L.P. Wang$^{20}$, L.Y. Wang$^{1,3}$, P.H. Wang$^{8}$, R. Wang$^{20}$, W. Wang$^{17}$, X.G. Wang$^{28}$, X.Y. Wang$^{9}$, Y. Wang$^{8}$, Y.D. Wang$^{1,3}$, Y.J. Wang$^{1,3}$, Z.H. Wang$^{25}$, Z.X. Wang$^{18}$, Zhen Wang$^{13}$, Zheng Wang$^{1,3,6}$, D.M. Wei$^{12}$, J.J. Wei$^{12}$, Y.J. Wei$^{1,2,3}$, T. Wen$^{18}$, C.Y. Wu$^{1,3}$, H.R. Wu$^{1,3}$, S. Wu$^{1,3}$, X.F. Wu$^{12}$, Y.S. Wu$^{7}$, S.Q. Xi$^{1,3}$, J. Xia$^{7,12}$, J.J. Xia$^{8}$, G.M. Xiang$^{2,14}$, D.X. Xiao$^{11}$, G. Xiao$^{1,3}$, G.G. Xin$^{1,3}$, Y.L. Xin$^{8}$, Y. Xing$^{14}$, Z. Xiong$^{1,2,3}$, D.L. Xu$^{13}$, R.F. Xu$^{1,2,3}$, R.X. Xu$^{27}$, W.L. Xu$^{25}$, L. Xue$^{20}$, D.H. Yan$^{18}$, J.Z. Yan$^{12}$, T. Yan$^{1,3}$, C.W. Yang$^{25}$, F. Yang$^{11}$, F.F. Yang$^{1,3,6}$, H.W. Yang$^{17}$, J.Y. Yang$^{17}$, L.L. Yang$^{17}$, M.J. Yang$^{1,3}$, R.Z. Yang$^{7}$, S.B. Yang$^{18}$, Y.H. Yao$^{25}$, Z.G. Yao$^{1,3}$, Y.M. Ye$^{22}$, L.Q. Yin$^{1,3}$, N. Yin$^{20}$ X.H. You$^{1,3}$, Z.Y. You$^{1,3}$, Y.H. Yu$^{7}$, Q. Yuan$^{12}$, H. Yue$^{1,2,3}$, H.D. Zeng$^{12}$, T.X. Zeng$^{1,3,6}$, W. Zeng$^{18}$, M. Zha$^{1,3}$, B.B. Zhang$^{9}$, F. Zhang$^{8}$, H.M. Zhang$^{9}$, H.Y. Zhang$^{1,3}$, J.L. Zhang$^{16}$, L.X. Zhang$^{10}$, Li Zhang$^{18}$, P.F. Zhang$^{18}$, P.P. Zhang$^{7,12}$, R. Zhang$^{7,12}$, S.B. Zhang$^{2,16}$, S.R. Zhang$^{11}$, S.S. Zhang$^{1,3}$, X. Zhang$^{9}$, X.P. Zhang$^{1,3}$, Y.F. Zhang$^{8}$, Yi Zhang$^{1,12}$, Yong Zhang$^{1,3}$, B. Zhao$^{8}$, J. Zhao$^{1,3}$, L. Zhao$^{6,7}$, L.Z. Zhao$^{11}$, S.P. Zhao$^{12,20}$, F. Zheng$^{32}$, J.H. Zheng$^{9}$, B. Zhou$^{1,3}$, H. Zhou$^{13}$, J.N. Zhou$^{14}$, M. Zhou$^{31}$, P. Zhou$^{9}$, R. Zhou$^{25}$, X.X. Zhou$^{8}$, C.G. Zhu$^{20}$, F.R. Zhu$^{8}$, H. Zhu$^{16}$, K.J. Zhu$^{1,2,3,6}$, X. Zuo$^{1,3}$, (The LHAASO Collaboration)

$^{1}$ Key Laboratory of Particle Astrophysics & Experimental Physics Division & Computing Center, Institute of High Energy Physics, Chinese Academy of Sciences, 100049 Beijing, China

$^{2}$ University of Chinese Academy of Sciences, 100049 Beijing, China

$^{3}$ TIANFU Cosmic Ray Research Center, Chengdu, Sichuan, China

$^{4}$ Dublin Institute for Advanced Studies, 31 Fitzwilliam Place, 2 Dublin, Ireland

$^{5}$ Max-Planck-Institut for Nuclear Physics, P.O. Box 103980, 69029 Heidelberg, Germany

$^{6}$ State Key Laboratory of Particle Detection and Electronics, China

$^{7}$ University of Science and Technology of China, 230026 Hefei, Anhui, China

$^{8}$ School of Physical Science and Technology & School of Information Science and Technology, Southwest Jiaotong University, 610031 Chengdu, Sichuan, China

$^{9}$ School of Astronomy and Space Science, Nanjing University, 210023 Nanjing, Jiangsu, China

$^{10}$ Center for Astrophysics, Guangzhou University, 510006 Guangzhou, Guangdong, China

$^{11}$ Hebei Normal University, 050024 Shijiazhuang, Hebei, China

$^{12}$ Key Laboratory of Dark Matter and Space Astronomy & Key Laboratory of Radio Astronomy, Purple Mountain Observatory, Chinese Academy of Sciences, 210023 Nanjing, Jiangsu, China

$^{13}$ Tsung-Dao Lee Institute & School of Physics and Astronomy, Shanghai Jiao Tong University, 200240 Shanghai, China

$^{14}$ Key Laboratory for Research in Galaxies and Cosmology, Shanghai Astronomical Observatory, Chinese Academy of Sciences, 200030 Shanghai, China

$^{15}$ Key Laboratory of Cosmic Rays (Tibet University), Ministry of Education, 850000 Lhasa, Tibet, China

$^{16}$ National Astronomical Observatories, Chinese Academy of Sciences, 100101 Beijing, China

$^{17}$ School of Physics and Astronomy (Zhuhai) & School of Physics (Guangzhou) & Sino-French Institute of Nuclear Engineering and Technology (Zhuhai), Sun Yat-sen University, 519000 Zhuhai & 510275 Guangzhou, Guangdong, China {sup}18` School of Physics and Astronomy, Yunnan University, 650091 Kunming, Yunnan, China

{sup}19` Département de Physique Nucléaire et Corpusculaire, Faculté de Sciences, Université de Genève, 24 Quai Ernest Ansermet, 1211 Geneva, Switzerland

{sup}20` Institute of Frontier and Interdisciplinary Science, Shandong University, 266237 Qingdao, Shandong, China

{sup}21` APC, Université Paris Cité, CNRS/IN2P3, CEA/IRFU, Observatoire de Paris, 119 75205 Paris, France

{sup}22` Department of Engineering Physics, Tsinghua University, 100084 Beijing, China

{sup}23` School of Physics and Microelectronics, Zhengzhou University, 450001 Zhengzhou, Henan, China

{sup}24` Yunnan Observatories, Chinese Academy of Sciences, 650216 Kunming, Yunnan, China

{sup}25` College of Physics, Sichuan University, 610065 Chengdu, Sichuan, China

{sup}26` Institute for Nuclear Research of Russian Academy of Sciences, 117312 Moscow, Russia

{sup}27` School of Physics, Peking University, 100871 Beijing, China

{sup}28` School of Physical Science and Technology, Guangxi University, 530004 Nanning, Guangxi, China

{sup}29` Department of Physics, Faculty of Science, Mahidol University, Bangkok 10400, Thailand

{sup}30` Moscow Institute of Physics and Technology, 141700 Moscow, Russia

{sup}31` Center for Relativistic Astrophysics and High Energy Physics, School of Physics and Materials Science & Institute of Space Science and Technology, Nanchang University, 330031 Nanchang, Jiangxi, China

{sup}32` National Space Science Center, Chinese Academy of Sciences, 100190 Beijing, China

Supplementary Materials

LHAASO detector

LHAASO [13] [15] consists of three detector arrays: square kilometer array (KM2A), Water Cherenkov Detector Array (WCDA), and Wide-Field-of-view Cherenkov Telescope Array (WFCTA). WFCTA is mainly for cosmic ray physics, while the two particle detector arrays KM2A and WCDA are mainly for gamma-ray physics. When a high-energy extraterrestrial particle, a gamma-ray or cosmic ray, enters Earth’s atmosphere, it initiates a cascade consisting of secondary hadrons, muons, leptons, and gamma-rays known as an air shower. The WCDA and KM2A detectors record different components of these air showers, which are used to reconstruct the type, energy, and arrival direction of the primary particles.

WCDA consists of three water ponds with a total area of $300\ \mathrm{m} \times260\ \mathrm{m}$ and 3120 detector units. Each detector unit is $5\ \mathrm{m} \times5\ \mathrm{m}$ and is separated by non-reflecting black plastic curtains and equipped with two upward-facing PMTs on the bottom at the center of the unit. Each pond is filled with purified water up to 4 m above the photo-cathodes of the PMTs. The whole LHAASO-WCDA detector has been operational since March 5th, 2021, and the duty cycle is about 98%. A trigger algorithm was implemented to record air showers by requiring at least 30 PMTs fired among a $12 \times12$ PMT array simultaneously within a window of 250 ns, and the trigger rate is around 35 kHz. The event reconstruction method and the corresponding performance of the array is described elsewhere (55).

KM2A is composed of 5216 electromagnetic particle detectors (EDs) and 1188 muon detectors (MDs), which are distributed in an area of $1.3\ \mathrm{km}^{2}$. Each ED consists of a $1\ \mathrm{m}^{2}$ plastic scintillator covered by a 0.5 cm thick lead plate and equipped with a 1.5-inch photomultiplier tube (PMT). Each MD consists of a cylindrical water tank, with a diameter of 6.8 m and a height of 1.2 m, and an 8-inch PMT, which is buried under 2.5 m of soil. The MDs are designed to detect the muon component of showers, which is used to discriminate between gamma-ray and hadron-induced showers. The whole KM2A detector was completed and operational on July 19th, 2021, and the duty cycle is about 99%. A trigger is generated when 20 EDs are fired within a 400 ns window, and the trigger rate is about 2.5 kHz. The performance, including angular resolution, energy resolution, and gamma-ray/cosmic-ray discrimination power, of KM2A for gamma-rays has been thoroughly tested using the observation of the Crab Nebula (16).

The detailed spectral information observed by KM2A

The resulting differential flux has been shown in Figure 2. The detailed information about these results, including the number of events from the source region ($N_{\rm on}$), the number of background events ($N_{\rm b}$), and the differential flux at the median energy ($E_{\rm LP}$ for LP spectrum and $E_{\rm PLEC}$ for PLEC spectrum) of the bin, are listed in Table S1.

The highest energy events observed by KM2A

During the period from $T_0+230\ \mathrm{s}$ to $T_0+900\ \mathrm{s}$, nine events with reconstructed energy above 10 TeV were observed by KM2A, adopting a primary reconstruction as listed in Table S1. Detailed information about these events is listed in Table S2, including the number of detected secondary electromagnetic particles ($N_{\mathrm{e}}$), the number of detected muons ($N_{\mu}$), the incident zenith angle ($\theta$), the distance of the shower core from the nearest edge of the active detector array ($D_{\mathrm{edge}}$), the space angle ($\Delta\psi$) between the event and the direction of the GRB, the arrival time ($T_{\mathrm{event}}$) since $T_{0}$, median energy ($E_{\mathrm{LP}}$) and its errors (using the LP spectrum shown in panel A of Figure 2, median energy ($E_{\mathrm{PLEC}}$) and its errors using the PLEC spectrum shown in panel B of Figure 2, and median energy ($E_{\mathrm{EBL}}$) and its errors using the EBL model and spectrum shown in panel B of Figure 3. The chance probability of each event due to background is estimated using the characteristics of each event. Firstly, for each event, the number of background (denoted as $b$) and signal (denoted as $s$) events, with arriving time $T_{\mathrm{event}}$ during the period from $T_{0}+230\mathrm{s}$ to $T_{0}+900\mathrm{s}$, ratio less than $\log((N_{\mu}+0.0001)/N_{\mathrm{e}})$, space angle with the GRB less than ($\Delta\psi$), and reconstructed energy above $E_{\mathrm{rec}}$ of the event, is estimated. Hereafter, the chance probability of the event is calculated using $b/(s+b)$. The chance probability (denoted as P) of each event is also listed in Table S2.

SED fitting using different EBL models

The gamma-ray flux from GRB 221009A is estimated using the number of excess events and the corresponding statistical uncertainty in each energy bin. To test different EBL models, such as Saldana-Lopez et al. 2021 (17), Gilmore et al. 2012 (18), Dominguez et al. 2011 (19), and Finke et al. 2010 (20), we adopt a log-parabolic form to characterize the intrinsic GRB spectrum (corrected for EBL absorption) and then fit an attenuated model of the form

to the data. The intrinsic spectrum using different EBL models is shown in Figure S1. The $\chi^{2}/ndf$ values for different EBL models are listed in Table S3. For the interval from $T_{0}+230\mathrm{s}$ to $T_{0}+300\mathrm{s}$, the minimum $\chi^{2}/ndf$ is achieved using the Gilmore et al. 2012 model. However, the minimum $\chi^{2}/ndf$ is achieved using the Saldana-Lopez et al. 2021 model for the interval from $T_{0}+300\mathrm{s}$ to $T_{0}+900\mathrm{s}$. The total $\chi^{2}/ndf$ after summing up the two intervals is comparable in the two cases.

To further test the EBL model, we divided the distribution of EBL adopted in the Saldana-Lopez et al. 2021 (17) model into three wavelength ranges, i.e., $<8\ \mu\mathrm{m}$, 8 to 28 $\mu\mathrm{m}$, and $>28\ \mu\mathrm{m}$. We also adopted a log-parabolic form to characterize the intrinsic GRB spectrum and then fit an attenuated model to the spectral data. We then tuned the scaling factor of each EBL range to minimize the total $\chi^{2}/ndf$ with the sum of the two intervals. According to our fitting result, the best fit factors for the three EBL ranges are $1.30^{+0.33}_{-0.20}$, $1.20^{+0.23}_{-0.20}$, and $0.40^{+0.44}_{-0.16}$, respectively.

The LHAASO constrained EBL model is plotted in Figure 4. For comparation the Saldana-Lopez EBL model and data points from different EBL measurements are also shown in the figure. The intrinsic SEDs are presented in Figure 3. The corresponding $\chi^2/ndf$ values using the best fitting factors for the two intervals are also listed in Table S3.

Origin of the $\sim10$ TeV gamma-rays

In the work on the WCDA result of GRB 221009A (15), the multi-wavelength data, including Swift-XRT, Fermi-LAT, and WCDA data of GRB 221009A, are modeled with the synchrotron plus SSC radiation within the framework of the afterglow emission from external forward shocks. In this model, a GRB jet drives a forward shock expanding into the ambient medium, which accelerates electrons into a power law distribution described by $dN/d\gamma_{\mathrm{e}} \propto\gamma_{\mathrm{e}}^{-p}$, where $\gamma_{\mathrm{e}}$ is the electron Lorentz factor. The modeling takes into account a full Klein-Nishina cross section for the inverse Compton scattering and the internal $\gamma\gamma$ absorption within the shock region.

The comparison between the SSC models used previously for WCDA data (15) and at present for KM2A data is given in Figure 7. It can be seen that the SSC emission spectra become increasingly softer at higher energy, thus deviating from the data considerably at the highest energy. This is because both the Klein-Nishina effect and internal $\gamma\gamma$ absorption become stronger at higher energies.

To solve the discrepancy between the SSC model and observed data at the high-energy end, we can assume an additional hard spectral component that becomes dominant at high energies. One possible component is the hadronic emission from the shock-accelerated relativistic protons. It has been proposed that the proton synchrotron emission in the external reverse shock can produce a hadronic component (31). The peak energy from proton synchrotron emission can reach $\sim10$ TeV if the magnetic field equipartition factor is sufficiently high in the reverse shock (31). For a proton spectrum $dN_{\mathrm{p}}/dE \sim E^{-\alpha_{\mathrm{p}}}$, the spectral index of the energy flux ($\nu f_{\nu}$) from proton synchrotron emission is $(3-\alpha_{\mathrm{p}})/2$. With $\alpha_{\mathrm{p}}=2$, the spectral index of the energy flux is $1/2$, much harder than that of the SSC emission.

Another possibility for the hard spectral component above several TeV is an intergalactic electromagnetic cascade due to the propagation of ultra-high-energy cosmic rays (UHECRs) that are accelerated by internal or external shocks of GRB 221009A (32, 56, 57). If these UHECR protons can escape from the source and propagate through the extragalactic medium from their sources to Earth, the interactions lead to the production of secondary cascade particles. These particles can initiate various energy loss processes for the electromagnetic cascade, such as the inverse-Compton scattering of background photons to higher energies. Because some interactions occur so close to us that the generated TeV gamma-rays do not suffer from significant EBL absorption, the observed cascade radiation could have a hard spectrum (32). In this scenario, the extragalactic magnetic field (EGMF) can deflect the UHECRs and cause a time delay. To reconcile with the temporal property of the KM2A emission above a few TeV, which arrived within hundreds of seconds after the GRB trigger, a suitable value of the EGMF is required (32, 56). An alternative explanation could be to invoke a new hard leptonic component. This could be realized in the multi-zone models where the magnetic field is inhomogeneous throughout the emitting volume (33). Synchrotron photons from the strong magnetic field zone provide the dominant target for the IC cooling of the electrons in the weak magnetic field zone. If the IC cooling is in the Klein-Nishina regime, a hard electron distribution will be formed (33). A hard electron spectral component could also be formed by the hydrodynamical turbulence that is excited in the GRB forward shock and stochastically accelerates protons and electrons (34). The stochastic acceleration can yield a hard electron spectrum with $p < 2$, though the maximum electron energy depends on the model parameters.

Axion-like particle estimation

An axion-like particle (ALP) is a very light pseudoscalar boson with a characteristic coupling to two photons described by the Lagrangian:

where $a$ represents the ALP field, $F_{\mu\nu}$ and $\widetilde{F}^{\mu\nu}$ are the electromagnetic tensor and its dual, respectively, $\mathbf{E}$ and $\mathbf{B}$ are the electric and magnetic components, respectively, and $g_{a\gamma\gamma}$ is the coupling constant. According to Eq. (7), in the presence of external magnetic fields, ALP and gamma-ray conversion $\gamma\leftrightarrow a$ takes place. Therefore, when $\gamma$-rays propagate, they oscillate with axions. As a consequence, the absorption of high-energy $\gamma$-rays by EBL is weakened, and the optical depth for high-energy $\gamma$-rays decreases.

The propagation of $\gamma$-rays is described by a Schrödinger-like equation (58). The conversion occurs in three different environments, namely, in the source region, in extragalactic space, and in the Milky Way. As a benchmark scenario, we consider the minimal astrophysical environment. In the source region, we assume that the conversion occurs when $\gamma$-rays are emitted from the GRB and propagate in the host galaxy. We assume the transverse magnetic field component to be 0.5 $\mu\mathrm{G}$ with a coherence length of 10 kpc and an electron density of approximately 0.04 $\mathrm{cm}^{-3}$ (59). As the magnetic field in extragalactic space is very weak and has large uncertainties, we ignore the conversion in this region. The magnetic field of the Milky Way is characterized by the regular component of the model in (60). The electron density of the Milky Way is described by the NE2001 model (61).

Given the circumstances established above and adopting the Saldana-Lopez et al. EBL model (17), we solve the propagation equation and obtain the gamma-ray survival probability at Earth. The left panel of Figure S2 shows the gamma-ray survival probability for an axion mass $m_a = 10^{-7}\,\mathrm{eV}$ and for different coupling constants. It is shown that the heavy absorption of high energy $\gamma$-rays by EBL is greatly alleviated. However, if the coupling constant $g_{a\gamma\gamma}$ is too large, there will be too many high energy gamma-rays detected, and the fitting to data becomes worse, as shown in the right panel of Figure S2. This gives us a constraint on the axion coupling, as shown in Figure 6.

Lorentz Invariance Violation

Lorentz Invariance Violation (LIV) modifies the energy-momentum dispersion relation for photons. As a consequence, the threshold energy for photon-photon pair production is changed in the presence of LIV. For the subluminal case of first-order LIV, the energy-momentum conservation yields the modified pair-creation threshold (62):

where $E_{\mathrm{LIV}}$ represents the LIV energy scale. From Eq. (8), it is obvious that the threshold energy is increased, leading to the suppression of the pair production process and more transparency for high energy $\gamma$-rays.

The $\gamma$-ray survival probability in LIV is calculated in the Saldana-Lopez et al. EBL model for LIV energy scales ranging from 1 $M_{\mathrm{Pl}}$ to 2.8 $M_{\mathrm{Pl}}$, as shown in the left panel of Figure S3. Similar to the axion case, if $E_{\mathrm{LIV}}$ is too low, there will be too many high energy gamma-rays, leading to conflicts with observation, as shown in the right panel of Figure S3. This leads us to obtain a lower bound of the LIV energy scale at about 1.5 $M_{\mathrm{Pl}}$.