On 9 October 2022 the Large High Altitude Air Shower Observatory LHAASO reported the observation of the very early TeV afterglow of the brightestofalltime gammaray burst 221009A recording the highest photon statistics in the TeV band ever obtained from a gammaray burst. We use this unique observation to place stringent constraints on the energy dependence of the speed of light in vacuum a manifestation of Lorentz invariance violation LIV predicted by some quantum gravity QG theories. Our results show that the 95 confidence level lower limits on the QG energy scales are ESUBQG 1SUBgt10 times the Planck energy ESUBPlSUB for the linear LIV effect and ESUBQG 2SUBgt6 10SUP8SUPESUBPlSUB for the quadratic LIV effect. Our limits on the quadratic LIV case improve previous best bounds by factors of 57.
['Astrophysics - High Energy Astrophysical Phenomena', 'General Relativity and Quantum Cosmology', 'High Energy Physics - Phenomenology', 'High Energy Physics - Theory']
Stringent Tests of Lorentz Invariance Violation from LHAASO Observations of GRB 221009A
{'Stringent Tests of Lorentz Invariance Violation from LHAASO Observations of GRB 221009A': "Zhen Cao, 1, 2, 3 F. Aharonian, 4, 5 Axikegu, 6 Y.X. Bai, 1, 3 Y.W. Bao, 7 D. Bastieri, 8 X.J. Bi, 1, 2, 3 Y.J. Bi, 1, 3 W. Bian, 9 A.V. Bukevich, 10 Q. Cao, 11 W.Y. Cao, 12 Zhe Cao, 13, 12 J. Chang, 14 J.F. Chang, 1, 3, 13 A.M. Chen, 9 E.S. Chen, 1, 2, 3 H.X. Chen, 15 Liang Chen, 16 Lin Chen, 6 Long Chen, 6 M.J. Chen, 1, 3 M.L. Chen, 1, 3, 13 Q.H. Chen, 6 S. Chen, 17 S.H. Chen, 1, 2, 3 S.Z. Chen, 1, 3 T.L. Chen, 18 Y. Chen, 7 N. Cheng, 1, 3 Y.D. Cheng, 1, 2, 3 M.Y. Cui, 14 S.W. Cui, 11 X.H. Cui, 19 Y.D. Cui, 20 B.Z. Dai, 17 H.L. Dai, 1, 3, 13 Z.G. Dai, 12 Danzengluobu, 18 X.Q. Dong, 1, 2, 3 K.K. Duan, 14 J.H. Fan, 8 Y.Z. Fan, 14 J. Fang, 17 J.H. Fang, 15 K. Fang, 1, 3 C.F. Feng, 21 H. Feng, 1 L. Feng, 14 S.H. Feng, 1, 3 X.T. Feng, 21 Y. Feng, 15 Y.L. Feng, 18 S. Gabici, 22 B. Gao, 1, 3 C.D. Gao, 21 Q. Gao, 18 W. Gao, 1, 3 W.K. Gao, 1, 2, 3 M.M. Ge, 17 L.S. Geng, 1, 3 G. Giacinti, 9 G.H. Gong, 23 Q.B. Gou, 1, 3 M.H. Gu, 1, 3, 13 F.L. Guo, 16 X.L. Guo, 6 Y.Q. Guo, 1, 3 Y.Y. Guo, 14 Y.A. Han, 24 M. Hasan, 1, 2, 3 H.H. He, 1, 2, 3 H.N. He, 14 J.Y. He, 14 Y. He, 6 Y.K. Hor, 20 B.W. Hou, 1, 2, 3 C. Hou, 1, 3 X. Hou, 25 H.B. Hu, 1, 2, 3 Q. Hu, 12, 14 S.C. Hu, 1, 3, 26 D.H. Huang, 6 T.Q. Huang, 1, 3 W.J. Huang, 20 X.T. Huang, 21 X.Y. Huang, 14 Y. Huang, 1, 2, 3 X.L. Ji, 1, 3, 13 H.Y. Jia, 6 K. Jia, 21 K. Jiang, 13, 12 X.W. Jiang, 1, 3 Z.J. Jiang, 17 M. Jin, 6 M.M. Kang, 27 I. Karpikov, 10 D. Kuleshov, 10 K. Kurinov, 10 B.B. Li, 11 C.M. Li, 7 Cheng Li, 13, 12 Cong Li, 1, 3 D. Li, 1, 2, 3 F. Li, 1, 3, 13 H.B. Li, 1, 3 H.C. Li, 1, 3 Jian Li, 12 Jie Li, 1, 3, 13 K. Li, 1, 3 S.D. Li, 16, 2 W.L. Li, 21 W.L. Li, 9 X.R. Li, 1, 3 Xin Li, 13, 12 Y.Z. Li, 1, 2, 3 Zhe Li, 1, 3 Zhuo Li, 28 E.W. Liang, 29 Y.F. Liang, 29 S.J. Lin, 20 B. Liu, 12 C. Liu, 1, 3 D. Liu, 21 D.B. Liu, 9 H. Liu, 6 H.D. Liu, 24 J. Liu, 1, 3 J.L. Liu, 1, 3 M.Y. Liu, 18 R.Y. Liu, 7 S.M. Liu, 6 W. Liu, 1, 3 Y. Liu, 8 Y.N. Liu, 23 Q. Luo, 20 Y. Luo, 9 H.K. Lv, 1, 3 B.Q. Ma, 28 L.L. Ma, 1, 3 X.H. Ma, 1, 3 J.R. Mao, 25 Z. Min, 1, 3 W. Mitthumsiri, 30 H.J. Mu, 24 Y.C. Nan, 1, 3 A. Neronov, 22 L.J. Ou, 8 P. Pattarakijwanich, 30 Z.Y. Pei, 8 J.C. Qi, 1, 2, 3 M.Y. Qi, 1, 3 B.Q. Qiao, 1, 3 J.J. Qin, 12 A. Raza, 1, 2, 3 D. Ru ff olo, 30 A. S'aiz, 30 M. Saeed, 1, 2, 3 D. Semikoz, 22 L. Shao, 11 O. Shchegolev, 10, 31 X.D. Sheng, 1, 3 F.W. Shu, 32 H.C. Song, 28 Yu.V. Stenkin, 10, 31 V. Stepanov, 10 Y. Su, 14 D.X. Sun, 12, 14 Q.N. Sun, 6 X.N. Sun, 29 Z.B. Sun, 33 J. Takata, 34 P.H.T. Tam, 20 Q.W. Tang, 32 R. Tang, 9 Z.B. Tang, 13, 12 W.W. Tian, 2, 19 C. Wang, 33 C.B. Wang, 6 G.W. Wang, 12 H.G. Wang, 8 H.H. Wang, 20 J.C. Wang, 25 Kai Wang, 7 Kai Wang, 34 L.P. Wang, 1, 2, 3 L.Y. Wang, 1, 3 P.H. Wang, 6 R. Wang, 21 W. Wang, 20 X.G. Wang, 29 X.Y. Wang, 7 Y. Wang, 6 Y.D. Wang, 1, 3 Y.J. Wang, 1, 3 Z.H. Wang, 27 Z.X. Wang, 17 Zhen Wang, 9 Zheng Wang, 1, 3, 13 D.M. Wei, 14 J.J. Wei, 14 Y.J. Wei, 1, 2, 3 T. Wen, 17 C.Y. Wu, 1, 3 H.R. Wu, 1, 3 Q.W. Wu, 34 S. Wu, 1, 3 X.F. Wu, 14 Y.S. Wu, 12 S.Q. Xi, 1, 3 J. Xia, 12, 14 G.M. Xiang, 16, 2 D.X. Xiao, 11 G. Xiao, 1, 3 Y.L. Xin, 6 Y. Xing, 16 D.R. Xiong, 25 Z. Xiong, 1, 2, 3 D.L. Xu, 9 R.F. Xu, 1, 2, 3 R.X. Xu, 28 W.L. Xu, 27 L. Xue, 21 D.H. Yan, 17 J.Z. Yan, 14 T. Yan, 1, 3 C.W. Yang, 27 C.Y. Yang, 25 F. Yang, 11 F.F. Yang, 1, 3, 13 L.L. Yang, 20 M.J. Yang, 1, 3 R.Z. Yang, 12 W.X. Yang, 8 Y.H. Yao, 1, 3 Z.G. Yao, 1, 3 L.Q. Yin, 1, 3 N. Yin, 21 X.H. You, 1, 3 Z.Y. You, 1, 3 Y.H. Yu, 12 Q. Yuan, 14 H. Yue, 1, 2, 3 H.D. Zeng, 14 T.X. Zeng, 1, 3, 13 W. Zeng, 17 M. Zha, 1, 3 B.B. Zhang, 7 F. Zhang, 6 H. Zhang, 9 H.M. Zhang, 7 H.Y. Zhang, 1, 3 J.L. Zhang, 19 Li Zhang, 17 P.F. Zhang, 17 P.P. Zhang, 12, 14 R. Zhang, 12, 14 S.B. Zhang, 2, 19 S.R. Zhang, 11 S.S. Zhang, 1, 3 X. Zhang, 7 X.P. Zhang, 1, 3 Y.F. Zhang, 6 Yi Zhang, 1, 14 Yong Zhang, 1, 3 B. Zhao, 6 J. Zhao, 1, 3 L. Zhao, 13, 12 L.Z. Zhao, 11 S.P. Zhao, 14 X.H. Zhao, 25 F. Zheng, 33 W.J. Zhong, 7 B. Zhou, 1, 3 H. Zhou, 9 J.N. Zhou, 16 M. Zhou, 32 P. Zhou, 7 R. Zhou, 27 X.X. Zhou, 1, 2, 3 X.X. Zhou, 6 B.Y. Zhu, 12, 14 C.G. Zhu, 21 F.R. Zhu, 6 H. Zhu, 19 K.J. Zhu, 1, 2, 3, 13 Y.C. Zou, 34 and X. Zuo 1, 3  \n(The LHAASO Collaboration) ∗  \n1 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, 610000 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 School of Physical Science and Technology & School of Information Science and Technology, Southwest Jiaotong University, 610031 Chengdu, Sichuan, China 7 School of Astronomy and Space Science, Nanjing University, 210023 Nanjing, Jiangsu, China 8 Center for Astrophysics, Guangzhou University, 510006 Guangzhou, Guangdong, China 9 Tsung-Dao Lee Institute & School of Physics and Astronomy, Shanghai Jiao Tong University, 200240 Shanghai, China 10 Institute for Nuclear Research of Russian Academy of Sciences, 117312 Moscow, Russia 11 Hebei Normal University, 050024 Shijiazhuang, Hebei, China 12 University of Science and Technology of China, 230026 Hefei, Anhui, China 13 State Key Laboratory of Particle Detection and Electronics, China 14 Key Laboratory of Dark Matter and Space Astronomy & Key Laboratory of Radio Astronomy, Purple Mountain Observatory, Chinese Academy of Sciences, 210023 Nanjing, Jiangsu, China 15 Research Center for Astronomical Computing, Zhejiang Laboratory, 311121 Hangzhou, Zhejiang, China 16 Key Laboratory for Research in Galaxies and Cosmology, Shanghai Astronomical Observatory, Chinese Academy of Sciences, 200030 Shanghai, China 17 School of Physics and Astronomy, Yunnan University, 650091 Kunming, Yunnan, China 18 Key Laboratory of Cosmic Rays (Tibet University), Ministry of Education, 850000 Lhasa, Tibet, China 19 National Astronomical Observatories, Chinese Academy of Sciences, 100101 Beijing, China  \n21  \n20 School of Physics and Astronomy (Zhuhai) & School of Physics (Guangzhou)  \n& Sino-French Institute of Nuclear Engineering and Technology (Zhuhai), Sun Yat-sen University, 519000 Zhuhai & 510275 Guangzhou, Guangdong, China Institute of Frontier and Interdisciplinary Science, Shandong University, 266237 Qingdao, Shandong, China 22 APC, Universit'e Paris Cit'e, CNRS / IN2P3, CEA / IRFU, Observatoire de Paris, 119 75205 Paris, France 23 Department of Engineering Physics, Tsinghua University, 100084 Beijing, China 24 School of Physics and Microelectronics, Zhengzhou University, 450001 Zhengzhou, Henan, China 25 Yunnan Observatories, Chinese Academy of Sciences, 650216 Kunming, Yunnan, China 26 China Center of Advanced Science and Technology, Beijing 100190, China 27 College of Physics, Sichuan University, 610065 Chengdu, Sichuan, China 28 School of Physics, Peking University, 100871 Beijing, China 29 Guangxi Key Laboratory for Relativistic Astrophysics, School of Physical Science and Technology, Guangxi University, 530004 Nanning, Guangxi, China 30 Department of Physics, Faculty of Science, Mahidol University, Bangkok 10400, Thailand 31 Moscow Institute of Physics and Technology, 141700 Moscow, Russia 32 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 33 National Space Science Center, Chinese Academy of Sciences, 100190 Beijing, China 34 School of Physics, Huazhong University of Science and Technology, Wuhan 430074, Hubei, China  \n(Dated: August 27, 2024)  \nOn October 9, 2022, the Large High Altitude Air Shower Observatory (LHAASO) reported the observation of the very early TeV afterglow of the brightest-of-all-time GRB 221009A, recording the highest photon statistics in the TeV band ever obtained from a Gamma-ray burst. We use this unique observation to place stringent constraints on an energy dependence of the speed of light in vacuum, a manifestation of Lorentz invariance violation (LIV) predicted by some quantum gravity (QG) theories. Our results show that the 95% confidence level lower limits on the QG energy scales are E QG , 1 > 10 times of the Planck energy E Pl for the linear, and E QG , 2 > 6 × 10 -8 E Pl for the quadratic LIV e ff ects, respectively. Our limits on the quadratic LIV case improve previous best bounds by factors of 5-7.", 'I. INTRODUCTION': 'Lorentz invariance, the fundamental symmetry of Einstein\'s relativity, has withstood strict tests over the past century [1]. However, deviations from Lorentz invariance at energies approaching the Planck scale E Pl = p ℏ c 5 / G ≃ 1 . 22 × 10 19 GeV are predicted in many quantum gravity (QG) theories seeking to unify quantum theory and general relativity [2-15]. Although any signals of Lorentz invariance violation (LIV) are expected to be very tiny at attainable energies ≪ E Pl, they can increase with energy and accumulate to detectable levels over large propagation distances. Astrophysical observations involving high-energy radiation and long distances are therefore suitable for performing sensitive tests of Lorentz invariance.  \nOne of the manifestations of LIV can be characterised as energy-dependent modifications to the photon dispersion relation in vacuum [16]:  \nE 2 ≃ p 2 c 2        1 -∞ X n = 1 s E E QG , n ! n        , (1)  \nwhere E and p are the energy and momentum of a photon, s = ± 1 is the \'sign\' of the LIV e ff ect, corresponding to the  \n\'subluminal\' or \'superluminal\' scenarios, and E QG , n denotes the hypothetical QG energy scale. Since the sum is dominated by the lowest-order term of the series at small energies E ≪ E QG , n , only the first two leading terms ( n = 1 or n = 2) are of interest for independently LIV tests. They are usually referred to as linear and quadratic LIV corrections, respectively. Taking into account only the leading LIV modification of order n , the photon group velocity is then given by  \nυ ( E ) = ∂ E ∂ p ≈ c " 1 -s n + 1 2 E E QG , n ! n # . (2)  \nBecause of the energy dependence of υ ( E ), two photons with di ff erent energies (denoted by Eh and El , where Eh > El ) emitted simultaneously from the same source at redshift z would reach us at di ff erent times. The energy-dependent time delay due to LIV e ff ects can be expressed as [17]  \n∆ t LIV = s n + 1 2 E n h -E n l E n QG , n Z z 0 (1 + z \' ) n H ( z \' ) dz \' , (3)  \nwhere H ( z ) = H 0 p Ω m(1 + z ) 3 + ΩΛ , assuming a flat Λ CDM cosmology with Hubble constant H 0 = 67 . 36 km s -1 Mpc -1 , matter density parameter Ω m = 0 . 315, and vacuum energy density ΩΛ = 1 -Ω m [18]. For convenience, in Eq. (3) we introduce the dimensionless LIV parameters  \nη 1 = sE Pl / E QG , 1 (4)  \nand  \nη 2 = 10 -15 × sE 2 Pl / E 2 QG , 2 (5)  \nfor linear ( n = 1) and quadratic ( n = 2) modifications, respectively, to replace E QG , 1 and E QG , 2.  \nIt is obvious from Eq. (3) that the greatest sensitivities on η n (or E QG , n ) can be expected from those astrophysical sources with rapid signal variability, large distances, and high-energy emission. As the most violent explosions occurring at cosmological distances, Gamma-ray bursts (GRBs) have been deemed as excellent probes for searching for the LIV-induced vacuum dispersion [16, 19-22]. Indeed, the most stringent limits to date on η n (or E QG , n ), resulting from vacuum dispersion time-of-flight studies, have been obtained using the GeV emission by GRB 090510 observed by the Fermi-LAT. The limits set for the subluminal (superluminal) scenario are η 1 < 0 . 13, or equivalently E QG , 1 > 9 . 3 × 10 19 GeV ( η 1 > -0 . 09, or equivalently E QG , 1 > 1 . 3 × 10 20 GeV) and η 2 < 8 . 8, or equivalently E QG , 2 > 1 . 3 × 10 11 GeV ( η 2 > -16 . 8, or equivalently E QG , 2 > 9 . 4 × 10 10 GeV) for linear and quadratic LIV e ff ects, respectively [21]. Based on the detection of sub-TeV emission from GRB 190114C, MAGIC Collaboration (2020) [23] obtained competitive lower limits on the quadratic LIV energy scale, i.e., η 2 < 37 . 0, or equivalently E QG , 2 > 6 . 3 × 10 10 GeV ( η 2 > -48 . 0, or equivalently E QG , 2 > 5 . 6 × 10 10 GeV) for the subluminal (superluminal) case.  \nThe Large High Altitude Air Shower Observatory (LHAASO) detected more than 64,000 photons in the energy range of 0.2-7 TeV from GRB 221009A within the first 4000 s after the MeV burst trigger [24]. This object is located at redshift z = 0 . 151 [25, 26]. In this Letter, we study Lorentzviolating e ff ects using the time-of-flight measurements of the unprecedentedly very-high-energy (VHE, > 100 GeV) photons from GRB 221009A.', 'II. LHAASO OBSERVATIONS OF GRB 221009A': 'At 13:16:59.99 UTC on 9 October 2022 (denoted as T 0), a long-duration GRB, numbered as GRB 221009A, triggered the Gamma-ray Burst Monitor (GBM) onboard the Fermi satellite [27, 28]. The subsequent detection with Fermi-LAT made clear that it is an extraordinarily bright burst [29, 30]. The Gamma-ray emission of GRB 221009A was also detected by several other space missions [31-45], and by the groundbased air shower detector LHAASO [46].  \nLHAASO [47] is a new generation Gamma-ray and cosmic-ray observatory situated in Daocheng, China, at an elevation of ∼ 4410 meters. Due to the large area, wide field of view, and broad energy coverage, the LHAASO detectors are meticulously designed to delve into new frontiers physics, including investigations into LIV, among other scientific objectives.  \nAt T 0, GRB 221009A was observed by LHAASO at a zenith angle of 28 . 1 · , and remained within LHAASO\'s field of view for the next 6000 seconds. In the initial 4000 seconds, the Water Cherenkov Detector Array (WCDA) of LHAASO captured over 64,000 photons in the 0.2-7 TeV energy range, both the light curve and energy spectrum of VHE photons were measured [24]. The intrinsic light curve reveals a rise to a peak from 231 to 244 seconds after T 0, followed by a  \ndecay lasting 650 seconds.  \nThe light curve of energy flux in the specified time range, as detected by LHAASO-WCDA, can be well described by a smoothly broken power-law function [24],  \nλ ( t ) ∝ " t t b ! -ωα 1 + t t b ! -ωα 2 # -1 /ω , (6)  \nwhere all time-related variables are relative to a reference time T ∗ = T 0 + 226 s. Here, α 1 = 1 . 82 and α 2 = -1 . 115 denote the power-law indices before and after the break time t b = 15 . 37 s, and ω = 1 . 07 represents the sharpness of the break. The intrinsic time-resolved spectrum can be fitted with a power-law function, and the positive power-law spectral index varies with time following the expression [24]  \nγ ( t ) = a log( t ) + b , (7)  \nwhere the unit of t is seconds, and a = -0 . 135 and b = 2 . 579. When a time delay due to LIV is introduced, this λ ( t ) will be modified to λ GLYPH<2> t -∆ t LIV( E , η n ) GLYPH<3> .  \nThe observed count rate light curve, characterizing the probability of observing a photon in a range of the number of hits ∆ N hit and the arrival time t from the GRB, can be converted from the energy flux light curve with  \nf ( t , ∆ N hit | η n ) = Z + ∞ 0 λ GLYPH<2> t -∆ t LIV( E , η n ) GLYPH<3> × ζ ( t ) E -γ ( t ) P EBL( E ) S ( E , t , ∆ N hit) d E . (8)  \nHere, ζ ( t ) = A / R E 2 E 1 E 1 -γ ( t ) d E serves as the conversion factor from energy flux (within the energy range from E 1 = 0 . 3 TeV to E 2 = 5 TeV) to flux coe ffi cient, where γ ( t ) represents the power-law index evolution as per Eq. (7), and A is a constant to be determined. The term P EBL( E ) denotes the survival probability of photons subject to extra-galactic background light (EBL) attenuation, adopting the model [48]. Lastly, S ( E , t , ∆ N hit) accounts for the e ff ective detection area of photons at energy E and time t , with the number of fired cells in a given segment ∆ N hit.', 'III. ANALYSIS METHODS AND RESULTS': 'We utilize two analysis methods to examine the LIV lags in the VHE Gamma-ray signals from GRB 221009A. The cross-correlation function (CCF) is employed to directly measure the time delays between di ff erent energy bands, while the maximum likelihood (ML) method is adopted to extract energy-dependent arrival time delays. These two methods are widely employed in similar LIV studies.', 'A. Cross-correlation Function Method': 'We segment the light curve of the count rate detected by LHAASO-WDA from GRB 221009A into ten intervals, covering the time span from 232 to 400 seconds after T 0. The segmentation is based on the number of fired cells ( N hit), with approximately the same number of events in each segment. The  \nFIG. 1: Count rate light curves of GRB 221009A, as detected by LHAASO-WCDA, presented in ten N hit segments. The time binning of the light curves used for analysis is 0.1 s; however, they are depicted here with 2 s intervals for clarity. The blue-marked range 232-400 s after the GBM trigger T 0 is selected for calculating the time lags.  \n<!-- image -->  \nN hit segments roughly correspond to di ff erent energy ranges, and the median energies are 0.354, 0.375, 0.395, 0.419, 0.457, 0.486, 0.556, 0.658, 0.843, and 1.601 TeV, respectively, considering the spectral index evolution from [24] for the interested time span. The energy-dependent light curves of GRB 221009A for the ten N hit segments (denoted by Seg0-Seg9) are displayed in Fig. 1.  \nIn our analysis, we use the the CCF method to calculate the time lags ∆ t between the lowest energy band (Seg0) and any of the other nine high energy bands (Seg1-Seg9) (see Section A of the Supplemental Material for details). Assuming that the observed time lags ∆ t are primarily caused by LIV e ff ects, we can establish a conservative limit on the LIV parameter η n ( n = 1 or 2). Utilizing the 9 pairs of CCF measurements, we conduct a global fit to constrain η n by minimizing  \nχ 2 ( η n ) = X j h ∆ t j -∆ t LIV( η n ) i 2 σ 2 ( ∆ t j ) , (9)  \nwhere σ ( ∆ t j ) is the uncertainty of ∆ t j , regarded as a statistical origin. This uncertainty is obtained by a bootstrapping method: mocking the data for this CCF pair 1000 times, and RMS of the obtained ∆ t j is set to σ ( ∆ t j ).  \nIn Eq. (9), the median energy of photons is used to calculate the LIV time delay ∆ t LIV( η n ) for each N hit segment. However,  \nFIG. 2: Relationship between the induced LIV value of η n (with η 1 in the left panel and η 2 in the right panel) and the input η n (denoted as η 1 and η 2 respectively), obtained through a simulation procedure.  \n<!-- image -->  \nFIG. 3: χ 2 as a function of the dimensionless LIV parameter η n for the linear case (on the left) and the quadratic case (on the right), utilizing the CCF method. The vertical dashed lines indicate the best fits, and the shaded areas represent the 95% confidence intervals.  \n<!-- image -->  \nthis method may introduce bias due to several factors: spectral index evolution, wide dispersion of photon energies within each N hit segment (mainly caused by air shower fluctuations), and significant overlap in energy ranges among adjacent segments (see Section B of the Supplemental Material).  \nThe biases are estimated as follows. Assuming a LIV parameter η 1 or η 2, the intrinsic light curve of energy flux as Eq. (6) can be simulated, and then translated by Eq. (8) into ten count rate light curves. Employing simulated light curves, the same procedure as with the data is executed, leading to the determination of the measured LIV parameters, as shown in Fig. 2. This involves iterating over various η 1 and η 2 values as input, enabling a polynomial fitting. Using the fitted function, the bias η n -η n for the analysis results on experimental data can be evaluated and subtracted.  \nThe χ 2 distribution as a function of unbiased η n is calculated and depicted in Fig. 3. The MINUIT package [49], along with its error estimation processor MINOS [50], is used for fitting and constructing the 95% confidence levels (by setting χ 2 value 3.84 up the minimum). MINOS accounts for parameter correlations and non-linearities, providing asymmetric er-  \nintervals with correct probability coverage for both χ 2 and likelihood fits. The best-fit values of η n and their uncertainties are translated into limits on the QG energy scale E QG at the 95% confidence level. These limits are detailed in Table I.', 'B. Maximum Likelihood Method': "Another approach for inferring η n is the ML method, using Eq. (8) as the function to describe the count rate light curve. A binned Poissonian likelihood method is employed for the analysis, with a time binning set to 0.1 s and N hit split into 10 segments as illustrated in Fig. 1. The background distribution, as a function of time for each N hit segment, is determined using data from the same transit as the GRB, encompassing two sidereal days before and after the burst. The Poisson probability is given by  \nP i , j = e -( µ b , i , j + µ s , i , j ) ( µ b , i , j + µ s , i , j ) N on , i , j N on , i , j ! , (10)  \nwhere index i represents time bins, and j represents N hit segments. The number of signals from the GRB µ s , i , j is calculated using Eq. (8), the number of background events µ b , i , j is evaluated from the polynomial fitting of the background as a function of time (see Section C of the Supplemental Material), and N on , i , j is the number of observed events. The logarithmic likelihood ratio is defined as  \nΛ = -2 ln L 0 L = -2 X i , j        GLYPH<16> µ ' s , i , j -µ s , i , j GLYPH<17> + N on , i , j ln µ b , i , j + µ s , i , j µ b , i , j + µ ' s , i , j        . (11)  \nThis ratio is utilized for minimization to determine the LIV parameter η n and all other light curve parameters. The reference time ( T ∗ , see the context of Eq. (6)) and two parameters for spectral index evolution ( a and b , see Eq. (7)) are fixed to the values provided in [24]. In this equation, L = Q i , j P i , j is the likelihood function, and L 0 is the one for a null hypothesis that no LIV delay exists, where the number of signals corresponds to µ ' s i j through setting ∆ t LIV ( E , η n ) = 0 in Eq. (8).  \n, , The MINUIT package and MINOS processor is utilized for the likelihood fitting and calculating errors. The best-fit values are η 1 = 0 . 066 and η 2 = 0 . 20.  \nThere could be biases in this approach for η 1 and η 2. For instance, Eq. (8) does not fit the data well due to the stochastic nature of the afterglow process and various disruptive phenomena involved. To address this, we analyzed 1000 shu ffl ed data sets by randomly exchanging the time and N hit information of events to decouple the correlation between time and energy. To preserve the behavior of the spectral index evolution, we locally shu ffl ed the events within each time bin of 6 seconds and reshu ffl ed them after shifting half of the time bin. The chosen time binning for shu ffl ing is su ffi ciently large, as the maximum time delay is less than 6 seconds within the energy coverage of the data for | η 1 | < 0 . 15 or | η 2 | < 0 . 35. The means of the 1000 best-fit values of η n from the shu ffl ed  \nFIG. 4: Distribution of best-fit values of η n from the shu ffl ed data, presented for both linear (on the left) and quadratic (on the right) cases. The vertical lines indicate the means.  \n<!-- image -->  \nFIG. 5: Profile likelihood distributions resulting from the analysis with the ML method, depicted for both linear (on the left) and quadratic (on the right) cases. The vertical dashed lines signify the best fits after bias subtraction, and the shaded areas correspond to the 95% confidence intervals.  \n<!-- image -->  \ndata are considered as the biases. As shown in Fig. 4, the mean values are η 1 , bias = 0 . 063 and η 2 , bias = 0 . 19. The profile likelihood curves after subtracting these biases are depicted in Fig. 5.  \nAs no significant LIV time delay from GRB 221009A are detected in this approach, we set limits for η n by constructing 95% confidence intervals with the assistance of MINOS. The obtained intervals are subsequently used to compute the limits on the QG energy scale E QG. The corresponding results are listed in Table I.  \nAn alternative approach via bootstrapping and shu ffl ing, as developed in [21, 23], is also employed to obtain the so-called 'calibrated' limits. As shown in the last rows of Table I, the results appear to be very similar.  \nThe EBL model could introduce systematic uncertainties in our analysis. We conducted another two rounds of analyses for the same data, considering two extreme cases of the EBL models (corresponding to the upper and lower boundaries of the uncertainty of the EBL model in [48]). For the linear LIV e ff ect, we observe that the EBL models would enlarge the η 1 limits by 18% (12%) in subluminal (superluminal) scenario. In the quadratic case, the η 2 limits would be reduced by 6% (5%) in subluminal (superluminal) scenario.  \nTABLE I: Values for the best fits (BF) and the 95% lower (LL) and upper (UL) limits, provided for the dimensionless LIV parameter η n using both the CCF and ML methods. Additionally, the 95% confidence level (CL) lower limits on the quantum gravity (QG) energy scale E QG for the linear ( n = 1) and quadratic ( n = 2) cases are listed. a", 'IV. SUMMARY': "LHAASO observed unprecedented large number of VHE photon events from the brightest GRB 221009A at the earliest epoch, marking the identification of the onset of a TeV GRB afterglow for the first time. These characteristics render this signal a unique opportunity to probe LIV in the photon sector. Utilizing both CCF and ML methods, we searched for LIVinduced lags in the arrival time of the energetic photons. In both methods, compatible limits on the LIV energy scale E QG are obtained.  \nOur limit on the linear modification of the photon dispersion relation is E QG , 1 > 1 . 0 × 10 20 GeV ( E QG , 1 > 1 . 1 × 10 20 GeV), considering a subluminal (superluminal) LIV effect. This is comparable to the most constraining lower limit on E QG , 1 obtained by the GeV emission of GRB 090510 [21]. In the quadratic case, our result on the energy scale E QG , 2 > 6 . 9 × 10 11 GeV ( E QG , 2 > 7 . 0 × 10 11 GeV) for a subluminal (superluminal) LIV e ff ect represents the best time-of-flight limit, improving previous bounds [21] by a factor of 5 (7).  \nMoreover, we emphasize that, thanks to the adoption of the true spectral time lags of bunches of high-energy photons, our constraints on E QG , n could be statistically more robust compared to previous results of GRB 090510, which were based on the rough time lags of sparse GeV-scale photons.  \n- [1] V. A. Kosteleck'y and N. Russell, Reviews of Modern Physics 83 , 11 (2011), arXiv:0801.0287 [hep-ph].  \n- [3] G. Amelino-Camelia, J. Ellis, N. E. Mavromatos, and D. V. Nanopoulos, International Journal of Modern Physics A 12 , 607 (1997), arXiv:hep-th / 9605211 [hep-th].  \n- [5] G. Amelino-Camelia, Nature (London) 418 , 34 (2002), arXiv:gr-qc / 0207049 [gr-qc].\n- [6] J. Magueijo and L. Smolin, Phys. Rev. Lett. 88 , 190403 (2002), arXiv:hep-th / 0112090 [hep-th].  \nFuture observations of VHE prompt emission instead of afterglow from GRBs would further enhance sensitivity to LIV e ff ects using vacuum dispersion (time-of-flight) tests.", 'ACKNOWLEDGMENTS': "We would like to thank all the members who constructed and operate the LHAASO detectors. This work is supported by the following grants: the National Natural Science Foundation of China (Nos. U1831208, 12005246, 12173039, 12321003, 12373053, and 12375108), the Strategic Priority Research Program of the Chinese Academy of Sciences (grant No. XDB0550400), the 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, the Key Research Program of Frontier Sciences of Chinese Academy of Sciences No. ZDBS-LY-7014, the Natural Science Foundation of Jiangsu Province No. BK20221562, and in Thailand by the National Science and Technology Development Agency (NSTDA) and National Research Council of Thailand (NRCT): High-Potential Research Team Grant Program (N42A650868).  \n- [7] J. Alfaro, H. A. Morales-T'ecotl, and L. F. Urrutia, Phys. Rev. D 65 , 103509 (2002), arXiv:hep-th / 0108061 [hep-th].\n- [8] D. Mattingly, Living Reviews in Relativity 8 , 5 (2005), arXiv:gr-qc / 0502097 [gr-qc].\n- [9] T. Li, N. E. Mavromatos, D. V. Nanopoulos, and D. Xie, Physics Letters B 679 , 407 (2009), arXiv:0903.1303 [hep-th].\n- [10] G. Amelino-Camelia, Living Reviews in Relativity 16 , 5 (2013), arXiv:0806.0339 [gr-qc].\n- [11] J. D. Tasson, Reports on Progress in Physics 77 , 062901 (2014), arXiv:1403.7785 [hep-ph].\n- [12] J.-J. Wei and X.-F. Wu, Frontiers of Physics 16 , 44300 (2021), arXiv:2102.03724 [astro-ph.HE].\n- [13] P. He and B.-Q. Ma, Universe 8 , 323 (2022), arXiv:2206.08180\n- [astro-ph.HE].\n- [14] A. Addazi, et al. , Progress in Particle and Nuclear Physics 125 , 103948 (2022), arXiv:2111.05659 [hep-ph].\n- [15] R. Alves Batista, et al. , arXiv e-prints , arXiv:2312.00409 (2023), arXiv:2312.00409 [gr-qc].\n- [16] G. Amelino-Camelia, J. Ellis, N. E. Mavromatos, D. V. Nanopoulos, and S. Sarkar, Nature (London) 393 , 763 (1998), arXiv:astro-ph / 9712103 [astro-ph].\n- [17] U. Jacob and T. Piran, JCAP 2008 , 031 (2008), arXiv:0712.2170 [astro-ph].\n- [18] Planck Collaboration, et al. , Astron. & Astrophys. 641 , A6 (2020), arXiv:1807.06209 [astro-ph.CO].\n- [19] A. A. Abdo, et al. , Science 323 , 1688 (2009).\n- [20] A. A. Abdo, et al. , Nature (London) 462 , 331 (2009), arXiv:0908.1832 [astro-ph.HE].\n- [21] V. Vasileiou, et al. , Phys. Rev. D 87 , 122001 (2013), arXiv:1305.3463 [astro-ph.HE].\n- [22] J. Ellis, et al. , Phys. Rev. D 99 , 083009 (2019), arXiv:1807.00189 [astro-ph.HE].\n- [23] V. A. Acciari, et al. (MAGIC Collaboration), Phys. Rev. Lett. 125 , 021301 (2020).\n- [24] LHAASO Collaboration, et al. , Science 380 , 1390 (2023), arXiv:2306.06372 [astro-ph.HE].\n- [25] A. de Ugarte Postigo, et al. , GRB Coordinates Network 32648 , 1 (2022).\n- [26] D. B. Malesani, et al. , arXiv e-prints , arXiv:2302.07891 (2023), arXiv:2302.07891 [astro-ph.HE].\n- [27] P. Veres, et al. , GRB Coordinates Network 32636 , 1 (2022).\n- [28] S. Lesage, et al. , Astrophys. J. Lett. 952 , L42 (2023), arXiv:2303.14172 [astro-ph.HE].\n- [29] E. Bissaldi, N. Omodei, M. Kerr, and Fermi-LAT Team, GRB Coordinates Network 32637 , 1 (2022).\n- [30] R. Pillera, et al. , GRB Coordinates Network 32658 , 1 (2022).\n- [31] S. Dichiara, et al. , GRB Coordinates Network 32632 , 1 (2022).\n- [32] H. A. Krimm, et al. , GRB Coordinates Network 32688 , 1 (2022).\n- [33] A. Ursi, et al. , GRB Coordinates Network 32650 , 1 (2022).\n- [34] G. Piano, et al. , GRB Coordinates Network 32657 , 1 (2022).\n- [35] D. Gotz, et al. , GRB Coordinates Network 32660 , 1 (2022).\n- [36] H. Xiao, S. Krucker, and R. Daniel, GRB Coordinates Network 32661 , 1 (2022).\n- [37] I. Lapshov, et al. , GRB Coordinates Network 32663 , 1 (2022).\n- [38] D. Frederiks, et al. , GRB Coordinates Network 32668 , 1 (2022).\n- [39] D. Frederiks, et al. , Astrophys. J. Lett. 949 , L7 (2023), arXiv:2302.13383 [astro-ph.HE].\n- [40] J. Ripa, et al. , GRB Coordinates Network 32685 , 1 (2022).\n- [41] J. ˇ R'ıpa, et al. , Astron. & Astrophys. 677 , L2 (2023), arXiv:2302.10047 [astro-ph.HE].\n- [42] L. J. Mitchell, B. F. Phlips, and W. N. Johnson, GRB Coordinates Network 32746 , 1 (2022).\n- [43] J. C. Liu, et al. , GRB Coordinates Network 32751 , 1 (2022).\n- [44] Z.-H. An, et al. , arXiv e-prints , arXiv:2303.01203 (2023), arXiv:2303.01203 [astro-ph.HE].\n- [45] K.-K. Duan, et al. , GRB Coordinates Network 32973 , 1 (2022).\n- [46] Y. Huang, et al. , GRB Coordinates Network 32677 , 1 (2022).\n- [47] X.-H. Ma, et al. , Chinese Physics C 46 , 030001 (2022).\n- [48] A. Saldana-Lopez, et al. , Mon. Not. R. Astron. Soc. 507 , 5144 (2021), arXiv:2012.03035 [astro-ph.CO].\n- [49] F. James and M. Roos, Computer Physics Communications 10 , 343 (1975).\n- [50] B. A. Murtagh and M. A. Saunders, (1983), https: // web.stanford.edu / group / SOL / guides / minos551.pdf.\n- [51] T. Piran and D. D. Ofengeim, Phys. Rev. D 109 , L081501", '(2024).': '- [52] T. N. Ukwatta, et al. , Astrophys. J. 711 , 1073 (2010), arXiv:0908.2370 [astro-ph.HE].\n- [53] R.-J. Lu, et al. , Astrophys. J. 865 , 153 (2018), arXiv:1808.00636 [astro-ph.HE].', 'Appendix A: The Cross-correlation Function and time Lag': 'For a given pair of light curves, we utilize the crosscorrelation function (CCF) to determine the time lag ∆ t . The CCF with a delay ∆ t is defined as  \nFj ( ∆ t ) = P i R 0( ti ) Rj ( ti + ∆ t ) q P i R 2 0 ( ti ) P i R 2 j ( ti ) , (S1)  \nwhere j = 1 , 2 , 3 , . . . , 9 represents the N hit segment number except the first one, i is the time bin number, and R 0( ti ) and Rj ( ti ) are the rates of segment 0 and j at time bin ti . A discrete analysis is performed with a time step 0.1 s. The corresponding CCF as a function of time delay for the two light curves between the lowest energy band (Seg0) and any of the other nine high energy bands (Seg1-Seg9) are shown in Figure S1. We define ∆ t as the time delay corresponding to the global maximum of the CCF. By fitting Fj ( ∆ t ) around the peak with a Gaussian function (see the orange dashed line), a more precise peak position ∆ t j can be obtained. From these plots we would say that the distribution around the peak behaves like a Gaussian, and the fitting is quite reasonable, where the obtained peak position reflects the mean behaviour. Some other similar analyses such as [52, 53] also took this kind of treatment.  \nWe employ the bootstrap method to estimate the variance of the time lag ∆ t , whose square root is taken as the error. The details are as follows. Firstly, we randomly sample these 10 light curves of data using the Monte Carlo (MC) method to obtain 10 mocked distributions, where the size of each MC sample is same as the data. Secondly, using the same CCF analysis method, we calculate the ∆ t for the mocked distributions. Lastly, we repeat above two procedures for 1000 times, and then obtain 1000 ∆ t values, whose variance is finally calculated. This kind of treatment is also adopted by other similar analysis, such as [52]. We assume the systematic uncertainty of ∆ t plays negligible role in the fitting, and only statistical fluctuation takes part. This is the prerequisites of the bootstrap method that we applied.  \nIn our CCF analysis method, we extract 10 energydependent light curves of GRB 221009A, and the time binning of the light curves used for analysis is 0.1 s. In order to understand the e ff ects of energy and time binning more fully, we have ever investigated the cases for three time-bin sizes of 0.1, 0.5, and 1.0 s, and two N hit-segment numbers of 5 and 10,  \nand found that the di ff erences in the best-fitted η n (the dimensionless LIV parameters) values of the CCF are less than 2%, and the di ff erences in the uncertainty of the time lag are less than 5%. Tests show that the choice of energy and time bins for the light curves has negligible impact on the results.', 'Appendix B: Bias from the Energy Overlapping of Di ff erent N hit Segments': 'The energy ranges of ten N hit segments are listed in Table S1. Due to large fluctuation of the cascade processes of low energy particles / photons in the atmosphere, the energy TABLE S1: Energy ranges of 10 N hit segments. The lower and upper limits represent the 16% and 84% percentiles of the energy distribution, respectively.  \nranges of neighbouring N hit segments are much overlapped. Evaluating the LIV time delay ∆ t LIV( η n ) based on the median energy may introduce bias, since photon energies exhibit wide dispersion within each N hit segment. The procedure of estimating this bias is presented in Section III of the main letter.', 'Appendix C: Polynomial Fitting of the Background': 'The polynomial fit to the background has also been adopted in the previous study that was published in [24]. Two examples for the background fitting of Segments 0 and 3 are shown in Figure S2.  \n(a) Seg1 N Seg0  \n<!-- image -->  \n(d) Seg4 N Seg0  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \nFigure S1: CCF as a function of time delay for the two light curves between the lowest energy band (Seg0) and any of the other nine high energy bands (Seg1-Seg9). The orange line represents the best Gaussian fitting result.  \n<!-- image -->  \n(b) Seg2 N Seg0  \n<!-- image -->  \n(e) Seg5 N Seg0  \n<!-- image -->  \n(c) Seg3 N Seg0  \n<!-- image -->  \nCCF  \n(f) Seg6 N Seg0  \n<!-- image -->  \n<!-- image -->  \nFigure S2: Background count rate as a function of time and its polynomial fitting (red line). The upper and lower panel shows the fitting results for N hit segment 0 and 3, respectively.  \n<!-- image -->'}