First Combined Study on Lorentz Invariance Violation from Observations of Energy-dependent Time Delays from Multiple-type Gamma-Ray Sources. I. Motivation, Method Description, and Validation through Simulations of H.E.S.S., MAGIC, and VERITAS Data Sets

Modern physics is based on two fundamental pillars: quantum mechanics (QM) and Einsteinian general relativity (GR). When taken separately, these theories can claim success in satisfactorily describing many physical phenomena, but all attempts to make them compatible with each other have failed so far. The goal of quantum gravity (QG) research is to find a common approach to coherently merge quantum theory and GR. The QG problem has remained unsolved for more than 80 years now and keeps challenging physicists, who, in the struggle to find a solution, have proposed a myriad of models (see e.g., Polyakov 1981; Bombelli et al. 1987; Oriti 2001; Niedermaier & Reuter 2006; Rovelli 2007; Ambjørn et al. 2012). However, none of these models can claim full success. One of the main obstructions to progress in this field is the lack of experimental guidance. However, in the last two decades, the situation has changed, and in recent years there have been important advances in the field of QG phenomenology (Mattingly 2005; Amelino-Camelia 2013; Liberati 2013).

It is notoriously difficult to extract observable predictions from full-fledged QG approaches. Different models usually start from different conceptual premises and use different mathematical formalisms in such a way that it is difficult to determine whether they make compatible predictions. In some cases, the formal complexity forbids observable outcomes from being produced at all. Then, to guide experimental efforts, bottom-up approaches have been proposed (Amelino-Camelia 2002; Kowalski-Glikma & Nowak 2002; Magueijo & Smolin 2004; Livine et al. 2011; Barrau et al. 2015; Brahma & Ronco 2018; Calcagni et al. 2019). They rely on somewhat simpler models, suitable for describing only a subset of expected QG features but have the advantage of producing opportunities for experimental testing.

In this regard, at the end of the 1990s, independent semiclassical analyses inspired by QG models brought to the attention of the QG community the fact that it is a highly nontrivial task to retain Lorentz symmetries when quantizing the spacetime geometry of GR. These models include, most notably, string theory (see, e.g., Mavromatos 2010 and references therein), loop QG (Gambini & Pullin 1999), noncommutative geometry (Carroll et al. 2001), and standard model extension (Kostelecký & Mewes 2008 and references therein). From then on, departures from Lorentz invariance have become one of the rare observable features we would expect in a QG theory, and, as we shall see briefly, different bottom-up models to implement them have been proposed. According to this view, Lorentz invariance could be an emergent symmetry that arises in the low-energy limit but is modified at higher energies approaching the Planck scale, i.e., the energy scale at which both GR and QM effects should play an important role.

A much-studied way to encode departures from Lorentz invariance, either violations (noted LIV for Lorentz invariance violation) or deformations, consists of modifying the energy–momentum dispersion relation of free relativistic particles as follows (Amelino-Camelia et al. 1998):

$$\begin{eqnarray}&&{E}^{2}\simeq {p}^{2}{c}^{2}\times \left[1\pm \sum _{n=1}^{\infty }{\left(\displaystyle \frac{E}{{E}_{\mathrm{QG}}}\right)}^{n}\right],\end{eqnarray} \tag{ 1 }$$

where c is (the low-energy limit of) the speed of light, and E_QG the energy scale of QG effects, which is usually expected to be around the Planck scale (${E}_{P}=\sqrt{{\hslash }{c}^{5}/G}\,\simeq {10}^{19}$ GeV). The "±" sign in Equation (1) takes into account the possibility of having subluminal or superluminal effects.

Published one year after the first redshift of a gamma-ray burst (GRB) was measured, the article by Amelino-Camelia et al. (1998) also proposed for the first time the use of transient, distant, and high-energy gamma-ray sources as a way to probe the quantum nature of spacetime by searching for energy-dependent delays. In the following, we will focus on this particular way to probe a dispersion relation, such as the one in Equation (1), in so-called "time-of-flight" studies. Since then, other possibilities have emerged to search for QG effects in gamma-ray astronomy. For example, astrophysical sources have been used to search for vacuum birefringence (Götz et al. 2014) and spacetime "fuzziness" (Vasileiou et al. 2015). Possible modifications of the cross section of γ γ interaction between high-energy photons and the extragalactic background light were also investigated (Biteau & Williams 2015; Abdalla et al. 2019). Some of the limits published in these papers exceed the Planck scale, sometimes even by several orders of magnitude, but there is also a possibility that LIV could occur only through energy-dependent delays. In principle, all these effects could also coexist, even if they have only been tested separately so far. A comprehensive review of different possible effects of LIV on gamma-rays, as well as on other messengers (cosmic rays, gravitational waves, neutrinos), is given by Addazi et al. (2021).

Heuristically, Equation (1) can be justified as follows: at Planckian distances (∼10⁻³³ cm), QG effects are believed to cause fluctuations of spacetime geometry, which, then, would behave as a dynamical medium characterized by a nontrivial refractive index. Consequently, photons with different energies would have different interactions with the "foamy" structure of spacetime (sometimes called "quantum spacetime") and, thus, they would propagate in vacuum at different velocities, thereby producing an effect of in-vacuo dispersion. This explains the dependence of Equation (1) on some power n of the energy E of the probe. For simplicity, n is generally assumed to be an integer, and we will keep that assumption in this paper. However, in so-called fractional or (multi)fractional models, the modifications of the dispersion relation depend on noninteger powers of the energy (Amelino-Camelia et al. 2017; Calcagni 2017).

Regardless of the model to be used, the expected scale of QG effects is typically several orders of magnitude higher than the energy of observed photons. For this reason, we can treat the anomaly induced by QG as a small correction to the photon group velocity and, in particular, only linear n = 1 or quadratic n = 2 modifications are of interest for experimental searches, taking into account the sensitivity of current detectors. It is important to stress that there are counterexamples where E_QG can be far from the Planck scale (being either above or below E_P ). Among others, we highlight two particular cases. In the approach of asymptotic safety, renormalization group techniques generate a running of the gravitational constant, thereby affecting the value of E_QG (Niedermaier & Reuter 2006). In string theory, the compactification of extra dimensions can produce testable effects at energies much lower than E_P , even of the order of tens of teraelectronvolts (TeV; Antoniadis et al. 1998). Some stringent constraints already exist on these models (Aad et al. 2016). Given that, different types of experiments play a crucial role in constraining the value of E_QG.

To compensate for the smallness of the effect (E/E_QG is typically of the order 10⁻¹⁹–10⁻¹⁴), it has been recognized that very distant astrophysical sources could be used to probe the properties of quantum spacetime (Amelino-Camelia et al. 1998; Alfaro 1999; Jacobson et al. 2006). Indeed, if the emitted photons travel over large distances, then even extremely tiny quantum-spacetime effects could accumulate, and eventually, the overall effect could become detectable in the form of energy-dependent time delays in the light curves. Variable or transient sources at cosmological distances such as GRBs and flaring active galactic nuclei (AGNs) are good candidates for looking for LIV, but it is important to stress that the involvement of cosmological distances forces us to face the problem of combining curvature with quantum-spacetime effects. In other words, the delays should depend on the redshift. On the other hand, fast-spinning pulsars (PSRs) detected at TeV energies are within our galaxy and, thus, their Euclidean distances can be used instead of the redshift.

Considering only the leading dominant term in Equation (1), either linear (n = 1) or quadratic (n = 2), it can be shown that the group velocity of photons acquires a dependence on their energies. In particular, the delay between two photons emitted at the same time by a source at redshift z with energies E_h > E_l is

$$\begin{eqnarray}&&{\rm{\Delta }}{t}_{n}\simeq \pm \,\displaystyle \frac{n+1}{2}\,\displaystyle \frac{{E}_{h}^{n}-{E}_{l}^{n}}{{{\rm{H}}}_{0}{E}_{\mathrm{QG}}^{n}}\ {\kappa }_{n}(z),\end{eqnarray} \tag{ 2 }$$

where κ_n (z) is a parameter depending on the distance of the source. The symbol ± allows us to take into account both a subluminal (sign +) or a superluminal (sign −) LIV effect. In this paper, two different expressions for κ_n (z) will be compared for the first time: one obtained in a pure LIV framework (Jacob & Piran 2008) and another obtained in the doubly (or deformed) special relativity (DSR) approach (Rosati et al. 2015). This will be discussed in more detail in Section 2.

The delay Δt_n only takes into account Lorentz violation effects, therefore ignoring any time lag originating from emission mechanisms, also referred to as "source-intrinsic" delay. A hint of such kind of delay is observed for GRBs (Ajello et al. 2019) and has been also reported once in the case of an AGN, for the flare of Mrk 501 in 2005 recorded by MAGIC ¹² (Albert et al. 2007). With only one source and with only a rough knowledge of how particles are emitted and accelerated, intrinsic delays cannot be separated from propagation effects. Modeling of astrophysical sources is an ongoing effort and a first study of source-intrinsic effects in connection with Lorentz violation searches has been published recently in the case of blazar flares (Perennes et al. 2020; Levy et al. 2022). On the other hand, when several sources are combined, it could be possible, at least in principle, to separate intrinsic and propagation effects. Indeed, it is reasonable to assume that intrinsic delays do not depend on the distance. It is therefore essential that these studies could be performed on a large population of objects.

From Equation (2), another parameter λ_n can be defined as

$$\begin{eqnarray}&&{\lambda }_{n}\equiv \displaystyle \frac{{\rm{\Delta }}{t}_{n}}{{\rm{\Delta }}{E}_{n}\ {\kappa }_{n}(z)}=\pm \displaystyle \frac{n+1}{2{H}_{0}\ {E}_{\mathrm{QG}}^{n}}\end{eqnarray} \tag{ 3 }$$

using the simplified notation ${\rm{\Delta }}{E}_{n}\equiv {E}_{h}^{n}-{E}_{l}^{n}$. This parameter λ_n , which will be used later, has the advantage of being independent of the distance of the source and is therefore suitable for multisource analysis.

Since the late 1990s, the field has rapidly expanded, with more and more sources being analyzed in the search for LIV effects. With the notable exception of the flare of Mrk 501 in 2005 already mentioned above, no significant delay has been reported so far when using only photons as the messenger. Constraints have been improving on a regular basis, even reaching the Planck scale in some cases in analyses of individual objects (see, e.g., Vasileiou et al. 2013). Table 1 gives a partial selection of the best limits available on E_QG,n for time-of-flight studies, with the three types of sources (AGNs, GRBs, and PSRs). This new notation E_QG,n reflects the fact that LIV analyses have different sensitivities for linear and quadratic effects.

Table 1. A Selection of Limits for a Subluminal Propagation Obtained with Various Instruments and Various Types of Objects

Source	Experiment	Year	Distance a	Lower Limit on EQG,1	Lower Limit on EQG,2	Reference	Note
				(95% CL, GeV)	(95% CL, GeV)
35 GRB	BATSE, HETE-2, Swift	⋯	⋯	1.4 × 1016	⋯	1	b
8 GRB	Fermi LAT	⋯	⋯	1.0 × 1017	⋯	2
GRB 090510	Fermi LAT	2009	0.903	9.3 × 1019	1.3 × 1011	3
GRB 190114C	MAGIC	2019	0.4245	0.6 × 1019	6.3 × 1010	4
Mrk 501	MAGIC	2005	0.034	0.3 × 1018	5.7 × 1010	5	c , e
Mrk 501	H.E.S.S.	2014	0.034	3.6 × 1017	8.5 × 1010	6
PKS 2155–304	H.E.S.S.	2006	0.116	2.1 × 1018	6.4 × 1010	7	e
PG 1553+113	H.E.S.S.	2012	0.49 ± 0.04	4.1 × 1017	2.1 × 1010	8	d , e
PSR B0531+21	VERITAS	2007-14	2.2 kpc	1.9 × 1017	⋯	9
PSR B0531+21	MAGIC	⋯	2.2 kpc	5.5 × 1017	5.9 × 1010	10	e
PSR B0833–45	H.E.S.S.	⋯	294 pc	4.0 × 1015	⋯	11	e

Notes.

^a Redshift is given for extragalactic objects. ^b The limits of Ellis et al. (2006) were corrected in Ellis et al. (2008) taking into account the factor $(1+z^{\prime} )$ in the numerator of the integral in Equation (4). Only the limit obtained for a linear correction is given. ^c These numbers are actually reported as best-fit values by Martínez & Errando (2009). ^d The redshift of this source was not measured but only estimated. ^e Sources used as benchmark in the present paper.

References. (1) Ellis et al. (2006, 2008), (2) Ellis et al. (2019), (3) Vasileiou et al. (2013), (4) Acciari et al. (2020), (5) Martínez & Errando (2009), (6) Abdalla et al. (2019), (7) Abramowski et al. (2011), (8) Abramowski et al. (2015), (9) Zitzer et al. (2013), (10) Ahnen et al. (2017), (11) Chretien et al. (2015).

Download table as:  ASCIITypeset image

The results of Ellis et al. (2006, 2008, 2019) are of particular interest because they were obtained from the analyses of several GRBs. This kind of analysis, repeated with different experiments (e.g., Bolmont et al. 2008; Bernardini et al. 2017) consists of two steps: first, the time lags are computed for each individual source and then the obtained data points are fitted with a function Δt = a z + b(1 + z) (for n = 1). The value of parameter a is subsequently used to constrain E_QG,1 while b represents source-intrinsic effects, assumed to be identical for all bursts. In the present paper, we describe and test a more advanced and presumably more sensitive method to perform such a population study, based on a likelihood technique.

For completeness, let us mention that recently it was suggested from the analysis of several GRBs that Lorentz invariance could be violated at a scale of ∼3.7 × 10¹⁷ GeV (Xu & Ma 2016a, 2016b, 2018). This result is contradictory with the best limits listed in Table 1 and still needs to be confirmed.

One of the main objectives of the first phase of describing QG phenomenology was the ability to prove that in-vacuo dispersion (or, more generally, Planck-scale effects) could be tested with current experiments. Now, the stringent limits established with GRB observations, together with the growing amount of relevant data and progress on the theory side, can bring us to a more mature phase where we can start constraining actual QG models in a robust manner. The present work can be considered a first step in this direction. We aim to combine, for the first time, the data obtained with the three major imaging atmospheric Cherenkov telescope (IACT) experiments, H.E.S.S., ¹³ MAGIC, and VERITAS ¹⁴ in order to constrain QG effects through the time-of-flight technique (see Terzić et al. 2021 for a recent review). This combination will extract the most information out of each type of source to produce robust constraints from existing data, while taking into account the redshift dependence of the LIV-induced time lag.

The paper will be divided in two parts. In the present article (part I), two possible ways to account for the dependence of time delays as a function of redshift will first be described (Section 2). Then, the method used to compute and combine the likelihoods to measure time-lag parameters λ_n and τ_n will be described in detail (Section 3). Several nuisance parameters are included in the computation to take into account various sources of systematic uncertainties. Then, in Section 4, the method is tested on simulated data sets, mimicking the data of several representative sources observed in the TeV domain. These simulations are used to evaluate statistical errors and study the impact of various sources of systematic errors in the lag measurement. The results, as well as the impact of redshift dependence, will be given and discussed in Section 5.

In the second part of the paper, to appear later, the method will be used with available data from H.E.S.S., MAGIC, and VERITAS, and possibly from other gamma-ray experiments, in order to produce a combined limit on E_QG,n .

It is rather natural to believe that curvature and quantum effects are deeply intertwined because curvature is a key characteristic of spacetime geometry. In light of this, a complete QG theory would be needed to tell us whether there is a phenomenon of in-vacuo dispersion and then compute its magnitude. However, in the absence of such a theory, simplified speculative approaches to model in-vacuo dispersion in curved spaces have been proposed. Among them, especially for reasons of simplicity, a model where Lorentz invariance is explicitly broken in a specific way proposed by Jacob and Piran (Jacob & Piran 2008, J&P for short) attracted particular interest and has been systematically used so far in experimental analyses constraining in-vacuo dispersion.

In this approach, parameter κ_n (z) is expressed as

$$\begin{eqnarray}{\kappa }_{n}^{{\rm{J}}\&{\rm{P}}}(z)\equiv {\int }_{0}^{z}\displaystyle \frac{{\left(1+z^{\prime} \right)}^{n}}{\sqrt{{{\rm{\Omega }}}_{m}{\left(1+z^{\prime} \right)}^{3}+{{\rm{\Omega }}}_{{\rm{\Lambda }}}}}\ {dz}^{\prime} ,\end{eqnarray} \tag{ 4 }$$

where the denominator relates to the Hubble parameter $H(z)={{\rm{H}}}_{0}\sqrt{{{\rm{\Omega }}}_{m}{\left(1+z\right)}^{3}+{{\rm{\Omega }}}_{{\rm{\Lambda }}}}$. In the following, values of cosmological parameters are taken from Planck results (Aghanim et al. 2020) as recommended by the Particle Data Group (Zyla et al. 2020): H₀ = 67.4 ± 0.5 km s⁻¹ Mpc⁻¹, Ω_m = 0.315 ± 0.007, and Ω_Λ = 0.685 ± 0.007.

As has been pointed out in recent literature (Rosati et al. 2015; Barcaroli et al. 2016; Pfeifer 2018), Equation (4) offers only one possible parameterization among many others. It has been shown that if, following the DSR approach, Poincaré symmetries are modified in order to preserve the invariance of Equation (1) under relativistic transformations, then one can obtain a different result for the distance parameter (Rosati et al. 2015):

$$\begin{eqnarray}&&{\kappa }_{n}^{\mathrm{DSR}}(z)\equiv {\int }_{0}^{z}\displaystyle \frac{{h}^{2n}(z^{\prime} ){dz}^{\prime} }{{\left(1+z^{\prime} \right)}^{n}\,\sqrt{{{\rm{\Omega }}}_{m}{\left(1+z^{\prime} \right)}^{3}+{{\rm{\Omega }}}_{{\rm{\Lambda }}}}},\end{eqnarray} \tag{ 5 }$$

with

$$\begin{eqnarray}\begin{array}{rcl}h(z^{\prime} ) & \equiv & 1+z^{\prime} -\sqrt{{{\rm{\Omega }}}_{m}{\left(1+z^{\prime} \right)}^{3}+{{\rm{\Omega }}}_{{\rm{\Lambda }}}}\\ & & \times \,{\int }_{0}^{z^{\prime} }\displaystyle \frac{{dz}^{\prime\prime} }{\sqrt{{{\rm{\Omega }}}_{m}{\left(1+z^{\prime\prime} \right)}^{3}+{{\rm{\Omega }}}_{{\rm{\Lambda }}}}}\,,\end{array}\end{eqnarray} \tag{ 6 }$$

and this can result in consistently different limits on E_QG,n . Note that Equation (5) is only one possible outcome of DSR, chosen here as a benchmark. Using observations of multiple sources at different redshifts, we establish for the first time limits on these two different models.

Figure 1 shows the functions κ^DSR and κ^J&P as a function of redshift for n = 1 and n = 2. κ^DSR is smaller than κ^J&P for both linear and quadratic cases. When no lag is measured, this leads to less stringent limits on E_QG,n . Both functions κ^J&P and κ^DSR increase with redshift, thus increasing the expected time delay. However, it has to be pointed out that ultimately, the distance at which sources can be detected at high energies is limited by the absorption by the extragalactic background light (EBL). This distance depends on the energy and does not exceed z ∼ 1 in the TeV range.

Zoom In Zoom Out Reset image size

To conclude this section, let us add that in the case of nearby sources, such as pulsars, the Euclidean approximation is valid, i.e., κ_n (z) = d H₀/c, where d is the Euclidean distance to the source. In addition, the ratio κ^DSR/κ^J&P converges to unity for low distances. As a result, a given pulsar will give the same constraints on E_QG,n for both J&P and DSR cases.

All observations considered for combination in this work are analyzed with a maximum likelihood (ML) method to search for linear or quadratic LIV delays and to extract limits on E_QG,n . Compared to alternative methods such as e.g., PairView (Vasileiou et al. 2013), dispersion correction (Barres de Almeida & Daniel 2012), and peak comparison (Ahnen et al. 2017), the ML method allows optimal use of the information in data and provides a relatively straightforward way to combine analyses of multiple sources and observatories. On the other hand, it relies on parameterization of intrinsic photon emission time and energy distributions, which are currently not fully understood at the theoretical level. The uncertainties related to these parameterizations are taken into account when deriving the limits on E_QG,n .

The code for likelihood computation as well as for simulations was developed using the ROOT ¹⁵ framework (Brun & Rademakers 1997).

3.1. Single-source Likelihood

First applied by Martínez & Errando (2009) to analyze the 2005 flare of Mrk 501 observed by MAGIC, the ML method relies on defining a probability density function (PDF) that describes the probability of observing a gamma-ray photon at a certain arrival time and with a certain reconstructed energy.

In its simplest form, the PDF for signal photons is defined as a function of time and energy with λ_n as the single parameter to be estimated. The PDF is obtained by convolving the spectrum of the source Γ_s (E_t ) with the light curve C_s , both as observed on Earth, i.e., after propagation:

$$\begin{eqnarray}&&{F}_{s}({E}_{t},t;{\lambda }_{n})=\displaystyle \frac{{{\rm{\Gamma }}}_{s}({E}_{t}){C}_{s}\left(t-D({E}_{t},{\lambda }_{n},z)\right)}{{N}_{s}},\end{eqnarray} \tag{ 7 }$$

where E_t is the true energy of the gamma-ray photon, t is the arrival time corrected by the factor

$$\begin{eqnarray}&&D({E}_{t},{\lambda }_{n},z)={\lambda }_{n}\times {\kappa }_{n}(z)\times {E}_{t}^{n},\end{eqnarray} \tag{ 8 }$$

which defines the propagation delay due to LIV, where λ_n is given by Equation (3), and N_s is a normalization term expressed as follows:

$$\begin{eqnarray}&&{N}_{s}=\iint {{\rm{\Gamma }}}_{s}({E}_{t}){C}_{s}\left(t-D({E}_{t},{\lambda }_{n},z)\right){{dE}}_{t}\,{dt}.\end{eqnarray} \tag{ 9 }$$

In Equations (7) and (9), function C_s is often called the template light curve. It is usually obtained by fitting a light curve at low energies, where LIV effects are assumed to be weak or negligible. Because there is no fully accepted model available that reproduces the shape of the light curves for GRBs, AGNs, or PSRs, Gaussian or Lorentzian functions, or the sum of several of these functions, are usually used. The function Γ_s is obtained from the data on the full energy range considered for the LIV analysis (see Section 4).

Background events of several different origins are also taken into account. They include hadrons misreconstructed as gamma-rays and baseline photons emitted either by an AGN in its quiescent state or by the nebula surrounding a PSR. The PDF for background events of type k (hadrons or baseline photons), which are not affected by LIV propagation effects, is written as

$$\begin{eqnarray}&&{F}_{b,k}({E}_{t},t)=\displaystyle \frac{{{\rm{\Gamma }}}_{b,k}({E}_{t}){C}_{b,k}}{{N}_{b,k}}.\end{eqnarray} \tag{ 10 }$$

Γ_b,k is the background spectrum taken as a power law. For hadrons, the index is set to 2.7 while the values for signal and baseline photons are given in Table 3. C_b,k is the time distribution of background events assumed to be a constant, and N_b,k the normalization term defined as

$$\begin{eqnarray}&&{N}_{b,k}=\iint {{\rm{\Gamma }}}_{b,k}({E}_{t}){C}_{b,k}\,{{dE}}_{t}{dt}.\end{eqnarray} \tag{ 11 }$$

From Equations (7) and (10), the complete definition of the PDF is obtained, accounting for detector performance assessed from instrument response functions (IRFs):

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{{dP}}{{{dE}}_{m}{dt}}={w}_{s}\ \displaystyle \frac{\int A({E}_{t},\vec{\varepsilon })M({E}_{t},{E}_{m})\times {F}_{s}({E}_{t},t;{\lambda }_{n}){{dE}}_{t}}{{N}_{s}^{{\prime} }}\\ +\,\displaystyle \sum _{k}{w}_{b,k}\ \displaystyle \frac{\int A({E}_{t},\vec{\varepsilon })M({E}_{t},{E}_{m})\times {F}_{b,k}({E}_{t},t){{dE}}_{t}}{{N}_{b,k}^{{\prime} }},\end{array}\end{eqnarray} \tag{ 12 }$$

where the source and background terms F_s and F_b,k are convoluted with the detector effective area $A({E}_{t},\vec{\varepsilon })$ and energy resolution M(E_t , E_m ). Source and background terms are weighted by w_s and w_b,k , respectively, with w_s + ∑_k w_b,k = 1. E_t still denotes the true energy while E_m is the corresponding measured energy. The parameters ${N}_{s}^{{\prime} }$ and ${N}_{b,k}^{{\prime} }$ are the normalization factors of the PDF. In addition to energy E_t , the effective area depends on a set of factors $\vec{\varepsilon }$, which vary with observation conditions and with the method used for event reconstruction and identification. The IRFs were kindly provided by the H.E.S.S., MAGIC, and VERITAS collaborations. Distinct IRFs are used for each source and each observation period.

The confidence levels (CLs) for either a measurement or the derivation of lower limits on λ_n (and E_QG,n ) can then be obtained summing the log-likelihood of all the events for a given source S:

$$\begin{eqnarray}&&{L}_{S}({\lambda }_{n})=-\sum _{{\rm{i}}}\mathrm{log}\left(\displaystyle \frac{{dP}}{{{dE}}_{m}{dt}}({E}_{m,i},{t}_{i});{\lambda }_{n}\right).\end{eqnarray} \tag{ 13 }$$

3.2. Combining Likelihoods

While each source may require a different analysis strategy, either using a single-parameter likelihood or a profile likelihood, the combination of multiple sources is straightforward. Once log-likelihood functions L_S (λ_n ) are obtained for all sources, the combined log-likelihood L_comb is simply given by their sum:

$$\begin{eqnarray}&&{L}_{\mathrm{comb}}({\lambda }_{n})=\sum _{\mathrm{all}\,\mathrm{sources}}{L}_{S}({\lambda }_{n}).\end{eqnarray} \tag{ 14 }$$

3.3. Statistical and Systematic Uncertainties

Statistical and systematic uncertainties are propagated in the final result through the use of the profile likelihood. The log-likelihood for each source is then written as

$$\begin{eqnarray}\begin{array}{rcl}L({\lambda }_{n},\vec{\theta }) & = & {L}_{{\rm{S}}}({\lambda }_{n},\vec{\theta })+{L}_{\mathrm{template}}({\vec{\theta }}_{{\rm{C}}})+{L}_{\gamma }({\theta }_{\gamma })\\ & & +\,{L}_{{\rm{B}}}({\vec{\theta }}_{{\rm{B}}})+{L}_{\mathrm{ES}}({\theta }_{\mathrm{ES}})+{L}_{{\rm{z}}}({\theta }_{{\rm{z}}}),\end{array}\end{eqnarray} \tag{ 15 }$$

where $\vec{\theta }$ is the vector of all nuisance parameters defined as:

• 1.  
  ${\vec{\theta }}_{{\rm{C}}}$, the parameters of the light-curve analytic parameterization,
• 2.  
  θ_γ , the power-law index of the signal events spectrum,
• 3.  
  ${\vec{\theta }}_{{\rm{B}}}$, the ratio of signal and of background event numbers to the total number of events,
• 4.  
  θ_ES, the energy scale,
• 5.  
  θ_z, the distance or redshift.

As already mentioned above, the template light curve C_s of Equation (7) is obtained by fitting a low-energy light curve, for which LIV is assumed to be negligible. From this parameterization, it is possible to evaluate errors directly, defining ${L}_{\mathrm{template}}({\vec{\theta }}_{{\rm{C}}})$ as the sum of the log-likelihoods of each event generated from the low-energy template parameterization:

$$\begin{eqnarray}&&{L}_{\mathrm{template}}({\vec{\theta }}_{{\rm{C}}})=-\sum _{i=1}^{{N}_{\mathrm{template}}}\mathrm{log}\left(\displaystyle \frac{{C}_{s}({t}_{i},{\vec{\theta }}_{{\rm{C}}})}{{N}_{c}}\right),\end{eqnarray} \tag{ 16 }$$

with C_s the light curve and N_c its normalization. In this equation, the new notation for the light curve ${C}_{s}({t}_{i},{\vec{\theta }}_{{\rm{C}}})$ denotes the fact the template is evaluated for zero lag (D(E_i , λ_n = 0, z) = 0) and explicitly shows the parameter vector ${\vec{\theta }}_{{\rm{C}}}$ of the template function. On the other hand, some other analyses use the template fit results as nuisance parameters. In that case, ${L}_{\mathrm{template}}({\vec{\theta }}_{{\rm{C}}})=0$ and the uncertainty on ${\vec{\theta }}_{{\rm{C}}}$ is then accounted for in the generated data sample log-likelihood ${L}_{{\rm{S}}}(\lambda ,\vec{\theta })$ defined in Equation (13).

L_γ (θ_γ ) is obtained from the statistical and systematical uncertainties of the spectral index as provided in the analyses of the different sources. The flux normalization and energy-scale uncertainties provided by the different observatories are taken into account by ${L}_{{\rm{B}}}({\vec{\theta }}_{{\rm{B}}})$ and L_ES(θ_ES), respectively. The energy-scale parameter is introduced in the data sample log-likelihood by a scale factor applied to the event energy. The uncertainties on redshift for extragalactic sources, or distance for galactic sources, are accounted for in L_z(θ_z).

For the power-law index, the ratio of signal and of background, energy-scale, and redshift uncertainties, a normal distribution is assumed, which allows to use a simple chi-square approach:

$$\begin{eqnarray}&&{L}_{{\rm{x}}}({\vec{\theta }}_{{\rm{x}}})=\sum _{i}\displaystyle \frac{{\left({\theta }_{{\rm{x}},i}-{\bar{\theta }}_{{\rm{x}},i}\right)}^{2}}{2{\sigma }_{{\theta }_{{\rm{x}},i}}^{2}},\end{eqnarray} \tag{ 17 }$$

where ${\sigma }_{\theta }^{2}$ is the uncertainty of the nuisance parameter θ and x is a different type of systematics. The full list of uncertainties assigned to each nuisance parameter for each source is shown in Table 2.

Table 2. Nuisance Parameter Uncertainties for the Individual Sources

Source	Energy Scale	Background Proportion	Spectral Index a	Distance/Redshift	References b
GRB 190114C	17%	11%	0.21	Δz = 1 × 10−3	1, 2, 3, 4
PG 1553+113	10%	20%	0.31	Δz = 4 × 10−2	5, 6
Mrk 501	17%	11%	0.04	Δz = 1 × 10−4	1, 2, 3, 6
PKS 2155–304	10%	20%	0.1	Δz = 1.7 × 10−2	5, 7
Crab (M)	17%	11%	0.07	Δd = 506 pc	1, 2, 3, 8
Crab (V)	20%	22%	0.5	Δd = 506 pc	9, 8
Vela	10%	20%	0.67	Δd = 76 pc	5, 10

Notes.

^a Uncertainty for the spectral index includes both statistical and systematic errors. ^b References are in the same order as the columns: energy scale, background proportion, spectral index (the three values are sometimes given in the same reference), and distance.

References. (1) Aleksić et al. (2012), (2) Aleksić et al. (2016), (3) Aleksić et al. (2015), (4) Acciari et al. (2019b), (5) Aharonian et al. (2006), (6) Mao (2011), (7) Ganguly et al. (2013), (8) Kaplan et al. (2008), (9) for the energy scale and spectral index, Pueschel (2019). The value for the background proportion was provided by the VERITAS collaboration. (10) Caraveo et al. (2001).

Download table as:  ASCIITypeset image

In order to illustrate the impact of the different sources of uncertainties, the uncertainty on lambda is derived by varying only one of the different nuisance parameters. The systematic errors are then derived assuming the total uncertainty is the squared sum of the statistical and systematical uncertainties. They are presented in Appendix A with Table A1 for the J&P case and Table A2 for the DSR case, for each source and each source combination. These results will be commented on further in Section 5.

4.1. Simulated Data Sets

4.2. Method Calibration and Performance

The normalization factor ${N}_{s}^{{\prime} }$ of the PDF of Equation (12) is a triple integral, the computation of which is particularly time consuming because it needs to be done for each minimization step and for each event of the sample. To decrease the computation time, the PDF is precalculated and stored in tables binned over the measured energy E_m , nondelayed arrival times t, and time delays D(E_t , λ_n , z). The same number of bins is used for each of these three variables. A trilinear interpolation is performed on these tables to extract PDF values for the likelihood computation.

The number of bins used in the tables has been chosen for each source to minimize the bias λ_rec − λ_inj between the injected (λ_inj) and reconstructed (λ_rec) time delays. An example is shown in Figure 2 for GRB 190114C. In this particular case, the plot shows that a minimum of ∼140 bins for each variable is required, and a conservative number of 200 was actually chosen. The optimal number of bins was found to be independent of the injected lag. Four sets of tables have been produced for each source, which accounts for the four configurations explored in this work: J&P or DSR formalism for distance and for linear and quadratic LIV effects.

Zoom In Zoom Out Reset image size

In order to assess the sensitivity and precision of the lag reconstruction, simulated data sets were produced with different values for λ_inj. For each value of the injected lag, 1000 realizations of the light curve were simulated. The distribution of reconstructed values of λ_n is shown in the central panel of Figure 3 for GRB 190114C, in the J&P case, and n = 1. The lower and upper limits of the confidence intervals are taken as the values of λ_n for which $2\,[{L}_{S}({\lambda }_{n})-\min ({L}_{S})]=1$ for the 68% CL and $2\,[{L}_{S}({\lambda }_{n})-\min ({L}_{S})]=3.84$ for the 95% CL (Equation (13)). Their distributions for the 68% CL are displayed in the left and right panels of Figure 3, respectively. All three distributions were fitted with asymmetric Gaussian functions providing three parameters: the average (λ_LL, λ_rec, and λ_UL) and standard deviations separately defined on the left and on the right of the maxima (σ_λ,l , σ_λ,r ). While the latter accounts for statistical uncertainties only, the lower and upper limits λ_LL and λ_UL account for both statistical and systematic uncertainties.

Zoom In Zoom Out Reset image size

For extragalactic sources, the range for λ_inj goes from −5σ₀ to +5σ₀, where ${\sigma }_{0}=\max ({\sigma }_{{\lambda }_{\mathrm{rec}},l},{\sigma }_{{\lambda }_{\mathrm{rec}},r})$ for λ_inj = 0. In the case of PSR, the range for λ_inj is chosen so that the highest-energy photons are not shifted out of the phase range given in Table 3.

The plots for λ_rec versus λ_inj were then produced for individual sources, as well as for combinations, for the two correction orders and the two lag-distance models. Figure 4 shows two examples of calibration plots for GRB 190114C alone (left) and for all sources combined (right) in the linear and J&P case. The reconstructed lag is fitted with a linear function λ_rec = a λ_inj + b. The plot for the GRB alone shows a clear decrease of the reconstruction error as the injected lag increases. This is a consequence of a peculiar shape of the light curve, which has a narrow peak, followed by a power-law decay. As the value of λ_inj increases, the peak of the light curve enters progressively the time window where the likelihood is computed, resulting in an improvement of the reconstruction precision. The plot on the right shows the same behavior, illustrating how the GRB dominates over the other sources. Other examples of calibration plots are shown in Appendix B. As none of them include the GRB, the consistency in the reconstruction error is maintained. All of the plots produced show a very good reconstruction of the injected lag, with slopes a very close to unity, while the bias b is found to be close to zero.

Zoom In Zoom Out Reset image size

5.1. Systematic Uncertainties

5.2. Limits

From Equation (3), limits on E_QG are given by

$$\begin{eqnarray}&&{\left[\displaystyle \frac{2}{n+1}\left({\lambda }_{n,\pm }+\sqrt{{\delta }_{\mathrm{stat}}^{2}+{\sigma }^{2}{\delta }_{\mathrm{syst}}^{2}}\right){{\rm{H}}}_{0}\right]}^{\tfrac{1}{n}},\end{eqnarray} \tag{ 18 }$$

where the subscript ± refers to subluminal and superluminal cases, δ_stat is the statistical error (standard deviation) on the normally distributed reconstructed value of λ_n , δ_syst is the overall systematic error obtained from the values listed in Tables A1 and A2 computed for a CL of 68%, and σ is a real number allowing for a shift in CL using the same systematic errors. Because systematic errors are computed for 68% CL and statistical errors for 95% CL, σ is set to 2 in the following.

Using Equation (18), E_QG limits were obtained for both subluminal and superluminal cases. Both approaches give comparable results, and only the subluminal limits are shown in Table 4. In Figure 5, results are given with and without accounting for systematic uncertainties, clearly demonstrating the importance of taking them into account. In some cases, systematics lead to upper limits smaller by a factor of ≳ 2 as compared to the case where they are not taken into account. For pulsars, DSR and J&P formalisms lead to the same limits so they are given only for the J&P case in the table.

Zoom In Zoom Out Reset image size

Table 4. 95% CL Limits Obtained for Individual Objects and Combinations

Source	EQG,1				EQG,2
	J&P		DSR		J&P		DSR
	(1018 GeV)		(1018 GeV)		(1010 GeV)		(1010 GeV)
	w/o syst.	w/ syst.	w/o syst.	w/ syst.	w/o syst.	w/ syst.	w/o syst.	w/ syst.
GRB 190114C	9.2	4.0	6.5	2.7	14.2	8.3	9.5	5.8
PKS 2155–304	2.8	1.0	2.6	0.9	8.2	6.2	7.2	5.5
Mrk 501	1.1	0.5	1.1	0.5	9.6	7.1	9.3	6.9
PG 1553+113	0.17	0.11	0.10	0.07	1.3	1.0	0.87	0.68
Crab (M)	0.80	0.65	⋯	⋯	3.0	2.5	⋯	⋯
Crab (V)	0.48	0.10	⋯	⋯	1.5	0.94	⋯	⋯
Vela	5.1 × 10−3	3.5 × 10−3	⋯	⋯	5.6 × 10−2	5.5 × 10−2	⋯	⋯
Crab (M+V)	1.0	0.28	⋯	⋯	3.3	2.6	⋯	⋯
PSR	1.0	0.28	⋯	⋯	3.3	2.8	⋯	⋯
AGN	3.0	1.1	2.8	1.0	10.8	8.3	10.5	7.9
AGN+PSR	3.2	1.2	3.0	1.1	10.6	8.5	10.1	8.3
GRB+PSR	9.2	4.1	6.6	2.8	14.3	9.2	9.1	7.0
GRB+AGN	9.5	4.1	6.9	3.0	14.5	9.7	11.4	8.2
All combined	9.5	4.1	7.0	2.9	14.4	9.7	11.1	8.4

Download table as:  ASCIITypeset image

Figure 6 shows a comparison between the already published results taken from the references listed in Table 1 and the ones obtained in the present study, for n = 1 (left) and n = 2 (right). Overall, the agreement between simulations and data is good, showing the simulated data sets represent the actual data well. The observed differences are most probably due to three different factors. First, the limits obtained in this work come from several hundreds realizations of the light curves while already published limits were derived from one (measured) light curve. Second, the systematic uncertainties in previous publications were evaluated with different methods. These methods can vary from one analysis to another, but most use a frequentist approach, while in the present article, nuisance parameters and the profile likelihood where used (Section 3.3). Finally, IRFs were fully taken into account in the present analysis while it was often approximated as a constant of energy in earlier articles. The latter point was fully justified at the time by the use of a somewhat reduced energy range, while we wanted to get rid of this restriction in the present analysis.

Zoom In Zoom Out Reset image size

As expected from previously published results, GRB 190114C is the most constraining source due to its high redshift, high variability, and large statistics, as well as the fact it has been observed on a wide energy range. Therefore, it dominates the final result whenever it is included in the combination.

When the GRB observation is not included, AGNs dominate, with a competition between PKS 2155–304 and Mrk 501. Due to its smaller number of events and limited energy range, the PG 1553+113 limit is less constraining, even though its redshift is the highest of all the sources included in this work. While PKS 2155–304 dominates over the other sources in the linear case due to its higher redshift and event statistics, Mrk 501 dominates the limit in the quadratic case due to its energy range extending twice as high as the one of PKS 2155–304.

PSR have only a marginal impact on the overall combination due to their closeness. The Crab pulsar dominates the combined PSR limit thanks to its higher statistics, wider energy range, and greater distance. However, it is important to note that the limits provided by pulsars are independent of the redshift-dependence model, providing model-free constraints.

Figure 7 shows the limits on E_QG,n as a function of redshift for both DSR and J&P models. Differences in the results from the two approaches start to become significant for high-redshift sources such as GRB 190114C or PG 1553+113, which start behaving in accordance with the κ_n parameter evolution shown in Figure 1. Due to the fact that κ^J&P > κ^DSR and κ^J&P increases faster than κ^DSR (Section 2), the J&P model emphasizes the impact of high-redshift sources on the limits. Therefore, the GRB dominates more in the J&P case than in the DSR case, where all source contributions are more balanced.

Zoom In Zoom Out Reset image size

In the present paper, we have described an implementation of likelihood analysis designed with the goal to combine data from different sources and experiments in the search for LIV-induced energy-dependent time delays. One of the most important benefits of the likelihood technique is its simplicity for such a combination. In order to check the method and evaluate its performance, simulated data sets were produced mimicking actual observations of one GRB, three flaring AGNs, and two pulsars by the H.E.S.S., MAGIC, and VERITAS experiments. We paid particular attention to the implementation of the algorithm, checking for any bias and carefully evaluating statistical and systematic errors, and their combination within the different experiments. For the first time, two different formalisms were studied concerning the way the distance is taken into account in the time-lag computation. Others could be added in the future (see, e.g., Amelino-Camelia et al. 2021). As a next step, the software developed for this work will be applied to all available data sets recorded so far by H.E.S.S., MAGIC, and VERITAS, and perhaps by other experiments, and the results will be published in the second part of this work.

Another important advantage of likelihood analysis is its adaptability. Indeed, nothing prevents, in principle, including other effects on the production or propagation of photons in the PDF. Two examples can be pointed out. First, as mentioned in the introduction, it is known that LIV could modify the absorption of VHE photons by the EBL changing the shape of high-energy spectra. Assuming that QG affects both the photon group velocity and photon interactions, the likelihood technique could be used to provide combined EBL and delay constraints on QG models. It is not clear, however, whether the different effects would manifest at the same energy scale or if a different energy scale is applicable for each effect. Second, it should be possible to include other types of delays in the probability function to probe both propagation and source-intrinsic time lags. Despite some recent exploratory work on that topic (see, e.g., Perennes et al. 2020, for the case of blazar flares), the latter are still poorly understood. In addition, intrinsic effects are most probably different from one type of source to another, and even from one subtype to another: short or long GRBs, blazars or flat-spectrum radio quasars. We therefore chose not to include them in the present study. Intrinsic effects are a critical aspect, and they will need to be addressed in the future.

The Cherenkov Telescope Array ¹⁶ (CTA) will start operating soon, superseding the current generation of IACTs in the years 2025–2030 (Acharya et al. 2019). Thanks to its better overall performance and dedicated observation strategies to maximize the number of transient-event detections, it is expected that both CTA arrays (one in each hemisphere) will be able to detect a large number of PSRs, AGN flares, and GRBs. Different subarray configurations will be used in order to optimize the observation program and combining data will therefore become very important. As a result, CTA will be much more sensitive to LIV effects than current-generation experiments. The tools developed in this work will be made publicly available concurrently with the publication of the second paper and adapted to be used in CTA analysis software architecture.

The authors would like to thank collaborations H.E.S.S., MAGIC, and VERITAS for their support in the making of this joint effort and for allowing the use of IRFs for the set of sources used in this paper. They also would like to acknowledge networking support by the COST Action CA18108 (https://qg-mm.unizar.es/). This project has received funding from the European Union's Horizon 2020 research and innovation program under the Marie Skłodowska-Curie grant agreement No. 754510, from the ERDF under the Spanish Ministerio de Ciencia e Innovación (MICINN), grants PID2019-107847RB-C41 and PID2019-107847RB-C42, the Centro de Excelencia "Severo Ochoa" (SEV-2016-0588), and from the CERCA program of the Generalitat de Catalunya. T.T. acknowledges funding from the University of Rijeka, project number uniri-prirod-18-48, and from the Croatian Science Foundation (HrZZ), project number IP-2016-06-9782.

The authors would like to thank G. D'Amico for his useful comments on the draft as well as G. Rosati and C. Pfeifer for insightful discussions on the lag redshift dependence. Finally, the authors express their gratitude to the anonymous referee who helped clarify some parts of the paper.

This paper is dedicated to the memory of our colleague and friend A. Jacholkowska, who initiated this work and put it on the best track toward successful completion.

Facilities: HESS - High Energy Stereoscopic System, MAGIC - Major Atmospheric Gamma Imaging Cherenkov Telescope, VERITAS - Very Energetic Radiation Imaging Telescope Array System.

Software: ROOT.

All systematic uncertainties are listed for individual sources and combinations in Table A1 (J&P case) and Table A2 (DSR case).

Figures 8 and 9 show the calibration plots λ_rec versus λ_inj for all AGNs combined and all PSRs combined, respectively. For PSRs, note that the scale is not the same for n = 1 and n = 2. This leads to an apparent higher value of uncertainty.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Initially created by A.J. and M.M., the task force was led by M.M. for the MAGIC Collaboration, A.N.O. for the VERITAS Collaboration, and J.B. for the H.E.S.S. Collaboration. J.B. acted as the main task force leader after the passing of A.J. in 2018. He was also the principal coordinator for the writing of the present article.

Software development activities were shared between S.C., A.G., D.K., C.L., T.L., L.N., C.P., and T.T. M.R. contributed in studying the J&P and DSR redshift dependence of the time delays, as well as in the writing of the introduction. J.B., S.C., M.G., A.G., D.K., C.L., T.L., T.T. provided the IRFs used in the paper. S.C. and C.L. were responsible for producing the final results and plots. Finally, all the authors had a significant contribution in writing and editing the draft.