Probing Quantum Gravity with Imaging Atmospheric Cherenkov Telescopes

High energy photons from astrophysical sources are unique probes for some predictions of candidate theories of Quantum Gravty (QG). In particular, imaging atmospheric Cherenkov telescopes (IACTs) are instruments optimised for astronomical observations in the energy range spanning from a few tens of GeV to ~100 TeV, which makes them excellent instruments to search for effects of QG. In this article, we will review QG effects which can be tested with IACTs, most notably the Lorentz invariance violation and its consequences. It is often represented and modelled with photon dispersion relation modified by introducing energy-dependent terms. We will describe the analysis methods employed in the different studies, allowing for careful discussion and comparison of the results obtained with IACTs for more than two decades. Loosely following historical development of the field, we will observe how the analysis methods were refined and improved over time, and analyse why some studies were more sensitive than others. Finally, we will discuss the future of the field, presenting ideas for improving the analysis sensitivity and directions in which the research could develop.


Introduction and Motivation
The general theory of relativity is a beautiful and elegant theory, which connects the local matter and energy content to the curvature of spacetime, thus giving a classical description of gravity. It has been heavily tested and scrutinised ever since Albert Einstein proposed it in 1915 [1]. Nevertheless, it breaks down in extreme circumstances such as singularities within black holes, or the early universe. Therefore, it is expected that there exists a more fundamental quantum theory of gravity, which can handle these extreme situations. Furthermore, a quantum theory of gravity would be a giant leap towards unification of all fundamental forces.
Theoretical endeavours in formulating the theory of QG have explored different avenues (see, e.g., [2][3][4][5][6][7]). However, despite significant efforts, a complete and consistent description of gravity on a quantum level remains unknown. In addition, many of the QG models include departures form the Lorentz symmetry (see, e.g., [8][9][10][11]). Performing measurements in the expected realm of QG would strongly hint in which direction theoretical research should proceed. Unfortunately, the expected domain of QG is the Planck scale 1 . Even if QG emerges at energies several orders of magnitude below E Pl , it is still vastly above the highest energies accessible in contemporary human-built accelerators. When technology falls short, we turn to nature's own accelerators: active galactic nuclei, gamma-ray bursts, supernova remnants, pulsars, etc. The most energetic particles detected up to date, a cosmic ray at ∼3.2 × 10 11 GeV [12,13], a neutrino at ∼6.3 × 10 6 GeV [14], and a gamma ray at ∼1.4 × 10 6 GeV [15], while reaching energies higher than those achievable in Earth-based accelerators, are still more than a little shy of E Pl . So what allows us to hope that we will measure an effect of QG? 1 Planck energy E Pl = √h c 5 /G ≈ 1.22 × 10 19 GeV, Planck length L Pl = √h G/c 3 ≈ 1.62 × 10 −35 m, Planck time t Pl = √h G/c 5 ≈ 5.39 × 10 −44 s.

A Proposal to Probe Quantum Gravity
In 1997, the distance of a gamma-ray burst (GRB) 2 was measured for the first time. Indeed, following the detection of GRB 970508 [16], an optical counterpart was observed [17], which allowed estimating of its redshift to 0.835 ≤ z 2.3 [18]. A strong flux of gamma rays from a quickly varying source detected at a cosmological distance incited Amelino-Camelia et al. [19] to suggest that the signal from GRBs could be used to probe the structure of spacetime. The proposal was based on the idea of spacetime as a dynamical medium, which experiences quantum fluctuations due to QG effects. While the scale of fluctuations was expected to be comparable to the Planck units, propagation of photons of substantially smaller energies could still be affected by it. Probing the fluctuations would result in an energy-dependent propagation speed, similar to what visible light experiences when propagating through a medium such as water or air.

Modified Photon Dispersion Relation
This behaviour can be modelled by modifying the standard photon dispersion relation in the following way: where E and p are respectively the energy and momentum of a photon, c is the Lorentz invariant speed of light, and the E QG,n are the energy scales at which effects of QG become significant. We will start discussing the values of E QG,n in a short while, for now let us just acknowledge that E/E QG,n 1 even for the most energetic gamma rays. Different modifying terms in the dispersion relation contribute less and less with increasing n. Therefore, usually, only the first two leading terms (n = 1 or n = 2) of the series are considered and independently tested for 3 . They are often referred to as linear and quadratic energy-dependent contributions, respectively. Letting E QG,n → ∞ for all n leads to the well-known Lorentz invariant photon dispersion relation. Parameters S n can take values ±1, and their role will become apparent immediately. From Equation (1) the energy-dependent photon group velocity can be easily derived as: ( Considering each modifying term independently, one can see that for S n = +1, the velocity becomes greater than c, while for S n = −1 it becomes smaller than c. These two behaviours are known as superluminal and subluminal, respectively. Once a modification is introduced in the dispersion relation, various effects (other than changes in the photon speed) are conceivable, such as the modification of the electromagnetic interaction. But whatever the effects of modifying the dispersion relation may be, they are minuscule because of the ratio E/E QG,n . The good news, however, is that the effects are cumulative. This is extremely important because gamma rays from some astrophysical sources take billions of years to reach Earth, allowing for these potential effects to accumulate, thus giving hope that we might be able to measure their consequences from Earth.
Additionally, the effects of modifying the dispersion relation are more pronounced for higher energy photons. Thus, searching for them with IACTs, which are instruments optimised for astronomical observations in the very high energy (VHE) gamma-ray band (100 GeV < E < 100 TeV), is a sensible thing to do. Given their large collection area and good sensitivity, IACTs are excellent instruments for testing effects of QG on gamma rays, and will be the main focus of this review. At lower energies, satellite-born detectors such as the Fermi-Large Area Telescope (LAT) benefit from more distant observations, but suffer from their lower effective area (see Sections 2.13 and 4 for a brief comparison with the IACT results). On the other hand, at higher energies water Cherenkov detectors such as High Altitude Water Cherenkov (HAWC) or Large High Altitude Air Shower Observatory (LHAASO) have an advantage of observing in a higher energy range than IACTs. However, due to the rapid decrease of the flux of gamma rays at these energies they are handicapped by smaller statistics, which makes them less sensitive to fast flux variations (see Sections 3.3 and 4 for a brief comparison with the IACT results).
Before taking a dive into the methods and results of probing QG with IACTs, let us acknowledge that the modified photon dispersion relation is the usual starting point of experimental tests of QG on gamma rays. While some QG models indeed do not preserve Lorentz symmetry, it is important to note that Equation (1) is not a direct consequence of any particular QG model. Given that there is no fully formulated theory of QG, it would be overambitious to expect exact predictions. Rather, the modified dispersion relation can be regarded as a simple way of parameterizing and modelling phenomena not predicted by the current physical theories and laws. It, therefore, enables us to experimentally search for effects of those phenomena.
That being said, there are two main ways of modifying the dispersion relation that are usually considered. LIV, the main focus of this review, implies the existence of a preferred inertial frame of reference, which breaks the Lorentz symmetry [21]. However, there are also ways of modifying the photon dispersion relation, while at the same time preserving the Lorentz symmetry. One example is the so-called Doubly Special Relativity (DSR) [24,25]. In this model, the symmetry is deformed, rather than broken, and there is no preferred inertial frame of reference. Moreover, in order to keep the conservation laws covariant with respect to deformed symmetries of DSR, the conservation laws themselves need to be modified. This fundamental difference between LIV and DSR becomes important when different possible effects of QG are discussed. In particular, the kinematics, and possibly the dynamics, of electromagnetic interactions in the LIV framework will differ from the Lorentz invariant ones. In the DSR framework, on the other hand, the descriptions of interactions will be the same as (or only slightly different from) the Lorentz invariant descriptions. DSR is a recently discussed promising avenue of research gaining attention and traction. However, there has been no published results from IACTs mentioning explicitly DSR effects thus far. Therefore, in the rest of the text, we will, refer to all effects as LIV, regardless of their true origin. We will, however, keep using E QG,n to note the energy scales at which the effects become relevant. The details of either of these models, and their differences are out of the scope of this work. An interested reader is referred to a review paper by the COST Action 18108 4 (in preparation) and references therein.
In this paper, we will focus on searches for signatures of LIV in measurements with IACTs. We will discuss various effects of modifying the photon dispersion relation and their respective probes, adopting a chronological course. However, it is our intention (instead of simply recalling the most important studies performed) to analyse the evolution of the field, with a particular focus on the development of the analysis methods. Hopefully, this approach will inspire the authors and the readers alike to formulate new ideas on how to search for the effects of QG, and pave the path for future research. Historically, the first effect to be tested was the energy dependence of the photon group velocity, so results of different measurements of the photon time of flight will be covered first, in Section 2. As stated above, LIV can affect the kinematics and dynamics of electromagnetic interactions. This other important class of effects will be discussed in Section 3. We might as well break the suspense and state right away that no effects of QG have been detected so far. Nevertheless, strong constraints have been set on the minimum value of the LIV energy scale. These are usually expressed as lower limits at the 95% confidence level. The results of different effects, obtained from various experiments and analysis methods will be mutually compared and their differences discussed in Section 4. Finally, we turn towards the future in Section 5, to discuss opportunities for development and progress of this field of research.

Testing Energy-Dependent Photon Group Velocity
Assuming energy-dependent propagation speeds, two photons of energies E 2 > E 1 emitted from a source at the same time will have different, energy-dependent, times of flight t 2 and t 1 respectively, finally reaching Earth with an energy-dependent time delay [26]: where t is the time needed for a photon travelling with speed c to reach the Earth 5 . The time delay is proportional to a source distance parameter: where z s is the source redshift, and H 0 , Ω m , and Ω Λ represent cosmological parameters, respectively: the Hubble constant, the matter density parameter, and the dark-energy density parameter 6 . The time delay expression was derived from comoving trajectories of particles, starting from their modified dispersion relations. More general and alternative expressions can be obtained by modifying the general relativistic dispersion relation as was done in [27], or by adopting that the spacetime translations are modified alongside with the modification of the dispersion relation [28].

The First Test with an Imaging Atmospheric Cherenkov Telescope
Soon after it was proposed that GRBs could be used to search for effects of LIV, the first test using data from IACTs was performed. In fact, researchers observing with the Whipple telescope 7 already had a suitable data set available. Albeit, the source was not a GRB, but the very first active galactic nucleus (AGN) ever detected in the VHE gamma-ray band, Markarian 421 (Mrk 421, redshift z = 0.031) [30]. On 15 May 1996, Whipple observed the most rapid flare from Mrk 421 up to that time, with the flux doubling time of less than 15 5 Note that the time delay can be both positive or negative, depending on whether the behaviour is subluminal or superluminal, respectively. According to the usual convention, a time delay is positive for subluminal behaviour, i.e., photon of a higher energy propagating at a lower speed than a lower energy photon. 6 Various LIV studies use different values for cosmological parameters. However, given the precision of these studies, their final results are not strongly affected by the differences in the values of the cosmological parameters used, so we will treat them equally in that respect. 7 The Whipple telescope (https://veritas.sao.arizona.edu/whipple, accessed 15 July 2021) was the first IACT. It consisted of a single 10 m reflector dish. In operation from 1968 until 2013, it detected the very first TeV gamma-ray source, the Crab nebula [29], and the first AGN detected in the same energy range Mrk 421 [30]. minutes and including photons of energies up to several TeV [31]. This groundbreaking study used a rather rudimentary analysis: the data set was split in two energy bands (E < 1 TeV and E > 2 TeV) [32]. In each energy band, the events were further subdivided in time bins of 280 s. The distribution of arrival times of photons with energies E > 2 TeV was compared to the distribution of arrival times of events below 1 TeV. The authors used the likelihood-ratio test 8 to compare the contents of time bins in the two energy ranges. In this study, no distinction was made between the subluminal and the superluminal behaviour. No delay in either direction was detected at the 95% confidence level. Combined with the distance to Mrk 421, this result was translated into a lower limit on the LIV energy scale E QG,1 > 4 × 10 16 GeV. Only the linear contribution was considered in this first study.

Fastest Variability in Blazars
The observation of Mrk 421 with the Whipple telescope drew attention to flaring blazars as possible probes of LIV. In the summer of 2005, the Major Atmospheric Gamma Imaging Cherenkov (MAGIC) telescopes 9 observed two flares from the AGN Markarian 501 (Mrk 501, redshift z = 0.034) [35]. The data analysis revealed that the flux doubled in only 2 minutes, which remains until today the fastest flux variability ever observed from a blazar in the VHE gamma-ray band. With the highest energies reaching ∼10 TeV, the flux varying by an order of magnitude, it was a chance not to be missed. Moreover, there was an indication of a 4 ± 1 min delay between the peaks in the light curves in the lowest (0.15-0.25 TeV) and the highest (1.2-10 TeV) energy bins on the 9th of July. A search for an energy-dependent photon time of flight in the Mrk 501 flare of this night was performed employing two distinct statistical analysis methods which we will now describe.
Energy cost function (ECF) method utilises the fact that a signal pulse propagating through a dispersive medium will be diluted, and its power (total energy per unit time), consequently, decreased (see, e.g., Section 7.9 in [36]). In the case of the Mrk 501 flare [37], the data sample was chosen by selecting the most active part of the flare, i.e., the time interval in which the temporal distribution of events differs the most from a uniform distribution. We will mark the beginning and the end of this time interval as t 1 and t 2 , respectively. The power of the signal was calculated as the sum of the energies of all the photons within the interval divided by the duration of the interval. Had the photons experienced any energy-dependent time delay, the power would have been smaller than without dispersion. One can then search for the maximal possible power by applying dispersion in the opposite direction, assuming different values of the LIV energy scale. Specifically, in order to compute a new signal power, a new arrival time t i was calculated for each photon in the sample for a particular value of E QG,n : where t i and E i are, respectively, the measured arrival time and the measured energy of the i-th photon. Given the large values of E QG,n , various parameters are often introduced to facilitate numerical computations. Here the parameters η n is defined from Equation (3) as: This parameter is introduced for computational reasons because η n is O(1). Usually expressed in units of [s/GeV] (for n = 1) or [s/GeV 2 ] (for n = 2), η n indicates how much a photon will be delayed in arrival compared to a photon propagating at c per every GeV of its energy. Limits on E QG,n are then derived by inverting Equation (6).
The arrival time recalculation as described in Equation (5) was performed for each individual gamma ray, and only photons whose recalculated arrival times fell in the [t 1 , t 2 ] time interval were retained. In this way, an alternative sample of photons was constituted, and its total energy calculated. The procedure was repeated for different values of η n (i.e., different values of E QG,n ). The ECF was defined as the total energy as a function of η n . The value of η n which maximises the ECF, would recover the maximal signal power. In other words, it would correspond to the measurement of the dispersion which the gamma rays experienced because of the LIV effects, assuming no other effects play a significant role. The sensitivity of the method and the confidence interval for parameters η n were estimated using Monte Carlo simulations of the observed signal. Next, 1000 of simulated data sets were generated, and the ECF method was applied to each of them. The most probable value of η n and its confidence interval were estimated from the distribution of η n , which were then translated into lower limit on the LIV energy scale at the 95% confidence level. This particular analysis yielded E QG,1 > 2.1 × 10 17 GeV, which was more constraining than the Whipple result on Mrk 421 flare from 1996 [32] by an order of magnitude. In addition, for the first time the quadratic contribution was constrained, setting E QG,2 > 2.6 × 10 10 GeV. Unlike the approach used for the Mrk 421 data analysis, the ECF allowed for testing of superluminal, as well as subluminal behaviours. Nevertheless, only the subluminal behaviour was investigated.
There are several methods based on the idea of removing the dispersion from the data: the sharpness maximisation method (SMM) [38,39], the dispersion cancellation (DisCan) [40], and the minimal dispersion (MD) [41]; the main difference between these approaches being the way the sharpness of the light curve is quantified. We will investigate in more details a variation of the DisCan method in Section 2.8 and the SMM in Section 2.13.

Introducing the Maximum Likelihood Method
Originally, the maximum likelihood (ML) method was proposed for the analysis of the Mrk 501 data set described in the previous section. However, as we shall soon see, it became the standard analysis method used for searches of energy-dependent gamma-ray group velocity in IACTs data. Therefore, we will dedicate a separate section to its description. Introduced by Martínez & Errando in [42], the authors argued that the analysis methods used should be unbinned (unlike the one previously employed in the case of Mrk 421 [32]), in order to fully exploit the information carried by a relatively small gamma-ray sample. The ECF method, used in [37], was indeed unbinned, however, it depended upon identifying and isolating the flares from the rest of the light curve. While this particular Mrk 501 light curve from 2005 had a relatively simple structure [35], it was already recognised that the ECF method would not be suited for the analysis of complex light curves, or segments of flares. Therefore, the unbinned ML method soon became a standard approach in searches for energy-dependent time delays, with every new study incorporating additional features and improvements. Here we will depart from the historical course, and describe the ML method in its present form.
In order to search for LIV, the ML method makes use of a profile likelihood ratio test: where η n is the LIV parameter of order n of interest, ν represents the nuisance parameters, η n andν are values that maximize the likelihood L,ν maximizes the likelihood L for a given η n , and D represents the observed data on which the analysis is performed. According to Wilks' theorem [43], the distribution of −2 ln λ p (η n |D) follows a χ 2 distribution with 1 degree of freedom for the true value of η n , i.e., the one we are looking for. The 95% confidence level one-sided upper limits are therefore derived by solving the following equation: while 95% confidence level two-sided upper limits are obtained using: In the case where the conditions for Wilks' theorem are not fulfilled, one can calibrate intervals using Monte Carlo simulated samples of the null hypothesis. The right value for any particular case can then be derived from the quantiles of the distribution of these simulations. For instance, the 95% two-sided confidence interval is delimited by the lower and upper 2.5% quantiles.
The likelihood function L, for an observed number of events N ON , can be written as: • G(E, E true ), contains the information about the energy resolution and the bias of the instrument. E true is the true energy of a particular event, and G(E, E true ) is the PDF of E true being measured as E.

•
The final ingredient, A eff (E true , t) represents the collection area (i.e., acceptance) of the instrument expressed in true energy E true . In the most general case, it can change with time, especially if the data were collected in different observation conditions.
The PDF for observing a background event of reconstructed energy E i at the moment t i has fundamentally the same form as the PDF for signal events, see Equation (11). However, the origin of background events is generally not known, so time of flights of individual events (whether affected by LIV or not) cannot be determined. Therefore, both temporal and energy distributions of events are taken as observed on Earth. Concretely, in Equation (11), when used for background events η n = 0, and F(t) and Φ obs (E) are the measured background light curve and spectrum on Earth. The final pieces of puzzle are probabilities for each event to be part of the signal or background, p i . In IACTs, the signal is estimated from a region around the source position in the field of view, usually referred to as the ON region. However, besides the signal, the ON region contains an irreducible contribution from the background. The background is estimated from the so-called OFF region, a region in the field of view which contains no sources of gamma rays, and is observed under the same conditions as the ON region. Usually, the probabilities for each event to be a part of the signal or background are calculated as follows: where N ON is the total number of events in the ON region, N OFF the total number of events in the OFF region, and α is the ratio of effective exposure times in the two: α = t ON /t OFF 10 . A legitimate objection to the ML method is that it relies on our knowledge of source-intrinsic processes, which is limited at best. In that sense, the ECF or similar methods like the SMM (see Section 2.13) have the advantage of not depending on our knowledge of source-intrinsic effects. However, it is quite imaginable that there are source-intrinsic dispersive processes, which could mimic effects of LIV 11 . These would not depend on the source redshift, and could be "filtered out" by considering sources at different redshifts. Nevertheless, combining the results of different analyses might prove tricky for ECF and related methods. The likelihood function, on the other hand, should tackle that task with relative ease, as we will discuss in Section 5.2.
Another possible source of systematic effects are secondary gamma rays, which can be produced through one of the following processes: (i) hadrons accelerated within a source interact with the surrounding electromagnetic fields to produce neutral pions, which decay into gamma rays, (ii) gamma rays emitted from the source interact with magnetic fields to produce electron-positron pairs, which can create secondary gamma rays either through annihilation, or by inverse-Compton scattering of lower-energy photons. Secondary gamma rays could create a false signal, especially in the analysis methods based on individual events, such as time of flight studies. However, since secondary gamma rays are not produced within the observed source, their origin is not necessarily on the line of sight. A significant rate of secondary gamma rays would manifest as an extended emission, so-called halo, around an otherwise point-like source. Indeed, several studies searching for gamma-ray halos have been performed, but have shown no evidence thereof (see, e.g., [46][47][48], see also [49] and references therein). Though an occasional secondary gamma ray might be mistakenly treated as the 10 For a more detailed definition and discussion of al pha, we refer the interested reader to [44]. 11 For more details on the modelling of these effects in the context of LIV, an interested reader can refer to [45]. signal, the effect should be minor, and it would diminish with an increasing size of a data sample.

Results from the Maximum Likelihood Method on the Mrk 501 Flare from 2005
In this first application of the ML method, the light curve of the Mrk 501 flare from 2005 was modelled with a Gaussian superimposed on top of a constant baseline emission from the source. A background contribution was not considered because of its negligible contribution in such a flare. The results showed η 1 (see Equation (6)) departing from zero by slightly more than 2σ, implying an energy-dependent time delay [37]. η 2 , on the other hand, was consistent with zero. Again, only the subluminal behaviour was tested for. The corresponding LIV energy scales were E QG,1 = 0.30 +0. 24 −0.10 × 10 18 GeV and E QG,2 = 0.57 +0.75 −0.19 × 10 11 GeV for the linear and quadratic contributions, respectively [42]. Martínez & Errando refrained from interpreting the results and focused on the description of the method. The nonzero time delay was instead discussed in [37]. A possibility of a bias in the ML analysis was investigated on a set of simulated Monte Carlo samples. An independent researcher simulated data sets with injected energy-dependent time delays. These data sets were blindly analysed using the ML method, which correctly reconstructed the injected delay values. It was concluded that the effect was real, although the statistical significance was too low to claim a discovery. Finally, it was concluded that the results obtained with the ECF and the ML methods were mutually consistent. Furthermore, some investigations of emission models suggest that the energy-dependent time delay could be a consequence of source intrinsic spectral variability in time, occurring either because of the acceleration of particles or the absorption of gamma rays [50,51]. In summary, the study on the Mrk 501 flare data from 2005, not only significantly tightened the constraints on the LIV energy scale, compared to the pioneering study of Whipple [32], but it also motivated the introduction of a novel analysis method, and served as a cross check between two fundamentally different analysis approaches. This data set was also the first one to be studied with two fundamentally different analysis approaches, allowing comparisons between their results.

Sensitivity to the Lorentz Invariance Violation Effects
After taking a look at the first searches for the possible signatures of LIV effects in IACTs data, this is a good place to analyse what properties a signal should have in order to be considered a good probe of such an effect. As we have already discussed in the beginning of Section 2, the more energetic photons will be more strongly affected by the LIV. Therefore, sources with spectra extending to higher energies, and with comparatively larger population of higher energy photons (colloquially called "harder spectra", as opposed to "softer spectra"), are more favourable. Furthermore, the farther the source, the more the effect will be accumulated. However, a large distance carries the caveat that VHE gamma rays are partially absorbed on the EBL (see the description of the ML method in page 8 et sec., and Section 3.1), which softens the spectra and depletes data samples of the most energetic photons. The time delay between two photons of different energies will also be more pronounced for smaller values of E QG,n , so smaller time delay means stronger constraint on the LIV energy scale. However, it is entirely possible that there are emission time delays present within sources, which could mimic or conceal LIV-induced arrival time delays. Considering our limited knowledge of emission mechanisms, the emission times cannot be precisely modelled. Instead, the emission time has to be constrained based on the flux variability timescale. Emission is more probable during periods of higher flux. However, high flux on its own is not enough. If the flux is constant, or changing monotonically, an application of a spectral dispersion, will not change the shape of the light curve. A variable light curve, on the other hand, will be smeared due to spectral dispersion. The effect will be more pronounced for stronger dispersion, and more detectable for faster changing flux. By inverting Equation (3), one can make a crude estimate of how LIV energy scale depends on the highest energies of detected photons (E max ), light curve variability timescale (t var ), and the redshift of the source (z s ). These dependencies are summarized in Table 1. Note that the power of the dependence on the redshift was numerically computed for z s up to ∼10, which is much further than what current IACTs can probe. Table 1. Dependence of E QG,n on the characteristics of the source and the sample. E max is the highest photon energy in the sample, t var is the shortest variability timescale in the light curve, and z s is the redshift of the source.
There is another parameter, not present in Equation (3), whose importance becomes apparent through the data analysis. That is the size of the sample. Its influence on E QG,n depends on the analysis method, and is difficult to estimate it the way it was done for other parameters in Table 1. The general rule, though, is simple: the more the better. More specific estimates will be discussed on particular cases.
Based on this simple analysis, three types of sources are considered to be suitable for testing of LIV on gamma rays: Pulsars can have rotation periods as short as a few milliseconds, although the ones detected with IACTs so far have periods of at least a few tens of milliseconds. The only four pulsars that have been detected with IACTs so far are the Crab pulsar [52], the Vela pulsar [53], the Geminga pulsar [54], and PSR B1706-44 [55]. Their pulsation is highly regular, which makes it predictable, and allows stacking of signal from different periods, thus increasing the detected statistics. Additionally, these four pulsars are located in the Milky Way. This relatively close proximity significantly impairs the sensitivity of LIV tests performed on pulsar data.
Gamma-ray busts are powerful transient cosmic explosions, usually associated with collapses of massive stars into black holes (long GRBs), or mergers of neutron stars (short GRBs). Their light curves are variable on timescales of a second. Unlike pulsars, GRBs are completely unpredictable. Satellite-borne detectors with a large field of view, such as Gammaray Burst Monitor (GBM) [56] and LAT [57] onboard satellite Fermi 12 on average detect one GRB almost every day [58]. However, IACTs with a rather small field of view (an order of few degrees) rely on alerts from satellite borne detectors to trigger observations. Furthermore, because of their large distances, VHE gamma rays are strongly absorbed on the EBL. For these reasons, GRBs are elusive and notoriously difficult to detect with IACTs, with only four detected to date ( [59][60][61][62]). However, due to their short variability timescales, combined with large distances, once a GRB is detected, the signal becomes a valuable asset for probing QG (see Section 2.12).
Active galactic nuclei are persistent sources at distances comparable to GRBs. During their flaring states, they emit signals abundant in VHE gamma rays, with flux variability timescales on the order of minutes. Although unpredictable, flares usually last longer than GRBs. In addition, they emit stronger fluxes, with the most energetic photons reaching higher energies. All of this makes flares from AGN easier to detect with IACTs compared to GRBs.

Lorentz Invariance Violation Study on the Most Variable Blazar Flare
While the Mrk 501 flare observed by MAGIC showed the fastest changing gamma-ray flux in blazars, it had a rather simple structure. Almost exactly one year later, on 28 July 2006, while the LIV data analysis on the Mrk 501 sample was still ongoing, another promising flare occurred. This time around it was the High Energy Stereoscopic System (H.E.S.S.) 13 that observed a flare from blazar PKS 2155-304 [64]. During an ∼85 min observation, flare with a quite complex structure was detected, variable on the scale of ∼200 s with several local minima and maxima, and with the signal to background ratio above 300. At the same time, no significant changes of spectrum were found. The highest flux reached more than 15 Crab units (C.U.) 14 above 200 GeV, and a total of more than eleven thousand gamma rays were detected, reaching the highest energies of ∼4 TeV. Moreover, PKS 2155-304 is located at a redshift of z = 0.116, more than three times larger than Mrk 421 and Mrk 501.
Several studies of energy-dependent time delay were performed using this signal. The first one, published soon after the flare was observed, used two different statistical methods, both estimating time lag between light curves in different energy ranges [65].
Modified cross correlation function (MCCF) was originally developed for timescale analysis of spectral lags, and it enables searches for time lags shorter than the temporal resolution of light curves [66]. In this case, the data were split in two energy bins: 200-800 GeV and >800 GeV, and the MCCF was used to estimate the time lag between the light curves in in these two energy ranges. The analysis resulted in the most stringent constraint on the linear contribution up to that time E QG,1 > 7.2 × 10 17 GeV; more than two times stronger limit than the one set by MAGIC on Mrk 501 data. The lower limit on the quadratic contribution, on the other hand, was set at E QG,2 > 1.4 × 10 9 GeV; more than 40 times lower than the one from Mrk 501 by MAGIC using the ML method, and almost 20 times lower than the one set with the ECF.
Continuous wavelet transform (CWT) method relies on identifying extrema in two energy bands and measuring their relative time delay. In this case, the chosen energy ranges were: 210-250 GeV and >600 GeV, and two pairs of extrema were identified. Only the constraint on the linear term was set at E QG,1 > 5.2 × 10 17 GeV, thus confirming the constraint obtained using the MCCF method.
Relying on the rule of thumb, laid out in Table 1, it was expected that the larger distance and faster flux variability of PKS 2155-304, compared to Mrk 501, would make this study more sensitive to the linear modifying term of the dispersion relation. The influence of these two variables to the quadratic contribution is somewhat smaller, because of the exponents, allowing a stronger influence of the highest gamma-ray energies. Nevertheless, it seems unlikely that a factor of 2.5 difference in the highest energies alone would result in a factor of forty difference between limits on the quadratic contribution. It is more likely that the MCCF and CWT methods do not fully exploit all of the potentials of the PKS 2155-304 data sample.
When proposing the ML method in [42], the authors were already aware of the PKS 2155-304 flare, and decided to test their method on that signal as well. Since they did not have the access to the actual data set, and the method relied on individual events, they generated Monte Carlo simulated data sets, based on the published information on the PKS 2155-304 flare. It was estimated that the application of the ML method on the PKS 2155-304 flare sample would be more than six times more sensitive to E QG,n compared to the Mrk 501 case. Moreover, the authors analysed where the sixfold improvement came from, and came up with similar conclusions as we have just discussed: (i) the higher redshift contributed a factor of three, (ii) larger sample of PKS 2155-304, albeit with the highest energies lower than in the Mrk 501 sample, added another factor of two, and (iii) more complex light curve shape was responsible for an additional factor. However, the authors also noted that it was in fact the fastest single change of flux, i.e., the fastest rise time or fall time in the entire light curve, which dominated the sensitivity.
Following the study by Martínez & Errando [42], the H.E.S.S. Collaboration performed another search for effects of LIV in the PKS 2155-304 flare data, this time fully adopting the ML method [67]. For this occasion, a particular H.E.S.S. data analysis was performed, focusing on the initial 4000 s of the observation, during which both the flux and its variability were the highest. Upon applying some additional cuts on the data, only 3526 events remained (out of more than 11,000 in the original data set) in the 0.25-4.0 TeV energy range. This resulted in a strong background suppression, and a very good fit of the light curve and the spectrum. Based on optimisation using Monte Carlo simulations, the data were finally separated in two energy bins: 0.25-0.28 TeV and 0.3-4.0 TeV. The lower bin was used to create the light curve template. The data were fitted with a sum of a constant baseline emission and five consecutive asymmetric Gaussian curves. The events from the higher energy bin were used to calculate the likelihood. The results were E QG,1 > 2.1 × 10 18 GeV and E QG,2 > 6.4 × 10 10 GeV for the linear and quadratic term, respectively, both significantly more constraining than the ones obtained in the previous analysis by H.E.S.S. using MCCF, demonstrating the dominance of the ML method on a concrete case. Furthermore, both results were in line with the assessments by Martínez & Errando, and finally, both were the most constraining lower limits on the LIV energy scale up to that time. Discussing their results, the authors reached similar conclusions as Martínez & Errando in their work. In particular, the higher sensitivity was due to the high flux variability and large data sample, while the lower maximal energies somewhat impaired the sensitivity. Furthermore, the uncertainty on the estimated parameter depended mostly on the width of the individual flux peaks, which was in agreement with the conclusion by Martínez & Errando that the sensitivity is dominated by the fastest single change of flux. Final important point was that the estimated parameter uncertainty only mildly depended on the number of events used to calculate the likelihood, meaning that robust results are obtainable even with small data sets.

Extending to Higher Redshifts
On 26 and 27 April 2012, the H.E.S.S. telescopes observed a flare from the blazar PG 1553+113 [68]. The flux was three times higher than the archival measurements, with an indication of intra-night variability. Interestingly, the redshift of the source had been only loosely constrained prior to this study. In order to estimate the redshift more precisely, the authors devised a method based on Bayesian statistics, which relies on accounting for the absorption of VHE gamma rays on the EBL 15 . This enabled them to estimate the redshift to be z = 0.49 ± 0.04. Though the flux showed only a hint of intra-night variability, the relatively large redshift encouraged the authors to perform a search for an energy-dependent time delay. Observations from the second day were used for that purpose. Unlike the flare from PKS 2155-304 (see Section 2.6) the signal to background ratio in this case was only 2. Due to this high background contamination, a PDF for the background had to be introduced into the likelihood function for the first time. Events from the energy range 300-789 GeV, the upper edge corresponding to the last significant bin, were used. The sample was separated into a lower energy bin used to create the light curve template, and a higher bin used for the ML calculation. The delimiter between these two bins was set at 400 GeV, approximately corresponding to the median of the sample. The results, E QG,1 > 4.1 × 10 17 GeV, E QG,2 > 2.1 × 10 10 GeV for the subluminal scenario, and E QG,1 > 2.8 × 10 17 GeV, E QG,1 > 1.7 × 10 10 GeV for the superluminal scenario, did not further constrain the LIV energy scale, but confirmed the already existing limits on the quadratic term. The bounds on the linear term were an order of magnitude below the ones set by H.E.S.S. on PKS 2155-304 flare [67]. The authors did not discuss the reasons for the lower sensitivity, however, referring again to our rule of thumb (Table 1), it seems safe to conclude that this study benefited from the high redshift of the source, while paying dues to the lower gamma ray energies detected, the modest sample size, and a marginal flux variability.

Exploring Lower Time Variability with the Crab Pulsar Observations by VERITAS
The idea of using pulsar emission to search for LIV was first applied to Crab pulsar observations by EGRET [69]. The Crab pulsar (PSR J0534+2200) is located at the center of the Crab nebula at 2.0 ± 0.5 kpc [70] from Earth, and has a period of rotation of ∼33 ms [71]. In 2011, the Very Energetic Radiation Imaging Telescope Array System (VERITAS) reported the observation of gamma-ray emission from the Crab pulsar above 100 GeV [72]. Its phaseogram, i.e., its emission as a function of the pulsar rotational phase φ, distinctly shows a main pulse (referred to as P1) and an inter-pulse (referred to as P2) at a phase φ∼0.4 from P1. For the LIV analysis [73,74], the authors made use of the peak comparison (PC) method. This method can be used to look for an average phase delay ∆φ between photons from two different energy bands with mean energies E 1 and E 2 for the lower and the higher energy band, respectively: where d Crab is the distance to the Crab pulsar, P Crab its period and c the Lorentz invariant in vaccuo speed of light. Note that the phase φ is a practical quantity when describing pulsar behavior, nevertheless since ∆t = ∆φP Crab one immediately recovers Equation (3) from Equation (13) under the assumption that D n (z s ) ≈ d which is true for such nearby sources as pulsars. The authors used this method to compare the mean fitted pulse position obtained with VERITAS above 120 GeV to the one obtained with Fermi-LAT above 100 MeV [75]. The peak positions agreed within statistical uncertainties, therefore a 95% confidence upper limit on their timing difference of 100 µs could be derived. This limit was then converted into limits on E QG n by reversing Equation (13): yielding E QG,1 > 3.0 × 10 17 GeV and E QG,2 > 7.0 × 10 9 GeV in the subluminal scenario (S = −1). Note that in Equations (13) and (14), the distance parameter D n (z s ) from Equation (  4) was replaced by a more standard distance d as pulsars are sources within the Milky Way, hence their distance is not properly described by the redshift. This means that the last column of Table 1 is different in the case of pulsars, indeed E QG,n will be proportional to d 1/n . A variation of the DisCan method was also used in this work [74]. It was first introduced in 2008 [40] and, as its name suggests, consists in looking for the LIV parameter that best cancels out any time dispersion in the data. As such this method is a variation of the ECF with a different cost function. The variation consisted in the use of the Z 2 m test [76] as a test statistic (with m = 20 resulting from Monte Carlo optimization for this particular case) applied to the phased data to look for the potential LIV effect. This Z 2 20 DisCan method yields a best value of η 1 = −0.49 µs/GeV. The calibration of the method using 1000 Monte Carlo simulations allowed the authors to establish that this value was only 1.4 σ away from the null hypothesis and therefore compatible with it. The 95% confidence level limits on η 1 reached −1.2 µs/GeV and 1.1 µs/GeV for the lower and upper limits, respectively. These results were then translated into the following limits: E QG,1 > 1.9 × 10 17 GeV and E QG,1 > 1.7 × 10 17 GeV for the subluminal and the superluminal scenario, respectively.

Applying the Maximum Likelihood Method to the Crab Pulsar with MAGIC
The Crab pulsar was also observed and detected by MAGIC. The LIV analysis performed on the Crab pulsar [77] focused on the events from the P2 pulse as they reach higher energies, which increases the sensitivity to a LIV effect. For this analysis, the authors used ∼326 h of excellent quality data. This dataset was analysed with two different methods. Three energy bands (mean energies ∼75 GeV, ∼465 GeV, and ∼770 GeV,) were defined for the analysis but the analysis focused primarily on the two highest. The reason for this choice is that the emission's mechanism is likely to be different between the lowest and the highest energies of the pulse. Therefore the comparison focused on the two high energy bands, which are more likely to arise from the same mechanism that will not affect the search for a LIV effect. The first method used was the PC, already introduced in Section 2.8, which yielded the following limits on the LIV energy scale: E QG,1 > 1.1 × 10 17 GeV and E QG,2 > 1.4 × 10 10 GeV for the subluminal scenario, and E QG,1 > 1.1 × 10 17 GeV and E QG,2 > 1.5 × 10 10 GeV for the superluminal scenario. The second method used by the authors was the ML method, here used for the first time to analyse data from a pulsar. The likelihood approach follows what we introduced in Section 2.3, adapted to the study of events describe by their phase φ instead of their absolute time t. In addition, the likelihood included terms to describe nuisance parameters among which the parameters used to fit the pulse profile and the background events. The former was used to evaluate systematic uncertainties in the analysis, while the latter was particularly important in the case of pulsar located in a Nebula, itself an important and steady source of gamma rays. An extended investigation of the possible origin of systematic uncertainties in this work was performed, including the uncertainty on the absolute energy and flux scale, the possible contribution from events outside the pulse region, and the relatively large uncertainty on the estimation of the distance to the Crab pulsar. In total, the authors estimated the systematic uncertainties on E QG,1 to be less than 42 % and on E QG,2 to be less than 36 %. The obtained limits reached E QG,1 > 5.5 × 10 17 GeV and E QG,2 > 5.9 × 10 10 GeV, including systematic uncertainties, in the subluminal scenario, and E QG,1 > 4.5 × 10 17 GeV and E QG,2 > 5.3 × 10 10 GeV, including systematic uncertainties, in the superluminal scenario. The ML method, thus, provided limits a factor 4-5 more stringent than the limits obtained with the PC method.

Lorentz Invariance Violation Study on a New Vela Pulsar
The Vela pulsar (PSR J0835-4510), located at 0.29 +0.08 −0.05 kpc [78] from Earth, was observed by H.E.S.S. from March 2013 to April 2015 [53]. In order to reach an energy threshold as low as possible, the analysis only used events recorded by the large 28 m telescope telescope at the centre of the array. The LIV analysis [79] made use of 24 h of good quality data from 2013 to 2014. In this period, the telescope recorded about 10,000 pulsed events above ∼20 GeV. The energy range considered for the analysis was 20 to 100 GeV, yielding a statistics of ∼9300 excess events associated to the pulsar for a signal to noise ratio of ∼0.025. The authors used the same ML method as described in Section 2.9. The ON phase region was defined as the interval [0.5,0.6]. The signal template was obtained from the fitting of the low energy (20-45 GeV) events from the ON phase region by an asymetrical Lorentzian function (for the signal) plus a constant (for the background). This constant is determined from the fitting of events from the OFF phase region chosen as [0.7,1]. The authors used dedicated toy Monte Carlo simulations to calibrate their analysis by simulating mock data reproducing Vela's sample characteristics and injecting different simulated phase delays, similar to what was presented in Section 2.4. The method exhibits an almost unbiased reconstruction of the LIV induced delay. Therefore the results of the distribution of the reconstructed delay, when no LIV effect was injected, was used to evaluate the statistical uncertainty of the measurements as well as the systematic uncertainty. Applied to the 45-100 GeV range, the ML analysis provided a measurement of the delay φ = (−2.0 ± 5.0 stat ± 3.0 sys ) × 10 −2 TeV −1 compatible with no delay. The results were, therefore, converted to 95% confidence level lower limits on the linear term E QG,1 , yielding E QG,1 > 4.0 × 10 15 GeV and E QG,1 > 3.7 × 10 15 GeV in the subluminal and superluminal cases, respectively.

First Parallel Study of Energy-Dependent Photon Group Velocity and Gamma-ray Absorption on the Same Data Sample
As discussed in Section 2.2, Mrk 501 was already observed and studied in the search of LIV after a flare detected by MAGIC in 2005. In 2014, another flare was detected during a monitoring campaign of First G-APD Cherenkov Telescope (FACT) [80]. The alert of this flare triggered observations by the full array of five telescopes of H.E.S.S. on the night of 23-24 June 2014. Observations were performed at high zenith angle (63 • to 65 • ) leading to a high energy threshold of 1 TeV. The LIV analysis on this flare [81] was done using the ML method presented in Section 2.3. The only noticeable difference was the use of the variable η n , defined in Equation (6), as the main likelihood parameter and the explicit mention of the normalization factor depending on η n . The sample of events was divided between the 733 events between 1.3 TeV and 3.25 TeV, which were used to compute the template and the 662 events above 3.25 TeV, which were used to compute the likelihood and the best values of η n . It is important to note that in this specific analysis, given the high energy threshold of the observations, the low energy template included a potential LIV effect. In practice, while the delay in the template is usually taken as null (D = η n E n ), here they modelled it as D = η n (E n − E T n ) where E T is the mean energy of the events in the energy range of the template. As the low energy events were used to built the template, only the 662 high energy events were used in the likelihood in the search for a LIV effect in this dataset. The best fitted value of η n were compatible with the Lorentz invariant scenario. 1000 Monte Carlo simulations with no LIV effect were used to derive calibrated intervals from which uncertainties were derived. Finally, limits on the energy scale of LIV were set to E QG,1 > 3.6 × 10 17 GeV (E QG,1 > 2.6 × 10 17 GeV) and E QG,2 > 8.5 × 10 10 GeV (E QG,2 > 7.3 × 10 10 GeV) for the subluminal (superluminal) scenario. These limits include systematic uncertainties, the main one being the determination of the template. Note that, to date, this dataset is the only one that has been used to perform a time of flight study, as described in this section, and a universe transparency study later described in Section 3.1.3.

First Lorentz Invariance Violation Study on a Gamma-ray Burst Observed with Imaging Atmospheric Cherenkov Telescopes
More that two decades after the proposal by Amelino-Camelia et al., an opportunity presented itself to test the LIV on a signal from a GRB observed with IACTs. The MAGIC Collaboration announced a discovery of a GRB with IACTs for the first time ever [82] 16 . A signal from GRB 190114C was detected at energies above 1 TeV [59]. The analysis of this signal for the purpose of testing LIV started immediately. The ML method was applied (in fact, Equations (10) and (11) were adopted from the LIV study on GRB 190114C [83]). The most troublesome issue about the analysis was the formulation of the light curve template. The MAGIC observations started 62 s after the burst, almost completely missing the prompt phase, and detecting gamma rays from almost only afterglow phase of the GRB. The signal in the TeV energy band was observable until ∼40 min after the burst [84], however it was estimated that only the first 20 min (the duration of a single observation run) would be relevant for the test of LIV. After that, the signal rate became comparable to the background rate, meaning that it would not have considerably improved the sensitivity of the analysis, while at the same time, the systematic effects would have increased. The MAGIC data analysis revealed that during the first 20 min of observation about 700 gamma rays were detected with the energies in range of 0.3− 2 TeV. The intrinsic spectral distribution of events was well fitted with a power law [84]. More interestingly, the light curve also demonstrated a monotonic, power law decay of the flux. A monotonic change of flux is no more useful in searches for a spectral dispersion than no change of flux at all would be. A spectral dispersion, applied to a monotonic temporal distribution, would change the rate of a change, but not the functional shape of the distribution. Thus, any effect of a spectral dispersion would be undetectable. Therefore, in order to perform the LIV test, the authors used the light curve model obtained from theoretical inference, and based on the observations performed with the MAGIC telescopes and other facilities observing in lower energy ranges [84]. This template was dubbed theoretical by the authors of the LIV study. All ∼700 events were used to calculate the likelihood. Before estimating the values and the confidence interval of LIV parameters, the sensitivity of the method was estimated. This was done by creating 1000 mock data sets and using them to calculate the likelihood. Each mock data set was produced from the original data set by shuffling the arrival times of detected events and then randomly selecting events from this reshuffled data set. In this way, the generated mock data sets consisted of the same number of events as the original data set, with the same energy and temporal distributions. Furthermore, bootstrapping, the procedure of randomly selecting events from the existing (reshuffled) data set, allowed both the energy and temporal distributions to vary in line with their statistical uncertainties. In this way, these uncertainties were propagated to final result. Reshuffling, on the other hand, had the role of removing any correlation between the energy and arrival time, if present in the first place. Therefore, if there was any energy-dependent time delay present in the original data set, it would have been washed out by the reshuffling. After calculating the likelihood for each of the mock data sets, a distribution of the results was made, revealing a bias in the method, for which the final data, obtained on the real data set, were corrected. The same mock data sets were used to calibrate the confidence interval, as described in Section 2.3.
As was already mentioned, the light curve model adopted from [84] was constructed based on observations in lower energies, and theoretical considerations. The power law decay, observed with the MAGIC telescopes, in the model is preceded by a rather sharp peak. The peak was before the MAGIC observation window, so neither confirmed, nor disproved. So even before obtaining the final results of the LIV test, there was a genuine concern that such fast change of the flux was introducing artificially high sensitivity to the LIV effects. As a sort of a sanity check, the LIV analysis was performed on another light curve template. This template, dubbed "minimal", was a step function, with zero value before the burst, and constant value afterwards. Translated to the signal PDF (Equation (11)), it means that there is zero probability of a gamma ray being emitted before the burst, and equal probability of emitting any gamma ray at any time after the burst. This very simple function is clearly not the correct description of the intrinsic light curve. Nevertheless, it avoids sharp peaks not confirmed by observations, consequently, in a sense, minimizing the influence of the light curve template on the sensitivity to the LIV effects. This light curve template will cause the likelihood profile to be minimal and flat for small and negative values of the LIV parameters, thus preventing the estimation of the bias, and only allowing setting constraints on the subluminal scenario. The results obtained using this minimal model (E QG,1 > 2.8 × 10 18 GeV and E QG,2 > 7.3 × 10 10 GeV) are compatible with the ones obtained using the theoretical model, meaning that the usage of the theoretical model did not introduce unreasonably high sensitivity into the analysis.
The bounds on the LIV energy scale obtained in this study, were comparable to the most constraining lower limits present at that time. However, more than confirming the constraints resulting from other studies, the importance of this work was particularly in the fact that it was the first one ever performed on a signal from a GRB observed with IACTs. Especially in the upcoming era of the Cherenkov Telescope Array (CTA) 17 , which carries a promise of observing a few GRBs each year with significantly larger data samples for every GRB [86] 18 , the test of LIV on GRB 190114C presents an important stepping stone for the future of LIV research.

Lorentz Invariance Violation on Fermi-LAT Gamma-ray Bursts
In previous sections we laid out analysis methods and results of different studies performed on the IACTs data, searching for the signatures of energy dependence in the photon velocity. Results of all these studies are usually compared to the results from a benchmark work by Vasileiou et al. [38], where the authors collected four GRBs observed with the Fermi-LAT instrument 19 . Vasileiou et al. analysed the Fermi-LAT data from four bright GRBs with well determined redshifts: GRB 080916C (z = 4.35 ± 0.15 [87]), GRB 090510 (z = 0.903 ± 0.003 [88]), GRB 090902B (z = 1.822 [89]), and GRB 090926A (z = 2.1071 ± 0.0001 [90]). All of these are much farther away than any source used for LIV tests with IACTs. In addition, and unlike the case of GRB 190114C observed with the MAGIC telescopes, a quickly variable prompt GRB phases were observed in all these four cases. An LIV test was performed on each of these sources individually, and three different analysis methods were used on each source.
The PairView (PV) method was developed for the purposes of this study. It calculates once the energy-dependent differences in the arrival times for each pair (i, j) of photons in the sample: The distribution of l i,j will be peaked at η n defined in Equation (6), giving the value of the LIV parameter.
Sharpness maximisation method (SMM) [38,39] is analogue to the ECF method, which was previously applied to the MAGIC sample of Mrk 501 and explained in Section 2.2. It employs the aforementioned fact that an application of a spectral dispersion to a data set will decreases sharpness of the light curve. While the ECF method maximizes the power in the selected time interval, the SMM measures the sharpness of the light curve, e.g., 17 CTA ( [85], https://www.cta-observatory.org/, accessed 15 July 2021) is an array of Cherenkov telescopes currently being built in two locations. The approved "Alpha Configuration" in the Southern Site in Paranal Observatory (Chile) will consist of 14 Medium-Sized Telescopes and 37 Small-Sized Telescopes, covering the area of ∼3 km 2 . The Northern Site will be located in the Roque de los Muchachos Observatory (Spain), consisting of four Large-Sized Telescopes and nine Medium-Sized Telescopes, which will cover the area of ∼0.25 km 2 . 18 Note that this estimate was obtained for a larger number of telescopes in each site compared to the Alpha Configuration. 19 Considering our focus on research performed the Vasileiou et al. work is strictly speaking out of the scope of this review. However, it derived some interesting results, and other LIV studies are often compared to it, so we will outline its main points.
after applying an opposite dispersion as described in Equation (5). ρ is a fixed parameter making sure that events which are very close together are not considered in the denominator, because that would dominate the function. The intrinsic light curve is expected to be the sharpest one. Therefore, η n for which the light curve is the sharpest, will be the measure of the spectral dispersion present in the data sample. The third and final method used was the ML, which we already described in Section 2.3. Final limits on the LIV energy scale were obtained for each source individually, by taking average value of results obtained from three different methods and after accounting for systematic effects (see Table 5 in [38]). The most constraining lower limits resulted from the GRB 090510: E QG,1 > 2.2 × 10 19 GeV, E QG,2 > 4.0 × 10 10 GeV for the subluminal, and E QG,1 > 3.9 × 10 19 GeV, E QG,1 > 3.0 × 10 10 GeV for the superluminal scenarios. The staggering lower limits on the linear term, surpassing the Planck energy are the reason why every other LIV study is compared to this one. Interestingly, GRB 090510 was the one with the smallest redshift in the sample, and the only one with z < 1. However, it is also the only one in the sample which was classified as a short GRB, while the other three were long GRBs, with emission spread over somewhat longer time. Moreover, the highest energies in all four data sets were detected from GRB 090510. While the authors at first considered combining the results from all four sources into one single bound on the LIV energy scale, they gave up on the idea because the result from GRB 090510 was so much more constraining than the other three GRBs that a combination would not significantly increase the lower limit.
It should be noted that the Fermi-LAT detector is sensitive in the energy range of 20 MeV-300 GeV. The higher energy part of this band partially overlaps with the IACTs sensitivity range (∼30 GeV-few tens of TeV). However, in the overlapping energy region, the Fermi-LAT sensitivity deteriorates with increasing energy, while the opposite is true for IACTs 20 . Lower energy gamma rays, detectable with Fermi-LAT, are not absorbed by the EBL, thus enabling Fermi-LAT to detect sources at significantly higher redshifts than IACTs, increasing the sensitivity to LIV effects. On the other hand, Fermi-LAT reaches significantly lower energies than IACTs, which limits its sensitivity to LIV effects. These characteristics will be important when we compare different results in Section 4.

Modified Photon Interactions
Very soon after the modified photon dispersion relation was introduced, it has been realised that it can have consequences on kinematics and dynamics of the processes (see e.g., [92][93][94]). In some quantum electrodynamics processes, modifications of dispersion relation may cause the change of the reaction energy threshold. On the other hand, some processes forbidden by energy-momentum conservation law in Lorentz invariant scenario, may become allowed if the Lorentz symmetry is broken. In this chapter, we will look more closely into several of these phenomena.

Testing Lorentz Invariance Violation with Universe Transparency
The universe is filled with low energy photon fields such as extragalactic background light (EBL), cosmic microwave background (CMB) and radio background (RB). Gamma rays traversing cosmological distances scatter off those photons creating electron-positron pairs. Consequently, their flux, observed from Earth, is attenuated [95][96][97][98]. The EBL is responsible for the attenuation of gamma rays in 10-10 5 GeV range, which roughly corresponds to observable energy range of current IACTs. Unfortunately, direct EBL measurements are obstructed by bright foreground emissions, mainly zodiacal light [99], which makes it hard to determine its precise spectrum (for more information about the EBL, photon-photon interactions and the opacity of the universe to gamma rays we referee the reader to [100] and references therein). To tackle this problem, different phenomenological approaches predicting overall EBL spectrum have been followed. Remarkably, EBL models obtained through different methodologies, such as Franceschini et al. [101], Domínguez et al. [102] and Gilmore et al. [103], are in a good agreement. These models were tested on VHE data from sets of AGN by current IACTs [104][105][106]. Those tests were done presuming Lorentz invariance.
The gamma-ray spectrum observed from Earth is usually written as a convolution of the source intrinsic spectrum and the EBL attenuation effect: where E is the observed gamma-ray energy and z s is the redshift of the observed source. τ(E, z s ) is the optical depth, dependent on the two aforementioned parameters and is given by 21 : In this expression, • denotes the energy of an EBL photon in the comoving frame, while n( , z) is the comoving number density of EBL photons per unit energy. • The probability of the interaction between a gamma ray and background photons is given by the cross section σ γγ (s), where s is the square of the center of mass energy. In the gamma-ray energy range relevant for IACTs, by far the most dominant channel is the Breit-Wheeler process of electron-positron pair creation [107]. • The angle of interaction between a gamma ray and EBL photons is indicated by θ . • th denotes the EBL energy reaction threshold for electron-positron pair creation, i.e., the minimal energy of an EBL photon, in the comoving frame, necessary for the reaction to take place. Derived from the kinematics laws of special relativity, it can be expressed as: The threshold energy, and its changes due to modifications of the special relativity kinematics, will play a vital role in constraining E QG . • The final integral accounts for the distance traveled by the gamma ray, assuming flat ΛCDM cosmology: Beyond the gamma-ray horizon (τ(E, z s ) = 1) the universe becomes progressively opaque for VHE gamma rays (for further readings on this topic in connection with IACTs we suggest [108] and references therein). For a sources at redshift 0.034, which is a redshift of Mrk 501, gamma-ray horizon is around 10 TeV [101]. When doing calculations, one must be careful to take into account cosmic expansion and notice that measurements are affected by a factor (1 + z). Namely, E and change along the line of sight inversely proportional to (1 + z); for example, = /(1 + z). 21 Henceforward, when denoted with prime physical quantity is written in the comoving frame at which the interaction occurs. When prime does not occur, it is written in the observer's frame of reference.

Influence of Lorentz Invariance Violation on Universe Transparency
Detected gamma-ray emission up to ∼22 TeV from Mrk 501 [109] in 1997 by High Energy Gamma Ray Astronomy (HEGRA) experiment 22 hinted that the universe is more transparent to VHE gamma rays than expected. One possible solution to this newly arisen problem was the aforementioned modification of photon dispersion relation. Added terms in the photon dispersion relation can cause a change in the energy threshold for pair creation, consequently leading to changes in the gamma-ray absorption. In this scenario, the new energy reaction threshold is [111]: Changes in energy reaction threshold are depicted in Figure 1 for a head-on collision and z = 0. As defined in Section 1.2, S = +1 for superluminal, and S = −1 for subluminal behaviour. This modified energy reaction threshold has been derived under two assumptions: (i) the standard energy-momentum conservation law is maintained in LIV scenario, and (ii) LIV affects only dispersion relation for photons, while electrons remain unaffected. Indeed, the effects of LIV on electrons were strongly constrained by independent studies 23 . Therefore, the majority of studies considering electromagnetic interaction rely on the assumption that only the photon dispersion relation is modified by LIV.
Modifications in the energy reaction threshold could lead to changes in the observed spectra of a distant source, depending on the LIV scale [92,97,[113][114][115][116]. In the superluminal behaviour (S = +1), modifications in the photon dispersion relation will cause lowering of the energy reaction threshold for the electron-positron pair creation. In that case, gamma rays would be absorbed by lower energy photon fields than in the Lorentz invariant case. For example, a 50 TeV gamma ray in Lorentz invariant scenario does not have enough energy to reach the reaction threshold with CMB photons. However, in a LIV superluminal scenario, for sufficiently low values of E QG , the reaction threshold will be reached (see Figure 1). This would lead to additional depletion of the most energetic photons resulting in a steeper observed spectrum. Still, so far no way was found to unambiguously disentangle the effects caused by the lowering of the photon-photon energy threshold due to LIV, from the effects arising due to Lorentz invariant EBL attenuation, or intrinsic properties of the source such as a spectral cut off. This is arguably the main reason why all experimentally set limits on E QG using the universe transparency to gamma rays were derived for the subluminal behaviour only. In the subluminal scenario (S = −1), modifications of the photon dispersion relation will lead to an increase of the energy reaction threshold, resulting in a reduced opacity of the universe to VHE gamma rays. Moreover, the reaction threshold as a function of the gamma-ray energy will have a global minimum [115] as can be seen in Figure 1. Note that there is no equivalent minimum in the Lorentz invariant nor LIV superluminal scenario, since th , as defined in Equation (19), is a monotonous function of the gamma ray energy. Contrary to the Lorentz invariant scenario and the LIV superluminal scenario, for which once the reaction energy threshold is reached the reaction is allowed for all gamma rays with energies above the reaction energy threshold, a pair creation in the LIV subluminal scenario is kinematically forbidden for gamma-ray energies higher than the reaction energy threshold. The existence of the global minimum implies that the energy domain of EBL photons, as targets for absorption of VHE gamma rays, would be reduced, regardless of the gamma-ray energy. Consequently, a certain number of gamma rays would evade absorption and thus reach the Earth. This most particularly holds for gamma rays with energies above the position of the energy threshold minimum.
The aforementioned Breit-Wheeler cross section, as a function of the gamma-ray energy, is shown in Figure 2. In the Lorentz invariant scenario, it is represented with a black dashed line. Once the reaction energy threshold is reached, the cross section rises quickly. At the gamma-ray energy roughly twice the threshold energy [117], the Lorentz invariant cross section reaches its maximal value of 1.70 × 10 −25 cm 2 [118]. Afterwards, as the gammaray energy increases, the cross section drops and asymptotically approaches zero. In the superluminal scenario, the energy reaction threshold is lower than in the Lorentz invariant scenario. Moreover, lower E QG,1 results in lower reaction threshold. The cross section shape remains the same as in the Lorentz invariant scenario, although, it becomes narrower as E QG,1 decreases and reaches its maximum at lower gamma-ray energies. A somewhat more interesting development of the Breit-Wheeler cross section occurs in a subluminal LIV scenario. There are three distinct cases: (i) As we saw in Figure 1, for E QG,1 low enough, the reaction energy threshold will never be reached. Consequently, the cross section will be zero for all gamma-ray energies (red full line in the bottom panel of Figure 2). (ii) For higher values of E QG,1 , the horizontal line will be crossed twice. Hence, there will be a lower and an upper reaction energy thresholds, and the reaction will be possible for gamma-ray energies between these thresholds. For relatively low E QG,1 , this interval will be narrow, and the cross section will never reach its maximum possible value of 1.70 × 10 −25 cm 2 (green and violet full lines in the bottom panel of Figure 2). (iii) For even higher values of E QG,1 , the gamma-ray energy interval between the reaction energy thresholds will be wide enough for the cross section to reach its maximum possible value. Moreover, the cross section will start to decrease with increasing gamma-ray energy, roughly following the shape of the Lorentz invariant cross section. However, as the gamma-ray energy rises, the cross section reaches a local minimum, and starts increasing to reach its maximum possible value once again, just below the upper reaction threshold. Once the threshold is reached, the cross section is cut off (light and dark blue full lines in the bottom panel of Figure 2). If one continued to increase the value of E QG,1 , the second peak would become sharper and move towards higher energies, and the intermediate part of the cross section would more closely follow the Lorentz invariant cross section. In addition to these three cases, there are borderline cases. Between cases (i) and (ii), for E QG,1 precisely such that the horizontal line in Figure 1 is a tangent to the energy threshold line, the reaction energy threshold will be reached at precisely one gamma-ray energy, and the cross section will be a vertical line at that energy, and zero elsewhere. Between cases (ii) and (iii), the cross section would reach its maximum possible value, and monotonically drop to zero.
After seeing how modifications of photon dispersion relation influence the reaction energy threshold and the cross section, now it is time to see how the optical depth changes due to those modifications, since all experimental limits on E QG , based on the universe transparency, were set only for the subluminal scenario, we will focus only on this scenario. In Figure 3 we depicted a hypothetical gamma-ray absorption for a source at redshift z s = 0.03 and gamma-ray energies up to 100 TeV 24 , assuming different values of E QG,1 . In the top panel of Figure 3 the gamma-ray horizon is denoted with a maroon line. As previously mentioned, beyond the gamma-ray horizon the universe becomes increasingly opaque for VHE gamma rays and thus the probability of their detection is diminishing. In the Lorentz invariant scenario, once the gamma-ray horizon is reached, the optical depth only increases. On the other hand, in the subluminal LIV scenario, the optical depth has a global maximum, different for different energies depending on the value of E QG , after which it decreases again. At some point it goes below the gamma-ray horizon allowing the gamma rays to evade absorption, which would lead to the recovery of the photon flux. The higher the E QG scale, the higher the energy of the gamma-ray at which the recovery would occur. The absorption coefficient (e −τ ), as a function of a gamma-ray energy (E) is depicted in the bottom panel of Figure 3 and shows how the survival probability of the photons behaves in this scenario.
It should be noted that there are other phenomena other than LIV which could leave imprints in the spectra of observed sources. Most notable are the axion-like particles, into which VHE gamma rays can oscillate in the presence of the external magnetic field. Nevertheless, imprints which axion-like particles and LIV would potentially leave could be mutually distinguished. Namely, axion-like particle imprints should be independent of the source redshift, while at the same time dependent on the magnetic field. The opposite holds for the LIV. For now, these two effects are being investigated separately, even in studies investigating both of them (see, e.g., [97]). For more information about axion-like particles and their searches with IACTs, we refer the interested reader to [119,120].

Testing Lorentz Invariance Violation on Universe Transparency
The first experimental test of LIV on the EBL absorption of gamma rays using data from IACTs was performed by Biteau & Williams [97]. The authors derived a simplified expression for optical depth using analytical methods. Furthermore, they constructed the EBL spectrum using 86 published gamma-ray spectra of 30 blazars with well-established redshifts. A total of ∼270,000 gamma rays constituted this gamma-ray sample. The EBL model was described through eight free parameters, denoted with A. In the LIV scenario, E QG was added as one additional free parameter to the analysis, using the optical depth dependency on LIV     described in the previous section. Furthermore, the EBL parameters were allowed to vary with E QG . In order to investigate the possible effects of LIV, Biteau & Williams compared spectra of the aforementioned blazars under the assumption of Lorentz invariance on the one hand, and under the assumption of LIV on the other hand. The effect of LIV was quantified using a test statistic (TS) defined as follows, with L = exp(−TS/2): The term χ 2 E QG , A QG represented the best fit in the LIV scenario. The χ 2 (∞, A ∞ ) represented the best fit in the Lorentz invariant scenario since, as mentioned in Section 1.2, letting E QG,1 → ∞ leads to the Lorentz invariant photon dispersion relation.
Biteau & Williams adopted the formalism of [115], which presumes that LIV affects both photons and leptons equally. Under that assumption, the last term in the energy threshold expression (Equation (21)) gains another factor of (1 − 2 −n ). One of the 86 spectra selected for this study was the spectrum of Mrk 501 historical flare detected by HEGRA in 1997 [109]. When performing the analysis on the originally published Mrk 501 spectrum, Biteau & Williams showed that E QG,1 was approximately E Pl at 4σ level. However, when a newly derived spectrum of the same data set (from [121]) was used, the significance decreased to 2.4 σ. The re-analysed spectrum had better energy resolution, therefore it excluded the initially obtained highest energy data point. This example demonstrates how a single spectrum, and the most energetic photons in it, can greatly influence the final result. The first experimentally set 95% confidence level lower limit, using universe transparency to VHE gamma rays, was found to be E QG,1 > 9.5 × 10 18 GeV. This value changed to E QG,1 > 8.6 × 10 18 GeV when a 10% systematic uncertainty on the energy scale (typical for IACTs) was accounted for.

The Most Constraining Limits Based on Single Source Analysis
Seventeen years after the historical flare from Mrk 501, which triggered discussions about the influence of LIV to universe's transparency to VHE gamma rays, another bright flare brought Mrk 501 in the spotlight once again [122]. This time it was used to constrain E QG by the H.E.S.S. Collaboration using two independent channels [81]. The test of energy-dependent photon group velocity was described in Section 2.11, while the spectral analysis will be discussed here. Only the possible subluminal behaviour for linear or quadratic contributions was investigated. The pair production cross-section was calculated according to [123]. In this approach, the modified expression for the square of the center of mass energy s can be written as: The dependence of the cross section on s is considered to be the same in LIV and Lorentz invariant scenarios. The cross section as a function of a gamma-ray energy for the background photon of 11 meV and a head-on collision is depicted in Figure 2.
To quantify the possible effects of LIV on a spectra of this bright flare, authors defined a TS similar to Equation (22): However, unlike the TS defined in Equation (22), this TS did not contain free-varying EBL parameters. The EBL model of Franceschini et al. [101] was used and the optical depth was calculated in a standard way, as described in Equation (18). The χ 2 values in Equation (24) were obtained by varying E QG logarithmically and fitting the measured spectrum of a flare by an assumed intrinsic spectrum convoluted with the EBL attenuation effect. The intrinsic spectrum was assumed to be a simple power law. From the TS profiles, lower limits on E QG at 95% confidence level were set to E QG,1 > 2.6 × 10 19 GeV and E QG,2 > 7.8 × 10 11 GeV for the linear and quadratic contributions, respectively.
A peculiarity of this work lies in the fact that the same data set was used to put the constraints on E QG via energy-dependent time delay and the universe transparency effects, which we will discuss in more detail in Section 4.

On How the Most Constraining Limits Were Obtained
At the moment of writing this review, the best lower limits on E QG obtained by testing the universe transparency were obtained by Lang et al. [124]. The starting data sample consisted of 111 spectra from 38 different sources available in the online catalogue TeVCat 25 , of which only a subset was used for setting constraints on LIV. In order to select only the relevant sources for probing LIV, the authors defined the attenuation a(E, z) as the ratio between the measured J meas (E) and the intrinsic J int (E, z) spectra of each source: Lang et al. calculated the ratio between the attenuation assuming LIV and the attenuation in the Lorentz invariant scenario, at the maximum energy (E max ) measured in a given spectrum. Only 18 spectra from 6 different sources for which the a LIV /a LI ratio differed by at least 10% and which could be used to further constrain E QG were selected for further analysis.
In general, the intrinsic spectrum is obtained via the process of so-called deabsorption which consists in reverting the Equation (17). In previous studies, deabsorption was done under the assumption of Lorentz invariance, not taking into account LIV in this step of the analysis. In order to rectify it, Lang et al. used an energy interval which they call fiducial region. It is defined as the energy range starting at the lowest measured spectral point. The highest spectral point is the last one at which the difference between the fluxes assuming LIV and assuming Lorentz invariance are indistinguishable, considering the measurement uncertainties. Therefore only measured spectral points from bins that satisfy the following condition were used to determine the intrinsic spectrum: Throughout their work, the authors assumed ρ = 1, implying the tightest energy interval and hence leading to the most conservative limits on E QG . Every intrinsic spectrum (J int ) was modeled as a power law with an exponential cut off. For each selected spectrum (J int ), the energy spectrum on Earth (J cal ) was computed for multiple E QG,n values using Subsequently, all computed spectra were compared with the complete measured spectra J meas using a log-likelihood method. Finally, the authors combined the likelihood results from all the sources to achieve the best possible sensitivity.
In their work, Lang et al. report 2σ confidence level lower limits on E QG obtained with three different EBL models, using the same procedure. The most conservative limits were derived using the EBL model by Domínguez et al. [102] and were set to be E QG,1 > 6.85 × 10 19 GeV and E QG,2 > 1.56 × 10 12 GeV.

Constraints on Violation of Lorentz Invariance from Atmospheric Showers Initiated by Multi-TeV Photons
The imaging technique of Cherenkov telescopes relies on recording flashes of Cherenkov light produced in the atmosphere by ultrarelativistic particles constituting extensive air showers. When a VHE gamma ray enters the Earth's atmosphere, it is absorbed in the Coulomb field of an atomic nucleus in the air, creating an electron-positron pair. Each created particle carries approximately one half of the primary gamma ray's energy. Leptons are emitting additional gamma rays through bremsstrahlung, each of which again go through the process of pair production. In that way, an electromagnetic cascade is created. The pair-creation is, fundamentally, the same process as the gamma-gamma interaction in which VHE gamma rays get absorbed by the EBL. Therefore, if the gamma-gamma interaction was affected by modifying the photon dispersion relation, it would also influence the development of particle showers in the atmosphere. This interesting notion was proposed by Rubtsov, Satunin & Sibiryakov in [125] and tested on data from HEGRA and H.E.S.S. in [126].
A shower development is governed by the Bethe-Heitler (B-H) process. In particular, the depth of the first interaction in the atmosphere is exponentially distributed with the mean value inversely proportional to the cross section [127] Here Z is the atomic number of the nucleus, α is the fine structure constant, and m e is the electron mass. In LIV scenario, the cross section will not change significantly for superluminal photons, unless the threshold for photon decay is reached. However, in that case, the photon decay will be the dominant process, making the LIV influence on the B-H process negligible.
As a consequence, the shower development in the LIV scenario will be impeded. The first gamma-gamma interaction will occur deeper in the atmosphere, and the effect will be more pronounced for higher gamma-ray energies. This will lead to showers reaching their maximal sizes also deeper in the atmosphere. Height of the shower maximum is an important parameter in IACTs data analysis (see, e.g., [128]). Depending on the experimental setup and the details of the data analysis, changes in the B-H cross section might lead to the showers induced by the most energetic gamma rays being misrepresented and excluded from further analysis. Ultimately, this will result in an apparent cut off in the spectrum at the high end.
Rubtsov, Satunin & Sibiryakov applied this method on two independent measurements of the Crab nebula spectrum. The first one was obtained by the HEGRA Collaboration, based on 385 h of observations performed between 1997 and 2002 [129]. The highest energy bin in the spectrum was centered at 75 GeV. The analysis method used to compare the spectra in Lorentz invariant versus LIV scenarios, and to determine limits to the LIV energy scale was based on ML method similar to the one used in [81]. They set a limit to the LIV energy scale to E QG,2 > 2.1 × 10 11 GeV. The second sample was the Crab nebula spectrum measured by the H.E.S.S. Collaboration, based on 4.4 h during the flaring episode in March 2013 [130]. In this case, the spectrum was determine up to ∼40 TeV. Because of the smaller data set and lower energies reached, the result was less constraining: E QG,2 > 1.3 × 10 11 GeV.
Note that this effect would be opposite to the one caused by modified absorption of gamma rays on the EBL (described in Section 3.1). The constraints on the LIV energy scale set in the work by Rubtsov, Satunin & Sibiryakov were lower than the ones obtained based on the universe transparency to gamma rays. However, the constraints based on the universe transparency were obtained assuming that shower development, and consequently measurement with IACTs, is not modified by LIV. The work by Rubtsov, Satunin & Sibiryakov tested and validated that assumption.

Constraints on Lorentz Invariance Violation Based on Photon Stability
A modifying term in the photon dispersion relation can be treated as the mass term in the (unmodified) dispersion relation of a massive particle. Assigning a photon a mass in the superluminal scenario renders it unstable and prone to decay. A superluminal photon of energy E γ can: • decay into an electron-positron pair

This reaction becomes possible under condition [131]
where m e is the electron mass.
While the photon splitting has no reaction threshold, and is kinematically allowed for every superluminal photon, the process rate (see [132]) is significantly smaller than the photon decay rate [131,133].
Both processes have a similar effect on the observed spectra from astronomical sources, i.e., spectra will be attenuated at higher energies. Since there is no reaction threshold for the photon splitting, the attenuation will happen gradually. The photon decay rate quickly increases with the gamma-ray energy once the reaction threshold is reached. Consequently, it will be manifested as a cut off in the spectrum. The effects are similar to spectral attenuation due to gamma-ray absorption on the EBL. Therefore, in order to test for photon decay or photon splitting, one needs to exclude the possibility of EBL absorption. An obvious choice of a VHE gamma-ray source for these studies is the Crab nebula. Its spectrum reaches energies well above 100 TeV [134,135], while at the same time, because of its small distance from the Earth, EBL absorption is virtually negligible. These effects were independently tested for on Crab nebula spectral measurements in several studies (see, e.g., [131,132,136]), however, the results setting substantially stronger constraints came not from any IACT experiment, but from the HAWC Collaboration 26 [133]. The photon decay was used to constrain both liner and quadratic terms, obtaining E QG,1 > 2.2 × 10 22 GeV and E QG,2 > 0.8 × 10 14 GeV, respectively. The photon splitting was used only for the quadratic term, resulting in a much stronger constraint than the photon decay, E QG,2 > 1.0 × 10 15 GeV. Sheer moments before concluding this review, a very exciting result was published by LHAASO 27 , announcing a detection of gamma rays with energies up to 1.4 PeV from 12 sources [15]. Based on these measurements the LHAASO Collaboration performed a similar study as HAWC. They searched for a cut of in spectra of the two sources with the highest energies LHAASO J0534+2202 (Crab nebula) and LHAASO J2032+4102 (the source which the 1.4 PeV event was associated with) [138]. Due to significantly higher spectral measurements by LHAASO, the resulting constraints on the LIV energy scale were also higher. Specifically, their most constraining limits were based on the analysis of LHAASO J2032+4102 spectrum. After including the systematic uncertainties, the limits were set to E QG,1 > 1.2 × 10 24 GeV and E QG,2 > 1.1 × 10 15 GeV, when the photon decay was considered. As in the HAWC analysis, only quadratic term was constrained based on photon splitting, resulting in E QG,2 > 2.0 × 10 16 GeV.
The scope of this review are studies performed with IACTs. Nevertheless, we included these results from other type of observatories for comparison purposes. In addition, similar studies could be performed with IACTs. With the prospect of the CTA to be commissioned in the next few years, and a recent result from the MAGIC Collaboration measuring the Crab nebula spectrum up to 100 TeV [139], feasibility of a similar study with IACTs does not seem so far-fetched.

Summary and Discussion
In this section we will review and compare the results presented in previous sections. The results are summarised and listed chronologically in Table 2. A quick glance already reveals that the constraints, both on linear and quadratic terms, are the strongest in the case of the photon stability measurements [133] by LHAASO, surpassing the Planck energy by four orders of magnitude in the case of the linear term. Apparently, that is the effect the most sensitive to LIV. However, as already stated in Section 3.3, photon decay and photon splitting are processes only allowed in superluminal scenario, and cannot be used to constrain E QG in the subluminal scenario. Nevertheless, it was important to constrain this LIV effect. Possible photon decay and photon splitting competes with modified absorption of gamma rays on EBL. Without a confirmation of the photon stability to these substantially high values of E QG , it would be virtually impossible to resolve and independently search for LIV effects of modified universe transparency to gamma rays. As it happens, the second most constraining bounds were set precisely on measurements of the universe transparency for gamma rays, or gamma-gamma interaction. The most notable result was obtained by Lang et al. [124] through a simultaneous analysis of spectra of several AGNs. As already argued, combining different sources, should wash out dependence on the properties of a given source and, therefore, help to lift the degeneracy between the LIV and source intrinsic effects. The main characteristics of a desirable source for this method are large source distance and the highest spectral energy measurement. Combining various sources in a single study enables the most pronounced characteristic of each source to be fully exploited. However, the method heavily relies on the EBL modeling and the assumptions made on the intrinsic source spectra. Uncertainties of EBL models, as well as discrepancies between the different models, are the predominant source of systematic effects. Lang et al. considered this uncertainty by presenting the results of their analysis using three different EBL models. In this review, we presented the most conservative result. Table 2. Chronological census of bounds on the LIV energy scale as reported in the respective publications. The columns show from left to right: (i) LIV effect used to probe the QG, (ii) name of the observed source (when more than one source was used, we give the numner of the sources in the sample), (iii) type of the observed source, (iv) source redshift for extragalactic sources, and distance in kpc for sources in Milky Way, (v) analysis method used to perform the test, (vi) lower limit on the linear term in the modified photon dispersion relation, (vii) lower limit on the quadratic term in the modified photon dispersion relation, (viii) the experiment which produced the data sample, (ix) reference for the respective publication, (x) section in which the result was discussed in this work. Markers (+) and (−) in columns (vi) and (vii) represent superluminal and subluminal behaviours, respectively. The lower limits are expressed on the 95 % confidence level. Symbol † represents the best fit result in case where the estimated parameter was > 2σ away from Lorentz invariant scenario. A second look at the Table 2 clearly shows that analyses based on photon interactions result in stronger constraints on E QG compared to the ones based on photon time of flight. A study by the H.E.S.S. Collaboration was the only one so far in which both tests were performed on the same data set [81]. Granted, the two analyses were performed independently, the time of flight test ignoring modified gamma-ray absorption and vice versa. Nevertheless, it allowed a more direct comparison of the two effects of LIV. The constraints based on absorption of gamma rays on the EBL were more stringent by one order of magnitude on the quadratic term, and two orders of magnitude on the linear term, compared to the constraints based on energy-dependent time delay. Apart from the source distance and the detected energy, when it comes to time delay studies, another important property comes into play. Indeed, fast variability of flux is crucial in constraining emission times of individual photons (see Table 1). It should be noted that, while changes in flux do not strongly interfere with analysis based on EBL absorption, it is extremely important that the spectrum remains constant. A change in the spectrum would introduce additional uncertainties and hence lower the analysis sensitivity. Therefore, one may argue that a faster flux variability is needed to increase the sensitivity of time of flight tests. However, even the most sensitive of the time of flight analyses ( [38] for the linear, and [81] for the quadratic term) are below the sensitivities of analyses based on photon interactions, suggesting that (assuming the same E QG ) LIV affects interactions more strongly than it does the photon group velocity. So, is there a point in testing energy-dependent photon group velocity, when we were able to set much stronger bounds on E QG through tests of modified photon interactions? As we pointed out in Section 1.2, the theory of QG has not been formulated yet. Consequently, we do not know what the effects of QG are. It is quite possible that the photon group velocity is affected by the LIV, while interactions remain unaffected, or the other way around. Another possibility is that both interactions and propagation are affected, but on different scales, effectively introducing separate and different values of E QG . It should be remembered that the modified dispersion relation (Equation (1)) is merely a mathematical model facilitating experimental tests of LIV, and, in the most general case, different values of E QG are applicable in different cases. Yet another imaginable scenario is that both subluminal and superluminal behaviours occur as different effects of QG. In that case, they might start to manifest at different energy scales, resulting in an even more complex expression for the photon group velocity (Equation (2)).
Considering the time delay studies alone, we observe a gradual increase of the lower limits on E QG . Of course, a study declaring a stricter limit (or a detection for that matter) is more likely to be published, however we would like to argue that this improvement is a result of several circumstances: (i) improvement on the performances of detectors allowed observations of sources at larger redshifts with better sampling, (ii) longer operations increased the probability of observing transient phenomena as flaring AGNs and GRBs, as well as increased statistics on observations of pulsars, and (iii) analysis techniques, the ML in particular, have been refined over time, allowing for higher analysis sensitivity. A notable exception is the result obtained on the observation of GRB 090510 with Fermi-LAT. Published in 2013, the study by Vasileiou et al. [38] still holds the record for the most stringent bound on E QG,1 . Let us analyse where this sensitivity came from. As conveyed in Table 1, the sensitivity of time of flight analyses increases with the energy of the gamma rays within the sample (E max ) and the redshift of the source (z s ), and decreases with the timescale of flux variations (t var ). GRB 090510 has the second and third criterion well satisfied. The data set used in the analysis was taken from a time interval of less than 5 s (see Figure 1 in [38]), significantly shorter than any other sample covered in this review. At the same time, the redshift of ∼0.9 is more than two times larger than the second furthest source considered in Table 2. Apparently, the downside of the low highest energy in the sample (31 GeV) is more than well compensated by these two advantageous characteristics. In the case of the quadratic term, the relationship between E max , t var , and z s is somewhat different. E QG is still inversely proportional to the variability timescale, however t var now enters through a square root, decreasing its influence on the analysis sensitivity. Furthermore, E QG depends on the source redshift as z ∼2/3 s (in contrast to z ∼1 s in case of the linear term). Given that all sources are at redshifts smaller than 1, the influence of redshift is weaker on the quadratic term. Weaker influences of variability and redshift make more room for the influence of E max . Now, the tables have turned. The bound on E QG,1 set using MAGIC observation of GRB 190114C [83] was only about a factor 4 (7) below the limit from GRB 090510 for the subluminal (superluminal) behaviour. In the quadratic term, the influence of the highest gamma-ray energy is more pronounced than in the linear term. So, the constraint on E QG,2 became more stringent because of 2 TeV photons in the GRB 190114C sample. Actually, regarding the quadratic case, the strongest constraint came from the H.E.S.S. analysis of AGN Mrk 501, due to ∼20 TeV photons in the sample [81].
Emission time from pulsars is excellently constrained, with the variability timescales down to t var ∼10 ms. Furthermore, gamma rays from pulses are detected up to ∼7 TeV. Unfortunately, their relatively close proximity renders LIV analyses on pulsars comparatively less sensitive, in particular in the linear scenario. On the other hand, in the quadratic scenario, their fast variability gives them a huge advantage, making them competitive sources for the search of LIV effect. Nonetheless, contrary to flaring AGNs or GRBs, new observations of pulsar such as the Crab one do not depend on luck but can be carefully planned. This continuous accumulation of new data can lead to a predictable increase of statistics and therefore improved sensitivity of the analysis. This is particularly interesting as the current ML analyses are still mainly limited by background fluctuations and systematics such as the pulse shape and its energy evolution in the case of the Crab pulsar. In [77] the authors stated that a total dataset of ∼2000 h on the Crab pulsar is within reach for the MAGIC collaboration alone given the regular observations performed for calibration purposes. And this data set could be further enlarged by taking into account the data accumulated on the Crab pulsar by H.E.S.S. and VERITAS (see also Section 5.2) with the potential of addressing some of the main limitations of the current analyses mentioned above and thus exploring QG scales beyond the current best limits, in particular in the quadratic scenario. Lastly, only pulsars with periods of a few tens of milliseconds have been detected with IACTs so far. Detection of VHE gamma rays from pulsars with periods down to a few milliseconds would allow to constrain the emission time more strongly. This would lead to an improvement of the current limits on E QG by an order of magnitude, further exploiting the potential of these sources in the search for LIV.

An Eye on the Future
In previous sections, we discussed the evolution of the search for LIV effects with IACTs. The studies, which we reported on, have set strong constraints on the LIV energy scale, and significantly restricted the parameter space. More importantly, in those works, diverse ideas were proposed, various analysis methods had been developed, and different effects investigated. Ever since the first LIV study with an IACT, numerous ameliorations have been brought to the field in the form of technical or analysis improvements and the future ahead of us is certainly no exception. The most important question at this point is where do we go from here. What can we do to accelerate the research and contribute to the understanding of QG? In this section we try to outline some ideas on what could be the next steps in that regard. Hopefully, this will motivate the reader to perform some of the proposed research. We certainly intend to take our part in this endeavour.

Refinement of the Analysis Technique
In Section 2.3 we presented the state of the art ML method for testing energy-dependent photon group velocity. Here, we will discuss some possibilities for improving the method and increasing the analysis sensitivity.
Firstly, let us return to the definition of the likelihood function (Equation (10)) and in particular how individual events are selected and weighted. Note that the probability for each event to be a part of the signal, as defined in Equation (12), is the same for every event. The same is true for the probability for being a part of the background. Moreover, in order to reduce the background, the IACT data analysis involves applying cuts on the parameters describing events properties. This approach inevitably leads to cutting out some of the signal events as well. As we have seen, most of the LIV analysis methods rely on individual events. Therefore, cutting out signal events reduces the analysis sensitivity. However, in a recent paper, D'Amico et al. [140] proposed an alternative IACT data analysis method, which considers all events in the ON region without applying any cut. Instead, PDFs for the signal and the background are calculated based on the parameters of individual events. This method could, therefore, be used to calculate p as PDFs, without discarding any events, whether of signal or background origin.
Secondly, the real strength of the ML method relies in its modularity, meaning that the components of the likelihood function listed in Section 2.3 can be refined, and additional terms describing nuisance parameters can be added without limitations. Furthermore, likelihoods from individual targets can easily be combined in a joint likelihood (see Section 5.2). For example, in case the intrinsic spectrum changes with time, the function Φ int (E) can be generalised to Φ int (E, t). Additionally, by taking the product F(t + η n E n )Φ int (E) in Equation ( 11) we assumed that the emission time t of individual photons does not depend on their energy. Indeed, in the LIV studies performed so far, there was no strong evidence of changes of spectral shape during AGN flaring episodes nor in GRBs in the VHE gamma-ray range. Nevertheless, waving this simplification might increase the sensitivity of the ML analyses. Present day state of the art models used to describe emission from astrophysical sources are nowhere near refined and accurate enough to predict the exact emission time of each particular photon. Hopefully, future emission models will be precise enough to allow creating emission light curve templates simultaneously depending on emission time and energy. E.g., instead of having independent temporal and spectral distributions as F(t + η n E n ) Φ int (E), we could take a two-dimensional distribution as F 2D (t + η n E n , E) to account for potential source intrinsic delays. Therefore, further progress in the modelling of the sources emission mechanism such as initiated in [45] will certainly play an important role in future LIV studies.

Combining Data from Different Sources and Instruments
We have already seen on the example of the gamma-ray absorption (see Sections 3.1.2-3.1.4) how a combination of several sources in a single study improved the sensitivity of the analyses. An equivalent approach is still to be fully applied when testing the photon group velocity. Indeed, a combination of sources observed at a wide range of redshifts is the key to disentangling of intrinsic source effects from a real LIV effect. A source-intrinsic energydependent photon emission time can mimic an effect of LIV, leading to a misinterpretation as energy-dependent time of flight. Alternatively, intrinsic effects could have the same magnitude in the opposite direction from the LIV effect, canceling it and preventing a detection; a scenario known as a conspiracy of nature. Nonetheless, a LIV effect, if it exists, should be present in all observational data, and directly depend on the distance of a source. A combination of sources at different redshifts would mitigate the potential effect of source intrinsic effects. Furthermore, emission from different types of sources is a result of different physical processes, and subject to different emission dynamics. Therefore, combining different types of sources could further limit the contribution of source-intrinsic effects. These two factors, combining data from different types of sources and from observations at different redshifts, are instrumental in the development of a significantly more robust LIV analysis.
The likelihood function, can be relatively easily extended to consider additional data sets from the same or other sources. Note that a joint likelihood can be constructed from the product of individual likelihoods, each of them described by Equation (10). This single analysis of multiple data sets (sharing the same LIV parameter η n ) is possible as long as each data set is accompanied by its set of functions (light curve template, spectral distribution, acceptance, energy resolution, etc.) describing its particularity and its condition of observations.
Recently, an inter-experiment collaboration has been formed. The so-called LIV Consortium assembles researchers from all three currently operating IACTs experiments (MAGIC, H.E.S.S., VERITAS) with the goal of combining data from different sources. In addition, the LIV Consortium is working on unifying observational data from different IACTs facilities. This will immediately give access to a notably larger pool of sources than what has been used so far in LIV studies with IACTs. Furthermore, combining different instruments in a single analysis, with a particular consideration of individual instrument response functions, is expected to decrease systematic uncertainties. Finally, such combination effort provides the necessary environment to harmonize the details of the analysis. As we have previously seen, analysis techniques can slightly differ from one experiment to the other. Combining individual best practices will further contribute to the research efficiency. Therefore, a combination of observations from all three currently operating IACTs will provide a major improvement in the constraints on the LIV energy scale, not specifically on the facial value of the limits, but more importantly on the robustness of the results. The work done by the LIV Consortium can also be regarded as preparatory activities for the CTA era, which we discuss in the following paragraph. Preliminary results have been presented in conferences [141], while the final results are expected soon.
Another improvement in the direction of combining instruments would be to extend data samples to lower and higher energies, e.g., Fermi-LAT for the MeV-GeV energy range and HAWC or LHAASO for energies above PeV. These experiments provide complementary information to the one of IACTs. It thus makes sense to combine the observations from these experiments to further increase the sensitivity of tests of LIV on gamma rays. Possible benefits are quite tantalizing, although such endeavor would be far from trivial, especially when it comes to treatment of different instrumental effects.
However, near future holds the prospect of the CTA, which will be an order of magnitude more sensitive than any of the existing Cherenkov telescopes [91]. Combined with a large number of telescopes, located in both hemispheres, the CTA will cover larger portions of the sky and observe more sources. This is particularly noteworthy for transient events, such as GRBs and AGN flares, and less bright sources, such as pulsars, all of which are essential for LIV studies. So far, only one LIV study was performed on a GRB observed with an IACT, and only a total of four GRBs have been observed with IACTs until now. The CTA is expected to improve on this statistics. In Section 2.13, we discussed how Vasileiou et al. considered four GRBs in their study. GRB 090510, although located at the lowest redshift, yielded constraints on E QG stronger by a factor of 2-20 in the linear and 3-15 in the quadratic scenarios, compared to results from other three GRBs [38]. This was achieved due to a combination of the highest energy in the sample and the fastest variability. In Section 4 we compared the results obtained on GRB 090510 to the ones from GRB 190114C. The first one was a short GRB, with more than double the redshift of the second one. The signal from GRB 190114C was detected at two orders of magnitude higher energies, however the MAGIC telescopes detected mostly the afterglow phase. With the CTA, we hope to detect prompt emission as well, which is expected to be a more variable phase of GRBs. Furthermore, the CTA will extend the range of accessible energies both to lower (down to 20 GeV) and higher bands (up to 300 TeV). The highest energies detectable with the CTA will be almost an order of magnitude higher than the the ones accessible with the currently operating IACTs. On the other hand, extending the observation window to lower energies, will grant access to gamma rays not absorbed by the EBL, enabling observations of sources at higher redshifts. Referring again to Table 1, higher redshift of a source improves the sensitivity to the LIV energy scale proportionally to the redshift in the linear, and somewhat less in the quadratic contribution. Alas, these two improvements cannot be combined. Gamma rays with the highest energies, emitted from the highest redshifts, will be absorbed by EBL before reaching Earth, so only one of these advantages will be accessible at a time. Nevertheless, there is an improvement coming from the wide energy range itself. Whichever effect of LIV is tested for (time of flight, universe transparency, etc.), the flux at the highest energies is compared to what is assumed to be the intrinsic emission. The latter is estimated from observations in the lower energy bands. LIV, if real, affects the gamma rays in the lower energy band as well the most energetic events. Granted, the effect is smaller for low energy band, but still present. Assuming that low energy photons are unaffected by the LIV induces a bias, and, ultimately, decreases the sensitivity of analyses. The bias will be smaller for a wider gap between the two energy bands used. Observations in the lower energy band can be obtained either using the same instrument, as in the case of e.g., PKS 2155-304 (see [67] and Section 2.6), or from other instruments combined with theoretical inferences, as in the case of GRB 190114C (see [83] and Section 2.12). Ideally, observations in both energy bands would be performed with the same instrument with a wide range of observable energies. The former would reduce possible systematic effects, while the latter would decrease the potential LIV-induced bias. The CTA, with the range of accessible energies spanning over more than four orders of magnitude, will answer this need. While it is difficult to predict the light curves and spectra that the CTA will observe from GRBs, we can draw some conclusions by extrapolating the case of GRB 090510 to higher energies. Assuming that gamma rays emitted during the GRBs prompt phases can reach energies as high as few hundred GeV to a few TeV, the sensitivity to LIV effects would increase by one to two orders of magnitude. Even if we have to settle with lower redshifts, thus somewhat reducing the sensitivity (see Table 1; e.g., for z = 0.3 the loss of sensitivity is at most a factor of ∼3 compared to z = 0.9, the redshift of GRB 090510), the gain would still be substantial. A similar reasoning can be applied to studies based on spectral analysis. By lowering the detection energy threshold, the intrinsic spectra will be more precisely determined at the lower energy end of spectra, leading to better spectral fitting, and thus decreasing the uncertainties on the LIV energy scale. The CTA Consortium has published a projection of the CTA's capabilities for probing fundamental physics, including LIV [120]. The study was limited to universe transparency to gamma rays. The authors estimated that the CTA will be able to probe the LIV energy scale a factor of two to three higher than the most sensitive studies published so far, whether based on a single source, or combining several sources.

Additional and Alternative Lorentz Invariance Violation Effects and Related Phenomena
Apart from the numerous interesting studies described in the previous sections, there are still a plethora of LIV-induced effects and related phenomena which have not been studied with IACTs data. We will briefly mention those here.
Firstly, the energy-dependent arrival time delay between two simultaneously emitted photons from the same source is calculated from Equations (3) and (4). These were derived taking a comoving trajectories of the photons and their respective energy-dependent velocities. As mentioned in Section 2, one could start the derivation from a modified general relativistic dispersion relation [27], or by modifying spacetime translations with the photon dispersion relation [28], and obtain different results for the photon time of flight. It would be interesting to investigate the differences between these models on IACTs data.
Secondly, all the tests of energy-dependent photon group velocity performed on IACTs data were based on Equation (2), which is deterministic in the sense that it assumes that all photons of the same energy will propagate with the same group velocity. However, there are models which propose fluctuations of the photon group velocity as a consequence of fluctuations of the spacetime foam [142]. This phenomenon, often referred to as stochastic LIV models photon group velocity as [143]: where the velocity modification δv γ (E) is randomly distributed according to normal distribution with mean in zero, and a standard deviation given with: The stochastic LIV was tested on Fermi-LAT observation of GRB 090510 [143]. Only linear term was constrained to E QG,1 > 3.4 × 10 19 GeV. So far, no similar study was performed on IACTs observations. As pointed out by Bolmont in [144], distant pulsars with sharp pulsation peaks would be excellent probes of stochastic LIV.
A substantial portion of this work was dedicated to effects LIV has on photon interactions. An important process in astrophysical sources of gamma rays is (inverse) Compton scattering. Abdalla & Böttcher analysed a possible influence of LIV on Compton scattering in [98] and concluded that LIV signatures were expected to be important only for incoming gammaray energies above ∼1 PeV. While in light of the recent results from LHAASO (see Section 3.3) this effect might draw some attention, Abdalla & Böttcher also concluded that even in the superluminal LIV scenario, the Klein-Nishina cross section would still be strongly dominated by the Thomson cross section, while in the subluminal scenario, the Klein-Nishina cross section would be even more strongly suppressed. Overall conclusion was that an LIV modified Compton scattering was not likely to be relevant in realistic astrophysical situations.
We would also like to point out that all investigations of universe transparency for gamma rays were performed on photons of energies up to 100 TeV. Remembering again the LHAASO detection of a photon of 1 PeV, we are strongly encouraged to extend our test to higher gamma-ray energies. As we have demonstrated in Section 3.1, and in particular in Figure 1, this will extend the photon target field to the CMB, which is substantially more dense than any of the EBL components. Moreover, in Figure 2, we have shown how the shape of the Breit-Wheeler cross section deforms in LIV subluminal scenario. These deformations were not significant for gamma-rays up to ∼100 TeV, but might play an important role for higher gamma-ray energies.
As previously mentioned, a study done by the H.E.S.S. Collaboration [81] is the only one in which time of flight and universe transparency tests were performed on the same data set (see Sections 2.11 and 3.1.3). However, these two effects were tested independently of each other. In fact, there has never been a study that combined two different LIV effects. We have argued in Section 4 why strong limits based on one LIV effect do not necessarily constrain other effects. For example, limits on E QG based on the universe transparency do not apply to the photon group velocity. Excluding one effect does not automatically exclude all LIV effects. Nonetheless, considering VHE gamma rays from astrophysical sources, it seems natural to wonder what would be the net observational result if several LIV effects were present. Namely, in Equation (11), term F(E) takes into account propagation effects, such as EBL absorption, but under the assumption that photon interactions were not LIV affected. A combined effect study would presume that both photon group velocity and photon interactions were affected by the LIV. Different effects could manifest on the same, or different energy scales. Since there is no fully formulated theory of QG, and we do not know what the effects of QG are, this test is well justified. Implementation of such a test could be realised by modifying terms in the signal PDF in Equation (11). This might present a significant challenge, especially a computational one when maximising the likelihood function, and in particular if an independent value of E QG is associated to each considered LIV effect.
Finally, there are possible effects of LIV, which have not been mentioned earlier, such as vacuum birefringence, vacuum dispersion, and vacuum anisotropy [145]. The vacuum birefringence implies that the polarization vector of a linearly polarized photon will rotate dependently on the energy of the photon. As a consequence, a linearly polarised signal from an astrophysical source will be depolarised by the time it reaches Earth. Therefore, measurements of polarisation degree in a signal would constrain the vacuum birefringence effect (see, e.g., refs. [146,147] for studies performed on optical observations, and refs. [148][149][150] for studies performed on hard X-ray -soft gamma-ray observations.). These results are substantially more constraining than the tests based on modified photon kinematics. The ratio of sensitivities of the birefringence tests and the time of flight tests, performed on the same data sample, is proportional to the energy of the photons in the sample, as discussed by Kostelecký & Mewes in [151]. As they put it descriptively, in order to achieve sensitivity in the time of flight tests comparable to the birefringence tests, one would need to achieve the timing resolution on the order of the inverse frequency of the photons. Even if this was instrumentally feasible (which is not), the measurement precision would be spoiled by the inability to constrain the photon emission times. Unfortunately, the IACT detection technique does not allow for measurements of gamma-ray polarisation. Recently, a novel satellite-based gamma-ray detector was proposed, which would be capable of measuring polarisation [152] of gamma rays in a lower energy band. Unfortunately, the proposal was not accepted, but there is hope it will be selected for some future space mission, or that a similar concept such as AMEGO 28 will be realised. X-ray polarimetry will certainly gain with soon to be launched IXPE 29 .
Given that IACTs cannot measure the polarisation of gamma rays, a broad discussion of birefringence would be out of the scope of this work. Nevertheless, considering strong constraints on LIV based on birefringence tests, we should note that these do not render time of flight tests useless. Indeed, already from theoretical considerations it is clear that in some cases the energy-dependent photon group velocity does not imply vacuum birefringence [151] (see also [144] for a critical discussion). Therefore, strong limits on the vacuum birefringence do not necessarily constrain the energy dependence of the photon group velocity. Furthermore, while IACTs are not convenient for tests on photon polarisation, one should keep in mind that, as remarked in [74], vacuum anisotropy could be constrained by performing LIV tests on sources in different directions. Hence, we should not be satisfied with setting stronger and stronger constraints on the LIV energy scale only. We should also strive towards building a rich statistics of sources located at different directions, as well as at different redshifts, as we already argued. Again, the CTA is expected to considerably contribute in this respect.

Conclusions
Most of investigators agree that there is a fundamental, quantum description of gravity. Though not formulated yet, the theory of QG is expected to resolve what happens in extreme gravitational potentials, such as singularities within black holes predicted by the general theory of relativity, or early universe, but also to push us in the direction of formulating the next unification theory describing all interactions. The expected realm of QG is the Planck scale, far above the reach of any physical laboratories present today, or any experiment envisaged in the near future. Nevertheless, some investigators believe even VHE gamma rays from astrophysical sources would feel minuscule effects of QG. These effects, consequences of the socalled LIV, would manifest as the photon group velocity deviating from the Lorentz invariant speed of light c, or as modified photon interactions. Given the cosmological distances gamma rays cover to reach Earth, there is a hope that these effects would accumulate sufficiently enough to be detected with gamma-ray detectors.
In this review we presented all experimental tests of LIV performed with IACTs. We followed the historical development of the field. However, in order not to create a confusion, we first covered tests of energy-dependent photon group velocity, and then tests of modified photon interactions. A strictly chronological overview of the results was given in Table  2. 24 years have past since the first proposal that gamma rays emitted from astrophysical objects could be used to search for effects of QG. In fact, the proposal singled out GRBs as gamma-ray sources to be used for these tests. Only two years later, the first study using IACT Whipple was published on data from a blazar Mrk 421. Meanwhile, twenty more years had to pass before we were able to test the LIV on GRB 190114C observed by MAGIC. So far, no significant violation of the Lorentz symmetry was detected. However, in the past 22 years since the first experimental result was published, quite strict bounds were set on the level of LIV. In particular, considering the linear term in modified photon dispersion relation (Equation (1)), the lower limit on E QG,1 has surpassed the Planck energy. This is especially true when photon interactions are considered. While the expected energy level of QG is indeed the Planck scale, there is no strictly defined interval for E QG , so the test will continue. When it comes to the quadratic term in Equation (1), there is still quite some parameter space of E QG,2 to be investigated, both considering the photon group velocity and photon interactions; and in Section 4 we discussed why it is important to test for all possible effects, as well as both for subluminal and superluminal scenarios.
We hope to have demonstrated that IACTs played an important role in search for LIV and setting strong constraints on its energy scale. More importantly, we tried to argue that IACTs still have much to say in this field, and that observations with these instruments will ultimately lead to either detecting LIV, or confirming that the universe remains Lorentz invariant up to trans-Planckian energies. In Section 5 we presented our vision of the future development of this field. It goes without saying that we have great expectations of future facilities. In particular, the CTA, which promises a detection of several GRBs each year. However, new instruments will not do on their own, and additional improvements and development of analysis techniques will be necessary. In addition, there are hypothesised phenomena (e.g., stochastic LIV) which have not been tested for yet. Obviously, there is quite some work cut out for us. Whether we detect some effect of LIV, or it turns out that the Lorentz symmetry is perfectly preserved, one thing is for sure: exciting times are ahead.

Acknowledgments:
We would like to thank G. Bonnoli, C. Levy, M. Martínez, H. Martínez-Huerta, D. Paneque, and F. Tavecchio for fruitful discussions, and verification of historical facts. The authors acknowledge networking support by the COST Action CA18108.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following acronyms and abbreviations are used in this manuscript: