Probing quantum gravity using photons from a flare of the active galactic nucleus Markarian 501 observed by the MAGIC telescope

We analyze the timing of photons observed by the MAGIC telescope during a flare of the active galactic nucleus Mkn 501 for a possible correlation with energy, as suggested by some models of quantum gravity (QG), which predict a vacuum refractive index \simeq 1 + (E/M_{QGn})^n, n = 1,2. Parametrizing the delay between gamma-rays of different energies as \Delta t =\pm\tau_l E or \Delta t =\pm\tau_q E^2, we find \tau_l=(0.030\pm0.012) s/GeV at the 2.5-sigma level, and \tau_q=(3.71\pm2.57)x10^{-6} s/GeV^2, respectively. We use these results to establish lower limits M_{QG1}>0.21x10^{18} GeV and M_{QG2}>0.26x10^{11} GeV at the 95% C.L. Monte Carlo studies confirm the MAGIC sensitivity to propagation effects at these levels. Thermal plasma effects in the source are negligible, but we cannot exclude the importance of some other source effect.


Introduction
It is widely speculated that space-time is a dynamical medium, subject to quantumgravitational (QG) effects that cause space-time to fluctuate on the Planck time and distance scales [1][2][3][4][5][6][7][8], for reviews see [9]. It has also been suggested that this 'foaming' of space-time might be reflected in modifications of the propagation of energetic particles, namely dispersive effects due to a non-trivial refractive index induced by the QG fluctuations in the space-time foam. There are microscopic string-inspired models [1][2][3] that predict only subluminal refraction, and only for photons [10], suppressed either linearly or quadratically by some QG mass scale: (1) One might guess that the scale M QG1 or M QG2 would be related toM P , wherê M P = 2.4 × 10 18 GeV is the reduced Planck mass, but smaller values might be possible in some string theories [2,3], or models with large extra dimensions [11]. Superluminal modes and birefringence effects are also allowed in some other models [4][5][6][7][8].
A favored way to search for such a non-trivial dispersion relation is to compare the arrival times of photons of different energies arriving on Earth from pulses of distant astrophysical sources [1,12]. The greatest sensitivities may be expected from sources with short pulses, at large distances or redshifts z, of photons observed over a large range of energies. In the past, studies have been made of emissions from pulsars [13], γ-ray bursts (GRBs) [1,11,12,[14][15][16][17] and active galactic nuclei (AGNs) [18,19]. In particular, a combined analysis of many GRBs at different redshifts made possible some separation between energy-and source-dependent effects, and yielded a robust lower limit M QG1 > 0.9 × 10 16 GeV [15]. Astrophysical sources that produce very high energy photons in the TeV range or higher could improve significantly the sensitivity to Lorentz violation, if one could distinguish source and propagation effects. Flaring AGNs are celestial objects with the desired properties, and a pioneering study of a flare of the AGN Mkn 421 yielded a sensitivity to M QG1 ∼ 4 × 10 16 GeV [18]. 4 In this Letter we analyze two flares of Mkn 501 (z = 0.034) observed by the Major Atmospheric Gamma-ray Imaging Cerenkov (MAGIC) telescope [22] between May and July 2005. After applying standard quality checks, data covering a total observation time of 31.6 h spread over 24 nights survived and were analysed [23]. The data were taken at zenith angles of 10 • −30 • , resulting in an energy threshold (defined as the peak of the differential event-rate spectrum after cuts) of ≈ 150 GeV. The air-shower events were subjected to the standard MAGIC analysis [24], which rejects about 99% of hadronic background events, while retaining 50−60% of the γ-ray induced showers. The γ-ray energies are, in a first approximation, proportional to the total amount of light recorded in the shower images; corrections are applied according to further image parameters [25] obtained from the analysis. The achieved energy resolution is ∼ 25% over the range 170 GeV to 10 TeV. The arrival time of each event is obtained with sub-ms precision.
During the observations, variations in the γ-ray flux by an order of magnitude were observed, with the maximum integrated flux above ≈ 150 GeV exceeding (11.0 ± 0.3)×10 −10 cm −2 s −1 . In the two nights with the highest flux, high-intensity outbursts of short duration (flares) were recorded with characteristic rise and fall times of 1−2 minutes. While the flare of July 9 was clearly visible over the full energy range 0.15−10 TeV and reached a peak flux more than a factor two higher than before and after the flare, that seen on June 30 was concentrated at low energies (0.25−1 TeV) and less significant. In the analysis below, we applied cuts on the image parameter ALPHA [25], describing the gamma shower arrival direction: |ALPHA| < 10 • , and on energy: E γ > 150 GeV.

Timing analysis
The spectral time properties of the most intense portions of the flares were quantified in [23] using four different energy bands with boundaries at 0.15, 0.25, 0.6 and 1.2 TeV, the fourth band extending to infinite energies. In the June 30 flare a signal above a uniform background appeared only in the energy band of 0.25−0.6 TeV, which did not permit any conclusion on the time-spectral properties of the signal. For the flare of July 9, a time lag of about 4 minutes was found for the maximum of the time profile envelope for photons in the 1.2−10 TeV energy band relative to those in the range 0.25−0.6 TeV. The difference between the mean energies in these two bands is ≈ 2 TeV, which would lead to a naive estimate of a time delay of about 0.12 s for each GeV of energy difference. However, this approach is too simplistic, since the energy range covered by the 1.2−10 TeV band is much larger than the energy difference between the two bands, so the binned estimator used in [23] is inadequate for constraining possible linear or quadratic energy dependences. In view of this and the limited number of photons, we improve here on the binned estimator by analyzing the complete information encoded in the time-energy distribution of individual photons in the flare, with the aim of probing possible systematic energy-dependent time lags induced by QG vacuum refraction during photon propagation to the Earth, or intrinsic to the source. The true shape of the time profile at the source is not known, so we choose the following analysis strategy. In general, the short pulse structure of any flare would be blurred by an energy-dependent effect on photon propagation. Conversely, one may correct for the effects of any given parametric model of photon dispersion, e.g., the linear or quadratic vacuum refractive index, by applying to each photon the appropriate time shift [15] corresponding to its propagation in a spatially-flat expanding uni- where H 0 is the Hubble expansion rate and h(z) = Ω Λ + Ω M (1 + z) 3 . If the correct energy-dependent QG shift is applied, the short pulse structure of the emission profile is restored. We implement this analysis strategy in two ways. In one analysis, we consider the most active part of the flare, that is distinguished clearly from the uniform background, and the QG shift is varied so as to maximize the total energy in this part.
In the other analysis we use the shape of the flare as extracted from untransformed low-energy data.

Energy cost function
It is well known [26] that a pulse of electromagnetic radiation propagating through a linearly-dispersive medium, as postulated here, becomes diluted so that its power (the energy per unit time) decreases. 5 Any transformation of a signal to reproduce the undispersed signal tends to recover the original power of the pulse. If the parameter M QGn is chosen correctly, the power of the recovered pulse is maximized. 5 The applicability of classical electrodynamics for estimating the low-energy behavior induced by space-time foam [1,2,[4][5][6][7][8][9]11] and the corresponding pulse-broadening effect have been discussed elsewhere (see [12] for details and an explicit example). The dilution effects for the linear or quadratic cases may easily be obtained as described in [ We implement this observation as follows. First, we choose a time interval (t 1 ; t 2 ) containing the most active part of the flare, as determined using a Kolmogorov-Smirnov (KS) statistic [27]. The KS statistic is calculated from the difference between the cumulative distribution function estimated from the unbinned data and that of a uniform distribution. The interval (t 1 ; t 2 ) covers the time range where the value of the KS difference varies from its maximum over the whole time support of the signal down to a negligible value. This procedure determines the proper time-width t 2 − t 1 of the most active (transient) part of the flare. 6 Having chosen this window, we scan over the whole support the time-distribution of all photons shifted by ∆t(E) and sum up the energies of photons in the window. For convenience, we re-parametrize the time shift as ∆t = ±τ l E or ∆t = ±τ q E 2 respectively, with τ l and τ q in s/GeV and s/GeV 2 units. The transformation is repeated for many values of τ l and τ q , chosen so that the shifts ∆t match the precision of the arrivaltime measurements, and for each τ l or τ q the scan is performed and the maximal summed energy in a window of width t 2 − t 1 is obtained. The maximal energies as a function of τ l or τ q define the 'energy cost function' (ECF). The position of the maximum of the ECF indicates the value of τ l or τ q that best recovers the signal, in the sense of maximizing its power. 7 This procedure is applied to 1000 Monte 6 The time interval chosen agrees very well with the spread of a Gaussian fit to the profile of the binned data, as well as the more complicated profile used in [23]. 7 Varying slightly the boundaries of the interval (t 1 ; t 2 ) has a negligible effect on the position of the maximum. We take into account the difference between the width at the Earth Carlo (MC) data samples generated by applying to the measured photon energies the (energy-dependent) Gaussian measurement errors. Fig. 1 shows the ECF for one such energy-smeared MC sample. It exhibits a clear maximum, whose position may be estimated by fitting it with a Gaussian profile in the peak vicinity. Fig. 2 shows the results of such fits to the ECFs with τ l for the 1000 energy-smeared realizations of the July 9 flare. From this distribution we derive the value τ l = (0.030±0.012) s/GeV, where M QG1 = 1.445×10 16 s/τ l , leading to a lower limit M QG1 > 0.21 × 10 18 GeV at the 95% C.L. 8 The same procedure applied to the ECF obtained using τ q leads to τ q = (3.71 ± 2.57) × 10 −6 s/GeV 2 , where M QG2 = 1.222 × 10 8 (s/τ q ) 1/2 , corresponding to M QG2 > 0.26 × 10 11 GeV at the 95% C.L. While our results for the June 30 flare have similar sensitivities and are compatible, they cannot be used to strengthen our results, as this flare is not very significant.

Likelihood function
We have confirmed this result using another technique to study the energy-dependent delay signal in the data. It is motivated by the initial time and energy-binned analysis performed in [23], which we used to check that the light-curve is well described by a simple Gaussian profile, superimposed on a time-independent background. We compute a likelihood function L based on the probability of a photon to be observed with energy E and arrival time t, using variables describing the energy spectrum at the source, the time distribution at emission obtained from the measured arrival times of the photons assuming an adjustable energy-dependent propagation delay, and the energy resolution of the detector, which is modelled as a Gaussian [28]. To describe the photon energy at the source a simple power law Γ(E s ) ∼ E −β s is taken, with β = 2.7 for the time-uniform part of the flare and 2.4 for the flaring part [23].
The likelihood function is fitted to the July 9 data minimizing − log L as a function of four parameters: (i) the energy-dependent delay parameterized in terms of M P /M QG1 , (ii) the position of the intrinsic maximum of the Gaussian flare, (iii) its width and (iv) the normalization of the time-independent background component in arbitrary units. The best fit yieldsM P /M QG1 = 8.2 +3.7 −3.4 , corresponding to τ l = (0.048 ± 0.021) s/GeV. The shape of the function χ 2 ≡ −2 log L + const around the minimum in these variables is quite parabolic almost up to the 2-σ level. In view of the correlations with these parameters, the sensitivity to τ l would be improved if they were known more precisely.
A similar procedure in the case of a quadratic dependency gives τ q = (4.60±5.46)× and at the source, also negligible.  Fig. 3. Comparison of the MAGIC measured lightcurve at low and high energies with the prediction given by the best set of parameters found using the likelihood method, and binning the data and the likelihood function in the same manner. 10 −6 s/GeV 2 . Fig. 3 shows that the L function gives a good overall fit to the data: binning in time and energy both the data and the L function, we find χ 2 /NDF ∼ 1.

Crosscheck with Monte Carlo data
To check the robustness of the ECF and likelihood analyses, we simulated several MC test samples with two components: (a) a time-independent background with the same energy spectrum as the measured data before the flare, and (b) a superposed signal generated at the source with an energy spectrum similar to that observed during the flare and an energy-independent Gaussian time distribution, each with the same numbers of photons as in the measurement. We then calculated the arrival times of all photons using various dispersion models and parameters, taking into account the MAGIC energy resolution. For each dispersion model and parameter, we generated 1000 incarnations, using different random seeds. These samples were then analyzed blindly, and the encoded effects were recovered successfully by the two estimators within the expected uncertainties. In addition, the analysis techniques were applied to MC samples with no energy-dependent dispersive sig-nal encoded, and found no effect, and both techniques also returned null results when applied to Mkn 501 data from outside a flare. These tests confirm the numerical sensitivities of the analyses and the estimates of the uncertainties given above. For the likelihood method, additional checks have been performed [28] assuming different flare energy spectra and shapes, besides the Gaussian one discussed here, which also fit reasonably well the binned data (c.f. Fig. 3).

Conclusions
The probability of the zero-delay assumption relative to the one obtained with the ECF estimator is P = 0.026. The observed energy-dependent delay thus is a likely observation, but does not constitute a statistically firm discovery. The results of the two independent analyses of the July 9 flare of Mkn 501 are quite consistent within the errors. Their results exhibit a delay between γ-rays of different energies, τ l = (0.030 ±0.012) s/GeV, corresponding to a lower limit M QG1 > 0.21 ×10 18 GeV at the 95% C.L. We also find a quadratic delay τ q = (3.71±2.57)×10 −6 s/GeV 2 , and M QG2 > 0.26 × 10 11 GeV at the 95% C.L., far beyond previous limits on a quadratic effect in photon propagation [11,14,18]. These numbers could turn into a real measurement of M QG1,2 , if the emission mechanism at the source were understood and the observed delays were mainly due to propagation. We cannot exclude, however, the possibility that the delay we find, which is significant beyond the 95% C.L., is due to some energy-dependent effect at the source. 9 However, we can exclude the possibility that the observed time delay may be due to a conventional QED plasma refraction effect induced as photons propagate through the source. This would induce [29] ∆t = D(α 2 T 2 /6q 2 ) ln 2 (qT/m 2 e ), where α is the fine-structure constant, q is the photon momentum, T is the plasma temperature, m e is the mass of electron, D is the size of the plasma, and we use natural units: c, = 1. Plausible numbers such as T ∼ 10 −2 MeV and D ∼ 10 9 km (for a review see [30]) yield a negligible effect for q ∼ 1 TeV. Exclusion of other source effects, such as time evolution in the mean emitted photon energy, might be possible with the observation of more flares, e.g., of different AGNs at varying redshifts. Observations of a single flare cannot distinguish the quantum-gravity scenarios considered here from modified synchrotronself-Compton mechanisms [23,31]. However, this pioneering study demonstrates clearly the potential scientific value of an analysis of multiple flares from different sources. The most promising candidate for applying the analyses proposed here is the flare from PKS 2155-304 detected recently by H.E.S.S. [32]. Unfortunately the occurrence of fast flares in AGNs is currently unpredictable, and since no correlation has yet been established with observations in other energy bands that could be used as a trigger signal, only serendipitous detections are currently possible.