© ESO 2020

1. Introduction

Energy-dependent time-lags in signals arriving from remote cosmic gamma-ray emitters are of particular interest both for understanding the physics of astrophysical sources and for investigating possible new phenomena impacting on photon propagation. Lorentz invariance violation (LIV) is an example of such a phenomenon. It appears as a striking outcome of some Quantum Gravity (QG) models in the form of a modified dispersion relation for photons in vacuum and is one of the most explored ways in QG phenomenology (Amelino-Camelia 2013). It is also included as a starting hypothesis in the standard model extension (SME), an effective field theory built from the Lagrangian of the standard model of particle physics including terms for LIV and charge parity time reversal (CPT) violation (Kostelecký & Mewes 2008, 2009).

Although the present paper focuses on spectral lags induced by a LIV effect in the photon sector, it is necessary to stress that these are not the only possible phenomenon arising from LIV that could be measured with astrophysical sources (see Mattingly 2005, for a review). A modified dispersion relation for photons in vacuum can indeed be interpreted as the photon taking a nonzero effective mass. In that case, normally forbidden processes such as photon decay (γ → e⁺e⁻) and Cherenkov radiation (e⁻ → γe⁻) would be allowed in vacuum, and cross-section for the γ_HE + γ_EBL absorption process of high-energy photons on the extragalactic background light (EBL) would be modified. The latter effect would result in the Universe being more transparent than expected with the standard EBL absorption with no LIV (Biteau & Williams 2015; Tavecchio & Bonnoli 2016; Abdalla et al. 2019). Depending on the QG model considered, it is also possible to obtain vacuum birefringence in addition to energy-dependent speed of light: photons could have different speeds depending on their polarization (Götz et al. 2014; Wei 2019). Several other types of effects are also actively analyzed in the electron sector (Altschul 2005, 2006) and in the gravitational sector (Le Poncin-Lafitte et al. 2016; Bourgoin et al. 2017; Kostelecký & Mewes 2017). Such LIV effects might further modify the launching mechanism of active galactic nucleus (AGN) jets and the radiative processes involved in the generation of blazar flares, but are neglected in the present work.

The first sources proposed to look for LIV spectral lags were gamma-ray bursts (GRBs) because they could be observed at large redshifts and in great numbers by satellites in the soft gamma-ray range (Amelino-Camelia et al. 1998), but flaring AGNs were used almost from the same time (Biller et al. 1999). As in GRBs, AGNs appear to be strong gamma-ray emitters and their active states are observed by detectors with a high effective area such as ground-based Imaging Atmospheric Cherenkov Telescopes (IACTs), with large enough sample sizes to measure fast variability. In addition, AGNs can be monitored regularly and since an AGN flare lasts longer than a GRB, the probability to catch a flare signal under alert is higher. In this regard, blazars emitting in the tera-electronvolt (TeV) range are especially interesting since they are the most variable population of gamma-ray loud AGNs. Their observation with the first generation of IACT, such as Whipple, H.E.S.S., MAGIC, and VERITAS, has already provided stringent limits on spectral time-lags and LIV parameters (Albert et al. 2008; Abramowski et al. 2011, 2015; Zitzer 2013).

The performance aimed for CTA, the Cherenkov Telescope Array (Acharya et al. 2019), will potentially allow for increasingly significant lag measurements thanks to a larger energy range (20 GeV–300 TeV), a higher sensitivity (×10), and a better temporal resolution with respect to the present generation of IACT, and observation strategies designed to optimize the number of transient or variable object detections. The measured lags, if significant, will have to be interpreted as propagation delays, as effects intrinsic to the sources, or as a superposition of both.

The present work is a first attempt to gain knowledge on source-intrinsic spectral lags of flaring AGNs at high and very high energies and on short timescales relevant for LIV searches, using leptonic AGN flare modeling. In this study, only the LIV effect on the propagation of photons in vacuum and the subsequent time-delays engendered are considered to keep a relatively simple AGN flare model. Section 2 presents the general context of the search for LIV signatures from the analysis of blazar gamma-ray flares. We briefly present and apply a standard Synchrotron-Self-Compton (SSC) scenario for such flares in Sect. 3, focusing on purpose on the dominant and unavoidable mechanisms needed to generate the burst. We explore the intrinsic SSC time delays which are induced in the gamma-ray range in Sect. 4, and further characterize them in Sect. 5. Section 6 focuses on the VHE domain (E >  100 GeV) to compare potential intrinsic and LIV time-lags, and further astrophysical issues are discussed. Conclusions and perspectives are presented in Sect. 7.

2. Search for Lorentz invariance violation from spectral lags in blazars

2.1. Spectral lags and Lorentz invariance violation

Focusing on spectral lags in the context of QG phenomenology, and neglecting birefringence effects, the modified dispersion relation is usually expressed as

$$\begin{array}{c}\hfill E^2\simeq p^2{c}^2\times [1\pm {\left(\frac{E}{E_{\mathrm{QG}}}\right)}^n],\end{array}$$(1)

where c is the low-energy limit of the speed of light, E_QG is the energy to be measured or constrained at which LIV effects should become non-negligible and the sign ± translates into the possibility to have an increasing (superluminal) or decreasing (subluminal) speed when photon energy increases. The value of E_QG is usually expected to be of the order of the Planck scale E_P ∼ 10¹⁹ GeV. Present day IACT experiments probe the linear effect (n = 1) while the quadratic effect (n = 2) is still far from reach. However, other orders are investigated at VHE with other experiments (see Kostelecký & Russell 2011, and references therein). Moreover, for odd values of n, there is a correspondence between the photon helicity and the subluminal or superluminal case. In this first work, we neglect this helicity effect and assume only either the subluminal or superluminal case for odd n values.

The use of variable or transient and distant astrophysical sources for LIV searches was first proposed in the late 90s (Amelino-Camelia et al. 1998). The modified dispersion relation of Eq. (1) naturally leads to an energy-dependent group velocity of light. Two photons of different energies (E_h and E_l, with E_h >  E_l), assumed to have been emitted at the same time from the same place by a source at redshift z would be detected with a spectral lag

$$\begin{array}{c}\hfill \mathrm{\Delta }{t}_{\mathrm{LIV},n}\simeq \pm \frac{n+1}{2}\frac{E_h^n-E_l^n}{E_{\mathrm{QG}}^n}{\kappa }_n(z),\end{array}$$(2)

where the distance parameter

$$\begin{array}{c}\hfill {\kappa }_n(z)={{\displaystyle \int }}_0^z\frac{{(1+z\prime )}^n}{H(z\prime )}\mathrm{d}z\prime \end{array}$$(3)

is an increasing function of redshift taking into account Universe expansion during photon propagation (Jacob & Piran 2008). Here, H(z) is the Hubble parameter. It is necessary to point out that this expression, although used in all LIV searches performed so far, was obtained under the implicit assumption that translations are not affected by Planck scale effects. Other expressions have been proposed, for example in the deformed special relativity (DSR) approach (see e.g., Rosati et al. 2015). From Eq. (2), it is common to express the “time-lag over energy difference” parameter τ_n such that

$$\begin{array}{c}\hfill {\tau }_n\equiv \frac{\mathrm{\Delta }{t}_{\mathrm{LIV},n}}{E_h^n-E_l^n}·\end{array}$$(4)

This parameter is constrained and limits on E_QG are derived from astrophysical source observations.

The lags induced by a LIV effect are expected to be small. Maximizing them requires observation of sources preferably with high redshifts and hard spectra on a wide energy range so that κ_n and $E_h^n-E_l^n$ are maximized. Fast variability is also required in order to measure the lags. The high-energy gamma-ray domain is therefore particularly suitable for such studies.

The expression of Eq. (2) was obtained assuming that high- and low-energy photons are emitted at the same time from the same place, that is, neglecting any source-intrinsic delays possibly resulting from emission mechanisms and source extent. In principle, the measured time-lag Δt_m should rather be expressed as the sum of delays with different origins:

$$\begin{array}{c}\hfill \mathrm{\Delta }{t}_{\mathrm{m}}=\mathrm{\Delta }{t}_{\mathrm{LIV}}+(1+z)\mathrm{\Delta }{t}_{\mathrm{s}}+{\displaystyle \sum _j}\mathrm{\Delta }{t}_j,\end{array}$$(5)

where Δt_s is the delay due to the emission processes at the source located at redshift z and Δt_j accounts for various additional effects which could affect Δt_m. These terms include for example the dispersion by free electrons along the line of sight, mostly important in the radio range, potential lags due to special relativistic effects in the case where photons have a nonzero rest mass, or lags caused by the gravitational potential integrated from the source to the Earth (see Gao et al. 2015; Wei et al. 2016, and references therein). These extra terms Δt_j are neglected in the present work.

2.2. Intrinsic time-lags from a blazar flare

In the gamma-ray domain, only one flare from Mrk 501 has shown an indication of a 4 ± 1 min time-lag between energy bands below 250 GeV and above 1.2 TeV (Albert et al. 2007). Barely significant spectral lags τ₁ = (0.030 ± 0.012) s GeV⁻¹ and τ₂ = (3.71 ± 2.57) × 10⁻⁶ s GeV⁻² were later reported from the same data set (Albert et al. 2008). This flare, recorded on July 9 2005 by the MAGIC Cherenkov telescopes, suggests that intrinsic delays can exist in AGN flares, while the fact that it is the only one detected implies that intrinsic effects are certainly different for each AGN and perhaps from one flare to another, even in the same source.

Apart from the previous example, no significant spectral lag has ever been measured at GeV and TeV energies from AGNs, and source intrinsic effects have been ignored when constraining the LIV energy scale. Neglecting intrinsic effects could be partially justified when the energy range considered in the analysis is restricted enough to ensure that observed photons can be considered as all emitted together at once in the cosmic source. However, using stringent energy selections has a drawback since it drastically decreases statistics.

Delays can indeed easily appear during the emission of photons from cosmic sources. Such intrinsic spectral time-lags have already been unambiguously detected in some GRBs, and have mostly been reproduced by models considering light-travel-time effects from extended and inhomogeneous emitting zones as expected in standard GRB scenarios (Dai et al. 2017, and references therein). In the case of blazar flares, intrinsic effects do not appear to be as important and have been poorly detected up to now. Nevertheless, they should be detected, either soon with current instruments during some extraordinary flare fully monitored over a large spectral gamma-ray range, or in the coming years with a new generation of instruments providing high-quality light curves and dynamical spectra at VHE. Therefore, in the context of LIV signature searches, intrinsic effects need to be further investigated, at least to help in constraining the QG energy scale in case of future detections of significant time delays in AGN flares (Perennes et al. 2017).

Several time-dependent emission models of nonthermal emission in AGNs have been proposed in the literature, such as for example by Blandford & Königl (1979), Marscher & Gear (1985), Celotti et al. (1991), Böttcher et al. (1997), Katarzyński et al. (2003), Joshi & Böttcher (2011), and Lewis et al. (2016). As a matter of fact, such scenarios usually induce possible spectral lags but such intrinsic lags have not been purposefully analyzed in the context of TeV blazar flares. We distinguish here two different types of time delays, the “macroscopic” and the “microscopic” ones, depending on their values and origins.

Macroscopic delays correspond to a variety of long lags induced by the global structure and properties of extended nonhomogeneous emitting zones. Different parts of the source emit in different spectral ranges, therefore inducing possibly important time-lags depending on the specific geometry and kinetics of the radiating plasma. This type of scenario was proposed for instance by Sokolov et al. (2004), while triggering an AGN flare by the collision of a stationary shock wave with a relativistic shock wave in the jet, later forming reverse and forward shocks which both accelerate particles. The complex jet geometry as well as acceleration and emission pattern which then develop can result in long time-lags. This model successfully reproduced a flare observed in 3C 273, in which X-rays were delayed by about one day with respect to their infrared counterparts (Sokolov et al. 2004). Another example of “macroscopic delay” has been obtained by Bednarek & Wagner (2008) who explained the time-lag of about four minutes measured in the 2005 Mrk 501 flare by considering a continuous increase of the global Lorentz factor of the emitting zone propagating along the jet. In such a model of an accelerating blob, with an increase of the Doppler boosting effect, lower-apparent-energy photons are, on average, emitted in a different location from higher-apparent-energy photons, creating the observed spectral time delay.

In the present study we focus on microscopic delays, which are barely discussed at TeV energies although they are easily produced in standard one-zone SSC models and naturally come from the temporal evolution of the distribution of emitting particles. Such microscopic delays have already been analyzed by Lewis et al. (2016) in the context of X-ray variability studies in blazar jets. Focusing on synchrotron emission, these latter authors produced flare light curves and their associated X-ray time-lags. Applying them to data from Mrk 421, they were able to characterize some acceleration parameters in that source. However, Lewis et al. (2016) entirely neglect the inverse-Compton losses which are mandatory to describe gamma-ray time delays. In the following section, we consider a minimal time-dependent leptonic model to generate blazar flares and study short intrinsic time delays at gamma-ray energies, above 1 MeV.

3. Time-dependent model

3.1. Description of the electron distribution

The time evolution of a flare is deduced from the description of the electron population N_e(γ, t) at the source in a single homogeneous emitting zone, with γ = E/(m_ec²) being the Lorentz factor of electrons. We adopt the standard SSC scenario presented in Katarzyński et al. (2001, 2003) but simplify the geometry to focus on gamma-ray flares arising from a small spherical blob of plasma. Further developments of the simulation codes were done to adapt them to a detailed analysis of spectral lags. The evolution of the electron density is expressed as:

$$\begin{array}{c}\hfill \frac{\partial N_{\mathrm{e}}^{\ast }(t,\gamma )}{\partial t}=\frac{\partial }{\partial \gamma }\left\{[{\gamma }^2{C}_{\mathrm{cool}}(t)-\gamma C_{\mathrm{acc}}(t)]N_{\mathrm{e}}^{\ast }(t,\gamma )\right\},\end{array}$$(6)

where C_cool accounts for radiative energy losses and C_acc for acceleration of electrons or energy gain (neglecting here adiabatic losses). This differential equation admits an analytical solution with an initial electron spectrum N_e(0, γ) defined as a power law function with a cut-off:

$$\begin{array}{c}\hfill N_{\mathrm{e}}(0,\gamma )=K_0{\gamma }^{-n}[1-{\left(\frac{\gamma }{{\gamma }_{\mathrm{c},0}}\right)}^{n+2}],\end{array}$$(7)

where K₀ is the electron density at γ = 1, n the electron spectrum index, and γ_c, 0 the cut-off energy.

The coefficient C_cool takes into account the SSC energy losses and can be expressed as follows:

$$\begin{array}{c}\hfill C_{\mathrm{cool}}(t)=\frac{4{\sigma }_T(U_B(t)+U_r(t))}{3m_{\mathrm{e}}c},\end{array}$$(8)

where the first term corresponds to synchrotron energy losses, with U_B(t) = B(t)²/8π being the magnetic field energy density, and the second term is related to the inverse-Compton (IC) energy losses on the synchrotron photon field, corresponding to the synchrotron energy density U_r(t). It is parameterized as:

$$\begin{array}{c}\hfill U_r(t)=U_B(t)/\eta .\end{array}$$(9)

This choice allows us to simplify the resolution of Eq. (6), which then admits an analytical solution, and to decrease the computation time significantly. The η parameter represents the relative importance between the synchrotron and IC radiative cooling. It has to be constant over time or to be large enough in order for the analytic solution to be valid. This limits the parameter space that can be explored to the one corresponding to synchrotron-dominated sources like high-frequency peaked BL Lac (HBL) objects.

We further assume a characteristic time t₀ for the temporal evolution of the main flare parameters. This time is typically related to the speed of sound in a blob of relativistic plasma characterized by a radius R₀ as $t_0=R_0/(c\sqrt{3})$, c being the speed of light in vacuum. The evolution of the magnetic field strength is then described as:

$$\begin{array}{c}\hfill B(t)=B_0{\left(\frac{t_0}{t}\right)}^{m_b},\end{array}$$(10)

with B₀ being the initial value of the magnetic field at t₀, and m_b its temporal evolution index.

Finally, C_acc corresponds to the acceleration processes and allows us to initiate the flare starting from low states of the electron spectrum distribution. It is expressed as:

$$\begin{array}{c}\hfill C_{\mathrm{acc}}(t)=A_0{\left(\frac{t_0}{t}\right)}^{m_{\mathrm{a}}},\end{array}$$(11)

where A₀ is the initial acceleration amplitude, and m_a the acceleration temporal index.

3.2. Spectral energy distribution

Spectral energy distributions (SEDs) are generated assuming a SSC emission model. To compute the SSC radiation, we use standard packages described in Katarzyński et al. (2001) which provide the synchrotron and IC emission from a given electron spectrum, taking into account synchrotron self absorption of the emitting zone. Absorption from the extragalactic background light (EBL) is included following the model of Franceschini et al. (2008). The evolution of the electron spectrum and of the SED is shown in Fig. 1 for the set of parameters given in Table 1. Such physical parameters correspond to standard values as previously deduced with the same SSC packages and others for the archetypal TeV sources Mrk 501, Mrk 421, and PKS 2155−304 (Katarzyński et al. 2001, 2003; Abramowski et al. 2012) and can be considered as typical values expected in blazar flares. This explains the small redshift value adopted here, leading to a low EBL absorption.

Fig. 1. Electron spectrum (top) and SED (bottom) evolution for the reference set of parameters given in Table 1, illustrating a typical TeV blazar flare. The two plots share the same color code. The full lines correspond to the rising phase of the flare and the dashed lines to its decay. The time tmax is defined as the time when the highest energy value γmax = max[γc(t)] is reached in the electron spectrum. The corresponding electron spectrum and SED are plotted in thick lines. The electron spectrum and SED at the maximum luminosity of the flare are not displayed because they almost coincide with the ones at tmax for this reference set of parameters.

3.3. Light curve determination

To compute the light curves at different energies, the SEDs are integrated over the required energy bands, each SED giving one flux point for each light curve. Examples of flare light curves are shown in Fig. 2 (left and center) for parameters of Table 1. In order to consider only realistic light curves that may be observed, a selection cut is applied on them. Only light curves with a maximum flux above 2 × 10⁻¹⁵ cm⁻² s⁻¹ are kept. The time t_max shown on the light curve is defined from Fig. 1 as the time when electrons reach their maximum energy γ_max = max[γ_c(t)]. This time indicates when the radiative cooling timescale becomes shorter than the acceleration timescale, for the most energetic electrons with an energy close to γ_max. For this set of parameters, the time t_max happens after all the light curves peak, meaning that the electron acceleration timescale is still shorter than the electron cooling timescale for all energies. As a result, the decaying light curves can only be explained by the decrease of the magnetic field. An opposite case can be defined when t_max happens before all the light curves peak as shown in Fig. 3 (left and center). This case is obtained from the first reference set of parameters increasing the radiative cooling power. This is achieved by changing the B₀ value from 65 mG to 90 mG, and provides us with our second reference set of parameters throughout this paper. Conversely to the previous case, the radiative cooling contributes to the light-curve decay, resulting in the fact that the most energetic light curves decay first.

Fig. 2. Normalized light curves (left), zoom close to the light-curve maxima (center), and time-delays between each light curve and the one in the range 1–1.8 MeV (right) obtained for the first reference set of parameters. The vertical dashed line on the light curves (left) corresponds to the time when electrons reach their maximum energy γmax and shows the moment when energy losses dominate over the acceleration at the energy γmax.

Fig. 3. Normalized light curves (left), zoom close to the light curve maxima (center), and time-delays between each light curve and the one in the range 1–1.8 MeV (right) obtained for the second reference set of parameters (i.e., with B = 90 mG). The vertical dashed line on the light curves (left and middle) corresponds to the time when electrons reach their maximum energy γmax and shows the moment when energy losses dominate over the acceleration at the energy γmax.

4. Study of time delays

4.1. Time delay determination

The time delay is determined by computing the time difference between the maximum of the light curve at the energy considered and the maximum of the lowest-energy light curve (1.0–1.8 MeV). A positive lag corresponds to the case when high-energy light curves peak after the lowest-energy ones. This method, later called the peak position method (PPM), is rather simplistic but provides a simple and robust estimation of the time delay, with an accuracy directly related to the time step chosen for the light curves. A cross-correlation function (CCF) from Edelson & Krolik (1988) was also considered but was found to incorrectly reconstruct the time delay in the case of light curves with varying widths. A comparison of the PPM and CCF is shown in Fig. 4. Two light curves are simulated with a similar shape as the light curves obtained from the model, approximated by an asymmetric Gaussian function. A delay of 800 s is applied to one of the light curves. The width of the lagged light curve is then varied with respect to the other one and the PPM and CCF methods are applied to reconstruct the lag. The CCF is not able to reconstruct the injected delay as the width difference increases, and more generally in case the shape of the two light curves is different. Concerning the light curves from the model, this point is further discussed in Sect. 6.3.

Fig. 4. Time delay measured with PPM and CCF between two simulated asymmetric Gaussian light curves varying the width of one of them. The injected delay value is fixed at Δtinj = 800 s. The measurement from the CCF shows a misreconstruction of the injected time delay as the width difference increases. The asymmetric Gaussian function used has σleft = 2500 s and σright = 20 000 s.

4.2. Main results: two time-delay regimes

The time delays obtained for the two reference sets of parameters shown in Figs. 2 (right) and 3 (right) highlight two different time-delay regimes found with the adopted SSC flare model. These two regimes can be better characterized by considering whether the time t_max is reached before or after the peaks of the light curves. The vertical error bars on the time-delay measurement is of 5 s and the horizontal ones for the energy correspond to the light-curve energy ranges.

In the first case (Fig. 2), all the light curves peak before t_max. Thus, the acceleration timescale is shorter than the cooling one when the flare decays. Electrons are still accelerated and their maximum energy still increases after the low-energy light curves peak. The decay of the flare is then induced by the decrease of the magnetic field B(t). Hence, the increase of the time delay above 1 GeV is explained by the time required for electrons to be accelerated up to γ_max and to emit the highest-energy photons. Conversely, below 1 GeV, the energy necessary to emit lower-energy photons is quickly reached by electrons. This leads to a decrease of the time delay, the low-energy light curves decaying due to the combined action of the energy-dependent radiative cooling and of the magnetic field decrease. In the following, this type of time delay evolution, with a decreasing time delay at low energies and an increasing one at high energies, is described as an “acceleration-driven regime” because of the origin of the high-energy delays.

In the second case (Fig. 3), t_max is reached before the peaks of the light curves, and the energy losses by radiative cooling are soon larger than the energy gains by acceleration. The decay of the flare is then mostly due to the radiative cooling. This is achieved by assuming a higher magnetic field strength initial value (B₀ = 90 mG) relative to the previous case, which enhances the electron radiative energy losses. As a consequence, the highest-energy light curves decay first due to a shorter cooling timescale, explaining the decreasing time delay above GeV energies. At lower energies, the cooling timescale of the electrons and the decreasing time delay appear similar to the previous case and due to the combined action of the magnetic field decrease and of the energy-dependent radiative cooling. In the following, such a case with time delays continuously decreasing are described as the “cooling driven regime” because of the main origin of the time delays above GeV energies.

The analysis presented in the following section confirms the existence of these two time-delay regimes over a large domain of parameters.

5. Influence of model parameters on time delays

To investigate the impact of the model parameters on the time delays, each of them is individually varied around the reference values defined in Table 1, namely B₀, m_b, A₀, m_a, δ, γ_c, 0 and n. All these values are chosen within a reasonable range for blazar modeling and also in the domain of validity of our model. Magnetic fields being a key driver for nonthermal blazar emission and SSC electron radiation, several values of B₀ are first explored and then each other parameter is investigated for B₀ = 65 mG and B₀ = 90 mG respectively.

5.1. Initial magnetic field strength variations

The magnetic field directly influences the cooling timescale of electrons. Therefore, it contributes to the maximum energy γ_max reached by electrons, meaning that B₀ is one of the main parameters acting on the intrinsic time delays.

Several values of B₀ are used to evaluate the influence of the magnetic field on the time-delay evolution, ranging from 50 mG to 110 mG. The resulting time delay is shown in Fig. 5 as a function of energy. In the low-energy domain, for E ≲ 1 GeV, the time delay does not qualitatively change with B₀ although it slightly decreases as B₀ increases due to a stronger radiative cooling. Conversely at high energies, for E ≳ 1 GeV, the variation of the initial magnetic field value induces large variations of the time delay and a significant change of its behavior, with a transition zone around 80 mG from an increasing to a decreasing phase when B₀ increases. Within this transition, some cases show an almost constant time delay relative to the MeV light curve, which means that they present just very small or even no relative time delay in the limited energy range from 1 GeV to 1 TeV.

Fig. 5. Time delay vs. energy for different B0 values. The two cases in bold in the legend correspond to the ones discussed in Sect. 4.2. All other parameters have the values given in Table 1.

5.2. Magnetic field temporal index

As defined in Eq. (10), the magnetic field strength is assumed to decrease over time with an index m_b. A value m_b = 1 was found adequate to reproduce flares observed in Mrk 421 within the SSC scenario adopted here (Katarzyński et al. 2003), and m_b = 2 corresponds for instance to the case of magnetic flux conservation in a blob propagating along a quasi-conical jet. Here, the index m_b is varied from 1 to 2 and time delays are shown for the two reference cases in Fig. 6.

Fig. 6. Time delay vs. energy for different magnetic field evolution index mb with B0 = 65 mG (left) and B0 = 90 mG (right). All other parameters are unchanged (Table 1). The two cases in bold in the legend correspond to the situation discussed in Sect. 4.2.

Starting from an acceleration-driven regime in the case where B₀ = 65 mG, an overall increase of the time delay value is observed. For E ≳ 1 GeV the time delay evolution does not qualitatively change and lags always increase. For E ≲ 1 GeV, the time delay decreases less and less as m_b increases and finally evolves towards only increasing time delay. This transition is explained by the rapid decrease of the magnetic field which induces low-energy flares quickly decaying while electrons are still accelerating. As a consequence, the highest-energy light curves peak at later times due to the time needed for the electrons to reach high γ values.

Starting from a cooling driven regime in the case where B₀  =  90 mG, increasing m_b induces a transition to the acceleration-driven regime. In those cases, the rapid decrease of B(t) reduces the radiative cooling power. Hence the acceleration time scale becomes shorter than the radiative cooling timescale, thus shifting t_max to later times. When m_b is larger, the highest-energy light curve peaks later while the faster magnetic field decrease induces the lowest-energy light curves to decay earlier.

5.3. Doppler factor variations

A modification of the Doppler factor δ does not change the temporal evolution of the source in its own frame and simply implies a change on the Doppler boosting effect for the observer leading to variations of the observed variability due to time contraction as well as variations of the observed flux and energy. The energy-dependent time delay for different δ ranging from 20 to 50 is shown in Fig. 7 for the two regimes. The dominant effect on the time-delay evolution is the time contraction, inducing smaller observed time delays as δ becomes larger. In reality, the time delays obtained for different δ values appear almost proportional to each other by the ratio of their Doppler factor, as shown for instance for B₀ = 90 mG by the delay at 1 TeV for δ = 40 and δ = 20 which gives Δt_δ = 40(1 TeV)≃293 s and Δt_δ = 20(1 TeV)≃600 s, with a ratio approximately equal to the Doppler factor ratio. However, it can also be noted that the variation of the maximum energy considered for the time-delay computation is an effect of the Doppler boosting on the apparent flux value. This is a consequence of our choice to neglect light curves with flux values below 2 × 10⁻¹⁵ cm⁻² s⁻¹. Also, the energy shift due to Doppler boosting can lead to a change of the time delay sign which is not expected from the time contraction alone as seen for B₀ = 65 mG at E ≈ 1 TeV.

Fig. 7. Time delay vs. energy for different Doppler factor δ values with B0 = 65 mG (left) and B0 = 90 mG (right). All other parameters are unchanged (Table 1). The two cases in bold in the legend correspond to the situation discussed in Sect. 4.2.

5.4. Acceleration parameter variations

The acceleration term, defined in Eq. (11), depends on two parameters, namely the initial amplitude A₀ and the evolution index m_a. Varying the acceleration parameters modifies the electron acceleration time scale and thus the time t_max with respect to the flare maxima. This can in turn induce a transition between the two time-delay regimes.

The time delay obtained for different A₀ values ranging from 4.0 × 10⁻⁵ s⁻¹ to 6.0 × 10⁻⁵ s⁻¹ is shown in Fig. 8. In the case where B₀ = 65 mG, decreasing A₀ increases the time delays for the highest energies. The acceleration power is weaker and high-energy electrons need more time to reach γ_max. In addition, the maximum electron energy γ_max is smaller, implying a weaker cooling effect which shifts t_max to later times. With larger A₀ values, a transition between the two regimes occurs, leading to a cooling-driven regime. The acceleration is then stronger involving a higher γ_max value, inducing a stronger cooling effect for the most energetic electrons. Hence, increasing A₀ brings electrons much faster to higher γ_max values which then quickly suffer from intense radiative energy losses, entering into the cooling-driven regime.

Fig. 8. Time delay vs. energy for different acceleration amplitude A0 values with B0 = 65 mG (left) and B0 = 90 mG (right). All other parameters are unchanged (Table 1). The two cases in bold in the legend correspond to the situation discussed in Sect. 4.2.

Conversely, starting from a cooling-driven regime with B₀  =  90 mG, the opposite situation occurs with a transition to the acceleration-driven regime when A₀ decreases. Small A₀ values induce longer acceleration timescales and lower γ_max values, the cooling power at γ_max becomes weaker, electrons need a long time to reach their maximum energy, the t_max is shifted to later times, and a transition occurs to the acceleration-driven regime when A₀ is small enough. On the other hand, larger A₀ values imply a shorter acceleration timescale, γ_max is reached at earlier times, and the cooling power becomes stronger. Then the highest-energy light curves reach their peaks much earlier than the reference case, which explains the negative time delays.

Similar reasoning applies to the acceleration evolution index parameter m_a (Fig. 9), varied from 4.5 to 5.9. Increasing m_a induces longer acceleration timescales. However, the variation of m_a does not lead to a significant change of regime. For B₀ = 90 mG, a hint of transition is observed for large m_a values when the acceleration is weaker. For B₀ = 65 mG, the overall delay decreases toward negative values for a stronger acceleration when m_a decreases.

Fig. 9. Time delay vs. energy for different acceleration evolution ma values with B0 = 65 mG (left) and B0 = 90 mG (right). All other parameters are unchanged (Table 1). The two cases in bold in the legend correspond to the situation discussed in Sect. 4.2.

5.5. Initial electron distribution index variations

The initial electron spectrum assumed in the present flare scenario follows a power law function with a high-energy cut-off Eq. (7). The initial electron density K₀ is only a scaling parameter and does not affect the time evolution. In the transfer equation, modifying n does not change the balance between electron acceleration and cooling effects and in practice γ_max and t_max remain at the same values when n changes. However, the initial electron spectrum index n impacts the ratio of low- to high-energy electrons, with a higher proportion of lower-energy electrons for high n values.

The time delays obtained with values from n = 2.2 to n = 2.8 are presented in Fig. 10. For B₀ = 90 mG, starting from a cooling-driven regime, a transition occurs to the acceleration-driven regime when n decreases. Such a transition can be easily explained because the flare is globally shorter when n is smaller while t_max remains the same whatever the value of n. This is the consequence of the fact that the electron population is on average more energetic for smaller n, therefore inducing light curves which peak at earlier times. The transition occurs when n is small enough to produce light curves peaking earlier than t_max. For B₀ = 65 mG, the variations of the time delay are small and do not really highlight the influence of the parameter n, however the evolution of the time delay follows the same behavior. A smaller n value induces an overall shorter flare, with light curves peaking at earlier times because they are produced by a more energetic electron population. The maxima of the light curves become shifted at earlier times than t_max thus inducing larger delays for the highest-energy light curves.

Fig. 10. Time delay vs. energy for different initial electron index n values with B0 = 65 mG (left) and B0 = 90 mG (right). All other parameters are unchanged (Table 1). The two cases in bold in the legend correspond to the situation discussed in Sect. 4.2.

5.6. Energy break of electron distribution evolution

The energy cut-off γ_c, 0 defines the maximum energy of electrons at the starting time t₀ Eq. (7). Using lower values of γ_c, 0 than 4 × 10⁴ increases the time needed for electrons to reach their highest energy γ_max. For higher γ_c, 0 values, electrons are quickly accelerated to high γ values and γ_max become larger leading to a shorter cooling timescale. The resulting time delays for γ_c, 0 ranging from 2 × 10⁴ to 8 × 10⁴ are shown in Fig. 11 for the acceleration- and cooling-driven cases.

Fig. 11. Time delay vs. energy for different initial maximum energy for electron γc, 0 values with B0 = 65 mG (left) and B0 = 90 mG (right). All other parameters are unchanged (Table 1). The two cases in bold in the legend correspond to the situation discussed in Sect. 4.2.

For B₀ = 65 mG, decreasing γ_c, 0 leads to an increase of the time delay values for all energies. Indeed, starting from less energetic electrons implies that they need more time to be accelerated up to γ_max, thus shifting the highest-energy light curves to later times. For larger γ_c, 0 values, the electron population is more energetic at t₀ and reaches a larger γ_max within a smaller time t_max. In this situation, the cooling timescale becomes shorter at γ_max leading to a transition from the acceleration- to the cooling-driven regime. For B₀ = 90 mG, a similar behavior appears.

6. Discussion and astrophysical issues

With the flare model presented here, all the cases investigated reveal the presence of an energy-dependent intrinsic time delay at gamma-ray energies. Two distinct regimes, referred to as “cooling driven” and “acceleration driven” are found for the time delays, corresponding to the mechanism driving the electron evolution when the light curves peak and the flare starts to decay. In addition, some specific cases corresponding to the transition between the two regimes show no delay between roughly 100 GeV and a few TeV. Adiabatic effects due to the expansion of the emitting zone can affect the quantitative values of time delays but do not qualitatively modify the global picture. When added in the differential equation describing the evolution of electrons Eq. (6), they contribute to reducing the acceleration term, therefore pushing the system towards the acceleration-driven regime, which will be reached or not depending on the given set of parameters.

6.1. Observational constraints on the model

The information on the energy-dependent time delay can be used in order to constrain the model parameters. Indeed, if the evolution of the measured time delay corresponds to one of the two regimes, the other one is obviously ruled out. Clearly, the time delay is a new observable which can be used to constrain the modeling of blazars.

For instance, the observation of a flare from the blazar Markarian 501 in 2005 (Albert et al. 2007) revealed a nonzero time delay in the VHE range, the unique case of detection of a time delay at VHE from a blazar. The authors found a time delay increasing with the energy which appears to correspond to an acceleration-driven regime. This suggests a qualitative scenario with a flare initiated through a sudden shock acceleration or magnetic reconnection in the emitting zone, immediately followed by a mechanism inducing the flare decay such as the decrease of the magnetic field (or adiabatic expansion). In parallel, acceleration processes are still more efficient than the radiative cooling at the highest energies and thus ensure the observed acceleration-driven regime.

However, most of the blazar flare observations do not show any significant time delays. One example is provided by the exceptional flare of PKS 2155-304 observed by H.E.S.S. in 2006. A CCF was applied to the light curves (Aharonian et al. 2008) between 400 and 800 GeV and above 800 GeV, and no significant delay was found in the data. The reason for such missing time-lags is not clear, because basic SSC flare models such as the one presented here predict significant intrinsic delays. Obviously this can be due to the limited number of observed blazar TeV flares, the poor time coverage and time accuracy, and the limited spectral range of the present data sets. The next generation of gamma-ray instruments should clarify this issue by measuring significant time delays and providing new and precise quantitative constraints on VHE flare models. Indeed, a flare which can only decay through radiative cooling leads inevitably to a cooling-driven regime. Therefore, for example, if future observations of blazar flares do not reveal the presence of time delays decreasing with the energy at VHE, the simplest scenarios with a fast acceleration (or injection) followed by radiative cooling in otherwise constant emitting zone and magnetic field could be excluded. Alternatively, it is also possible that basic SSC scenarios do not describe blazar flares in full details. If time delays are not at all observed at the VHE, then the majority of flares may occur in a specific domain of parameters corresponding to the transition zone between the two time-delay regimes we have identified. This would suggest a physical link and a fine-tuning between acceleration and cooling processes in the global evolution of the flares.

6.2. Focus on time delays at very high energies

For the search of LIV signatures with IACT, studies are performed in the VHE range only. The intrinsic delays obtained in Sects. 4 and 5 are therefore re-calculated here relative to a reference light curve at higher energies, in the range from 42 to 74 GeV which is now used as the zero origin for the delays. Comparing the energy dependency of the intrinsic time delays with the LIV ones can then provide direct constraints on specific QG models since for instance some models produce only one specific type of time delay such as Amelino-Camelia et al. (1997) or Ellis et al. (2000) with positive, linear, energy-dependent delays.

To quantify the energy dependency of the intrinsic delays, they are adjusted with a power law function similar to the one used for LIV studies, namely

$$\begin{array}{c}\hfill \mathrm{\Delta }t=\xi \times (E_i^{\alpha }-E_0^{\alpha }),\end{array}$$(12)

where α is the energy dependency index, ξ the amplitude of the delay in s TeV^−α, and E₀ the midpoint of the energy range of the reference light curve (58.5 GeV). The α and ξ parameters obtained from fitting the time delay for the cases with different B₀ values are shown in Table 2 and for all other cases in Table A.1. At the transition between the two regimes, α cannot be evaluated since there is no significant time delay. In addition, some cases are not able to produce a significant flare above 250 GeV, implying a poor energy coverage and preventing any estimation of the α parameter. Overall, the energy evolution index α for the cases producing significant time delays is found to be in the range [0.45−0.85].

As a consequence, energy-dependent intrinsic time delays obtained from the basic SSC flare model described here evolve with an index α different from the QG model predictions from Amelino-Camelia et al. (1997) and Ellis et al. (2000) where α_LIV = 1. This result illustrates how time-dependent blazar flare scenarios can be used to test these two QG models or any other model predicting a LIV delay with an energy dependency outside the range of values found for the α parameter. Indeed, if a specific QG model predicts a LIV energy-dependent delay with an index α_LIV outside this range, the two delays can be discriminated and the QG model can be constrained. Otherwise, the LIV and intrinsic delays are mixed together and remain difficult to disentangle.

To summarize, Fig. 12 shows the domain of intrinsic time delays generated in the two regimes and a linear LIV delay at the Planck energy scale (E_QG ∼ 10¹⁹ GeV) for subluminal and superluminal effects at different redshifts adopting the redshift evolution from Jacob & Piran (2008). A direct comparison can be done since the redshift does not affect the time-delay evolution with energy but only the observed flux due to the distance and EBL attenuation. Generally speaking, intrinsic time delays appear to be much larger than expected LIV delays in the linear case, except for very large redshifts where observations at VHE remain difficult due to EBL absorption. For higher values of E_QG, the LIV delay will simply become smaller. In the quadratic case, LIV delays will always be much smaller than intrinsic ones and very difficult to disentangle.

Fig. 12. Typical time delays obtained from the SSC flare model (shaded area) and expected from a linear LIV effect with an energy scale at the Planck energy (full lines) with the redshift evolution found in Jacob & Piran (2008). Other formalism such as DSR would lead to different values but of the same order of magnitude. The dashed line corresponds to Δt = 0. Positive and increasing LIV delays correspond to subluminal LIV effects, negative and decreasing ones to superluminal LIV effects.

6.3. Temporal evolution

A characteristic feature of intrinsic delays is that their magnitude can vary in time during a flare while LIV delays stay constant. This is another observational signature which can provide important information on the origin of the delays. In their study of X-ray variability of blazar flares, Lewis et al. (2016) performed a Fourier transform analysis of the delay between two light curves at different energies and obtained the time delay as a function of the Fourier frequency, inversely proportional to the time. Their results show a break in the Fourier transform of the time delay, which occurs at the Fourier frequency corresponding to the time when low-energy photons start to arrive before the high-energy ones. Such a temporal evolution of the delay is a consequence of the mechanisms generating the X-ray flare considered in Lewis et al. (2016). Conversely, the LIV delay is expected to be constant along the flare as it is a cumulative effect over the propagation of all photons from the cosmic source to the Earth. It is entirely determined by the distance and the photon energies emitted by the source.

To find out the temporal evolution of the delays induced by the SSC model presented here, a simpler method is applied. The evolution of the delay is evaluated by comparing the time difference Δt_evol(t) between two light curves reaching the same normalized flux value and by computing this time difference all along the flare. In other words, we noted the normalized fluxes of the two light curves F₁ and F₂, and considering two times t and t₂ such as F₁(t) = F₂(t₂), Δt_evol(t) = t − t₂. To compare this study with that of Lewis et al. (2016), we have to adopt their opposite sign convention for the time delays, meaning that in this section, a positive delay corresponds to low energies arriving after the high energies. The two light curves chosen for the comparison are integrated over the energy ranges 200−400 GeV and 2.6−4.7 TeV. The set of parameters considered here is the one presented in Table 1, corresponding to an acceleration-driven regime, but results can be obtained in the same way for the cooling-driven regime. The normalized light curves and the temporal evolution of the delay are shown in Fig. 13. The result obtained here is similar to that shown in Fig. 1 of Lewis et al. (2016). At large times (small Fourier frequencies), the delay is large because the high-energy light curve decays faster than the low-energy one. A break occurs at the specific Fourier frequency corresponding to the time when a given flux is reached earlier for the high-energy light curve than for the low-energy one. At small times (large Fourier frequencies), the delay is negative because the high-energy flare rises after the low-energy one due to the time needed for electrons to be accelerated. Such results are in agreement with the analysis by Lewis et al. (2016) and confirm that intrinsic time delays can be significantly variable along the flares and show characteristic time profiles. This variation of the intrinsic time delays explains the result obtained in Sect. 4.1, concerning the CCF which cannot correctly reconstruct the injected delay. This is directly due to the fact that the CCF uses the full light curves to evaluate time delays.

Fig. 13. Light curves for a flare with parameters of Table 1 with B = 65 mG for two energy bands (left) and evolution of the simulated time delays between them as a function of time (right). The LIV-like delay overlaid for comparison has been obtained by simulating two asymmetric Gaussian light curves with an injected constant delay of Δt = 500 s. Only positive delays are shown here in logarithmic scale, which can be directly compared to the Lewis et al. (2016) description.

To describe in the same way the situation when only a LIV time delay is present, two light curves were simulated, following the same asymmetric Gaussian shapes as the light curves coming from the SSC model. A constant time delay of Δt = 500 s is then injected between the two simulated light curves. The resulting temporal evolution of the delay is given in Fig. 13, which illustrates the constant LIV time delay obtained as a function of the Fourier frequency. Such results clearly show that temporal evolution of time delays is a direct marker for the presence of intrinsic effects in flares observed from blazars.

7. Conclusion

The blazar flare model considered in this paper describes the evolution of a population of relativistic electrons radiating in a single compact zone through their synchrotron-self-Compton emission. In order to explore basic intrinsic energy-dependent time delays expected in such sources, we focus on a minimal scenario taking into account only the dominant processes needed to generate flares, namely a generic acceleration mechanism and radiative cooling in a slowly varying magnetic field. Under reasonable assumptions, an analytical solution can outline the electron spectrum evolution and applies for instance when synchrotron losses dominate over inverse-Compton losses. Such a scenario clearly emphasizes the likely presence of significant intrinsic time delays. It reveals the existence of two main time-delay regimes, referred to here as acceleration-driven and cooling-driven regimes, over a large domain of parameters. The system evolves in one of these two regimes depending on the mechanism driving the evolution of the most energetic electrons at the time when the light curves reach their peak and the flare starts to decay. The cooling-driven regime typically corresponds to cases where the decay of the flare is mostly dominated by the radiative cooling effects. Conversely, the acceleration-driven regime corresponds to cases where the acceleration of emitting VHE electrons goes on while the flare starts to decay under the influence of some loss mechanisms other than radiative cooling, such as decrease of the ambient magnetic field or adiabatic losses. The confirmed detection of one of these regimes during blazar monitoring would provide precious information on detailed processes generating flares and significant measurements of time delays would further constrain source parameters in situ.

However, only upper limits on time delays have been firmly confirmed so far in TeV blazar flares. The detection of possible time delays in a flare of Mrk 501 has remained unique and was observed by only one instrument. Such a situation could be due to the lack of high-quality data on blazar flares. It could also directly put basic flare scenarios into question, which will require further investigation. In the most simple scenarios with fast initial injection or acceleration of particles immediately followed by flare decay due to radiative losses, a cooling-driven regime could have already been observed during bright flaring events since the cooling time over the VHE domain is longer than the time resolution reached. Radiative cooling alone should typically induce intrinsic time delays of several minutes. The fact that such lags have not yet been detected suggest that flares could mainly occur in a specific range of parameters, corresponding to an intermediate zone between the two time-delay regimes identified in this paper. As a consequence, there should be a physical link between acceleration and cooling processes with fine-tuning effects during the global evolution of the flare. A possible qualitative scenario would be to consider launching the flare by sudden shock acceleration or magnetic reconnection in the emitting zone, with subsequent mechanisms which induce the flare decay by adiabatic expansion and/or magnetic field decrease, while acceleration processes are still efficient enough at the highest energies. Future observations will be necessary to constrain scenarios and improve time-dependent modeling of blazar flares.

Moreover, the intrinsic delays obtained within the SSC scenario show specific characteristics which could help to constrain QG models or new physics involving time delays. The temporal dependency of the intrinsic delays appears to be different from the one expected by the current description adopted for LIV effects, which may provide a characteristic signature for the presence of intrinsic effects. Indeed, LIV delays are not expected to show any kind of evolution with time since they affect photons in the same way throughout their propagation. On the contrary, intrinsic delays evolve with time due to the different energy-dependent mechanisms involved in the generation of blazar flares. In addition, the energy dependency of intrinsic time delays at GeV–TeV energies was found to present typical power index α in the range [0.45−0.85]. This property can be explored in order to test specific QG models which predict energy-dependent LIV delays with an index different from the typical intrinsic one. Nevertheless, further study of QG models involving LIV effects is necessary to fully exploit the global time-delay information when it becomes available and to distinguish between the various effects. Another important feature to exploit is that LIV delays depend strongly on the propagation distance while intrinsic delays should not.

Briefly, tools and results presented in this paper contribute to the scientific preparation of the new gamma-ray instruments of the coming decade which should provide higher sensitivity and a much larger number of blazar flare detections than current IACTs. Future data will potentially lead to the detection of significant time delays. Flare scenarios should be further developed in order to explain the new observables on intrinsic time delays, and to help disentangle intrinsic and extrinsic effects, opening a way for time-of-flight LIV searches.

References

Appendix A: Table for high-energy time delay

All Tables

All Figures

Fig. 1. Electron spectrum (top) and SED (bottom) evolution for the reference set of parameters given in Table 1, illustrating a typical TeV blazar flare. The two plots share the same color code. The full lines correspond to the rising phase of the flare and the dashed lines to its decay. The time tmax is defined as the time when the highest energy value γmax = max[γc(t)] is reached in the electron spectrum. The corresponding electron spectrum and SED are plotted in thick lines. The electron spectrum and SED at the maximum luminosity of the flare are not displayed because they almost coincide with the ones at tmax for this reference set of parameters.

In the text

Fig. 2. Normalized light curves (left), zoom close to the light-curve maxima (center), and time-delays between each light curve and the one in the range 1–1.8 MeV (right) obtained for the first reference set of parameters. The vertical dashed line on the light curves (left) corresponds to the time when electrons reach their maximum energy γmax and shows the moment when energy losses dominate over the acceleration at the energy γmax.

In the text

Fig. 3. Normalized light curves (left), zoom close to the light curve maxima (center), and time-delays between each light curve and the one in the range 1–1.8 MeV (right) obtained for the second reference set of parameters (i.e., with B = 90 mG). The vertical dashed line on the light curves (left and middle) corresponds to the time when electrons reach their maximum energy γmax and shows the moment when energy losses dominate over the acceleration at the energy γmax.

In the text

	Fig. 4. Time delay measured with PPM and CCF between two simulated asymmetric Gaussian light curves varying the width of one of them. The injected delay value is fixed at Δtinj = 800 s. The measurement from the CCF shows a misreconstruction of the injected time delay as the width difference increases. The asymmetric Gaussian function used has σleft = 2500 s and σright = 20 000 s.
In the text

	Fig. 5. Time delay vs. energy for different B0 values. The two cases in bold in the legend correspond to the ones discussed in Sect. 4.2. All other parameters have the values given in Table 1.
In the text

	Fig. 6. Time delay vs. energy for different magnetic field evolution index mb with B0 = 65 mG (left) and B0 = 90 mG (right). All other parameters are unchanged (Table 1). The two cases in bold in the legend correspond to the situation discussed in Sect. 4.2.
In the text

	Fig. 7. Time delay vs. energy for different Doppler factor δ values with B0 = 65 mG (left) and B0 = 90 mG (right). All other parameters are unchanged (Table 1). The two cases in bold in the legend correspond to the situation discussed in Sect. 4.2.
In the text

	Fig. 8. Time delay vs. energy for different acceleration amplitude A0 values with B0 = 65 mG (left) and B0 = 90 mG (right). All other parameters are unchanged (Table 1). The two cases in bold in the legend correspond to the situation discussed in Sect. 4.2.
In the text

	Fig. 9. Time delay vs. energy for different acceleration evolution ma values with B0 = 65 mG (left) and B0 = 90 mG (right). All other parameters are unchanged (Table 1). The two cases in bold in the legend correspond to the situation discussed in Sect. 4.2.
In the text

	Fig. 10. Time delay vs. energy for different initial electron index n values with B0 = 65 mG (left) and B0 = 90 mG (right). All other parameters are unchanged (Table 1). The two cases in bold in the legend correspond to the situation discussed in Sect. 4.2.
In the text

	Fig. 11. Time delay vs. energy for different initial maximum energy for electron γc, 0 values with B0 = 65 mG (left) and B0 = 90 mG (right). All other parameters are unchanged (Table 1). The two cases in bold in the legend correspond to the situation discussed in Sect. 4.2.
In the text

Fig. 12. Typical time delays obtained from the SSC flare model (shaded area) and expected from a linear LIV effect with an energy scale at the Planck energy (full lines) with the redshift evolution found in Jacob & Piran (2008). Other formalism such as DSR would lead to different values but of the same order of magnitude. The dashed line corresponds to Δt = 0. Positive and increasing LIV delays correspond to subluminal LIV effects, negative and decreasing ones to superluminal LIV effects.

In the text

Fig. 13. Light curves for a flare with parameters of Table 1 with B = 65 mG for two energy bands (left) and evolution of the simulated time delays between them as a function of time (right). The LIV-like delay overlaid for comparison has been obtained by simulating two asymmetric Gaussian light curves with an injected constant delay of Δt = 500 s. Only positive delays are shown here in logarithmic scale, which can be directly compared to the Lewis et al. (2016) description.

In the text