Robust Constraint on Lorentz Violation Using Fermi-LAT Gamma-Ray Burst Data

Models of quantum gravity suggest that the vacuum should be regarded as a medium with quantum structure that may have non-trivial effects on photon propagation, including the violation of Lorentz invariance. Fermi Large Area Telescope (LAT) observations of gamma-ray bursts (GRBs) are sensitive probes of Lorentz invariance, via studies of energy-dependent timing shifts in their rapidly-varying photon emissions. In this paper we analyze the Fermi-LAT measurements of high-energy gamma rays from GRBs with known redshifts, allowing for the possibility of energy-dependent variations in emission times at the sources as well as a possible non-trivial refractive index in vacuo for photons. We use statistical estimators based on the irregularity, kurtosis and skewness of bursts that are relatively bright in the 100 MeV to multi-GeV energy band to constrain possible dispersion effects during propagation. We find that the energy scale characterizing a linear energy dependence of the refractive index should exceed a few $\times 10^{17}$ GeV, and we estimate the sensitivity attainable with additional future sources to be detected by Fermi-LAT.


Introduction and Summary
The idea that the space-time vacuum should be regarded as a non-trivial medium -baptized 'space-time foam' by Wheeler [1] -is based on very general intuition. This intuition arises from the feature of quantum mechanics that on time scales ∆t any physical system must exhibit virtual energy fluctuations ∆E with magnitudes ∆E ∼ /∆t. Wheeler [1] argued that on time scales ∆t ∼ 1/M P , where M P ∼ 10 19 GeV is the Planck mass: M P = c/G N and G N is the Newton constant of classical gravity, there would appear quantum-gravitational fluctuations in the space-time continuum with ∆E ∼ M P , resulting in a 'foamy' structure on short time scales ∆t ∼ /M P . This observation led Wheeler to argue that space-time would no longer appear smooth at distance scales ∆x ∼ /M P , and that it might exhibit both non-topological irregularities and topological fluctuations. This intuitive picture suggests the appearance of a refractive index η for particles such as photons propagating through 'empty' space, corresponding to a phase velocity v ph = p/E = c/η [2]. It can be argued on general grounds that photons should not travel faster than c 1 , because otherwise they would emit gravitationalČerenkov radiation and lose energy unacceptably quickly [3] 2 . Hence, the photon refractive index η ≥ 1, corresponding to subluminal propagation of energetic photons, as predicted in simple models [2,[5][6][7]. One would also expect that the refractive index should increase with energy, because gravitational interactions are proportional to some negative power of M P , in general, and therefore should increase as some positive power of the energy 3 . Models [2,[5][6][7] suggest that the photon group velocity might deviate from that of light linearly in photon energy E: where one might expect that M 1 = O(M P ). However, the Lorentz-violating (LV) scale M 1 would depend on unknown parameters of the microscopic theory, including the string scale, which = M P , in general. In the D-foam model discussed below [2], M 1 would depend also on couplings to D-particles, which would depend on the particle species, and on the local density of D-particles 4 . Moreover, other energy dependences: η −1 ∼ (E/M n ) n should also be considered, such as the case n = 2 [9,10].
The proposal that there might be observable effects on the propagation of particles such as photons [11] was made originally in the context of concrete models of space-time foam motivated by (non-critical) string/brane theories [2,9]. These models go beyond conventional local quantum field theories, and contain the necessary ingredients for discussing the interaction between a propagating matter particle and a quantum-gravitational 'environment'. The former is described as an (open) string excitation representing some observable (Standard-Model-like) particle of matter or radiation, moving through a (3+1)-dimensional brane Universe [12]. The 'environment' of quantum-gravitational fluctuations is provided in this context by ensembles of quantum space-time defects described as D-particles [2], which move in the higher-dimensional bulk. They consist of branes that are compactified in such a way that, from the point of view of a low-energy four-dimensional observer living on the (3+1)-dimensional brane Universe, they look approximately point-like [8]. When D-particles cross this D-brane world, they are perceived in our Universe as space-time events localized at specific locations x and specific times t, which we call 'D-foam'.
There are non-trivial interactions between bosonic open strings and such D-particles, consisting of splitting of the open string and emission of other open string excitations stretching between the D-particle and the brane world [8]. These interactions must respect the gauge symmetries on the brane world, such as the U(1) of electromagnetism. The consequent charge conservation implies that 'D-foam' appears transparent to charged bosons, but not to neutral ones such as photons and gravitons. Moreover, in the absence of low-energy supersymmetry, fermions such as right-handed neutrinos would have suppressed interactions with the Dfoam [13]. Conventional left-handed neutrinos are doublets under the electroweak SU(2) group of the Standard Model, so their interactions with the D-particles are further suppressed.
We have studied the possible observable consequences for photons propagating through our D-brane Universe in several previous papers [2,[8][9][10][11]14]. In general, if a photon encounters a D-particle, its interaction with it may resemble that of a photon propagating through a transparent material medium such as glass. In that case it may interact with the electrons that it contains, via absorption and subsequent re-emission, with the net effect of slowing down the photon. Thus, light travelling through glass has a refractive index η > 1. Moreover, the value of η varies with the colour of the light, i.e., the energy of the associated photon. Similarly, we expect in general that light travelling through the quantum-gravitational vacuum would acquire an energy-dependent refractive index η > 1, that we may model via interactions with D-particles [8]. On the other hand, because of the absence of interactions between charged particles and D-foam the deviation of the refractive index of the electron from unity would be suppressed [8,13], as required phenomenologically.
Many other models of Lorentz violation have been proposed. These include purely phenomenological models motivated by aspects of cosmic-ray physics [15] and other considerations [16]. A more theoretical suggestion -the 'Standard Model Extension' -is the possibility that Lorentz invariance is broken spontaneously [6,17,18], and it has been argued in the context of some models of loop quantum gravity [5] that the vacuum might exhibit non-trivial optical properties. Moreover, quantum field theories of the Lifshitz type [19], in which the space and time coordinates scale differently, have attracted renewed interest in the context of quantum gravity. In Lifshitz theories Lorentz invariance can be violated at high energies, but is restored in the low-energy limit. Another approach is that of doubly (or deformed) special relativity [7], in which Lorentz invariance is fundamentally deformed, rather than violated 5 .
How can one probe such ideas, in particular the modification (1) of photon propagation in vacuo? It was suggested in [11] that variable astrophysical sources would provide the most sensitive probes of such Lorentz violation, in view of their large distances and and non-trivial time structures, as well as their emissions of high-energy photons. Examples of such sources that were suggested in [11] include pulsars, active galactic nuclei and gamma-ray bursts (GRBs).
The Fermi-LAT sample of the latter is the subject of the analysis in this paper 6 . There have been many previous analyses of such variable astrophysical emissions. The first systematic study of possible Lorentz violation using the light curves of a number GRBs with emissions in the sub-MeV energy range distributed over a range of red shifts was presented in [9] 7 . The study was extended subsequently in [10,14], incorporating the light curves of substantially larger samples of GRBs, applying advanced time series analysis techniques (wavelets), and scrutinizing the possible systematic uncertainties inherent to such kind of analyses. Studies exhibiting similar levels of sensitivity have also been performed elsewhere [28,29]. Recent searches for Lorentz violation using samples of sub-MeV light curves of short GRBs [30] detected by the Swift satellite have not reported an improvement of sensitivity, compared with the analysis of [10]. On the other hand, observations by Fermi-LAT [31] of high-energy emission from GRBs where the energies of some individual gamma rays exceeded 10 GeV made possible a substantial increase in the sensitivity to M 1 , approaching the Planck scale. For example, an 5 See [27] for a discussion of the definition of measurable momenta in such models. 6 In field-theoretical models such as the Standard Model Extension [6,17,18] in which Lorentz invariance is broken spontaneously, there are also birefringence effects. Probes of the rotation of the polarization of light from distant astrophysical sources [20,21] constrain this effect very strongly [22][23][24][25][26]. However, birefringence is not expected in the models of space-time foam studied here. 7 We recall, however, that the effective quantum-gravity scales would depend on the density of D-brane defects in D-foam models [8], which could vary with the cosmological epoch.
analysis of time differences in the arrival times of individual gamma rays from the single source GRB080916C suggested a limit on M 1 that was about 2 orders of magnitude [32] stronger than that in [10]. Another source GRB090510A detected by Fermi-LAT was used to give a trans-Planckian lower limit on this scale of quantum gravity [33][34][35]. Moreover, it was argued in [36] that the assumption of a particular "rhythm" in the arrival of multi-GeV gamma rays from GRB090510A could even push the lower limit on the quantum gravity scale up to two orders higher than the Planck scale. Another assumption was made in [37][38][39][40], where it was suggested that the source frame time offset of the individual highest-energy gamma rays in emissions of Fermi-LAT objects should coincide to very high precision with the time offset of the peak emission of the sub-MeV energy light curves of the objects. This analysis led to a claim of a signal for the violation of Lorentz invariance, rather than a lower bound, corresponding to a quantum-gravity scales M 1 ∼ 10 17 GeV.
However, reports of sensitivities to Lorentz-violating effects with M 1 ∼ M P and, a fortiori, claims of signals, are beset with systematic uncertainties associated with our ignorance of the energy dependence of the times at which photons are emitted at the source. In particular, the literature also contains considerable discussion of the possibility that some higher-energy photons may be emitted later than prompt lower-energy photons, see [41], powered by a relativistic blast wave in the circum-GRB medium. However, we do not enter this discussion here.
Instead, our aim is to develop statistical techniques that minimize the impact of such source effects, which is the central point of our paper 8 . In the current study, we consider three distinct statistical measures of GRB emissions that mitigate source effects, which we use in an analysis of Fermi-LAT data in an attempt to obtain the most robust constraints on Lorentz violation associated with modified dispersion relations of photons during their propagation in a quantum foamy space-time medium. Although our analysis is motivated by one particular framework for Lorentz violation, the statistical techniques developed and the results obtained here are applicable to a wide class of such models.
The structure of the article is as follows: in Section 2 we review the basic features of propagation of a pulse in a dispersive quantum-gravity medium, while in Section 3 we present the Fermi-LAT data to be used in our analysis. In Section 4, we discuss methods to recover properties of source timing that will play an important rôle in our attempts to extract robust constraints on Lorentz violation. In Section 5, we embark on our main part of the analysis, by describing the various statistical measures of GRB emissions that we use to mitigate source effects. The first, the Irregularity Estimator is based on the observation that dispersion due to propagation through the space-time medium would tend to 'dilute' any burst-like 8 Some attempts to combine such source effects with propagation effects due to a potential quantum gravity medium, based on a particular "magnetic-jet" model for GRB emission [42], can be found in [43]. feature, leading asymptotically to a distribution that shows no time-dependent features above the background. One may then constrain the Lorentz violation parameters M n by minimizing this dilution. The second statistical measure, the Kurtosis Estimator is based on the related observation that the kurtosis, namely the height of a distribution relative to its standard deviation, would also be reduced by the effects of propagation through the space-time medium.
Finally, the Skewness Estimator exploits the fact that an energy-dependent reduction in photon velocity would increase the skewness, or asymmetry, of burst-like features in the emissions. Uncertainties in these estimators are discussed in Section 6, and in Section 7 we apply these methods to the ensemble of Fermi-LAT data on emissions from GRBs with bright emissions in the 100 MeV to multi-GeV energy band. Our analysis leads to a lower limit on M 1 in the range 2.4 to 8.4 × 10 17 GeV, which we consider to be the most robust constraint to date on this type of Lorentz violation induced by a dispersive quantum-gravitational medium. A brief discussion of these results, the associated uncertainties and ways to improve them, is presented in Section 8. Finally conclusions and outlook are presented in Section 9.

Pulse Propagation in a Quantum-Gravity Medium
In this Section we review basic features of the deformation of the envelope of an electromagnetic wave packet during its propagation in a quantum-gravity dispersive medium [9] that leads to the refractive index effect (1).
The basic solution of the wave equation is a plane wave of the form where k is the momentum and ω the frequency, and the superposition principle leads to the general solution Conversely, the amplitude A(k) can be expressed as Here, for simplicity, we consider a normalized Gaussian wave packet, with variance a 2 : The amplitude in momentum space of such a wave is If the distribution A(k) is sharply peaked around some value k 0 , the group velocity of a traveling pulse is given by Provided that the quantum-gravity-induced deviation of the propagation velocity from the speed of light would then imply that the dispersion relation has the form for β n k n << 1, the group velocity is given v g ≈ 1 + (n + 1)β n k n 0 , which yields a correction of the form (1) if β 1 is negative and n = 1.
In the case of a Gaussian packet formed by superposing traveling plane waves of momentum k 0 that locates at x = 0 at t = 0 and propagates along the x direction, the corresponding Fourier amplitude is given by At a later time t > 0, the pulse (11) will evolve as which for n = 1 readily reduces to where ω 0 = k 0 (1 + β 1 k 0 ). Evaluating the integral in (13), one arrives at the following final expression for the amplitude of the wave packet  (8), convoluted at two different points in time t 2 > t 1 > 0 with a powerlaw energy spectrum (17). The (green) solid line represents the profile of the wave envelope at the injection time (t = 0). The (orange) dashed and (magenta) dash-dotted lines represent the profiles modified by the propagation effects at times t 1 and t 2 , respectively. The energy range spans two orders of magnitude in dimensionless units, so as to reproduce a typical spectral range of the high-energy GRB emissions measured by Fermi-LAT used in the analysis.
Eventually, one also obtains the intensity of the wave packet: It is easy to see from (15) that, as time evolves, the peak of the amplitude gets shifted to x + v g t and the amplitude of the envelope reduces. The packet becomes wider in such a way that an initially narrower packet spreads much faster compared to one that has a larger initial width a.
Let us now assume that a signal with Gaussian profile and some spectral content Φ(k) is located within a band spanning a certain range from k 1 to k 2 . In this case, since the group velocity in a dispersive medium (9) depends on the wave number k, the overall intensity should be calculated as the convolution The Fermi-LAT high-energy emission spectra of GRBs can be well approximated by a power-law model [44] with positive spectral index α 9 : and results of the numerical convolution (16) at two points in time, using to the profile (15) and the spectral model (17) are presented in Fig. 1.
As can be seen in Fig. 1, the burst-like feature of a GRB intensity profile modelled by a Gaussian envelope is deformed during propagation in a quantum-gravity dispersive medium of the type considered in this work. Three features of this deformation can be distinguished, which we exploit subsequently in our analysis.
(i) As the signal of the GRB propagates, irregularities in its intensity profile, superposed upon a background, get diluted, causing any pulsing intensity profile to approach an almost featureless the background like time profile at large times. One can invert this possible quantumgravity propagation effect, converting the timings and energies of photons arriving in highenergy emissions from distant GRBs detected by Fermi-LAT back to the intensity distribution that would have been injected at the source, compensating the signal timings for the propagation delays (see Section 3 for details). One can then estimate the amount of Lorentz violation affecting the propagation of the signal by calculating the amount of compensation that maximizes the irregularities (spikes) of the GRB intensity distribution injected at the source, denoted by t = 0 in Fig. 1. The qualitative and quantitative picture of smoothing of irregularities in the initial intensity distribution of a burst-like signal by Lorentz violation during its propagation holds for any shapes of initial spikes, which may be quite irregular, unlike the Gaussian example shown in Fig. 1.
(ii) One such effect on the intensity profile is that the heights of peaks in the signal are reduced during its propagation in a dispersive medium. Therefore, one can quantify the relative sizes of the peaks of the compensated intensity distribution and estimate the amount of Lorentz violation by calculating the amount of compensation that would maximize the peaking of the distribution at the time of emission. In making this estimate, we use in Section 5.2 a method that does not depend on the particular shapes of the peaks.
(iii) The intensity time profile becomes more asymmetric with time. The degree of asymmetry of a time profile can be estimated by comparing it with a symmetric reference distribution, which we take for convenience to be a normal distribution. One can use as an estimate of the effect of propagation through a quantum-gravitational medium the amount of Lorentz violation that maximizes the symmetry of the initial distribution. We note that the estimator we use in Section 5.3 to quantify this kind of deformation of the signal is not restricted by any assumption on the shape of the initial peaks in the intensity distribution.
In later Sections we present in detail the methods used to estimate the deformation features outlined above.

Fermi-LAT Data
The data analyzed in this work are taken from the Fermi-LAT Pass 8 transient event class P8R2 TRANSIENT010 [45]. The background contamination in this set of data is calibrated to the best-fit power-law parametrization of the Isotropic Diffuse Gamma-Ray Background (IGRB) emission from [46]. The loosest selection criteria for this TRANSIENT class are designed for shortduration events, such as GRBs, that benefit from increased photon statistics while tolerating a higher background fraction and the broader point spread function (PSF) of LAT. This class has a background rate that is equal to the IGRB reference spectrum and requires the presence of a signal in both the tracker and the calorimeter.
The required data set is extracted from the publicly available archive [47] at Fermi Science Support Center.
The Fermi mission provides a suite of tools, called the Fermi Science Tools, for the analysis of LAT data, and the tool to perform selection cuts is called gtselect. A lower energy limit of 100 MeV on photon energies is chosen to reject photons with poorly reconstructed energies and directions. No maximum energy cut is applied, since photon energies can reach a few tens of GeV. We select photons reconstructed within a circular region of interest (ROI) centered on the best available GRB localization with a radius corresponding to the 95% containment radius of the transient LAT PSF estimated for a 100 MeV photon. The GRB directions used to specify the center of the ROI are obtained by follow-up ground-based observations, and can be assumed for our purposes to coincide with the true direction of the GRB.
The data for the 8 GRBs with measured red shifts and relatively high numbers of photons with energies above 100 MeV detected by Fermi-LAT that are used in our analysis are presented in Table 1. Fig. 2 shows scatter plots of the photon energies versus arrival times for two GRBs in this data sample. The data for all the GRBs in Table 1 are processed similarly, using the various estimation procedures described below.

Recovery of the Source Timing Properties
If β n in (8) is set to where M QGn represents the scale at which Lorentz-violating quantum gravity effects set in, the group velocity (9) acquires the form where h(z) = Ω Λ + Ω M (1 + z) 3 and H 0 = 68 km/s/Mpc is the present Hubble expansion rate (see for example [48]). Thus, the difference in proper distance covered by two photons emitted at red shift z src with velocity difference ∆v g is It follows from (19) that the velocity difference of two photons of energies E 2 > E 1 is given by Therefore, the difference in the arrival times of these photons is where the factors in (23) are a QGn = n+1 In the following, the earliest arrival time of a photon from a given GRB is always set to zero 10 .
As discussed above, the expectation that photons with lower energies (longer wavelengths) may be delayed less than photons with higher energies (shorter wavelengths) implies that the temporal pattern of photons arriving from a given GRB should be modified compared to the pattern when emitted by the source. However, in order to elucidate the possible magnitude of the quantum gravity dispersion effect, we need statistical estimators that enable us to discriminate between source and propagation effects.
For the purpose of our analysis, we allow the arrival time of every detected photon to incorporate an a priori unknown source-related time-lag b sf as well as the energy-dependent time delay (23) accumulated in the course of propagation because of quantum gravity dispersion [10].
In the case of linear energy dependence, n = 1 (1), one finds the following expression for the arrival times of individual photons: where is a "compensation parameter" that quantifies the amount of the linearly energy-dependent propagation effect encoded in the signal from a given source. "Compensation" is understood here in the sense of recovery of the original pattern of the intrinsic emission times by removing a possible propagation effect related to quantum-gravity dispersion. In practice, instead of the source frame intrinsic timings b sf (E i ) we use the detector frame intrinsic timings given by The correct value of the compensation parameter applied in (26) recovers the intrinsic pattern of the timings b df (E i ) before being dispersed by quantum-gravity effects. Because of unknown details of the source activity due to our imperfect knowledge of the radiation mechanism of GRBs, as well as of potential stochasticity during the burst, the intrinsic source distribution b df (E i ) is expected to be different for different GRBs. This is the main challenge for inferring a common propagation effect from samples of high-energy gamma rays emitted by different GRBs.
As we have demonstrated in Section 2, a Gaussian emission envelope would be deformed during its propagation through a dispersive quantum-gravity medium. We may assume that there would be similar deformations in the shape of an emission envelope of arbitrary profile with burst-like features in its temporal intensity distribution. Therefore, one may estimate the compensation parameter using statistical quantifications of the deformations in the intensity profile of an envelope propagating in such a medium.
For a given source, the data are represented by N points, each one associated with the arrival time and energy of a photon reconstructed by the Fermi-LAT (for two examples of sources used in our analysis, see Fig. 2). Describing the pattern of intrinsic timings of individual photons in the detector frame by a probability distribution function F(b df (E i , τ )), the shape of the intensity distribution is given by where is the energy weight of a given photon with energy E i normalized by the energy E min of the softest photon measured by the Fermi-LAT within the sample from a given GRB. The compensation parameter τ in (27) is defined in such a way that the profile I(b df (E i , τ = 0)) coincides with that measured in the detector, after the deformation by propagation effects.
This deformation can be characterized by non-parametric estimators whose optimization, using appropriate criteria, allows one to estimate the correct value of the compensation parameter. In practice, we calculate the estimators for trial values of the compensation parameter τ j applied in (26), where j indicates one of a set of values taken either from a regular grid or random values distributed uniformly over a pre-defined range of values of τ . The correct value of τ j should generate an intensity distribution I(b df (E i , τ j )) that satisfies in the best possible way the criteria for recovery of the genuine pattern of the detector-frame intrinsic timing for a given estimation technique, as we discuss in the next Section.

Estimation Techniques
We present three estimators in this Section, one sensitive to each type of deformation of a wave envelope with burst-like features described in Section 2. They are used subsequently in our analysis of the Fermi-LAT data presented in Section 3. As we demonstrate below, these techniques lead to robust constraints on the photon refractive index potentially induced by Lorentz-violating quantum-gravity effects.

Irregularity Estimator
We have shown in Section 2 that a burst-like signal with a power-law spectrum, as modelled by a Gaussian shape, gets "diluted" with time as it propagates in a dispersive medium. In general, a qualitatively similar result is expected for a signal possessing any kind of burst-like activity superposed on a "regular" background distribution. Thus we expect that during dispersive propagation any irregular signal with burst-like features degenerates in shape, approaching this background distribution 11 . Conversely, application of the procedure for compensating quantum gravity propagation effects described above should recover the intrinsic irregularities at the source. The degree of irregularity can be estimated by comparing a compensated intensity distribution with a reference one, the latter being an a priori featureless (smooth) distribution with the same statistical strength. It is clear that , in the absence of any insight into the physics of the non-variable part of the high-energy emission by the GRB engine and hence any assumptions about the shape of the background, the best reference featureless distribution is the uniform one. Since the shape of the intensity PDF at the source is unknown, we make the comparison on the basis of non-parametric statistics. This ensures that the analysis is as independent as possible of assumptions on the shape of the irregularities an hence the dynamics of the variability of of GRB engines in the high-energy band.
We use the Kolmogorov-Smirnov (KS) statistic (see, for example, ref. [49]) to estimate the degree of irregularity of an intensity distribution. We define the KS difference between two distribution functions Ξ T (t) and U T (t) as where t 0 and t N −1 , which themselves are functions of τ , represent the timings of the first and the last photon within a compensated distribution, respectively.
Following (27), the Fermi-LAT list of photons is converted as follows into the distribution function Ξ T (τ, t). First, for a given source, we associate every photon with its energy weight (28). Then, the function Ξ T (τ, t) is constructed as the fraction of those energy weights Examples of the D(τ ) distributions for two GRBs (GRB090510A and GRB080916C) are presented in Fig. 3. The values of τ at the maxima are the best estimates of the values of the compensation parameter that recover the initial irregularities of the intensity distribution at the source. Since the D(τ ) curves in Fig. 3 exhibit significant fluctuations, in particular around the peaks, we use three methods to localize the maxima of the D(τ ) distributions, which are described in more detail in Section 6 below. As an example, the positions of the maxima shown in Fig. 3 were obtained by averaging the distribution of τ γ values for which D(τ γ ) exceeds 95% of its peak value. The corresponding detector-frame intrinsic times derived from the positions of the maxima of the D(τ ) curves are plotted in Fig. 2 as open triangles, which can be compared with the detected arrival times (green squares).

Kurtosis Estimator
As already discussed in Section 2, as a signal becomes more "diluted" during its propagation in a dispersive medium, the peaks in its burst-like features degrade. Therefore, the relative sizes of the peaks can serve as another measure of signal deformation by quantum gravity effects. To quantify this effect on the intensity distribution we use a measure of kurtosis, which provides information on the height of the peak of a distribution relative to the value of its standard deviation. For the intensity defined in (27), the kurtosis formula for a compensated distribution is where, as above, every photon is associated with a weight W i given by (28), while b df (τ ) stands for average of intrinsic time of a given signal in the detector frame, and N W is a normalization factor 12 . Expression (30) gives the excess of the kurtosis of the intensity distribution relative to that of a normal distribution. Whatever the shape of the time profile at the source one expects that the kurtosis of the intensity distribution always changes in a certain way in the course of propagation. Namely, a burst-like signal evolves towards a platykurtic (flattened) type GeV × (s τ of intensity distribution (27) as it propagates in a dispersive medium 13 . Therefore, we use the compensation procedure described in Section 4, see (26), to return the shape of the intensity distribution towards the maximally leptokurtic (peaky) type 14

Skewness Estimator
Skewness is a measure of the degree to which a distribution is asymmetrical. In Section 2 we used a symmetric distribution, namely a Gaussian, to model a burst-like feature at a GRB source, see the solid line in Fig. 1. The dashed and dashed dotted lines in Fig. 1 are asymmetric distributions showing how this Gaussian envelope evolves when propagating trough a dispersive 13 It is convenient to calculate the kurtosis (30) of a distribution relative to that of a normal distribution. However, this does not entail any assumption about the actual initial shape of the time profile of a GRB. 14 We note that in other applications (see for example [50]) the kurtosis is used as a measure of the tail of a distribution, rather than the shape of the peak mediim with a power-law energy spectrum. The asymmetry may be measured using skewness (see for example [49]), which takes the following form for the intensity distribution defined in (27): where a weight W i given by (28) is assigned to every data point, b df (τ ) stands for the average of the detector-frame intrinsic timing of a given signal, and N W is a normalization parameter 15 .
Mathematically, the skewness is the ratio of the third moment of a distribution to its second moment raised to the power 3/2. The dispersed distribution shown in Fig. 1 is described as negatively skewed (or skewed to the left). In general, dispersion of the form (1) causes the skewness of a signal with a burst-like feature to become more negative. Therefore, whatever the shape of the time profile at the source and its degree of symmetry, one expects that the skewness of the intensity distribution always changes towards more negative values in course of propagation of the signal due to dispersion. Conversely, the compensation procedure of Section 4 tends to increase the skewness of the intensity distribution towards positive values, and we consider the optimal value of the compensation parameter in (26) to be that maximizing the skewness.
Examples of the values of S(τ ) calculated for different points τ j in a grid of values for the compensation parameter are shown in Fig. 5. The values of τ maximizing the skewness of the intensity distributions are considered to be those that best recover the original signal at the source. The estimates are made for the same GRBs as in cases of the irregularity and kutosis estimators discussed previously, namely GRB090510A and GRB080916C, and we find values of τ that are similar to those found previously The skewness estimator we utilize here is fully non-parametric, and does not rely on any assumption about the shape of the time profile at the source.

Uncertainties in the Estimators
In this Section we quantify the stability of the estimators described above with respect to the performance of Fermi-LAT [31] and account for the bias-induced systematic uncertainty in the overall measurement of the compensation parameter.
We comment first that the shortest timing shift in our studies is expected to be at the level of the smallest estimated |τ | 0.1 s/GeV multiplied by the lower energy cut, 100 MeV, which amounts to about 1 ms. Since the time resolution of the instrument is better than 10µs, we may assume that our results are insensitive to the timing accuracy. However, the evolution of the timing patterns during the propagation of the signals depends upon the energies of individual photons. Thus, the energy resolution of the instrument can influence the stability of the estimation of the correct value of the compensation parameter.
In the following we apply the estimators discussed in the previous Section to toy data sets generated by smearing the energies of the individual photons using a model resolution function, so as to assess the instability of the estimated compensation parameters. For this purpose we use one of the P8R2 V6 energy resolution performance plots from [51], and parameterize empirically the energy resolution for 68% half-width containment of the reconstructed incoming photon energy as: where x ≡ log 10 (E/MeV). Inaccuracy of the energy measurements superimposes an instability into the estimations of the correct value of the compensation parameter. We assign to every photon within an emission episode an energy generated randomly using a normal distribution defined by the mean value of its observed energy E and the standard deviation derived from (32). To maximize confidence in the accuracy of estimates, we have analyzed ∼ 100 thousand toys for every individual source. Examples of distributions of the corresponding values of the compensation parameters obtained for the toy data sets using the different estimators are shown in Fig. 6. For every source we apply five different estimation procedures based on the estimators  Three distributions in Fig. 6 are obtained from the irregularity estimator described in Section 5.1, using different methods to analyze the KS difference curve. The issue is that KS difference curves like those in Fig. 3 are quite irregular, in particular near the peak. This is due to the fact that, in order to avoid an unwanted systematic, we utilize different uniformly-distributed reference samples for every choice of τ . These irregularities can introduce an ambiguity in the estimation of the position of the maximum of the D(τ ) function derived as in (29). We utilize three methods to estimate the positions of the maxima of D(τ ) curve. The first is simply to define the position of a maximum as the centre of a segment formed by a horizontal line cutting the D(τ ) curve at a certain fraction of the total height of the D(τ ) curve 16 . The resulting distributions for two GRBs are shown by (magenta) dotted lines in Fig. 6. Another method is to define the maximum by the weighted average of the top part of the D(τ ) curve after being cut by the same horizontal line, which is shown by the (green) short-dashed three-dotted lines in Fig. 6. Finally, we also used the kernel density estimation (KDE) technique [52], which provides an estimate of D(τ ) within its whole support. The maxima of the KDE curves are then used to estimate the correct values of the compensation parameter, with the results shown by (orange) dashed lines in Fig. 6.
Unlike the KS difference curves, the kurtosis curves (see Fig. 4) and the skewness curves (see Fig. 5) are quite regular over the whole support, so that the positions of the maxima for K(τ ) and S(τ ) are unambiguous. The resulting distributions for kurtosis and skewness are presented in Fig. 6 by (black) long-dashed and (blue) short-dashed-dotted lines, respectively.
Another measure of the robustness of an estimator is its bias. In our case, this is expressed by the deviation of the estimate of the correct value τ recov of the compensation parameter from its true value τ true : where the average τ recov stands for the expected recovered value over a large number of repeated experiments.
We study the bias by analyzing a variety of realizations of the emission from a given source generated by sampling the timing and energy distributions with respect to their distributions in the data. The method of generating realizations used here is akin to the "flux randomization" procedure initially prescribed in [53] and later applied in [54] for simulations of Fermi-LAT GRB light curves. Following the general idea of [53], we simulate realizations of a given signal in such a way that the average temporal and energetic characteristics of each realization are equal, at some level of accuracy, to a specific timing-energy distribution obtained from the data. The timing-energy distributions for different realizations are obtained as random numbers along two axes, distributed according to the cell-contents of two-dimensional progenitor histograms with N 2 cells. (We recall that N is the number of events arriving at Fermi-LAT from a given source.) The total energy of the photons composing a generated realization is set to be equal to the total energy of the original signal with a certain accuracy. The progenitor histograms are constructed from the detected patterns with timings compensated for the propagation effect Fig. 7. The amount of the compensation is defined by the "correct" value of the compensation parameter found in any given estimation procedure. Examples of such realizations generated from the data from GRB090510A and GRB080916C compensated for the results of the kurtosis estimation procedure are presented in Fig. 7.
The realizations produced by the compensated data are regarded as being unaffected by the dispersion effect. In general, since the compensation values are different for the five estimation procedures used, one should study five different sets of realizations of the detected emission episode for every source. However, in practice, in order to reduce CPU time, for each GRB we use only one progenitor histogram compensated with respect to the result of one particular estimation technique, and offsets of the compensations for other estimators are taken into account in the calculations of the final uncertainties attributed to the bias corrections. Once the progenitor histogram is obtained we produce a number of realizations with a common degree of dispersion corresponding to a particular injected value of τ true . The total energy of the generated photons is required to be the same, to within 15%, as that measured in the data.
We then apply our estimation procedures to the set of realizations with a given injected dispersion signal, and calculate the average over the set of toy realizations of the estimated correct value of the compensation parameter. This average represents the expected recovered value of the compensation parameter τ recov in (33). Several reference values of τ true , each injected into separately-generated sets of realizations, are tested for every estimation procedure.
Parameters of straight line fits to τ recov versus τ true , like that ones shown in Fig. 8 and Fig. 9, are used for the determination of the uncertainties related to bias in the estimates of the correct values of the compensation parameters obtained using our estimation procedures. The required value of B(τ ) given by (33) is given by the difference between a given value of τ and the value of the linear fitted function calculated at the same τ . Finally, the bias calculated in this way for every estimation procedure is included in the uncertainties of the estimates of the compensation parameter that we present. The impacts of the bias corrections are illustrated in the lower row of Fig. 6, to be compared with the upper row, where no bias correction has been applied.
Almost identical work flows were used to estimate bias uncertainties for all estimators in the data from all the sources analyzed. We use for illustration results for the progenitor histogram for GRB080916C, shown in Fig. 7 by stars, compensated for the result obtained using the  Figure 9: Same as Fig. 8, but for the irregularity estimator with its maximum values calculated by averaging of the tops of the KS difference curves.
kutosis estimation procedure, namely τ = 0.85 s/GeV (the pattern originally detected is shown by the solid squares in Fig. 7). For the bias study, we use five reference values of τ true , each injected into a separate set of 15 realizations of the detected emission generated from the progenitor histogram. Thus, a total of 75 realizations has been generated for the bias study of GRB080916C. One of the realizations (modified with τ true = 0.4 s/GeV) is shown in Fig. 7 by downward-(upward-)pointing triangles. For every realization we apply the kurtosis estimation procedure to a set of 16 thousand toys generated with the energy smearing procedure described earlier, to obtain the optimal values of the compensation parameters 17 . The final distribution of the recovered value of the compensation parameter τ recov for a particular amount of injected dispersion τ true is obtained as the average over all 15 realizations within a single set. The result of the processing for this specific object using the kurtosis estimator is shown in the right panel of Fig. 8. The difference between the outcome of the kurtosis estimator for the data of GRB090816C, τ γ0 = 0.8 s/GeV, and the value of the function obtained from the straight line fit with errors related to the fit added in quadrature yields an uncertainty of ±0.069 s/GeV, to be compared with ±0.048 s/GeV when the data are used directly. The irregularity estimator is more affected by bias, as seen in the right panel of Fig. 9. The estimators are least biased in the case of GBR090510A (see the left panels of Fig. 8 and Fig. 9).
As can be seen from the examples in Fig. 6, the precisions of the different estimation techniques are quite similar to each other, although differences appear at the 1σ level, in particular when the bias is not taken into account (upper row of Fig. 6). We attribute these differences to an unidentified systematic that is probably related to the fact that different estimators deal with different kinds of deformation of the signal envelope. To be conservative, instead of giving a preference to any particular estimator, we simply average the results of the five estimation procedures for each energy smearing toy. In this way, we obtain the overall distributions shown as the solid lines in the upper panels of Fig. 6. The results obtained from the bias-corrected distributions shown in the lower panels of Fig. 6 as solid lines are used for combination studies in the next Section.

Consolidated Distribution and the Robust Limit
Our final goal is to infer the common degree of quantum-gravity-induced dispersion that is most compatible with the estimates obtained for all the sources we have analyzed. We note that the relation (25) implies that correct value of the compensation parameter τ (z) obtained for sources at different red shifts can be adjusted to the value at a reference red shift z 0 via the scaling τ (z 0 ) = τ (z) We apply this adjustment to every toy model generated by energy smearing for each source. For simplicity, we choose z 0 = 0.8944, which corresponds to K 1 (z 0 ) = 1. In this case the compensation parameter can be converted trivially into the main parameter of interest, namely, the scale of violation of Lorentz invariance We have transformed the overall distributions, point by point, into a combined distribution of values (34), which we call the K-reduced distribution. The overall distributions of all the sources entered in the analysis together with their K-reduced versions are presented in Fig. 10.
We now address the problem of consolidating our measurements of the compensation parameter obtained for different sources. In general, we would need to minimize a likelihood function to give an estimate for the distribution of the τ (z 0 ) that combines the information of the individual K-reduced overall probability distribution functions (PDFs). However, given that our individual K-reduced distributions are very close to Gaussian as seen in the left panel of Fig. 11 18 , we can use the minimum χ 2 method to obtain the consolidated PDF.
Assuming that there are no correlations between the measurements of different sources, one can construct a common χ 2 function:   Table 1), and their K-reduced normalized versions (right panel). Every individual distribution is obtained as an average of probability distributions obtained using the five estimation procedures described in Section 6. Lower row: the same distributions as in the top row, corrected for the bias systematic as described in Section 6. where τ i and σ τ i are the means and the standard deviations of the individual Gaussians (see Table I for details) and τ is the mean of the consolidated distribution that minimizes (36). It is well known that the solution is given by the weighted average and the standard deviation of the consolidated distribution is given by The results (37) and (38) can be proved as theorems in the framework of conflation [55], which provides a recipe for combining the PDFs of different measurements on a point-by-point basis.
Consolidating the K-reduced PDFs of our sources, using the prescriptions (37) and (38), one arrives at a large raw χ 2 value of 261, which implies a negligible probability for the individual distributions entering in the combination to be compatible with each other, implying that the sources are not identical. Each emission episode is affected by an intrinsic process which might introduce ether stochastic or systematic scatter of the results for individual sources.
The nature of the radiative processes and energy dissipation mechanisms of GRBs have not been clearly identified yet, which limits our ability to model the temporal spectral properties of the emitting region. Without additional inputs on the physics of the processes responsible for high energy emission of GRBs [41], one can only assume that there are some source-dependent contributions to the spectral evolution of individual sources that could be responsible for the mistuning we find in the K-reduced distributions of the compensation parameters obtained for different sources.
In our ignorance, we estimate the possible uncertainties that might be introduced by such unknown effects, in two ways, namely, using two different scalings of the individual distributions to render them compatible with a single overall consolidated distribution. The first possibility is to rescale the standard deviation of the individual distributions by a factor χ 2 raw , so that the resulting χ 2 scaled becomes unity [48]. The corresponding re-scaled K-reduced distributions together with the consolidated one are presented in the left panel of Fig. 11. One can see that the region of 95% incompatibility with a zero result for the correct value of the compensation parameter lies beyond the line τ (z 0 )[95%CL] = 0.54 s/GeV, which corresponds to the following lower limit on the scale of linear Lorentz violation: An alternative way of taking into account unknown source-related intrinsic temporal spectral variations would be to allow for an additional universal stochastic spread of the PDFs. This may be achieved by adding in quadrature, for all the PDFs of the sources entered in the analysis, a universal variation in the τ (z 0 ) distributions whose normalization is fixed so that the overall χ 2 1. We estimate this standard deviation to be 2.30 s/GeV, and the corresponding rescaled K-reduced distributions together with the consolidated one are presented in the right panel of Fig. 11. In this case, the region of 95% incompatibility with zero result exceeds τ (z 0 )[95%CL] = 1.86 s/GeV, which corresponds to the following lower limit on the scale of linear Lorentz violation: which is significantly weaker than (39).
We note that in the case of the χ 2 raw re-scaling method the consolidated result is most affected by sources with lower variations in the K-reduced distributions. On the other hand, in the case of the method of adding of a universal stochastic spread the result is depends more equally on the different sources, since those with narrower distributions are expanded more substantially than in the re-scaling case.

Discussion of the Results
The mean values of the compensation parameter τ and its standard deviation ∆τ encoded in the K-reduced distributions define the sensitivity of the source sample to propagation effects  Fig. 12, and the data measurements are indicated by crosses. In the following, we perform some simple simulation exercises to assess the sensitivity which would be achieved if the statistical realization of the measurements was different, but assuming the same pattern of (τ, ∆τ ) distribution from high-energy GRBs with known red shifts as has been measured by

Fermi-LAT.
We first assess what would be expected if we had another realization of the current data set. We generate 100 thousand realizations of the pattern of 8 sources, as shown by triangles in the left panel of Fig. 12. Every realization containing only simulated measurements (triangles in the left panel of Fig. 12) has then been processed using the χ 2 raw rescaling method described in the previous section. The resulting distribution of 95% CL limits obtained for realizations with χ 2 raw rescaled measurements is presented in the left panel of Fig. 13. One can see that in 95% of the cases the limits fall below 1.4 × 10 18 GeV (indicated by the vertical dashed dotted line), which we interpret as an effective limit on the sensitivity of the pattern of measurements we have in our disposal. In the other words, whatever the red-shift distribution, the spectral and temporal content found for 8 emissions leading to a (τ, ∆τ ) distribution similar to the current data, this is the best sensitivity one could reasonably expect to achieve, which we term the sensitivity end-point.
If instead we process the realizations of eight measurements by adding a universal stochastic spread, as described previously, we find that the sensitivity end-point is at 3.3 × 10 17 GeV, as one can see in the right panel of Fig. 13.
The most probable values of the 95% CL constraints are quite similar in the two cases, namely 3.2 × 10 17 GeV in the case of χ 2 raw rescaled measurements and 2.1 × 10 17 GeV for processing with the universal stochastic spread. Doubling the number of simulated measurements in the realizations (see right panel of Fig. 12), however, we find a sensitivity end-point of 1.0 × 10 18 GeV for χ 2 raw rescaling and 2.6 × 10 17 GeV for universal stochastic spreading, although the most probable values of the 95% CL limit stay unchanged.
In practice, working with actual data, it is important not to underestimate the uncertainties at each step in the analysis flow. In particular, since the temporal distributions of the highenergy emissions of GRB engines are still poorly understood (see [41]), one has to be careful when cross-correlating directly the Fermi-LAT multi-GeV events with sub-MeV light curves detected by the the Gamma-ray Burst Monitor (GBM) [56]. In general, the paucity of multi-GeV photons makes it difficult to assess the importance of variability and temporal correlations with the emissions at lower energies. The latter implies that common features of signals in the sub-MeV and multi-GeV spectral bands could be established within some time intervals [57] that exceed substantially the time resolution of the detectors. This ultimately implies an uncertainty whose neglect can lead to an overstated assessment of the significance of the measurement obtained on the basis of an analysis [37][38][39][40] cross-correlating sub-MeV and multi-GeV photons.
One can also study the potential impact of accumulating more GRBs with K-reduced compensation parameter measurements that agree with the pattern of the eight sources we have analyzed. As a first exercise, we assume that the existing statistics are doubled so that eight additional measurements, like those indicated by the triangles in the left panel of Fig. 12, are available to be processed along with the current eight measurements indicated by the crosses  Fig. 14). However, the distribution of the 95% CL for the χ 2 raw re-scaling method looks rather smooth and symmetric, which slightly decreases the most probable value of the constraint to 6.0 × 10 18 GeV.
It is also instructive to perform two extreme exercises. One is to add just one simulated source to the present eight sources, processing the realizations with the χ 2 raw rescaling method. The resulting distribution, presented in the left panel of Fig. 15, clearly indicates that the most probable value of the 2σ limit is substantially boosted towards to higher values, namely to 7.3×10 17 GeV, getting quite close to the limit (39). On the other hand, a substantial increase of statistics, modelled by adding 28 sources to the eight present sources would lead to a distribution rather similar to that one shown in the left panel of Fig. 13, with the most probable value of the limit approaching that one obtained from the statistics of the present data alone. We recall that χ 2 raw re-scaling method of obtaining limit weights mostly the measurements with lower variances. Therefore, simulating one additional source provides, in most of the realizations, one additional measurement with low variance that increases substantially the constraint obtained.
In this sense, the example with one additional simulated source is evidence of an instability in conclusions about quantum-gravity effects on photon propagation drawn from analysis of a single GRB [32][33][34]36].
In closing this discussion, we emphasize that our analysis was performed in the context  Figure 14: Lefts panels: Distribution of the 95% CL limits obtained in 100 thousand realizations with eight (upper) and 16 (lower) simulated measurements added to the current data, processed with the χ 2 raw re-scaling method. We note that 95% of the entries do not exceed the value indicated by the vertical dotted dashed line. Right panels: The same as in the left panels, but processing the real and simulated data by adding a universal stochastic spread. of a model for space-time foam that does not predict birefringence, so that the strong constraints [22][23][24][25][26] are inapplicable. That said, the statistical techniques developed here could be applied to a wide class of models for Lorentz violation that predict anomalous dispersion in vacuo, providing model-independent constraints that are complementary to those from searches for birefringence.

Conclusions
We have developed in this paper three distinct statistical non-parametric measures of GRB emissions, which we have used in an analysis of Fermi-LAT data to search for the possible effect of quantum gravity on the propagation of high-energy gamma rays from GRBs. The measures utilize different types of deformation of the intensity profile of an envelope of electromagnetic radiation with a burst-like feature that would arise from propagation through a dispersive quantum-gravity medium. Applying five different estimation procedures developed on the basis of these statistical measures to the eight observed GRBs that are relatively bright in multi-GeV energies detected by Fermi-LAT, we constrain the possibility of a non-trivial vacuum refractive index for photons. Depending on the method of consolidation of the results for individual sources, we find that the energy scale M 1 characterizing a linear energy dependence of the refractive index should exceed either 8.4 × 10 17 GeV or 2.4 × 10 17 GeV. We have also made simple numerical exercises to explore the possible sensitivity of the current statistics of Fermi-LAT sources with measured red shifts together with sources that might be detected in the future, finding that the sensitivity would probably not exceed significantly M 1 ≈ 10 18 GeV.