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

The LHAASO Collaboration111Corresponding authors: S.Z. Chen
(chensz@ihep.ac.cn), X.J. Bi (bixj@ihep.ac.cn),
S.C. Hu (hushicong@ihep.ac.cn), X.Y. Wang (xywang@nju.edu.cn)222The LHAASO
Collaboration authors and affiliations are listed in the supplementary
materials

> The highest energy photons from gamma-ray bursts (GRBs) have important
> implications for the physics of particle acceleration and constraints on the
> extragalactic background light (EBL). Here we report for the first time the
> detection of photons beyond 10 TeV from the unprecedentedly brightest GRB
> 221009A by the Large High Altitude Air-shower Observatory (LHAASO). The
> LHAASO-KM2A detector registered more than 140 photons with energies above 3
> TeV from GRB 221009A during 230$-$900s after the trigger. The intrinsic
> energy spectrum of gamma-rays can be described by a power-law after
> correcting for EBL absorption, with the highest energy photon of 13 TeV.
> Such a hard spectrum challenges the synchrotron self-Compton (SSC) scenario
> of relativistic electrons for the afterglow emission above several TeV,
> because the Klein-Nishina effect and the internal $\gamma\gamma$ absorption
> within the GRB generally lead to a spectral softening towards high energies.
> Observations of photons up to 13 TeV from a source with a measured redshift
> of z=0.151 requires more transparency in intergalactic space than previously
> expected, in order to avoid an unusual pile-up at the end of the spectrum.
> Alternatively, one may invoke new physics such as Lorentz Invariance
> Violation (LIV) or assume an axion origin of very high energy (VHE) signals.

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

##  1 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 TeV gamma-rays during the afterglow
period (?, ?, ?, ?). The synchrotron self-Compton (SSC) process of
relativistic electrons in the afterglow has been proposed to explain the
origin of the TeV emission (?, ?, ?). However, challenges (?) 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 (RA=288.264∘, Dec=19.768∘) was detected
by Fermi-GBM on October 9, 2022, at 13:16:59.99 UT (denoted as T0 hereafter)
(?), about one hour earlier than the Swift trigger (?). The GBM light curve
consists of two emission episodes: a single isolated peak from T0 to T0+20s,
followed by a longer, extremely bright, multi-pulsed emission episode from
about T0+220s to T0+550s (?). It was also detected by the Large Area Telescope
(Fermi-LAT) with the highest-energy photon of 99.3 GeV observed 240 seconds
after the GBM trigger (?). This GRB was also clearly observed by many other
detectors, such as Insight-HXMT, GECAM-C (?), and Konus-Wind (?). Follow-up
optical observations revealed that the redshift of this GRB is $z$=0.151
(corresponding to a distance $\sim$753 Mpc) (?). The estimated isotropic
energy release is about $10^{55}$ ergs (?, ?, ?), which denotes it as an
extremely energetic GRB.

##  2 LHAASO observations of GRB 221009A

LHAASO observed GRB 221009A at the highest energy. LHAASO (?) consists of
three interconnected detectors located at 4,410 m above sea level in Sichuan
Province, China. The sub-arrays, a 78,000 m2 Water Cherenkov Detector Array
(WCDA) and a 1.3 km2 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 T0-20000s, culminating at a zenith angle of 9.5∘ at T0-7000s, and left
the FOV at T0+6000s. The observations cover both the prompt period and the
afterglow period, implemented at TeV for a detectable GRB for the first time.
The zenith angles were 28∘ at T0 and 31.5∘ at T0+1000s, which are favorable
positions for LHAASO observation. Within 2000s after T0, more than 5000
photons at $>$0.5 TeV were detected by LHAASO-WCDA from GRB 221009A with
significance above 100 s.d. (?). A previous report has focused on the temporal
characteristics of the TeV emission, based on LHAASO-WCDA data at lower
energies, which revealed the earliest TeV afterglow from GRB 221009A (?).
Here, we focus on the spectral results extending to energies above 10 TeV, as
measured by LHAASO-KM2A.

LHAASO-KM2A is optimized for gamma-rays with energies from 10 TeV to 10 PeV
(?). It also has sensitivity to detect gamma-rays below 10 TeV with an
effective area of 10,000 m2 at 4 TeV and 100,000 m2 at 7 TeV, which overlaps
with LHAASO-WCDA. Below 10 TeV, KM2A can still remove 98% of the cosmic ray
background using a “muon-less” content criterion with an angular resolution of
about 1∘ and an energy resolution of around 40%. The count-rate light curve
observed by LHAASO-KM2A is shown in Figure 1. The emission started at T0+230s
with two flux peaks around T0+245s and T0+285s. The emission gradually
decreased since T0+300s and faded out after T0+900s. During the period from
T0+230s to T0+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 RA=288.26∘
$\pm$0.07∘, Dec=19.83${}^{\circ}\pm$0.07∘ using the KM2A data, which is
consistent with that of Fermi-LAT. The significance map can be found in the
Materials and methods section (see Figure S4).

Figure 1:  The photon-count light curve of GRB 221009A obtained by LHAASO-
KM2A. Each time-bin is 10s. The black curve indicates the events from the
angular cone centered on the GRB, and the blue curve indicates the number of
events due to cosmic ray background estimated from 20 similar angular cones at
off-source directions with the same zenith angle. The gray dashed lines
indicate the peak times of the multi-pulsed emission observed by GECAM-C (?)
in the MeV band. The green dashed lines indicate the times of T0+230s,
T0+300s, and T0+900s. The pink points indicate the energy marked by the right
label and the arrival time of each event.

The reconstruction of the energy for each event depends on the assumption of
the spectral function due to the limitation of energy resolution of KM2A at
the low-energy band. Assuming a log-parabola (denoted as LP) function (as
shown in panel A of Figure 2), eight events with energy above 10 TeV were
observed during the period from T0+230s to T0+900s. The event with the maximum
energy is at 17.8${}_{-5.1}^{+7.4}$ TeV, as reported in the GCN (?). However,
if assuming a power-law with exponential cutoff (denoted as PLEC) function (as
shown in panel B of Figure 2), the energy of the maximum energy event is
12.2${}_{-2.4}^{+3.5}$ TeV. The energies of each event shown in Figure 1 were
obtained under the PLEC assumption. If taking into account the absorption of
EBL, the energies will also have a model-dependent change, as will be
discussed later. The energies of the nine gamma-ray-like events with the
highest energies have been estimated under three assumptions, i.e., LP, PLEC,
and an EBL model. The chance probability due to the CR background has also
been estimated using the characteristics of each event (see the details in the
Materials and methods section). Gamma-ray events with energy above 100 TeV
were searched for during a long period from T0 to T0+6000s; however, no event
was detected.

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 T0+230s to T0+300s
and interval 2 from T0+300s to T0+900s. We combined the WCDA and KM2A data to
determine the spectral energy distribution (SED) of GRB 221009A, which is
shown in Figure 2, assuming an LP or 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. More details about the analysis of WCDA data can be found elsewhere (?).
The two measurements are consistent with each other in the overlapping
energies from 3 TeV to 7 TeV. Assuming an LP function, the two independent
measurements are dN/dE = (3.57$\pm$0.13)$\times 10^{-7}$
(E/1TeV)-(3.46±0.03)-(0.58±0.05)log(E/1TeV) TeV-1 cm-2 s-1 for interval 1, and
dN/dE = (6.37$\pm$0.16)$\times 10^{-8}$
(E/1TeV)-(3.38±0.02)-(0.73±0.05)log(E/1TeV) TeV-1 cm-2 s-1 for interval 2. The
corresponding $\chi^{2}/ndf$, where $ndf$ is the number of degrees of freedom,
are 14.1/9 and 37.3/10, respectively. If assuming an PLEC function, the two
independent measurements are dN/dE = (5.19$\pm$0.33)$\times 10^{-7}$
(E/1TeV)-2.87±0.05 e-E/(2.62±0.26TeV) TeV-1 cm-2 s-1 for interval 1, and dN/dE
= (1.05$\pm$0.06) $\times 10^{-7}$(E/1TeV)-2.65±0.05 e-E/(1.99±0.15TeV) TeV-1
cm-2 s-1 for interval 2. The corresponding $\chi^{2}/ndf$ are 10.1/9 and
18.3/10, respectively. The fitting using the PLEC function is better than that
using the LP function. The spectrum of interval 2 is slightly harder than that
of interval 1; however, the cut-off is slightly lower. 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.

Figure 2: Observed VHE spectra of GRB 221009A by LHAASO for the two intervals,
i.e., interval 1 from T0+230s to T0+300s (red points) and interval 2 from
T0+300s to T0+900s (blue points). The solid lines indicate the best-fitting
results, and the shaded regions indicate the 1-sigma error region. (A) The
log-parabola function is used to fit the observational data. (B) The power-law
with exponential cutoff function is adopted to fit the observational data.

##  3 EBL absorption of emission from GRB 221009A

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. 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 S5). Therefore, the spectra of GRB 221009A measured by LHAASO can be
adopted to test the 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 =
J0E${}^{a+b.log(E)}e^{-\tau(E)}$ to the data. The fitting results are given in
the Materials and methods section. The goodness of fit combining the two
intervals is comparable for different EBL models, as shown in Table S3. The
intrinsic GRB spectrum obtained using the recent EBL model of Saldana-Lopez et
al. 2021 (?) is shown in panel A of Figure 3. 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. Although the deviation is not statistically significant,
the measured flux is 11 times higher than the fit line, which is remarkable
(panel A of Figure 3). This situation is similar to that obtained for the
blazar Mrk501 during its huge flare in 1997 (?). Even though some hypotheses
were proposed to produce a non-classical sharp pile-up spectrum (?), 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 model 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.

Figure 3:  (A) Filled points show the intrinsic spectrum of GRB 221009A
corrected for EBL absorption using the model of Saldana-Lopez et al. 2021 (?).
The red points are for interval 1 from T0+230s to T0+300s, and the blue points
are for interval 2 from T0+300s to T0+900s. The solid lines indicate the best-
fitting results using the power-law function, and the shaded regions indicate
the 1-sigma error region. The unfilled points and shaded regions are
corresponding observed spectra. (B) Filled points show the intrinsic spectrum
of GRB 221009A corrected for EBL absorption using the LHAASO-constrained EBL
model. The red solid line indicates the best-fitting result for interval 1,
which is a power-law function, and the blue solid line indicates the best-
fitting result for interval 2, which is a log-parabolic function. The points
and shaded regions are similar to those in panel A.

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\pm 1/2)\lambda$
centered on $\lambda^{\approx}1.5(E/1TeV)\mu m$ (?). Therefore, the measured
gamma-ray spectra from GRB 221009A can be used to constrain the EBL
intensities at different wavelengths. The EBL model of Saldana-Lopez et al.
2021 (?) is independently re-scaled at three wavelength ranges, i.e.,
$\lambda<8\mu m$, 8$\mu m<\lambda<28\mu m$, and $\lambda>28\mu 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.20}^{+0.33}$,
1.20${}_{-0.20}^{+0.23}$, and 0.40${}_{-0.16}^{+0.44}$, respectively. The
constrained EBL intensity is shown in Figure S6. It should be emphasized that
the EBL intensity at large wavelengths is now close to the lower bound set by
galaxy count observations (?). The intrinsic spectra using the re-scaled best-
fit EBL model are shown in panel B of Figure 3. For interval 1, the intrinsic
SED follows a single power-law model with an index of -(2.12$\pm$0.03), with a
$\chi^{2}$/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$\pm$0.02)-(0.15$\pm$0.06).log(E/TeV). The corresponding $\chi^{2}$/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 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${}_{-2.4}^{+3.2}$ TeV.

##  4 Implication for GRB physics

The light curve of GRB 221009A in 0.2-5 TeV, as observed by LHAASO-WCDA, is
characterized by a rapid rise to a peak, followed by a decay that persists for
at least 3000 s after the peak (?). The smooth temporal profile suggests an
external shock origin for the emission, caused by the interaction of the GRB
ejecta with the ambient medium. The light curve observed by LHAASO-KM2A is
overall similar to that of WCDA, indicating that the bulk of KM2A emission
should have the same origin. 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.

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 (?, ?,
?). 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 (?). This can also be seen from the
comparison between the model that explains the WCDA data (?) and the KM2A data
(see Figure S7). Such a discrepancy would suggest more complicated processes
occurring during the early afterglow phase.

One possible solution to this discrepancy is to assume that an additional hard
spectral component emerges at the highest energy end. 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
(?). 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 (?). For a
proton spectrum $dN_{p}/dE\sim E^{-2}$, the photon 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 (?) can also produce a hard spectral component at the highest energy,
because $>10$ TeV photons could be produced by UHECRs that have propagated to
a nearby distance to us. To explain the observed flux of $>10$ TeV photons in
the hadronic scenario, both models require efficient production of UHECRs from
GRBs (?, ?).

##  5 Discussion on the highest energy photon

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 require 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 photon 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 (?), the expected number of events
beyond the highest energy during the interval from T0+300s to T0+900s is
0.009. 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$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, further discussions about these
scenarios will be plagued with uncertainties of EBL models. Therefore, instead
of adopting these “exotic” physics models to interpret the observation of the
highest energy photon, we discuss constraints on new physics using the spectra
of the GRB measured by LHAASO.

##  6 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 (?) here.

Figure 4:  Constraints on the axion-photon coupling constant by LHAASO
observation of the GRB 221009A. It improves the constraint from CAST (?) and
is comparable with that derived from observations by HESS (?) and of Mrk 421
(?).

The axion is a hypothetical particle that solves the strong CP problem of the
Standard Model by a spontaneously broken Peccei-Quinn symmetry. Axion-like
particles (ALPs) exist extensively in theoretical models developed beyond the
SM. Both axions and ALPs are among the most attractive cold dark matter
candidates. The most distinctive feature of the axion is its coupling vertex
with two photons, leading to axion-photon conversion in external magnetic
fields for astrophysical detection. The axion-photon conversion leads to a
suppression of EBL absorption of high-energy photons and higher transparency
of space. We still utilize a log-parabolic form to describe the intrinsic GRB
spectrum in this scenario. The goodness of the best fitting for the spectra
during the two intervals is improved by only about $1.4\sigma$ compared to
that without axion-photon 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 4,
the constraint on the axion-photon 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 (?) (CERN Axion Solar
Telescope) and is comparable to that derived from the observation of PKS
2155-304 energy spectrum by HESS (?) and a long-term monitoring of Mrk 421 by
ARGO-YBJ (?).

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 MPl = $1.22\times
10^{28}$eV. At low energies, LIV terms are also induced, which are suppressed
by MPl and lead to modification of the particles’ dispersion relation. This
may allow many interesting phenomena which are otherwise forbidden processes.
For example, an ultra-high energy free gamma-ray particle can decay into an
electron-positron pair in the vacuum automatically, which has been adopted by
LHAASO to explore the LIV effect with the observation of PeV photons (?).
Another consequence of LIV 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\rightarrow e^{-}e^{+}$ is
suppressed, and the universe becomes more transparent. In such a scenario, the
detection of a gamma-ray beyond 10 TeV event fits the expectation very well.
However, the significance of the LIV effect is only about $1.6\sigma$, 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.5 MPl.
This constraint is comparable to that derived by exploring the energy
dependence of the propagation speed of photons from GRB090510 by Fermi-LAT
(?).

In summary, we report the first detection of photons beyond 10 TeV from a
gamma-ray burst, from GRB 221009A during the early afterglow period by LHAASO.
The measured energy spectrum has a clear cutoff at the high-energy end, which
is broadly consistent with the standard picture of gamma-ray propagating and
interaction with EBL. However, to avoid a nonstandard shape of the intrinsic
spectrum obtained when correcting the observed spectrum for the absorption in
EBL, the EBL flux at mid-IR wavelengths should be smaller than current models’
prediction. 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 variation of
the spectrum during different periods suggests a complex process during the
early afterglow phase. These observations also place strong constraints on new
physics parameters, such as an axion origin of the signal or LIV.

## References

## References

  * 1. MAGIC Collaboration, et al., Nature 575, 455 (2019).
  * 2. MAGIC Collaboration, et al., Nature 575, 459 (2019).
  * 3. H. Abdalla, et al., Nature 575, 464 (2019).
  * 4. H. E. S. S. Collaboration, et al., Science 372, 1081 (2021).
  * 5. E. Derishev, T. Piran, Astrophys. J. 880, L27 (2019).
  * 6. X.-Y. Wang, R.-Y. Liu, H.-M. Zhang, S.-Q. Xi, B. Zhang, Astrophys. J. 884, 117 (2019).
  * 7. S. Lesage, et al., arXiv e-prints p. arXiv:2303.14172 (2023).
  * 8. J. A. Kennea, M. Williams, Swift Team, GRB Coordinates Network 32635, 1 (2022).
  * 9. R. Pillera, et al., GRB Coordinates Network 32658, 1 (2022).
  * 10. Z.-H. An, et al., arXiv e-prints p. arXiv:2303.01203 (2023).
  * 11. D. Frederiks, et al., arXiv e-prints p. arXiv:2302.13383 (2023).
  * 12. A. de Ugarte Postigo, et al., GRB Coordinates Network 32648, 1 (2022).
  * 13. X.-H. Ma, et al., Chinese Physics C 46, 030001 (2022).
  * 14. Y. Huang, et al., GRB Coordinates Network 32677, 1 (2022).
  * 15. LHAASO Collaboration, et al., Science 380, 1390 (2023).
  * 16. F. Aharonian, et al., Chinese Physics C 45, 025002 (2021).
  * 17. A. Saldana-Lopez, et al., Mon. Not. R. Astron. Soc. 507, 5144 (2021).
  * 18. F. A. Aharonian, et al., Astron. Astrophys. 349, 11 (1999).
  * 19. F. A. Aharonian, A. N. Timokhin, A. V. Plyasheshnikov, Astron. Astrophys. 384, 834 (2002).
  * 20. F. A. Aharonian, 27th International Cosmic Ray Conference (ICRC27) (2001), vol. 27 of International Cosmic Ray Conference, p. 250.
  * 21. S. Berta, et al., Astron. Astrophys. 518, L30 (2010).
  * 22. E. Nakar, S. Ando, R. Sari, Astrophys. J. 703, 675 (2009).
  * 23. B. T. Zhang, et al., Astrophys. J. 947, L14 (2023).
  * 24. S. Das, S. Razzaque, Astron. Astrophys. 670, L12 (2023).
  * 25. V. Anastassopoulos, et al., Nature Phys. 13, 584 (2017).
  * 26. A. Abramowski, et al., Phys. Rev. D 88, 102003 (2013).
  * 27. H.-J. Li, J.-G. Guo, X.-J. Bi, S.-J. Lin, P.-F. Yin, Phys. Rev. D 103, 083003 (2021).
  * 28. Z. Cao, et al., Phys. Rev. Lett. 128, 051102 (2022).
  * 29. M. Ackermann, et al., Nature 462, 331 (2009).
  * 30. F. Aharonian, et al., Chinese Physics C 45, 085002 (2021).
  * 31. R. C. Gilmore, R. S. Somerville, J. R. Primack, A. Domínguez, Mon. Not. R. Astron. Soc. 422, 3189 (2012).
  * 32. A. Domínguez, et al., Mon. Not. R. Astron. Soc. 410, 2556 (2011).
  * 33. J. D. Finke, S. Razzaque, C. D. Dermer, Astrophys. J. 712, 238 (2010).
  * 34. M. G. Hauser, et al., Astrophys. J. 508, 25 (1998).
  * 35. G. Lagache, L. M. Haffner, R. J. Reynolds, S. L. Tufte, Astron. Astrophys. 354, 247 (2000).
  * 36. J. P. Gardner, et al., Astron. J. 119, 486 (2000).
  * 37. D. Elbaz, et al., Astron. Astrophys. 384, 848 (2002).
  * 38. G. G. Fazio, et al., Astrophys. J. Suppl. 154, 39 (2004).
  * 39. C. K. Xu, et al., Astrophys. J. Lett. 619, L11 (2005).
  * 40. M. Bethermin, H. Dole, A. Beelen, H. Aussel, Astron. Astrophys. 512, A78 (2010).
  * 41. S. Matsuura, et al., Astrophys. J. 737, 2 (2011).
  * 42. E. N. Voyer, J. P. Gardner, H. I. Teplitz, B. D. Siana, D. F. de Mello, Astrophys. J. 736, 80 (2011).
  * 43. M. Zemcov, et al., Science 346, 732 (2014).
  * 44. S. P. Driver, et al., Astrophys. J. 827, 108 (2016).
  * 45. D. P. Finkbeiner, M. Davis, D. J. Schlegel, Astrophys. J. 544, 81 (2000).
  * 46. P. Madau, L. Pozzetti, Mon. Not. Roy. Astron. Soc. 312, L9 (2000).
  * 47. L. Metcalfe, et al., Astron. Astrophys. 407, 791 (2003).
  * 48. C. Papovich, et al., Astrophys. J. Suppl. 154, 70 (2004).
  * 49. D. T. Frayer, et al., Astrophys. J. Lett. 647, L9 (2006).
  * 50. R. C. Keenan, A. J. Barger, L. L. Cowie, W.-H. Wang, Astrophys. J. 723, 40 (2010).
  * 51. Y. Matsuoka, N. Ienaka, K. Kawara, S. Oyabu, Astrophys. J. 736, 119 (2011).
  * 52. A. Penin, et al., Astron. Astrophys. 543, A123 (2012).
  * 53. K. Mattila, P. Väisänen, K. Lehtinen, G. von Appen-Schnur, C. Leinert, Monthly Notices of the Royal Astronomical Society 470, 2152 (2017).
  * 54. R. Alves Batista, arXiv e-prints p. arXiv:2210.12855 (2022).
  * 55. A. Rudolph, M. Petropoulou, W. Winter, Ž. Bošnjak, Astrophys. J. 944, L34 (2023).
  * 56. G. Raffelt, L. Stodolsky, Phys. Rev. D 37, 1237 (1988).
  * 57. A. Fletcher, ASP Conf. Ser. 438, 197 (2011).
  * 58. R. Jansson, G. R. Farrar, Astrophys. J. 757, 14 (2012).
  * 59. J. M. Cordes, T. J. W. Lazio (2002).
  * 60. A. Saldana-Lopez, et al., Mon. Not. Roy. Astron. Soc. 507, 5144 (2021).
  * 61. J. Biteau, D. A. Williams, Astrophys. J. 812, 60 (2015).

## 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, 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 through the LHAASO web page
(http://english.ihep.cas.cn/lhaaso/index.html) in the section ‘Public Data ?
All data are in ASCII, and the code used to produce the figures is publicly
accessible and listed in the Public Data section. The data are also available
on the public data web page at https://www.nhepsdc.cn/resource/astro/lhaaso/.

##

Supplementary Materials for

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

The LHAASO Collaboration∗

∗ Corresponding authors: S.Z. Chen(chensz@ihep.ac.cn), X.J. Bi
(bixj@ihep.ac.cn), S.C. Hu(hushicong@ihep.ac.cn), X.Y. Wang
<EMAIL_ADDRESS>

This PDF file includes:

LHAASO Collaboration Author List
Materials and Methods
Figures S1 to S9
Tables S1 to S3

## LHAASO Collaboration authors and affiliations

Zhen Cao1,2,3, F. Aharonian4,5, Q. An6,7, Axikegu8, Y.X. Bai1,3, Y.W. Bao9, D.
Bastieri10, X.J. Bi1,2,3, Y.J. Bi1,3, J.T. Cai10, Q. Cao11, W.Y. Cao7, Zhe
Cao6,7, J. Chang12, J.F. Chang1,3,6, A.M. Chen13, E.S. Chen1,2,3, Liang
Chen14, Lin Chen8, Long Chen8, M.J. Chen1,3, M.L. Chen1,3,6, Q.H. Chen8, S.H.
Chen1,2,3, S.Z. Chen1,3, T.L. Chen15, Y. Chen9, N. Cheng1,3, Y.D. Cheng1,3,
M.Y. Cui12, S.W. Cui11, X.H. Cui16, Y.D. Cui17, B.Z. Dai18, H.L. Dai1,3,6,
Z.G. Dai7, Danzengluobu15, D. della Volpe19, X.Q. Dong1,2,3, K.K. Duan12, J.H.
Fan10, Y.Z. Fan12, J. Fang18, K. Fang1,3, C.F. Feng20, L. Feng12, S.H.
Feng1,3, X.T. Feng20, Y.L. Feng15, S. Gabici21, B. Gao1,3, C.D. Gao20, L.Q.
Gao1,2,3, Q. Gao15, W. Gao1,3, W.K. Gao1,2,3, M.M. Ge18, L.S. Geng1,3, G.
Giacinti13, G.H. Gong22, Q.B. Gou1,3, M.H. Gu1,3,6, F.L. Guo14, X.L. Guo8,
Y.Q. Guo1,3, Y.Y. Guo12, Y.A. Han23, H.H. He1,2,3, H.N. He12, J.Y. He12, X.B.
He17, Y. He8, M. Heller19, Y.K. Hor17, B.W. Hou1,2,3, C. Hou1,3, X. Hou24,
H.B. Hu1,2,3, Q. Hu7,12, S.C. Hu1,2,3, D.H. Huang8, T.Q. Huang1,3, W.J.
Huang17, X.T. Huang20, X.Y. Huang12, Y. Huang1,2,3, Z.C. Huang8, X.L. Ji1,3,6,
H.Y. Jia8, K. Jia20, K. Jiang6,7, X.W. Jiang1,3, Z.J. Jiang18, M. Jin8, M.M.
Kang25, T. Ke1,3, D. Kuleshov26, K. Kurinov26, B.B. Li11, Cheng Li6,7, Cong
Li1,3, D. Li1,2,3, F. Li1,3,6, H.B. Li1,3, H.C. Li1,3, H.Y. Li7,12, J. Li7,12,
Jian Li7, Jie Li1,3,6, K. Li1,3, W.L. Li20, W.L. Li13, X.R. Li1,3, Xin Li6,7,
Y.Z. Li1,2,3, Zhe Li1,3, Zhuo Li27, E.W. Liang28, Y.F. Liang28, S.J. Lin17, B.
Liu7, C. Liu1,3, D. Liu20, H. Liu8, H.D. Liu23, J. Liu1,3, J.L. Liu1,3, J.Y.
Liu1,3, M.Y. Liu15, R.Y. Liu9, S.M. Liu8, W. Liu1,3, Y. Liu10, Y.N. Liu22, R.
Lu18, Q. Luo17, H.K. Lv1,3, B.Q. Ma27, L.L. Ma1,3, X.H. Ma1,3, J.R. Mao24, Z.
Min1,3, W. Mitthumsiri29, H.J. Mu23, Y.C. Nan1,3, A. Neronov21, Z.W. Ou17,
B.Y. Pang8, P. Pattarakijwanich29, Z.Y. Pei10, M.Y. Qi1,3, Y.Q. Qi11, B.Q.
Qiao1,3, J.J. Qin7, D. Ruffolo29, A. Sáiz29, D. Semikoz21, C.Y. Shao17, L.
Shao11, O. Shchegolev26,30, X.D. Sheng1,3, F.W. Shu31, H.C. Song27, Yu.V.
Stenkin26,30, V. Stepanov26, Y. Su12, Q.N. Sun8, X.N. Sun28, Z.B. Sun32,
P.H.T. Tam17, Q.W. Tang31, Z.B. Tang6,7, W.W. Tian2,16, C. Wang32, C.B. Wang8,
G.W. Wang7, H.G. Wang10, H.H. Wang17, J.C. Wang24, K. Wang9, L.P. Wang20, L.Y.
Wang1,3, P.H. Wang8, R. Wang20, W. Wang17, X.G. Wang28, X.Y. Wang9, Y. Wang8,
Y.D. Wang1,3, Y.J. Wang1,3, Z.H. Wang25, Z.X. Wang18, Zhen Wang13, Zheng
Wang1,3,6, D.M. Wei12, J.J. Wei12, Y.J. Wei1,2,3, T. Wen18, C.Y. Wu1,3, H.R.
Wu1,3, S. Wu1,3, X.F. Wu12, Y.S. Wu7, S.Q. Xi1,3, J. Xia7,12, J.J. Xia8, G.M.
Xiang2,14, D.X. Xiao11, G. Xiao1,3, G.G. Xin1,3, Y.L. Xin8, Y. Xing14, Z.
Xiong1,2,3, D.L. Xu13, R.F. Xu1,2,3, R.X. Xu27, W.L. Xu25, L. Xue20, D.H.
Yan18, J.Z. Yan12, T. Yan1,3, C.W. Yang25, F. Yang11, F.F. Yang1,3,6, H.W.
Yang17, J.Y. Yang17, L.L. Yang17, M.J. Yang1,3, R.Z. Yang7, S.B. Yang18, Y.H.
Yao25, Z.G. Yao1,3, Y.M. Ye22, L.Q. Yin1,3, N. Yin20, X.H. You1,3, Z.Y.
You1,3, Y.H. Yu7, Q. Yuan12, H. Yue1,2,3, H.D. Zeng12, T.X. Zeng1,3,6, W.
Zeng18, M. Zha1,3, B.B. Zhang9, F. Zhang8, H.M. Zhang9, H.Y. Zhang1,3, J.L.
Zhang16, L.X. Zhang10, Li Zhang18, P.F. Zhang18, P.P. Zhang7,12, R. Zhang7,12,
S.B. Zhang2,16, S.R. Zhang11, S.S. Zhang1,3, X. Zhang9, X.P. Zhang1,3, Y.F.
Zhang8, Yi Zhang1,12, Yong Zhang1,3, B. Zhao8, J. Zhao1,3, L. Zhao6,7, L.Z.
Zhao11, S.P. Zhao12,20, F. Zheng32, J.H. Zheng9, B. Zhou1,3, H. Zhou13, J.N.
Zhou14, M. Zhou31, P. Zhou9, R. Zhou25, X.X. Zhou8, C.G. Zhu20, F.R. Zhu8, H.
Zhu16, K.J. Zhu1,2,3,6, X. Zuo1,3, (The LHAASO Collaboration)
1 Key Laboratory of Particle Astrophyics & 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
18 School of Physics and Astronomy, Yunnan University, 650091 Kunming, Yunnan,
China
19 Département de Physique Nucléaire et Corpusculaire, Faculté de Sciences,
Université de Genève, 24 Quai Ernest Ansermet, 1211 Geneva, Switzerland
20 Institute of Frontier and Interdisciplinary Science, Shandong University,
266237 Qingdao, Shandong, China
21 APC, Université Paris Cité, CNRS/IN2P3, CEA/IRFU, Observatoire de Paris,
119 75205 Paris, France
22 Department of Engineering Physics, Tsinghua University, 100084 Beijing,
China
23 School of Physics and Microelectronics, Zhengzhou University, 450001
Zhengzhou, Henan, China
24 Yunnan Observatories, Chinese Academy of Sciences, 650216 Kunming, Yunnan,
China
25 College of Physics, Sichuan University, 610065 Chengdu, Sichuan, China
26 Institute for Nuclear Research of Russian Academy of Sciences, 117312
Moscow, Russia
27 School of Physics, Peking University, 100871 Beijing, China
28 School of Physical Science and Technology, Guangxi University, 530004
Nanning, Guangxi, China
29 Department of Physics, Faculty of Science, Mahidol University, Bangkok
10400, Thailand
30 Moscow Institute of Physics and Technology, 141700 Moscow, Russia
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
32 National Space Science Center, Chinese Academy of Sciences, 100190 Beijing,
China

##  S1 Materials and methods

###  S1.1 LHAASO detector

LHAASO (?) consists of three detector arrays: KM2A, WCDA, and 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 photons 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 is formed with three water ponds with a total area of 300 m $\times$ 260
m and 3120 detector units. Each detector unit is 5 m $\times$ 5 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 $\times$ 12 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 can be found elsewhere (?).

KM2A is composed of 5216 electromagnetic particle detectors (EDs) and 1188
muon detectors (MDs), which are distributed in an area of 1.3 km2. Each ED
consists of a 1 m2 plastic scintillaton covered by a 0.5 cm thick lead plate
and 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 (?).

###  S1.2 LHAASO-KM2A observations of GRB 221009A

Both WCDA and KM2A arrays can monitor the overhead sky with a zenith angle
less than 50∘. For the LHAASO observatory, the zenith angle of GRB 221009A as
a function of time is shown in Figure S1. It entered into the FOV of LHAASO at
about T0-20000s, culminating at a zenith angle of 9.5∘ at T0-7000s, and left
the FOV at T0+6000s. The duty cycle of KM2A and WCDA were both 100% during
this period. The zenith angles were 28.8∘ at T0+230s, 31.2∘ at T0+900s, and
35.1∘ at T0+2000s. The effective area of LHAASO-KM2A for gamma-rays with
different energies at these zenith angles is shown in Figure S1. At a zenith
angle of 31.2∘, the effective area is 10,000 m2 at 4 TeV and 100,000 m2 at 7
TeV, reaching a roughly constant value of 900,000 m2 around 20 TeV. For
comparison, the detector area of WCDA is 78,000 m2.

Figure S1:  Panel A shows the zenith angle of GRB 221009A within the FOV of
LHAASO as a function of the observation time. The dotted green lines mark the
times of T0+230s and T0+900s, respectively. Panel B shows the effective area
of KM2A for gamma-rays as a function of the gamma-ray energy at three zenith
angles, 28.8∘, 31.2∘, and 35.1∘, corresponding to the observation times
T0+230s, T0+900s, T0+2000s, respectively.

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. Taking into account the derived spectrum for GRB 221009A
using the KM2A data, the relation function between $\rho_{50}$ and energy
(denoted as Erec) was renewed, assuming a power-law energy spectrum with an
index of $\alpha$=-4.7. In this work, events with reconstructed energy (Erec)
above 2.5 TeV are divided into five bins per decade. The distributions of
primary true energy (Etrue) of the simulated gamma-ray events in each
reconstructed energy bin are shown in Figure S2, which presents the energy
resolution of KM2A at energies below 25 TeV.

Figure S2:  The distributions of simulated gamma-ray primary true energy
(Etrue) in each reconstructed energy (Erec) bin in the energy range from 2.5
TeV to 25 TeV. The number of events has been rescaled to make the peak value
1. This figure presents the energy resolution of KM2A at energies below 25
TeV. The bumpy structure shown in the low energy tail of the lowest energy bin
is mainly due to statistical fluctuation.

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=$\mathrm{\log((N_{\mu}+0.0001)/N_{e})}$ between the
measured muons and electrons is used to discriminate primary gamma-rays from
cosmic nuclei. Figure S3 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.

Figure S3:  The distribution of the ratio
R=$\mathrm{\log((N_{\mu}+0.0001)/N_{e})}$, where Nu is the number of muons
measured by MDs and Ne is the number of electromagnetic particles measured by
EDs. The red curve represents the events from the GRB source region, while the
blue curve represents the events from the 20 off-source regions, with a weight
of 1/20 used for each event. The dotted green line marks the position of
R=-1.9.

The significance sky maps around GRB 221009A during the period from T0+230s to
T0+900s are shown in Figure S4, which are smoothed using the point-spread
function (PSF) of KM2A in the corresponding energy range. The significance is
estimated using a likelihood ratio test method, and details can be found in
our previous reference. Clear signals from GRB 221009A are detected with a
statistical significance of 20.6 $\sigma$. To estimate the position of the
gamma-ray signal around the GRB direction, a 2-dimensional Gaussian is used to
fit the event excess map. The positions obtained, in right ascension (RA) and
declination (Dec), are RA=288.26∘ $\pm$0.07∘ and Dec=19.83${}^{\circ}\pm$0.07∘
in J2000.0 epoch, which is consistent with that of Fermi-LAT within
statistical errors.

Figure S4:  The significance map around GRB 221009A as observed by KM2A. The
plus sign and corresponding length denote the position and error determined by
LHAASO-KM2A. The black circle denotes the position of the GRB reported by
Fermi-LAT. The white circle shows the size of the PSF that contains 68% of the
events.

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
best-fit values of the spectrum are obtained using a forward-folded method to
minimize a $\chi^{2}$ function for all energy bins:

$\chi^{2}=\sum_{i=1}^{n}\begin{pmatrix}\frac{N_{s_{i}}-N_{MC_{i}}(f(E))}{\sigma_{Ns_{i}}}\end{pmatrix}^{2},$
(S1)

where $N_{s_{i}}$ is the number of excess events, $N_{MC_{i}}(f(E))$ is the
number of signals expected by the MC simulation using the spectrum $f(E)$, and
$\sigma_{Ns_{i}}$ is the uncertainty of $N_{s_{i}}$ in the $i$th 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 S2
are taken into account. 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 (Non), the number of background events (Nb),
and the differential flux at the median energy (ELP for LP spectrum and EPLEC
for PLEC spectrum) of the bin, are listed in Table S1. According to our
observation of the Crab Nebula, the overall systematic uncertainty of the flux
is about 7%, and that of the spectral index is 0.02 (?).

Table S1: Number of events, background, median energy and corresponding differential flux from GRB 221009A Time after T0 | log(E${}_{\rm rec}/TeV$) | Non | Nb | ELP | Flux (LP/10-10) | EPLEC | Flux (PLEC/10-10)
---|---|---|---|---|---|---|---
|  |  |  | (TeV) | (ergs cm-2 s-1) | (TeV) | (ergs cm-2 s-1)
230-300s | 0.4$-$0.6 | 8 | 0.40 | 4.84 | 549${}_{-195}^{+240}$ | 4.32 | 548${}_{-194}^{+239}$
0.6$-$0.8 | 24 | 1.20 | 6.10 | 179${}_{-37}^{+42}$ | 5.31 | 207${}_{-42}^{+48}$
0.8$-$1.0 | 29 | 1.55 | 7.85 | 111${}_{-20}^{+23}$ | 6.53 | 156${}_{-28}^{+32}$
1.0$-$1.2 | 3 | 0.15 | 11.6 | 12.8${}_{-8.2}^{+10.3}$ | 8.61 | 28.5${}_{-18.4}^{+23.1}$
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${}_{-38}^{+45}$ | 4.12 | 100${}_{-37}^{+44}$
0.6$-$0.8 | 34 | 6.05 | 5.82 | 35.2${}_{-7.2}^{+8.0}$ | 4.95 | 41.0${}_{-8.3}^{+9.3}$
0.8$-$1.0 | 28 | 5.45 | 7.50 | 12.9${}_{-3.0}^{+3.3}$ | 5.96 | 20.4${}_{-4.7}^{+5.3}$
1.0$-$1.2 | 5 | 1.40 | 10.8 | 2.42${}_{-1.48}^{+1.89}$ | 8.04 | 5.46${}_{-3.35}^{+4.26}$
1.2$-$1.4 | 1 | 0.15 | 17.6 | 0.59${}_{-0.42}^{+1.22}$ | 11.9 | 3.01${}_{-2.16}^{+6.22}$
1.4$-$1.6 | 0 | 0 | 27.2 |  $<$2.2 | 17.2 |  $<$29

###  S1.3 The highest energy events observed by KM2A

During the period from T0+230s to T0+900s, 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
(Ne), the number of detected muons (Nμ), the incident zenith angle ($\theta$),
the distance of the shower core from the nearest edge of the active detector
array (Dedge), the space angle ($\Delta\psi$) between the event and the
direction of the GRB, the arrival time (Tevent) since T0, reconstructed energy
(ELP) and its errors using the LP spectrum shown in panel A of Figure 2,
reconstructed energy (EPLEC) and its errors using the PLEC spectrum shown in
panel B of Figure 2, and reconstructed energy (EEBL) 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 Tevent during the period
from T0+230s to T0+900s, ratio less than $log((N_{\mu}+0.0001)/N_{e})$, space
angle with the GRB less than ($\Delta\psi$), and reconstructed energy above
Erec 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.

Table S2: Detail information of the nine events with the highest energy from GRB 221009A Tevent(s) | ELP (TeV) | EPLEC (TeV) | EEBL (TeV) | Ne | Nμ |  $\theta$ (∘) |  $\Delta\psi$ (∘) | Dedge (m) | P (%)
---|---|---|---|---|---|---|---|---|---
236.6 | 12.7${}_{-3.8}^{+6.2}$ | 9.7${}_{-2.1}^{+3.3}$ | 9.8${}_{-2.3}^{+3.1}$ | 60.6 | 0 | 28.5 | 0.46 | 77 | 7.0
242.5 | 10.5${}_{-3.2}^{+5.0}$ | 8.3${}_{-2.1}^{+3.0}$ | 8.4${}_{-2.2}^{+3.2}$ | 57.4 | 0 | 28.8 | 0.45 | 111 | 10
262.4 | 12.6${}_{-3.8}^{+5.5}$ | 9.5${}_{-2.3}^{+3.4}$ | 9.6${}_{-2.4}^{+3.3}$ | 57.3 | 0 | 28.6 | 0.53 | 180 | 5.7
358.1 | 10.0${}_{-3.2}^{+4.8}$ | 7.4${}_{-1.8}^{+3.1}$ | 7.9${}_{-2.2}^{+3.3}$ | 46.0 | 0 | 28.7 | 0.54 | 119 | 6.0
571.1 | 9.4${}_{-3.0}^{+5.1}$ | 7.4${}_{-2.5}^{+2.6}$ | 7.7${}_{-2.5}^{+3.0}$ | 45.7 | 0 | 29.5 | 0.52 | 99 | 7.8
643.0 | 17.8${}_{-5.1}^{+7.4}$ | 12.2${}_{-2.4}^{+3.5}$ | 12.5${}_{-2.4}^{+3.2}$ | 81.8 | 0.3 | 29.7 | 0.62 | 181 | 4.5
812.4 | 11.1${}_{-4.3}^{+5.9}$ | 7.4${}_{-2.8}^{+3.6}$ | 7.6${}_{-3.0}^{+3.9}$ | 68.0 | 0 | 30.3 | 0.66 | 112 | 11
863.8 | 12.9${}_{-3.9}^{+6.1}$ | 9.2${}_{-2.8}^{+3.0}$ | 9.7${}_{-3.1}^{+3.2}$ | 100.2 | 0.8 | 30.1 | 1.07 | 81 | 17
894.1 | 13.6${}_{-4.2}^{+6.1}$ | 9.7${}_{-2.5}^{+3.4}$ | 10.4${}_{-3.0}^{+3.3}$ | 60.5 | 0 | 31.8 | 0.83 | 214 | 16

###  S1.4 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 (?), Gilmore et
al. 2012 (?), Dominguez et al. 2011 (?), and Finke et al. 2010 (?), 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 dN/dE =
J0Ea+b.log(E)e-τ(E) to the data. The intrinsic spectrum using different EBL
models is shown in Figure S5. The $\chi^{2}/ndf$ values for different EBL
models are listed in Table S3. For the interval from T0+230s to T0+300s, 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 T0+300s to T0+900s. The total $\chi^{2}/ndf$
after summing up the two intervals is comparable in the two cases.

Table S3: $\chi^{2}/ndf$ using different EBL models 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/9 | 5.49/10 | 11.42/19

Figure S5:  Panel (A) shows the EBL absorption models for VHE gamma-rays from
a redshift of z=0.151. Different lines represent different EBL models, i.e.,
Saldana-Lopez et al. 2021 (?), Gilmore et al. 2012 (?), Dominguez et al. 2011
(?), and Finke et al. 2010 (?). Panel (B) shows the intrinsic spectrum of GRB
221009A corrected for EBL absorption using the Gilmore et al. 2012 model. The
red points are for the interval from T0+230s to T0+300s, while the blue points
are for the interval from T0+300s to T0+900s. The filled points are obtained
using LHAASO-KM2A data, while the unfilled points are obtained using LHAASO-
WCDA data. The solid lines are the fitting result using the log-parabolic
function. Panels (C) and (D) show the intrinsic spectrum of GRB 221009A
corrected for EBL absorption using the EBL models of Dominguez et al. 2011 and
Finke et al. 2010, respectively.

To further test the EBL model, we divided the distribution of EBL adopted in
the Saldana-Lopez et al. 2021 (?) model into three wavelength ranges, i.e.,
$<$8 $\mu$m, 8 to 28 $\mu$m, and $>$28 $\mu$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.20}^{+0.33}$, 1.20${}_{-0.20}^{+0.23}$, and
0.40${}_{-0.16}^{+0.44}$, respectively. The LHAASO constrained EBL model is
plotted in Fig. S6. 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.

Figure S6:  The best fit EBL from the GRB 221009A measurement (red), relative
to the Saldana-Lopez et al. 2021 EBL model at z = 0 (blue), in three
wavelength bins. The data points represent direct EBL measurements, taken from
(?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?). Downward
triangles correspond to galaxy counts (?), which are taken as the lower limits
of EBL.

###  S1.5 Origin of the $\sim 10$ TeV photons

In the work on the LHAASO-WCDA result of GRB 221009A (?), the multi-wavelength
data, including Swift-XRT, Fermi-LAT, and LHAASO-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_{\rm
e}\propto\gamma_{\rm e}^{-p}$, where $\gamma_{\rm 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.

Figure S7:  Comparison between the SSC emission model (?) and the LHAASO data
for different EBL models. The LHAASO data for the two time intervals, obtained
using different EBL models, are same as those shown in Figure 3 and Figure S5.

The comparison between the SSC models used previously for WCDA data (?) and at
present for KM2A data is given in Figure S7. 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 (?). The peak energy from proton synchrotron emission can
reach $\sim 10$ TeV if the magnetic field equipartition factor is sufficiently
high in the reverse shock (?). For a proton spectrum $dN_{p}/dE\sim
E^{-\alpha_{p}}$, the spectral index of the energy flux ($\nu f_{\nu}$) from
proton synchrotron emission is $(3-\alpha_{p})/2$. With $\alpha_{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 (?, ?, ?). 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 photons do not suffer from significant EBL absorption,
the observed cascade radiation could have a hard spectrum (?). 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 (?, ?).

###  S1.6 Axion-like particle estimation

Figure S8:  Panel A shows EBL absorption models for very high-energy gamma-
rays from a redshift of $z=0.151$, taking into account the oscillation between
photons and axion-like particles (ALPs) assuming $m_{a}=10^{-7}$ eV and
$g_{a\gamma}$=(1 to 6)$\times 10^{-11}$ GeV-1. The EBL model used is Saldana
et al. 2021. Panel B shows $\Delta\chi^{2}$ relative to the minimum that fits
the spectral energy distribution data as a function of the ALP $g_{a\gamma}$
for $m_{a}=10^{-7}$ eV. The line indicates $\Delta\chi^{2}=2.71$ used to
define the upper limit on $g_{a\gamma}$.

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

$\displaystyle\mathcal{L}_{a\gamma\gamma}$
$\displaystyle=-\frac{1}{4}g_{a\gamma\gamma}aF_{\mu\nu}\tilde{F}^{\mu\nu}\ $
$\displaystyle=g_{a\gamma\gamma}a\bm{E}\cdot\bm{B}\ ,$ (S2)

where $a$ represents the ALP field, $F_{\mu\nu}$ and $\tilde{F}^{\mu\nu}$ are
the electromagnetic tensor and its dual, respectively, $\bm{E}$ and $\bm{B}$
are the electric and magnetic components, respectively, and
$g_{a\gamma\gamma}$ is the coupling constant. According to Eq. (S2), in the
presence of external magnetic fields, ALP and photon 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$\ddot{\rm o}$dinger-
like equation (?). 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~{}\mathrm{\mu G}$ with a
coherence length of $10~{}\mathrm{kpc}$ and an electron density of
approximately $0.04~{}\mathrm{cm}^{-3}$ (?). 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 (?). The electron
density of the Milky Way is described by the NE2001 model (?).

Given the circumstances established above and adopting the Saldana-Lopez et
al. EBL model (?), we solve the propagation equation and obtain the photon
survival probability at Earth. The left panel of Figure S8 shows the photon
survival probability for an axion mass $m_{a}=10^{-7}$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 photons
detected, and the fitting to data becomes worse, as shown in the right panel
of Figure S8. This gives us a constraint on the axion coupling, as shown in
Figure 4.

###  S1.7 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 (?):

$\epsilon_{\mathrm{thr}}=\frac{m_{e}^{2}c^{4}}{E}+\frac{1}{8}\left(\frac{E}{E_{\mathrm{LIV}}}\right)E\
,$ (S3)

where ELIV represents the LIV energy scale. From Eq. (S3), 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 MPl to 2.8 MPl, as
shown in the left panel of Figure S9. Similar to the axion case, if
$E_{\mathrm{LIV}}$ is too low, there will be too many high energy photons,
leading to conflicts with observation, as shown in the right panel of Figure
S9. This leads us to obtain a lower bound of the LIV energy scale at about 1.5
MPl.

Figure S9:  (A) EBL absorption models for VHE gamma-rays from a redshift of
z=0.151 taking into account the LIV assuming ELIV= (1 to 2.8)$\times$MPl. The
EBL model is from Saldana-Lopez et al. 2021. (B) $\Delta\chi^{2}$ relative to
the minimum from fitting the SED data as a function of ELIV. The line
indicates $\Delta\chi^{2}$=2.71 used to achieve the lower limit of
E${}_{\mathrm{LIV}}>$1.5MPl.