Search for Lorentz Invariance breaking with a likelihood fit of the PKS 2155-304 flare data taken on MJD 53944

Several models of Quantum Gravity predict Lorentz Symmetry breaking at energy scales approaching the Planck scale (∼10 19 GeV). With present photon data from the observations of distant astrophysical sources, it is possible to constrain the Lorentz Symmetry breaking linear term in the standard photon dispersion relations. Gamma Ray Bursts (GRB) and flaring Active Galactic Nuclei (AGN) are complementary to each other for this purpose, since they are observed at different distances in different energy ranges and with different levels of variability. Following a previous publication of the High Energy Stereoscopic System (H.E.S.S.) collaboration [1] , a more sensitive event-by-event method consisting of a likelihood fit is applied to PKS 2155-304 flare data of MJD 53944 (July 28, 2006) as used in the previous publication. The previous limit on the linear term is improved by a factor of ∼3 up to M QG l > 2.1 × 10 18 GeV and is currently the best result obtained with blazars. The sensitivity to the quadratic term is lower and provides a limit of M QG q > 6.4 × 10 10 GeV, which is the best value obtained so far with an AGN and similar to the best limits obtained with GRB. Keywords Active galaxies PKS 2155-304 H.E.S.S. Quantum Gravity Lorentz Invariance breaking 1 Introduction 1.1 Lorentz Invariance and Quantum Gravity A development of new space-borne and ground-based experiments in the last decade, covering a large domain of gamma-ray astronomy, opened a new window for fundamental physics tests. In particular, the confrontation of the results from astrophysical sources with challenging predictions of various models in the frame of Quantum Gravity theories became of general interest [2] . The theory of Quantum Gravity (see e.g. [3] ), which is still incomplete, would give a unified picture based on both Quantum Mechanics and General Relativity, thus leading to a common description of the four fundamental forces. Quantum Gravity effects in the framework of String Theory [4] , where the gravitation is considered as a gauge interaction, result from a graviton-like exchange in a classical space–time. In most String Theory models involving large extra-dimensions, these effects would take place at the Planck scale, thus not leading to “spontaneous” Lorentz Symmetry breaking, as may happen in models with “foamy” structure [5,4] of quantum space–time. In this second class of models, photons propagate in a vacuum which may exhibit a non-negligible refractive index due to its foamy structure on a characteristic scale approaching the Planck length or equivalently the Planck energy E P = 1.22 × 10 19 GeV. This implies a group velocity of light changing as a function of energy of the photons, in analogy to the dispersion effects in any theoretical description of plasma media. On the other hand, in models based on General Relativity with Loop Quantum Gravity [6,7] which postulate discrete (cellular) space–time in the Planckian regime, the fluctuations would introduce perturbations to the propagation of photons. The light wave going through the discrete space matrix would feel an induced perturbation which increases with decreasing wavelength or equivalently with increasing energy of photons. In consequence, photons with different wavelengths propagate with different velocity. As a result, one may expect a spontaneous violation of Lorentz Symmetry at high energies to be the generic signature of the Quantum Gravity. In four dimensions, the Quantum Gravity scale is presumed to be close to the Planck mass and the standard photon dispersion relations up-to second order corrections in energy can be written as: (1) c 2 p 2 = E 2 ( 1 ± ξ ( E / M ) ± ζ ( E / M ) 2 ± ⋯ ) , where ξ and ζ are positive parameters. This implies an energy-dependent speed-of-light with (2) v = c ( 1 - ξ ( E / M ) ) , considering only the linear term in the expansion. Alternatively assuming the linear correction equal to zero, the quadratic term leads to: (3) v = c ( 1 - ζ ( E / M ) 2 ) . High-energy photons could propagate either slower (the sub-luminal case) or faster (the super-luminal case) than low-energy photons. In this study, a causality conserving scheme was adopted following to the theoretical arguments of the non-critical string models, developed in [8] . The sign minus in Eqs. (2) and (3) is required to avoid Cherenkov radiation in vacuum [9] . As results of this analysis and as in case of most of other quoted results in Table 1 , only limits on Quantum Gravity scale for the sub-luminal photons will be given. As suggested by Amelino-Camelia [5] , the tiny effects due to Quantum Gravity can add up to measurable time delays for photons from cosmological sources. The energy dispersion is best observed in sources that show fast flux variability, that are at cosmological distances and are observed over a wide energy range. This is the case of Gamma Ray Bursts (GRB) and Very High Energy (VHE) flares of Active Galactic Nuclei (AGN). Both types of sources are the preferred targets of these time-of-flight studies, which provide the least model-dependent tests of the Lorentz Symmetry. The case of pulsed emission by Pulsars has also been considered, and provided valuable results [10] . It should be underlined that studying both GRB and AGN is of fundamental interest. GRB are observed with satellites at very large distances (up to z ∼ 8) but with very limited statistics above a few tens of GeV. On the contrary, AGN flares can be detected by ground based instruments with large statistics up to a few tens of TeV. Due to the absorption of the high energy photons by Extragalactic Background Light (EBL), TeV observations are limited to sources with low redshifts. Thus, GRB and AGN are complementary since they allow to test different energy and redshift ranges. It should be noted that time-of-flight measurements are subject to a bias related to a potential dispersion introduced from the intrinsic source effects, which could cancel-out or enhance the dispersion due to modifications to the speed of light. When considering sources at cosmological distances, the analysis of time lags as a function of redshift requires a correction due to the expansion of the universe [11] , which depends on the cosmological model. Following the analysis of the BATSE data and recently of the HETE -2 and Swift GRB data [12–15] and within a framework of the Standard Cosmological Model [16] with a flat expanding universe and a cosmological constant, the difference in the arrival time of two photons emitted at the same time is given by (4) Δ t Δ E ≈ ξ E P H 0 ∫ 0 z ( 1 + z ′ ) dz ′ Ω m ( 1 + z ′ ) 3 + Ω Λ , for a linear effect, and by (5) Δ t Δ E 2 ≈ 3 ζ 2 E P 2 H 0 ∫ 0 z ( 1 + z ′ ) 2 dz ′ Ω m ( 1 + z ′ ) 3 + Ω Λ , for a quadratic effect. In the present study the cosmological parameters were set to Ω m = 0.3, Ω Λ = 0.7 and H 0 = 2.3 × 10 −18 s −1 . Δ E and Δ E 2 represent linear and quadratic energy ranges respectively. 1.2 Present results The most important results obtained so far with transient astrophysical sources are given in Table 1 . The most robust limits to date have been obtained by Ellis et al. [13,14] and Bolmont et al. [15] , with the use of several GRB at different redshifts, resulting in limits of the order of M QG l > 10 16 GeV. In case of a significant detection, using several sources at different redshifts allows in principle to take into account the possible time-lags originating from source effects. Fermi is now in operation and has provided the best limits so far [17,18] . However, for the moment the results are obtained with no study of the variation of the time-lag with redshift and with limited statistics in the GeV range. The recent studies with flares of Mrk 501 [19,20] and the former result of H.E.S.S. with PKS 2155-304 [1] provide results setting a lower limit on the linear Quantum Gravity scale at a few 10 17 GeV for each source. Table 1 also shows the different analysis techniques used to measure the energy-depend time-lags. The analysis methods can be divided in two categories, analyses applied on binned data on the one hand and unbinned analyses on the other hand. Binned analyses are applied on binned light curves and consist of finding the position of local extrema in two different energy bands and deducing the value of the lag. A simple fit [21] or a more sophisticated technique such as Continuous Wavelet Transform (CWT) can be used to localize extrema. CWT was used for GRB [12–15] and PKS 2155-304 [1] . The Modified Cross Correlation Function (MCCF) was also used for PKS 2155-304 to determine directly the lag from two light curves in two different energy bands [1] . Unbinned analyses have also been used for both GRB and AGN. Likelihood fits require a parameterization of the light curves and this parameterization can be taken from a binned [20] or unbinned fit [22] to the light curve. Another method, the Energy Cost Function (ECF), was recently used by Albert et al. [19] . In this method one introduces a time lag τ for each photon and adds up their energies if they fall into a given time interval. The maximum of ECF ( τ ) gives the value of the measured lag. The recently developed method called DisCan [23] was used for the analysis of Fermi data [18] . As discussed in [23] , this method applies successfully to the analysis of very sharp peaks as can be found in GRB light curves. It is important to point out here that the use of at least two methods is recommended since different methods can probe different aspects of the light curve. The present analysis completes previously published studies [1] where two different methods were employed. As no significant time-lag was detected so far, careful error calibration studies using Monte Carlo simulations are mandatory for the extraction of the limits, as will be further discussed in the following. 1.3 Overview of the paper In this paper, a new analysis of PKS 2155-304 data is presented, using a likelihood fit to determine both linear and quadratic corrections to the dispersion relations. In Section 2 , the method is described. The data in use as well as the selections applied are detailed in Section 3 . Then, in Section 4 , the procedure used to calibrate the method and to evaluate systematics is presented. The results are given in Section 5 and finally discussed in Section 6 . 2 Method description To study the correlations between the arrival time and the energy of the photons, a method relying on a probability density function (p.d.f.) on an event-by-event basis is used. Following Martinez and Errando [20] , the p.d.f. is defined by: (6) P ( t , E ) = N ∫ 0 ∞ A ( E S ) Λ ( E S ) G ( E - E S , σ ( E S ) ) F S ( t - τ l E S ) dE S , where Λ ( E S ) is the emission spectrum, G ( E − E S , σ ( E S )) is the smearing function in energy (a usual point spread function for an energy resolution of 10%), A ( E S ) is the acceptance of the H.E.S.S. array, F S is a parameterization of the emission time distribution (or light curve) and N is a normalization factor. The time shift was re-parameterized as Δ t = − τ l E and Δ t = − τ q E 2 for linear and quadratic effects respectively, where τ l and τ q are expressed in s TeV −1 and s TeV −2 . The value of N was computed in view of unbiased τ -parameter estimation. It guarantees the correct normalization of the p.d.f. and takes into account its model dependence. Here, integration over measured domain of energy and time has been performed. The function F S is a function of time at the source t S = t − τE S . The probability P ( t , E ) depends only on the relative delay of photons for different energies. As the light curve at the emission time is not known, the best estimate of the time structure at the source can be taken from the time distribution of the lowest energy photons. In the propagation models, the change of the QG scale induces basically a sliding effect in time of the light curve. This point has been checked to be valid by the Monte Carlo simulations described below. An important assumption of this study should be underlined: the intrinsic light curves have been considered to have the same shape in different energy bands. As no lag is found, this can be partially justified by comparing the light curve at low energy ( Fig. 2 ) and in the full energy range of the instrument where the statistics is five times higher [24] . In both cases the light curve shows a five peak structure with maxima located roughly at the same times. The likelihood function is given by the product of probability density functions over all photons: (7) L = ∏ i P i ( t , E ) . The maximum of L provides the time-lag parameter τ l and τ q for linear and quadratic models respectively. The likelihood method is in principle very sensitive and allows a lag to be measured even if the number of selected photons is limited, which is not the case for wavelet or MCCF analyses. As a consequence of this high statistical sensitivity, it would be possible to probe a light curve locally, introducing a selection in time. The main aspect of the method to be taken into account is a need for a very precise parameterization of the photon distributions. 3 Data sample used in the analysis PKS 2155-304 is a high-peaked BL Lac object located at redshift z = 0.116. On July 28, 2006 (MJD 53944), an extreme flare of this source was observed by H.E.S.S. [24] . The full data sample was obtained using a combination of two analysis methods [25] : the standard Hillas analysis [26] and the Model analysis [27] . This technique greatly improves the signal to background ratio and results in 8148 on-source events recorded in ∼85 min (three data-taking runs) with an energy threshold of ∼120 GeV. The signal-to-background ratio in the flare data exceeds 300. The effect of the background contribution will be investigated in the next section. The average zenith angle of selected events is θ z ≈ 11°. Table 2 summarizes the different cuts applied to the data. In this study, only the first 4000 s of data (7713 events) are considered, where the flux and variability are the highest. To further improve the signal-to-background ratio, a cut θ 2 < 0.005 deg. 2 is applied. This cut reduces the number of photons to 5974 events but allows a very good fit of the light curve and of the energy spectrum. A total of 3526 events pass the cut on θ 2 in the first 4000 s of data in the range 0.25–4.0 TeV. All events in the data sample fulfilling the selection were used in the analysis. As shown in Fig. 1 , the reconstructed spectrum is compatible with a power law of index Γ = 3.46 in the range of 0.25–4.0 TeV (with χ 2 /dof = 16.7/23). This result is compatible with the high energy part of the broken power law spectrum obtained for the same flare by Aharonian et al. [24] with a different event selection and energy reconstruction method, and also with the spectrum obtained for the whole 2005–2007 period [28] . The dependence of the spectral index as a function of time has been studied extensively in [24] and no significant variation was found during the flare. Fig. 2 shows the light curve in the range 0.25–0.28 TeV with a binning of 61 s. The shape of this light curve will be used later as a template for the time-lag determination (Section 5 ). To estimate the main parameters of the light curve (widths and asymmetry of the pulses, distance between consecutive spikes), a parameterization was performed with a sum of asymmetric Gaussian functions. For this, the number of spikes was first determined using a peak-finding procedure based on a Markov chain algorithm [29] . Five pulses were found and their positions were used as initial parameters for the light curve fit. The fit function F is a sum of asymmetric Gaussian curves, defined as follows: (8) F n ( t ) = ∑ i = 1 n f i ( t , A i , μ i , α i , β i ) , where n is the number of spikes, and where (9) f ( t , A , μ , α , β ) = y = Ae - ( t - μ ) 2 2 α 2 , if t < μ , y = Ae - ( t - μ ) 2 2 β 2 , if t ⩾ μ . In this definition, A and μ are the normalization factor and the position of the pulse and α and β are the left-hand and right-hand widths of the pulse respectively. The asymmetry of the curve can be quantified by the parameter ν = β / α . When ν = 1, f becomes a Gaussian function with a width α = β = σ 0 . The parameterization was carried out with a function F 5 plus a constant term B . Values of the parameters are given in Table 3 , leading to χ 2 /dof = 48.8/38. This fit gives the following information about the light curve: • width of the pulses in the range 250–600 s, • time difference between consecutive spikes in the range 500–1000 s, and • ratio ν in the range 0.6–7. These intervals were investigated for the calibration procedure described in the next section. 4 Calibration of the method A calibration procedure was developed using a toy Monte-Carlo program in order to reflect as closely as possible the real data analysis and to provide the best strategy to minimize the systematic errors on the estimated time lag. As discussed in the previous section, the light curve of PKS 2155-304 is well fitted by a superposition of five asymmetric Gaussian spikes. The energy spectrum follows a power law distribution. These parameters have been used as initial conditions for the simulation. The zenith angle was considered to be constant at θ z = 10°. A simulation procedure comprises the following steps: • the template F S is chosen to be either a Gaussian function or derived from the measured light curve in the range 0.25–0.28 TeV; • photons are generated from a parameterizations F S and Λ in the range 0.25–4 TeV; • for each photon, a lag Δ t is introduced depending on the energy; • a weight is associated to each photon taking into account the acceptance of the detector during data-taking in 2006; • the energy is smeared taking into account the H.E.S.S. energy resolution (Δ E / E ∼ 10%); • when a simple Gaussian light curve is generated, histograms obtained for the light curve and the spectrum are fitted. As a result, the functions used for the fits reproduce well those used for the generation. This point was checked to estimate the effects of the statistical smearing, which were found to be negligible; • finally, the likelihood is computed and the minimum of −2Δln( L ) is searched for, first using a Gaussian template parameterization for the detailed studies or the PKS 2155-304 template. As shown later, the distribution of the minima has always been found to be compatible with a Gaussian distribution (with a mean τ ¯ and a width σ τ ). For each simulation run, these steps are repeated 500 times with ∼2000 photons after selections in each realization. The injected lag ( τ 0 ) range was chosen between −100 and 100 s TeV −1 (respectively s TeV −2 ) with steps of 10 s TeV −1 (s TeV −2 ). The results of the simulations can be summarized studying the variation of the reconstructed lag τ ¯ as a function of the injected lag. Fig. 3 (left) shows such a plot for simulation settings which will be described later (Gaussian spike). The shaded area corresponds to a 68% CL range of the reconstructed lag. A fit within ± σ τ yields a slope parameter value close to unity. A small systematic shift of the slope is observed; however, the bias on the reconstructed time lag is much smaller than the statistical uncertainty σ τ for time lags between −50 s TeV −1 and 50 s TeV −1 . In the following sub-sections, different features of the error calibration procedure are detailed. First, the influence of the pulse shape (width and asymmetry) and the superposition of pulses is investigated. Then, the dependence on the spectral index is evaluated. In the end, the effect of background events is studied. The simulation of light curves with five pulses provide the value of the calibrated statistical error which is used in Section 5 to compute the limits on E QG . 4.1 Pulse shape and superposition of pulses Fig. 3 (left) shows the reconstructed lag as a function of the injected lag for a Gaussian pulse of a mean width σ 0 of 300 s for a linear model. The calibration slope is 0.97 ± 0.02, showing the likelihood fit provides efficient means for the reconstruction of a lag. No deviation from linearity is observed in this case. The same figure shows an example of the distribution of the minima of −2Δln( L ) fitted with a Gaussian curve, for an injected lag of 0 s TeV −1 . This distribution is fitted by a Gaussian curve with τ ¯ = 0.6 ± 0.2 s TeV −1 and σ τ = 5.1 ± 0.2 s TeV −1 . For a given value of the pulse width σ 0 , the spread σ τ is very stable. On the other hand, as shown in Fig. 3 (right), σ τ increases rapidly with the width of the pulse. For a single pulse width of 600 s, σ τ rises to ∼13 s TeV −1 . When assuming asymmetry in the pulse shape, the calibration slope remains very stable and close to unity. The maximal asymmetry factor observed for PKS 2155-304 is ν = 7, leading to σ τ ≈ 10 s TeV −1 for σ 0 of 300 s. The superposition of spikes requires detailed studies concerning the evolution of the calibration slope and the error in case of complex light curves such as the one measured with PKS 2155-304 flare. In order to check the effect of spike superposition, light curves were generated from the sum of several asymmetric Gaussian curves. Various configurations were tested and the resulting calibration curves were compared to those obtained when injecting the real PKS 2155-304 light curve. The plots of Fig. 4 were obtained with light curves generated from real data in the range 0.25–0.28 TeV for the linear and quadratic modelling. These simulations were used to investigate the effect of superimposing two to five pulses with realistic parameters leading to similar calibrated error values. In the quadratic case, the slope is somewhat smaller than 1. The spread σ τ is smaller than in the linear case close to τ 0 = 0 but increases at larger values of ∣ τ 0 ∣ where the calibration curve becomes non-linear. In the linear case, a systematic shift of the calibration values of the order of 1/2 σ τ is observed. This effect will be included in the overall calculation of the systematic errors (see Table 4 ). The error values obtained from Fig. 4 were used in the following to derive the limits for linear and quadratic corrections. The error σ τ was found to depend mainly on individual pulse widths rather than on the addition of the widths of the five pulses, leading to 10.9 s TeV −1 for the linear and 6.3 s TeV −2 for the quadratic case respectively. Anticipating the results from data, these values were taken in the vicinity of the injected time-lag τ 0 = 0. In addition, other functional forms were used when fitting template light curve: Norris parameterization as also used in [24] and Landau function. No detectable effect was found when deriving the results. Moreover, the effect of the different functional forms was evaluated with MC simulations. 4.2 Dependence on spectral index Varying the spectral index of the generated light curves from 1 to 4.5 resulted in no significant change in the calibration slope (below 5%) for a single Gaussian pulse. The systematic error due to the uncertainty on the spectral index determination can then be considered as negligible. In contrast, the error value is a strong function of the spectral index as also shown in Fig. 5 . A harder spectrum gives a higher precision, due to the fact that the photons are spread more uniformly between high and low energies. This enhances the potential of studies with closer AGN such as Mrk 421 and Mrk 501. 4.3 Effect of background contribution Even though the signal to background ratio is very high for the studied data sample, it is possible that a small contribution of mis-identified photons introduces a small bias in the measurement of the time-lag. To study the impact of this effect, light curves were simulated with a fraction of particles with no energy-dependant time-lag as is expected for charged cosmic-rays. A realistic power law spectrum with an index of 2.7 was used to generate these particles. When 1% of the particles are considered as protons, no significant change in either slope or parameter σ τ is observed. 4.4 Summary of calibration studies In this section, the most important features of the calibrations with the likelihood method are summarized. As the emission light curve at the source is not known, the function F S was derived from real data in the range 0.25–0.28 TeV. In this case, it is shown that the likelihood fit can reconstruct the injected lag extremely well. In the future, if a satisfying emission model is proposed, the likelihood approach can be very useful to test propagation and emission models as well. The error on the measured lag ( σ τ ) depends strongly on the width of individual pulses.This result is very important and allows competitive results to be obtained with the PKS 2155-304 light curve. On the other hand, the asymmetry of the pulses and the superposition of several pulses has been shown to have less important effect on the results. The parameter σ τ also depends on the spectral index, which is another source of uncertainty. However, for PKS 2155-304 flare in 2006, the spectral index can be determined with a very good precision. The systematic uncertainties in the method employed were quantified by varying selections and cuts and treating the light curve shape parameterization in different ways. The various contributions are detailed in Table 4 . The main uncertainty was found to be related to event selections and to the light curve parameterizations. In order to get a precise estimation of the systematic uncertainty related to the F S parameterization, the covariance matrix of the fit presented in Fig. 2 was studied. The values of the fit parameters were varied following to the gaussian distributions with parameters derived from the covariance matrix after its diagonalization. The values in Table 4 were obtained with this procedure for the linear and the quadratic case separately. The overall systematic error is estimated to be <10.3 s TeV −1 and <6.6 s TeV −2 for the linear and quadratic effects respectively. These values will be combined with the statistical calibrated error given in Section 4.1 when determining the Quantum Gravity scales. 5 Results with PKS 2155-304 flare of MJD 53944 The final result was obtained with a parameterization of the light curve in the range 0.25–0.28 TeV ( E 1 ¯ = 0.26 TeV, see light curve in Fig. 2 ) and a likelihood computed in 0.3–4.0 TeV ( E 2 ¯ = 0.49 TeV). The choice of these ranges was dictated by the balance found between lever-arm in energy and the statistics used in the likelihood fit. These choices were studied and optimized with the simulations described above. Fig. 6 shows the −2Δln( L ) curve as a function of the injected τ for the linear and quadratic cases, leading to the value for the minimum of τ 0 l ′ = - 5.5 ± 10.9 ( stat ) ± 10.3 ( sys ) s TeV - 1 , for the linear and τ 0 q ′ = 1.7 ± 6.3 ( stat ) ± 6.6 ( sys ) s TeV - 2 , for the quadratic effect. These results imply that there is no significant time-lag measured. Consequently, the 95% CL limits on the Quantum Gravity energy scales are derived: M QG l > 2.1 × 10 18 GeV ( ξ < 5.7 ) and M QG q > 6.4 × 10 10 GeV ( ζ < 3.6 × 10 16 ) . The more conservative 99% CL limits are M QG 99 % l > 1.3 × 10 18 GeV ( ξ < 9.2 ) and M QG 99 % q > 5.3 × 10 10 GeV ( ζ < 5.2 × 10 16 ) . These results were found to be very stable when considering different energy domains and time selections. The stability studies of the estimated systematic effects allows the quoted values to be considered as conservative. The limit on the linear term, while lower by a factor of ten in comparison with the one obtained by Fermi with GRB 090510 [18] , shows that with AGN too, the conclusion about the models predicting a linear effect with energy [4] tends to be confirmed. As far as the limit on the quadratic term is concerned, it is the best one obtained so far with AGN and GRB, even if it still remains far from the Planck scale value. As stated in the introductory part, the presented paper focused on models proposed in [12,13] where predictions of linear Lorentz Invariance Violation effects were expected and the causality condition was conserved. Moreover, it may be noticed that because of negligible negative values found both for the linear and quadratic terms, the super-luminal case limits will differ only slightly from those provided for the sub-luminal case. 6 Summary and prospects The limits obtained in this study are the most constraining ever obtained with an AGN. They are almost ten times higher than those obtained by Martinez and Errando [20] with the same method for the Mrk 501 flare recorded by MAGIC, and consistent with their rough estimate for the PKS 2155-304 flare observed by H.E.S.S. They are also a factor of ∼3 higher than the previous H.E.S.S. result [1] . The increase in sensitivity as compared to the MAGIC result is mainly due to excellent parameters of the data taken on July 28, 2006: low zenith angle (∼10°) which leads to a low energy threshold and high statistics with high variability. On the other hand, the steepness of the PKS 2155-304 energy spectrum was one of the penalizing factors in this study. Given the complexity of the PKS 2155-304 light curve, a detailed calibration study of the method was necessary. For this purpose, a toy Monte-Carlo program was developed. The simulations allowed the influence of different parameters on the likelihood performance to be investigated, such as the shape of the light curve or spike superposition. In particular, a very interesting conclusion for a light curve with multiple spikes is that the error on the lag depends mainly on the width of the individual pulses. This feature of the likelihood fit allows a smaller error on the lag measurement to be obtained, yielding higher limits on M QG . It should still be noticed that the superposition of different spikes leads to a decrease in sensitivity due to a limited number of “useful events” which contribute to the likelihood fit. Despite this loss in statistical power, the method allows robust and precise results on the fitted time-lag to be obtained. Another important conclusion about the method is the possibility to obtain robust results even with low photon statistics. The use of different time selections showed that a large decrease in the number of selected events entering the fit had only a small influence on the error value of the time-lag measurement. The main difficulty of time of flight studies with photons from AGN and GRB is that the origin of the measured lag is unknown. If the measured time-lag is not related to Quantum Gravity induced propagation effect, it can be due to an emission effect at the source. The measured lag can also be a superposition of source and propagation effects. As no lag was found in the studied sample, either the source effects compensate exactly the propagation effects (“conspiracy scenario”), or there are no source or propagation induced effects detected within present experimental sensitivity. In this study the second hypothesis was retained, which leads to a conclusion that any measured lag would be entirely due to Lorentz Invariance Violation effects. In future, modelling of the emission light curve could bring new insight in the propagation and emission time studies. Population studies of AGN examining the effect of the distance on the lag measurement could allow this assumption to be avoided but would require that intrinsic effects are similar for different AGN. Better results on propagation effects will then come in parallel with a better understanding of the emission models. More and more results on Lorentz Symmetry breaking are published which give limits for the linear correction that are very close to the Planck energy scale. In particular, the latest results from Fermi observations of distant GRB provide exclusion limits for a large class of linear models above the Planck energy [17,18] . However, as pointed out by Ellis et al. [30] , ‘extraordinary claims require extraordinary evidence’, so more studies will be needed in the future to give more robust results and to be able to definitively reject or validate proposed models. The time-lag measurements with photons are far from being sensitive enough to constrain the quadratic term of the dispersion relations and do not exclude any model in this case. Multi-messenger and multi-wavelength campaigns will be necessary in future to constrain this quadratic term. Combining results obtained with GRB and AGN detected by Fermi and ground-based detectors would allow considering common scheme with photon energies ranging from the sub-GeV range to the multi-TeV range and from sources with higher redshift values. Future observations with neutrino telescopes of GRB and AGN flares would further extend the energy range to extreme, multi-hundred TeV domain. However, at present, the multi-messenger and multi-wavelength procedures are still to be developed. Concerning the increase of the experimental potential in observations of the AGN, a key step towards a better understanding of the emission and the propagation effects will be achieved with a new generation of instruments, which will allow the observation of a larger number of AGN flares. HESS-II, with an energy threshold of ∼30 GeV and later the large Cherenkov array CTA 4 http://www.cta-observatory.org/ . 4 will increase dramatically the sensitivity to energy-dependent time-lags. A rough estimation obtained by extrapolating results from the PKS 2155-304 flare to the extended domain in energy from a few GeV to a hundred TeV, and assuming a gain in sensitivity of a factor of ten, would rise the limit on the quadratic term by two orders of magnitude. This type of results may be further confronted with present and future theoretical models. It may also be underlined here that this will open a new unexplored domain for Lorentz Invariance Violation searches. Acknowledgments The support of the Namibian authorities and of the University of Namibia in facilitating the construction and operation of H.E.S.S. is gratefully acknowledged, as is the support by the German Ministry for Education and Research (BMBF), the Max Planck Society, the French Ministry for Research, the CNRS-IN2P3 and the Astroparticle Interdisciplinary Programme of the CNRS, the UK Science and Technology Facilities Council (STFC), the IPNP of the Charles University, the Polish Ministry of Science and Higher Education, the South African Department of Science and Technology and National Research Foundation, and by the University of Namibia. We appreciate the excellent work of the technical support staff in Berlin, Durham, Hamburg, Heidelberg, Palaiseau, Paris, Saclay, and in Namibia in the construction and operation of the equipment. References [1] F. Aharonian H.E.S.S. Collaboration Phys. Rev. Lett. 101 2008 170402 [2] G. Amelino-Camelia, arXiv:0806.0339 , 2008. [3] L. Smolin, Three Roads to Quantum Gravity (Basic Books), 2001. [4] J. Ellis N.E. Mavromatos D.V. Nanopoulos Phys. Rev. D 61 1999 027503 [5] G. Amelino-Camelia Nature 395 1998 525 [6] J. Alfaro Phys. Rev. D 65 2002 103509 [7] R. Gambini J. Pullin Phys. Rev. D 59 1999 124021 [8] S. Sarkar Mod. Phys. Lett. A 17 2002 1025 [9] G. Amelino-Camelia J. Ellis N.E. Mavromatos D.V. Nanopoulos Int. J. Mod. Phys. A 12 1997 607 [10] P. Kaaret A&A 345 1999 L32 [11] U. Jacob T. Piran JCAP 01 2008 031 [12] J. Ellis N.E. Mavromatos D.V. Nanopoulos A.S. Sakharov A&A 402 2003 409 [13] J. Ellis Astropart. Phys. 25 2006 402 [14] J. Ellis Astropart. Phys. 29 2008 158 [15] J. Bolmont A. Jacholkowska J.-L. Atteia ApJ 676 2008 532 [16] N.A. Bahcall J.P. Ostriker S. Perlmutter P.J. Steinhardt Science 284 1999 1481 [17] A.A. Abdo, et al., Science Express, 02/19/2009, 2009a. [18] A.A. Abdo Nature 462 2009 331 including supplementary information [19] J. Albert MAGIC Collaboration Phys. Lett. B 668 2008 253 [20] M. Martinez M. Errando Astropart. Phys. 31 2009 226 [21] S. Boggs C. Wunderer K. Hurley W. Coburn ApJ 611 2004 L77 [22] P. Lamon INTEGRAL Collaboration Gen. Rel. Grav. 40 2008 1731 [23] J.D. Scargle J.P. Norris J.T. Bonnell ApJ 673 2008 972 [24] F. Aharonian H.E.S.S. Collaboration ApJL 664 2007 L71 [25] M. De Naurois, et al., in: Proceedings Towards a Network of Atmospheric Cherenkov Detectors VII (Palaiseau), 149, 2005. [26] A. Hillas, in: Proceedings of 19th I.C.R.C. (La Jolla), vol. 3, 445, 1985. [27] M. De Naurois, et al., in: Proceedings of 28th I.C.R.C. (Tsukuba), vol. 5, 2907, 2003. [28] A. Abramowski et al., H.E.S.S. Collaboration 2010, A&A, submitted. [29] M. Morhac NIM A 443 2000 108 [30] J. Ellis Phys. Lett. B 674 2009 83 [31] S.D. Biller Whipple Collaboration Phys. Rev. Lett. 83 1999 2108 [32] A. Pélangeon J.-L. Atteia D.Q. Lamb Gamma-Ray Bursts in the Swift Era. AIP Conf. Proc. 836 2006 149