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

Julien Bolmont Author contributions are detailed in Appendix C.  Sorbonne
Université, CNRS/IN2P3,
Laboratoire de Physique Nucléaire et de Hautes Energies, LPNHE,
4 Place Jussieu, F-75005 Paris, France  Sami Caroff Present address:
Laboratoire d’Annecy de Physique des Particules, Université Grenoble Alpes,
Université Savoie Mont Blanc, CNRS, LAPP, F-74000 Annecy, France  Sorbonne
Université, CNRS/IN2P3,
Laboratoire de Physique Nucléaire et de Hautes Energies, LPNHE,
4 Place Jussieu, F-75005 Paris, France  Markus Gaug Universitat Autònoma de
Barcelona and CERES-IEEC, E-08193 Bellaterra, Spain  Alasdair Gent School of
Physics & Center for Relativistic Astrophysics,
Georgia Institute of Technology,
837 State Street NW, Atlanta, GA 30332-0430  Agnieszka Jacholkowska  Deceased.
Sorbonne Université, CNRS/IN2P3,
Laboratoire de Physique Nucléaire et de Hautes Energies, LPNHE,
4 Place Jussieu, F-75005 Paris, France  Daniel Kerszberg Institut de Física
d’Altes Energies (IFAE),
The Barcelona Institute of Science and Technology (BIST), E-08193 Bellaterra
(Barcelona), Spain  Christelle Levy  Sorbonne Université, CNRS/IN2P3,
Laboratoire de Physique Nucléaire et de Hautes Energies, LPNHE,
4 Place Jussieu, F-75005 Paris, France  LUTH, Observatoire de Paris, PSL
Research University,
CNRS, Université de Paris,
5 Place Jules Janssen, F-92190 Meudon, France  Tony Lin  Physics Department,
McGill University,
Montreal, QC H3A 2T8, Canada  Manel Martinez Institut de Física d’Altes
Energies (IFAE),
The Barcelona Institute of Science and Technology (BIST), E-08193 Bellaterra
(Barcelona), Spain  Leyre Nogués  Institut de Física d’Altes Energies (IFAE),
The Barcelona Institute of Science and Technology (BIST), E-08193 Bellaterra
(Barcelona), Spain  A. Nepomuk Otte School of Physics & Center for
Relativistic Astrophysics,
Georgia Institute of Technology,
837 State Street NW, Atlanta, GA 30332-0430  Cédric Perennes  Università di
Padova and INFN,
I-35131 Padova, Italy  Michele Ronco Sorbonne Université, CNRS/IN2P3,
Laboratoire de Physique Nucléaire et de Hautes Energies, LPNHE,
4 Place Jussieu, F-75005 Paris, France  Tomislav Terzić University of Rijeka,
Department of Physics, 51000 Rijeka, Croatia

###### Abstract

Gamma-ray astronomy has become one of the main experimental ways to test the
modified dispersion relations (MDRs) of photons in vacuum, obtained in some
attempts to formulate a theory of Quantum Gravity. The MDRs in use imply time
delays which depend on the energy, and which increase with distance following
some function of redshift. The use of transient, or variable, distant and
highly energetic sources, already allows us to set stringent limits on the
energy scale related to this phenomenon, usually thought to be of the order of
the Planck energy, but robust conclusions on the existence of MDR-related
propagation effects still require the analysis of a large population of
sources.

In order to gather the biggest sample of sources possible for MDR searches at
teraelectronvolt energies, the H.E.S.S., MAGIC and VERITAS collaborations
enacted a joint task force to combine all their relevant data to constrain the
Quantum Gravity energy scale. In the present article, the likelihood method
used, to combine the data and provide a common limit, is described in detail
and tested through simulations of recorded data sets for a gamma-ray burst,
three flaring active galactic nuclei and two pulsars. Statistical and
systematic errors are assessed and included in the likelihood as nuisance
parameters. In addition, a comparison of two different formalisms for distance
dependence of the time lags is performed for the first time. In a second
article, to appear later, the method will be applied on all relevant data from
the three experiments.

gamma rays: general, methods: statistical, active galactic nuclei, gamma-ray
bursts, pulsars, quantum gravity, Lorentz invariance violation

††journal: ApJ††facilities: HESS,MAGIC,VERITAS††software: ROOT

##  1 Introduction

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 eighty 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 et al., 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 recent years have held important advances in the field of QG
phenomenology (Mattingly, 2005; Amelino-Camelia, 2013; Liberati, 2013).

It is notoriously difficult to extract observable predictions from fully-
fledged QG approaches. Different models usually start from different
conceptual premises and use different mathematical formalisms in such a way it
is difficult to determine whether they make compatible predictions. In some
cases, the formal complexity forbids producing observable outcomes at all.
Then, to guide experimental efforts, bottom-up approaches have been proposed
(Amelino-Camelia, 2002; Kowalski-Glikman & 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 90s, independent semi-classical analyses
inspired by QG models brought to the attention of the QG community the fact
that it is a highly non-trivial task to retain Lorentz symmetries when
quantizing the space-time geometry of GR. These models include, most notably,
String Theory (see e.g. Mavromatos, 2010, and references therein), Loop
Quantum Gravity (Gambini & Pullin, 1999), Non-commutative 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 in modifying the energy-momentum dispersion relation of free
relativistic particles as follows (Amelino-Camelia et al., 1998):

$E^{2}\simeq
p^{2}c^{2}\times\left[1\pm\sum_{n=1}^{\infty}\left(\frac{E}{E_{QG}}\right)^{n}\right],$
(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{\hbar c^{5}/G}\,\simeq 10^{19}$ GeV). The sign $\pm$ in
Equation (1) takes into account the possibility to have 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 space-time 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 of Equation (1) in the 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 space-time ‘fuzziness’ (Vasileiou et al., 2015). Possible modifications of
the cross section of $\gamma\gamma$ interaction between high-energy photons
and the extra-galactic 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 ($\sim 10^{-33}$ cm), QG effects are believed to cause fluctuations
of space-time geometry which, then, would behave as a dynamical medium
characterised by a non-trivial refractive index. Consequently, photons with
different energies would have different interactions with the ‘foamy’
structure of space-time (sometimes called ‘quantum space-time’) 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 non-integer 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 away from the
Planck scale (being either above or below $E_{P}$). Among others, let us
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 et al., 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^{-19}-10^{-14}$), it has been recognized that very distant
astrophysical sources could be used to probe properties of quantum space-time
(Amelino-Camelia et al., 1998; Alfaro et al., 1999; Jacobson et al., 2006).
Indeed, if the emitted photons travel over large distances, then even
extremely tiny quantum-space-time 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 GRB and flaring active galactic nuclei (AGN) are good
candidates 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-space-time 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.

Table 1: A selection of limits for a sub-luminal propagation obtained with
various instruments and various types of objects.

Source | Experiment | Year |  Distancea |  Lower limit on $E_{QG,1}$ |  Lower limit on $E_{QG,2}$ | Reference | Note
---|---|---|---|---|---|---|---
|  |  |  | (95% CL, GeV) | (95% CL, GeV) |  |
35 GRB | BATSE, HETE-2, Swift | - | - | $1.4\times 10^{16}$ | - | 1 | b
8 GRB | Fermi LAT | - | - | $1.0\times 10^{17}$ | - | 2 |
GRB 090510 (catalog ) | Fermi LAT | 2009 | 0.903 | $9.3\times 10^{19}$ | $1.3\times 10^{11}$ | 3 |
GRB 190114C (catalog ) | MAGIC | 2019 | 0.4245 | $0.6\times 10^{19}$ | $6.3\times 10^{10}$ | 4 |
Mrk 501 (catalog ) | MAGIC | 2005 | 0.034 | $0.3\times 10^{18}$ | $5.7\times 10^{10}$ | 5 | c,⋆
Mrk 501 (catalog ) | H.E.S.S. | 2014 | 0.034 | $3.6\times 10^{17}$ | $8.5\times 10^{10}$ | 6 |
PKS 2155-304 (catalog ) | H.E.S.S. | 2006 | 0.116 | $2.1\times 10^{18}$ | $6.4\times 10^{10}$ | 7 | ⋆
PG 1553+113 (catalog ) | H.E.S.S. | 2012 | $0.49\pm 0.04$ | $4.1\times 10^{17}$ | $2.1\times 10^{10}$ | 8 | d,⋆
PSR B0531+21 (catalog ) | VERITAS | 2007-14 | 2.2 kpc | $1.9\times 10^{17}$ | - | 9 |
PSR B0531+21 (catalog ) | MAGIC | - | 2.2 kpc | $5.5\times 10^{17}$ | $5.9\times 10^{10}$ | 10 | ⋆
PSR B0833-45 (catalog ) | H.E.S.S. | - | 294 pc | $4.0\times 10^{15}$ | - | 11 | ⋆

Notes.
a Redshift is given for extra galactic 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 integral in Eq. 4. Only the limit
obtained for a linear correction is given. c These numbers are actually
reported as best fit values by Martinez & Errando (2009). d The redshift of
this source was not measured but only estimated. ⋆ 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) Martinez & 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).

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:

$\Delta
t_{n}\simeq\pm\,\frac{n+1}{2}\,\frac{E_{h}^{n}-E_{l}^{n}}{\mathrm{H}_{\mathrm{0}}E_{QG}^{n}}\
\kappa_{n}(z),$ (2)

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

The delay $\Delta t_{n}$ takes only into account Lorentz violation effects,
therefore neglecting any time lag originating from emission mechanisms, also
referred to as ‘source intrinsic’ delay. A hint of such kind of delay are
observed for GRB (Ajello et al., 2019) and has been also reported once in the
case of an AGN, for the flare of Mrk 501 (catalog ) in 2005 recorded by MAGIC
111Major Atmospheric Gamma Imaging Cherenkov, https://magic.mpp.mpg.de (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 on-going
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., 2021). 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 $\lambda_{n}$ can be defined as

$\lambda_{n}\equiv\frac{\Delta t_{n}}{\Delta E_{n}\
\kappa_{n}(z)}=\pm\frac{n+1}{2\mathrm{H}_{\mathrm{0}}\ E^{n}_{QG}},$ (3)

using the simplified notation $\Delta E_{n}\equiv E_{h}^{n}-E_{l}^{n}$. This
parameter $\lambda_{n}$, which will be used later, has the advantage to be
independent of the distance of the source and is therefore suitable for a
multi-source analysis.

Since the late 90s, the field has rapidly expanded, with more and more sources
being analyzed in the search for Lorentz invariance violation effects. With
the notable exception of the flare of Mrk 501 (catalog ) in 2005 already
mentioned above, no significant delay has been reported so far when using only
photon 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 (AGN, GRB and PSR). This new notation $E_{QG,n}$ reflects the
fact that LIV analyses have different sensitivities for linear and quadratic
effects.

The results of Ellis et al. (2006, 2008, 2019) are of particular interest
since they were obtained from the analyses of several GRB. This kind of
analysis, repeated with different experiments (e.g. Bolmont et al., 2008;
Bernardini et al., 2017) consists in two steps: first, the time-lags are
computed for each individual source, and then the obtained data points are
fitted with a function $\Delta 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 GRB that Lorentz invariance could be violated at a scale
of $\sim 3.7\times 10^{17}$ GeV (Xu & Ma, 2016a, b, 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. 222High Energy
Stereoscopic System, https://www.mpi-hd.mpg.de/hfm/HESS/, MAGIC and VERITAS
333Very Energetic Radiation Imaging Telescope Array System,
https://veritas.sao.arizona.edu 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 $\lambda_{n}$ and
$\tau_{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}$.

##  2 Redshift dependence of time delays

It is rather natural to believe that curvature and quantum effects are deeply
intertwined since curvature is a key characteristic of space-time 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 a particular interest and has been
systematically used so far in experimental analyses constraining in-vacuo
dispersion.

Figure 1: Parameter $\kappa$ for $n=1$ (black) and $n=2$ (gray) in the J&P
case (solid line) and in the DSR case (dashed line).

In this approach, parameter $\kappa_{n}(z)$ is expressed as:

$\kappa^{\mathrm{J\&P}}_{n}(z)\equiv\int_{0}^{z}\frac{(1+z^{\prime})^{n}}{\sqrt{\Omega_{m}\,(1+z^{\prime})^{3}+\Omega_{\Lambda}}}\
dz^{\prime},$ (4)

where the denominator relates to the Hubble parameter
$H(z)=\mathrm{H}_{\mathrm{0}}\sqrt{\Omega_{m}\,(1+z)^{3}+\Omega_{\Lambda}}$.
In the following, cosmological parameters values are taken from Planck results
(Aghanim et al., 2018) as recommended by the Particle Data Group (Zyla et al.,
2020): $\mathrm{H}_{\mathrm{0}}=67.4\pm 0.5\
\mathrm{km\,s}^{-1}\,\mathrm{Mpc}^{-1}$, $\Omega_{m}=0.315\pm 0.007$ and
$\Omega_{\Lambda}=0.685\pm 0.007$.

As recent literature has pointed out (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 Deformed Special
Relativity (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):

$\kappa^{\mathrm{DSR}}_{n}(z)\equiv\int_{0}^{z}\frac{h^{2n}(z^{\prime})dz^{\prime}}{(1+z^{\prime})^{n}\,\sqrt{\Omega_{m}\,(1+z^{\prime})^{3}+\Omega_{\Lambda}}},$
(5)

with

$\begin{split}h(z^{\prime})\equiv
1+z^{\prime}-&\sqrt{\Omega_{m}\,(1+z^{\prime})^{3}+\Omega_{\Lambda}}\\\
&\times\int_{0}^{z^{\prime}}\frac{dz^{\prime\prime}}{\sqrt{\Omega_{m}\,(1+z^{\prime\prime})^{3}+\Omega_{\Lambda}}}\,\end{split}$
(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 functions $\kappa^{\mathrm{DSR}}$ and $\kappa^{\mathrm{J\&P}}$
as a function of redshift for $n=1$ and $n=2$. $\kappa^{\mathrm{DSR}}$ is
smaller than $\kappa^{\mathrm{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 $\kappa^{\mathrm{J\&P}}$ and $\kappa^{\mathrm{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\sim
1$ in the TeV range.

To conclude this section, let us add that in case of nearby sources, such as
pulsars, the euclidean approximation is valid, i.e.
$\kappa_{n}(z)=d\,\mathrm{H}_{\mathrm{0}}/c$ where $d$ is the euclidean
distance to the source. In addition, the ratio
$\kappa^{\mathrm{DSR}}/\kappa^{\mathrm{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.

##  3 Methodology

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), 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 444https://root.cern.ch framework (Brun & Rademakers, 1997).

###  3.1 Single Source Likelihood

First applied by Martinez & Errando (2009) for analyzing the 2005 flare of Mrk
501 (catalog ) 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 $\lambda_{n}$ as the single parameter to be estimated.
The PDF is obtained convolving the spectrum of the source $\Gamma_{s}(E_{t})$
with the light curve $C_{s}$, both as observed on Earth, i.e. after
propagation:

$F_{s}(E_{t},t;\lambda_{n})=\frac{\Gamma_{s}(E_{t})\,C_{s}\left(t-D(E_{t},\lambda_{n},z)\right)}{N_{s}},$
(7)

where $E_{t}$ is the true energy of the gamma-ray photon, $t$ is the arrival
time corrected by the factor:

$D(E_{t},\lambda_{n},z)=\lambda_{n}\times\kappa_{n}(z)\times E_{t}^{n},$ (8)

which defines the propagation delay due to LIV, where $\lambda_{n}$ is given
by Equation (3), and $N_{s}$ is a normalization term expressed as follows:

$N_{s}=\iint\Gamma_{s}(E_{t})\,C_{s}\left(t-D(E_{t},\lambda_{n},z)\right)\,dE_{t}\,dt.$
(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. Since there is no fully
accepted model available which reproduces the shape of the light curves for
GRB, AGN or PSR, Gaussian or Lorentzian functions, or the sum of several of
these functions are usually used. The function $\Gamma_{s}$ is obtained from
the data on the full energy range considered for the LIV analysis (see Section
4).

Table 2: Nuisance parameter uncertainties for the individual sources.

Source | Energy Scale | Background proportion |  Spectral indexa | Distance/redshift |  Referencesb
---|---|---|---|---|---
GRB 190114C |  $17$% |  $11$% | $0.21$ | $\Delta z=1\times 10^{-3}$ | 1, 2, 3, 4
PG 1553+113 |  $10$% |  $20$% | $0.31$ | $\Delta z=4\times 10^{-2}$ | 5, 6
Mrk 501 |  $17$% |  $11$% | $0.04$ | $\Delta z=1\times 10^{-4}$ | 1, 2, 3, 6
PKS 2155-304 |  $10$% |  $20$% | $0.1$ | $\Delta z=1.7\times 10^{-2}$ | 5, 7
Crab (M) |  $17$% |  $11$% | $0.07$ |  $\Delta d=506$ pc | 1, 2, 3, 8
Crab (V) |  $20$% |  $22$% | $0.5$ |  $\Delta d=506$ pc | 9, 8
Vela |  $10$% |  $20$% | $0.67$ |  $\Delta d=76$ pc | 5, 10

Note.
a Uncertainty for 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 energy
scale and spectral index, Pueschel (2019). The value for background proportion
was provided by the VERITAS collaboration, (10) Caraveo et al. (2001)

Background events of several different origins are also taken into account.
They include hadrons mis-reconstructed 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:

$F_{b,k}(E_{t},t)=\frac{\Gamma_{b,k}(E_{t})\,C_{b,k}}{N_{b,k}}.$ (10)

$\Gamma_{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:

$N_{b,k}=\iint\Gamma_{b,k}(E_{t})\,C_{b,k}\,dE_{t}dt.$ (11)

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

$\frac{dP}{dE_{m}dt}=w_{s}\ \frac{\int
A(E_{t},\vec{\varepsilon})M(E_{t},E_{m})\times
F_{s}(E_{t},t;\lambda_{n})dE_{t}}{N^{\prime}_{s}}\\\ +\sum_{k}w_{b,k}\
\frac{\int A(E_{t},\vec{\varepsilon})M(E_{t},E_{m})\times
F_{b,k}(E_{t},t)dE_{t}}{N^{\prime}_{b,k}},$ (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}+\sum_{k}w_{b,k}=1$. $E_{t}$
still denotes the true energy while $E_{m}$ is the corresponding measured
energy. Parameters $N^{\prime}_{s}$ and $N^{\prime}_{b,k}$ are the
normalization factors of the PDF. In addition to energy $E_{t}$, 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 for either a measurement or the derivation of lower
limits on $\lambda_{n}$ (and $E_{QG,n}$) can then be obtained summing the log-
likelihood of all the events for a given source $S$:

$L_{S}(\lambda_{n})=-\sum_{\mathrm{i}}\log\left(\frac{dP}{dE_{m}dt}(E_{m,i},t_{i});\lambda_{n}\right).$
(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}(\lambda_{n})$ are obtained for all sources, the combined log-likelihood
$L_{comb}$ is simply given by their sum:

$L_{comb}(\lambda_{n})=\sum_{\mathrm{all\ sources}}L_{S}(\lambda_{n}).$ (14)

###  3.3 Statistical and systematic uncertainties

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

$L(\lambda_{n},\vec{\theta})=L_{\mathrm{S}}(\lambda_{n},\vec{\theta})+L_{\mathrm{template}}(\vec{\theta}_{\mathrm{C}})+L_{\mathrm{\gamma}}(\theta_{\mathrm{\gamma}})+\\\
L_{\mathrm{B}}(\vec{\theta}_{\mathrm{B}})+L_{\mathrm{ES}}(\theta_{\mathrm{ES}})+L_{\mathrm{z}}(\theta_{\mathrm{z}}),$
(15)

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

  1. a.

$\vec{\theta}_{\mathrm{C}}$, the parameters of the light curve analytic
parameterization,

  2. b.

$\theta_{\mathrm{\gamma}}$, the power law index of signal events spectrum,

  3. c.

$\vec{\theta}_{\mathrm{B}}$, the ratio of signal and of background event
numbers to the total number of events,

  4. d.

$\theta_{\mathrm{ES}}$, the energy scale,

  5. e.

$\theta_{\mathrm{z}}$, the distance or redshift.

Table 3: Simulation settings for the individual sources.

Source | Energy Range |  Time Rangea | Spectral index | Lightcurve shape | Number of events | Background proportion
---|---|---|---|---|---|---
| (TeV) |  |  $\Gamma_{s}$, $\Gamma_{b}$ |  |  likelihoodb, templatec | hadronic, baseline
GRB 190114C |  $0.3$ \- $2$ |  $60$ \- $1200$ s |  $5.43$, - | curved power law |  $726$, - |  $0.055$, $0.$
PG 1553+113 |  $0.4$ \- $0.8$ |  $0$ \- $8000$ s |  $4.8$, $4.8$ | double Gauss |  $72$, $82$ |  $0.29$, $0.15$
Mrk 501 |  $0.25$ \- $11$ |  0 - $1531$ s |  $2.2$, $2.2$ | single Gauss |  $1800$, - |  $0.39$, $0.$
PKS 2155-304 |  $0.28$ \- $4$ |  $0$ \- $4000$ s |  $3.46$, $3.32$ | 5 asymmetric Gauss |  $2965$, $561$ |  $0.$, $0.02$
Crab (M) |  $0.4$ \- $7$ |  $0.36$ \- $0.45$ |  $2.81$, $2.47$ | single Gauss + Baseline |  $14869$, - |  $0.$, $0.961$
Crab (V) |  $0.2$ \- $10$ |  $0.37$ \- $0.43$ |  $3.25$, $2.47$ | single Gauss + Baseline |  $22764$, - |  $0.$, $0.964$
Vela |  $0.06$ \- $0.15$ |  $0.50$ \- $0.60$ |  $3.9$, $1.75$ | asymmetric Lorentzian |  $330820$, - |  $0.$, $0.998$

Notes.
a For pulsars, the phase range is given, i.e. the time range normalized with
respect to the rotation period. b Number of photons considered when computing
the likelihood, i.e. excluding the ones used for template determination. c A
sign ’-’ means no template was used (see Section 3.3 for details).

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}_{\mathrm{C}})$ as the
sum of the log-likelihoods of each event generated from the low energy
template parameterization:

$L_{\mathrm{template}}(\vec{\theta}_{\mathrm{C}})=-\sum^{N_{\mathrm{template}}}_{i=1}\log\left(\frac{C_{s}(t_{i},\vec{\theta}_{\mathrm{C}})}{N_{c}}\right),$
(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}_{\mathrm{C}})$
denotes the fact the template is evaluated for a zero-lag
($D(E_{i},\lambda_{n}=0,z)=0$), and explicitly shows the parameter vector
$\vec{\theta}_{\mathrm{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}_{\mathrm{C}})=0$ and the uncertainty
on $\vec{\theta}_{\mathrm{C}}$ is then accounted for in the generated data
sample log-likelihood $L_{\mathrm{S}}(\lambda,\vec{\theta})$ defined in
Equation (13).

$L_{\mathrm{\gamma}}(\theta_{\mathrm{\gamma}})$ 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_{\mathrm{B}}(\vec{\theta}_{\mathrm{B}})$ and
$L_{\mathrm{ES}}(\theta_{\mathrm{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_{\mathrm{z}}(\theta_{\mathrm{z}})$.

For the power law index, 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:

$L_{\mathrm{x}}(\vec{\theta}_{\mathrm{x}})=\sum_{i}\frac{(\theta_{\mathrm{x},i}-\bar{\theta}_{\mathrm{x},i})^{2}}{2\sigma^{2}_{\theta_{\mathrm{x},i}}},$
(17)

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

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 statistical and systematical uncertainty.
They are presented in the Appendix with Table 5 for the J&P case and Table 6
for the DSR case, for each source and each source combination. These results
will be commented further in Section 5.

##  4 Simulations

###  4.1 Simulated data sets

####  4.1.1 Data sets choice criteria

The sources used in this study are listed in Table 3. They have all been
detected by the three experiments H.E.S.S., MAGIC and VERITAS, and have been
selected to gather a representative sample. Namely, the three types of source
were selected: one GRB, three flaring AGN and two PSR, with LIV results
already published, and with the following additional criteria:

  1. a.

The three flaring AGN show different signal to background ratios: negligible
background for PKS 2155-304 (catalog ) and Mrk 501 (catalog ) and substantial
background level for PG 1553+113 (catalog ),

  2. b.

The sources show very different light curve shapes, from a single Gaussian
pulse for Mrk 501 (catalog ) flare of 2005 to multiple asymmetric spikes for
PKS 2155-304 (catalog ) flare of 2006,

  3. c.

The sources selected cover a wide range in distance, from 2 kpc for the Crab
PSR to a redshift of 0.49 for PG 1553+113 (catalog ),

  4. d.

The two PSR have different distances and were observed on very different time
scales,

  5. e.

In addition, PG 1553+113 (catalog ) has a large uncertainty on the distance
which was taken into account in the analysis.

In the following sub-section, the most important characteristics of the
sources as taken from the references listed in Table 1 are briefly summarized.
The numbers given in Table 3 were extracted from these references or provided
by the authors in private communications. Then, the use of simulated data-sets
to assess the performance of the method is described in Section 4.2.

Except specified otherwise, a spectral index $\Gamma_{k}=2.7$ was used for
hadrons.

####  4.1.2 Source description

GRB 190114C (catalog ) is a gamma-ray burst detected on 2019 January 14 at
20:57:03 Universal Time (UT) and located at redshift $z=0.4245\pm 0.005$
(Selsing et al., 2019; Castro-Tirado et al., 2019). Following the alert sent
by Swift (Gropp et al., 2019), MAGIC observed GRB 190114C (catalog ),
detecting a strong VHE $\gamma$-ray signal (Acciari et al., 2019a). The
observations started 62 seconds after the beginning of the burst. A total of
$\sim$700 events with energy ranging from 300 GeV to $\sim$2 TeV were recorded
during the first 19 minutes of observations. The intrinsic energy distribution
of the signal was fitted with a power law of index $2.5\pm 0.2$, leading to an
index of $\Gamma_{s}=5.43\pm 0.22$ (statistical error only) when EBL
absorption is taken into account. The time distribution of the events recorded
by MAGIC follows a power law with index $1.51\pm 0.04$. MAGIC did not observe
the peak of the burst. Therefore, the light curve of the full burst, including
the sharp rise to the peak flux followed by a power-law decay was modeled
based on multiwavelength observations of the event and theoretical inference
(Acciari et al., 2019b). The prompt emission of GRB 190114C inferred from the
keV–MeV light curves and spectra lasted no more than $25$ seconds after the
onset of the GRB. This indicates that the emission observed by MAGIC is
associated with the afterglow phase, rather than with the prompt phase, which
typically shows irregular variability. However, as reported in Acciari et al.
(2020), a sub-dominant contribution from the prompt phase (at most 20%) at
early times of the afterglow ($t\lesssim 100$ s) cannot be entirely excluded.
The lower bound of 60 s chosen for the present study was chosen to minimize
this contribution while retaining statistics as high as possible.

Mrk 501 (catalog ) is a BL Lac object at redshift $z=0.03364$. The flare of
2005 July 9 was detected by the MAGIC telescope, at the time operating in
monoscopic configuration (Albert et al., 2007). The flux of this flare reached
a peak more than a factor of two higher than before and after the flare. A
total of $\sim$1800 events with energy from 0.15 to 10 TeV were recorded
during the flare among which $\sim$700 could be associated to the background.
The energy distribution of the signal and baseline events is well described by
a power law of index $\Gamma_{s,b}=2.2$ while the time distribution was
parameterized by a single Gaussian spanning over $1600$ seconds.

PKS 2155-304 (catalog ) is another BL Lac object at higher redshift $z=0.116$.
The flare of 2006 July 28 detected by H.E.S.S. telescopes is seemingly one of
the brightest flares detected by the experiment so far with a signal to noise
ratio exceeding 300 (Aharonian et al., 2007). The lightcurve is parameterized
by five asymmetric Gaussians with $2\%$ background over $4000$ seconds with a
total of $3526$ photons. The energy distribution is described by a power law
of index $\Gamma_{s}=3.46$ ranging from 0.25 to 4 TeV during the flare while
the quiescent state leads to an index of $\Gamma_{b}=3.32$.

PG 1553+113 (catalog ), yet another BL Lac object, is the furthest source of
this list with an estimated redshift $z=0.49\pm 0.04$. The flare of 2012 April
26-27 was detected by H.E.S.S. telescopes where its flux increased three-fold
as compared to its quiescent state (Abramowski et al., 2015). The time
distribution was parameterized by two Gaussians with 154 photons over $8000$
seconds, where background accounts for $44\%$ of the events with $30\%$ gamma-
like hadrons and $14\%$ baseline photons. The energy distribution spreading
between $0.3$ and $0.8$ TeV is described by a power law of index
$\Gamma_{s,b}=4.8$ for signal and baseline photons.

The Vela Pulsar (PSR B0833-45 (catalog )) located at $294\pm 76$ pc rotates
with a periodicity of $89$ ms. The data simulated in this work is from a
compilation of observations with H.E.S.S. large telescope from March 2013 to
April 2014, for which a LIV analysis was performed (Chretien et al., 2015;
Chretien, 2015). 330,820 pulsed events between 60 and 150 GeV were recorded
with a signal to noise ratio of $0.012$. The phase distribution is
parameterized by an asymmetric Lorentzian between $0.5$ and $0.6$. Background
accounts for $98.8\%$ of the events with only baseline photons. The energy
distribution is described by a power law of index $\Gamma_{s}=3.9$ for signal
and $\Gamma_{b}=1.75$ for baseline photons.

Figure 2: Bias $\lambda_{rec}-\lambda_{inj}$ vs. number of bins for GRB
190114C in the linear case and J&P formalism. The number of bins in the table
is chosen so that the bias (black line) is compatible with zero within its
1$\sigma$ uncertainty range (gray envelope). The same number of bins is used
for measured energy, non-delayed arrival times, and time delays.

The Crab Pulsar (PSR B0531+21 (catalog )) has a $33.7$ ms period and is
located at $2.0\pm 0.5$ kpc. One of the data sets used in this work, noted
“Crab M” hereafter, is a compilation of observations with MAGIC telescopes
from 2005 to 2017 (Ahnen et al., 2017). $3080\pm 460$ excess events from the
P2 region of the phase were recorded, from which $544\pm 92$ have a
reconstructed energy above 400 GeV and are used in the LIV analysis. The phase
distribution of the P2 peak was parameterized by a Gaussian. A profiling of
the nuisance parameters yielded a mean of $\Phi=0.403$ (respectively $0.401$)
and standard deviation $0.015$ ($0.011$) for $n=1$ ($n=2$). Background
accounts for 96% of the events with only baseline photons. The energy
distribution was described by a power law of index $\Gamma_{s}=2.81$ for
signal and $\Gamma_{b,k}=2.47$ for combined background events and baseline
photons.

The other data set, noted “Crab V”, is a compilation of high quality data
taken with VERITAS telescopes between 2007 and 2011. 22,764 pulsed events were
recorded from the P2 region and its baseline where background account for
$96.4\%$ of the events with again only baseline photons (Zitzer et al., 2013).
The phase distribution was also parameterized with a Gaussian centered on
$0.398$ with a standard deviation of $0.0116$. The energy distribution was
again described by a power law of index $\Gamma_{s}=3.25$ for signal and
$\Gamma_{b}=2.47$ for baseline photons.

###  4.2 Method calibration and performance

The normalization factor $N^{\prime}_{s}$ of the PDF of Equation (12) is a
triple integral, computation of which is particularly time consuming since it
needs to be done for each minimization step and for each event of the sample.
To decrease the computation time, the PDF is pre-calculated and stored in
tables binned over measured energy $E_{m}$, non-delayed arrival times $t$, and
time-delays $D(E_{t},\lambda_{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.

Figure 3: The center plot shows the distribution of the reconstructed lag in
the case of GRB 190114C, J&P formalism for the linear case. The plot on the
left (respectively on the right) shows the distribution of the lower (upper)
limits of the confidence interval for 68% CL. The three distributions are
obtained with a zero injected lag. The histograms are fitted with asymmetric
Gaussian functions, parameters of which are used in turn to produce
calibration plots (see the text for details).

Figure 4: Calibration plot showing $\lambda_{rec}$ vs. $\lambda_{inj}$ for GRB
199114C (left) and all sources combined (right) in the linear case and J&P
formalism. The light gray area corresponds to the standard deviation of the
$\lambda_{rec}$ distribution while the dark gray region shows the statistical
uncertainty. For both plots, a function $a\ \lambda_{inj}+b$ is fitted (black
line).

The number of bins used in the tables has been chosen for each source to
minimize the bias $\lambda_{rec}-\lambda_{inj}$ between the injected
($\lambda_{inj}$) and reconstructed ($\lambda_{rec}$) time delays. An example
is shown in Figure 2 for GRB 190114C (catalog ). In this particular case, the
plot shows that a minimum of $\sim$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, for linear and quadratic LIV
effects.

In order to assess the sensitivity and precision of the lag reconstruction,
simulated data sets were produced with different values for $\lambda_{inj}$.
For each value of the injected lag, one thousand realizations of the light
curve were simulated. The distribution of reconstructed values of
$\lambda_{n}$ is shown in the central panel of Figure 3 for GRB 190114C
(catalog ), in the J&P case and $n=1$. Lower and upper limits of the
confidence intervals are taken as the values of $\lambda_{n}$ for which
$2\,[L_{S}(\lambda_{n})-\mathrm{min}(L_{S})]=1$ for 68% CL and
$2\,[L_{S}(\lambda_{n})-\mathrm{min}(L_{S})]=3.84$ for 95% CL (Equation 13).
Their distributions for 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 ($\lambda_{LL}$,
$\lambda_{rec}$ and $\lambda_{UL}$) and standard deviations separately defined
on the left and on the right of the maxima ($\sigma_{\lambda,l}$,
$\sigma_{\lambda,r}$). While the latter accounts for statistical uncertainties
only, the lower and upper limits $\lambda_{LL}$ and $\lambda_{UL}$ account for
both statistical and systematic uncertainties.

Figure 5: Limits obtained for all individual sources and combinations for the
linear (top) and quadratic (bottom) cases, for the J&P (dots) and DSR
(crosses) redshift dependence. Blue markers correspond to the case where only
statistical errors are taken into account (“stat only”) while orange markers
correspond to the case where both statistical and systematic errors are
included (“stat+syst”).

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

The plots for $\lambda_{rec}$ versus $\lambda_{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 (catalog ) alone (left) and for all sources combined
(right) in the linear and J&P case. The reconstructed lag is fitted with a
linear function $\lambda_{rec}=a\ \lambda_{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
$\lambda_{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 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.

##  5 Results and discussion

###  5.1 Systematic uncertainties

All systematic uncertainties are listed for individual sources and
combinations in Table 5 (J&P case) and Table 6 (DSR case) in the Appendix. For
most of the individual sources, the dominant systematic is the statistical
uncertainty of the light curve template. Since the time lag intensifies as the
correction order $n$ gets larger, the template uncertainties contribute
comparatively less in the quadratic case than the linear one. For other
individual sources, the precision of the energy distribution of the events
prevails. Indeed, the energy scale uncertainty is found to be the most
important source of systematics for the Crab pulsar observed by MAGIC, the
Vela pulsar and Mrk 501 (catalog ), for the quadratic case. This is expected
since the time delay depends on the energy squared. A similar behavior is
observed for GRB 190114C (catalog ) and the Crab pulsar observed by VERITAS,
where the uncertainty on the spectral slope dominates.

Table 4: 95% CL limits obtained for individual objects and combinations. Source | $E_{QG,1}$ | $E_{QG,2}$
---|---|---
| J&P | DSR | J&P | DSR
|  ($10^{18}$ GeV) |  ($10^{18}$ GeV) |  ($10^{10}$ GeV) |  ($10^{10}$ 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\times 10^{-3}$ | $3.5\times 10^{-3}$ | - | - | $5.6\times 10^{-2}$ | $5.5\times 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

For the combinations, dominant systematic uncertainties are the ones of
sources that dominate the sample. The pulsar combinations are dominated by
template statistics, while the combination of AGN shows a predominance of
template statistics for $n=1$, and a domination of energy scale for $n=2$,
confirming the importance of the energy distribution uncertainty for the
quadratic case. The combinations that include GRB 190114C (catalog ) follow a
very similar trend due to the dominance of the GRB over the other sources.
They show a clear ascendancy of the uncertainty of the power law index, which
is the main source of systematic uncertainty for the combination of all the
sources.

Figure 6: Limits obtained from the simulated data sets for all individual
sources for the linear (left) and quadratic (right) cases for the J&P redshift
dependence. Blue dots show the limits obtained taking into account statistical
errors only (“stat only”) while yellow dots show the limits including both
statistical and systematic errors (“stat+syst”). Blue crosses give the limits
published from actual data sets (Table 1).

###  5.2 Limits

Figure 7: Comparison between the limits obtained in the J&P framework (red
dots) and the limits obtained in the DSR formalism (green crosses) in the
linear case (left) and quadratic case (right). The limits shown include both
statistical and systematic errors.

From equation (3), limits on $E_{QG}$ are given by:

$\left[{\frac{2}{n+1}\left(\lambda_{n,\pm}+\sqrt{\delta_{stat}^{2}+\sigma^{2}\delta_{syst}^{2}}\right)\mathrm{H}_{\mathrm{0}}}\right]^{\frac{1}{n}},$
(18)

where the subscript $\pm$ refers to subluminal and superluminal cases,
$\delta_{stat}$ is the statistical error (standard deviation) on the normally
distributed reconstructed value of $\lambda_{n}$, $\delta_{syst}$ is the
overall systematic error obtained from the values listed in Tables 5 and 6
computed for a confidence level of 68%, and $\sigma$ is a real number allowing
for a shift in confidence level using the same systematic errors. Since
systematic errors are computed for 68% CL, and statistical errors for 95% CL,
$\sigma$ is set to two 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 $\gtrsim 2$ as compared to the case 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.

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 paper nuisance parameters and 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.

As expected from previously published results, GRB 190114C is the most
constraining source due to its high redshift, high variability, 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 (catalog ) and Mrk 501 (catalog ). Due to its smaller
number of events and limited energy range, PG 1553+113 (catalog ) limit is
less constraining, even though its redshift is the highest of all the sources
included in this work. While PKS 2155-304 (catalog ) dominates over the other
sources in the linear case due to its higher redshift and event statistics,
Mrk 501 (catalog ) dominates the limit in the quadratic case due to its energy
range extending twice as high as the one of PKS 2155-304 (catalog ).

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 the 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 (catalog )
or PG 1553+113 (catalog ), hence behaving in accordance with the $\kappa_{n}$
parameter evolution shown in Figure 1. Due to the facts that
$\kappa^{\mathrm{J\&P}}>\kappa^{\mathrm{DSR}}$ and $\kappa^{\mathrm{J\&P}}$
increases faster than $\kappa^{\mathrm{DSR}}$ (Section 2), the J&P model
emphasizes the impact of large redshift sources on the limits. Therefore, the
GRB dominates more in the J&P case than for the DSR case, where all sources
contributions are more balanced.

##  6 Conclusions

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 AGN 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 the 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 production
or propagation of photons in the probability density function. 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 sub-type 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 555https://www.cta-observatory.org (CTA) will
start operating soon, superseding the current generation 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
detection, it is expected that both CTA arrays (one in each hemisphere) will
be able to detect a large number of PSR, AGN flares and GRB. Different sub-
array 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 programme under the Marie Skłodowska-
Curie grant agreement No. 754510, from the ERDF under the Spanish Ministerio
de Ciencia e Innovación (MICINN), grant PID2019-107847RB-C41, 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 lag-redshift dependence. Finally, the
authors express their gratitude to the anonymous referee who helped clarifying
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 tracks towards a successful completion.

##  Appendix A Systematic uncertainties

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

Table 5: Summary of systematic uncertainties for all sources and combinations simulated for the J&P case. Source | Correction | Template | Energy | Background | Uncertainty on | Distance/redshift | Reconstruction | All syst.
---|---|---|---|---|---|---|---|---
| order | statistics | scale | normalization | power law index | uncertainty | uncertainty | combined
|  |  (s.TeV-n) |  (s.TeV-n) |  (s.TeV-n) |  (s.TeV-n) |  (s.TeV-n) |  (s.TeV-n) |  (s.TeV-n)
GRB 190114C | $n=1$ | $17.8$ | $6.9$ | $8.0$ | $9.4$ | $<7.7$ | $3.0$ | $25.6$
$n=2$ | $9.4$ | $12.4$ | $1.7$ | $15.4$ | $<9$ | $4.2$ | $24.1$
PKS 2155-304 | $n=1$ | $101$ | $11.7$ | $<20$ | $<22$ | $17.8$ | $<3.3$ | $107$
$n=2$ | $21.8$ | $19.3$ | $0.7$ | $8.1$ | $12.0$ | $<2.2$ | $37.4$
Mrk 501 | $n=1$ | $155$ | $56$ | $<51$ | $49$ | $1.$ | $<8.5$ | $197$
$n=2$ | $11.2$ | $18.3$ | $<10.3$ | $9.3$ | $0.19$ | $<1.6$ | $28.8$
PG1553+113 | $n=1$ | $631$ | $150$ | $324$ | $<361$ | $112$ | $<64$ | $727$
$n=2$ | $916$ | $638$ | $537$ | $<552$ | $338$ | $<112$ | $1282$
Crab V | $n=1$ | $897$ | $137$ | $<73$ | $142$ | $145$ | $<25$ | $1135$
$n=2$ | $1141$ | $410$ | $<264$ | $694$ | $265$ | $<174$ | $1820$
Crab M | $n=1$ | $371$ | $66$ | $7$ | $23$ | $74$ | $<11$ | $416$
$n=2$ | $167$ | $64.5$ | $61$ | $24$ | $48$ | $<72$ | $190$
Vela | $n=1$ | $1.36\times 10^{4}$ | $1.03\times 10^{4}$ | $0.46\times 10^{4}$ | $<1.3\times 10^{4}$ | $1.30\times 10^{3}$ | $<5.87\times 10^{3}$ | $2.28\times 10^{4}$
$n=2$ | $1.0\times 10^{5}$ | $2.05\times 10^{5}$ | $0.48\times 10^{5}$ | $<1.5\times 10^{5}$ | $1.57\times 10^{5}$ | $<0.95\times 10^{5}$ | $3.05\times 10^{5}$
Crab (M+V) | $n=1$ | $357$ | $49$ | $<56$ | $32$ | $61$ | $<32$ | $398$
$n=2$ | $161$ | $59$ | $45$ | $59$ | $38$ | $<83$ | $197$
PSR | $n=1$ | $355$ | $52$ | $<58$ | $38$ | $58$ | $<11$ | $394$
$n=2$ | $90$ | $71$ | $49$ | $24$ | $62$ | $<55$ | $138$
AGN | $n=1$ | $89.5$ | $12$ | $<15$ | $3.7$ | $15.8$ | $<2.9$ | $94.9$
$n=2$ | $10.1$ | $11.1$ | $<6$ | $6.2$ | $3.4$ | $<1.3$ | $19.7$
AGN+PSR | $n=1$ | $85$ | $11$ | $<18$ | $5$ | $15$ | $<2.9$ | $91$
$n=2$ | $9.6$ | $10.9$ | $<8$ | $5.9$ | $4.5$ | $<1.1$ | $17.8$
GRB+AGN | $n=1$ | $17.8$ | $5.8$ | $6.8$ | $8.3$ | $1.4$ | $3.3$ | $24.5$
$n=2$ | $6.8$ | $7.8$ | $<6.6$ | $9.0$ | $1.7$ | $1.4$ | $16.2$
GRB+PSR | $n=1$ | $17.5$ | $6.7$ | $7.9$ | $9.1$ | $1.0$ | $3.2$ | $24.9$
$n=2$ | $8.1$ | $11.3$ | $1.6$ | $12.7$ | $2.8$ | $<1.1$ | $19.4$
All | $n=1$ | $18.0$ | $5.8$ | $6.7$ | $8.2$ | $1.5$ | $4.1$ | $24.8$
$n=2$ | $7.5$ | $7.7$ | $<6.2$ | $8.2$ | $2.4$ | $4.8$ | $16.4$
Table 6: Summary of systematic uncertainties for all sources and combinations simulated for DSR case. Source | Correction | Template | Energy | Background | Uncertainty on | Distance/redshift | Reconstruction | All syst.
---|---|---|---|---|---|---|---|---
| order | statistics | scale | normalization | power law index | uncertainty | uncertainty | combined
|  |  (s.TeV-n) |  (s.TeV-n) |  (s.TeV-n) |  (s.TeV-n) |  (s.TeV-n) |  (s.TeV-n) |  (s.TeV-n)
GRB 190114C | $n=1$ | $26.2$ | $10.2$ | $11.9$ | $13.9$ | $<11.2$ | $5.5$ | $38.0$
$n=2$ | $18.0$ | $25.5$ | $3.8$ | $30.0$ | $6.2$ | $10.6$ | $47.8$
PKS 2155-304 | $n=1$ | $113$ | $12.7$ | $<22.5$ | $<24.2$ | $17.3$ | $<3.6$ | $119$
$n=2$ | $25.8$ | $23.7$ | $3.4$ | $7.0$ | $14.8$ | $<2.9$ | $45.6$
Mrk 501 | $n=1$ | $160$ | $58$ | $<53$ | $51$ | $1.$ | $<8.0$ | $204$
$n=2$ | $12.0$ | $19.6$ | $<11$ | $10.0$ | $0.2$ | $<1.8$ | $30.9$
PG1553+113 | $n=1$ | $968$ | $311$ | $545$ | $<555$ | $<522$ | $<104$ | $1131$
$n=2$ | $2200$ | $1545$ | $1259$ | $<1377$ | $295$ | $<250$ | $2965$
AGN | $n=1$ | $98.4$ | $12.9$ | $<17$ | $4.2$ | $14.8$ | $<3.2$ | $103$
$n=2$ | $11.1$ | $13.0$ | $<6.6$ | $7.3$ | $2.1$ | $<1.5$ | $22.5$
AGN+PSR | $n=1$ | $94$ | $12$ | $<19$ | $4.3$ | $15$ | $<3.0$ | $99$
$n=2$ | $9.1$ | $11.9$ | $<8.2$ | $6.1$ | $3.9$ | $<1.2$ | $19.1$
GRB+AGN | $n=1$ | $26.2$ | $7.7$ | $9.1$ | $11.2$ | $2.4$ | $1.7$ | $34.7$
$n=2$ | $10.1$ | $11.2$ | $<8.5$ | $9.8$ | $1.7$ | $4.3$ | $21.7$
GRB+PSR | $n=1$ | $26.0$ | $9.7$ | $11.3$ | $13.3$ | $1.8$ | $3.9$ | $37.4$
$n=2$ | $8.0$ | $18.0$ | $<15.4$ | $18.5$ | $6.5$ | $<2.5$ | $28.7$
All | $n=1$ | $27.0$ | $7.7$ | $8.7$ | $10.9$ | $2.8$ | $<4.5$ | $35.6$
$n=2$ | $10.1$ | $11.0$ | $<0.96$ | $8.3$ | $3.2$ | $<4.2$ | $19.8$

##  Appendix B Calibration plots for combined AGN and combined PSR

Figures 8 and 9 show the calibration plots $\lambda_{rec}$ vs. $\lambda_{inj}$
for all AGN combined and all PSR combined respectively. For PSR, note that the
scale is not the same for $n=1$ and $n=2$. This leads to an apparent higher
value of the uncertainty.

Figure 8: Calibration plots showing $\lambda_{rec}$ vs. $\lambda_{inj}$ for
all AGN combined for the linear case (left) and the quadratic case (right)
with the J&P formalism. The light gray area corresponds to the standard
deviation of the $\lambda_{rec}$ distribution while the dark gray region shows
the statistical uncertainty. For both plots, a function $a\ \lambda_{inj}+b$
is fitted (black line).

Figure 9: Calibration plots showing $\lambda_{rec}$ vs. $\lambda_{inj}$ for
all PSR combined for the linear case (left) and the quadratic case (right)
with the J&P formalism. The light gray area corresponds to the standard
deviation of the $\lambda_{rec}$ distribution while the dark gray region shows
the statistical uncertainty. For both plots, a function $a\ \lambda_{inj}+b$
is fitted (black line). Note that the scales are not the same for $n=1$ and
$n=2$.

##  Appendix C Author contributions

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

Software development activities were shared between S. Caroff, A. Gent, D.
Kerszberg, C. Levy, T. Lin, L. Nogués, C. Perennes and T. Terzić. M. Ronco
contributed studying J&P and DSR redshift dependence of the time delays, as
well as in the writing of the introduction. J. Bolmont, S. Caroff, M. Gaug, A.
Gent, D. Kerszberg, C. Levy, T. Lin, T. Terzić provided the IRFs used in the
paper. S. Caroff and C. Levy were responsible for producing the final results
and plots. Finally, all the authors had a significant contribution writing and
editing the draft.

## References

  * Aad et al. (2016) Aad, G., Abbott, B., Abdallah, J., et al. (ATLAS Collaboration) 2016, J. High Energ. Phys. 2016, 26, https://arxiv.org/abs/1512.02586
  * Abdalla et al. (2019) Abdalla, H., et al. (H.E.S.S. Collaboration) 2019, ApJ, 870, 93, https://arxiv.org/abs/1901.05209
  * Abramowski et al. (2011) Abramowski, A., et al. (H.E.S.S. Collaboration) 2011, Astropart. Phys., 34, 738, https://arxiv.org/abs/1101.3650
  * Abramowski et al. (2015) Abramowski, A., et al. (H.E.S.S. Collaboration) 2015, ApJ, 802, 65, https://arxiv.org/abs/1501.05087
  * Acciari et al. (2019a) Acciari, V. A., et al. (MAGIC Collaboration) 2019a, Nature, 575, 455, https://arxiv.org/abs/2006.07249
  * Acciari et al. (2019b) Acciari, V. A., et al. (MAGIC Collaboration) 2019b, Nature, 575, 459, https://arxiv.org/abs/2006.07251
  * Acciari et al. (2020) Acciari, V. A., et al. (MAGIC Collaboration) 2020, Phys. Rev. Lett., 125, 021301, https://arxiv.org/abs/2001.09728
  * Acharya et al. (2019) Acharya, B. S., et al. (CTA Consortium) 2019, World Scientific Publishing Co. Pte. Ltd., https://arxiv.org/abs/1709.07997
  * Addazi et al. (2021) Addazi, A., et al. (COST Action CA18108), Prog. Part. Nucl. Phys., in press, https://arxiv.org/abs/2111.05659
  * Ade et al. (2016) Ade, P. A. R., et al. (Planck Collaboration) 2016, A&A, 594, A13, https://arxiv.org/abs/1502.01589
  * Aghanim et al. (2018) Aghanim, N., et al. (Planck Collaboration) 2020, A&A, 641, A6, https://arxiv.org/abs/1807.06209, [erratum: A&A, 652 (2021), C4]
  * Aharonian et al. (2006) Aharonian, F., et al. (H.E.S.S. Collaboration) 2006, A&A, 457, 899, https://arxiv.org/abs/astro-ph/0607333
  * Aharonian et al. (2007) Aharonian, F., et al. (H.E.S.S. Collaboration) 2007, ApJ, 664, 71, https://arxiv.org/abs/0706.0797
  * Ahnen et al. (2017) Ahnen, M. L., et al. (MAGIC Collaboration) 2017, ApJS, 232, 9, https://arxiv.org/abs/1709.00346
  * Ajello et al. (2019) Ajello, M., et al. (Fermi-LAT Collaboration) 2019, ApJ, 878, 52, https://arxiv.org/abs/1906.11403
  * Albert et al. (2007) Albert, J., et al. (MAGIC Collaboration) 2007, ApJ, 669, 862, https://arxiv.org/abs/astro-ph/0702008
  * Aleksić et al. (2012) Aleksić J., et al. (MAGIC Collaboration) 2012, Astropart. Phys., 35, 435, https://arxiv.org/abs/1108.1477
  * Aleksić et al. (2016) Aleksić J., et al. (MAGIC Collaboration) 2016, Astropart. Phys., 72, 76, https://arxiv.org/abs/1409.5594
  * Aleksić et al. (2015) Aleksić J., et al. (MAGIC Collaboration) 2015, JHEAp, 5, 30, https://arxiv.org/abs/1406.6892
  * Alfaro et al. (1999) Alfaro, J. 1999, Phys. Rev. Lett., 84, 2318, https://arxiv.org/abs/hep-th/0412295
  * Ambjørn et al. (2012) Ambjørn, J., et al. 2012, Phys. Rep., 519, 127, https://arxiv.org/abs/1203.3591
  * Amelino-Camelia et al. (1998) Amelino-Camelia, G., et al. 1998, Nature, 393, 763, https://arxiv.org/abs/astro-ph/9712103
  * Amelino-Camelia & Majid (2000) Amelino-Camelia, G., & Majid, S. 2000, Int. J. Mod. Phys. A, 15, 4301, https://arxiv.org/abs/hep-th/9907110
  * Amelino-Camelia (2002) Amelino-Camelia, G. 2002, Nature, 418, 34, https://arxiv.org/abs/gr-qc/0207049
  * Amelino-Camelia (2013) Amelino-Camelia, G. 2013, Living Rev. Relativ., 16, 5, https://link.springer.com/article/10.12942/lrr-2013-5
  * Amelino-Camelia et al. (2017) Amelino-Camelia, G., Calcagni, G., Ronco, M. 2017, Phys. Lett. B, 774, 630, https://arxiv.org/abs/1705.04876
  * Amelino-Camelia et al. (2021) Amelino-Camelia, G., Rosati, G., Bedić, S. 2021, Phys. Lett. B, 820, 136595, https://arxiv.org/abs/2012.07790
  * Antoniadis et al. (1998) Antoniadis, I., et al. 1998, Phys. Lett. B, 436, 257, https://arxiv.org/abs/hep-ph/9804398
  * Barcaroli et al. (2016) Barcaroli, L., et al. 2016, Phys. Rev. D, 95, 024036, https://arxiv.org/abs/1612.01390
  * Barrau et al. (2015) Barrau, A., et al. 2015, JCAP, 05, 051, https://arxiv.org/abs/1404.1018
  * Barres de Almeida & Daniel (2012) Barres de Almeida, U., & Daniel, M. K. 2012, Astropart. Phys., 35, 850, https://arxiv.org/abs/1204.2205
  * Bernardini et al. (2017) Bernardini, M. G., et al. 2017, A&A, 607, A121, https://arxiv.org/abs/1710.08432
  * Biteau & Williams (2015) Biteau, J., & Williams, D. A. 2015, ApJ, 812, 60, https://arxiv.org/abs/1502.04166
  * Bolmont et al. (2008) Bolmont, J., et al. 2008, ApJ, 676, 532, https://arxiv.org/abs/astro-ph/0603725
  * Bombelli et al. (1987) Bombelli, L., et al. 1987, Phys. Rev. Lett., 59, 521
  * Brahma & Ronco (2018) Brahma, S., & Ronco, M. 2018, Phys. Lett. B, 778, 184, https://arxiv.org/abs/1801.09417
  * Brun & Rademakers (1997) Brun, R., & Rademakers, F. 1997, Nucl. Inst. & Meth. in Phys. Res. A, 389, 81
  * Calcagni (2017) Calcagni, G. 2017, J. High Energ. Phys. 2017, 138, https://arxiv.org/abs/1612.05632
  * Calcagni et al. (2019) Calcagni, G., et al. 2019, JCAP, 10, 012, https://arxiv.org/abs/1907.02489
  * Caraveo et al. (2001) Caraveo, P.A., et al. 2001, ApJ, 561, 930, https://arxiv.org/abs/astro-ph/0107282
  * Carroll et al. (2001) Carroll, S. M., et al. 2001, Phys. Rev. Lett., 87, 141601, https://arxiv.org/abs/hep-th/0105082
  * Castro-Tirado et al. (2019) Castro-Tirado, A. J., et al. 2019, GCN Circulars 23708, https://gcn.gsfc.nasa.gov/gcn3/23708.gcn3
  * Chretien et al. (2015) Chretien, M., et al. 2015, Proceedings of ICRC 2015, https://arxiv.org/abs/1509.03545
  * Chretien (2015) Chretien, M. 2015, PhD thesis, https://tel.archives-ouvertes.fr/tel-01271357
  * Ellis et al. (2006) Ellis, J., et al. 2006, Astropart. Phys., 25, 402, https://arxiv.org/abs/astro-ph/0510172
  * Ellis et al. (2008) Ellis, J., et al. 2008, Astropart. Phys., 29, 158, https://arxiv.org/abs/0712.2781
  * Ellis et al. (2019) Ellis, J., et al. 2019, Phys. Rev. D, 99, 083009, https://arxiv.org/abs/1807.00189
  * Gambini & Pullin (1999) Gambini, R., & Pullin, J. 1999, Phys. Rev. D, 59, 124021, https://arxiv.org/abs/gr-qc/9809038
  * Ganguly et al. (2013) Ganguly, R., et al. 2013, MNRAS, 435, 1233, https://arxiv.org/abs/1307.6799
  * Götz et al. (2014) Götz, D., et al. 2014, MNRAS, 444, 2776, https://arxiv.org/abs/1408.4121
  * Gropp et al. (2019) Gropp, J. D., et al. 2019, GCN Circulars 23688, https://gcn.gsfc.nasa.gov/gcn3/23688.gcn3
  * Jacob & Piran (2008) Jacob, U., & Piran, T. 2008, JCAP, 01, 031, https://arxiv.org/abs/0712.2170
  * Jacobson et al. (2006) Jacobson, T., et al. 2006, Annals Phys., 321, 150, https://arxiv.org/abs/astro-ph/0505267
  * Kaplan et al. (2008) Kaplan, D.L., et al., 2008, ApJ, 677, 1201, https://arxiv.org/abs/0801.1142
  * Kostelecký & Mewes (2008) Kostelecký, V.A., & Mewes, M. 2008, ApJ, 689, L1–L4, https://arxiv.org/abs/0809.2846
  * Kowalski-Glikman & Nowak (2002) Kowalski-Glikma, J., & Nowak, S. 2002, Phys. Lett. B, 539, 126, https://arxiv.org/abs/hep-th/0203040
  * Levy et al. (2021) Levy, C., Sol, H., Bolmont, J. 2021, Proceedings of ICRC 2021, https://arxiv.org/abs/2110.06734
  * Liberati (2013) Liberati, S. 2013, Class. Quant. Grav., 30, 133001, https://arxiv.org/abs/1304.5795
  * Livine et al. (2011) Livine, E., et al. 2011, Class. Quant. Grav., 28, 245010, https://arxiv.org/abs/1104.5509
  * Magueijo & Smolin (2004) Magueijo , J., & Smolin, L. 2004, Class. Quant. Grav., 21, 1725, https://arxiv.org/abs/gr-qc/0305055
  * Mao (2011) Mao, L.S., 2011, New Astronomy, 16, 503
  * Martinez & Errando (2009) Martínez, M. & Errando, M. 2009, Astropart. Phys., 31, 226, https://arxiv.org/abs/0803.2120
  * Mattingly (2005) Mattingly, D. 2005, Living Rev. Rel., 8, 5, https://link.springer.com/article/10.12942/lrr-2005-5
  * Mavromatos (2010) Mavromatos, N. E. 2010, IJMPA, 25, 5409, https://arxiv.org/abs/1010.5354
  * Niedermaier et al. (2006) Niedermaier, M., et al. 2006, Living Rev. Relat., 9, 5, https://link.springer.com/article/10.12942/lrr-2006-5
  * Oriti (2001) Oriti, D. 2001, Rept. Prog. Phys., 64 1703, https://arxiv.org/abs/gr-qc/0106091
  * Perennes et al. (2020) Perennes, C. , Sol, H. and Bolmont, J. 2020, A&A, 633, A143, https://arxiv.org/abs/1911.10377
  * Pfeifer (2018) Pfeifer, C. 2018, Phys. Lett. B, 780, 246, https://arxiv.org/abs/1802.00058
  * Polyakov (1981) Polyakov, A. M. 1981, Phys. Lett. B, 103, 207
  * Pueschel (2019) Pueschel, E. 2019, Proceedings of ICRC 2019, https://arxiv.org/abs/1908.04163
  * Rosati et al. (2015) Rosati, G., et al. 2015, Phys. Rev. D, 92, 124042, https://arxiv.org/abs/1507.02056
  * Rovelli (2007) Rovelli, C. 2007, Cambridge University Press, Cambridge U.K.
  * Selsing et al. (2019) Selsing, J., et al. 2019, GCN Circulars 23695, https://gcn.gsfc.nasa.gov/gcn3/23695.gcn3
  * Terzić et al. (2021) Terzić. T., Kerszberg, D., Strišković, J. 2021, Universe, 7, 345, https://arxiv.org/abs/2109.09072
  * Vasileiou et al. (2013) Vasileiou, V., et al. 2013, Phys. Rev. D, 87, 122001, https://arxiv.org/abs/1305.3463
  * Vasileiou et al. (2015) Vasileiou, V., et al. 2015, Nature, 11, 344
  * Xu & Ma (2016a) Xu, H., & Ma, B.-Q. 2016, Astropart. Phys., 82, 72, https://arxiv.org/abs/1607.03203
  * Xu & Ma (2016b) Xu, H., & Ma, B.-Q. 2016, Phys. Lett. B, 760, 602, https://arxiv.org/abs/1607.08043
  * Xu & Ma (2018) Xu, H., & Ma, B.-Q. 2018, JCAP, 1801, 050, https://arxiv.org/abs/1801.08084
  * Zitzer et al. (2013) Zitzer, B., et al. 2013, Proceedings of ICRC 2013, https://arxiv.org/abs/1307.8382
  * Zyla et al. (2020) Zyla, P.A., et al. (Particle Data Group) 2020, Prog. Theor. Exp. Phys. 2020, 083C01, https://pdg.lbl.gov/rpp-archive/index-2021.html


* [52] Hadžić M, Jang J. Nonlinear Stability of Expanding Star Solutions of the Radially Symmetric Mass‐Critical Euler‐Poisson System. Communications on Pure and Applied Mathematics. 2018 May;71(5):827-91.
  * [53]
  * [54] Makino T. Blowing up solutions of the Euler-Poisson equation for the evolution of gaseous stars. Transport Theory and Statistical Physics. 1992 Aug 1;21(4-6):615-24.
  * [55]
  * [56] Ducomet, B., 1996. Hydrodynamical models of gaseous stars. Reviews in Mathematical Physics, 8(07), pp.957-1000.
  * [57] item Deng Y, Liu TP, Yang T, Yao ZA. Solutions of Euler-Poisson Equations¶ for Gaseous Stars. Archive for rational mechanics and analysis. 2002 Sep;164(3):261-85.
  * [58]
  * [59] Guo Y, Hadžić M, Jang J. Continued gravitational collapse for Newtonian stars. Archive for Rational Mechanics and Analysis. 2021 Jan;239(1):431-552.
  * [60]
  * [61] Hadžić M, Jang J. Nonlinear Stability of Expanding Star Solutions of the Radially Symmetric Mass‐Critical Euler‐Poisson System. Communications on Pure and Applied Mathematics. 2018 May;71(5):827-91.
  * [62]
  * [63] Yuen M. Stabilities for Euler–Poisson equations in some special dimensions. Journal of mathematical analysis and applications. 2008 Aug 1;344(1):145-56.
  * [64]
  * [65] Deng Y, Xiang J, Yang T. Blowup phenomena of solutions to Euler–Poisson equations. Journal of mathematical analysis and applications. 2003 Oct
  * [66]
  * [67] Shapiro SL, Teukolsky SA. Black holes, white dwarfs, and neutron stars: The physics of compact objects. John Wiley and Sons; 2008 Nov 20.
  * [68]
  * [69] Oppenheimer JR, Volkoff GM. On massive neutron cores. Physical Review. 1939 Feb 15;55(4):374.
  * [70]
  * [71] Makino T. On spherically symmetric stellar models in general relativity. Journal of Mathematics of Kyoto University. 1998 Jan;38(1):55-69.
  * [72]
  * [73] Boonserm P, Visser M, Weinfurtner S. Solution generating theorems for the Tolman-Oppenheimer-Volkov equation. Physical Review D. 2007 Aug 30;76(4):044024.
  * [74]
  * [75] Wald RM. General relativity. University of Chicago press; 2010 May 15.
  * [76]
  * [77] Bardeen JM, Carter B, Hawking SW. The four laws of black hole mechanics. Communications in mathematical physics. 1973 Jun;31(2):161-70.
  * [78]
  * [79] Tolman RC. Relativity, thermodynamics, and cosmology. Courier Corporation; 1987.
  * [80]
  * [81] Cocke WJ. A maximum entropy principle in general relativity and the stability of fluid spheres. In Annales de l’IHP Physique théorique 1965 (Vol. 2, No. 4, pp. 283-306).
  * [82]
  * [83] Gao S. General maximum entropy principle for self-gravitating perfect fluid. Physical Review D. 2011 Nov 9;84(10):104023.
  * [84]
  * [85] Sorkin RD, Wald RM, Jiu ZZ. Entropy of self-gravitating radiation. General Relativity and Gravitation. 1981 Dec;13(12):1127-46.
  * [86]
  * [87] Roupas Z. Thermodynamical instabilities of perfect fluid spheres in General Relativity. Classical and Quantum Gravity. 2013 May 13;30(11):115018.
  * [88]