Giovanni Amelino-Camelia

Arache Djannati-Atai

Giles Henri

Frédérique Piron D'

Andrew

Bruno

Dmitry David

François Fabian

Jean-François Floriana

Markus Matthias

Monica Nukri

Anyssa

Changqiao, Dilia, Émilie, François, Ilaria Joao Audrey

Simon Robert Louis

Julien Martine Guillaume

Zakaria Merci

les remarques qu'ils m'ont suggèré

à la soutenance. Je voudrais aussi remercier toutes les personnes de la collaboration H.

Introduction

The Quantum Gravity (QG) theory aims at describing Physics at the Planck scale (E P ∼ 10 19 GeV). This theory, still under construction, would allow for the description of systems that cannot be fully understood with current theories or by combining quantum physics and general relativity. Several tentative models exist trying to achieve such a universal theory, like string theory or loop quantum gravity. However, due to the nature of the QG theory and the extreme conditions it intends to describe, there is no earth based experiment able to constrain or reject these tentative models.

Nevertheless, these models can be indirectly tested using some of their properties. One of which is the possible violation of the Lorentz invariance principle. Indeed, some categories of models predict a Lorentz Invariance Violation (LIV) close to the Planck scale. Despite the fact that Lorentz invariance is found to be exact at low energy, faint manifestations of LIV are hoped to be observed at γ-ray energies with Imaging Atmospheric Cherenkov Telescopes such as the H.E.S.S. (High Energy Stereoscopic System) experiment. One of these manifestations is the alteration of the propagation of photons in vacuum, which can lead to energy dependent velocities for photons and thus delaying the arrival time of photons of different energies. With a sufficiently large energy range and a very long photon propagation distance, these delays are expected to be detectable, using variable astrophysical sources such as flaring Active Galactic Nuclei (AGN).

The most promising AGN to search for LIV signatures are blazars which are the ones showing the most energetic and variable emission. These astrophysical objects exhibit strong flaring activities leading to emission observed up to tens of TeV, with a variability down to the minute time scale. However, if an energy dependent time delay is found from the observation of a blazar flare, this delay will have to be clearly understood. Indeed, another type of delay can be induced in situ during generation of a flare at the level of emission mechanisms. Such intrinsic effects will then come in addition to the LIV propagation effect. As a consequence, intrinsic time delays from blazars need to be studied in the energy range where LIV signatures are searched for, in order to evaluate their contribution and hopefully help distinguishing them from LIV time delays. This thesis attempts, for the first time, to study the search of LIV signatures with the modeling of blazar flares in order to investigate energy dependent time delays. On the one hand, a time dependent model of blazar flares, focusing on the emission at γ-ray energies, is developed in order to characterize intrinsic time delays. This study aims at understanding the origin of intrinsic delays and at finding a way to differentiate them from the LIV ones. On the other hand, a search for energy dependent time delays is carried out, for a flare from the Markarian 501 blazar which occurred in the night of June 23-24 2014. This flare, observed by H.E.S.S., shows very high energy photons, up to 20 TeV, and a variability of a few minutes making it a good candidate for the search of LIV signatures. This manuscript begins with a short introduction on γ-ray astronomy with the H.E.S.S. experiment in Chapter 1, Quantum Gravity and LIV in Chapter 2 and the non-thermal emission and acceleration mechanisms at work in blazars in Chapter 3. The first part of the thesis work describes the time dependent model developed to study intrinsic time delay in Chapter 4. This model is based on a temporal evolution of the particle population responsible for the γ-ray emission. Then, the delays produced by the model are fully characterized by evaluating the influence of each parameter on them and the results from the model are discussed to help the search of LIV signatures in Chapter 5. The second part of the thesis work presents the main procedure used to measure the LIV energy dependent time delay, based on the maximum likelihood method which is applied to Monte-Carlo simulations in Chapter 6. These simulations allow the introduction of a new correction term for the method to deal with very high energy data sets like the flare of Markarian 501. This flare is then analyzed and used to search for a LIV signature in order to provide constraints on the energy scale of LIV in Chapter 7. Finally, the conclusion puts in perspective the modeling results with the search of LIV signature, as well as future developments to establish more robust constraints on QG models and blazar flare scenarios.

Chapter 1

Gamma-ray astronomy with H.E.S.S. Gamma-ray astronomy aims at observing photons of astrophysical origin and with energies ranging from hundreds of keV up to several TeV. These photons can be detected by two main processes related to their interaction with matter. Below 10 MeV, γ rays mainly interact via Compton scattering. Detection techniques in this regime have a major drawback: a poor accuracy on the reconstruction of the γ-ray photon direction. As a consequence, γ-ray astronomy at such energies is very difficult. For γ rays above 10 MeV, the dominant process when they interact with matter is electron-positron pair creation, which triggers the development of an electromagnetic shower. The particles created during the shower development can be used in order to reconstruct the initial γ ray properties.

The direct measurement of secondary particles requires the use of space telescope such as Fermi [START_REF] Atwood | The Large Area Telescope on the Fermi Gamma-Ray Space Telescope Mission[END_REF]. Indeed, the electromagnetic shower has to develop entirely inside the detector or secondary particles can be lost, preventing an accurate determination of the initial γ-ray photon properties. However, due to the universal nature of the power law photon energy spectrum from γ-ray sources, γ rays are rarer as their energy increases which prevents satellite experiments to be sensitive at energies higher than ∼ 300 GeV due to their small collection area. To measure a significant amount of γ-ray photons at higher energies, ground-based telescopes can be used to perform observations of electromagnetic showers. When a shower develops in the atmosphere, secondary particles emits Cherenkov radiation because their velocities are higher than the speed of light in the air. This Cherenkov radiation can thus be used to image the electromagnetic shower and deduce the properties of the primary γ ray.

In this chapter, the techniques used to detect γ rays with ground based telescopes are introduced. Then, the H.E.S.S. (High Energy Stereoscopic System) experiment is described. Finally, the full analysis chain is briefly presented, from the event reconstruction to the spectrum and light curve determination.

Detection of gamma rays with ground based telescopes 1.Extensive air shower

When a γ ray penetrates the Earth atmosphere, it quickly interacts with the Coulomb field of atmospheric nuclei leading to a pair creation process: γ → e + e -. Then, each of the two created particles interacts with the same Coulomb field and emits a deceleration radiation called Bremsstrahlung emission, producing several γ-ray photons during their propagation in the atmosphere. These secondary γ rays also interact with atmospheric nuclei leading to pair creation while the electrons (positrons) continue emitting Bremsstrahlung until these particles reach a critical energy value E C 84 MeV, which corresponds to the energy when ionization losses become more important than the ones induced by Bremsstrahlung. The propagation of electrons and in the atmosphere is represented by the atmospheric depth X (in g cm -2 ) which is the quantity of atmosphere crossed depending on its density. With such a physical quantity, one describes the Bremsstrahlung emission using its typical radiation length X 0 , corresponding to the average atmospheric depth crossed by an electron before emitting a Bremsstrahlung photon, and losing a fraction 1/e of its energy:

1 X 0 = 4dα A (r e Z) 2 log 183 Z 1/3 , (1.1)
with α the fine-structure constant, r e the electron radius, and Z, A and d respectively the atomic number, mass and density of the medium. For air, X 0 ≈ 36.7g cm -2 while the total depth of atmosphere is about 27 X 0 .

Detection of gamma rays with ground based telescopes

All the particles created from the primary γ ray constitute the electromagnetic shower. The charged particles in the shower are produced with velocities higher than the speed of light in air, defined as c air = c/n air with n air the refractive index of air. As a consequence, these particles produce a Cherenkov radiation near the ultraviolet (λ ∼ 300 -600 nm). All the Cherenkov photons from the shower are emitted extremely fast, from 5 to 20 ns. The angle θ at which Cherenkov photons are emitted depends on the refractive index n of the medium and on the particle velocity. For particles with a velocity v ∼ c, it is expressed as: .2) In this atmosphere, with a refractive index very close to 1, a shower with an altitude of 10 km emits a Cherenkov light pool illuminating a surface of about 10 5 m 2 at the ground level. This allows for ground based telescope with a large instrumented area to obtain an extremely large collection area in order to deal with the low flux of γ-ray photons above hundreds of GeV.

cos θ = c n × v ∼ 1 n . ( 1
Unfortunately, not only γ rays produce atmospheric showers when arriving on Earth. For one γ-ray photon, there is ∼ 10 5 hadronic particles (proton and nuclei). These hadronic particles trigger the development of hadronic showers when interacting with an atmospheric nuclei. The result of these interactions leads mainly to the creation of charged and neutral pions. Neutral pions immediately decay to produce a pair of γ-ray photons. In opposition, charged pions have a longer life time, allowing them to interact and produce more pions, until their reach the critical energy E π C ∼ 20 GeV when they decay and produce muons and neutrinos. These hadronic showers are a background for γ-ray detection and have to be discriminated, as detailed in Section 1.3.

Heitler models

A simple model to represent the electromagnetic shower development is the Heitler model [START_REF] Heitler | Quantum theory of radiation[END_REF]. It allows for a description of all the particles created during the shower development with only two parameters: the initial energy E 0 of the γ ray and the radiation length X 0 between each interaction. For this model, the radiation length of pair creation, λ γ is supposed to be the same as the one of Bremsstrahlung (actually λ γ = 9/7X 0 ). It allows for a deduction of the particle energy at an atmospheric depth X:

E(X) = E 0 exp - X X 0 . (1.
3)

The model describes the electromagnetic shower using a step R = X 0 ln(2) corresponding to the atmospheric depth that a particle has to cross in order to lose half of its energy. At each step R, a particle in the shower interacts either via Bremsstrahlung emission for electrons or pair creation for γ rays during which they loose half of their energy, split into the two created particles. The scheme of shower development from the description of the Heitler model is shown in Figure 1.1a. The shower ends when the particles at the last stage n C reach the energy E C corresponding to the domination of ionization losses for electrons and the beginning of the Compton regime for γ rays.

Hadronic showers can also be described with an extended Heitler model [START_REF] Matthews | A Heitler model of extensive air showers[END_REF]. The interaction of the initial hadron leads to the creation of several pions. About one third of the pions are neutral while the rest are charged. The characteristic radiation length for pions is X π 0 = 120 g cm -2 . After one X π 0 , π 0 decays into two γ rays which trigger the development of an electromagnetic shower and π ± interact with atmosphere leading to more pions (∼ 1/3π 0 and ∼ 2/3π ± ). The end of the shower is reached when charged pions reach their critical energy E π C and decay into muons and when the secondary γ-ray photons from neutral pions end their electromagnetic showers. The scheme of this hadronic shower model is shown in Figure 1.1b. Hadronic showers lead to wider and asymmetric atmospheric shower than electromagnetic ones and so is their Cherenkov emission. Such a property can be used to distinguish them in order to reject them for the data analysis. As a remark, such a simple model is not sufficient to fully describe hadronic showers and more complex models are needed such as numerical models. These numerical models use precise descriptions of the different interaction during the shower development, based on data collected by particle accelerators like the Large Hadron Collider (LHC). Some of the main models used to describe hadronic shower development are [START_REF] Fletcher | Extragalactic optical-infrared background radiation, its time evolution and the cosmic photon-photon opacity[END_REF], [START_REF] Kalmykov | Quark-Gluon-String Model and EAS Simulation Problems at Ultra-High Energies[END_REF], Pierog et al. (2015). (right). For hadronic showers, π 0 decay creating two γ rays leading to electromagnetic showers described by the left shower. R is the radiation length after which a particle loses half its energy.

Detection techniques

To assess the initial particle properties (nature, energy and direction), the Cherenkov light emitted by all the charged particles composing the atmospheric shower is measured and used to determine the properties of the shower which then can be compared to a shower model such as Heitler models. Actually, there are two main techniques using the Cherenkov light of the shower: the Cherenkov imaging technique and the air-shower particle sampling technique.

Imaging Atmospheric Cherenkov Telescopes (IACT) use the imaging technique with Cherenkov light which is the most used detection technique at the moment. Measuring the Cherenkov image of an atmospheric shower allows for the determination of its size, intensity and orientation which is related to the energy and the direction of the initial particle. This technique was widely used and allowed for the development of TeV γ-ray astronomy with experiments such as Whipple [START_REF] Cawley | A high resolution imaging detector for TeV gamma-ray astronomy[END_REF], HEGRA [START_REF] Daum | A high performance likelihood reconstruction of γ-rays for imaging atmospheric Cherenkov telescopes[END_REF] and CAT [START_REF] Barrau | The CAT imaging telescope for very-high-energy gamma-ray astronomy[END_REF]. While CAT (Cherenkov Array at Themis) and Whipple had a single telescope, HEGRA (High Energy Gamma-Ray Astronomy) and now the current generation of IACT use multiple telescopes to have multiple Cherenkov images of the showers from different direction to better characterize the primary particle using stereoscopic mode. The instruments currently in operation are:

• H.E.S.S. : The High Energy Stereoscopic System, is an hybrid array composed of five telescopes located in Namibia [START_REF] Aharonian | Observations of the Crab nebula with HESS[END_REF],

• MAGIC: The Major Atmospheric Gamma-ray Imaging Cherenkov is an array composed of two telescopes located in Canary Islands [START_REF] Aleksić | The major upgrade of the MAGIC telescopes, Part II: A performance study using observations of the Crab Nebula[END_REF],

• VERITAS: The Very Energetic Radiation Imaging Telescope Array System is an array composed of four telescopes located in Arizona (USA) (Park & VERITAS Collaboration, 2015),

• FACT : The First Geiger-mode Avalanche photodiode Cherenkov Telescope, is an automatized telescope located in Canary Islands [START_REF] Biland | Calibration and performance of the photon sensor response of FACTthe first G-APD Cherenkov telescope[END_REF].

These experiments use large mirrors to collect the Cherenkov light and to focus it on a camera with a fast electronic to be able to record such short Cherenkov light flashes. They have a relatively small Field of View (FoV), of the order of a few degrees and observe during moonless night time for optimal efficiency. MAGIC [START_REF] Ahnen | Variable Very High Energy γ-Ray Emission from Markarian 501[END_REF], VERITAS [START_REF] Archambault | Gamma-ray observations under bright moonlight with VERI-TAS[END_REF] and FACT [START_REF] Knoetig | FACT -Long-term stability and observations during strong Moon light[END_REF] also observed during moon time but with degraded performances.

The particle sampling technique consists in measuring directly the secondary particles at a given stage of the atmospheric shower development using water tanks. The use of water is justified by a larger refractive index of n w = 1.333, inducing a higher Cherenkov light yield. However, this technique requires the instrument to be located at high altitude to record airshowers during their development. As examples, the experiment Milagro [START_REF] Atkins | TeV Gamma-Ray Survey of the Northern Hemisphere Sky Using the Milagro Observatory[END_REF] and its successor HAWC (Abeysekara et al., 2017) use this air-shower particle sampling technique. Such instruments measure less events in average, due to a higher energy threshold and a smaller collection area but they are able to observe at any time. Since photo-detectors are inside a water tank, they are not affected by sunlight or moonlight allowing the instrument to monitor a very large fraction of the sky. Such experiment are therefore quite complementary to IACT.

H.E.S.S. instrument

The H.E.S.S. experiment is located in Namibia at a hundred kilometers south of the capital Windhoek in the Khomas region and at an altitude of 1800 m above sea level. The site was chosen because of the good atmosphere quality (for an optimal shower development and to reduce the Cherenkov light absorption) and of the desert climate of the region with small cloud coverage and the absence of light pollution. In addition, H.E.S.S. is the only IACT located in the southern hemisphere allowing for a good coverage of the galactic center and the Magellanic clouds. The experiment observes only during night time when the moon is not present in the sky, which represents 26 nights of observations every lunar month (29 nights). The effective area of the instrument is ∼ 50000 m 2 averaging over energy and under optimal observation conditions.

H.E.S.S. telescopes

The instrument is an hybrid array of Cherenkov telescopes composed of four small telescopes and a larger one (Figure 1.2). The array is a square of 120 m side length with one small telescope at each corner and the big one in the center. The diagonals of the array are oriented North-South and East-West. All the telescopes have an alt-azimuth mount to be able to point any position in the sky.

Small telescopes (CT1-4) were designed for the first H.E.S.S. phase (HESS-I). They have been operational since December 2003 and work only in stereoscopic mode. They have a 12 m diameter segmented mirror. This mirror, composed of 382 round mirrors of 60 cm each, is a spherical mirror with a surface of 108 m 2 . They are able to reconstruct γ-ray photons with energies from ∼ 100 GeV up to tens of TeV.

The big telescope (CT5) has started observing since July 2012, initiating the second H.E.S.S. phase (HESS-II). This telescope has a much larger segmented mirror of 28 m long, composed of 876 small mirrors of 90 cm. This represents a surface of 614.5 m 2 . In opposition to the small telescopes, this large mirror has a parabolic shape, in order to reduce the temporal dispersion of Cherenkov photons on the camera. Having a larger surface to collect Cherenkov light allows this telescope to detect γ-ray photons at lower energies down to 20 GeV (H.E.S.S. Collaboration et al., 2018). Indeed, a γ ray with such an energy induces a very faint Cherenkov light pool which then requires such a large mirror to be collect enough Cherenkov light for the detection. CT5 can work together with CT1-4 in stereoscopic mode or operate alone in monoscopic mode.

H.E.S.S. cameras

The cameras of CT1-4 (Aharonian et al., 2004a) are composed of 960 photo-multipliers (PM). They are arranged in 60 drawers, which are a group of 16 PM sharing the same front-end electronics (see one example on Figure 1.3 with a CT5 drawer). In order to ensure an optimal Cherenkov light collection, each PM has a Winston cone which guides the photons to the PM entrance. Winston cones also reduce the non-instrumented area between PM. Each PM has a FoV of 0.16 • in the sky, with a total FoV of 5 • for the whole camera. With the whole electronic system, these cameras were initially designed with a dead time, i.e. the time needed for the electronic to read an event, of about 460 µs. A recent upgrade of the camera electronics was done to reduce this dead time to 7 µs [START_REF] Klepser | Hardware and software architecture of the upgraded H.E.S.S. cameras[END_REF].

The CT5 camera (Bolmont et al., 2014) is composed of 2048 PM arranged in 128 drawers. Each pixel has a FoV of 0.07 • in the sky, for a total FoV of 3.2 • for the whole camera. The total dead time of the camera is about 20 µs. For an average trigger rate of 1 kHz, the dead time represents 1% of the total acquisition time of the camera.

Trigger system

In order for the whole array to trigger and record an event, several conditions must be fulfilled, based on individual telescopes and the whole array. For each individual telescope, there are two levels of trigger. The first level is defined at the pixel level and requires for one pixel to exceed a collected charge of four photo-electrons (p.e.). The second level requires at least three pixels passing the first level in a group of four neighbor drawers (defined as a "sector") within a temporal window of ∼ 1.3 ns. If the second level is passed, a trigger decision is sent to the central trigger of the array. If a minimum of two telescopes send a trigger decision, within a coincidence window of 80 ns, the event is recorded with the full array. The use of such a stereoscopic mode allows for the rejection of most of the Night Sky Background (NSB) noise which is supposed to produce a random signal in each camera.

CT5 can also record events on its own in the monoscopic mode. It is mostly the case for events with an energy below 100 GeV for which CT1-4 are not sensitive enough. If the central trigger receives a CT5 signal, it automatically records the event. This leads for CT5 to a larger fraction of events recorded being background noise than for the whole array.

Data calibration

In order to precisely reconstruct γ-ray photon, there are two main calibration parameters to determine. The first one is the conversion coefficient between the charge (in p.e.) and the numerical signal produced by an Analog-to-Digital Converter (ADC). For one PM, the amplitudes measured in two channels with a high (HG) and a low gain (LG) amplifier are expressed as:

A HG = F F ADC HG -P HG γ HG , ( 1.4
) .5) where A is the amplitude in p.e. in each channel, ADC is the amplitude in ADC counts in each channel, P is the mean value of the electronic noise (pedestal) in ADC counts in each channel, γ HG is the nominal ADC count in the HG channel, HG/LG is the amplification gain ratio between the two channels and F F is the Flat Field coefficient which corresponds to the response of the camera to a uniform field. Each of these parameters can be determined with specific observation runs.

A LG = F F × HG LG ADC LG -P LG γ HG , ( 1
To determine P HG and P LG , a so called pedestal run is done. It consists of recording the pedestal of the camera, in absence of any kind of signal.

The F F coefficient is computed from Flat Field runs. A LED (Light-Emitting Diode) is used to uniformly illuminate the camera and the intensity is recorded on each pixel. Then, the Flat Field coefficient can be deduced as:

F F = I cam I P M (1.6)
where I cam is the average intensity by the whole camera and I P M by the one of the PM. F F takes into account the variations of the quantum efficiency of a PM as well as the efficiency of its Winston cone Finally, γ HG is determined from single photo-electron runs. Using a pulsing LED, flashes are sent to the camera with an intensity adjusted to obtain on average one p.e. at the photo-cathode of each PM. The camera is triggered by a pulse synchronized with the light flash. Figure 1.4 shows the results of such a run for one pixel. The first peak on the figure corresponds to the pedestal and the second peak corresponds to the pixel response to a single photon electron. The γ HG value corresponds to the distance in ADC counts between the two peaks. The second calibration parameter is the optical transmission from the top of the atmosphere down to the cameras. This includes the atmosphere absorption, the collection efficiency, the optical efficiency of the mirror and of the Winston cones and the quantum efficiency of each PM. To compute it, the Cherenkov emission from muons (which are emitted during hadronic showers development) is used. The Cherenkov image of muons is unique and gives a ring which is easily recognizable. Muon rings are compared to analytical models of muon Cherenkov emission (Figure 1.5) in order to deduce the optical efficiency of the instrument based on the measured intensity [START_REF] Chalme-Calvet | Muon efficiency of the H.E.S.S. telescope[END_REF].

Event reconstruction

From the shower images recorded by the cameras, as the examples in Figure 1.6, one needs to determine if an event is a γ-like or a hadron-like event. Then, after selecting all γ-like events, the shape of the image is used in order to determine the size and the orientation of the shower, to finally deduce the energy and the direction of the γ ray.

The Hillas model

This model was developed by Hillas (1985) in order to discriminate between electromagnetic and hadronic showers using the image shape. A shower induced by a γ-ray photon produces a Cherenkov light pool giving elliptical shapes in each camera (Figure 1.6a). In opposition, the development of hadronic showers contains a lot of secondary particles which can themselves decay and induce secondary electromagnetic showers and muons, leading to a wider and irregular Cherenkov image (Figure 1.6b). The Hillas model uses the elliptical shape property of electromagnetic showers to identify γ-like events with several parameters (the Hillas parameters):

• the ellipse half length L These parameters are determined from images recorded by the telescopes after a cleaning procedure to remove most of the NSB. The direction of the shower image is defined as the prolongation of the major axis. Using more than one telescope allows for a better reconstruction of the source position. Then the angle θ is defined as the error on the source position in the sky. Using the Hillas parameters L and W , two discriminant variables are defined for each telescope which are the Scaled Width SW and the Scaled Length SL:

SW = W -W σ W ; SL = L -L σ L , ( 1.7)
where W and L are the mean ellipse half length and half width in all telescopes and σ W and σ L are their variance. Then, each SW and SL values of the different telescopes are combined to respectively produce the Mean Scale Width M SW and Mean Scale Length M SL defined for the whole array:

M SW = 1 √ N tel N tel i=1 SW i ; M SL = 1 √ N tel N tel i=1 SL i , (1.8)
where N tel is the number of telescopes having a shower image recorded. The distributions of these variables do not follow the same behavior for hadron-like and γ-like events and can be used to reject a part of the hadron-like event of the data. An example of the distribution of M SW for simulated γ-like and hadron-like showers as well as for an observation with no source in the FoV is shown on Figure 1.8. While the simulated γ-like events present a sharp distribution, the simulated hadron-like events have a wider distribution. Thus, some selection cut can be applied to the variable M SW and M SL in order to keep the majority of the γ-like events while rejecting most of the hadron-like event. However, some hadron-like events still remain after selection cut but are removed during the analysis when performing a background subtraction. To determine the energy of the initial γ ray, the amplitudes of the images from the different telescopes are used. For each telescopes, simulated shower images are produced with different initial energies, following the Hillas parameters obtained from the shower image. Then, the comparison between simulated and data images allows for the determination of the energies E i of an event for each telescope. The values are then averaged to produce the reconstructed event energy E rec from the array:

MRSW

E rec = 1 N tel N tel i=1 E i .
(1.9)

The semi-analytical model reconstruction

This event reconstruction method, called Model, is based on a semi-analytical model (de Naurois & Rolland, 2009) which consists of comparing directly the images recorded by the cameras to simulated images. These simulated images are produced by computing the Cherenkov emission from an analytical model of shower development, taking into account the efficiency of the instrument to measure the Cherenkov light. The distribution of Cherenkov light from the simulated showers is derived from the longitudinal, lateral and angular distributions of particles inside the shower. These distributions are parameterized from showers produced by the analytical development model. This produces a semi-analytical model for the Cherenkov emission coming from the shower induced by a primary γ-ray photon.

Then, images from data and simulations are compared using a maximum likelihood (ML) method at the pixel level. Using ML has several benefits such as taking into account the NSB level on each pixel and directly comparing raw images without any cleaning. The log-likelihood function used to compare the data image and the simulated one is defined as: (1.10) where P i (s i |µ i , σ p , σ γ ) is the probability density function to observe in a pixel i a signal of intensity s i where a value µ i is expected. σ p is the pedestal width, σ γ the photo-electron pulse width and N pixels is the total number of pixels participating to the image. Different simulated images are compared to the data depending on several parameters such as the γ ray energy, its direction, the impact distance from the telescope and the altitude of the first interaction (or primary depth). The log-likelihood is maximized when the simulated image matches the data one and thus the parameters used to simulate the image are attributed to the data event.

log(L) = -2 N pixels i=1 log P i s i | µ i , σ p , σ γ ,
From this method, several discriminant variables are defined for γ-hadron separation based on the log-likelihood function. They are the Goodness of the fit G, the Shower Goodness SG and the Background Goodness BG:

G = 1 2(N pixel -6) N pixels i=1 log L(s i | µ i ) -log(L) | µ i (1.11) SG = 1 √ N in N in i=1 log L(s i | µ i ) -log(L) | µ i (1.12) BG = 1 √ N out Nout i=1 log L(s i | µ i ) -log(L) | µ i (1.13)
where L(s i | µ i ) is the log-likelihood value for each pixel, log(L) | µ i is the expected log-likelihood value at µ i , N in is the number of pixels associated to the shower image and N out the other pixels (N pixels = N in + N out ). G computes how likely the best simulated image obtained with the semi-analytical model matches the image recorded from the data. SG computes how likely the pixel participating to the shower image intensity in the camera match between the simulated and the data image. Finally BG compare how likely the rest of the pixels match between the two images. For each event, each telescope which has a shower image recorded provides the values G i , SG i and BG i and allows for the computation of the Mean Scale Goodness M SG, the Mean Shower Scale Goodness M SSG and the Mean Scale Background Goodness M SBG:

M SG = 1 √ N tel N tel i=1 G i -G σ G , (1.14) M SSG = 1 √ N tel N tel i=1 SG i -SG σ SG , (1.15) M SBG = 1 √ N tel N tel i=1 BG i -BG σ BG . (1.16)
The Model analysis presents better performances than the Hillas model (de Naurois & Rolland, 2009) and is currently in use for the H.E.S.S. data analysis.

Data analysis 1.4.1 Analysis configuration

To analyze a data set, several analysis cuts have to be set as well as the desired telescope configuration. In H.E.S.S. there are four types of configurations, or profiles, available for the analysis:

• The Stereo profile is the most standard profile for the analysis and uses all the telescopes in Stereoscopic mode to reconstruct an event.

• The Mono profile is specifically designed for CT5. It only considers the events recorded by the large telescope alone. This profile aims at analyzing low energy data set.

• The Stereo HESS1 profile is a variant of the Stereo profile which only uses the small telescopes. This profile is mostly designed for observations prior the installation of CT5 or if CT5 is not available. This profile focuses more on high energy events.

• Combined3 : The Combined3 uses all the previous profiles to provide the best reconstruction. Each event is reconstructed with Stereo, Mono and Stereo HESS1 profile, and then the best reconstruction is chosen to determine the energy and direction of the event, based on the reconstruction errors [START_REF] Holler | Photon Reconstruction for H.E.S.S. Using a Semi-Analytical Shower Model[END_REF].

For all these profiles, different cut configurations are defined depending on the source of interest. These cuts are tuned to reject more or less hadron-like events which also impact the number of γ-like events passing the analysis cuts. The most to the less selective cuts configuration are called Safe, Standard (Std), Loose and VeryLoose. As an example, two sets of analysis cuts are presented (Combined3 Std and Stereo VeryLoose) in Table 1.1 which are going to be used in Chapter 7. These cuts concerns the camera images with a minimum charge value in p.e. recorded, a minimum number of pixels in the image and a distance between the shower and the camera center. Some other cuts are made during the event reconstruction on M SSG, the primary depth, the level of NSB and the error on the reconstructed direction of the γ ray. The last cut on θ 2 is to keep only events coming from a region centered around the source called the ON region.

Acceptance and energy resolution determination

With a specific analysis configuration, the response of the instrument can be determined. The Instrument Response Functions (IRF) are composed of the acceptance of the instrument, the energy resolution and the energy bias. They are all computed using Monte-Carlo simulations of γ-like events, using the program KASCADE [START_REF] Kertzman | Computer simulation methods for investigating the detection characteristics of TeV air Cherenkov telescopes[END_REF] for the numerical simulation of atmospheric showers and a simulation of the detector measurement. Then, these simulations, with known parameters, are analyzed and their results allow for the determination of the IRF.

The acceptance is computed from the ratio between the number of events passing the analysis cuts and the total number of simulated events. The acceptance depends on the energy of the events as well as on the observation conditions such as the Zenith Angle (ZA), the position of the shower image in the camera and the atmospheric conditions. From the acceptance, the energy threshold of the analysis can be defined as the energy corresponding to 15% of the maximum acceptance value is reached. The energy threshold corresponds to the minimum energy when the analysis is considered as optimal. The energy resolution and bias are determined by looking at the distribution reconstructed energy for a single value of true energy. This distribution follows in average a Gaussian distribution in logarithmic scale. The width of this distribution is taken as the energy resolution and the shift between the true energy and the mean of the Gaussian distribution is the energy bias.

Background subtraction

A last step for the data analysis is to remove the contribution of the γ-ray background as well as the hadron-like events which passed the selection cuts. To estimate this background, OFF regions outside the ON region are used. From the events reconstructed in these regions, the number of excess events N ex in the ON region (with respect to the OFF regions) is defined as: .17) where N ON and N OFF are the number of events in the ON and OFF regions, and β is a normalization factor corresponding to the ratio of the acceptance between the two regions. It is expressed as:

N ex = N ON -βN OFF , ( 1
β = ON A ON (Φ x Φ y Φ z )dΦ x dΦ y dΦ z dEdt OFF A OFF (Φ x Φ y Φ z )dΦ x dΦ y dΦ z dEdt (1.18)
where A ON and A OFF are the acceptance in the ON and OFF regions, Φ x and Φ y correspond to the event position in the FoV and Φ z is the zenith angle. The main procedure to define the OFF region for the H.E.S.S. analysis is the Reflected Background method. It consists on defining multiple OFF regions of the same size than the ON region. ON region and OFF regions are all defined with the same offset with respect to the camera center. The main advantage of this method is the computation of the normalization β for the computation of N ex . Indeed, as the acceptance has a radial symmetry around the camera center, the disposition of the ON and OFF regions as well as their size guarantee that they have the same acceptance. β then simply becomes the inverse of the number of OFF regions. Besides, in case of the presence of another source in the FoV some OFF regions can be excluded if they are on this source for the determination of the background contribution.

Combined3 Std Stereo VeryLoose

(deg) [0, 2] [0, 1.4] [0, 2] M SSG [-4, 0.9] [-4, 0.6] [-4, 0.9] PrimaryDepth (X 0 ) [-1.1, 3.4] [-1.1, 3.4] [-1.1, 3.4] NSBGoodness > 32 >32 - N SBG > 32 >32 - Direction Error [0, 0.2] [0, 0.3] [0, 0.2] θ 2 cut 0.0015 • 0.0015 • 0.01 •

Spectral analysis

In order to deduce the source energy spectrum and light curve, methods based on loglikelihood are used from the event list taking into account the instrument IRF.

Forward folding technique

The forward folding technique (Piron, 2000) requires an assumption on the source spectral shape Φ(E true , -→ Λ ) with some parameters -→ Λ and with E true being the true energy of the incoming γ ray. The most common function for the spectral shape are:

• Power-law spectrum, with Φ 0 being the energy spectrum normalization at the energy E 0 and Γ the spectral index:

Φ(E true | Γ, E 0 ) = Φ 0 E true E 0 -Γ , ( 1.19)
• Curved power-law spectrum, with β being the curvature index

Φ(E true | Γ, β, E 0 ) = Φ 0 E true E 0 -Γ-β log(Etrue/E 0 ) , ( 1.20)
• Broken power-law spectrum, with Γ 1 and Γ 2 the two spectral indices and E cut the break energy:

Φ(E true | Γ 1 , Γ 2 , E cut , E 0 ) = Φ 0      Etrue E 0 Γ 1 forE true ≤ E cut Ecut E 0 Γ 1 -Γ 2 Etrue E 0 -Γ 2 forE true ≥ E cut (1.21)
• Power-law spectrum with exponential cut off, E cut being the cut off energy:

Φ(E true | Γ, E cut , E 0 ) = Φ 0 E true E 0 -Γ exp E true E cut . (1.22)
After selecting a spectral shape, a log-likelihood comparison is done between the data events and the expected number of events computed from the spectral shape with a given set of parameters. The formula to compute the expected number of γ-like events n γ in the energy bin [E rec,1 ,E rec,2 ] taking into account the instrument IRF is: (1.23) where A(E true | Φ x Φ y Φ z ) is the acceptance, and R(E rec , E true | Φ x Φ y Φ z ) is the energy resolution, which depends on the observation conditions and on the energy of the event.

n γ = E rec,2 E rec,1 ∞ 0 R(E rec , E true | Φ x Φ y Φ z )A(E true | Φ x Φ y Φ z ) × Φ(E true | - → Λ )dE true dE rec ,

Spectral analysis

The probability to detect N ON and N OF F events, while n γ and n bck background events are expected from the spectral shape, is given in Poisson statistics by:

P (N ON , N OF F | n γ , n bck ) = n γ + τ n N ON bck N ON ! e -nγ + τ n bck × n N OF F bck N OF F ! e -n bck , (1.24)
where τ is the live time normalization, i.e. the ratio between the observing time in the ON region and the sum of the observing times of all the OFF regions. The expected number of background events is evaluated by maximizing the log-likelihood function L = log(P ) over n bck which gives the results:

n bck = C + ∆ 2τ (τ + 1)
, (1.25)

C = τ (N ON + N OF F ) -(1 + τ )n γ , (1.26) ∆ 2 = C 2 + 4τ (τ + 1) × N OF F × n γ . (1.27)
Finally, the spectral parameters are evaluated by maximizing again the log-likelihood function over the -→ Λ parameters to provide the best fit parameters -→ Λ 0 with their uncertainties as well as the covariance matrix (giving the correlation between the spectral parameters), the expected number of events in each spectral bin and the goodness of the fit defined as the log-likelihood function value at

- → Λ = - → Λ 0 .
The expected number of events in each bin obtained from the ML is then compared to the data to deduce the residuals.

Light curve determination

The standard determination of light curve in the H.E.S.S. spectral analysis relies on a similar approach as the forward folding method. It computes the evolution of the integrated flux, above the energy E int . For each temporal bin, n γ and n bck above E int are obtained from the reconstructed spectrum and compared to the number of events observed N ON and N OF F above the same energy value. However, this method does not converge when the number of events in either the ON or OFF regions is zero.

Another approach is the Simple Statistics method which allows for a determination of the flux if either N ON or N OF F is zero. It's a simple method using the excess number of events N ex in each temporal bin. The errors are determined using the Li&Ma statistics [START_REF] Li | Constraining Lorentz Invariance Violation Using the Crab Pulsar Emission Observed up to TeV Energies by MAGIC[END_REF]. It is expressed for each bin as: .28) where N γ is the expected number of events computed from the spectrum, and Φ(E, -→ Λ ) is the energy spectrum.

Φ(t, E > E int ) = N ON -βN OF F N γ ∞ E int Φ(E, - → Λ )dE. ( 1

Summary

In this chapter, the development of atmospheric showers induced by γ-ray photons entering the atmosphere was briefly presented, as well as the different techniques used to detect them from the Cherenkov photons they produced while crossing the atmosphere. Then, the H.E.S.S. experiment was described. The Model analysis, with the description of the log-likelihood comparison pixel by pixel from a semi-analytical shower model was then detailed. This method allows for the discrimination between γ-like and hadron-like events and also the determination of the event energies and directions. Finally, the forward folding technique was presented which determines the source spectra from the photon list and allows for the determination of the light curves.

The H.E.S.S. experiment allows for the investigation of various astrophysical sources such as active galactic nuclei or supernovae remnants but can also probe fundamental physics. As an example, the measurement of γ-ray photon from distant sources can provide tests for Lorentz invariance in the framework of Quantum Gravity models.

Chapter 2

Quantum gravity and Lorentz invariance violation Contents

At the Planck scale, Quantum Gravity (QG) is hoped to described the four fundamental interactions in the same unified theory. Currently, on one hand, electromagnetic, weak and strong interactions are described together in the Standard Model (SM) using Quantum Field Theory (QFT). On the other hand, Gravity is described alone by the General Relativity (GR) theory. The scale of this unified theory implies physical quantities using both gravity and quantum mechanic constants. This scale is given by the Planck length l P , the Planck mass m P or the Planck energy E P :

l P = hG c 2 ≈ 1.6 × 10 -35 m, (2.1
)

m P = hc G ≈ 2.18 × 10 -8 kg, (2.
2)

E P = m P c 2 ≈ 1.22 × 10 19 GeV. (2.3)
with h the reduced Planck constant, G the universal gravitational constant and c the speed of light. E P is believed to be the scale where all coupling constants of fundamental interactions converge at the same value. Actually, [START_REF] Wheeler | Geons[END_REF] has shown that the Planck length is the quantum limit for GR and the scale of QG theory. Some attempts to elaborate a full theory of QG were made but, due to the high energy scale of the theory, there is no laboratory experiment sensitive enough to test them. Nevertheless, some classes of models share a common property which can be used indirectly to probe QG: Lorentz Invariance Violation (LIV). In those models, some effects impacting the propagation of very high energy photons (i.e. GeV and beyond) in vacuum leads to a Modified Dispersion Relation (MDR) inducing a LIV. As a result, looking for LIV with astrophysical experimental tests has emerged as a promising way to indirectly constrain or reject some QG models.

At first in this chapter, the main motivations to establish a QG theory are discussed. Secondly, different models of QG are briefly introduced with their specificity and also their limitations. Then, two specific experimental tests are presented, based on MDR to search for LIV effect: the photon time of flight and the γ-γ absorption interaction modification. Finally, focusing on the time of flight approach, the astrophysical sources of interest to constrain LIV are described as well as the best results available in the literature.

Current theory issues

General relativity

General Relativity is the theory of the gravitation interaction. It describes the geometry of space-time, deformed by the presence of matter or radiation. GR was established by Einstein (1915a). It allows for the resolution of many issues of Newtonian gravity like the anomalous advance of Mercury's perihelion. The observed value was not in agreement with the Newtonian prediction but compatible with the GR one, computed by Einstein (1915b). In the same paper, he also computed the bending of star light by gravity, obtaining a different value compared to Newton predictions. Observations done by [START_REF] Dyson | A Determination of the Deflection of Light by the Sun's Gravitational Field, from Observations Made at the Total Eclipse of May 29, 1919[END_REF] to measure the deviation of light by the Sun during the eclipse on 29 th May 1919 were in agreement with the value computed by Einstein, providing one of the first tests of GR. The main effect of light deviation by masses is the gravitational lensing which causes multiple images of the same light emitter if there is a large mass distribution in the light path to the observer (example on Figure 2.1 with extragalactic sources). Finally, one of the most recent and remarkable success of GR is the first detection of gravitational waves by [START_REF] Abbott | Observation of Gravitational Waves from a Binary Black Hole Merger[END_REF]. They were suggested by [START_REF] Poincaré | Sur la dynamique de l'électron[END_REF] and predicted by Einstein (1916) thanks to GR. Nonetheless, GR has some limitations. First of all, the theory does not include naturally dark matter and dark energy which are required for the theory to be in agreement with observations. On the one hand, dark matter is needed to explain the shape of rotational velocities of galaxies (Begeman et al., 1991). One would expect the rotational velocity to decrease as the distance to the center of the galaxy increases because of the distribution of visible matter in the galaxy. But this is in disagreement with actual measurements, presented in Figure 2.2. Also, dark matter explain another inconsistency related to the apparent mass of galaxy clusters. The mass of a cluster estimated using the virial theorem or the gravitational lens effects do not lead to the same value as the one deduced with the luminosity of the cluster. This made [START_REF] Zwicky | Die Rotverschiebung von extragalaktischen Nebeln[END_REF] postulates the existence of a non-luminous mass in the Coma cluster coming from a "dark matter". A last example for the requirement of dark matter is related to the evolution of the Universe, and the formation of large scale structures. From the Cosmic Microwave Background (CMB), a remnant emission coming from the epoch of recombination between proton and electron forming the first atoms, the primordial fluctuations of the Universe can be deduced. The CMB photons are the first able to propagate freely in the early Universe. The primordial fluctuations trace the origin of the formation of structure during the cosmological history. If only baryonic matter is considered, the primordial fluctuations become washed out by the interaction of radiation on matter during the radiation-dominated era leading to a universe without any galaxy cluster or galaxy. Dark matter, which is not affected by any kind of radiation (photons and neutrinos) is required to preserve these fluctuations in order to form the structures visible today. On the other hand, dark energy is necessary in the theory to explain the acceleration of the Universe expansion, discovered by [START_REF] Riess | Observational Evidence from Supernovae for an Accelerating Universe and a Cosmological Constant[END_REF] and [START_REF] Perlmutter | EPOS LHC: Test of collective hadronization with data measured at the CERN Large Hadron Collider[END_REF]. Those two ingredients are keys to explain the evolution of the universe, in the ΛCDM (Lambda Cold Dark Matter) model. But their nature is still unknown whereas dark matter is believed to be non-baryonic matter [START_REF] Peebles | Dark matter and the origin of galaxies and globular star clusters[END_REF][START_REF] White | Expected signatures from hadronic emission processes in the TeV spectra of BL Lacertae objects[END_REF]. Investigations are ongoing to look for the nature of these two constituents. Some other issues are related to Black Holes (BH). In the GR theory, BH have singularities inside them, implying that the theory diverges. The GR breaks down and is not sufficient to describe the physics inside of a BH, at such high curvatures of space-time. Another problem arises when looking at the Hawking radiation mechanism [START_REF] Hawking | Particle creation by black holes[END_REF] which allows for the evaporation of BH by mass losses. After some time, a BH will ultimately reach the size l P with mass m P which is the limit of the Hawking radiation calculation. In other words, GR is not able to describe the fate of BH after complete evaporation.

Quantum Field Theory

Quantum Field Theory is the theory merging together quantum mechanics and special relativity. QFT led to Quantum Electrodynamics, describing electro-magnetism and the weak interaction, and Quantum Chromodynamics describing the strong interaction. These interactions, together with all the particles define the Standard Model. SM has shown an excellent capability to predict the existence of several particles with a precise value of their masses. Few examples can be mentioned such as the discovery of the W boson [START_REF] Arnison | Experimental observation of isolated large transverse energy electrons with associated missing energy at sqrt(s)=540 GeV[END_REF] and the top quark [START_REF] Abe | Observation of Top Quark Production in Pbar-P Collisions with the Collider Detector at Fermilab[END_REF]. More recently the first detection of Higgs bosons [START_REF] Aad | Observation of a new particle in the search for the Standard Model Higgs boson with the ATLAS detector at the LHC[END_REF][START_REF] Chatrchyan | Observation of a new boson at a mass of 125 GeV with the CMS experiment at the LHC[END_REF] allows for the confirmation of the Brout-Englert-Higgs mechanism (Higgs, 1964) which gives the masses of all particles in the theory. Without the Higgs, the theory is either massless or violates energy conservation.

Despite the success of the SM, there are some issues remaining. One of them is QFT being only applicable on a flat space-time. This precludes the use of QFT for a full description of the Hawking radiation. In addition, the speculative graviton which is believed to be the particle carrying the gravitational interaction in the frame of the SM, has not been found so far. Only limits are known about the physical properties of the graviton, such as the one obtained with the detection of gravitational waves [START_REF] Abbott | Observation of Gravitational Waves from a Binary Black Hole Merger[END_REF]. Besides, compared to other interactions, gravity has an extremely low coupling constant, a factor 10 32 weaker than the weak interaction. This is not clearly understood and known as the "hierarchy problem". In addition, SM has around 20 extra parameters which have to be included in the theory in order to match observations. One of them is the ratio between the masses of electron and muon, which has to be set up at a value about 207.

Another issue is about neutrinos. The SM predicts three neutrino flavors (ν e , ν µ and ν τ ) which do not have masses. This is a consequence of the theory which predicts neutrinos to be only left-handed, preventing a Yukawa interaction to give them masses. In 1998, the Homestake Solar Neutrino Detector measured the neutrino flux coming from the Sun [START_REF] Cleveland | Reaction rates and energy distributions for elementary processes in relativistic pair plasmas[END_REF]. They reported some inconsistencies on electron neutrino flux coming from the 8 B and 7 Be nuclear reactions compared to solar models. Few months later, the Kamiokande experiment also reported a deficit of neutrino coming from the Sun [START_REF] Fukuda | Evidence for Oscillation of Atmospheric Neutrinos[END_REF]. A neutrino oscillation mechanism between the different flavors was thus proposed to explain this deficit of neutrino flux, as suggested by [START_REF] Pontecorvo | Mesonium and Antimesonium[END_REF]. These oscillations imply that neutrinos are massive which is not consistent with SM and requires an extra mechanism to generate their masses.

A last example is given by the proton size problem within Quantum Electrodynamics. This problem is not directly related to QG but can motivate the necessity of a model beyond the SM. The proton size defines the Rydberg constant as its inverse value. Before 2010, the proton size was well known and measured from two different techniques giving compatible values. The first method relies on the Rutherford experiment. When looking at the scattering of electrons fired on a proton target, the deflection angles of electrons give access to the proton size. The second method makes use of the energy levels of the electron in the Hydrogen atom. Via spectroscopy, the energy of the emitted photon depends on the Rydberg constant and so on the proton size. In 2010, [START_REF] Pohl | The size of the proton[END_REF] used this second technique but with a muonic hydrogen, i.e. the electron of hydrogen is replaced by a negative muon. Because of the heavier mass of the muon, the precision on the measurement is higher. The result presents a 4% lower value than the accepted r p = 0.84184(67) fm [START_REF] Mohr | CODATA recommended values of the fundamental physical constants: 2006[END_REF], about a 4 σ deviation. More recently, Beyer et al. (2017) made the measurement again using a hydrogen atom but using another transition frequency. This new result is 5% lower than the value by [START_REF] Mohr | CODATA recommended values of the fundamental physical constants: 2006[END_REF] but is compatible with [START_REF] Pohl | The size of the proton[END_REF]. In this regard, one has either to assume that something is wrong in the theory, or to find a way to explain the discrepancy between the measurements.

Effective field theory

To get a common framework to describe all fundamental interactions, one would like to simply combine GR and QFT. However, it turns out to be impossible in practice. One naive way to point it out is that GR is a deterministic theory while QFT is a probabilistic one. These two approaches in that sense are not compatible with each other. In a more rigorous way, when trying to quantize gravity to introduce it in QFT, divergences appear and cannot be removed resulting in a non-renormalizable theory. This comes mostly from the coupling constant of gravity G being dimensionful (Equation 2.2 with hc = 1 and G = 1/m 2 P ). When applying perturbation theory, two powers of the momentum come to dimensionally compensate the two powers from the Planck mass. This worsen ultraviolet divergences and would require an infinite number of counter-terms to remove them compared to the usual SM which needs a finite number of counter-terms. Nevertheless, this approach is still of interest and is called Effective Field Theory (EFT). It consists of using the low-energy predictions of this theory, where it converges, to investigate Physics beyond the currents theory. As an example, the use of an EFT allowed for the discovery of the Hawking radiation [START_REF] Hawking | Particle creation by black holes[END_REF].

The approach to Quantum Gravity

The previous section shows that both GR and QFT have limitations and inconsistencies. Given these points, a theory beyond GR and QFT is required, possibly capable of solving those issues. QG is hoped to be the answer. Currently there is no complete QG theory, but rather several models that may lead to the QG theory in the future. In the following, the two most popular of them will be shortly explained with few of their accomplishments and open questions.

Loop Quantum Gravity

This part benefits from the reading of Rovelli (2004) and [START_REF] Gambini | A First Course in Loop Quantum Gravity[END_REF] As seen in the previous section, quantizing gravity leads to a non-normalizable theory in the QFT framework. This mainly comes from the requirement of QFT to describe a phenomenon on a background: the flat Minkowski space-time. So, as gravity is defined by the space time itself in GR, trying to describe gravity in QFT removes the background on which QFT is applicable. This is directly connected to the non-renormalizability of the EFT approach. In other words, the background-dependence of the QFT leads to the failure of quantizing gravity in this framework. Loop Quantum Gravity (LQG) is a tentative to define a background-independent QFT. On the contrary of QFT which uses fields, LQG uses Wilson loops. The reason for this is that fields require a background space-time to be used which is not compatible with LQG.

The first success of LQG is that it provides a general-relativistic QFT, with a background independent formulation. The goal of quantizing gravity is achieved through spin networks defining localization and spatial relations. The ultraviolet divergences, present in QFT and removed by renormalization, are non-existent in LQG which gives a complete description of quantized gravity. Furthermore, LQG provides a quantum cosmology without singularities, giving a minimal size to the initial state of the Universe. Finally, LQG has the benefit of having precise physical predictions. But these predictions are currently out of reach of current detection capabilities and so there is no way to directly test if LQG is the final QG theory. This theory has also some issues not resolved yet. One of them is that the classical world is difficult to recover from LQG: computing solutions of Einstein equations at low energy is quite complex. Some work is being done to achieve that goal but it is still ongoing.

String theory

This part benefits from the reading of [START_REF] Zwiebach | A First Course in String Theory[END_REF] A more ambitious approach to QG is String Theory (ST). It was initially developed in the 1960s, to describe the strong interaction but was disfavored over Quantum Chromodynamics. It is only when QG started to be investigated that theorists realized that ST is well suited for QG. One of its advantage is that it is a truly unified theory. While at present, there are 29 fundamental particles in SM (12 particles, 12 anti-particles and 5 boson vectors), ST offers only one fundamental particle: the strings. Each particle of the SM is a vibrational mode of a string at a specific frequency. On the contrary to the SM, ST has also the beauty to only have one adjustable parameter: the string length l s . The lack of adjustable parameters implies that the theory cannot be easily changed or deformed which prevents the theory to provide an infinite number of predictions. ST does not describe the Universe in 4 dimensions (3 spatial and 1 temporal dimension) but rather 11 or 24 depending on the type of string theory considered (bosonic string theories or superstring theories). In those theories, our Universe exists on a 4-D hypersurface. Hypersurfaces are called D-branes which are physical objects in ST.

As LQG, ST suffers from the problem of having no predictive power at low energy scales. The extra spatial dimensions of the theory are hoped to be detected by gravitational experiments. If D-branes exist, some deviation of the Newton law about the attraction of two masses are expected to show up but only at very small distances. With the low gravity coupling constant G, only distances down to very few nanometers are investigated [START_REF] Haddock | Search for deviations from the inverse square law of gravity at nm range using a pulsed neutron beam[END_REF]. Besides, the parameter l s is mostly believed to be of the order of l p which is out of reach for any current instrument. It also strongly relies on supersymmetry, searched for many decades in accelerators but not yet found. Some theorists think that string theory already provided one successful prediction: the unavoidable existence of gravity which arises naturally from relativistic strings. Indeed, none of the classical vibration of strings is able to produce the particle of gravity. Only the quantum vibrations of relativistic strings corresponds to the graviton. Nevertheless, ST, as LQG, offers a good candidate for QG despite the lack of experimental evidence.

Probing Quantum gravity via Lorentz invariance violation

Lorentz invariance can be used to test some classes of model of QG. LIV is not a general prediction of all QG models but a possibility for some of them. As seen in the last section, EFT allows for the establishment of predictions which can be of interest for LIV studies. A work by Colladay & Kostelecký (1998) includes the gravitational force in the SM with a single Lagrangian. All the terms which can violate the Charge-Parity-Time (CPT) symmetry and Lorentz invariance appear and allow for comprehensive framework to search for such violations. This define the Standard Model Extension (SME). An annually updated table of parameters is published by [START_REF] Kostelecký | Data tables for Lorentz and CPT violation[END_REF]. In the following, only the photon sector is investigated. However, a simpler test theory is favored in this thesis which is described in this section.

Modified dispersion relation

LIV is not a requirement in LQG or ST models but can appear in some them. LIV can emerge from a Modified Dispersion Relation (MDR) for photons in vacuum with for example interactions between photons and D-branes [START_REF] Amelino-Camelia | Distance Measurement and Wave Dispersion in a Liouville-String Approach to Quantum Gravity[END_REF][START_REF] Ellis | Quantum-Gravitational Diffusion and Stochastic Fluctuations in the Velocity of Light[END_REF] in some ST model, or because of the discrete nature of space-time [START_REF] Gambini | Nonstandard optics from quantum space-time[END_REF] in some LQG model. This MDR can lead to promising observables. To simplify experimental investigations, all those specific models are represented with a unique toy-model based on the standard dispersion relation E 2 = p 2 c 2 and supposing a LIV effect to occur at an energy scale E QG , expected to be close to the Planck energy. This MDR is expressed as: (2.4) where the sign ± allows for superluminal and subluminal effects. If the LIV effect is linearly or quadratically suppressed with the energy, then experimental observations are believed to be sensitive enough to detect it. Derived from this MDR, two main observables inducing LIV will be discussed in the following.

E 2 p 2 c 2 ×   1 ± ∞ n=1 E E QG n   ,

Time of flight

From the MDR of Equation 2.4, one can derive the resulting velocity of photons, v γ (E), which become energy-dependent (See Appendix 7.5). It is expressed, only keeping the leading order term, as:

v γ (E) = ∂E ∂p = c ×   1 ± n + 1 2 E E QG n   . (2.5)
With such a velocity, two photons emitted at the same time with different energies should travel at different velocities and thus arrive with a time lag at a distant observer location.

As E QG is expected to be of the order of E P , such an experimental test cannot be done on Earth because the effect would be extremely small. But one can make the most of the Universe as a laboratory, using astrophysical sources as photon emitters and accelerators. Those sources are part of the most energetic accelerators known up to date, and are also extremely distant making them very well suited to measure the time difference induced by this energy dependent velocity. In this time of flight approach, describe by Equation 2.5, ST models such as [START_REF] Amelino-Camelia | Distance Measurement and Wave Dispersion in a Liouville-String Approach to Quantum Gravity[END_REF], [START_REF] Ellis | Quantum-Gravitational Diffusion and Stochastic Fluctuations in the Velocity of Light[END_REF] produce only a subluminal (+) linear (n = 1) effect whereas with the LQG model from [START_REF] Gambini | Nonstandard optics from quantum space-time[END_REF], both subluminal and superluminal (-) can be produced, each of them with a linear or quadratic (n=2) effect with the energy.

However, one has to be careful and to take into account the expansion of the Universe in the computation. Travelling at slightly different speeds than c, photons will not travel the same distance because the Universe is expanding during their propagation. Considering a source at a redshift z, Jacob & Piran (2008) have computed the time delay between two photons of energies E 1 , E 2 from Equation 2.5, including the cosmological expansion: .6) In this expression, Ω m and Ω Λ are the cosmological parameters of the ΛCDM (Lambda -Cold Dark Matter) model, respectively matter and dark energy densities, and H 0 is the Hubble parameter. Their values are taken from the latest Planck mission results [START_REF] Ade | Planck 2015 results. XIII. Cosmological parameters[END_REF]:

∆t = ± 1 + n 2H 0 E n 2 -E n 1 E n QG z 0 (1 + z ) n dz Ω m (1 + z ) 3 + Ω Λ . ( 2
Ω m = 0.3089 ± 0.0062, (2.7
)

Ω Λ = 0.6911 ± 0.0062, (2.8) H 0 = 67.74 ± 0.46 km s -1 Mpc -1 .
(2.9)

Cross section modification

A MDR can also induce a modified photon kinematics, which can lead to modifications of the γ-γ interaction cross-section. This effect can be probed when observing distant extragalactic objects. TeV energy photons can be absorbed because of the interaction between these photons with the Extragalactic Background Light (EBL). EBL is the optical/UV light coming from all the galaxies in the Universe through its cosmological evolution including its diffusion by dust. It suppresses the TeV γ-ray flux above a certain energy, depending on the redshift of the source. A LIV effect can modify the interaction between γ rays and EBL photons, leading to a change of absorption. To constrain such an effect, one needs both at the same time a distant source to have significant EBL absorption, and hard power spectrum from the source. This condition greatly limits the number of sources usable for such a study.

Sources of interest for time of flight studies

This work focuses on the time of flight approach to probe LIV. There are several astrophysical sources of interest for this method. The few crucial properties required are the following:

• Distance: the farther away the source is, the more important the time delay difference will be. Due to the expansion of the Universe, a photon with higher velocity will travel less distance than a slower photon, thus inducing a larger LIV time delay.

• High energy range: the energy is the key ingredient in order to be able to detect a time delay. The photons energy difference comes in the equation, and requires to be able to observe emitter on the widest energy range possible. High energies are also essential for investigating the second order term in the MDR.

• Variability: one needs to have sharp variability in order to be able to measure a time delay.

If the source emits a constant flux over time, no delay can be measured. The sharper the variability of the source will be, the more precise will be the measurement of a time lag.

The different kinds of sources used for time of flight studies are going to be presented below as well as the most important results available in the literature.

Gamma-ray bursts

Gamma-Ray Bursts (GRB) are transients which release a colossal amount of energy, about 10 51-54 ergs, in a fractions of seconds up to a few minutes. They were first detected by a set of military satellites called Vela, launched by the U.S. during the 1960s. Vela primary mission was to monitor nuclear activity during the Cold War. They measured γ-ray flashes of unknown origin which lead to the discovery of GRB (Klebesadel et al., 1973) after public release of the data. Then, there was a debate about the origin of GRB: galactic or extragalactic. In 1991, NASA launched a dedicated satellite, the Compton Gamma-Ray Observatory, for γray observations with 4 different instruments on-board. One of them, the Burst and Transient Source Experiment, allowed for the measurement of the distribution of GRB in the sky, revealing an isotropic distribution favoring the extragalactic origin [START_REF] Meegan | Spatial distribution of gamma-ray bursts observed by BATSE[END_REF].

GRB usually show a fast γ-ray prompt emission, observed up to tens of GeV, followed by an afterglow spreading in energy from radio to X-rays. There are two classes of GRB depending on their duration. Long GRB last a few seconds to a few minutes and represent the vast majority of GRB. They originate from collapse of massive stars, as they can be linked to star-forming galaxies and some of them coincide with supernovae. Short GRB are typically less than 2 seconds long. Their origin is more difficult to determine, as they are rarer. The observation of one of them, GRB 130603B shows flux values compatible with a kilonovae [START_REF] Tanvir | A 'kilonova' associated with the short-duration γ-ray burst GRB 130603B[END_REF]. More recently, the observation of GRB 170817A (Abbott et al., 2017a) associated with a gravitational wave event (Abbott et al., 2017b) allowed for the determination of this GRB origin coming from a kilonovae. The gravitational wave signal which identifies the two merging objects as neutron stars.

GRB are promising candidates for LIV search thanks to their fast variability (see Figure 2.3 for example) and their high energy emission. They are also observed at large distances, from the closest at z = 0.0085 to z = 9.4 for the farthest. However, distance comes with a cost which is the maximum energy observed due to EBL absorption. The further photons will travel, the more probable they will be absorbed during propagation leading to an attenuation, or even a suppression of the high energy flux. Moreover GRB are random events which make their detection quite difficult and require wide-field detectors, mostly satellites which limit the number of detected photons due to their small collection area. Up to now, ground based telescopes have not been able to catch a GRB due to their small fields of view and large reaction times, between the reception of an alert and the repositioning of the telescopes. Thus the observed energy range of GRB is restricted by the combination of the small collection area of satellites and the very small high energy photon flux. The red vertical line corresponds to the trigger time of the observation, while the green ones are used for the spectral analysis on the paper from which the plot is extracted [START_REF] Ackermann | Fermi Observations of GRB 090510: A Short-Hard Gamma-ray Burst with an Additional, Hard Power-law Component from 10 keV TO GeV Energies[END_REF].

Pulsars

A pulsar is a fast rotating magnetized neutron star. Neutron stars originate from massive stars which have consumed all their nuclear fuel. During the collapse, the star goes from a typical size of 10 6 km to tens of kilometers at the final state. By conservation of angular momentum, this results in a very fast rotation. The magnetic field is also significantly increased during the subsequent collapse, up to typically 10 12 G. Pulsars emit narrow beams of light ranging from radio to γ rays, and due to their rotation, they are seen as pulsating. The period ranges from tens of seconds down to a few milliseconds. The first pulsar, detected by J. Bell during her PhD [START_REF] Hewish | Observation of a Rapidly Pulsating Radio Source[END_REF], was PSR B1919+21, originally named "Little Green Man 1" because of the unknown and extremely regularity nature of the pulsating emission at that time. Variability for pulsars is related to the rotation period. A phasogram, which is a folded light curve into the pulsar period (see example Figure 2.4), is used instead of the full light curve. The phase is defined as the arrival time expressed as a multiple of the pulsar period. Since they are persisting sources in the sky, they can be observed for very long times in order to accumulate statistics. This allows for the determination of a very detailed phasogram, sometimes showing millisecond variability. Also, the periodicity allows for a great improvement of the signal over background ratio, by only considering photon arriving on phase with respect to the pulse shape. Up to now, only two pulsars were detected at TeV energy: the Crab pulsar (d = 2 kpc) and Vela (d = 0.29 kpc).

Active galactic nuclei

Among galaxies, some exhibit a very energetic core, more than what stars within them can produce. These extremely luminous cores are called Active Galactic Nuclei (AGN). The core itself is composed of a super-massive BH, with strong accretion power leading to a disk of matter surrounding the BH. It exists a huge number of AGN categories (more than 50). This large number of categories came from the misunderstanding of their stucture, and the huge variety of different properties that AGN show. Some of them were already reported in optical catalogs like the "Catalogue des nébuleuses et des amas d'étoiles" [START_REF] Messier | Catalogue des Nébuleuses et des Amas d'Étoiles (Catalog of Nebulae and Star Clusters[END_REF], but later on they were recognized to emit unusual bright and large spectral emission lines. They were named as "Seyfert galaxies" due to the classication by [START_REF] Seyfert | Nuclear Emission in Spiral Nebulae[END_REF] using emission line widths. Similarly, when first radio catalogs were published (Edge et al., 1959), some very compact objects associated with optical point-like sources were detected with strong variability, strong emission in radio and spectral emission lines. They were named quasars, standing for "quasi-stellar objects". A lot of different objects, all known as AGN today, were defined when first discovered due to the variety of attributes that AGN can have: presence or absence of radio jets, strong variability in optical, presence or absence of spectral emission lines, broad or narrow emission lines and many more. It's finally [START_REF] Antonucci | Unified models for active galactic nuclei and quasars[END_REF] for radio-quiet sources and [START_REF] Urry | Unified Schemes for Radio-Loud Active Galactic Nuclei[END_REF] for radio-loud ones who proposed a unification model to merge all those different objects in a single one able to explain the full zoology of AGN. This unification is based on the angle on which the object is seen. This is illustrated in Figure 2.5. Only one important dichotomy remains, which is the presence or absence of radio jets in respectively radio-loud and radio-quiet objects. As a remark, there is another sort of unification model, based on the evolution of the host galaxy of the AGN [START_REF] Hickox | Host Galaxies, Clustering, Eddington Ratios, and Evolution of Radio, X-Ray, and Infrared-Selected AGNs[END_REF][START_REF] Goulding | Tracing the Evolution of Active Galactic Nuclei Host Galaxies over the Last 9 Gyr of Cosmic Time[END_REF], Azadi et al., 2015).

Only blazars are used for LIV studies because of their better γ-ray flux, energy and variability with respect to other classes of AGN. More details about blazars will be discussed in Chapter 3. Blazars are most of the time seen in their low emission state but they can occasionally go in high states and even show some flaring activity. Flares can have observed variability time scales ranging from hours down to a few minutes at γ-ray energies (see example Figure 2.6) which is needed for LIV searches. They are observed at energies up to tens of TeV for the closest ones and hundreds of GeV for the farthest, up to z ∼ 1, due to EBL absorption. Unfortunately, flares are random transients and can only be seen during night time with ground based telescopes. Besides, the redshift is not always easy to determine, especially when emission lines are hidden by the jet emission itself. Blazars appear to be at the moment the best compromise for deep LIV studies with the time of flight approach and are the main sources used for this work.

Results on LIV studies

Different measurements were done using the different sources presented above which all have their specificities and interests. When measuring time delays, only linear and quadratic terms in the MDR are investigated. Higher orders are not accessible with current instruments in use at the moment. A summary of the best current limits is presented in Table 2.1. One can note that GRB 090510 outclasses by far all other sources on the limit on E QG,1 which exceeds the Planck scale for a linear model. This observation is the only one able to put such a stringent limit by its high number of detected events up to 30 GeV coming from a redshift z = 0.903. But this GRB is the only observed with these properties so far and more of them are required to confirm this limit.

For AGN, PKS 2155-304 is dominating the limits with its exceptional characteristics. Its light curve is composed of five consecutive short flares and ∼ 10000 photons were recorded by H.E.S.S. in 90 minutes. On the contrary, PG 1553+113 suffers from low statistics (about a hundred photons detected) but benefits from a high redshift value (z ∼ 0.5). The last AGN, Markarian 501, has low redshift (z = 0.034) explaining the modest linear limits but it shows photon being detected above 10 TeV which results in a competitive limit on E QG for n = 2 compared to the one established with the PKS 2155-304 flare.

Finally, pulsars provide in overall weaker limit values on the QG energy scale due to their close distance. However, they greatly benefit from their fast variability down to milliseconds, thus providing good constraints on E QG,2 . As a example, the last result from MAGIC Collaboration et al. ( 2017) uses very long observation time of the Crab pulsar (1000 hours distributed over 12 years) which allows for the establishment of a very competitive limit despite the drawbacks of the Crab pulsar being a local source in our Galaxy.

Summary

The current theories of General Relativity and Quantum Field Theory have been very successful and allowed for the determination of predictions and new discoveries. They are able to describe the four fundamental interactions and all the particles of the Standard Model. However, a new theory is necessary to truly unify gravity with quantum mechanics, named as the Quantum Gravity theory. Several theoretical models, as String Theory or Loop Quantum Gravity, attempt to describe the full Quantum Gravity theory, but none of them predicts testable observables, preventing their validation. Nevertheless, some classes of models presents properties like Lorentz Invariance Violation that can still be tested in order to constrain those models.

The time of flight approach to constrain LIV is the one used in this work, from a Modified Dispersion Relation affecting the velocity of photon in vacuum. This velocity becomes energy dependent which leads to time delays between photons with different energies, especially coming from distant astrophysical sources. From the different cosmic sources presented here, blazars are the ones considered in this thesis. The method used to measure energy dependent delays in the arrival time of photons as well as an application on data will be presented in Chapter 6. But one needs to assume that during the emission of high energy photons within the source, there is no delay introduced by the source itself. By doing so, any measured time delay could be considered only coming from LIV effect which would then allows for the establishment of QG model constraints.

The next chapter focuses on blazars which are used for this work to search for LIV signatures. These object are going to be described as well as the mechanisms which allows for particle acceleration and photon emission. The term "blazar" comes from the merging of two categories of AGN: BL Lac objects and Optically Violent Variability quasars. The term BL Lac object itself comes from the misidentification of the AGN BL Lacertae as a variable star. Blazars are the most promising objects for LIV study among the AGN. As already discussed, they emit very high energetic photon flux, as observed up to tens of TeV, and show variability down to the minute scale. In this chapter, the main characteristics of blazars are introduced. Then the processes responsible of the Spectral Energy Distribution (SED) observed in blazar are described, followed by the particle acceleration mechanisms, allowing for the generation of particles with a sufficient energy to emit the observed SED. Finally, the main emission models using the different non-thermal emission processes are presented.

Characteristics of blazars

Jets and superluminal motion

Blazars belong to the family of radio-loud AGNs. Radio-loud AGNs present radio jets as shown on Figure 3.1 for the radio-galaxy M87. Inside jets, denser zones can be distinguished, called "knots". In the case of blazars, detailed observations of the jets are difficult due to their specific orientation, pointing towards the Earth. Nevertheless, knots can be detected and tracked in some cases for blazars as shown in Figure 3.2. Actually, in some radio-loud AGNs, these knots can be seen moving at superluminal velocity (v obs > c).

These superluminal apparent velocities are explained by a simple geometrical argument, as illustrated in Figure 3.3. As the jet moves toward the observer at velocity v ∼ c, the time difference ∆t obs between the observed photons does not reflect the emission time difference ∆t. With a small viewing angle θ, the observed distance traveled by the knot is ∆L obs = v∆t sin θ while the observed time is ∆t obs ∆t -v c ∆t cos θ. In overall, the observed velocity is:

v obs = ∆L obs ∆t obs = v sin θ 1 -β cos θ , ( 3.1)
with β = v/c. For a knot with v ∼ c, the apparent velocity v obs can exceed c when the observing angle θ is small enough. The maximum velocity is reached for cos θ = β which gives:

v obs = v sin θ 1 -β 2 = v √ 1 -cos 2 θ 1 -β 2 = v 1 1 -β 2 = Γv (3.2)
where Γ = 1/ 1 -β 2 is the Lorentz factor of the knot. The observation of superluminal velocity is a direct evidence of relativistic motion inside jets. For blazar, one of the jet is pointing toward the observer and photons benefit from a blue-shift coming from a Doppler boosting effect.

Doppler boosting

The Doppler boosting effect is the consequence of the time dilatation induced by Lorentz transformation. For an object moving at velocity v = β × c and observed at a small angle θ, the characteristic Doppler factor δ is expressed as: The time dilatation and its subsequent effects imply that:

δ = 1 Γ (1 -β cos θ) , (3.3)
t obs = t s δ ; ν obs = δν s ; E obs = hν obs = δhν s = δE s . (3.4)
Consequently, from the Doppler boosting, the observed variability time scale is reduced and the observed frequencies are increased. Besides, the Doppler boosting also affects the observed photon flux. As noted in [START_REF] Rybicki | Radiative processes in astrophysics[END_REF], I(ν)/ν 3 is a Lorentz invariant as it is simply the number of particle within a volume, with I(ν) being the specific intensity at frequency ν. I(ν). Therefore, one can write:

I obs (ν obs ) ν 3 obs = I s (ν s ) ν 3 s , (3.5) I obs (ν obs ) = I s (ν s ) ν 3 obs ν 3 s = I s (ν s )δ 3 . (3.6)
The observed specific intensity is boosted by a factor δ 3 . Finally, concerning the total intensity integrated over all the energies, the effect is further enhanced with:

I obs = ∞ 0 I obs (ν obs ) dν obs = ∞ 0 δ 3 I s (ν s ) δdν s = δ 4 I s . (3.7)
On the whole, Doppler boosting effect increases simultaneously the observed flux and the observed energy as well as the observed variability. All these observables benefit to the time of flight analysis for LIV, enhancing the major physical properties required from the source for these studies.

Blazar sequence

Within the unified scheme of AGN, there are several types of blazars: Flat Spectrum Radio Quasar (FSRQ) and BL Lac objects. These categories follow the so called "blazar sequence" based on their radio luminosity, originally defined by Fossati et al. (1998). A total of 126 blazars, split into 5 GHz radio luminosity bins, were used to produce an average SED for each band. The resulting SEDs are shown on Figure 3.4 (right). As a result, a specific trend called the blazar sequence is observed: when the luminosity increases, the SED peaks at lower energy. The blazar sequence goes from bright and low frequency sources, FSRQ, to weaker and higher frequency ones, High BL Lac object (HBL), with in between Low and Intermediate BL Lac object.

A recent study from [START_REF] Ghisellini | The Fermi blazar sequence[END_REF] attempts to update the blazar sequence. A total of 747 blazars were used, increasing significantly the statistics from the previous sequence. The luminosity bins used to construct the SED are also changed, and defined at γ-ray energies instead of radio. This updated sequence is shown on Figure 3.4 (left) and shows similar behavior as the previous one.

In the initial and revised blazar sequences, FSRQs have the highest luminosity while HBL have the lowest. This translates into HBL being generally observed with a lower redshift value than FSRQ. Nevertheless, the sequence is still under debate and some exceptions to the sequence have been reported [START_REF] Padovani | The discovery of high-power high synchrotron peak blazars[END_REF][START_REF] Cerruti | Luminous and high-frequency peaked blazars: the origin of the γ-ray emission from PKS 1424+240[END_REF].

Non-thermal radiation processes

For blazars, several mechanisms exist which can produce photons. As shown in Figure 3.4, typical SEDs of TeV blazars present a two-bump structure which suggest at least two main emission mechanisms. This section presents the main interactions able to produce photon which depict the blazar SEDs.

Synchrotron emission

The detailed synchrotron emission description is inspired by [START_REF] Rybicki | Radiative processes in astrophysics[END_REF].

Accelerated charged particles radiate electromagnetic waves as described by the Lienart-Wieckart potential [START_REF] Rybicki | Radiative processes in astrophysics[END_REF]. For classical particles moving within a magnetic field -→ B , the process is called cyclotron emission, while for relativistic particles, it is called synchrotron emission.

Considering a charged particle accelerated in a homogeneous magnetic field of strength B, with a pitch angle θ between its trajectory and the magnetic field lines, with Lorentz factor of γ = E/mc 2 , of charge q, the emitted power at frequency ν by the particle is expressed as:

P (γ, θ, ν) = √ 3q 3 B mc 2 sin θ ν ν c ∞ ν νc K 5 3 (x)dx, (3.8)
with ν c = 3qB 4πmc 2 γ 2 sin θ, (3.9)

where K 5 3 (x) is the modified Bessel function of second kind of order 5/3. The general formula for K n (x) is:

K n (x) = π 2x e -x n -1 2 ! ∞ 0 e -y y n-1/2 1 - y 2x n-1/2
dy.

(3.10)

The synchrotron emissivity j s (ν) is then deduced from P (γ, θ, ν). It corresponds to the synchrotron photon at frequency ν. Considering a distribution of particles N (γ) with γ min < γ < γ max and a isotropic distribution of angle θ:

j s (ν) = 1 8π 2π θ=0 sin θ γmax γ min N (γ)P (γ, θ, ν)dγdθ. (3.11)
The synchrotron photons, in the medium where they are emitted, suffer from absorption due to the same population of emitting particles. The synchrotron self-absorption coefficient is written as:

k s (ν) = 1 16πmν 2 π θ=0 sin θ γmax γ min γ 2 d dγ N (γ) γ 2 P (γ, θ, ν)dγdθ. (3.12)
Using the expression of j s (ν) and k s (ν), one can obtain the synchrotron intensity I rad (ν) which propagates outside the emission zone, called blob. Considering a spherical blob of radius R b , filled with a particle distribution N (γ), I rad (ν) is the solution of the transfer equation in a spherical geometry in case of synchrotron dominated emission:

I rad (ν) = j s (ν) k s (ν) 1 - 1 τ 2 1 -e -τ (τ + 1) , (3.13) with τ = 2R b k(ν). (3.14)
Finally, the average emitted power can be deduced from Equation 3.8. It corresponds to the average energy loss by an electron via synchrotron radiation. Considering particles with an isotropic distribution of velocity, the mean synchrotron emitted power P synch (γ) of one particle, integrated over all emitted frequencies and averaged over the pitch angles θ, is:

dE dt syn (γ) = P synch (γ) = 4 3 σ T cβ 2 γ 2 U B (3.15) with σ T = 8πr 2 p 3 , U B = B 2 8π , (3.16)
where σ T is the Thomson cross-section for a particle of radius r p , and U B is the magnetic energy density.

Inverse-Compton interaction

Particle frame The Inverse Compton (IC) scattering is the interaction between a photon and a relativistic electron. During the interaction, an energy transfer occurs from the electron to the photon, resulting in a boost of the photon energy. A scheme of the interaction is presented in Figure 3.5. A photon of energy E i in the observer frame has an energy E i in the frame of the target electron of Lorentz factor γ such that:

E i = E i γ(1 -cos θ), (3.17)
with θ being the angle of approach of the photon in the observer frame. After the interaction, the final energy E f of the photon in the electron frame is:

E f = E i 1 + E i mec 2 (1 -cos θ 1 ) , (3.18)
Blazars with θ 1 being the angle of exit of the photon in the electron frame and m e the electron mass.

Then, the energy in the observer frame can be deduced with the inverse Doppler formula:

E f = E f γ(1 + cos θ 1 ). (3.19)
In the case where E i m e c 2 , the photon energy does not change in the electron frame:

E i = E f .
However, in the case where θ ≈ π and θ 1 ≈ 0, the final energy of the photon in the observer frame is:

E f = 2E f γ = 2E i γ = 4E i γ 2 .
(3.20)

The energy gain from IC interaction is proportional to γ 2 which is a significant increase which for instance allows for the boost of an initial keV energy photon up to TeV energies. The requirement of E i m e c 2 is defined as the Thomson regime, where the cross-section of the interaction is simply σ T (Equation 3.16). Conversely, the Klein-Nishina regime is defined when E i > m e c 2 . In that case the general cross-section is:

σ KN (x) = σ T 3 4   1 + x x 3 2x(1 + x) 1 + 2x -log(1 + 2x) + 1 2x log(1 + 2x) - 1 + 3x (1 + 2x) 2 )   , (3.21)
where x = hν/m e c 2 is the adimensional photon energy. In this regime, the boost of the photon towards higher energies becomes less efficient compared to the Thomson regime. The Thomson cross-section is recovered when x 1.

The IC emissivity j c (ν) at final frequency ν, corresponding to the adimensional photon energy x, is computed by considering a distribution of target electrons N e (γ), and a density of initial photons n γ ( 0 ) at initial frequency ν 0 (or adimensional energy x 0 ). It is given by the formula:

j c (ν) = h 4π x x 0 n γ (x 0 ) γ N (γ)C(x 0 , x, γ)dγdx 0 , (3.22)
where C(x 0 , x, γ) is defined as the Compton kernel. Computed by [START_REF] Jones | Time-dependent Radiation Transfer in the Internal Shock Model Scenario for Blazar Jets[END_REF], the Compton Kernel is expressed as:

C(x 0 , x, γ) = 2πr 2 e c γ 2 x 0 = 2κ log κ + (1 + 2κ)(1 -κ) + (4x 0 γκ) 2 2(1 + 4x 0 γκ) (1 -κ) , (3.23)
where κ =

x 4xγ(γ-x) and r e = e 2 mec 2 is the electron radius. The kinematics of the interaction constrain the resulting energy x to range:

x 0 ≤ x ≤ γ 4x 0 γ 1 + 4x 0 γ . (3.24)
Finally, the average first order IC emitted power P IC (γ) for an electron of energy γ in the Thomson regime, which corresponds to the IC energy losses by the electron, is obtained by integration of Equation 3.22 and gives:

dE dt IC (γ) = P IC (γ) = 4 3 σ T cβ 2 γ 2 U rad (3.25)
where U rad is the energy density of the initial photons involved in the IC interaction.

Hadronic emission processes

Hadronic interactions can lead to the production of γ-ray photons. The main processes considered in blazars involve only energetic proton. For the interactions of protons with photons, there are three different channels:

p + γ → p + π 0 ; p + γ → n + π + ; p + γ → p + e + + e -. (3.26)
Then, the produced pions decay through different channels, depending on their charge. Their main decay channel are:

π 0 → 2γ ; π + → µ + + ν µ ; π -→ µ -+ ν µ . (3.27)
For π 0 , there is a production of photons. For π ± , only muons are produced which can either emit synchrotron photon, or decay into electrons:

µ + → e + ν e ν µ ; µ -→ e -ν e ν µ .
(3.28)

The electrons originating from the muon decay can then either emit synchrotron photon or interact with IC.

Another important way to produce photons is through p-p interactions such as:

p + p → p + p + k(π + + π -) + lπ 0 , (3.29) p + p → p + n + k(π + + π -) + lπ 0 , (3.30) p + p → n + n + 2π + + k(π + + π -) + lπ 0 , (3.31)
with k and l natural numbers which depend on the energy of the initial protons by the energy conservation principle. As for the previous case, the production of π 0 is the only way to directly produce γ-ray photons. The π ± produce secondary particles, i.e. muons and electrons/positrons, allowing for synchrotron and IC emission. One specificity of hadronic processes is the production of secondary neutrinos which are characteristic of hadronic processes. More details on such hadronic processes can be found in [START_REF] Mannheim | The proton blazar[END_REF], [START_REF] Protheroe | High Energy Neutrinos from Blazars[END_REF].

γ -γ interaction absorption

Photons of high energy can be absorbed when encountering lower energy photons via the interaction γ + γ → e -+ e + . The associated cross-section for the interaction is [START_REF] Heitler | Quantum theory of radiation[END_REF]:

σ γγ (E, , θ) = 3σ T 16 (1 -s 2 ) 2s(s 2 -2) + (3 -s 4 ) log 1 + s 1 -s , (3.32)
where E and are respectively the energies of the high and low energy photons, θ is the pitch angle between the photons, s = (1 -4m 2 e c 4 /ζ), ζ = 2E w and w = (1-cos θ). From Equation 3.32, one can compute the energy max and wavelength λ max which maximize the absorption:

max (E) 2(m e c 2 ) 2 E 0.5 1TeV E eV, (3.33) λ max (E) 1.24 E 1TeV µm. (3.34)
Then considering a high energy photon propagating through a medium filled with low energy photons of density n γ ( , r), one can derive the associated optical depth:

τ (E, D) = c D 0 2 0 w 2 ∞ 2m 2 e c 4 ζ dn γ ( , r) d σ γγ (s)d dw dr, (3.35)
with D being the total distance traveled by high energy photons in the medium.

Extragalactic Background Light absorption

The EBL light is UV and optical light emitted by galaxies in the Universe considering also diffusion by dust which can also induces γ-γ absorption. When considering EBL absorption of γ rays emitted by remote blazars, the optical depth formula needs to be modified. Indeed, the distances involved are cosmological and the expansion of the Universe has to be taken into account. In that case, the optical depth for a source at redshift z s is written as:

τ EBL (E, z s ) = c zs 0 dt dz 2 0 w 2 ∞ 2m 2 e c 4 ζ(1+z) dn γ ( , z) d σ γγ (s)d dw dz, (3.36) with dt dz = 1 H 0 (1 + z) (1 + z) 2 (1 + zΩ m ) -z(z + 2)Ω λ -1 2 .
(3.37)

One difficulty arises in order to compute the EBL optical depth which is the evaluation of n γ ( , z). At TeV energies, the wavelength of the low energy photon field is in the infrared range

Particle acceleration processes

The emission mechanisms reviewed in Section 3.2 require high energy particles, in order to produce high energy photons. The main acceleration mechanisms to produce such population of high energy particles are presented in this section.

Second order Fermi mechanism

This mechanism was proposed by [START_REF] Fermi | On the Origin of the Cosmic Radiation[END_REF] to explain the energies observed in the measured cosmic-ray spectrum. It relies on the deflection of charged particles by magnetized clouds. When entering in the cloud, particles travel within magnetic fields inducing an electric field which can then accelerate particles. The energy increase through this process can be computed in two ways: via the induced electric field acceleration or by referential frame changes [START_REF] Parizot | Accélération des particules : les mécanismes de Fermi[END_REF]. These approaches are equivalent and the reference frame change is chosen for simplicity.

As illustrated in Figure 3.7, a particle in a observer frame entering a cloud at energy E in is seen in the cloud frame with energy E in expressed as:

E in = γE in (1 -β c cos θ in ) , (3.38)
with v c = β c c being the cloud velocity in the observer frame, γ its Lorentz factor and θ in the pitch angle between the particle and cloud trajectories. Then, the particle exits the cloud with an energy E out , seen in the observer frame as E out , which is expressed as:

E out = γE out 1 + β c cos θ out , (3.39)
with θ out the outgoing pitch angle after deflection in the cloud. In the observer frame, E out can then be expressed with respect to E in as:

E out = γ 2 E in (1 -β c cos θ in ) 1 + β c cos θ out , (3.40)
which allows for the computation of the energy gain of the particle:

∆E E = β c cos θ out -cos θ in + β 2 c 1 -cos θ int cos θ out 1 -β 2 c .
(3.41)

To compute the average energy gain ∆E E , one needs to determine the incoming and outgoing angles. Assuming an isotropic distribution of the outgoing angle results in cos θ out = 0. For the incoming angles θ in , the probability of having a collision P (θ in ) is proportional to the relative velocity between the particle and the cloud, giving thus P (θ in ) ∝ v p -v c cos θ in which gives:

cos θ in = 1 -1 P (θ in ) cos θ in d (cos θ in ) 1 -1 P (θ in )d (cos θ in ) = -2v c /3 2v p - β c 3 . (3.42)
where v p is the particle velocity, assumed to be relativistic so that v p c. Finally, the average energy gain per particle is: As a result, this mechanism provides on average a second order energy gain in β c . This leads to a rather slow acceleration process requiring long acceleration time for particle to be able to emit up to TeV energies. In case of cosmic rays interacting with interstellar clouds of typical velocity of tens of km/s and considering a typical distance between two clouds of the order of a parsec, the acceleration time is of the order of billions of years in order to double the particle initial energy. Besides, it strongly relies on the presence of magnetic clouds in the path of the particles and on specific values of the initial pitch angle allowing for the increase of the particle energy. Indeed, even if the average energy gain is positive, not all interactions between particles and clouds induce a gain of energy. The second order Fermi mechanism appears as a stochastic process which can appear in turbulent plasmas.

∆E E = 4 3 β 2 c . ( 3

First order Fermi mechanism

A more efficient acceleration mechanism, based on the same principle, would ideally require only head-on collisions between particles and clouds in order to get cos θ in = -1 and ∆E E ∝ β. The presence of shock waves and the transfer of a particle between the non-shocked and shocked media allows for the obtention of such conditions. In the shock wave referential frame, the shocked medium (index 1) follows the shock wave meaning v 1 = v s , while the non-shocked medium (index 2) moves at a velocity v 2 = vs r , where r is the compression ratio. For a particle in the shocked medium (Figure 3.8, left), the non-shocked medium goes towards the particle at velocity ∆v = v 1 -v 2 = r-1 r v s . Then, when the particle crosses the shock wave and enters the non-shocked medium (Figure 3.8, right), the particle "sees" the shocked medium going towards it at the same velocity ∆v. Consequently, the particle can be pushed back to the shocked medium and do round trips between the two media.

With this effect and the double frame change presented in the previous section, an efficient looped acceleration mechanism is obtained. Most of the time head-on collisions occurs between particles and the shock wave front giving average incoming and outgoing angles cos θ in = -cos θ out = 2 3 β s where β s = ∆v/c. From Equation 3.41, the resulting energy gain per particle is:

∆E E = 4 3 β s . (3.44)
This results in a first order energy gain in β s . For a particle of initial energy E 0 , and doing one round trip between the two media, its final energy E 1 is:

E 1 = E 0   1 + ∆E E   = E 0 4 3 r -1 r v s c , (3.45)
After k round trips, its final energy E k is: (3.46) where u = 4 3 r-1 r vs c . From this equation, the number of cycles k(E) needed to reach the energy E for a particle of initial energy E 0 can be deduced:

E k = E 0 1 + 4 3 r -1 r v s c k = E 0 (1 + u) k ,
k(E) = log E E 0 log (1 + u) . (3.47)
However, this mechanism cannot imply an infinite number of cycles k because particles can escape from the shock. Particles can only escape from the shock when they are in the shocked medium. Indeed, as shown in Figure 3.8, the shock always moves towards the non-shocked medium so the particle will always be caught up by the shock since they are isotropic in their medium frame (v p = 0). On the contrary, in the shocked medium, the particle "sees" the shock moving away, which bring the possibility of escaping the acceleration loop. For an initial particle density n 0 , the total flux of particle crossing the shock is computed considering particles with incoming angle of which π/2 < θ in < π:

Φ 0 = n 0 c 1 0 cos θ in d (cos θ in ) 1 0 d (cos θ in ) = n 0 c 4 . (3.48)
The flux of escaping particles is computed considering the particle drift with respect to the shock which happen at velocity v 2 , giving:

Φ e = n 0 v 2 = n 0 v s /r (3.49)
The escaping probability for a particle is then simply deduced by:

P e = Φ e Φ 0 = n 0 vs r n 0 c 4 = 4v s rc . (3.50)
From an initial number of particles N 0 at energy E 0 , the number of particles accelerated at least at the energy E, i.e. particles remaining in the shock loop after k(E) cycles, is1 :

N (≥ E) = N 0 (1 -P ) k(E) = N 0 (1 -P e ) log(E/E 0 ) log(1+u) (3.51) = N 0 E E 0 log(1-Pe) log(1+u) (3.52) = N 0 E E 0 -r+2 r-1 . (3.53)
Finally, the resulting particle spectrum, n(E), produced by the first order Fermi mechanism can be deduced:

n(E) = dN (≥ E) dE = r + 2 r -1 -1 N 0 E E 0 -r+2 r-1 . (3.54)
Hence, this acceleration mechanism produces particle spectra which follow power law distribution with an index only depending on the compression ratio r. In the case of strong shocks (when the velocity v s is greater than the velocity of sound in the shocked medium) and assuming an ideal gas, the compression ratio is r ≈ 4 which results in a particle spectrum index of -2. This result is important for cosmic rays studies and blazar modeling and justify at least partially that spectra assumed for the radiating particles are typically power law functions with index close to -2. This seems to give a kind of universal explanation for the frequent observation of power law electromagnetic spectra in the case of non-thermal comic sources, as synchrotron emission from a power law particle distribution gives a power law electromagnetic spectrum. For the case of the relativistic shocks, several works based on Monte Carlo simulation [START_REF] Stecker | Blazar γ-Rays, Shock Acceleration, and the Extragalactic Background Light[END_REF][START_REF] Spitkovsky | Particle Acceleration in Relativistic Collisionless Shocks: Fermi Process at Last?[END_REF] show the generation of power law particle distributions with indexes ranging from 1.8 to 4.

Magnetic reconnection

Another efficient mechanism to accelerate particles is the magnetic reconnection. In highly conductive astrophysical plasma, magnetic fields are frozen in the plasma and so carried by the plasma motion. So depending on the velocity field, magnetic field lines of opposite directions can be spliced together, leading to a reconnection mechanism and changing their topological pattern. A simple configuration of magnetic reconnection is given in Figure 3.9.

The evolution of the magnetic field can be deduced from Maxwell equations:

∂ - → B ∂t = - → ∇ × - → v × - → B - - → ∇ × η - → ∇ × - → B (3.55)
where η is the magnetic diffusivity. In the equation, when the magnetic field is frozen in the plasma, the diffusive term, characterized by η, is negligible which defines the frozen-in description. At a stagnation point, where two magnetic fields frozen in two opposite plasma flows connect (Figure 3.9), the velocity field -→ v → -→ 0 and the diffusion term dominates the magnetic field evolution. Thus, the magnetic field energy can be converted to electric and kinetic energy leading to an efficient acceleration site for charged particles.

Magnetic reconnection can arise in magnetic turbulences which can also accelerate particles. A detailed discussion about turbulences, as well as magnetic reconnection can be found in [START_REF] Lazarian | Turbulence, Magnetic Reconnection in Turbulent Fluids and Energetic Particle Acceleration[END_REF] and references therein. The particle spectrum produced follows power law functions with index ranging from -2.5 to -1 for this mechanism [START_REF] Kowal | Particle Acceleration in Turbulence and Weakly Stochastic Reconnection[END_REF][START_REF] Sikora | Comptonization of diffuse ambient radiation by a relativistic jet: The source of gamma rays from blazars?[END_REF]. As for the first order Fermi mechanism, magnetic reconnection produces power law particle distribution with however a different index which favored the production of higher energy particles.

Modeling blazars spectral energy distribution

The origin of blazar emission is usually attributed to a high density region, moving at relativistic velocities along a larger scale jet (as seen in Section 3.1) and filled with magnetic fields and high energy particles. The SED is characterized by two distinct bumps, involving at least two different emission processes. A general consensus on the origin of the low energy bump, ranging from radio to X-ray energies, is that it originates from electrons and positrons synchrotron emission (Alfvén & Herlofson, 1950). For the second one, different models explain the high energy emission, depending on the nature of the particles responsible of the blazar emission.

Leptonic models

Leptonic models are considering electrons and positrons as the dominant particles responsible for blazar emission. The high energy radiation, from hard X-rays to γ-ray energies, is explained by IC interactions, with various photon field targets. For low luminosity blazars, the emission is mainly produced via Synchrotron Self-Compton (SSC) emission: the synchrotron photons are boosted by the same leptons which produced them. An example of leptonic SSC model applied to the blazar Markarian 421 is shown on Figure 3.10a. Different models of SSC emission were developed by [START_REF] Dermer | Model for the High-Energy Emission from Blazars[END_REF], [START_REF] Ghisellini | On the origin of the gamma-ray emission in blazars[END_REF], [START_REF] Bednarek | Testing the homogeneous synchrotron self-Compton model for gamma-ray production in MRK 421[END_REF][START_REF] Tanvir | A 'kilonova' associated with the short-duration γ-ray burst GRB 130603B[END_REF], Katarzyński et al. (2001).

For some blazars, additional IC processes called External IC emission (EIC) are required to reproduce the high flux observed. External source of photons contribute to enhance the total flux compared to the SSC process only. These photon fields can be generated by the accretion disk around the BH, the broad and narrow line regions which are clouds located respectively close and distant to the central BH (inducing a broadening of emission line via Doppler boosting), the dust torus, the jet synchrotron emission or the X-ray corona (see Figure 2.5). A SED considering EIC as well as SSC on the blazar AP Librae is shown Figure 3.11. The EIC emission is usually involved in the emission from FSRQ while a simple SSC scenario is sufficient to model HBL emission thus partially explaining the blazar sequence discussed in Section 3.1.3. More details about EIC process can be found in the models developed by [START_REF] Sikora | Comptonization of diffuse ambient radiation by a relativistic jet: The source of gamma rays from blazars?[END_REF], [START_REF] Ghisellini | Canonical high-power blazars[END_REF], [START_REF] Hervet | Linking radio and gamma-ray emission in Ap Librae[END_REF].

Hadronic and lepto-hadronic models

Hadronic models are considering protons to produce the high energy emission. Hadronic models are mainly based on proton synchrotron radiation and p-γ interactions. On the contrary to leptonic synchrotron, protons can emit synchrotron photons with much higher energies, due to the mass of the proton m p /m e ≈ 1836 but it requires a much higher magnetic field strength in order for protons to emit synchrotron radiation. The p-γ interaction involves synchrotron photons coming from electrons which provide an external photon field for EIC interactions. These interactions also produce muons which can emit synchrotron photons contributing to high energy emission. The p-p interaction is usually less important as it requires a very high proton density which is not common in hadronic models [START_REF] Reynoso | A lepto-hadronic model for highenergy emission from FR I radiogalaxies[END_REF]. An example of hadronic SED modeling is shown in Figure 3.10b. A few number of hadronic models can be found in [START_REF] Mücke | First combined studies on Lorentz Invariance Violation from observations of astrophysical sources[END_REF], [START_REF] Cerruti | A hadronic origin for ultra-high-frequency-peaked BL Lac objects[END_REF].

One can also consider both leptons and protons to produce the high energy photons in leptohadronic models. These models have the advantage of requiring lower magnetic field strength than pure hadronic models making them more realistic. It also naturally brings a mechanism for the injection of electrons for leptonic processes, from the protons interactions. For examples of such a model called lepto-hadronic model, see [START_REF] Cerruti | A hadronic origin for ultra-high-frequency-peaked BL Lac objects[END_REF], [START_REF] Ansoldi | The Blazar TXS 0506+056 Associated with a High-energy Neutrino: Insights into Extragalactic Jets and Cosmic-Ray Acceleration[END_REF], [START_REF] Cerruti | Lepto-hadronic single-zone models for the electromagnetic and neutrino emission of TXS 0506+056[END_REF], [START_REF] Keivani | A Multimessenger Picture of the Flaring Blazar TXS 0506+056: Implications for High-energy Neutrino Emission and Cosmic-Ray Acceleration[END_REF].

As a remark, hadronic models rely on several processes in order to produces γ-ray emission. The addition of these different contribution imprint a specific signatures in the energy spectrum of blazar at TeV energies which can be used as indication for hadronic processes. Up to now, current IACT are not sensitive enough to measure such signatures but CTA will provide such a performance (Zech et al., 2017). Besides, neutrinos can also be used to reveal the presence of hadronic processess, as they are not expected from leptonic models. A recent observation from IceCube Collaboration & al. (2018) shows the first coincidence between an astrophysical neutrinos, observed by the neutrino observatory IceCube, and a blazar flare observed by MAGIC. This result provides the first hints of hadronic emission processes for blazar emission.

Summary

The specificities of blazars, the mechanisms at play for particle acceleration and high energy photons production were presented in this section as well as different models based on the nature of the particle population responsible for the emission. In this thesis, leptonic models will be favored over the hadronic ones since fast variability seems easier to produce in the frame of leptonic models, based on the simple fact of the mass difference between electron and proton leading to larger radiative time scales for protons [START_REF] Catanese | Very High Energy Gamma-Ray Astronomy[END_REF]. Besides, pure hadronic models requires higher magnetic field values in average to emit enough proton synchrotron emission, and lepto-hadronic models, due to the poor conversion efficiency of proton into γ-rays (via pions decays), need a lot of protons and so very high power in the jet, above the Eddington limit. As a result, leptonic models are more favored compared to hadronic ones because of more reasonable physical parameter values.

In the following, a specific SSC scenario is presented and developed to build a time-dependent model. This model will be described and used to investigate and characterize intrinsic time delay.

Searching for energy dependent time delays in the arrival time of γ-ray photons from remote cosmic sources is a way to study LIV effects and investigate QG models. Such measurements rely on the hypothesis that LIV is the only effect producing energy dependent time delays. However, for a single source, delays can also appear at the emission of these photons. Such intrinsic effects are already observed with GRB and reproduced by models describing GRB emission. A list of GRB models can be found in Dai et al. (2017) and references therein. In the case of blazar flares, intrinsic effects are not as important because no significant time delays has been found up to now for LIV studies. Nevertheless, to study intrinsic time delays in blazar flare within the LIV framework, time dependent modeling is investigated in this thesis, to help constraining the QG energy scale in case of future significant detections of time delays.

In this chapter, after a short review of different intrinsic delay origins, a minimal timedependent leptonic model to generate blazar flares is described, explicitly developed to study intrinsic time delays at γ-ray energies, above 1 MeV. This model relies on the description of the electron population through a transfer equation and an instantaneous SSC model to compute the emission coming from these electrons. To determine the time delay, two methods are compared and the more robust of them is chosen and applied to provide the first results on intrinsic effect from the model.

Time dependent models

Several time dependent emission models of non thermal emission in AGN have been proposed in the literature, as for instance [START_REF] Blandford | Relativistic jets as compact radio sources[END_REF], [START_REF] Marscher | Models for high-frequency radio outbursts in extragalactic sources, with application to the early 1983 millimeter-to-infrared flare of 3C 273[END_REF], [START_REF] Celotti | A model for the spectral variability of BL Lacertae objects at high frequencies[END_REF], [START_REF] Boettcher | γ-ray emission and spectral evolution of pair plasmas in AGN jets. I. General theory and a prediction for the GeV -TeV emission from ultrarelativistic jets[END_REF], [START_REF] Katarzyński | The multifrequency variability of Mrk 421[END_REF], Joshi & Böttcher (2011), [START_REF] Lewis | Time-dependent Electron Acceleration in Blazar Transients: X-Ray Time Lags and Spectral Formation[END_REF]. They usually include frequency dependent time delays. However, there is no complete dedicated study of such intrinsic time delays available for TeV blazar flares. Time delays can be classified in two categories depending on their origins: "microscopic" and "macroscopic" delays. The microscopic delays originate from the temporal evolution of particles responsible for the emission. They mostly come from acceleration and cooling processes for particles. The macroscopic delays origin comes from the evolution of the macroscopic parameters such as the magnetic field strength or the global geometry of the blob responsible of the emission and the jet. Following this classification, a few examples of published models are presented which appear to produce time delays in their light curves results.

"Macroscopic" delays

In the models presented here, the discussion only deals with non-thermal high energy emission of blazars. A model proposed by Bednarek & Wagner (2008) suggests that the time delay originate from a continuous increase of the global Lorentz factor of the emitting zone. They consider that emission at energies E 1 and E 2 come from the same population of electrons. They parameterize the acceleration of the blob via an increase of its Lorentz factor, Γ blob = A×t, where A is the acceleration parameter expressed in (Γ b ) s -1 . This increase of the blob Lorentz factor leads to an increase of the Doppler boosting effect. So, at the earliest stage, when the Lorentz factor is small, electrons emit photons at energies E 1 . At later times, when the acceleration have reached a high enough Lorentz factor, the emission is observed at energy E 2 and thus creating a delay for the observer between energies E 1 and E 2 . In this model the time delay comes from the intensification of the Doppler boosting effect during the evolution of the emission zone. It can be also interpreted as the lower energies being emitted at a different location than the higher 4.1 Time dependent models energies, since the blob travels a large distance before reaching high Lorentz factor values. In this model of accelerating blob, high energy photons always arrive later than lower energy ones.

Another model from [START_REF] Sokolov | Synchrotron Self-Compton Model for Rapid Nonthermal Flares in Blazars with Frequency-dependent Time Lags[END_REF] describes the emission with a complex jet geometry and shock acceleration. The model considers the collision between a stationary jet and a relativistic shock wave which results in the production of a flare. The collision leads to the creation of a reverse and a forward shock defining several regions inside the jet depending on if one of the shock already crossed the region or not. Each shock provides particle acceleration. Finally, they also consider the photon internal travel time in the jet, which depends, for an observer, on the viewing angle in the plasma reference frame. In total the sum of the SSC emission coming from these different regions leads to complex light curve shapes. Focusing on emission coming from the acceleration of the forward shock, the IC emission shows time delays between light curves of different energies (Figure 4.1). This delay is explained by the cooling time of electrons. Electrons of high energies with short cooling time are distributed in a thin region while less energetic electrons are spread in a larger region, because of a longer cooling time. Thus, the shock leaves the high energy electrons region shutting down the high energy light curves earlier, while less energetic electrons are still accelerated by the shock, inducing a time delay between their respective emission. In addition, depending on the observer viewing angle, the photon internal travel time can become significant and induce light curve peaking after the shock left the region, which modify the observed time delay.

"Microscopic" delays

Whereas macroscopic delays have been explicitly studied as presented above, microscopic delays appear in many models but are barely discussed at TeV energies. In the X-ray range, [START_REF] Lewis | Time-dependent Electron Acceleration in Blazar Transients: X-Ray Time Lags and Spectral Formation[END_REF] produce delayed light curves with a transfer equation describing the evolution of the electron population. The flare is produced on top of a stationary emission described by a second transfer equation. A shock injects relativistic electrons responsible for the flares at a level of 10% of the steady emission. Their model takes into account first and second order Fermi accelerations, synchrotron and adiabatic energy losses, particle injection and escape. From these consideration, the authors produce light curves presented in Figure 4.2 which model the flare from the blazar Markarian 421 in 1998 [START_REF] Catanese | Hour-Scale Multiwavelength Variability in Markarian 421[END_REF]. The resulting light curves at two different energies (1.05 keV and 6 keV) present a time delay where the low energy light curve peaks at later time than the high energy one. This delay between the two light curves arises because of their difference of flux increase and decrease phase. In their work, the authors describe the time delay using Fourier transform which gives the evolution of the time delay with time between two light curves. This method will be discussed later on in Section 5.3.3. This model only considers synchrotron emission and thus is not in the relevant energy range considered for this work. Indeed, [START_REF] Lewis | Time-dependent Electron Acceleration in Blazar Transients: X-Ray Time Lags and Spectral Formation[END_REF] entirely neglect the IC losses because they cannot be modeled using non-linear transfer equation without approximation. For this work, the IC processes are investigated on a similar electron evolution via the elaboration of a timedependent model. This model describes the electron evolution using a transfer equation which is then used to produce the SED evolution allowing for the deduction of γ-ray light curves.

Electron spectrum time evolution

General equation

The general transfer equation in a plasma for isotropic particles of a kind i at a time t, energy γ and at a point r, is given by Ginzburg & Syrovatskii (1969). This equation describes the temporal evolution and takes into account energy losses, acceleration, particles interaction and spatial diffusion. It is expressed as:

∂N i (t, γ, r) ∂t -div D i (t, γ, r)∇N i (t, γ, r) A + ∂ ∂γ b i (t, γ)N i (t, γ, r) B - 1 2
∂ 2 ∂γ 2 d i (t, γ)N i (t, γ, r) C = Q i (γ, t, r) D -p i (γ)N i (t, γ) E + k γ P k i (γ , γ)N k (γ , t)dγ F , ( 4.1)
where N i (t, γ, r) is the density of particle i. The different terms are detailed and discussed in the following in the context of the description adopted here for the particle spectrum evolution involved in blazar flares.

The second term on the left side (A) represents the diffusion of particles, with the diffusion coefficient D i (t, γ, r). For blazar flares, the particle distribution function is usually assumed homogeneous, leading to ∇N (t, γ) = 0 inside the blob. At the edge, where there is a discontinuity between the blob and the surround medium, this term is neglected as it represents a small fraction of the total particle population. Thus this term is not considered for the blazar flare model.

The third and fourth terms (B and C) on the left side correspond to energy variations of the particle spectrum, which corresponds mostly in acceleration process and radiative cooling but can also include adiabatic losses. The third one (B) corresponds to the systematic energy variations, with the coefficient b i (t, γ) being the sum of mean energy increment of each process j happening in the source:

b i (γ, t) = j dγ dt j = j 1 mc 2 dE dt j , ( 4.2)
For synchrotron and IC processes, this coefficient was already expressed in Section 3.2. The fourth term (C) complements the third one, accounting for fluctuations of the systematic energy variations. The coefficient d i (t, γ) can be expressed in the general case as:

d i (t, γ) = d dt (∆γ) 2 , (4.3)
where (∆γ) 2 is the mean squared energy gain. This term can represent a particle acceleration process like the second order Fermi acceleration or energy losses like a second order IC process (photons which interact a second time with electrons by IC process). Theses processes are not considered in the blazar flare model which focuses on the systematic energy variations processes.

The three terms on the right side (D, E and F) represent injection and loss of particles. In the term D, Q i (γ, t) is a source term dealing with particle injections, corresponding to the number of particles supplied by unit of time dt in the energy range [γ, γ + dγ] in the particle population N i (t, γ). For the term E, p i (γ)N i (t, γ) is related to the "catastrophic" losses of particles, either by escaping the medium or by "catastrophic" interaction leading to the disappearance of particles and is characterized by the probability p i (γ). In case of disappearance of the particle due to an interaction, p i (γ) is expressed as:

p i (γ) = 1 T i = nvσ i , ( 4.4)
where T i and σ i are respectively the mean life time for the particle and the cross section of the "catastrophic" interaction, n is the density of particles and v the velocity of particles. Finally, the last term (F) corresponds to the creation of secondary particles from "catastrophic" interactions of primary particles of kind k with P k i (γ , γ) being the probability per unit of time of the creation of particle i of energy γ in a "catastrophic" interaction of a particle of kind k with an energy γ .

Simplified equation

In the Equation 4.1, the most important term for modeling a blazar flare is the systematic energy variation (term B) represented by the b i (t, γ) coefficient. It allows for the generation of flares via acceleration processes and energy losses like for instance emission of photons. Particle injection processes (term D and F) can also generate flares but in order to obtain a minimal scenario for this model, with the smallest number of parameter possible, they are not considered as well as statistical energy variations (term C) and particle escape (term E). This allows for the highlight of the most dominant and fundamental processes needed to generate a flare. Considering an homogeneous population of electrons, the Equation 4.1 becomes:

∂N i (t, γ) ∂t + ∂ ∂γ b i (t, γ)N i (t, γ) = 0. ( 4
b syn (t, γ) = γ 2 4σ T c 3m e c U B (t), (4.6) b IC (t, γ) = γ 2 4σ T c 3m e c U syn (t), (4.7)
where U syn (t) is the total synchrotron field energy density and U B defined by Equation 3.16. Since the synchrotron energy density is the energy lost by electrons via synchrotron emission, U syn (t) is expressed as:

U syn (t) = ∞ 0 b syn (t, γ)N e (t, γ)dγ. (4.8)
The expression of U syn (t) complicates the resolution of the Equation 4.5, preventing a simple determination of an analytical solution without an approximation. To be able to use the exact expression of U syn (t), numerical resolution of the equation is required. Nevertheless, with a judicious choice of approximations, [START_REF] Katarzyński | The multifrequency variability of Mrk 421[END_REF] determined an analytic solution to Equation 4.5. The complete formula of the analytic solution is detailed in [START_REF] Katarzyński | The multifrequency variability of Mrk 421[END_REF]. This solution takes into account the SSC emission, adiabatic expansion and a parameterized acceleration. The main approximation to be able to derive an analytical solution relies on a parameterization of the synchrotron photon density in terms of the magnetic energy:

U syn (t) = U B (t) η , (4.9)
where η is the ratio between the synchrotron and the IC emission power. This parameterization implies that the energy lost by IC emission evolves as the magnetic energy density. Thus, the addition of synchrotron and IC losses coefficient can be rewritten as a single coefficient:

b rad (t, γ) = γ 2 4σ T c 3m e c U B (t) 1 + 1 η = γ 2 C rad (t). ( 4

.10)

A second hypothesis is that the physical parameters B(t) and R(t) respectively decrease and increase as power laws of time:

B(t) = B 0 t 0 t m b , (4.11) R(t) = R 0 t 0 t -mr , ( 4.12)
where B 0 and R 0 are the initial magnetic field strength and initial radius in the source frame, m b and m r are respectively their evolution index and t 0 is the evolution time scale in the blob frame, related to adiabatic expansion. This time t 0 is computed as:

t 0 = R 0 V exp , (4.13)
which corresponds to the time needed to cross a blob of size R 0 with a wave velocity of the order of the velocity of sound in a relativistic plasma: V exp = c/ √ 3. As a result, since the total number of particles is conserved (no particle injection or loss), the adiabatic expansion coefficient is:

b adi (t, γ) = γ m r t = γC adi (t). (4.14)
The time t 0 , and thus the radius R 0 govern the temporal evolution of the other parameters and the whole evolution of the model. In particular, t 0 is related to the duration of the flare in the source frame and so to a factor δ for the observed time.

The acceleration coefficient b acc allows for the launch of the flare. In the model, the acceleration process is not specified because of the variety of process able to accelerate electrons, but is rather parameterized with an energy independent acceleration rate decreasing with time. In this regard, it is expressed as:

b acc (t, γ) = γA 0 t 0 t ma = γC acc (t), (4.15)
with A 0 (in s -1 ) is the initial acceleration amplitude at t = t 0 and m a its evolution index. The value of A 0 , corresponding to the inverse of the characteristic acceleration time, mainly determines the maximum energy reached by electrons during the flare while the index m a sets the duration of the acceleration phase, before the energy losses overcome the acceleration. The energy dependency of b acc (t, γ) is the same as b adi (t, γ) which result in an attenuated acceleration when adiabatic expansion is considered. The resulting differential equation with all terms is written as:

∂N e (γ, t) ∂t = ∂ ∂γ C rad (t)γ 2 -C acc (t) -C adi (t) γ N e (γ, t) . (4.16)
However, due to adiabatic expansion, one has to be careful with the electron density evolution. Indeed, in that case, the radius of the emitting blob R(t) increases with time (Equation 4.12) and so does its volume. This implies that the density needs to be corrected at each time of the radius increase. Thus, in case of an adiabatic expansion, the effective electron density N * e (γ, t) becomes:

N * e (γ, t) = N e (γ, t)

t 0 t mr . (4.17)
Finally, the last ingredient of the analytic solution is the initial electron spectrum N e (t 0 , γ). It is taken as a power law function with a high energy cut-off:

N * e (t 0 , γ) = K 0 γ -n   1 - γ γ c,0 n+2   , (4.18)
with K 0 the initial electron density in cm -3 , n the power law index and γ c,0 the initial energy cut-off of the electron spectrum.

Blob in jet emission model

From the analytic solution of the electron spectrum evolution, a set of electron spectra N * e (t i , γ) is obtained corresponding to different times t i . Each electron spectrum of the set is used to deduce the corresponding SED at time t i by applying an instantaneous emission model on each electron spectrum.

Presentation of the emission model

The emission model computing the SSC emission is taken from Katarzyński et al. (2001). This model considers the synchrotron emission, taking into account the synchrotron self-absorption and the associated IC emission considering the γ -γ absorption inside the blob. The EBL absorption during the propagation of photon from the blob to the Earth is also considered from the model of Franceschini et al. (2008). Moreover, Doppler boosting is applied to obtain the observed energy and photon flux.

The model considers an homogeneous spherical blob of radius R b filled with a magnetic field of strength B. This blob has a bulk relativistic motion toward the observer characterized by its Doppler factor δ b . The magnetic field strength and the blob radius at time t i are taken as B = B(t i ) and R b = R(t i ) for each electron spectrum.

Additional simplifying assumptions

The model uses several computation approximations widely used in the literature. For the synchrotron radiation computation, an integral over all the pitch angles between the particle and the magnetic field lines is required. This can induce complex computation requiring a lot of time, especially when multiple computations of synchrotron emission have to be done for each electron spectrum describing the time evolution. To ease the computation, a simple analytical 4.3 Blob in jet emission model formula approximating the integration over the pitch angle is used:

1 2 π θ=0 sin θ ∞ 2t sin θ K 5 3 (x)dxdθ ≈ C 1 t C 2 e -C 3 t , ( 4.19)
where t = ν/(3γ 2 ν B ) and ν B = eB/(2πmc). The coefficient values are: C 1 = 0.78, C 2 = 0.25 and C 3 = 2.175. This integration formula is required twice, for the synchrotron computation. It is used once for the synchrotron emissivity, after injecting Equation 3.8 in Equation 3.11:

j(ν) = e 3 B 8πmc 2 π θ=0 sin θ γmax γ min N e (γ)P (γ, θ, ν)dγdθ, (4.20) j(ν) = e 3 B 8πmc 2 ν ν B π θ=0 sin θ tmax t min N e ν 3tν B t -1/2 ∞ 2t sin θ K 5/3 (x)dxdγdθ. (4.21)
Similarly, it is used for the synchrotron self-absorption, in Equation 3.12 which gives:

k(ν) = √ 3e 3 B 24πm 2 c 2 ν B 1 ν γmax γmin d dγ N e (γ) γ 2 π θ=0 sin θ ∞ 2t sin θ dθdγ. (4.22)
For the computation of IC emission, the number density of synchrotron photons n(ν, r) is required to compute its emissivity and the γ-γ absorption effect. Assuming a uniform synchrotron intensity in blob n(ν, r) = n(ν), it is expressed as:

n(ν) = 4π hcE j s (ν) k s (ν) 1 -exp -k(ν)R b , (4.23)
where j s (ν) and k s (ν) are the synchrotron emissivity and self-absorption coefficient defined in Section 3.2.1. In reality, it is known that the photon field is not strictly uniform and decreases from the center to the edges of the blob [START_REF] Gould | Compton and synchrotron processes in spherically-symmetric non-thermal sources[END_REF]. To account for it, a correction factor 3/4 is added to U syn (ν), allowing for the scaling down of the synchrotron intensity. An additional approximation concerns the γ -γ absorption of IC photons on synchrotron photons. Using the approximated formula from Coppi & [START_REF] Cleveland | Reaction rates and energy distributions for elementary processes in relativistic pair plasmas[END_REF] and after integration over the synchrotron density and the cross section of the interactions, an optical depth corresponding to the γ -γ process is deduced as:

τ γ-γ (x) = 0.2σ T 1 x n 1 x R b . (4.24)
where x = hν/m e c 2 is the adimensional photon energy.

The emission model considers a homogeneous blob filled with electrons responsible of an isotropic photon emission. During the emission, photons are more likely emitted close the center of the blob rather than through the entire emission zone. Under such circumstances, local delays induced by the propagation times inside the emission zone cancel on average for an observer at the infinite and thus are not considered here.

Domain of validity: η parameter

The parameter η used to parameterize the total synchrotron energy U syn (t) in function of U B (t) can be computed at the emission level. Indeed, from the computation of synchrotron emission, I rad (ν), defined in Equation 3.13, U syn (t i ) can be computed with the electron spectrum N e (t i , γ):

U syn (t i ) = 4π c ∞ 0 I rad (ν)dν. (4.25)
One can compute η(t) = U B (t)/U syn (t) at the emission stage to compare it to the constant value η assumed for the electron spectrum evolution. Figure 4.3 shows the typical evolution of η(t) from the model, which reveals an important variation, spanned over several orders of magnitude during the flare. The decrease of η(t) is the result of the acceleration of electrons which enhances the total synchrotron energy density U syn (t) while at the same time U B (t) decreases due to the magnetic field evolution. After the flare reaches its maximum, electrons starts to lose their energy, inducing a diminution of the synchrotron energy density leading to an increase of η(t).

As η is assumed to be constant for the computation of the electron spectra, the time dependent model only allows for the investigation of cases where the variation of η(t) can be neglected. In the analytic solution, the parameter appears as η -1 (Equation ??) and if η(t) 1, it implies that its variations does not affect the electron spectrum evolution. All the scenarios investigated in this work satisfy this condition η(t) 1 so that its variation can be neglected. This corresponds to a blazar SED for which the IC peak remains smaller than the synchrotron one. This is usually the case for HBL sources (Figure 3 The decrease of η(t) corresponds to the fast acceleration of electrons. When electrons reach a large energy, the radiative cooling becomes more important than the electron acceleration leading to an increasing η(t).

From the time dependent model to time delays 4.4.1 Light curve determination

A typical case is investigated at first as illustration of the model outputs. The set of parameters considered for this reference case are listed in Table 4.1. Parameter values related to SSC emission are taken from typical values obtained for blazars (Lenain, 2009). The parameters related to the time evolution are deduced to obtain a flare lasting few hours. From these parameters, the electron spectrum evolution is obtained, as shown on Figure 4.4a with a few electron spectra. The electron population, starting from a low state, gains energy via the acceleration processes until the the maximum of the evolution is reached. Then, the radiative cooling starts to drives the evolution of the most energetic electrons as it become more important than acceleration due to the energy dependency of the electron cooling time scale t cool ∝ γ -1 . In the mean time, adiabatic expansion affects all energies decreasing the overall density via the increase of the radius R(t) (Equation 4.17).

The SEDs are obtained applying the instantaneous emission model on each electron spectrum of the evolution. They are presented in Figure 4.4b. The SED reflects the evolution of the electron spectrum, with a fast acceleration phase where the overall flux increases and the IC bump becomes more important, followed by the cooling phase with a slower decrease of flux.

Using these SEDs, the corresponding light curves in different energy bands are computed. The light curve in an energy range [E 1 , E 2 ] is obtained by integrating the flux in this energy range for all time t i . The integral is done using a Gauss-Legendre algorithm. 30 light curves are computed, spanned in the γ-ray energy range from 1 MeV to 50 TeV with a constant energy bin width in logarithm scale. Then only the light curves with a flux of at least 2 × 10 -15 cm -2 .s -1 are kept to avoid computational precision errors and consider the one with a flux measurable with current instruments. A few light curves are presented in Figure 4.5, showing a clear shift in time of their maxima. In order to measure precisely this time delay, a study is carried out to identify the appropriate method to use.

Parameter

Time delay measurement methods

The comparison of two simple methods to measure time delays is done to evaluate their respective accuracy: a simple Cross Correlation Function (CCF) and the so called "Peak Position Method" (PPM). For both methods, the time delay is defined as the delay between the lowest energy light curve compared to the other ones. The CCF method is the one published by [START_REF] White | Comments on cross-correlation methodology in variability studies of active galactic nuclei[END_REF]. It compares the full light curve while the PPM is a simple method which measures the delay as the time difference between the maxima of the two light curves considered. To test the two methods, simulations are done to reproduce as close as possible the normalized light curves from the model.

The light curves from the model have a fast flux increase followed by a slower decrease. Consequently, the shape used to reproduce them is an asymmetric Gaussian function. Besides, the model light curves width decreases when the energy increases. This effect acts significantly only on the decreasing flux phase of the flare. As a result, the right width of the simulated light curve is decreased as the energy increases. The resulting function reproducing the model light curves is: (4.26) where µ i = µ 0 + ∆t i is the position of the pulse i with a delay ∆t i and σ 1 , σ 2 are the left and right widths, the latter depending on the light curve number with ∆σ as the width evolution parameter. Each index i corresponds to a particular energy range. Finally, a progressive time delay is injected in the simulated light curve, with a parameter δt, with respect to the light curve for which i = 0. The delay is injected to the light curve of number i as:

F i (t) =          exp -(t-µ i ) 2 2σ 2 1 if t ≤ µ, exp - (t-µ) 2 i 2(σ 2 +i * ∆σ) 2 if t > µ,
∆t i = δt × i.
(4.27)

The two methods are applied to the simulated lagged light curve.The results of the time delay reconstruction with the CCF and the PPM are presented in Figure 4.6. For the two methods, the errors on the energy correspond to the energy range of each light curve and the ones on the time delay depends on the time step in the light curve. The test shows that the CCF is not able to reconstruct the injected time delay while the PPM obviously is. To investigate the discrepancy between the two methods, a second test is done. The width variations is removed and gives the results shown in Figure 4.7. In that case, both methods reconstruct correctly the injected time delay, indicating that the simple CCF method is not able to reconstruct a lag between light curves of changing widths. As a result, the PPM is preferred to reconstruct the time delay for this work.

Time delay from the reference scenario

Using the PPM, the time delay is determined on the light curves presented in Figure 4.5. The precision of PPM on the time delay reconstruction is limited at the level of the time step used for the model which is dt = 5 s. The lowest energy light curve, in the range 1 -1.8 MeV, is used as a reference, and the time delay is always measured with respect to this reference. A positive time delay corresponds to a light curve peaking later than the reference one. Figure 4.8 presents the results of the PPM and shows an increasing evolution of the delay with respect to the energy. This energy dependent time delay is understood by looking at the temporal evolution of the electron spectrum and the SED shown in Figure 4.4. The electron spectrum energy cut-off γ c (t) reaches its maximum energy γ max at the time t = 4440 s, while the density at lower γ values has already decreased with respect to earlier time. This translates in the SED evolution where the IC bump reaches the highest flux value at t = 3555 s whereas the maximum energy of the emission is reached at the same time the electron spectrum reaches γ max , i.e. t = 4440 s. This behavior is the result of the adiabatic expansion and the dilution of the electron density inside the expanding emission zone which shut down the lowest energy flares even if the electrons has not reached their maximum energy. In the mean time, γ c (t) is still increasing towards γ max for electrons inducing later high energy emission and increasing time delay with the energy.

Energy [TeV]

Conclusion

This chapter presented several blazar models which can produce different types of time delay depending on their origins: from macroscopic properties of the source or micro-physics of the emission processes. This work focused on the microscopic time delay origin. It investigated and described the evolution of the high energy electron population in order to produce a γray flare. The general transfer equation describing such an evolution was presented. Several reasonable assumptions were done and resulted in a simplified equation with an analytic solution which allows for the production of a "minimal flare scenario". This scenario takes into account only dominant and unavoidable processes to generate a flare for the electrons. This choice was motivated to limit the number of free parameters for the time dependent model in order to simplify the interpretation of the time delay origin.

For this minimal blazar flare model, the electron population evolution includes namely a generic acceleration, radiative cooling and adiabatic expansion processes in addition to a magnetic field decreasing with time. The resulting radiation computed with the instantaneous Synchrotron-Self-Compton emission model allows for the obtention of the high energy flux temporal evolution coming from the electron population. By using a Peak Position Method on the light curves produced by the model, the theoretical energy dependent intrinsic time delay an be found and estimated During this thesis, some attempts on a full numerical resolution for the transfer equations were done in order to extend the study to more general time dependent models for blazar flares. However, such investigation was found to be too ambitious for the time devoted to the present work. Further developments will allow for the confirmation that the analytical solution used in this thesis provides a first relevant global characterization of the intrinsic time delays in blazar flares and also the exploration of alternative and more complex scenarios.

The next chapter is devoted to the study of the time delays. An evaluation of the impact of each parameter on the delay is investigated to better understand its origin and to look into other possible evolution of the time delay with respect to the energy. This study is expected to provide constraints on intrinsic time delays. It should contribute to better modeling of blazar flare and help the search for LIV signatures in the arrival time of γ-ray photons.

The previous chapter presents the time dependent model, based on the description of the electron population evolution. From a transfer equation taking into account acceleration, cooling and adiabatic processes, the evolution of the electron spectrum allows for the deduction of the evolution of SED and light curves via a simple SSC emission model. With the reference set of parameters for the model, the resulting light curves in different energy bands show an energy dependent time delay. A time delay in the light curves at γ-ray energies can have strong implications, both for blazar modeling and LIV search. A study of intrinsic time delays induced by the models with different sets of parameters is essential to be better understood them.

In this chapter, the influence of the model parameters on the intrinsic time delay is investigated. Starting from the reference set of parameters, each parameter is varied one by one and the resulting variations of the energy dependent time delays are interpreted to understand its influence. Finally, the presence of energy dependent time delays is discussed in terms of blazar modeling as well as in the context of search for a LIV signature with a dedicated study of the intrinsic delay focusing on the energy range of IACT such as the H.E.S.S. experiment.

On the origin of intrinsic time delay

Generally speaking, the time delay can be interpreted as the addition of several effects inducing spectral evolution impacting the light curves. The main effects involved in the present study are the electron spectrum evolution with acceleration, adiabatic and cooling processes as well as the magnetic field decrease. The electrons acceleration is the only effect inducing a flux increase effect while radiative cooling, magnetic field evolution and adiabatic effect results in a flux decrease. Each effect does not have the same strength depending on the energy, thus leading to different transition times when flux decrease effects start to dominate over the acceleration process.

Scenario without adiabatic expansion

The first energy dependent time delay obtained with the model in the previous chapter is fully driven by the adiabatic expansion. To evaluate the influence of the other parameters, the adiabatic expansion is removed from the model. This corresponds to a scenario where the blob is entirely confined. The resulting transfer equation for electrons is:

∂N e (γ, t) ∂t = ∂ ∂γ C rad (t)γ 2 -C acc (t)γ N e (γ, t) . (5.1)
Besides, removing adiabatic expansion corresponds to a scenario with a constant blob radius: R(t) = R 0 . The evolution time t 0 , related to adiabatic expansion and required for the magnetic field and acceleration evolution, is kept as defined in Equation 4.13. This assumption is reasonable, since t 0 is related to the time needed for a perturbation to propagate to the entire emission zone.

Removing adiabatic losses results in a more efficient acceleration, leading to a more energetic electron spectrum and SED during the evolution (Figures 5.1 and 5.2). The SED present stronger emission compared to the scenario with adiabatic expansion (Figure 5.2). The light curves, shown in Figure 5.3, become wider and peak at later times, due to a longer acceleration phase and a reduced electron energy loss compare to the ones without adiabatic expansion (Figure 4.5). The resulting time delay for this scenario is presented in Figure 5.4 together with the previous scenario with adiabatic expansion. The resulting energy dependent time delay shows a decreasing phase at low energies as well as a lower time delay for the highest ones. To understand this new behavior, the light curve evolution is studied in detail.

Light curve evolution

The energy dependent time delay comes directly from the evolution of light curves since their maxima are used to determine the delay. This evolution is composed of two phases, namely a flux increase phase powered by the acceleration of electrons, and a flux decrease phase resulting from the energy losses of electrons and the decrease of B(t). To determine when the flare maximum occurs, one needs to take a look at when the electron cooling starts to dominate the electron evolution, i.e. when the electron spectrum reaches its maximum state. t max is defined as the time when maximum energy cut-off value γ max is reached. This time does not correspond to the time when light curves reach their maxima, but should happen close to them.

To determine t max , each electron spectrum N (t i , γ) is fitted with a power law function with a high energy cut-off, similar to the initial state defined by Equation 4.18. The reconstructed high energy cut-off values γ c (t i ) for each time are shown in Figure 5.5a. The time t max is determined by the time t i when γ c (t i ) is at maximum, noted γ max = γ c (t max ). As illustrated by Figure 5.5b, light curves peak before the time t max . This indicate that the decay of light curves happen while the electron are still accelerating and have not reached γ max yet. Thus, the radiative cooling cannot explain alone the decay of light curves and the magnetic field decrease is the only remaining flux deacreasing effect able to explaining it.

Time delay energy evolution

The information given by t max , about the magnetic field decrease inducing the light curve decay, is a crucial information to understand the energy dependent time delay. Light curves peaking before the time t max indicates that the electron spectrum is still under acceleration around the maximum of the flare. Thus, at E 10 -3 TeV the increasing time delay in Figure 5.4 is easily explained by the time needed by electrons to reach a sufficient γ c (t i ) to produce the highest energy photons.

For E 10 -3 TeV, the electrons emitting such energies reach quickly a sufficient energy, in opposition to more energetic electrons. Thus, the combination of the decrease of the acceleration power and of the magnetic field induces the light curve decay, inducing flares at E 10 -3 TeV earlier than the more energetic ones. The decrease of the delay is explained by the relative cooling effect being stronger for electrons emitting photons at a few GeV than for the ones emitting photons at MeV. This cooling power, which is not driving the electron evolution when the light curves peak, is sufficient enough to weaken the acceleration depending on the electron energy and induce earlier GeV flares compared to MeV ones.

Model parameters study

Until now, the energy dependent time delay is obtained from the reference case defined by the set of parameters shown in Table 4.1. This section investigates the influence of those parameters on the energy dependent time delay. Starting from the reference set, each parameter is varied one by one to assess its individual influence on the delay. For each case, the evolution of η(t) is checked in order to ensure the current is valid (η 1), which can limit the range for the variation of some parameters.

Initial magnetic field strength

The initial magnetic field strength is defined in the model by Equation 4.11 as B(t 0 ) = B 0 . The magnetic field is one of the main component of the radiative cooling for synchrotron losses and its decrease contributes to the low energy regime of the energy dependent time delay for the reference case without adiabatic expansion. Thus it makes B 0 one of the main parameter acting on the intrinsic time delay.

Several values of B 0 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 as a function of energy is shown in Figure 5.6. The variation of the initial magnetic field value induces large variations of the time delay and a significant change of its behavior at E 10 -3 TeV with a transition from an increasing to a decreasing phase when B 0 increases. Within this transition, some cases show a very small or even no time delay in this energy range. However, the delay behavior at E 10 -3 TeV does not change whereas the delay value slightly decreases as B 0 increases.

To understand this new decreasing phase of the energy dependent time delay, the case with B 0 = 90 mG is taken as an example. The light curves associated with this case are presented in Figure 5.7 which shows a time t max happening before the flare maxima in opposition to the reference case. This implies that the electron evolution, around light curve maxima, starts to be driven by radiative cooling over the acceleration for the more energetic electrons. The electron radiative cooling is an energy dependent process, with a cooling time scale inversely proportional to γ. The most energetic electrons within this case, close to γ max ≈ 10 6 , have a much smaller cooling time scale than less energetic electrons. Hence, the time delay at energies E > 10 -3 TeV is explained by the radiative cooling of electrons. The highest energy light curves peak first (negative ∆t) because the electrons emitting the highest energy photons lose their energy much faster then less energetic electrons. For the delay at E < 10 -3 TeV, the cooling time scale of electrons emitting these energies is larger than the total duration of the flare. Hence, the delay in this energy range behave as in the acceleration driven regime with a decreasing delay with respect to energy due to the magnetic field decrease and the relative cooling power.

Energy [TeV]

Time delay regimes

From all the cases in Figure 5.6, the time delay at E 10 -3 TeV is similar and originates from the same mechanism. Two different regimes emerge from the variation of B 0 characterized by either an increasing or a decreasing time delay as a function of energy above E 10 -3 TeV. As a remark, the cases in between the two regimes shows a time delay with respect to the energy with a plateau at E 10 -3 TeV which correspond to no time delay within this energy range. These specific cases can be of interest for constraining blazar modeling as discussed in Section 5.3.1. These origin of the two regimes are defined and described using the information on the time t max with respect to the light curve maxima.

Cooling driven regime: For the cases presenting a decreasing time delay phase at E > 10 -3 TeV, the study of their light curves shows that they all peak after t max . Thus, the evolution of the most energetic electrons is driven by the cooling effect, inducing the highest energy light curves to decay first. The energy dependent time delay is driven by the electron radiative cooling.

Acceleration driven regime:

For cases with increasing time delay at E > 10 -3 TeV, all light curves peak before t max indicating that electrons are still accelerating when light curves decay. In order to emit the highest energy photons, electrons require a long acceleration time to obtain enough energy and induces the highest energy light curves to peak later than the lower ones. The energy dependent time delay is driven by electrons acceleration.

As a result, the other parameters are scanned starting from each regimes. This provides more information on the influence of the others parameters and helps to identify if others transition between the two regimes are possible. The acceleration driven regime reference case is defined with B 0 = 65 mG while for the cooling driven regime it is the case with B 0 = 90 mG.

Magnetic field temporal index

As defined in Equation 4.11, the magnetic field strength decreases over time with an index m b . A value m b = 1 can be interpreted as a magnetic field with turbulent structure while m b = 2 can be interpreted as magnetic flux conservation assumption, in case of adiabatic expansion with a linear increase of the blob radius with time. A large m b value would shut down quickly the magnetic field strength leading to a shorter flare. Indeed, magnetic fields are one the main drive of blazar emission in order for electrons to emit SSC radiation. Besides, a smaller magnetic field strength induces a smaller radiative cooling for electrons.

The index m b is varied between 1 and 2 and is shown for the two cases representing the acceleration and cooling driven regime in Figure 5.8. For the case with B 0 = 65 mG, (Figure 5.8a) an overall increase of the time delay value is observed. While for E 10 -3 TeV the time delay evolution does not significantly change, the lower one becomes less decreasing as m b increases and evolves towards a regime like the adiabatic driven one. This transition is explained by the magnetic field decreasing so rapidly that the low energy flares quickly end up while electrons are still accelerating. As a consequence, the highest energy light curves peak at later due the time needed for the electrons to reach high γ values.

Starting from the cooling driven regime with B 0 = 90 mG (Figure 5.8b), increasing m b involves a transition to acceleration driven regime. In those cases, the rapid decrease of B(t) reduces the radiative cooling power which becomes not sufficient to drive the energy dependent delay. With a weaken cooling effect, electrons are able to be accelerated longer thus shifting t max to later time. As a consequence, when m b is larger, the highest energy light curve peaks at later time while the faster magnetic field decrease induces the lowest energy light curves to decay earlier.

Doppler factor

A modification of the Doppler factor δ implies a change in the Doppler boosting effect for the observer. It does not change the temporal evolution of the source in its own frame. It only leads to variations of the observed variability due to time contraction as well as a variations of the observed flux and energy. The energy dependent time delay for different δ ranging from 20 to 50 and for the two regimes is shown in Figure 5.9. For both regimes, the main effect on the time delay evolution is the time contraction due to Doppler boosting. As δ become larger, the time contraction effect is bigger leading to smaller observed time delays. Actually, the time delays obtained for different δ values are proportional to each other by the ratio of their Doppler factor. As example for the case with B 0 = 90 mG, the comparison of the delay at 1 TeV for δ values 40 and 20 gives ∆t δ=40 (1 TeV) 293 s and ∆t δ=20 (1 TeV) 600 s which correspond approximately to the ratio of the Doppler factors.

As a remark, the variation of the maximum energy considered for the time delay computation is an effect of the Doppler boosting on the flux value. This is a consequence of the choice to neglect light curves with flux value below 2 × 10 -15 cm -2 s -1 as mentioned in Section 4.4.1. Also, the energy shift due to Doppler boosting can lead to a change of the time delay sign which is not expected with a time contraction alone as seen for the case B 0 = 65 mG at E ≈ 1 TeV.

Acceleration parameters

The acceleration term, defined in Equation 4.15, depends on two parameters: the initial amplitude A 0 and the evolution index m a . A variation of the acceleration parameters leads to a change of the time t max with respect to the flare maxima which could imply a transition between the two time delay regimes.

The time delay obtained for different A 0 values, ranging from 4.0 × 10 -5 to 6.0 × 10 -5 s -1 is shown in Figure 5.10. For cases with B 0 = 65 mG (Figure 5.10a), a decreasing A 0 increases the time delay values for the highest energies. Indeed, as acceleration power is weaker, high energy electrons need more time to reach γ max which also has a smaller value because of the unchanged cooling effect. Besides, with a smaller γ max value the associated cooling effect is also weaker which shifts t max to later time. When A 0 increases, a transition between the two regimes occurs to a cooling driven regime. A case with a stronger acceleration implies a higher γ max value which is reached at an earlier time t max compared to weaker acceleration cases. As a result, a larger γ max induces a stronger cooling effect for the most energetic electrons and an earlier t max value implies a larger B(t max ) and an even stronger cooling effect. Hence, increasing the A 0 values brings electrons much faster to higher γ max values which then suffer more from radiative energy losses inducing a transition to cooling driven regime.

Starting from a cooling driven regime with B 0 = 90 mG (Figure 5.10b), the opposite situation occurs with a transition to an acceleration driven regime happening when A 0 decreases. Based on the previous explanation, this transition is easily interpreted as a smaller A 0 value induces a less efficient acceleration and lower γ max value. Hence, the cooling power at γ max become weaker and electrons have a longer acceleration time to reach their maximum energy. Therefore, t max is shifted to later time, and with a sufficient lower A 0 value after the light curve maxima thus leading to acceleration driven regime. In opposition, a larger A 0 value implies a more efficient acceleration with larger γ max reached at earlier time. In consequence, the cooling power become stronger and dominates earlier the most energetic electrons than cases with lower A 0 . This induces the highest energy light curves to peak at much earlier than the reference and therefore explain the more negative time delay. The same reasoning can be applied to the scans of the acceleration evolution index m a for the two initial magnetic field values (Figure 5.11). The parameter m a is varied between 4.5 and 5.9. Increasing m a induces a weaker acceleration which decrease faster. However, the scan of m a does not lead to a significant transition of regime. For the cases with B 0 = 90 mG (Figure 5.11b), a hint of transition is observed for large m a values when the acceleration becomes weaker and for B 0 = 65 mG, the overall delay decreases toward negative values when m a decreases corresponding to a stronger acceleration. Nevertheless, from the effect of acceleration deduced with variations of the A 0 parameter, a small enough m a value should improve the acceleration enough to induce a transition to the cooling driven regime. More generally, a strong acceleration leads to a cooling driven time delay due to more energetic electrons which suffer more from radiative energy losses. In opposition, a weak acceleration implies a longer acceleration time for electrons to reach a high energy enough to emit high energy photons thus leading to acceleration driven time delay.

Electron spectrum parameters

The electron spectrum used for the model follows a power law function with high energy cut-off (Equation 4.18). As a remark, the initial electron density K 0 is not considered for the scan because it does not affect the time evolution and is only a scaling parameter. Changing K 0 only change the number of electrons considered for the emission. However, when the electron spectrum index n changes, the proportion between low and high energies changes in the electron population, towards lower energies for high n values. In the transfer equation, changing n does not change the balance between acceleration and cooling effect for electrons. In practice, when n changes, γ max and t max remains at the same values..

The energy dependent time delays obtained with values from n = 2.2 to n = 2.8 are presented in Figure 5.12. For the cases with, B 0 = 65 mG (Figure 5.12a), the variations of the time delay are small when n changes and does not allow for the highlight of the influence of the parameter. Hence, its influence is evaluated only for cases with B 0 = 90 mG. In Figure 5.12b, starting from a cooling driven regime, a transition occurs to an acceleration driven regime when n is smaller. The transition is explained by the flare being shorter in overall when n is smaller while t max remains the same for all n values. This is the consequence of the electron population being in average more energetic when n is smaller leading to shorter duration. Hence, the light curves peak at earlier time when n become smaller. The transition occurs when n is small enough to produce light curve peaking earlier than t max . When looking back at cases with B 0 = 65 mG, the evolution of the time delay follows the same behavior. A larger n value induces an overall shorter flare, with light curves peaking at earlier time because they are produced by a more energetic electron population. The light curves maxima become shifted at earlier time than t max thus inducing larger delays for the highest energy light curve.

Energy [TeV]

Initial electron energy cut-off

The energy cut-off γ c,0 defines the maximum energy of electrons at the initial time t 0 of the time evolution (Equation 4.18). Using lower values of γ c,0 than 4 × 10 4 increases the time needed to reach high γ values, and to produce the highest energy photons. For higher γ c,0 values, electrons reach quickly high γ values and the γ max value is larger leading to a stronger cooling effect. The resulting time delays for γ c,0 ranging from 2 × 10 4 to 8 × 10 4 are shown in Figure 5.13 for the acceleration and cooling driven cases.

For the cases with B 0 = 65 mG (Figure 5.13a), decreasing γ c,0 lead to an increase of the time delay values for all energies. Indeed, starting from electrons with lower energies requires more time for them to be accelerated up to γ max . This leads for the highest energy light curve to peak later than the MeV-GeV. For larger γ c,0 values, more energetic electrons are considered at t 0 and so higher energies are achieved in a shorter time resulting in a shift of t max to earlier times. In this situation, the cooling effect becomes more important than acceleration and induces a transition from acceleration to cooling driven regime. For the cases with B 0 = 90 mG (Figure 5.13b), a similar behavior is observed. When decreasing γ c,0 , the acceleration of electrons lasts longer because the initial electrons population is less energetic. Electrons need more time to be accelerated up to γ max inducing a transition to acceleration driven regime. With larger γ c,0 , the radiative cooling becomes even more dominant over the electron acceleration especially for high γ values. Besides, γ max is reached at much earlier time when the magnetic field strength is larger leading to a larger cooling value at t max which becomes shifted to earlier time than cases with a smaller γ c,0 .

Discussion on adiabatic expansion

In this section, the adiabatic expansion is not considered and the energy dependent time delay appears to follow two regimes driven by either the acceleration or the radiative cooling of electrons. For both regimes, the evolution of the magnetic field is also important, explaining the origin of the low energy time-delay as well as the mechanism responsible for the highest energy light curve decay in the acceleration driven regime. For all the cases, time delays are investigated by comparing each phenomenon inducing either an increase (acceleration process) or a decrease (magnetic field and radiative cooling) of the flux for each energy value. The adiabatic expansion was removed from the model in order to simplify the interpretation on the origin of the time delay for each case.

When the adiabatic expansion is considered in the transfer equation, it dominates the evolution of electrons and explain the increasing time delay in function of energy shown in Figure 5.4 and explained in Section 4.4.3. On this specific scenario with adiabatic expansion, the determination of t max indicate that the electrons are still accelerating when the light curves decay (t max later than the light curve maxima). Hence, the main effect of adiabatic expansion for this scenario is the dilution of the electron population due to the increasing blob radius (Equation 4.17). The flux increase due to the acceleration of electrons of which γ < γ c t is not able to counterbalance the flux decrease induced by the decrease of the electron density. In the mean time, electrons are still increasing their maximum energy, towards γ max , explaining the increasing time delay as in the acceleration driven regime. However, there is no decreasing time delay at the lowest energies which is the consequence of adiabatic expansion which induces very early light curves decay when electrons does not have a high γ c (t).

As an example, a second case with adiabatic expansion is considered with B 0 = 90 mG, which presents a cooling driven regime. The resulting time delay with respect to energy is shown in Figure 5.14. As for the case B 0 = 65 mG, the time delay with B 0 = 90 mG is only increasing with energy indicating that the impact of the adiabatic expansion dominate the flare evolution. A larger B 0 induces an enhanced emission and so more energy loss for electrons. The acceleration becomes less efficient and induces shorter delay values. As a result, the two cases presented with adiabatic expansion produce an acceleration driven time delay. This regime is the consequence of a strong adiabatic effect which dominate the generation of a flare. Considering a smaller R 0 value should reduce the adiabatic expansion and allows for a transition towards a case with a cooling driven regime for the time delay.

As a final remark concerning the study on the model parameters, R 0 was not investigated because it is related to the definition of the initial time t 0 (Equation 4.13). Varying R 0 would change the t 0 value and thus the evolution of the magnetic field strength and the acceleration process making difficult any interpretation of the time delay. Besides, a different t 0 while the others parameters are unchanged could lead to cases where no flare can be generated. In the hypothesis where t 0 is fixed, variations of the initial radius does not affect the electrons evolution as the transfer equation is defined in density. R 0 as K 0 is simply a scaling parameter, by changing the total number of electrons at the emission stage of the model.

Discussion on intrinsic and LIV effects

Constraining the model with observations

With the model presented in this thesis, all cases reveal the presence of an energy dependent time delay at γ-ray energies. In total, two distinct regimes are found for the time delay, corre- sponding to the mechanism driving the electron evolution when the light curves peak. These regimes are the cooling driven and the acceleration driven regimes. In addition, some specific cases corresponding to the transition between the acceleration and the cooling driven regimes show no delay between GeV and TeV emission. Thus, even for observations limited to narrower spectral range around TeV energies, the information on the energy dependent time delay can be used in order to constrain the model parameters. Indeed, if a data set shows an energy dependent delay corresponding to the acceleration driven regime, all the cases producing cooling driven regimes cannot reproduce the emission of the source and have to be rejected. Clearly, the time delay is a new observable which can be used to constrain the modeling of blazars.

As an example, the observation of a flare from the blazar Markarian 501 in 2005 (Albert et al., 2007) revealed a non-zero time delay in the GeV-TeV range. This remains up to now the unique case reporting the detection of a time delay at γ-ray energies with a blazar. The authors found a positive delay (high energy photons arriving after the low energy ones) which indicates a time delay following the acceleration driven regime. This corresponds to an acceleration process which has to last long enough to counterbalance the radiative losses in order to produce the highest energy photons at the latest times.

Another example is given by the exceptional flare of PKS 2155-304 observed in 2006[START_REF] Aharonian | An Exceptional Very High Energy Gamma-Ray Flare of PKS 2155-304[END_REF] where a CCF method was applied to the light curves [START_REF] Aharonian | Limits on an Energy Dependence of the Speed of Light from a Flare of the Active Galaxy PKS 2155-304[END_REF] between 400-800 GeV and above 800 GeV. No significant delay was found in the data. Thus, the modeling of this flare should give parameter producing a negligible energy dependent delay. This corresponds to cases with parameters producing a time delay at the transition between acceleration and cooling driven regimes. For such case, the constraints of generating no delay to reproduce the data greatly constrain the available range for the parameters of the source and its giant flare.

Comparison with LIV delays: Focus on high energy

For the search of LIV signatures with IACT, studies are done in the GeV-TeV range. To be able to compare intrinsic delays from the model with the expected lags coming from the LIV effect, the intrinsic delays are re-calculated using a reference light curve at higher energy, in the range 42 -74 GeV, corresponding to the zero origin of the delay. Then, comparing the energy dependency of the intrinsic time delay in this energy range where LIV is hoped to be significant can provide direct constrains on specific QG models. Indeed as mentioned in Section 2.3.2, some models produce only one specific time delay such as [START_REF] Amelino-Camelia | Distance Measurement and Wave Dispersion in a Liouville-String Approach to Quantum Gravity[END_REF], [START_REF] Ellis | Quantum-Gravitational Diffusion and Stochastic Fluctuations in the Velocity of Light[END_REF] with a positive linear energy dependent delays. To look for the energy dependency of the intrinsic delay, intrinsic delays are fitted with a function inspired from Equation 2.6 which gives:

∆t = ξ × (E α i -E α 0 ) (5.2)
where α is the energy dependency index, ξ the amplitude of the delay in s TeV -α and E 0 the mean energy of the reference light curve (58.5 GeV). The time delays in the GeV-TeV range obtained with the reference case for B 0 = 65 mG and B 0 = 90 mG, with respectively increasing and decreasing time delays at high energies, are presented in Figure 5.15 and fitted with the function of Equation 5.2.

Energy [TeV] For both cases, the energy dependent delays evolve with an index α which is not linear, thus different from the QG model predictions from [START_REF] Amelino-Camelia | Distance Measurement and Wave Dispersion in a Liouville-String Approach to Quantum Gravity[END_REF] and [START_REF] Ellis | Quantum-Gravitational Diffusion and Stochastic Fluctuations in the Velocity of Light[END_REF]. As a result, a blazar flare following one of these cases could be used to efficiently test these two QG models. A more general study is done to identify the range α can take on all cases investigated in Section 5.2. For the special cases in between the two regimes, the energy dependent delay has a plateau in the GeV-TeV range corresponding to no delay at these energies. Hence, α cannot be evaluated since there is no time delay. Besides, some cases are not able to produce flare above 250 GeV, which does not provide a large enough energy coverage to estimate α. The resulting of this determination of the energy evolution index of the delay is shown on Table 5.1 and 5.2. The resulting energy evolution index of the delay is found within the range α ∈ [0.45 -0.85] for cases where a significant increasing or decreasing delay is observed (flat delay and low energetic cases excluded). As a general results from this study, the use of blazar observations can be done to constrain specific QG models predicting a LIV delay which is not in the range of α. Indeed, if a specific model predicts a LIV energy dependent delay with an index in the α range, this delay will probably be mixed with intrinsic delays On the contrary, if this index is outside the α range, the two delays can be discriminated and this model can be tested. -181 ± 23 0.68 ± 0.1 Table 5.1: Energy dependent time delay amplitude ξ and power index α for various initial magnetic field strengths in the GeV-TeV energy range. The missing values (shown with a dash) could not be evaluated because there was no significant delay in the considered energy range.

B 0 (mG) ξ [s.TeV -α ]
In order to compare the intrinsic time delay obtained with the model with a linear LIV delay defined in Equation 2.6, a E QG value is assumed at the Planck scale. Multiple values of redshift are investigated for the LIV delay while the intrinsic delay obtained with redshift z = 0.03 is kept identical for illustration purpose. Indeed, if z is varied in the model, the resulting energy dependent delay does not significantly change whereas the observed flux and the EBL absorption reduce the maximum energy used for the delay computation and prevent a decent comparison of the two effects. As a remark, using the scaling parameters of the model such as R 0 (without changing t 0 ) and K 0 could allow for the production of a sufficient emission to mimic the results at z = 0.03 with larger z values. The expected linear LIV delay with a QG energy E QG at the Planck scale is presented in Figure 5.16 for the two cases representing the two regimes (acceleration and cooling driven delays). For B 0 = 65 mG, the expected LIV effect is supposed subluminal to produce an positive time delay and is found to be always negligible compared to the intrinsic one. A superluminal effect LIV is assumed for the case B 0 = 90 mG and shows that with a high redshift value, the LIV effect become sufficient to produce a delay as large as intrinsic delays. From this comparison, if a delay is found at more than hundreds of seconds, a linear LIV effect at the Planck scale is expected to produce a negligible delay compare to the intrinsic one. In this situation, measuring the LIV delay to constrain QG models seems difficult unless the two effect can be efficiently disentangled. for all investigated parameters for the two initial magnetic field strength values taken as benchmark values. Bold lines correspond to these benchmark cases, given by the parameters in Table 4.1. The missing values (shown with a dash) could not be evaluated because there was no significant delay in the considered energy range or because the maximum energy emitted was below 250 GeV.

B 0 = 65 mG B 0 = 90 mG Parameter ξ [s.TeV -α ] α ξ [s.TeV -α ] α m b 1.

Temporal characteristic of intrinsic delay

From the Fourier transform analysis of the delay between two light curves at different energies, [START_REF] Lewis | Time-dependent Electron Acceleration in Blazar Transients: X-Ray Time Lags and Spectral Formation[END_REF] obtain the time delay as function of the Fourier frequency which is inversely proportional to the time. From there results, shown in Figure 4.2, a break is observed in the Fourier transform of the time delay. This break occurs at a frequency corresponding to the time when low energy photons starts to arrive before the high energy ones. This temporal evolution of the delay along the flare is a consequence of the mechanisms generating the flare, while the LIV delay is expected to induce a constant delay along the flare as it affect all photons during their propagation from the source to the Earth. To investigate this temporal evolution of the delay with the model presented in this thesis, the delay is studied along the whole duration of the flare.

To determine this temporal evolution, the delay is estimated in a simpler way with no Fourier transform. The delay temporal evolution is evaluated by comparing the time difference for two light curves to reach the same normalized flux value. The temporal evolution is deduced by computing the delay along all the flux values taken by one of the two light curve. To compare this work with the one by [START_REF] Lewis | Time-dependent Electron Acceleration in Blazar Transients: X-Ray Time Lags and Spectral Formation[END_REF], their definition of time delay is adopted where the high energy light curve is taken as the reference. This definition is in opposition to the definition considered in this thesis until now. With this definition, a positive delay corresponds to low energies arriving after the high energies.

For this study, the two light curves chosen for the comparison are in the energy ranges 200 -400 GeV and 2.6 -4.7 TeV. Two cases are considered following the reference set of parameter (Table 4.1) for the two time delay regimes represented by B 0 = 65 mG and B 0 = 90 mG. Their light curves as well as their time delay temporal evolution are shown respectively in Figure 5.17 and Figure 5.18. The temporal evolution of the delays obtained for both cases are similar to the one obtained by [START_REF] Lewis | Time-dependent Electron Acceleration in Blazar Transients: X-Ray Time Lags and Spectral Formation[END_REF]: at large time (small frequency), the delay is large because the high energy light curve decays faster than the low one. A break occurs at the specific frequency corresponding to the time when a given flux is reached earlier for the high energy light curve flux than for the low energy one. At small time (large frequency), the delay is negative because the high energy flare rises after the low one, due to the time needed for electrons to be accelerated. Negative delays cannot be represented in the figure because of logarithm scale. When the opposite definition of the time delay is considered, an opposite temporal evolution of the time delay is obtained. These results obtained from the model are in agreement with [START_REF] Lewis | Time-dependent Electron Acceleration in Blazar Transients: X-Ray Time Lags and Spectral Formation[END_REF] and seems to confirm this characteristic of the delay which is variable all along the flare. As a remark, this temporal evolution explains the result obtained in Section 4.4.2, where the CCF is not able to measure the injected delay. This is the consequence of the method using the full light curves to determine one delay value whereas light curves present a variable delay with time. To evaluate what a LIV effect would give with such a representation of the delay, two light curves are simulated, following the same asymmetric Gaussian shapes than the model light curves and no intrinsic effect. An average LIV delay between the two simulated light curves is injected of ∆t LIV = 500 s. The resulting temporal evolution of the delay is shown in Figure 5.19 and shows a constant delay along the whole flare. As a result, the temporal evolution of the delay seems to be a property that can trace the presence of an intrinsic effect if there is a variation of the delay with respect to the Fourier frequency or time.

Conclusion

The model developed for this work revealed the presence of energy dependent time delays. These delays at the highest energies were found to follow two regimes, driven by two processes which govern the electron population evolution: acceleration or radiative cooling. In presence of adiabatic expansion, the time delay evolution was found to follow the acceleration driven regime, in the case considered for this work. The absence of adiabatic expansion allowed for the highlight of the competition between the acceleration and radiative cooling. The magnetic field decrease was found to be crucial for the acceleration driven regime, as its effect induce the decay of low energy light curves while electrons are still accelerating, inducing high energy light curves peaking at later time than low energy ones. The cooling driven regime appears when a strong acceleration bring quickly electrons to very high energies which then suffer from more energy losses by radiative cooling than energy gain by acceleration. In those cases, the high energy light curves peak earlier as compared to lower energy ones which decay at later because of the magnetic field decrease and a longer cooling time.

The presence of energy dependent time delays in blazar flares leads to several implications, both for blazar modeling or fundamental physics. The time delay evolution with energy can inform about in situ physical parameters leading to the generation of a flare. Indeed, the parameters used to model a flare must be able to reproduce the same time delay regime, thus providing constraints for blazar modeling. For fundamental physics, the time delay is sometimes used to search for Lorentz Invariance Violation signature or for the measurement of the equivalence principle [START_REF] Wei | Tests of the Einstein Equivalence Principle Using TeV Blazars[END_REF]. For these studies, intrinsic time delays were usually neglected because until now, no significant delay has been found at TeV energies. However, if a significant delay is once measured one day, it will be necessary to be able to disentangle between an intrinsic effect and an effect related to new physics. The intrinsic delay energy power index α from the model was found between 0.45 and 0.85. As a result, specific Quantum Gravity models can be tested with blazar flares if their predictions on LIV delays depend on the energy with an index different from this range. Another property was also investigated which is the temporal evolution of the time delay. From the definition of the delay using a Fourier transform, proposed by [START_REF] Lewis | Time-dependent Electron Acceleration in Blazar Transients: X-Ray Time Lags and Spectral Formation[END_REF], the model presented here shows the same behavior with respect to time as the one they obtain at X-ray energies. This behavior is not expected from a LIV effect which brings tool to identify the presence of intrinsic effect.

A publication about the time dependent blazar flare model, presented in Chapter 4, and its results presented in this Chapter is in preparation.

However, such studies rely on a precise determination of light curves in multiple energy bands. This precision is not currently achievable, especially at TeV energies where current IACT suffer from a lack of statistics due to the power law energy spectrum of blazar flare as well as EBL absorption. The precise measurement of flux with time is limited by the number of photons detected, which is generally of the order of several hundreds for a blazar flare. In order to perform accurate energy dependent time delay measurements for the search of LIV signatures on TeV data, the use of individual photon information provides better results than Cross-Correlation Functions or the Peak Position Method, which use light curves. Such a method, the Maximum Likelihood method, is presented in the next chapter.

In LIV searches through the time of flight approach, energy-dependent time delays are investigated with astrophysical sources (Chapter 2) at γ-ray energies. Different methods were used to measure LIV time delays, like CCF or Maximum Likelihood (ML), in order to constrain the QG energy scale. [START_REF] Vasileiou | Constraints on Lorentz invariance violation from Fermi-Large Area Telescope observations of gamma-ray bursts[END_REF] applied these method to GRB observations and obtain the best constraints on E QG with ML. This method uses the individual photon information to measure the LIV delay and allows for the analysis of small sample data set. Thus this method, which was initially proposed by [START_REF] Martínez | A new approach to study energy-dependent arrival delays on photons from astrophysical sources[END_REF] is used in this work. From a mock function, which describes a given data set, the ML compares how likely the data set fits this mock function depending on one or several parameter. For the LIV study, in a simplified way, the ML compares the arrival time of high energy photons to a non LIV delayed light curve pattern which define the mock function. The ML gives then the most probable LIV time-delay in order for the photon arrival times to match the mock function.

In the present work, a complete toy MC simulation tool and an analysis tool based on the ML have been developed, allowing for the search of LIV signature with the time of flight approach. This chapter presents the analysis procedure and investigates the different hypothesis made for the analysis with a toy Monte Carlo (MC) simulation tool. After a description of the ML analysis, the MC procedure is going to be described and used on an illustrative simple case. Then to test the behavior and the precision of the ML, several scenarios are investigated with the MC simulation and the analysis tool. Finally a new development on the ML for LIV search will be presented in the view of real data analysis.

General view on maximum likelihood method for LIV

Parameter of interest τ n

To measure the LIV time delay, it is convenient to define the parameter τ n which is the time delay per unit of energy and the parameter of interest in the ML. It is expressed from the Equation 2.6 as:

τ n = ∆t E n = 1 + n 2H 0 1 E QG n κ n (z), (6.1)
with

κ n (z) = z 0 (1 + z ) n dz Ω m (1 + z ) 3 + Ω Λ . (6.2)
where n is the order of the LIV effect: linear (n=1) or quadratic (n=2), and κ n (z) takes into account the expansion of the Universe during the propagation of photons, depending on the redshift of the source considered. In this expression, τ n is the time delay between a photon of energy E compare to a photon of negligible energy, i.e. non delayed photon. Thus, it allows to easily compute such LIV time delay for a photon with energy E:

∆t LIV = τ n × E n . (6.
3)

The parameter τ n is measured with ML from data and allows for the deduction of the QG energy scale E QG as follows:

E QG,n = n 1 + n 2H 0 1 τ n κ n (z). (6.4)

Standard Probability Density Function

The mock function is a Probability Density Function (PDF). It corresponds to the probability of detecting a photon arriving at time t and measured with the energy E. To easily defined the PDF, the presence of intrinsic time delay is neglected. The complete PDF takes into account the source properties, the IRF of the instrument and the LIV effect:

P (E, t, τ n ) = 1 N (τ n ) ∞ E =0 A(E , t) R(E, E ) Γ(E ) F (t -E n × τ n , E )dE , ( 6.5)
where:

• F (t) corresponds to the photon time distribution non-delayed by LIV effect. This is the key ingredient, defining the light curve shape, and determining the position of the different pulses which are crucial for the time delay measurement. In Equation 6.5, the function argument is the measured time, corrected by the expected LIV effect given by E n × τ n .

• Γ(E ) is the energy distribution of the observed photons. This term allows for the introduction of a weight of high energy photons with respect to lower ones in the PDF. Indeed the high energy photons are rare but also more efficient to constrain the LIV parameter as they are more delayed so more constraining. The energy distribution at high energy can often be approximated by a power law function.

• R(E, E ) is the instrument energy resolution. It corresponds to the probability distribution of measuring a photon of true energy E at the energy E. It follows a Gaussian function in natural logarithm of E and can be written as:

R(E, E ) = A exp   log(E) -log(E ) -b(E ) 2 2r(E ) 2   , ( 6.6)
where A is an amplitude, b(E ) is the energy bias and r(E ) is the energy resolution, also noted as ∆E/E. At high energy, bias is usually negligible and the energy resolution is usually of the order of ∆E/E ≈ 10% for IACT. The energy resolution can impact the LIV time delay estimation since it is energy dependent.

• A(E , t) is the collection area of the instrument and is expressed in m 2 . This factor takes into account the change in the observing conditions during data taking. These changes can lead to significant change in the detection probability of photons.

• N (τ n ) is the normalization factor. It is expressed as :

N (τ n ) = tmax 0 E 1 E 0 P (E, t, τ n )dE dt, ( 6.7)
where [E 0 , E 1 ] is the energy range of the events considered for the ML method. The normalization factor is especially important when evaluating the PDF at large values of τ n . In case of events measured close to one edge of the observations windows, they can be evaluated at a time when the time distribution cannot be determined, i.e. outside the windows. If not taken into account this lead to a significant bias in τ n reconstructed values.

Likelihood function

The PDF is used in order to evaluate the Likelihood function L(τ n ), expressed as:

L(τ n ) = N i=1 P (E i , t i , τ n ), (6.8)
with N being the number of events used in the ML. For practical reason it is easier to compute the log-likelihood function, leading to the expression:

log L(τ n ) = N i=1 log P (E i , t i , τ n ) . (6.9)
The log-likelihood function computes the probability of a given data set to match the mock function, via the PDF, considering the value of τ n . It presents a maximum at the most probable value τ best n , considered as the best estimate. This function has the property to have an easy way to estimate the uncertainty on the measured τ best n value. Following Cowan (1997), one can rewrite the log-likelihood function like:

log L(τ n ) = log L(τ best n ) - (τ n -τ best n ) 2 2σ τ , ( 6.10)
where σ τ is the standard deviation of τ n . From this equation, the function can be estimated at τ n = τ best n ± kσ τ , and gives: (6.13) From Equation 6.13 the uncertainties on τ best n at kσ τ confidence level can be deduced at the τ n values for which -2∆ log L(τ n ) = k 2 . This allows for a simple determination of statistical uncertainties, with a graphical method, directly from the likelihood function itself. An example is shown on Figure 6.2.

log L(τ best n ± kσ τ ) = log L(τ best n ) - (τ best n ± kσ -τ best n ) 2 2σ τ , (6.11) log L(τ best n ) -log L(τ n ) = k 2 2 , (6.12) -2∆ log L(τ n ) = k 2 .

Simplifications and approximations

Determination of F(t)

Equation 6.5 gives the complete formula to define the mock up function needed for the ML. But some approximations can be done to be able to estimate all the factors. Also, depending on the case of interest, some terms in the PDF can be neglected or approximated. The function F (t, E ) is a priori not accessible since it would require to observe at the source level (i.e. without propagation effects). Besides, the presence of intrinsic effects are neglected, meaning that F (t, E ) = F (t). To determine this function, the lowest energy events are used to estimate F (t) but one assumption is made which consists of neglecting the LIV effect in the energy range used to determine F (t). This has the benefit of having an easy way to define the time function. The main drawback is that the photons used to estimate F (t) cannot be used again to estimate precision and special cases inducing a bad reconstruction of the lag in order to improve the method.

Simulation inputs

The MC generator requires a list of inputs in order to generate and analyze the data. These inputs can be chosen to be either as close as possible to real data or to investigate specific cases. These inputs are:

• The time distribution F sim (t) used for the random generation of photon detection time.

This function is built using Gaussian functions, either one or the sum of several ones.

• The energy distribution Γ sim (E). It generally follows a power law function but can also be either a curved power law, a broken power law or a power law with a cut off.

• The acceptance of the instrument. As discussed in Section 6.2.2, this input can be neglected if the acceptance is time independent and assumed constant over the whole energy range considered for the simulations.

• The energy resolution of the instrument. The smearing of photon energy can either strictly follows the true resolution of the instrument, or can be approximated with a fixed value which is typically ∆E E ∼ 10% for IACT. • The injected LIV parameter τ inj n , associated with n taken as 1 or 2 for respectively a linear or quadratic LIV effect.

• The energy ranges used for the ML which are: the "template region" at low energy for the determination of F (t) and the "likelihood region" at high energy for the ML calculations. These energy ranges are defined with an energy cut value E cut which is usually defined with the median energy of the data set to have the same number of photon in each range.

• The total number of events. This number can be given for the whole simulated data set which induces a fluctuation on the number of events in the likelihood and the template regions due to the random generation. Fixed numbers of event for each region can also be given, ensuring the same number of events for all the simulated data sets using the same input parameters.

• The definition of functions F (t) and Γ(E) used for the mock-up function. In MC simulations, the functions used for the generation of photons, F sim (t) and Γ sim (E), can be used for the ML, ensuring that the PDF describes the true origin of photons. For real data, this cannot be applied so the hypothesis of using low energy photons to determine F (t), and high energy photons for Γ(E) have to be considered. The differences of the two possibilities will be studied in Section 6.4.

Toy Monte-Carlo procedure

From these inputs, the MC tool randomly generates photons defined by their true time and energy which follow respectively the F sim (t) and Γ sim (E) distributions. The LIV time delay is added to the true time using the true energy and then IRF are applied to accept (or reject) photons and smear the true energy into the measured one taking into account the energy resolution. Finally, photons are split into the template and likelihood regions and the log-likelihood function is computed to determine the reconstructed value τ best n of the data set. A large number data sets can be generated and analyzed, allowing for the obtention of the distribution of the τ best n . The distribution usually has a Gaussian shape. The position and the standard deviation of the Gaussian give respectively the bias and the 68% statistical uncertainty of reconstructed LIV parameter. In the case where the τ best n does not follow a Gaussian shape, the median of the distribution can be used as the bias. The 68% statistical uncertainty can be obtained with the median of the τ best n ± 1σ distribution obtained with the graphical method discussed in Section 6.1.3.

Application on a simple case

To illustrate the MC tool and the ML analysis, 1000 data sets are generated considering a simple scenario. The time distribution is chosen as a simple Gaussian pulse with a mean µ = 2000 s and a width σ = 750 s. The energy distribution is taken as a power law function E -α with index α = 2.2. The acceptance is considered constant for the whole energy range. A perfect energy resolution is considered, meaning that E = E . Events are generated between 0.4 and 10 TeV, split at 0.8 TeV to define the template and likelihood regions. A fixed number of 500 events are generated for each region, which is similar to the number of photons obtained in the data set presented in Chapter 7. At first no LIV time lag is injected in order to determine the precision of the method and only linear LIV (n=1) is investigated with the ML. The input parameters are summarized in Table 6.1. Control plots for the distributions F sim and Γ sim are shown in Figure 6.1.

For the analysis, the simulated functions are used in the PDF. The log-likelihood function from a single realization is shown in Figure 6.2 as an example. The overall distribution of reconstructed τ best 1 from the 1000 simulations data sets is presented in Figure 6.3 and follows a Gaussian shape. The results of a Gaussian function fit on the τ best 1 distribution yield no significant bias, and a statistical uncertainties σ stat 11 s.TeV -1 . while the left and right plots correspond to lower (left) and upper (right) limits of the 68% confidence level interval, determined via the graphical method from the log-likelihood function. They are fitted with a Gaussian function to determine the bias and accuracy of the method via the mean and the width of the function.

Calibration plot

To evaluate the performance on the LIV parameter reconstruction for linear and quadratic effects, the procedure is redone using different values of injected time lag τ inj n . For each injected values: -20, -10, -5, 0, 5, 10, 20 s.TeV -n , a 1000 data sets are generated and analyzed to obtain the τ best n value associated with it statistical uncertainty σ stat from the distribution of reconstructed LIV parameter. By comparing the best reconstructed values with the injected ones, the efficiency of the method can be deduced. The resulting calibration plots where τ best n is plotted versus τ inj n are presented on Figure 6.4 for linear (n = 1) and quadratic (n = 2) LIV effects. Both calibration plots show perfect reconstructions with a slope 1 and a y-intercept, or bias, negligible compared to the statistical uncertainties (grey shaded area).

One can notice that the precision on the LIV parameter τ n is better for the quadratic effect than on the linear one. This is induced by the high energy photons of the simulations, up to 10 TeV. To give an example, let's assumed a value τ 2 = 1 s.TeV -2 for a quadratic effect. Equation 6.3 gives a LIV time delay for a 10 TeV photon of ∆t 2 ∼ 10 2 s whereas for a linear effect with τ 1 = 1 s.TeV -1 , the delay is ∆t 1 ∼ 10 2 s. Considering this facts, the results on quadratic LIV parameter are more constraining than for linear one due to the large LIV delay involved.

As a remark, for quadratic LIV effect the photon maximum energy is the most relevant properties in order to obtain the most precise evaluation of τ 2 . Indeed, from the previous example with τ 2 = 1 s.TeV -2 , a 0.5 TeV would experienced a the LIV delay of ∆t LIV ∼ 0.5 2 s. For this case, such a delay would be difficult to measured and thus worsening the precision on the LIV parameter. A quadratic LIV effect is better constrained if the energy is typically well above 1 TeV.

Monte Carlo simulation studies

In order to check the validity of some approximations given in Section 6.2, the toy MC tool is used to generate different scenario. From the simple scenario described in Section 6.3.3 and its parameter on Table 6.1, different inputs are modified and their effect is compared using calibration plots. All the scenarios investigated in this section are summarized in Table 6

Energy resolution scenario

For this scenario, the energy resolution assumption discussed in Section 6.2.3 of removing the energy resolution integral in the PDF is studied. In the simple scenario, a perfect energy resolution is considered while this scenario consider an energy resolution of ∆E/E = 10%. If the assumption on the energy resolution in the PDF is correct, the calibration plot should not change significantly. The other inputs are the same as the ones defined in Table 6.1, and the MC procedure follows the one defined in section 6.3. The resulting calibration plots are given on Figure 6.5 for linear and quadratic LIV effect and have to be compared with Figure 6.4. Adding energy resolution in the simulation does not significantly change the calibration plots, validating the hypothesis of neglecting the energy resolution term in the PDF.

Fitted template scenario

In case of real data, the true F (t) is not accessible in opposition to MC simulation where F sim (t) is known. One way to determine it is to fit the time distribution of the template region, usually with Gaussian functions. This scenario investigates the impact of parameterizing the template distribution on the τ best n value. The simulated events in the template region are used to determine F (t). The mathematical function used to fit the template region time distribution is the same as F sim (t): a simple Gaussian function. In case of real data, this mathematical function has also to be determined, comparing different fitted function via a χ 2 /ndf comparison. For this scenario, the energy resolution is applied and the LIV effect is not applied for the generation of template photons (as it will be discussed in the next subsection). Figure 6.6 presents calibration plots for linear and quadratic LIV effects. Fitting the template photons to deduce F (t) marginally decreases the calibration performance and also increases the uncertainties on τ best n as compared to Figure 6.5. This loss of sensitivity is explained by errors induced on F (t) parameters as compared to the ones of F sim (t) during the parameterization, mostly the position of the peak with is crucial to determine the LIV time delay. As a remark, the quadratic calibration plot is less impacted than the linear one. This is explain by the uncertainties on the position of the peak of the order of tens of seconds (Figure 6.1). A quadratic LIV time delay, expected to be larger than this uncertainties whereas a linear LIV delay is of the same order for this scenario.

Lagged template scenario

For the fitted template scenario, LIV effect was ignored in the template region. But in case of real data, all energies are affected by LIV. This can lead to an error in the reconstructed LIV parameter as mentioned in Section 6.4.2. In this lagged template scenario, the energy resolution is applied, F (t) is determined by fitting the photon time distribution in the template energy range and the LIV effect is applied on template photons. For this scenario, only linear LIV effect is shown but similar results are obtained with quadratic LIV effect. The Figure 6.7a shows the resulting calibration plot. The reconstruction power becomes worse than previous scenarios with a slope of 0.82, inducing a τ best 1 value systematically smaller than τ inj 1 . The origin of this lagged template effect, comes from the hypothesis of neglecting any LIV delay in the template region. But for this scenario, the energy range of the template region (0.3 -0.8 TeV) is not negligible which induces a significant delay in this range. As F (t) is determined from the template photons, the time function is lagged while the LIV delay was neglected in the template energy range. Considering this, expressing τ best 1 with respect to τ inj 1 gives: (6.15) (6.16) .17) where ∆t temp LIV is the averaged LIV delay in the template region, and E n T the averaged template energy to the power n. This equation demonstrates that the reconstructed LIV parameter is systematically smaller than the injected one.

τ best n × E n + ∆t temp LIV = τ inj n × E n ,
τ best n × E n + τ inj n × E n T = τ inj n × E n ,
τ best n = τ inj n × E n -E n T E n , ( 6
As a remark, this new template correction should not significantly affect previous studies on LIV signatures such as with the flares of PKS 2155-304 (H.E.S.S. Collaboration, 2011) ) or PG 1553+113 [START_REF] Abramowski | The 2012 Flare of PG 1553+113 Seen with H.E.S.S. and Fermi-LAT[END_REF] as they both have low energy photon (respectively 250 GeV and 300 GeV) and show no LIV effect. However it will be used for the data analysis presented in Chapter 7 with the flare of Mrk 501 which shows low energy photon starting from 1.3 TeV where LIV effect cannot be neglected for the template region.

Template correction

In order to deal with a lagged template, a correction was developed to solve the misreconstruction of τ n . It consists in an evaluation of the LIV time delay in the template region using E n T . This template LIV delay estimation is then subtracted to the argument of the function F (t) in the PDF. The PDF on Equation 6.5 becomes then: As a result, since the template correction depends on the LIV parameter, this correction induces an additional source of uncertainties in the log-likelihood function evaluation.

P cor (E, t, τ n ) = 1 N (τ n ) ∞ E =0 A(E , t) R(E, E ) Γ(E ) F t -E n × τ n -E n T × τ n dE . ( 6
To determine the efficiency of this correction, another scenario is investigated based on the lagged template scenario using the PDF of Equation 6.18. However, to generate proper simulations, one has to inject the correct LIV time-lag ∆t LIV in case of the template correction approach to reconstruct the same τ n as τ inj n . From Equation 6.16 the right injected LIV time delay in case of template correction is:

∆t LIV = τ inj n × E n + τ inj n × E n T (6.19)
To evaluate E n T at the simulation stage, the injected energy distribution Γ inj is used to evaluate the mean energy in the template energy range [E 0 ,E 1 ] by computing its first moment: (6.20) which corresponds to the average value of all the different E n T obtained in the simulated data set. For the analysis, E n T is estimated for each data set by averaging the measured energy of the photons to the power n in the template range. Figure 6.7b shows the calibration plot for linear LIV effect, and a perfect calibration is recovered with a slope close to 1 and the bias remains unchanged. As expected, the error bars become larger compared to the case without correction due to the addition of the template correction.

E n T = E 1 E 0 E n × Γ(E)dE,
Mrk 501 is a TeV blazar, located at redshift z = 0.034 [START_REF] Moles | Two elliptical galaxies with active nuclei -NGC 6212 and MKN 501[END_REF] and discovered in gamma rays by [START_REF] Quinn | Detection of Gamma Rays with E > 300 GeV from Markarian 501[END_REF]. It is one of the brightest sources in gamma rays, with flaring activities. In 1997, a flare from the source was observed by different γ-ray instruments like Whipple [START_REF] Catanese | Multiwavelength Observations of a Flare from Markarian 501[END_REF] and HEGRA (Barrau et al., 1997), showing an integral flux up to 4 times the one of the Crab Nebula. These observations provide a measurement of the source spectrum extending up to 20 TeV [START_REF] Djannati-Atai | Very High Energy Gamma-ray spectral properties of MKN 501 from CAT Čerenkov telescope observations in 1997[END_REF][START_REF] Aharonian | Reanalysis of the high energy cutoff of the 1997 Mkn 501 TeV energy spectrum[END_REF]. A second flare was measured by MAGIC in 2005 (Albert et al., 2007), with a very high flux of 3 Crabs which allowed for the search of LIV signatures [START_REF] Albert | Probing quantum gravity using photons from a flare of the active galactic nucleus Markarian 501 observed by the MAGIC telescope[END_REF][START_REF] Martínez | A new approach to study energy-dependent arrival delays on photons from astrophysical sources[END_REF].

During the night of June 23-24 2014, Mrk 501 went into another flaring state, observed by H.E.S.S. telescopes (Stegmann, 2014). After a presentation of the data set and its analysis, the ML is going to be applied on the data set to search for a LIV-induced energy-dependent delay. A dedicated study with the toy MC tool follows in order to estimate the systematic and statistical uncertainties related to the data set. Finally, the QG energy scale E QG,n for linear and quadratic LIV effect are assessed from the ML analysis and compared to other result from others flares observed by H.E.S.S.

Markarian 501 data analysis

H.E.S.S. observation conditions

The observation of this flare was triggered by the telescope FACT (Stegmann, 2014). Four runs were taken (∼2 hours of data) with the full H.E.S.S. array. Due to the position of the source in the sky, the observation were done with large ZA [63 • -65 • ]. With such an incoming angle for γ rays, the electromagnetic showers produced cross a larger depth of atmosphere than γ rays arriving with smaller ZA. This leads to an enhanced absorption of Cherenkov light through atmosphere. Besides, due to the geometry of the shower, the Cherenkov light pool image is larger at ground level resulting in a lower Cherenkov light density. This effect is illustrated on Figure 7.1. Consequently, low energy γ rays produce Cherenkov light pools at high zenith which are too faint to be detected, limiting the lowest energy observable. On the contrary high energy events which are intense enough to be detected benefit from the larger Cherenkov light pool on the ground, resulting in a larger collection area. Furthermore, the enhanced atmosphere absorption also affects the background noise. Hadronic showers, which are generally less intense than electromagnetic ones are even fainter ensuing a better γ-hadron rejection.

H.E.S.S. analysis results

To analyze the data, the so called Combined3 analysis, described in Section 1.4.1 is used. Because of the high zenith angle observations limiting the minimum energy observable, the use of Combined3 analysis is well suited and benefit both from the low energy performances of mono reconstruction with CT5 and the high energy performances of the stereo reconstruction. This maximization of the energy coverage is crucial for the search of a LIV signature. The detailed cuts of the analysis are shown in Table 1.1.

The control plots of the analysis are the θ 2 distribution (Figure 7.2) and the event map distribution (Figure 7.3). The map shows a point like source with almost no background, as also showed in the θ 2 histogram. The significance of the detection exceeds 80σ.

The total number of events from the source is 1930, with a high signal over background ratio of 46.5 in average. The individual runs event statistics are presented in Table 7.1. Table 7.1: Individual runs characteristics with Combined3 analysis profile for Mrk 501 flare analysis with the number of events in the ON and OFF regions and the observing time corrected for instrument dead time for each run. The full observation numbers are given on the last line.

Spectrum and light curve

For the energy spectrum reconstruction, the forward folding technique presented in Section 1.5.1 is used. The best spectral hypothesis for the Mrk 501 data set is a curved power law function defined as:

Φ(E) = Φ 0 E E 0 -α-β log E E 0 (7.1)
where E 0 is the reference energy, Φ 0 is the differential flux at the energy E 0 , α is the spectral index and β the curvature parameter. The result of the fit is presented in Figure 7.4. The threshold energy is defined as the minimum energy to get at least 85% of the maximum acceptance value for the corresponding analysis. With the ZA observation [63 • -65 • ] and the Combined3 analysis, the threshold energy is 1.3 TeV. So, to ensure an optimal quality of the data, only events above this threshold energy are considered. For the highest energy only events below 20 TeV are used, at the last significant point in the spectrum.

For the light curve, a time binning of 7 minutes (420 s) is adopted which gives 4 bins per runs for a total of 28 minutes. This binning is chosen to avoid having a bin partially outside the observation time window. If not taken into account, this can lead to bins with very low flux value due to a lack of events. The usual light curve reconstruction in H.E.S.S. is a log-likelihood method as detailed in Section 1.5.2, but this method cannot be applied to this data set. Indeed, because of the very low number of background events, the log-likelihood method does not have enough OFF events to evaluate the flux in a given time bin. The simple statistics method, described in Section 1.5.2 is used instead to circumvent the problem. The light curve for the full energy range is shown in Figure 7.5.

True energy (TeV

Search of LIV signatures

In order to search for an LIV effect in the data, the ML is used to measure the LIV time delay, described in Chapter 6. However, to apply the ML method, one assumption is required which consists of neglecting the presence of any intrinsic time delay. The general PDF expressed in Equation 6.5 does not take into account intrinsic delays and would require a complete modeling of the Mrk 501 flare in order to have a precise determination of intrinsic delay and include it in the PDF for the ML.

Template and fit region definition

To apply the ML method, an energy cut has to be applied in order to divide the data between the template and the likelihood region. The total number of events used for the analysis, between the minimum energy (1.3 TeV) and the maximum (20 TeV), is 1435. The energy cut E cut is taken at 3.25 TeV splitting the data with 773 events in the template region (1.3 -3.25 TeV) and 662 events in the likelihood region (3.25 -20 TeV).

For the Γ(E) function, the energy spectrum is fitted in the likelihood region. To simplify the computation of the PDF and remove the energy resolution, a first attempt is done with a power law function. The resulting parameterization is presented on Figure 7.6 and shows a decent representation of the high energy spectrum. To determine F (t) for the PDF, the events in the template region are used. Because of the ZA variations between the runs (Table 7.1 last column), the acceptance variations cannot be neglected. In order to take them into account easily, the light curve is fitted instead of the time distribution. Different shapes were tested like a simple Gaussian function (Figure 7.7a) or a double Gaussian function (Figure 7.7b). The latter result is found to be the best parameterization function to represent the low energy events, based on a χ 2 /ndf comparison. As a remark, the E cut value is usually defined as the median energy which is for this data set E med = 3 TeV. But for the determination of F (t), the parameterization result were not stable enough. Indeed, different parameterizations obtained with E cut = 3 TeV and E cut = 3.25 TeV in different energy bands are presented in Figure 7.8 and show that E cut = 3.25 TeV provide a more robust estimation of F s (t). Indeed, the presence of the second pulse is not clear using the median energy and seems to appear only if E > 1.5 TeV. For this reason, the higher E cut value is taken in order to increase the number of events used for the parameterization and provide a more robust estimation of F (t). Such variations of the light curve parameterization obtained with E cut = 3 TeV can indicate the presence of intrinsic effect at energies E < 3.25 TeV.

Application of the Maximum Likelihood

As the data set has a template region at TeV energies, the template correction is used. The different terms used in the PDF, defined in Equation 6.18, for the ML are summarized in Table 7.2. The log-likelihood function, defined in Equation 6.8, is computed using the events time and energy in the likelihood region. The functions -2∆ log(L) are presented for the linear (Figure 7.9a) and quadratic (Figure 7.9b) LIV effect. The ML method results on the LIV parameter are:

τ best 1 = -8.2 +20 -19 s TeV -1 , τ best 2 = -0.6 +1.4 -1.2 s TeV -2 .
Both linear and quadratic results are compatible with no LIV time delay. As a result, upper limits on the Quantum Energy scale are derived. However, the associated errors evaluated through the likelihood function originate from the PDF definition. In that sense, the errors obtained are model dependent. Beside, the data only allows for one realization and thus, to obtain a more robust evaluations of the statistical uncertainties, the toy MC tools presented in Chapter 6 is used to generate more realization of the data set.

Monte Carlo simulations and statistical uncertainties

With the toy MC tool, a large number of simulated data set are generated, reproducing the Mrk 501 observation. Analyzed as the real data, these simulations allow for the determination of statistical uncertainties on the reconstructed of the LIV parameter. The simulation parameters used for the simulation are defined in Table 7.2. To determine the statistical uncertainties, simulations with no injected lag are done, which is in agreement with the measurement on data. As the τ best n distribution are not perfectly Gaussian, the statistical uncertainties are then estimated by the difference between the τ best n median value and the upper and lower limits median values. The upper and lower statistical uncertainties for linear and quadratic effect yield:

σ stat,1 ≈ +22 -20 s TeV -1 , σ stat,2 ≈ +1.

Systematic uncertainties studies

To evaluate the overall systematic uncertainties, individual estimation of each contribution is done to deduce the average uncertainties. For one specific systematic effect, a dedicated scenario is generated with the MC tool. The associated systematic uncertainties ∆τ n,i are then estimated by looking at the shift of the median value of the upper and lower limits at 68% confidence level compared to the ones from the scenario investigated for statistical uncertainties. Some systematic effect are estimated directly from the data.

Maximum likelihood calibration

The ML calibration is determined by generating MC simulation as close as possible to the data parameters (Table 7.2), injecting different values of the LIV parameter. The procedure to produce the calibrations plots is the same as the one in Section 6.3.4. Figure 7.11 presents the two calibration curves for linear and quadratic LIV effect. The systematic uncertainties related to the ML calibration are estimated with the variations of the median of the upper and lower limits at 68% of τ inj n from Figure 7.10 deduced by calibration curve. The results yield:

∆τ 1,calib < +5.5 -2.8 s TeV 1 , ∆τ 2,calib < +0.4 -0.5 s TeV 2 . ] - 1
[s.TeV

Parameterization procedure

As H.E.S.S. Collaboration (2011) pointed out, one of the most important source of systematic uncertainties is the determination of F (t) and Γ(E), used for the ML, via a fitting procedure.

To investigate the associated error on the fitting procedure, simulated data sets are generated as close as possible to the data using the function F (t) and Γ(t) but for the analysis, in the PDF, each parameter of these two functions are randomly chosen, from Gaussian distributions. These distributions are centered on the fitted parameters values and their widths correspond to the standard deviation. This randomization procedure correspond to the errors made on the parameterization of the light curve and the spectrum.

The distributions of τ best n following this procedure are presented in Figure 7.12 for linear and quadratic LIV effect. These distributions are wider than the standard one in Figure 7.10. The resulting systematic uncertainties estimated via the median values of the upper and lower limit at 68% are: ∆τ 1,fit < +4.

Photon list determination

The determination of the photon list is done via the parameters of the H.E.S.S. analysis. Changing these parameters can lead to a different photon list which can be used to estimate the systematic uncertainty related to the H.E.S.S analysis. To do so, a second photon list is obtained using a different analysis configuration. The shift on variations of the confidence level interval of the reconstructed LIV parameter give then the estimation of the systematic error. The Stereo VeryLoose (SVL) analysis is chosen, which considers only stereoscopic events from at least two out of the five telescopes. In particular, this configuration removes the events reconstructed with only CT5 thus leading to an increase of the threshold energy and considers different the selection cut for the data analysis.

For the ML method, since the spectrum and the light curve are used in the PDF, they can also be used for the SVL analysis. Indeed, the light curve and the spectrum are properties of the source which do not depend on the analysis configuration in opposition to the photon list. So, only events in the likelihood region from the SVL event list are used for the ML procedure, resulting in the log-likelihood functions in Figure 7.13. Both results for linear and quadratic LIV effect are marginally compatible with zero lag at 68% confidence levels. The systematic uncertainties associated with the photon list determination are evaluated by comparing the upper and lower limits at 68% obtained with log-likelihood functions between the main (Combined3 ) and the SVL analysis and give: ∆τ 1,list < +10. 6 -8.2 s TeV -1 , ∆τ 2,list < +0. 7 -0.6 s TeV

Instrument response Function

The acceptance, the energy resolution and energy bias for the Mrk 501 observation conditions with the Combined3 analysis are shown in Figure 7.14. Their evaluation is handled by the H.E.S.S. analysis software and they are determined with the best accuracy possible leading to small uncertainties on their values. Acceptance and energy resolution are included in the light curve and the spectrum in the PDF. An error on the evaluation of the acceptance would lead to a wrong determination of the light curve and the spectrum with respect to the true ones. Thus the systematic uncertainties related to the acceptance is already taken into account via the systematic on the parameterization procedure. The energy resolution leads to errors already investigated in Section 6.4.1, and turned out to be relatively small on the LIV parameter reconstruction, especially in case of no delay in the data. Indeed, in case of no energy dependent delay, even if the energy is misreconstructed, the result is still going to be compatible with no energy-dependent delay. As a result, since the data shows no significant delay, the energy resolution systematic uncertainties is neglected.

For the energy bias, a dedicated study is done to evaluate its contribution to systematic uncertainties. Using the photon list, the energy of each events is unbiased by the values obtained by the H.E.S.S. analysis software (Figure 7.14c). This new photon list is then used in the ML analysis. The Figure 7.15 shows the resulting likelihood functions for linear and quadratic LIV effects. From theses results, the bias contribution to systematic uncertainties is taken as the shift of the 68% confidence level upper and lower limits from the main analysis and yield:

Systematic studies: Background contamination

For the LIV analysis, the background was neglected due to the high signal over background ratio. However, a small contamination is possible which can induce a misreconstruction of the LIV parameter. In order to evaluate how a small background contamination affect the result, MC simulations with an injection of background photons are done. The signal over background ratio in the likelihood region [3.25 -20] TeV in the data is 34. With a number of 663 events used for the log-likelihood function, 663/34 ∼ 20 background events are injected. These events, on the contrary to signal event, are generated following a flat light curve and a power law spectrum with index α bck = 2.7.

Since background photons do not show any variability, they should induce a worse estimation of the LIV parameter than without these photons. Figure 7.16 presents the distributions of reconstructed τ best n obtained with those simulations for linear and quadratic LIV effect. As expected, the background injection do not significantly impacts the LIV parameter measurement, compared to the results of Figure 7.10, due to the small number of background events. The associated systematic uncertainties are estimated as : ∆τ 1,bck < +0.8 -0.1 s TeV -1 , ∆τ 2,bck < +0.

Total systematic uncertainties

All the individual systematic uncertainties contributions ∆τ n,i are used to determine the upper and lower systematic error at 68% confidence level. Its value is computed by summing quadratically all the contributions:

σ syst,n = i ∆τ 2 n,i . (7.2)
All values obtained are summarized in Table 7.3 for n = 1 and n = 2. The resulting systematic uncertainties are then: σ syst,1 < +13 -11 s TeV -1 , σ syst,2 < +0.9 -1.0 s TeV -2 .

Compared with the statistical uncertainties obtained in Section 7.2.3, the systematic uncertainties correspond to 67% and 38% of the statistical ones for respectively the linear and quadratic LIV effect. ∆τ 1,i (s TeV -1 ) ∆τ

Results on the Quantum Gravity

Limits on the LIV time lag

The Maximum Likelihood method applied on data together with the statistical and systematic studies give the following results on the LIV parameter τ n with Mrk 501 data: with the errors corresponding to the 68% confidence level for statistical and systematic uncertainties. The measurement is compatible with an absence of any linear or quadratic LIV effect and only limits on the energy scale of QG can be deduced.

The limits on τ n at a 95% confidence level are computed considering the two LIV cases presented in Chapter 2: subluminal and superluminal. The subluminal case corresponds to a positive value of τ n with high energy photon arriving after low energy ones (positive errors are taken). The superluminal case corresponds to the negative τ n value (negative errors are considered). The two cases are tagged as τ n, + for subluminal and as τ n, -for superluminal. The resulting τ 95% n, ± limit is computed as:

τ 95% n, ± < τ best n ± 2 σ 2 stat + σ 2 syst . (7.3)
For the Mrk 501 data set, the limits on the LIV parameter for the linear and quadratic effects in both subluminal and superluminal cases, are: τ 95% 1, + < 42.9 s TeV -1 subluminal limit τ 95% 1, -< 53.9 s TeV -1 superluminal limit τ 95% 2, + < 3.1 s TeV -2 subluminal limit τ 95% 2, -< 4.0 s TeV -2 superluminal limit

These values can finally be used to compute the associated lower limits on E 95% QG,n at 95% CL for the subluminal and superluminal cases.

Limits on E QG

From to Equation 6.2, the lower limit on E 95% QG,n is expressed as a function of τ 95% n, ± as:

E 95% QG,n = 1 + n 2H 0   1 τ 95% n, ±   n
κ n (z), (7.4) With a redshift of z = 0.034, the final limits on the QG energy scale E 95% QG,n obtained with the Mrk 501 data set for a linear and quadratic LIV effect in both the subluminal and superluminal cases are:

E 95% QG,1 =    3.

Discussion on the results

The results of the present LIV study on Mrk 501 provide competitive lower limits on the QG energy scale with respect to previous results. Like other AGN lower limits on E QG,n , these results are still one order of magnitude below E P for the linear lower limits whereas for the quadratic ones, E P is still out of reach due the sensibility of current instruments. As previous results (listed in Table 2.1) were only given for the subluminal case, only a comparison between subliminal lower limit from AGN will be discussed. A direct comparison of all the subluminal lower limits obtained with IACT experiment on AGN flare is given in Figure 7.17 and the different characteristics of all the sources considered are presented Table 7.4.

The lower limit on linear LIV effect E 95% QG,1 with Mrk 501 is at the level of the one obtained with PG 1553+113 explained by the lower redshift of Mrk 501 while it has a larger number of photons and at higher energy. The source results obtained with H.E.S.S. are slightly better than the MAGIC one mostly explained by the overall better characteristics of the H.E.S.S. data. However, the result is less constraining than the one of PKS 2155-304 one because of the much smaller number of events even if the energetic are in favor or Mrk 501. The quadratic lower limit benefits more from the high energy photons up to 20 TeV providing from the Mrk 501 flare the best constrain on E 95% QG,2 even higher than lower limit obtain from the data of the exceptional flare of PKS 2155-304 in 2006.

The results obtained with this data set could be improved with a more precise determination of the systematic uncertainties. The use of MC simulations to estimate each individual contribution on systematic uncertainty to deduce the total systematic errors leads to an overestimation of its value. Indeed, the uncertainties of each contribution were taken conservatively, and also some effects could have been counted twice in the process. Besides, the ML procedure can be improved, such as taking into account the energy resolution in the PDF which allows for a better estimation of the spectrum in the whole energy range for the function Γ(E) instead of using a single power law in the high energy range. Another potential improvement concern the ML procedure for the LIV parameter reconstruction in order to potentially reduces one of the main uncertainties concerning the parameterization of the light curve. It consists of adding the positions of the pulses in the ML method which are the most important parameter for the τ n evaluation. In one hand, it allows for the use of all events in the ML, instead of splitting the data sample in two parts for the determination of F s (t). On the second hand, it provides a better estimation of the error of the position of the peaks thanks to the likelihood function, but also to take into account the correlation between the peak positions and the LIV parameter. However, such method induces a more complex evaluation of the likelihood function, which for the case of the Mrk 501 data set involves a 3 parameter likelihood function. It also requires to assess the global shape of the light curves, and the other parameter of F s (t) (width of each pulse and relative amplitude between the two pulses) if they are not also included in the likelihood function.

Conclusion

In this chapter, the search for Lorentz Invariance Violation signatures using the time of flight approach was done on the observation of the Markarian 501 flare which occurred in the night of June 23-24 2014. The source was observed for roughly 2 hours, with high zenith angle which results in a high energy threshold of 1.3 TeV for the data set. The associated energy spectrum is extending up to 20 TeV. The light curve presents a short variability of a few minutes, and revealed the presence of two Gaussian pulses above 1.5 TeV. The 1435 events recorded at energies between 1.3 and 20 TeV were used to search the presence of any Lorentz Invariance Violation signature, i.e. the presence of a linear or quadratic time delay in function of energy.

The method of Maximum Likelihood used revealed the absence of any significant energydependent time-delay for both linear and quadratic Lorentz Invariance Violation effects. Thus, this result allowed for the establishment of lower limits on the Quantum Gravity energy scale. A dedicated study was performed using the toy Monte-Carlo simulations to determine statistical and evaluate systematic errors. Ultimately, the lower limits at a 95% confidence level were computed for subluminal and superluminal cases which provide competitive lower limits compared to others obtained with Active Galactic Nucleus. In addition, the quadratic lower limits is found to be the best one obtained with Active Galactic Nucleus flare up to now. This analysis is part of a paper concerning the observation of this Mrk 501 flare and the search for Lorentz Invariance Violation signatures for both time-of-flight and cross-section modification approach. However, in order to fully exploit this data set, one should take a careful look on the light curve behavior below 1.5 TeV, as seen in Figure 7.8, which hints for potential intrinsic effect.

Conclusion

In this thesis, the topic of spectral time delays was investigated in the framework of blazar modeling and of the search of Lorentz Invariance Violation (LIV) signatures. The combination of these two studies turns out to be necessary if one wants to constrain Quantum Gravity (QG) models using LIV in the time-of-flight approach. Indeed, in case of measurement of a significant time delay in a blazar flare, one needs to ensure either that this delay is only coming from LIV in order to provide a constraint on the energy scale E QG or to estimate the contribution from intrinsic effects using blazar modeling.

The time dependent model developed and used for this work was based on the evolution of the electrons responsible for the flare through their synchrotron-self-Compton emission, described with a transfer equation. Using reasonable assumptions, an analytical solution of this equation was used in order to generate flares, within a specific domain of validity. The model showed the presence of intrinsic time delays. After studying its parameters, it revealed two main regimes depending on the mechanism driving the evolution of the high energy electron population when the light curves peak (acceleration or cooling effects). The presence of one of the regimes in a blazar flare observation can provide precious information on the mechanisms generating a blazar flare. This will be useful in case of time delay detection and will provide additional constraints on the in situ parameters of the source. Besides, the delays obtained with the model showed specific characteristics which can be used to constrain QG models. The temporal dependency of the intrinsic delays was found to be different from the one expected by a LIV effect and can be used as an indication of the presence of intrinsic effects. Indeed, a LIV delay is not expected to show any kind of evolution with time since it affects photons in the same way throughout their propagation. In opposition, the intrinsic delay evolves with time due to the different energy dependent mechanisms resulting in the generation of a blazar flare. Besides, the energy dependency of the intrinsic time delay at GeV-TeV energies was found to present a power index α ∈ [0.45 -0.85] for cases with a significant change of the time delay with the energy. This property can be used in order to test specific QG models which predict an energy dependent LIV delay with an index different from the intrinsic ones.

The flares investigated in this thesis were based on a minimal scenario able to generate a flare, emphasizing the dominant processes. It provided a comprehensive understanding of the mechanisms leading to intrinsic time delays. However, the analytic solution used was not able to describe alternative complex scenarios such as electrons dominated by inverse Compton cooling or the presence of particle escape. On one hand, the model can be improved by extending the analytic solution to scenarios considering external inverse Compton emission. Such a process is compatible with the analytic solution and would allow the model to generate stronger emission at γ-ray energies. On the other hand, the development of a full numerical resolution for the transfer equation would allow for the improvement of the domain of application of the model and also to include more processes, more specific acceleration processes and additional cooling effects. Such a resolution of the transfer equation would provide a wider description of time dependent emission of blazars to model all kinds of flares.

Also in this thesis, the measurement of LIV energy dependent time delay was carried out for the flare of Markarian 501 which occurred during the night of June 23-24 2014. The maximum likelihood method widely used for the search of LIV signatures, was presented and tested with Monte-Carlo simulations. These simulations allowed for the development of a new template correction for the LIV analysis of the flare. This correction also applies in general to all data sets presenting a high energy threshold, where the LIV effect cannot be neglected in the low energy range of the data to determine the light curve shape. Without this correction, the method is not well calibrated and the estimation of the LIV parameter is biased. This is generally the case for high zenith angle observations such as the one on the Markarian 501 flare. The analysis of this flare revealed an energy spectrum extending from 1.3 TeV up to 20 TeV and a light curve with variability of a few minutes showing a two peak structure above 1.5 TeV. After a selection of 1430 events in the energy range 1.3-20 TeV, the maximum likelihood was applied for the search of a linear and a quadratic LIV effect and showed no significant energy dependent time delay. These results allowed for the establishment limits at 95% confidence level on the QG energy scale E QG for subluminal and superluminal cases. For subluminal case, the limits are E 95% QG,1 > 3.63 × 10 17 GeV and E 95% QG,2 > 8.79 × 10 10 GeV and for the superluminal case, E 95% QG,1 > 2.89 × 10 17 GeV and E 95% QG,2 > 7.66 × 10 10 GeV. The limits on E QG,2 were the best that were obtained with AGN up to now.

However, during the analysis of the Markarian 501 data, a hint of intrinsic effects at energies E < 1.5 TeV was found with a variation of the light curve shape with different energy selections. This offers a potential use of the time dependent model to estimate a contribution from intrinsic effects and deserves further analysis. This estimation can be injected in the mock-up function used for the maximum likelihood in a similar way than the template correction was. This would allow the Probability Density Function to account for an intrinsic time delay contribution in the estimation of the LIV delay. One can also directly use the light curve shape produced by the model in the mock-up function for the LIV parameter measurement. However, to ensure an unbiased determination of the γ-ray light curves, one should not use γ-ray data to model the source. Such a procedure would require multiple simultaneous observations of a single flare at different energies. For example, one could try to constrain the magnetic field strength from very high angular resolution radio data, and the acceleration parameters could be determined to reproduce the flux increase on the MeV-GeV light curves. Then the resulting intrinsic light curves at GeV-TeV could be used in the maximum likelihood method to search for LIV with H.E.S.S. data.

In addition, the maximum likelihood method used to measure a linear and a quadratic LIV delay could be adapted to analyze different energy dependencies. One way would be to assume specific α values obtained from the modeling. Another way would be to estimate α directly with the maximum likelihood. This would provide the most probable energy dependent time delay corresponding to the data. The reconstructed energy dependent time delay index α can then appear as an indicator to test some specific QG models and could be compared with the blazar modeling prediction to determine a possible contribution from intrinsic effect for this data set. If no delay is measured with any kind of energy evolution, strong constraints could be put on the flare model. Such a situation should take place between the two time delay regimes where acceleration and cooling effects counterbalance themselves.

The tools and results presented in this thesis contribute to the scientific preparation of the next generation of Cherenkov telescopes. The future Cherenkov Telescope Array (CTA) observatory [START_REF] Acharya | Introducing the CTA concept[END_REF] will provide a higher sensitivity and a much larger number of blazar flare detections than current instruments. CTA will provide better quality data compare to current generation of instrument in order to search for LIV signatures but also to study intrinsic time delays. The larger number of flare observations expected with CTA will increase the probability of detecting significant time delays and probably lead to several positive detections. The tool provided in this work to describe intrinsic delays in the context of LIV study using photon time of flight will then become a necessity to be able to get robust constraints on the QG energy scale as well as to get new constraints on the modeling of blazar flares using the time delay information.

As a final remark, modeling blazar is not the only solution to deal with intrinsic delays. With the hopefully increased number of detection with CTA, another approach to search for LIV signatures is to exploit the redshift dependency of the LIV delay with a population study. Indeed, the LIV delay increases when the propagation distance for photons becomes larger. Thus, using a population of several AGN at different redshifts, a combined estimation of E QG should not be affected by intrinsic effects. Indeed, at first glance, these effects are not expected to depend significantly on the redshift. In this regard, during this thesis, a small contribution was done on a joint effort developed between the three major imaging atmospheric Cherenkov telescope experiments (H.E.S.S., MAGIC and VERITAS) to share data in order to operate such a population study and prepare for CTA. This combination aims at including blazar flare observations but also pulsars and GRB observations. Using simulations, the combination of a few blazar flares and one pulsar has shown to provide better limits on E QG than individual ones obtained from each source (Nogués et al., 2017).

Appendix A From the modified dispersion relation to the energy dependent velocity of photons

From the MDR shown in Equation 2.4, one can deduce the resulting group velocity of photon v g (E). This velocity is expressed as:

v g (E) = ∂E ∂p . ( 7.5)
To do such computation, the MDR first needs to be simplified. Only considering the leading order n, it can be written as: (7.6)

E 2 p 2 c 2 ×   1 ± E E QG n   ,
E pc   1 ± E E QG n   1 2 . (7.7)
With a first order Taylor expansion1 considering (E/E QG ) n → 0, the MDR becomes:

E = pc   1 ± 1 2 E E QG n   . (7.8)
Then v g (E) can be derived from Equation 7.8: (7.10) Considering that E E QG , one can deduce pc ≈ E which provides: (7.11)

v g (E) = ∂E ∂p = c   1 ± 1 2 E E QG n   ± pc 2E n QG ∂E n ∂E ∂E ∂p , (7.9)   1 ∓ pc 2E n QG E n-1   ∂E ∂p = c   1 ± 1 2 E E QG n   .
  1 ∓ n 2 E E QG n   ∂E ∂p = c   1 ± 1 2 E E QG n   ,
∂E ∂p = c   1 ± 1 2 E E QG n     1 ∓ n 2 E E QG n   -1
.

(7.12)

Using another Taylor expansion of the last factor of Equation 7.12, considering 1 2 n E/E QG n → 0, gives: (7.13) (7.14) Finally, neglecting the higher order term E/E QG 2n provides the expression for the energy dependent velocity of photons: 7.15)

∂E ∂p = c   1 ± 1 2 E E QG n     1 ± n 2 E E QG n   ,
∂E ∂p = c   1 ± n + 1 2 E E QG n + n 4 E E QG 2n   .
v g (E) = ∂E ∂p = c   1 ± n + 1 2 E E QG n   (

xii List of figures 1 . 1

 11 Heitler models of atmospheric showers. . . . . . . . . . . . . . . . . . . . . . . . . 1.2 The H.E.S.S. site in 2012. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3 A CT5 drawer with 16 photo-multiplier tubes. . . . . . . . . . . . . . . . . . . . 1.4 Charge distribution in one pixel during a single photo-electron run. . . . . . . . . 1.5 Comparison between muon rings observed by H.E.S.S. and an analytical model. . 1.6 Shower image in the H.E.S.S. camera for γ and hadron-like event. . . . . . . . . . 1.7 Scheme presenting the Hillas parameters. . . . . . . . . . . . . . . . . . . . . . . 1.8 Distribution of mean scale width with a blank field data set and simulated γ and hadron-like events. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.9 Definition of the OFF regions for the Reflected Background method. . . . . . . . 2.1 Infrared observations in H band of lensed Active Galactic Nuclei from CASTLES survey. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Photography and rotation curve of Andromeda galaxy. . . . . . . . . . . . . . . . 2.3 Photon distributions of GRB 090510. . . . . . . . . . . . . . . . . . . . . . . . . . 2.4 Phasograms of Vela pulsar at γ-ray energies from H.E.S.S. and Fermi-LAT. . . . 2.5 Unification scheme of AGN. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.6 Light curve of the 2006 PKS 2155-304 Big flare. . . . . . . . . . . . . . . . . . . . 3.1 Radio and X-ray observation of M87 inner jet. . . . . . . . . . . . . . . . . . . . 3.2 Observations of superluminal motion in the inner jet of PKS 1510-089. . . . . . . 3.3 Scheme illustrating the observed superluminal motion of the jets in blazars. . . . 3.4 Revised and initial blazar sequences. . . . . . . . . . . . . . . . . . . . . . . . . . 3.5 Scheme illustrating inverse Compton scattering geometry. . . . . . . . . . . . . . 3.6 EBL spectral energy distribution estimations from various studies. . . . . . . . . 3.7 Schematic view of the second order Fermi acceleration mechanism. . . . . . . . . 3.8 Schematic view of the first order Fermi acceleration mechanism. . . . . . . . . . . 3.9 Schematic view of a magnetic reconnection site example. . . . . . . . . . . . . . . 3.10 Leptonic and hadronic SED modeling of the blazar Markarian 421. . . . . . . . . 3.11 Synchrotron-self-Compton and external inverse-Compton spectral energy distribution modeling for the blazar AP Librae. . . . . . . . . . . . . . . . . . . . . . . 4.1 Inverse Compton light curves from the model of Sokolov et al. (2004). . . . . . . 4.2 Delayed X-ray light curves from time-dependent model of Mrk 421 1998 flare. . . 4.3 Time evolution of the parameter η(t) from computational of the emission for a single flare. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4 Electron spectrum and spectral energy distribution evolution for the reference case. 4.5 Normalized light curves for the reference scenario. . . . . . . . . . . . . . . . . . xiii 4.6 Simulated light curves with non-constant width and time delay measurement with CCF and PPM. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.7 Simulated light curves with constant width and time delay measurement with CCF and PPM. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.8 Time delay with respect to the energy for the reference scenario. . . . . . . . . . 5.1 Electron spectrum evolution for the reference case without adiabatic expansion. . 5.2 Spectral energy distribution for the reference case without adiabatic expansion. . 5.3 Normalized light curves for the reference case without adiabatic expansion. . . . 5.4 Time delay with respect to the energy for the reference case with and without adiabatic expansion. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.5 Evolution of the electron spectrum energy cut-off for the reference case without adiabatic expansion. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.6 Evolution of the energy dependent time delay for various initial magnetic field strength values. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.7 Normalized light curves for the reference case with B 0 = 90 mG and without adiabatic losses. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.8 Evolution of the energy dependent time delay with various magnetic field evolution index for two initial magnetic field values. . . . . . . . . . . . . . . . . . . . 5.9 Evolution of the energy dependent time delay with various Doppler factor for two initial magnetic field values. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.10 Evolution of the energy dependent time delay with various initial acceleration amplitude for two initial magnetic field values. . . . . . . . . . . . . . . . . . . . 5.11 Evolution of the energy dependent time delay with various acceleration evolution index for two initial magnetic field values. . . . . . . . . . . . . . . . . . . . . . . 5.12 Evolution of the energy dependent time delay with various electron spectrum index for two initial magnetic field values. . . . . . . . . . . . . . . . . . . . . . . 5.13 Evolution of the energy dependent time delay with various initial electron spectrum energy cut-off for two initial magnetic field values. . . . . . . . . . . . . . . 5.14 Time delay with respect to the energy for the reference case with adiabatic expansion. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.15 Times delay in the GeV-TeV energy range for the two initial magnetic field values. 5.16 Comparison between LIV and intrinsic time delay in the GeV-TeV energy range for the two initial magnetic fiend values. . . . . . . . . . . . . . . . . . . . . . . . 5.17 Times delay temporal evolution for the reference case without adiabatic expansion and B 0 = 65 mG. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.18 Times delay temporal evolution for the reference case without adiabatic expansion B 0 = 90 mG. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.19 Times delay temporal evolution for simulated light curves delays by an average LIV effect. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.1 Time and energy distributions for one realization of the simple scenario. . . . . . 6.2 Log-likelihood function for one realization of the simple scenario. . . . . . . . . . 6.3 Reconstructed LIV parameter distribution for a 1000 simulations of the simple example. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.4 Calibration plots for the simple scenario. . . . . . . . . . . . . . . . . . . . . . . . 6.5 Calibrations plots for the energy resolution scenario. . . . . . . . . . . . . . . . . 6.6 Calibration plots for the fitted template scenario. . . . . . . . . . . . . . . . . . . xiv 6.7 Calibrations plots for the lagged template and the corrected template scenarios for linear LIV effect. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105 7.1 Illustration of the high zenith angle observation. . . . . . . . . . . . . . . . . . . 109 7.2 Θ 2 histogram of Mrk 501 flare observations. . . . . . . . . . . . . . . . . . . . . . 109 7.3 Significance map of Mrk 501 flare observations. . . . . . . . . . . . . . . . . . . . 109 7.4 Energy spectrum of Mrk 501 flare observations. . . . . . . . . . . . . . . . . . . . 111 7.5 Light curve of Mrk 501 flare observations. . . . . . . . . . . . . . . . . . . . . . . 112 7.6 Energy Spectrum of Mrk 501 above 3.25 TeV. . . . . . . . . . . . . . . . . . . . . 113 7.7 Template light curve of Mrk 501 parameterized with a simple and double Gaussian functions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114 7.8 Template light curve parameterization results with different energy ranges . . . . 114 7.9 Normalized log-likelihood functions obtained from Mrk 501 data. . . . . . . . . . 115 7.10 Monte Carlo τ best n distributions with simulations of Mrk 501 data set. . . . . . . 116 7.11 Calibration plots with Mrk 501 Monte Carlo simulations. . . . . . . . . . . . . . 117 7.12 Monte Carlo τ best n distributions with Mrk 501 Monte Carlo simulations for the parameterization systematic uncertainties study. . . . . . . . . . . . . . . . . . . 118 7.13 Normalized log-likelihood functions obtained from Mrk 501 data with StereoVery-Loose analysis configuration. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119 7.14 Instrument response functions for the Mrk 501 observations conditions with Com-bined3 analysis configuration. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120 7.15 Normalized log-likelihood functions obtained from Mrk 501 data for the energy bias systematic uncertainties study. . . . . . . . . . . . . . . . . . . . . . . . . . . 121 7.16 Monte Carlo τ best n distributions with Mrk 501 Monte Carlo simulations for the background systematic uncertainties study. . . . . . . . . . . . . . . . . . . . . . 122 7.17E QG,n subluminal limits obtained different blazar flares. . . . . . . . . . . . . . . 125 xv List of tables 1.1 Analysis cuts for the Combined3 and StereoVeryLoose analyses . . . . . . . . . . 2.1 List of QG limits obtained with various types of sources and various instruments 4.1 Simulation input parameters for the reference case. . . . . . . . . . . . . . . . . . 5.1 Energy dependent time delay amplitude ξ and power index α for various initial magnetic field strengths in the GeV-TeV energy range. . . . . . . . . . . . . . . . 5.2 Energy dependent time delay amplitude ξ and power index α for all investigated parameters for the two initial magnetic field strength values. . . . . . . . . . . . . 6.1 Simulation input parameters of the simple scenario. . . . . . . . . . . . . . . . . . 6.2 Scenario list for the MC simulation studies. . . . . . . . . . . . . . . . . . . . . . 7.1 Individual run characteristics with Combined3 analysis for Mrk 501 flare analysis. 7.2 Analysis parameters for the LIV analysis of Mrk 501. . . . . . . . . . . . . . . . . 7.3 Systematic uncertainties for each contributions on the τ n measurement. . . . . . 7.4 Active Galactic Nuclei characteristics. . . . . . . . . . . . . . . . . . . . . . . . . xvii

Contents 1 . 1 1 . 4 1 . 5

 111415 Detection of gamma rays with ground based telescopes . . . . . . . 4 1.1.1 Extensive air shower . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 1.1.2 Heitler models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.1.3 Detection techniques . . . . . . . . . . . . . . . . . . . . . . . . . . .E.S.S. telescopes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 1.2.2 H.E.S.S. cameras . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 1.2.3 Trigger system . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 1.2.4 Data calibration . . . . . . . . . . . . . . . . . . . . . . . . . . . . Hillas model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 1.3.2 The semi-analytical model reconstruction . . . . . . . . . . . . . . . . . 14 Data analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 1.4.1 Analysis configuration . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 1.4.2 Acceptance and energy resolution determination . . . . . . . . . . . . . 16 1.4.3 Background subtraction . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 Spectral analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 1.5.1 Forward folding technique . . . . . . . . . . . . . . . . . . . . . . . . . . 18 1.5.2 Light curve determination . . . . . . . . . . . . . . . . . . . . . . . .

Figure 1 . 1 :

 11 Figure1.1: Heitler models of atmospheric showers in case of electromagnetic shower (left) and hadronic shower (right). For hadronic showers, π 0 decay creating two γ rays leading to electromagnetic showers described by the left shower. R is the radiation length after which a particle loses half its energy.

Figure 1 . 2 :

 12 Figure 1.2: The H.E.S.S. site in 2012. Credits from H.E.S.S. Collaboration, Frikkie van Greunen.

Figure 1

 1 Photomultiplier tubePhotonis XP-29600

Figure 1 . 4 :

 14 Figure 1.4: Charge distribution in one pixel during a single photo-electron run. The figure is taken from Aharonian et al. (2004a).

Figure 1 . 5 :

 15 Figure 1.5: Comparison between muon rings observed by H.E.S.S. (left) and an analytical model (right). The figure is taken from Chalme-Calvet et al. (2014).

  (a) γ-like event. (b) Hadron line event.

Figure 1 . 6 :

 16 Figure 1.6: Shower image in the H.E.S.S. camera after cleaning, for γ (top) and hadron (bottom) like events.

•

  the ellipse half width W • the total amplitude of the shower image • the center of gravity of the shower image d • the distance between the center of gravity of the shower image and the expected source position in the focal plane of the camera • the direction of the shower in the camera defined with the angle ϕ and the angular distance α as shown in Figure 1.7.

Figure 1 . 7 :

 17 Figure 1.7: Scheme presenting the Hillas parameters.Figure taken from Garrigoux (2015).

Figure 1 . 8 :

 18 Figure 1.8: Distribution of mean scale width with a blank field data set and simulated γ and hadron-like events. The data set used does not present any known source in the FoV. The variable M RSW stands here for the M SW . The vertical dashed lines represent the selection cuts for the γ-hadron separation at values -2 and 1. Figure taken from Aharonian et al. (2004b).

Figure 1 .

 1 Figure 1.9: Definition of the OFF regions for the Reflected Background method. The figure shows two observation positions corresponding to the camera center (yellow) with the same ON region (white) and their resulting OFF regions (red). Figure taken from Berge et al. (2007).

  Figure 1.9: Definition of the OFF regions for the Reflected Background method. The figure shows two observation positions corresponding to the camera center (yellow) with the same ON region (white) and their resulting OFF regions (red). Figure taken from Berge et al. (2007).

Figure 2 . 1 :

 21 Figure 2.1: Infrared observations in H band of lensed Active Galactic Nuclei from CAS-TLES survey. On each figure, the central spot corresponds to the foreground galaxy which acts as a gravitational lens. The other spots are the images of the lensed AGN located in the background. Images (a) and (c) also shows Einstein rings. The images are taken from https://www.cfa.harvard.edu/castles/.

Figure 2 . 2 :

 22 Figure 2.2: Photography of Andromeda galaxy (left) and its corresponding rotation curve (right). Each point on the galaxy photography corresponds to a measurement of rotation velocity at a given distance to the galaxy center on the rotation curve plot. Figures taken from Rubin & Ford (1970).

Figure 2 . 3 :

 23 Figure 2.3: Photons distribution of GRB 090510 in different energy bands, as measured by Fermi GBM and Fermi-LAT. The top four panels shows the time distributions of photons with a bin size of 0.01 s. The last panel represent individual photons energy versus time.The red vertical line corresponds to the trigger time of the observation, while the green ones are used for the spectral analysis on the paper from which the plot is extracted[START_REF] Ackermann | Fermi Observations of GRB 090510: A Short-Hard Gamma-ray Burst with an Additional, Hard Power-law Component from 10 keV TO GeV Energies[END_REF].

Figure 2 . 4 :

 24 Figure 2.4: Phasograms of Vela pulsar at γ-ray energies from H.E.S.S. (top) and Fermi-LAT (bottom) experiments. The grey shaded areas, bounded by the vertical lines represent the different pulse phases determined with the Fermi-LAT phasogram. The vertical dashed line on the top panel represents the level of background. The period of Vela pulsar is 0.089 s. Figure taken from H.E.S.S. Collaboration et al. (2018).

Figure 2 . 5 :Figure 2 . 6 :

 2526 Figure 2.5: Unification scheme of AGN based on the viewing angle for radio-loud (with jets) and radio-quiet sources (without jets). Figure taken from Beckmann & Shrader (2012).

  34

1 : 1 1 1

 1111 List of QG limits obtained with various types of sources and various instruments. The limits presented are the 95% confidence level for the subluminal case. The limits are shown for linear and quadratic term for each kind of sources. Characteristics of blazars . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1.1 Jets and superluminal motion . . . . . . . . . . . . . . . . . . . . . . . . 3.1.2 Doppler boosting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1.3 Blazar sequence . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Synchrotron emission . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.2 Inverse-Compton interaction . . . . . . . . . . . . . . . . . . . . . . . . 3.2.3 Hadronic emission processes . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.4 γ -γ interaction absorption . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.5 Extragalactic Background Light absorption . . . . . . . . . . . . . . . . Second order Fermi mechanism . . . . . . . . . . . . . . . . . . . . . . . 3.3.2 First order Fermi mechanism . . . . . . . . . . . . . . . . . . . . . . . . 3.3.3 Magnetic reconnection . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4 Modeling blazars spectral energy distribution . . . . . . . . . . . . . 3.4.1 Leptonic models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4.2 Hadronic and lepto-hadronic models . . . . . . . . . . . . . . . . . . . .

Figure 3 . 1 :

 31 Figure 3.1: Images of the M87 inner jet in radio at frequency 43 GHz (left) and in X-ray (right). The radio data were obtained with the Very Long Baseline Array and the X-ray with the Chandra X-ray Observatory. Figures taken from Wagner et al. (2009).

Figure 3 . 2 :

 32 Figure 3.2: Radio observations of superluminal motion in the inner jets of the TeV blazar PKS 1510-089 (z = 0.361) done with the Very Long Baseline Array. The contours represent different levels of intensity at frequency 43 GHz. The two panels correspond to two different knots observed at different epochs. The first knot is observed at an apparent velocity 24 ± 2 c while the second is at 21.6 ± 0.6 c (for the y-axis, 0.10 mas = 0.5 parsec). Figure taken from Marscher et al. (2010).

  Figure 3.3: Scheme illustrating the observed superluminal motion of the jets in blazars A knot is emitted with a relativistic speed from the central engine at a small viewing angle θ relatively to the observer. Due to the relativistic motion of the knot, the two emitted photon are observed with a smaller time difference with respect to the difference of their emission time leading to observation of superluminal motion.

Figure 3 . 4 :

 34 Figure 3.4: Revised (left) and initial (right) blazar sequences. The difference between the two sequences comes from the SEDs categorization in luminosity: the initial one uses radio luminosity while the revised one uses γ-ray luminosity. The equivalent energy range for the x-axis goes from 4.7 10 -8 eV to 470 TeV. Figure taken from Ghisellini et al. (2017).

Figure 3 . 5 :

 35 Figure 3.5: Scheme illustrating inverse Compton scattering geometry seen by an observer (left) and in the frame of the electron (right).

  Figure 3.6: EBL spectral energy distribution estimations from various studies. The red points corresponds to the results by H.E.S.S. Collaboration et al. (2017) obtained from the deconvolution of blazar spectra by EBL absorption.

Figure 3 . 7 :Figure 3

 373 Figure 3.7: Schematic view of the second order Fermi acceleration mechanism. A particle enters a magnetized cloud with a pitch angle θ in and an energy E in , is deflected and accelerated in the cloud and exits the cloud with a pitch angle θ out and an energy E out .

Figure 3

 3 Figure 3.9: Schematic view of a magnetic reconnection site example. The top and bottom magnetic field lines, carried by the motion of the plasma (illustrated by wide yellow arrows) reconnect together at the central X-point, and rearrange themselves with the magnetic lines on the left and right sides. Particles are accelerated by electric fields generated during the process.

Figure 3 .

 3 Figure 3.10: Leptonic and hadronic modeling of the blazar Markarian 421. The different contributions to the SED are specified in the legend. Figures taken from Zech et al. (2017).

Figure 3

 3 Figure 3.11: Synchrotron-self-Compton and external inverse-Compton spectral energy distribution modeling for the blazar AP Librae. The different components of the SED are specified in the legend. The figure is taken from Hervet et al. (2015).

  (a) Delay from cooling time scale. (b) Delay from internal light travel.

Figure 4 . 1 :

 41 Figure 4.1: Inverse Compton light curves from the model of[START_REF] Sokolov | Synchrotron Self-Compton Model for Rapid Nonthermal Flares in Blazars with Frequency-dependent Time Lags[END_REF] in the source frame. The time t ac corresponds to the time when the shock leaves the emission zone. The light curves are at frequencies 10 16 Hz (dashed line), 10 18 Hz (solid line) and 10 20 Hz (dot-dashed line) or respectively at energies 47 eV, 4.7 keV and 470 keV in the source frame.

Figure 4 . 2 :

 42 Figure 4.2: Delayed X-ray light curves (left) and Fourier transform of the delay (right) from time-dependent model of Mrk 421 1998 flare. The Fourier frequency is equivalent to the inverse time for the light curve. Positive lags correspond to low energies arriving after high energies. The break on the lag Fourier transform corresponds to the time at a given flux value, the low energy light curve reach this flux earlier than the high energy one. Figure taken from Lewis et al. (2016).

. 5 )

 5 Such equations keeps the dominant term to produce a flare in situ, corresponding to a minimal radiative scenario for an AGN γ-ray flare. For SSC emission, Equations 3.15 and 3.25 give the b syn and b IC coefficients, associated respectively with synchrotron and IC energy losses in the Thomson regime:

Figure 4 . 3 :

 43 Figure 4.3: Time evolution of the parameter η(t) from computation of the emission a single flare.The decrease of η(t) corresponds to the fast acceleration of electrons. When electrons reach a large energy, the radiative cooling becomes more important than the electron acceleration leading to an increasing η(t).

Figure 4 . 4 :

 44 Figure 4.4: Electron spectrum (left) and SED (right) evolution for the reference case. The color code is shared and indicated on the left panel. The spectra as the SED from 0 to 4440 s are shown during the electron acceleration phase while later times correspond to electron cooling phase.

Figure 4 . 5 :

 45 Figure 4.5: Left: normalized light curves for the reference case. The maximum flux value of each light curve is mentioned in the top-left legend. Right: Zoom around the light curve maxima.

Figure 4 . 6 :

 46 Figure 4.6: Simulated light curves with a non-constant width and time delay measurement with CCF and PPM. The CCF method appears does not reconstruct the injected delay.

Figure 4 . 7 :

 47 Figure 4.7: Simulated light curves with constant width and time delay measurement with CCF and PPM. The PPM and CCF are perfectly superimposed and reconstruct well the injected lag.

Figure 4 . 8 :

 48 Figure 4.8: Time delay with respect to the energy for the standard scenario. The time delay is defined with respect to the lowest energy light curve (1 -1.8 MeV) with no delay indicated with the horizontal dashed line.A positive time delay correspond to a light curve peaking later than the lowest energy one.

Figure 5 . 1 :Figure 5 . 2 :

 5152 Figure 5.1: Electron spectrum evolution for the reference case without adiabatic expansion. The spectra from 0 to 8450 s are shown during the electron acceleration phase while later time correspond to the electron cooling phase. The reference set of parameters is given in Table 4.1.

  around light curve maxima.

Figure 5 . 3 :

 53 Figure 5.3: Normalized light curves for the reference case without adiabatic expansion associated with the electron spectrum and SED evolution of Figure 5.1 and Figure 5.2.

Figure 5 . 4 :

 54 Figure 5.4: Time delay with respect to the energy for the reference case with and without adiabatic expansion. The time delay is defined with respect to the lowest energy light curve (1 -1.8 MeV) with no delay indicated with the horizontal dashed line. A positive time delay correspond to a light curve peaking later than the lowest energy one.

  Zoom around light curve maxima.

Figure 5 . 5 :

 55 Figure 5.5: Electron spectrum energy cut-off γ c (t) evolution (left) and associated light curves (right) for the reference case without adiabatic expansion. The time t max is defined when γ c (t) reaches its maximum value γ max and is indicated with the vertical dashed line for both figures. The energy ranges for each light curve are the same than Figure 5.3.

Figure 5 . 6 :

 56 Figure 5.6: Evolution of the energy dependent time delay for various initial magnetic field strength value B 0 . The delay is defined with respect to the lowest energy light curve. A constant time delay with respect to energy corresponds to an absence of delay between two adjacent energy light curves. The two B 0 values chosen to illustrate the two regimes are represented with thick crosses.

  Zoom around light curve maxima.

Figure 5 . 7 :

 57 Figure 5.7: Normalized light curves for the reference case with B 0 = 90 mG and without adiabatic losses. The time t max is indicated with a vertical dashed line. The delay associated with this cases is shown in Figure 5.6 with thick orange crosses.

Figure 5 . 8 :

 58 Figure 5.8: Evolution of the energy dependent time delay with various magnetic field evolution index m b for two initial magnetic field values.

Figure 5 . 9 :

 59 Figure 5.9: Evolution of the energy dependent time delay with various Doppler factor δ for two initial magnetic field values.

  B 0 = 90 mG.

Figure 5 .

 5 Figure 5.10: Evolution of the energy dependent time delay with various initial acceleration amplitude A 0 for two initial magnetic field values.

Figure 5 .

 5 Figure 5.11: Evolution of the energy dependent time delay with various acceleration evolution index m a for two initial magnetic field values.

Figure 5 .

 5 Figure 5.13: Evolution of the energy dependent time delay with various initial electron spectrum energy cut-off γ c,0 for two initial magnetic field values. Same description as Figure 5.6.

Figure 5 .

 5 Figure 5.14: Time delay with respect to the energy for two cases with adiabatic expansion with different magnetic field values. The time delay is the delay between the lowest energy light curve (1 -1.8 MeV) and the considered light curve. The dashed line correspond to no time delay (∆t = 0).

  (b) B 0 = 90 mG.

Figure 5 .

 5 Figure 5.15: Times delay in the GeV-TeV energy range for two initial magnetic field values. The dashed line correspond to ∆t = 0 while the full line corresponds to a power law function parameterization of the delay with respect to the energy.

  B 0 = 90 mG.

Figure 5 .

 5 Figure 5.16: Comparison between LIV and intrinsic time delay in the GeV-TeV energy range for two initial magnetic field values. The time delay, represented with black crosses, for B 0 = 65 mG (left) is represented with subluminal LIV delay while for B 0 = 90 mG (right), superluminal LIV effect is chosen. The dashed line correspond to ∆t = 0 while the full lines correspond to the expected LIV delays for a QG energy scale at the Planck energy for various redshifts. A redshift of z = 0.03 is assumed for the model.

  Model light curves. (b) Time delay temporal evolution.

Figure 5 .

 5 Figure 5.17: Times delay temporal evolution for the reference case with B 0 = 65 mG. Only positive delays are shown because of the y-axis scale, corresponding to low energy flare reaching a normalized flux after the high one.

Figure 5 .

 5 Figure 5.18: Times delay temporal evolution for the reference case with B 0 = 90 mG. Only positive delays are shown because of the y-axis scale, corresponding to low energy light curve reaching a normalized flux value after the high one.

]

Figure 5 .

 5 Figure 5.19: Times delay temporal evolution for simulated light curves delays by an average LIV effect.

Figure 6 . 1 :

 61 Figure 6.1: Time (left) and energy (right) distribution of the simple scenario. Both distributions are extracted from a simulated data set randomly selected. They are fitted with the functions used to generate the events. 99

Figure 6 . 2 : 1 .

 621 Figure 6.2: Log-likelihood function for one realization of the simple scenario with a injected parameter τ inj 1 = 0 s.TeV -1 . The log-likelihood function has a parabolic shape, with a single minimum corresponding the reconstructed value τ best 1. The red line corresponds to the threshold defined by the likelihood graphical method estimation of the 68% confidence level interval. The reconstructed value and the associated 1σ errors are reported on the top of the plot.

]

-

Figure 6 .

 6 Figure 6.3: LIV time lag distributions for a 1000 simulations of the simple example. The middle plot corresponds to the reconstructed lag τ best 1

]

Figure 6 . 4 :

 64 Figure 6.4: Calibration plots for the simple scenario with inputs defined in Table6.1. The full line is a parameterization of τ best n versus τ inj n with the fitted parameters indicated on bottom right. The dashed line represents the perfect reconstruction with a line of slope 1 and y-intercept of 0. The grey shaded area represents the 1σ confidence level on τ best n .

Figure 6 . 5 :

 65 Figure 6.5: Same as Figure 6.4 for the energy resolution scenario.

Figure 6 . 6 :

 66 Figure 6.6: Same as Figure 6.4 for the fitted template scenario.

  With lagged template (n=1).

  Corrected for lagged template (n=1).

Figure 6 . 7 :

 67 Figure 6.7: Same as Figure 6.4 for the lagged template (left) and corrected template (right) scenarios for a linear LIV effect (n=1).

Figure 7 . 1 :

 71 Figure 7.1: Illustration of the high zenith angle observation. The incoming γ ray with high ZA produce a larger and fainter Cherenkov light pool than the one with small ZA angle.

Figure 7 . 2 :

 72 Figure 7.2: θ 2 histogram of Mrk 501 flare observations. This represents the distribution of the observed events, depending on their distance to the center of the ON region. The θ distribution follows a Gaussian one centered in θ 2 = 0 and a standard deviation corresponding to the instrument PSF as Mrk 501 is a point-like source. At larger θ 2 , the flat distribution corresponds to the background events. In order reject most of the background events, a cut in θ 2 is done at 0.015 • .

Figure 7 . 3 :

 73 Figure 7.3: Significance map of Mrk 501 flare observations. The significance reflects the probability of the map being compatible with Gaussian noise only. The very high significance means that the data is not compatible with the hypothesis of being a Gaussian noise at a confidence level above 83σ.

Figure 7 . 4 :

 74 Figure 7.4: Energy spectrum of Mrk 501 flare observations, obtained with a curved power law spectral hypothesis. The top panel shows spectral reconstruction results with the black line, associated with the green area for the 68% confidence level contour (statistical errors only). The bottom panel shows the residual (difference between the observed number of event in the data and the expected number of event from the spectral shape hypothesis). The back data point are then derived from the residual, requiring a minimum significance of 3σ. For smaller significance, the values are taken as upper limits and shown as arrows. The detailed results of the fit are shown on the top right box.

Figure 7 . 5 :

 75 Figure 7.5: Light curve of Mrk 501 flare observations with a time binning of 420 seconds. The light curve gives the flux above the threshold energy of 1.3 TeV. The dashed line corresponds to the averaged flux on the whole observation at 74.5 × 10 -12 cm -2 .s -1 .

Figure 7 . 6 :

 76 Figure 7.6: Energy Spectrum of Mrk 501 above 3.25 TeV obtained with a power law spectral hypothesis. Same description as Figure 7.4.

Figure 7 . 7 :

 77 Figure 7.7: Template light curve of Mrk 501 parameterized with a simple and double Gaussian functions. The flux data points are drawn with triangle, associated with their error bars and red line is the best estimation of the light curve. The fitted parameters are indicated on the top right box. The double Gaussian function better represents of the data than the simple Gaussian one, as indicated by their respective χ 2 /ndf values.

Figure 7 . 8 :

 78 Figure 7.8: Template light curve parameterization results with different energy ranges. The function used to fit the different light curves is the sum of 2 Gaussian. The first raw of plots takes E cut = 3 TeV, the median of data energy distribution. The second one takes E cut = 3.25 TeV, the value used for the LIV data analysis. The second E cut values shows a more robust parameterization functions while the energy range changes.

Figure 7 .

 7 Figure 7.9: -2∆ log(L) functions obtained from Mrk 501 data. The horizontal line correspond to the condition -2∆ log(L) = 1 which shows the 68% two sided interval. The reconstructed values are written on the top of the plots. Both results are compatible with 0.

Figure 7 .

 7 10 shows the τ best n distribution. Both MC simulations show on average τ best n values compatible with the injected LIV parameter.

Figure 7 .

 7 Figure 7.10: Monte Carlo τ best n distributions with simulations of Mrk 501 data set. The middle plot corresponds to the distribution of reconstructed lag τ best n while the left and right plots correspond to the 68% confidence level lower (left) and upper (right) values, determined via the graphical method from the log-likelihood function. They are fitted with Gaussian functions and the median of the distributions is represented with a vertical dashed line.

Figure 7 .

 7 Figure 7.11: Calibration plots with Mrk 501 Monte Carlo simulations. The full line is a parameterization of the reconstructed τ best n versus τ inj n with fitted parameters indicated on bottom right. The dashed line represents the perfect reconstruction with a line of slope 1 and y-intercept of 0. The grey shaded area represents the 68% confidence level on τ best n .

  Quadratic LIV effect.

Figure 7 .

 7 Figure 7.12: Monte Carlo τ best n distributions with Mrk 501 Monte Carlo simulations for the parameterization systematic uncertainties study. Same description as Figure 7.10.

Figure 7 .

 7 Figure 7.13: -2∆ log(L) functions obtained from Mrk 501 data with SVL analysis configuration. Same description as Figure7.9.

Figure 7 .

 7 Figure 7.14: Instrument response functions for the Mrk 501 observations conditions with Com-bined3 analysis configuration.

Figure 7 .

 7 Figure 7.15: -2∆ log(L) functions obtained from Mrk 501 data for the energy bias systematic uncertainties study. Same description as Figure 7.9.

Figure 7 .

 7 Figure 7.16: Monte Carlo τ best n distributions with Mrk 501 Monte Carlo simulations for the background systematic uncertainties study. Same description as Figure 7.10.

  Quadratic subluminal E QG effect.

Figure 7 .

 7 Figure 7.17: E QG,n subluminal limits obtained with different blazar flares. The two limits on Mrk 501 were slightly shifted apart on redshift for better clarity. On the left plot, the dashed line corresponds to the Planck scale at the energy E P = 1.22 × 10 19 GeV. The red points correspond to the results of this work. The reference of each limits are mentioned inTable 7.4.

Table 1 .

 1

1: Analysis cut for the Combined3 and StereoVeryLoose analyses. These cuts are used for the analysis presented in Chapter 7.

Table 4 .

 4

	Value	Unit

1: Simulation input parameters of the reference case. This set of parameter is extensively used in Chapter 5.

Table 5 .

 5 2: Energy dependent time delay index

		0	175 ± 23	0.72 ± 0.1	-29 ± 8	0.57 ± 0.2
		1.25	300 ± 39	0.61 ± 0.1	154 ± 21	0.65 ± 0.1
		1.5	369 ± 59	0.49 ± 0.2	226 ± 29	0.63 ± 0.1
		2.0	391 ± 62	0.52 ± 0.2	365 ± 54	0.53 ± 0.2
	δ	20	413 ± 90	0.68 ± 0.5	-97 ± 15	0.63 ± 0.3
		30	227 ± 32	0.59 ± 0.2	-46 ± 10	0.62 ± 0.3
		40	175 ± 23	0.72 ± 0.1	-29 ± 8	0.57 ± 0.2
		50	96 ± 14	0.47 ± 0.1	-26 ± 10	0.49 ± 0.3
	A 0	4.0 × 10 -5	264 ± 36	0.61 ± 0.2	69 ± 12	0.60 ± 0.2
	(s -1 ) 4.5 × 10 -5	175 ± 23	0.72 ± 0.1	-29 ± 8	0.57 ± 0.2
		5.0 × 10 -5	81 ± 14	0.50 ± 0.1	-125 ± 18	0.63 ± 0.1
		5.5 × 10 -5	-	-	-202 ± 28	0.57 ± 0.1
		6.0 × 10 -5	-72 ± 11	0.58 ± 0.1	-225 ± 27	0.80 ± 0.1
	m a	4.7	51 ± 9	0.62 ± 0.1	-205 ± 28	0.59 ± 0.1
		5.0	94 ± 13	0.65 ± 0.1	-106 ± 15	0.66 ± 0.1
		5.3	158 ± 24	0.48 ± 0.2	-75 ± 12	0.67 ± 0.2
		5.6	175 ± 23	0.72 ± 0.1	-29 ± 8	0.57 ± 0.2
		5.9	190 ± 26	0.61 ± 0.2	-	-
	n	2.2	-	-	52 ± 11	0.43 ± 0.1
		2.3	173 ± 22	0.80 ± 0.1	-	-
		2.4	175 ± 23	0.72 ± 0.1	-29 ± 8	0.57 ± 0.2
		2.5	277 ± 153	0.26 ± 0.3	-72 ± 20	0.44 ± 0.2
		2.6	113 ± 26	0.66 ± 0.6	-104 ± 37	0.81 ± 0.7
		2.7	-	-	-	-
		2.8	-	-	-	-
	γ c,0	2 × 10 4	401 ± 51	0.64 ± 0.2	241 ± 31	0.66 ± 0.2
		3 × 10 4	265 ± 37	0.55 ± 0.2	57 ± 9	0.85 ± 0.2
		4 × 10 4	175 ± 23	0.72 ± 0.1	-29 ± 8	0.57 ± 0.2
		5 × 10 4	84 ± 13	0.57 ± 0.1	-82 ± 12	0.78 ± 0.2
		6 × 10 4	61 ± 14	0.42 ± 0.1	-131 ± 19	0.61 ± 0.1
		7 × 10 4	-	-	-151 ± 21	0.63 ± 0.1
		8 × 10 4	-	-	-205 ± 30	0.54 ± 0.1

Table 6 .

 6 2. 1: Simulation input parameters of the simple scenario.

	Inputs	Definition	Parameter values
	F sim (t)	Gaussian function	µ = 2000 s, σ = 750 s
	Γ sim (E)	Power law function	α = 2.2
	Acceptance	Constant acceptance assumed	None
	Energy Resolution Perfect energy resolution assumed None
	Template region	0.4 < E < 0.8 TeV	500 events
	Likelihood region	0.8 < E < 10 TeV	500 events

-9.4 s.TeV Median

	Number of sets 250			-1			/ ndf Prob 2 χ A LL τ LL tau σ	-	11.16 / 8 0.1926 9.6 ± 247.2 0.347 ± 9.683 0.24 ± 10.64	Number of sets 250	1.0 s.TeV Median	-1			/ ndf Prob 2 χ A best τ τ σ / ndf 2 χ Prob A best τ τ σ	8.355 / 9 0.4988 9.5 ± 247.7 0.339 ± 1.008 0.23 ± 10.65 8.355 / 9 0.4988 9.5 ± 247.7 0.339 ± 1.008 0.23 ± 10.65	Number of sets 250						/ ndf Prob 2 χ A UL τ UL tau σ	10.44 / 8 0.2355 9.5 244.7 ± 0.34 ± 11.19 0.24 ± 10.76
	200												200												200
	150												150												150
	100												100												100
	50												50												50
	100 -0	80 -	60 -	40 -	20 -	0	20	40		60	80	100	100 -0	80 -	60 -	40 -	20 -	0	20	40	60	80	100	100 -0	80 -	60 -	40 -	20 -	0	20	40	60	80 100
									LL τ	[s.TeV	1										best τ	[s.TeV	] -1									UL τ	[s.TeV	] -1

-1 11.3 s.TeV Median

Table 6 .

 6 .2. 2: Scenario list for the MC simulation studies.

	Scenario name	Description	Subsection
	Simple scenario	Reference scenario	6.3.3
	Energy resolution scenario Energy resolution included in event gen-	6.4.1
		erations
	Fitted template scenario	Energy resolution and non lagged tem-	6.4.2
		plate parameterization
	Lagged template scenario	Energy resolution and lagged template pa-	6.4.3
		rameterization
	Corrected template scenario Energy resolution, lagged template pa-	6.4.4
		rameterization and correction in the PDF

Table 7 .

 7 Double Gaussian function A 1 = 80.58 cm -2 .s -1 , µ 1 = 2361 s, σ 1 = 2153 s A 2 = 60.51 cm -2 .s -1 , µ 2 = 6564 s, σ 2 = 676 s 2: Analysis parameters for the LIV analysis of Mrk 501.

	Inputs	Definition	Parameter values
	F (t)
	Γ(E)	Power law function	α = 3.1
	Acceptance	Included in F (t) and Γ(E) None
	Energy Resolution Fixed energy resolution	∆E E = 15%
	Template region	1.3 < E < 3.25 TeV	773 events, E T = 2.14TeV
	Likelihood region	3.25 < E < 20 TeV	662 events

-2 -1.3 s.TeV Median

  6 -1.4 s TeV -2 .

	Number of sets 180 200	-18.4 s.TeV Median	-1		/ ndf Prob 2 χ A LL τ LL tau σ	21.85 / 13 0.05771 8.4 ± 195.6 0.64 ± 18.87 -0.57 ± 19.96	Number of sets 180 200	1.8 s.TeV Median	-1			/ ndf Prob 2 χ A best τ τ σ / ndf 2 χ Prob A best τ τ σ	25.88 / 15 0.03934 7.8 ± 181 0.701 ± 1.346 0.62 ± 21.48 25.88 / 15 0.03934 7.8 ± 181 0.701 ± 1.346 0.62 ± 21.48	Number of sets 180 200	24.1 s.TeV -1 Median	/ ndf Prob 2 χ A UL τ UL tau σ	27.51 / 15 0.02486 7.1 ± 0.77 ± 0.63 ± 23.01 23.88 168.8
	160												160									160
	140												140									140
	120												120									120
	100												100									100
	80												80									80
	60												60									60
	40												40									40
	20												20									20
	150 -0	100 -	50 -	0	50	100		150	150 -0		100 -		50 -	0		50	100	150	150 -0	100 -	50 -	0	50	100	150
							LL τ	[s.TeV	] -1								best τ	[s.TeV	] -1	UL τ	[s.TeV	] -1
															(a) Linear LIV effect.
	Number of sets 250 300					/ ndf Prob 2 χ A LL τ LL tau σ	-	31.96 / 9 0.0002024 14.0 ± 315.6 0.053 ± 1.365 0.050 ± 1.632	300 Number of sets 250							/ ndf Prob 2 χ A best τ τ σ / ndf 2 χ Prob A best τ τ σ	24.79 / 10 0.00575 12.6 ± 295.8 0.0595 ± 0.1233 0.050 ± 1.754 24.79 / 10 0.00575 12.6 ± 295.8 0.0595 ± 0.1233 0.050 ± 1.754
	200												200
	150												150
	100												100
	50												50
	20 -0	15 -	10 -	5 -	0	5	10			15		20	20 -0	15 -	10 -	5 -	0	5	10	15	20
							LL τ	[s.TeV	] -2									best τ	[s.TeV	] -2

-2 0.1 s.TeV Median

	Number of sets 250 300	1.7 s.TeV Median	-2			/ ndf Prob 2 χ A UL τ UL tau σ	29.04 / 9 0.0006383 11.8 ± 0.06 0.056 ± 1.902 ± 1.84 271.8
	200
	150
	100
	50
	20 -0	15 -	10 -	5 -	0	5	10	15	20
								UL τ	[s.TeV	] -2
	(b) Quadratic LIV effect.

  6 -3.5 s TeV -1 , ∆τ 2,fit < +0.4 -0.2 s TeV -2 .

	Number of sets 160 180 200	-14.9 s.TeV Median	-1			/ ndf Prob 2 χ A LL τ LL tau σ	45.66 / 18 0.000333 7.0 ± 159.9 1.35 ± 13.78 -1.18 ± 39.69	Number of sets 160 180 200	6.9 s.TeV Median	-1			/ ndf Prob 2 χ A best τ τ σ / ndf 2 χ Prob A best τ τ σ	60.06 / 19 06 -3.789e 6.8 ± 157.4 1.401 ± 6.627 1.14 ± 39.72 60.06 / 19 06 -3.789e 6.8 ± 157.4 1.401 ± 6.627 1.14 ± 39.72	200 Number of sets 160 180	28.7 s.TeV -1 Median	/ ndf Prob 2 χ A UL τ UL tau σ	53.28 / 18 05 -2.371e 6.3 ± 1.54 ± 1.23 ± 42.96 30.05 146.4
	140															140								140
	120															120								120
	100															100								100
	80															80								80
	60															60								60
	40															40								40
	20															20								20
	250 -0	200 -	150 -	100 -	50 -	0	50	100 150 200 250	250 -0	200 -	150 -	100 -	50 -	0	50 100 150 200 250	250 -0	200 -	150 -	100 -	50 -	0	50	100 150 200 250
										LL τ	[s.TeV	] -1									best τ	[s.TeV	] -1	UL τ	[s.TeV	] -1
																	(a) Linear LIV effect.
	Number of sets 180 200 220	-2 -1.1 s.TeV Median			/ ndf Prob 2 χ A LL τ LL tau σ	-	85.58 / 20 10 -4.339e 8.9 ± 176.9 0.095 ± 1.028 0.104 ± 2.749	180 220 Number of sets 200	0.3 s.TeV Median	-2			/ ndf Prob 2 χ A best τ τ σ / ndf 2 χ Prob τ σ A best τ	76.94 / 18 09 -2.924e 8.3 ± 173.4 0.1041 ± 0.4184 0.098 ± 2.833 76.94 / 18 09 -2.924e 0.098 ± 2.833 8.3 ± 0.1041 ± 0.4184 173.4
	160															160
	140															140
	120															120
	100															100
	80															80
	60															60
	40															40
	20															20
	20 -0	15 -	10 -		5 -	0	5		10			15		20	20 -0	15 -	10 -	5 -	0	5	10	15	20
										LL τ	[s.TeV	] -2									best τ	[s.TeV	] -2	UL	[s.TeV	] -2

  1 -0.1 s TeV -2 .

	Number of sets 160 180 200	-18.4 s.TeV Median	-1		/ ndf Prob 2 χ A LL τ LL tau σ	13.73 / 12 0.318 8.1 ± 201.9 0.63 ± 19.68 -0.5 ± 19.5	Number of sets 160 180 200	1.5 s.TeV Median	-1			/ ndf Prob 2 χ A best τ τ σ / ndf 2 χ Prob A best τ τ σ	15.07 / 13 0.3033 7.6 ± 185.1 0.684 ± 1.376 0.56 ± 21.27 15.07 / 13 0.3033 7.6 ± 185.1 0.684 ± 1.376 0.56 ± 21.27	200 Number of sets 160 180	23.3 s.TeV -1 Median		/ ndf Prob 2 χ A UL τ UL tau σ	23.56 / 14 0.05179 7.1 ± 0.75 ± 0.63 ± 22.91 23.42 170.5
	140											140											140
	120											120											120
	100											100											100
	80											80											80
	60											60											60
	40											40											40
	20											20											20
	150 -0	100 -	50 -	0	50		100		150	150 -0	100 -		50 -	0	50	100		150	150 -0	100 -		50 -	0	50	100	150
							LL τ	[s.TeV	] -1								best τ	[s.TeV	] -1						UL τ	[s.TeV	] -1
													(a) Linear LIV effect
	Number of sets 300	-2 -1.4 s.TeV Median		/ ndf Prob 2 χ A LL τ LL tau σ	-	26.45 / 9 0.001726 14.2 ± 327.8 0.051 ± 1.419 0.05 ± 1.58	Number of sets 300	0.0 s.TeV Median	-2			/ ndf Prob 2 χ A best τ τ σ / ndf 2 χ Prob A best τ τ σ	25.47 / 9 0.002496 12.6 ± 298.4 0.05580 ± 0.02535 0.048 ± 1.737 25.47 / 9 0.002496 12.6 ± 298.4 0.05580 ± 0.02535 0.048 ± 1.737	Number of sets 300	1.8 s.TeV Median	-2		/ ndf Prob 2 χ A UL τ UL tau σ	36.87 / 9 05 -2.776e 12.1 ± 0.063 0.056 ± 1.865 ± 1.728 275
	250											250											250
	200											200											200
	150											150											150
	100											100											100
	50											50											50
	20 -0	15 -	10 -	5 -	0	5	10		15		20	20 -0	15 -	10 -	5 -	0	5	10	15		20	20 -0	15 -	10 -	5 -	0	5	10	15	20
							LL τ	[s.TeV	] -2								best τ	[s.TeV	] -2						UL τ	[s.TeV	] -2
													(b) Quadratic LIV effect

  2,i (s TeV -2 ) Table7.3: Systematic uncertainties for each contributions on the τ n measurement. The total systematic is a quadratic sum of the individual ones.

	Likelihood method calibration	+5.5 -2.8	+0.4 -0.5
	F (t) and Γ(E)	+4.6 -3.5	+0.4 -0.2
	Photon list	+10.6 -8.2	+0.7 -0.6
	Energy calibration (bias)	+2.3 -5.2	+0.1 -0.6
	Background	+0.8 -0.1	+0.1 -0.1
	Total	+13 -11	+0.9 -1.0

  63 × 10 17 GeV, subluminal case 2.89 × 10 17 GeV, superluminal case 10 10 GeV, subluminal case 7.66 × 10 10 GeV, superluminal case

	E 95% QG,2 =	  8.79 ×
		

Table 7 .

 7

	Source	Redshift	Number of	Energy range Variability	Reference
	name		events in ML	(TeV)	(min)
	PKS 2155-304	0.113	3526	0.28 -4	1	H.E.S.S. Collaboration (2011)
	PG 1553+113	∼ 0.5	154	0.3 -0.8	10	Abramowski et al. (2015)
	Mrk 501	0.034	1500	0.15 -10	4	Martínez & Errando (2009)
		0.034	1430	1.3 -20	7	This thesis

4: Active Galactic Nuclei characteristics.

The mathematical relation a log b = b log a is used.

-16 -14 -12 -10 -8 -6 -4 -2 -

-14 -13 -12 -11 -10 -9 -8

-10 5 -10 4 -10 3 -10 2 -

-15 -14 -13 -12 -

-200 -100 -

-350 -300 -250 -200 -150 -100 -50

-10 3 -10 2 -10

-10 10 -10 8 -10 6 -10 4 -10

(1 + x) α ≈ 1 + αx when x → 0

Violation de l'Invariance de Lorentz (LIV en anglais) à l'échelle de Planck. Une des signatures de cette violation est la modification de la propagation des photons dans le vide, induisant des décalages temporels dépendant de l'énergie des photons observés sur Terre. De tels décalages dans le temps d'arrivée de rayons γ sont recherchés avec l'expérience H.E.S.S. (High Energy Stereoscopic System), grâce aux émissions de très hautes énergies en provenance de sources lointaines telles que les blazars. Néanmoins, l'origine du décalage temporel doit être comprise en détails. En effet, un décalage intrinsèque à la source pourrait venir biaiser les contraintes sur les modèles de gravitation quantique.Cette thèse propose dans un premier temps de s'intéresser à la modélisation temporelle des éruptions de blazars, pour étudier les possibles décalages intrinsèques liés aux processus d'émissions de ces éruptions. Grâce à l'élaboration d'un modèle simple, cette étude met en relief les différentes caractéristiques de ces décalages intrinsèques sur les scénarios d'éruptions de blazar afin d'essayer de les distinguer des décalages potentiellement dus à un effet de LIV et aussi de proposer de nouvelles contraintes basées sur ces décalages temporels. Dans un deuxième temps, la méthode de recherche de décalages temporels dépendant de l'énergie avec H.E.S.S. est présentée ainsi qu'une application sur l'éruption du blazar Markarian 501 ayant eu lieu en Juillet 2014. Cette analyse a permis d'établir la meilleure limite obtenue sur le terme quadratique de la signature de la LIV avec l'utilisation d'éruption de blazars.Mots-clés : éruption de blazar, décalage temporel, Violation de l'invariance de Lorentz, modélisation dépendent du temps, rayons gamma, H.E.S.S., astrophysique de haute énergie

Remerciements

Acknowledgments laboratoire avec qui j'ai pu discuter et partager une biennale à Pise. Merci à Ahmed, Andreas,

List of abbreviations

Chapter 6

Likelihood method for LIV search Contents the time delay, thus reducing the number of events for the ML. The assumption of neglecting LIV effect can be considered as being valid only at the lowest energy, e.g. typically below 1 TeV.

If LIV delay cannot be neglected, this can lead to a wrong estimation of the LIV parameter. A new correction to deal with those cases will be introduced in more details in Section 6.4.3.

Treatment of the acceptance term

The acceptance factor can be treated with two different approaches. During short observations, acceptance does not significantly change in time and thus time dependency can be neglected. Besides, as the acceptance above the threshold energy does not significantly vary, it can be considered constant and all the A(E i ) terms can be treated as a constant factor and can be removed from the integral since they do not depend on the parameter of interest.

On the contrary, for long observations, acceptance variations have to be taken into account. They can lead to energy dependent time delay which will bias the ML estimation of the LIV parameter. An efficient way to include acceptance in the PDF is to use the light curve and energy spectrum instead of the time and energy distributions. Light curves and spectra are already deconvoluted from acceptance variation in contrary to the time and energy distributions. In that case, the acceptance term should be removed from the PDF.

Treatment of energy resolution

The energy resolution involves in the PDF formula an integration to take into account all possible true energy values with respect to the available measured energy. It requires a lot of calculations when evaluating the log-likelihood function, increasing the computation time to evaluate the log-likelihood function. To reduce the calculation time, the energy resolution integration is removed from the PDF. Actually, the performances of the ML without energy resolution are not significantly degraded as it will be shown in Section 6.4.1.

But, this assumption is valid only if Γ(E) follows a power law function. Indeed, the energy resolution affects the true energy spectrum Γ(E ) and gives a measured spectrum by a convolution product between Γ(E ) and the energy resolution R(E, E ). In case of a true energy spectrum following power law functions, the convolution product gives a measured spectrum following a power law function not significantly different from the true one. Whereas, in case of more complex true spectra like curved power law or broken power law functions, the convolution product result in a measured spectrum significantly different from the initial one.

Simplified probability density function

After applying all these simplification and approximation, the PDF in Equation 6.5 become: .14) if Γ(E) is a power law function.

The toy Monte-Carlo simulation tool

In order to evaluate the performances of the method, a MC tool has been built and used to generate data sets with a lag. These simulated sets can then be analyzed via the ML described in the previous sections. This tool is used to investigate the behavior of the ML and also its

Conclusion

In this section, the Maximum Likelihood method was presented in order to measure energy dependent time delay coming from a LIV effect. This method relies on the definition of the Probability Density Function describing the data set considered, and the log-likelihood function. Different simplifications and approximations were discussed and even tested with the use a toy Monte-Carlo. The procedure used for this tool was described and applied to a simple case considering a Gaussian shape light curve and an energy spectrum following a power law function. Then multiple scenarios were investigated to check the validity of the different hypothesis that can be made in the Maximum Likelihood analysis. The two last scenarios showed that in case of very high energy data sets, the Lorentz Invariance Violation delay cannot be neglected anymore for the determination of F (t) which leads to a wrong reconstruction of the parameter of interest.

To solve this issue, a correction was designed and introduced in the Probability Density Function. This correction was tried with the toy Monte-Carlo tool and showed a recovery of the perfect reconstruction. As a general result, this correction should be always considered in the case of high energy threshold for the considered data set which correspond to case with a minimum energy around TeV energies. Actually, this situation correspond to the Markarian 501 2014 flare, observed by the H.E.S.S. experiment which is used in next chapter to search for Lorentz Invariance Violation signatures.

Chapter 7

Searching for Lorentz invariance violation with Markarian 501 Contents Abstract: Specific models of quantum gravity suggest the existence of a Lorentz Invariance Violation (LIV) at the Planck scale. One signature of that violation is a modification the propagation of photons in vacuum which induces energy dependent delays in the arrival time of photons on Earth. The H.E.S.S. (High Energy Stereoscopic System) experiment can search for such delays in the arrival time of gamma rays, thanks to the very high energy emission coming from distant blazars. However, the time delay origin have to be fully understood. Indeed, an intrinsic time delay coming from the source can bias the constraints made on quantum gravity models.

In the first part of this thesis, a time dependent blazar flare model is considered to search for the presence of intrinsic time delays related to the emission mechanisms of flares. With the elaboration of a simple scenario, this study highlights the different characteristics of intrinsic time delays in order to investigate how to disentangle them from delays due to LIV as well as to provide new constraints on blazar modeling. In the second part of this thesis, the method used to search for LIV signatures in blazar light curves at very high energy is presented as well as an application to the flare of Markarian 501 which occurred in July 2014. This analysis provides in particular the best upper limit on the quadratic term of LIV signature.

Key words: blazar flare, time delay, Lorentz invariance violation, time dependent modeling, gamma rays, H.E.S.S., high energy astrophysics Résumé : Des modèles spécifiques de gravitation quantique suggèrent l'existence d'une