Constraints on Lorentz Invariance Violation from GRB 221009A Using the DisCan Method

Lorentz symmetry represents a fundamental principle in modern physics, serving as a cornerstone of both special relativity and quantum field theory. However, despite the high-precision testing of symmetry across various contexts, there exists significant motivation to explore potential violations. A key motivation is the quest to resolve a fundamental challenge in modern physics: the reconciliation of General Relativity (GR) with Quantum Field Theory (QFT). While GR and QFT have achieved notable empirical success, they leave crucial questions unanswered, such as the emergence of spacetime singularities in GR. Addressing these issues necessitates a comprehensive quantum theory of gravity. Several Quantum Gravity (QG) models suggest the potential violation of Lorentz symmetry at scales where QG effects become significant, as evidenced by various studies [1–8]. This provides a compelling rationale for investigating Lorentz Invariance Violation (LIV) (for an overview, refer to [9]).

Recent advancements in high-energy astrophysical observations provide extra motivations of testing LIV beyond theoretical conjectures. Notably, the subtle spectral deviations between the Greisen-Zatsepin-Kuzmin (GZK) suppression features and theoretical predictions [10–12], alongside the anomalous attenuation deficit of TeV gamma rays from Active Galactic Nuclei (AGN), termed the TeV gamma-ray crisis [13], serve as pivotal observational probes due to their acute sensitivity to potential LIV [14]. These astrophysical phenomena serve as pivotal avenues for detecting potential violations of Lorentz symmetry. They offer essential insights into exploring deviations from standard dispersion relations in high-energy contexts. Pioneering work by researchers like collaborations such as MAGIC [15–19] have established increasingly stringent limits on potential Lorentz symmetry violations, particularly in high-energy astrophysical contexts. Our work is inspired by the MAGIC collaboration's efforts [15], aiming to verify whether similar constraints can be obtained in different source scenarios.

The DisCan (Dispersion Cancellation) method [20] emerges as a powerful technique for investigating LIV, does not require a pre-specified light curve model, in contrast to the maximum likelihood method. Utilizing the comprehensive photon dataset from the Large High Altitude Air Shower Observatory (LHAASO) [21], specifically the KM2A (Kilometer Square Array) and WCDA (Water Cherenkov Detector Array) detector. KM2A is a ground-based particle detector array covering approximately one square kilometer, composed of electromagnetic particle detectors and muon detectors. Its primary strength lies in high-sensitivity detection of gamma rays and cosmic rays in the energy range from tens of TeV to PeV (peta-electronvolts), enabling effective discrimination between cosmic-ray backgrounds and gamma-ray signals. WCDA comprises three large water pools with a total area of 78,000 square meters, utilizing Cherenkov light radiation to cover the 100 GeV to tens of TeV energy range. Key advantages of WCDA include its wide field of view, nanosecond-level timing resolution, and all-weather monitoring capability, making it particularly adept at capturing rapid variability in high-energy astrophysical phenomena. This method offers a novel approach to constraining potential Lorentz symmetry breaking. The motivation stems from the need to probe fundamental physics at energy scales approaching quantum gravity thresholds.

In this paper we focuse on analyzing the time-of-flight measurements of photons from GRB 221009A [22], one of the brightest gamma-ray bursts ever detected. By applying the DisCan method in conjunction with various information entropies as cost functions, and utilizing the Knuth binning method [23] as an optimization of the original fixed equal-bin distribution, aiming to provide improved constraints on potential LIV. Through the comparison of different information entropies, such as Shannon entropy, Rényi entropy, and Tsallis entropy, we derive effective limits on LIV that are consistent with previous NASA findings [20], indicating that Shannon entropy as a cost function yields higher credibility. More significantly, we extrapolate this approach to long-duration bursts and high-energy regimes, not just with simulated data, thus expanding the application domain of the DisCan method.

The structure of this paper is as follows : Section II presents the theoretical background and methodological approach. Section III details the DisCan method and its implementation. Section IV discusses the results and their implications. Finally, Section V offers our discussions conclusions

Quantum gravity theories propose that space-time exhibits quantum fluctuations, which may lead to a non-trivial refractive index that influences the propagation of particles [24]. Different quantum gravity models suggest that the speed of particle propagation can vary with energy, a phenomenon referred to as LIV.

A fundamental expectation associated with quantum gravity is a modified photon dispersion relation [25] :

where $ E_{Pl} \approx 1.22 \times 10^{19} $ GeV denotes the Planck energy. The parameter $ \varepsilon_{QG,n}^{(\sigma)} $ quantifies the LIV energy scale relative to $ E_{Pl} $. Typically, the value of n takes 1 or 2, corresponding to first-order or second-order breaking energy scales, respectively. The index $ \sigma = \pm1 $ indicates the direction of dispersion modification.

This dispersion relation introduces an energy-dependent photon velocity. In the regime where $ E \ll \varepsilon_{QG,n}E_{Pl} $, we can express the relationship as :

Photon propagation can exhibit either subluminal or superluminal characteristics, which correspond to $ \sigma = -1 $ and $ \sigma = +1 $, respectively. This energy-dependent velocity results in a time of flight (TOF) for a photon with energy E that differs from what would be expected if it traveled at the conventional speed of light [26].

The modified TOF can be mathematically represented as [26] :

where z denotes the redshift of the source. The Hubble-Lemaître constant, denoted as $ H_0 $, is taken to be approximately 67.5 km s^-1 Mpc^-1. The current fraction of matter in the Universe is represented by $ \Omega_M = 0.315 $, while the cosmological constant fraction is $ \Omega_\Lambda = 0.685 $ [27]. Specifically, for the case of GRB 221009A, which has a redshift of $ z = 0.151 $ [28], the relevant values are utilized in the expression.

In this work, which also makes use of the publicly accessible data from LHAASO for an initial exploration, we took a different approach from Piran's phenomenological average correction constant. Specifically, we utilized the complete time-energy information data of photons detected by KM2A [28] during GRB 221009A and applied a DisCan method to impose observational constraints on the potential Lorentz violation energy scale. The results obtained from our analysis mutually validate the findings of Piran's work.

The DisCan method is an algorithm for detecting quantum gravity dispersion in photons. The algorithm adjusts the arrival times of photons to cancel out the dispersion caused by quantum gravity effects, thereby detecting the energy-dependent time delays. Assuming the energy-dependent time delay is given by [20] :

where $ t_{\text{obs}} $ is the observed arrival time, and $ t_{\text{true}} $ denotes the time at which the photon would have arrived had it suffered no quantum gravity time shift [20]. $ \theta_1 $ is the dispersion coefficient when $ n=1 $ and E is the photon energy, with the time unit in s and the dispersion coefficient unit in s/TeV. Similarly, we can also assume a quadratic energy-dependent time delay as :

where $ \theta_2 $ is the dispersion coefficient when n=2, with the dispersion coefficient unit in s/TeV².

For both the first-order and second-order delay cases, we construct a time profile using the adjusted arrival times. Then, we optimize the sharpness of the time profile to find the best dispersion coefficient $ \theta_i $, using Shannon entropy, Rényi entropy and tsallis entropy.

Based on the previous formulas (3) (4) (5) and numerical values, and $ \delta t_{LIV}(E) = t_{\text{obs}} - t_{\text{true}}, $ we can derive :

Since $ \varepsilon_{QG} $ is always positive, we can determine σ through the parameter $ \theta_i $, thereby determining whether it is superluminal or subluminal.

The advantage of the DisCan method is its applicability to arbitrary dispersion models and its high sensitivity. By adjusting the photon arrival times, the DisCan method can effectively detect the small time delays caused by quantum gravity effects. We investigate time-of-flight constraints on LIV derived from KM2A’s observations of TeV photons from GRB 221009A, the brightest gamma-ray burst, and appropriately incorporated the low-energy photon data observed by WCDA.

The DisCan method is applicable to any dispersion model, as well as capable of detecting small time delays caused by quantum gravity effects. For this paper, after constructing the linear dispersion model, we first choose a cost function for the overview of time after binning, in order to optimize the selection of parameter $ \theta_1 $, $ \theta_2 $ and obtain the corresponding LIV parameters $ \varepsilon_{QG,1}^{(\sigma)} $, $ \varepsilon_{QG,2}^{(\sigma)} $. Finally, the breaking energy scales for $ E_{\text{QG},1} $ and $ E_{\text{QG},2} $ are obtained.

The LIV parameters $ \theta_i $ ($ i=1,2 $) are estimated by maximizing the information content of the reconstructed time profiles using transformed photon arrival times through dispersion relations (4) and (5). We employ three information-theoretic measures as cost functions:

● Shannon entropy

The Shannon entropy is defined [29] :

where $ p_n \equiv x_n/\sum x_n $ represents the normalized photon count probability [20], with $ x_n $ denoting the photon intensity in temporal bins. As the foundational measure of information uncertainty, it achieves maximal sensitivity for peaked distributions. Subsequent entropys also continue to employ these two concepts.

● Rényi entropy

The Rényi entropy is given [30] :

This generalized entropy reduces to Shannon entropy when $ \alpha \to 1 $, with $ \alpha > 1 $ emphasizing high-probability bins and $ \alpha < 1 $ enhancing sensitivity to rare events.

● Tsallis entropy

The expression for Tsallis entropy is as follows [31] :

It is characterized by non-extensive parameter $ q \neq 1 $, and has the same regulatory effect as α. Originally developed for non-equilibrium thermodynamics, this measure demonstrates superior performance in systems with long-range correlations.

When entropy is minimized, the information at the source is most certain, and it is most likely to obtain the true time profile, thereby obtaining the broken parameter.

We utilize the publicly available KM2A observational data from the LHAASO collaboration [28] to investigate the values of the parameters $ \theta_i $ obtained under first and second order linear models, thereby deriving the corresponding energy scale $ E_{\text{QG},n} $ for LIV. Initially, we utilized only the photon time-energy data from KM2A and referenced the 10-second binning provided by LHAASO to maintain the original equal-bin distribution [28]. Subsequently, we incorporated additional photon data from WCDA and replaced the original equal-bin distribution with the Knuth binning method. The Knuth binning method is a type of equal-bin method that employs a piecewise constant density model with equal bin widths. It utilizes the maximum Bayesian posterior probability to directly quantify the likelihood of binning. Its primary advantage is its adaptability to any distribution, including multi-peaked, asymmetric, and complex structures, enabling it to effectively capture complex temporal structures. Notably, this method requires little to no prior information [23]. Our approach to handling WCDA photon data also referenced the methodology of the LHAASO collaboration [22].

Since WCDA photons are initially assumed to possess identical energy, with any subsequent variations in energy attributable solely to statistical fluctuations, it is posited that these photons do not enhance the clarity of the source reconstruction. Instead, they substantially diminish the influence of KM2A photons. Consequently, it is deemed appropriate to subtract the information entropy generated by WCDA photons following the binning process. This leads to a natural adoption of a refined approach, wherein the information entropy, originally defined by the combined contributions of KM2A and WCDA photons, is adjusted by subtracting the information entropy derived exclusively from WCDA photons. This adjustment effectively replaces the traditional notion of information entropy with the concept of relative information entropy.

To account for photon energy uncertainties, we apply a Monte Carlo-based blurring technique to the photon energies, assuming Gaussian errors. In our Monte Carlo simulations ($ n_{\rm{sim}} = 1000 $), we model energy measurement uncertainties using asymmetric Gaussian distributions. For each data point i, positive errors follow $ \mathcal{N}(\mu = 0, \sigma = \sigma_{+i}) $ and negative errors follow $ \mathcal{N}(\mu = 0, \sigma = \sigma_{-i}) $, where $ \sigma_{+i} $ and $ \sigma_{-i} $ are the observed upper and lower $ 1\sigma $ error bounds respectively. For each realization, we generate random variates $ z \sim \mathcal{N}(0,1) $ and compute the error term as:

Synthetic energy samples are then generated through the transformation:

We construct the histogram distribution of the parameter $ \theta_i $. The 95% confidence interval is indicated by a dashed line in the plot.

Figs. 1 and 2 show the histogram distributions of the parameters $ \theta_1 $ and $ \theta_2 $ obtained from 1000 Monte Carlo samples, with dashed lines indicating the 95% confidence intervals. It can be observed that, when using Shannon entropy as the cost function, the resulting parameters $ \theta_i $ exhibit shorter confidence intervals and higher concentration in the same order, without the introduction of additional parameters. Therefore, we recommend the results obtained from Shannon entropy as the summary. Nevertheless, the results obtained using Rényi and Tsallis entropy remain significant. They offer insights that complement findings from other entropy measures, such as Shannon entropy. Also, as both our linear energy-dependent time delay formulas Eq. (4) and (5) and the use of information entropy as a cost function ignore some deeper physical considerations, such as the omission of extragalactic background light (EBL) absorption in the propagation model, source-intrinsic energy-dependent photon emission delays, and other unmodeled systematic uncertainties. Nevertheless, they still hold significance as observable tests for quantum gravity theories and comparative studies with other Lorentz symmetry violations.

Figure 1.  Probability density distributions of the first-order Lorentz invariance violation parameter $ \theta_1 $ estimated through: (a) Shannon entropy (baseline case), (b) Rényi entropy with $ \alpha=2 $, and (c) Tsallis entropy with $ q=2 $. Vertical dashed lines indicate 95% confidence intervals. All histograms contain 1,000 Monte Carlo realizations with identical color mapping and axis scaling.

Figure 2.  Probability density distributions of the second-order Lorentz invariance violation parameter $ \theta_2 $ estimated through: (a) Shannon entropy (baseline case), (b) Rényi entropy with $ \alpha=2 $, and (c) Tsallis entropy with $ q=2 $. Vertical dashed lines indicate 95% confidence intervals. All histograms contain 1,000 Monte Carlo realizations with identical color mapping and axis scaling.

Figs. 3, 4, and 5 present the histogram of the $ \theta_i $ distribution obtained after incorporating WCDA photons and applying the Knuth binning method.It can be observed that, upon employing the Knuth binning method, we directly exclude the scenario where relative Rényi entropy is used as the cost function. Additionally, a comparison is conducted for the parameter (q) in relative Tsallis entropy, revealing that the case with $ q=0.5 $ is superior to that with $ q=2 $. The physical interpretation of this finding is that our system exhibits significant non-extensivity, meaning that the total entropy of the system is not merely the sum of the entropies of its subsystems---a characteristic notably associated with Rényi entropy. In the analysis of gamma-ray burst (GRB) photon time series, the arrival times of photons may exhibit long-range correlations influenced by quantum gravity effects. Tsallis entropy with $ q<1 $ is more adept at capturing such correlations. Moreover, our adoption of relative entropy analysis reveals that the current task necessitates a focus on low-probability, high-value tail information—specifically, the KM2A photon data. This implies that the system may involve complex interactions that conventional extensive entropy measures cannot fully characterize. However, it is critical to note that the generality of this feature remains unverified. Future work should employ the DisCan method to conduct similar analyses on additional gamma-ray burst sources to validate these findings.

Figure 3.  The distribution histograms of parameters $ \theta_1 $ and $ \theta_2 $ derived from the application of Knuth's optimal binning algorithm with relative Shannon entropy as the cost function

Figure 4.  The distribution histograms of parameters $ \theta_1 $ and $ \theta_2 $ derived from the application of Knuth's optimal binning algorithm with relative Tsallis entropy ($ q=2 $) as the cost function

Figure 5.  The distribution histograms of parameters $ \theta_1 $ and $ \theta_2 $ derived from the application of Knuth's optimal binning algorithm with relative Tsallis entropy ($ q=0.5 $) as the cost function

Our initial motivation for introducing Tsallis entropy was its non-extensivity, which naturally describes non-local propagation effects and, when the dispersion relation is non-linear, the non-extensive framework of Tsallis entropy is more easily extended to multi-parameter optimization. This inspires us to consider modifying the dispersion relation and incorporating more quantum gravity spacetime characteristics into the model in future research.

Tables 1, 2, and 3 list the 95% confidence intervals for the parameters $ \theta_i $, $ \varepsilon_{{\rm{QG}},n} $, and $ E_{{\rm{QG}},n} $, where $ \sigma = +1 $ represents the superluminal case and $ \sigma = -1 $ represents the subluminal case. It can be observed that for the parameter θ, regardless of whether it is first-order or second-order, and regardless of the entropy measure used as the cost function, the absolute value of the lower bound is always greater than that of the upper bound. We have obtained conclusions consistent with those of the LHAASO collaboration group [22]. This indicates that the distribution is not symmetric, but rather biased. An additional feature observed in Figures 3–5 is that, after incorporating WCDA photons (which exhibit lower energy contrast compared to KM2A photons) and applying the Knuth binning method, the photon distribution shows a preference for subluminal propagation. This trend contrasts with the results obtained from KM2A data alone.

	$ \theta_1 $ (s/TeV)			$ \theta_2 $ (s/TeV2)
	Lower	Upper		Lower	Upper
$ H_{\rm{S}} $	−2.545	1.290		−0.380	0.090
$ H_{\rm{R}}(\alpha=2) $	−5.435	2.380		−0.380	0.240
$ H_{\rm{T}}(q=2) $	−6.007	2.770		−0.450	0.270

Table 1.  95% CIs for $ \theta_1 $ and $ \theta_2 $

	$ \varepsilon_{QG,1} $			$ \varepsilon_{QG,2} $ (10−7)
	$ \sigma=+1 $	$ \sigma=-1 $		$ \sigma=+1 $	$ \sigma=-1 $
$ H_{\rm{S}} $	2.2	4.4		2.0	8.3
$ H_{\rm{R}}(\alpha=2) $	1.0	2.4		2.0	3.1
$ H_{\rm{T}}(q=2) $	0.95	2.1		1.7	2.8

Table 2.  95% CIs for $ \varepsilon_{QG,1} $ and $ \varepsilon_{QG,2} $

	$ E_{\text{QG},1} $ (1019 GeV)			$ E_{QG,2} $ (1012 GeV)
	$ \sigma=+1 $	$ \sigma=-1 $		$ \sigma=+1 $	$ \sigma=-1 $
$ H_{\rm{S}} $	2.7	5.4		2.4	10.0
$ H_{\rm{R}}(\alpha=2) $	1.3	2.9		2.4	3.8
$ H_{\rm{T}}(q=2) $	1.2	2.5		2.0	3.4

Table 3.  95% CIs for $ E_{\text{QG},1} $¹⁾ and $ E_{\text{QG},2} $

In Table 4, it can be observed that for the case of relative Tsallis entropy with $ q=2 $, in the second-order scenario, both lower and upper limits for the second-order Lorentz violation energy scale are obtained. When employing the Knuth binning method and relative $ H_{\rm{T}}(q=2) $ entropy, the following results are derived: in the superluminal case, the lower limit for $ E_{\text{QG},1} $ is $ 2.0 \times 10^{19} \, \text{GeV} $; in the subluminal case, the lower limit for $ E_{\text{QG},1} $ is $ 1.2 \times 10^{19} \, \text{GeV} $; additionally, in the subluminal case, the lower limit for $ E_{\text{QG},2} $ is $ 5.2 \times 10^{12} \, \text{GeV} $ and the upper limit is $ 6.1 \times 10^{13} \, \text{GeV} $. These limits are determined based on the 95% confidence interval.The results in Tables 4 and 5 demonstrate that we have further considered the influence of the afterglow of GRB 221009A and the effect of parameter adjustment in Tsallis entropy.

	$ \theta_1 $ (s/TeV)			$ \theta_2 $ (s/TeV2)
	Lower	Upper		Lower	Upper
$ H_{\rm{S}} $	−7.476	1.586		−0.130	0.135
$ H_{\rm{T}}(q=2) $	−3.556	5.577		0.015	0.175
$ H_{\rm{T}}(q=0.5) $	−10.676	1.080		−0.330	0.110

Table 4.  95% CIs for $ \theta_1 $ and $ \theta_2 $ with Knuth's binning method and relative entropy

	$ E_{\text{QG},1} $ (1019 GeV)			$ E_{QG,2} $ (1012 GeV)
	$ \sigma=+1 $	$ \sigma=-1 $		$ \sigma=+1 $	$ \sigma=-1 $
$ H_{\rm{S}} $	0.93	4.4		7.0	6.8
$ H_{\rm{T}}(q=0.5) $	0.65	6.4		2.8	8.3

Table 5.  95% CIs for $ E_{\text{QG},1} $ and $ E_{\text{QG},2} $ with Knuth's binning method and relative entropy

In the process of deriving $ E_{{\rm{QG}},n} $ from $ \varepsilon_{{\rm{QG}},n} $, we find that, regardless of whether the first-order or second-order linear delay model is used, the subluminal constraint is always stronger than the superluminal one. Specifically, the first-order energy scale is on the order of $ 10^{19} $ GeV, and the second-order energy scale is on the order of $ 10^{12} $ GeV. The constraints obtained using Shannon entropy as the cost function are stricter, indicating that it provides more comprehensive results in the absence of additional physical considerations.

In this paper, we employed the DisCan method with different information entropy measures to investigate the constraints on LIV parameters using the LHAASO's KM2A observational data of GRB 221009A. Similar constraints were obtained using different entropy measures. Specifically, using the DisCan method with Shannon entropy, the 95% confidence level (CL) lower limits for the subluminal (superluminal) scenarios were determined to be $ E_{\text{QG},1} > 5.4 (2.7) \times 10^{19} \, \text{GeV} $ for the linear case ($ n = 1 $), and $ E_{\text{QG},2} > 10.0 (2.4) \times 10^{12} \, \text{GeV} $ for the quadratic case ($ n = 2 $). By using the alternative Rényi entropy with $ \alpha = 2 $, the lower limits were found to be $ E_{\text{QG},1} > 2.9 (1.3) \times 10^{19} \, \text{GeV} $ and $ E_{\text{QG},2} > 3.8 (2.4) \times 10^{12} \, \text{GeV} $, while Tsallis entropy with $ q = 2 $ produced similar results to Rényi entropy. The Lorentz-violating energy scale obtained in our study is approximately of the same order of magnitude as that derived from the recent work by LHAASO [32]. These results indicate that while the exact choice of entropy measure affects the precise values, the order of magnitude of the LIV energy scale remains consistent. As an optimization and supplement to the initial methods, we introduced WCDA photons and the Knuth binning method to enhance the original approach. This not only makes the results more diverse and engaging but also highlights the differences when various information entropies are employed as cost functions. Additionally, it underscores the further screening effect of the Knuth binning method on information entropy. Our findings indicate that Rényi entropy, being an additive entropy, is less suitable than Tsallis entropy for systems characterized by long-range interactions. Consequently, this insight directs our future research efforts toward a more in-depth investigation of the parameter (q). It should be emphasized that, at present, our findings do not definitively suggest the existence of long-range interactions in the observed gamma-ray burst spectra. To obtain more definitive signals, it is necessary to conduct similar analyses on other events, such as GRB 130427A and GRB 190114C. This aspect of the research will be undertaken in future work.

The derived constraints for the LIV parameters $ \theta_1 $ and $ \theta_2 $, as well as the corresponding energy scales $ E_{\text{QG},1} $ and $ E_{\text{QG},2} $, reveal a clear preference for subluminal scenarios, with stronger limits in the subluminal regime. Although we only used the KM2A data and lacked more detailed data from the LHAASO collaboration, this conclusion is consistent with theirs [22]. Moreover, the energy scales corresponding to the first-order LIV ($ E_{\text{QG},1} $) are on the order of $ 10^{19} \, \text{GeV} $, while the second-order LIV scales ($ E_{\text{QG},2} $) are on the order of $ 10^{12} \, \text{GeV} $. The use of Shannon entropy as a cost function yields reliable constraints, highlighting its advantage in providing more reliable results when additional physical factors are not incorporated into the model.

Based on Table 1 of [35], the lower limits on the energy scale of LIV are initially constrained by subluminal propagation scenarios. However, as photon energies increase to the TeV range, a significant number of lower limits are derived from superluminal propagation cases. Regarding the gamma-ray burst GRB 221009A, our analysis results are consistent with the conclusions of the LHAASO collaboration [22]. The underlying mechanism responsible for superluminal propagation, however, remains unclear, this will be the objective of our next research phase.

The LHAASO collaboration elucidates that the rapid rising and gradual decaying phases of the GRB afterglow's light curve are particularly sensitive to subluminal LIV effects [22]. This asymmetry arises because the time-delay signatures induced by subluminal dispersion (σ = -1) align more prominently with the steep flux variations in the early afterglow phase, thereby yielding tighter constraints compared to superluminal scenarios (σ = +1) where dispersion effects may be diluted by the smoother temporal evolution. Due to the rapid rise and slow decay characteristics of the afterglow, high-energy photons are predominantly concentrated in the early time periods, resulting in a higher accumulation of photons within individual bins. Our method appears to favor an even distribution of photons across each bin, thus preferring a dispersion coefficient θ that shifts the densely packed high-energy photons on the left toward the right. This likely explains why our constraints on subluminal LIV effects are stricter than those on superluminal LIV effects.

Another possible issue is that we neglect the potential effect of Extragalactic Background Light (EBL) in the present work. However, as found by the LHAASO collaboration in [21], EBL absorption in the present case exerts minimal influence on the light curves of both high-energy and low-energy photons, rendering our methodology largely insensitive to EBL effects.

It is noteworthy that our methodology is grounded in photon time series counting analysis, which differs from statistical methods that rely on parametric fitting. A promising avenue for future research lies in integrating these two approaches. Significantly, existing literature has accounted for photon emission time effects by employing data from multiple GRBs at different redshifts to conduct statistical joint analyses [36–38], thereby mitigating source-intrinsic effects. It is of great interest to combine our method to perform such multiple GRB analysis in the future, so as to effectively reduce the uncertainty of source effects.

In addition, to validate whether our statistical procedure is robust or not, it is of great significance to establish that the confidence intervals (CIs) produced by the method reliably meet the desired coverage criteria, as illustrated in [39]. This is an another important direction of future work.

In conclusion, our study provides stringent constraints on LIV using the LHAASO observations of GRB 221009A. The results show that the subluminal scenario offers stronger constraints, and Shannon entropy provides the most reliable estimates. However, further refinement of these constraints will require consideration of additional physical factors, such as EBL absorption and intrinsic energy-dependent emission time effects, as well as the exploration of other quantum gravity models and dispersion relations in future studies.

We extend our gratitude to Yuming Yang for useful discussions, and particularly thank Xiao-Jun Bi for his insightful discussions and valuable comments.