THE 2012 FLARE OF PG 1553+113 SEEN WITH H.E.S.S. AND FERMI-LAT

H.E.S.S. is an array of five imaging atmospheric Cherenkov telescopes located in the Khomas highland in Namibia ($23{}^\circ 16^{\prime} 18^{\prime\prime}$ S, 16°30'01'' E), at an altitude of 1800 m above sea level (Hinton & the H.E.S.S. Collaboration 2004). The fifth H.E.S.S. telescope was added to the system in 2012 July and is not used in this work, reporting only on observations prior to that time.

PG 1553+113 was observed with H.E.S.S. in 2005 and 2006 (Aharonian et al. 2008). No variability was found in these observations, which will be referred to as the "pre-flare" data set in the following. New observations were carried out in 2012 April after flaring activity at VHE was reported by the MAGIC collaboration ("flare" data set, Cortina 2012b).

The pre-flare data set is composed of 26.4 live time hours of good-quality data (Aharonian et al. 2006b). For the flare period, eight runs of ∼28 minutes each were taken during the nights of 2012 April 26 and 27, corresponding to 3.5 hr of live time. All the data were taken in wobble mode, for which the source is observed with an offset of $0\buildrel{\circ}\over{.} 5$ with respect to the center of the instrument's field of view, yielding an acceptance-corrected live time of 24.7 and 3.2 hr for the pre-flare and flare data sets, respectively.

Data were analyzed using the Model analysis (de Naurois & Rolland 2009) with Loose cuts. This method–based on the comparison of detected shower images with a pre-calculated model–achieves a better rejection of hadronic air showers and a better sensitivity at lower energies than analysis methods based on Hillas parameters. The chosen cuts, best suited for sources with steep spectra such as PG 1553+113,⁴⁵ require a minimum image charge of 40 photoelectrons, which provides an energy threshold of $\sim 217$ GeV for the pre-flare and $\sim 240$ GeV for the flare data set.⁴⁶ All the results presented in this paper were cross-checked with the independent analysis chain described in Becherini et al. (2011).

Events in a circular region (ON region) centered on the radio position of the source, ${{\alpha }_{{\rm J}2000}}={{15}^{{\rm h}}}{{55}^{{\rm m}}}43\buildrel{\rm{s}}\over{.}04,$ ${{\delta }_{{\rm J}2000}}=11{}^\circ 11^{\prime} 24\buildrel{\prime\prime}\over{.} 4$ (Green et al. 1986), with a maximum squared angular distance of 0.0125 deg², are used for the analysis. In order to estimate the background in this region, the reflected background method (Berge et al. 2007) is used to define the OFF regions. The excess of γ-rays in the ON region is statistically highly significant (Li & Ma 1983): 21.5σ for the pre-flare period and 22.0σ for the flare. Statistics are summarized in Table 1.

Table 1.  Summary of the Statistics for Both Data Sets (First Column)

Data Set	ON	OFF	r	Excess	Significance	${{E}_{{\rm th}}}({\rm GeV})$	Zenith Angle	$P_{{{\chi }^{2}}}^{{\rm cst}}$
Pre-flare	2205	13033	0.100	901.7	21.5	217	34°	0.77
Flare	559	1593	0.105	391.2	22.0	240	52°	3.3 × 10−3

Note. The second and third columns give the number of ON and OFF events. The fourth column gives the ratio between ON and OFF exposures (r). The excess and the corresponding significance are given, as well as the energy threshold and the mean zenith angle of the source during the observations. The last column presents the probability of the flux to be constant within the observations (see text).

Download table as:  ASCIITypeset image

The differential energy spectrum of the VHE γ-ray emission has been derived using a forward-folding method (Piron et al. 2001). For the observations prior to 2012 April, a power law (PWL) model fitted to the data gives a χ² of 51.7 for 40 degrees of freedom (dof, corresponding to a χ² probability of ${{P}_{{{\chi }^{2}}}}=0.10$). The values of the spectral parameters (see Table 2) are compatible with previous analyses by H.E.S.S. covering the same period (Aharonian et al. 2008). A log-parabola (LP) model,⁴⁷ with a χ² of 37.5 for 39 dof (${{P}_{{{\chi }^{2}}}}=0.54$), is found to be preferred over the PWL model at a level of 4.3σ using the log-likelihood ratio test. Note that systematic uncertainties, presented in Table 2, have been evaluated by Aharonian et al. (2006b) for the PWL model and using the jack-knife method for the LP model. The jack-knife method consists in removing one run and redoing the analysis. This process is repeated for all runs.

Table 2.  Summary of the Fitted Spectral Parameters for the Pre-flare and the Flare Data Sets and the Corresponding Integral Flux I Calculated Above 300 GeV

Data Set (Model)	Spectral Parameters	I (E > 300 GeV)	Edec
		(10−12 ph cm−2 s−1)	(GeV)
Pre-flare (PWL)	Γ = 4.8 ± 0.2stat ± 0.2sys	4.4 ± 0.4stat ± 0.9sys	306
Pre-flare (LP)	$a\;=\;5.4\pm {{0.4}_{{\rm stat}}}\pm {{0.1}_{{\rm sys}}}$	5.0 ± 0.6stat ± 1.0sys	⋯
	b = 4.0 ± 1.4stat ± 0.2sys
Flare (PWL)	Γ = 4.9 ± 0.3stat ± 0.2sys	15.1 ± 1.3stat ± 3.0sys	327

Note. The Last Column gives the Decorrelation Energy.

Download table as:  ASCIITypeset image

For the flare data set, the LP model does not significantly improve the fit and the simple PWL model describes the data well, with a χ² of 33.0 for 23 dof (${{P}_{{{\chi }^{2}}}}=0.08$). Table 2 contains the integral fluxes above the reference energy of 300 GeV. The flux increased by a factor of ∼3 in the flare data set compared to the pre-flare one with no sign of spectral variations (when comparing power law fits for both data sets). The derived spectra and error contours for each data set are presented in Figure 1, where the spectral points obtained from the cross-check analysis are also plotted.

Zoom In Zoom Out Reset image size

To compute the light curves, the integrated flux above 300 GeV for each observation run was extracted using the corresponding (pre-flare or flare) best fit spectral model. A fit with a constant of the run-wise light curve of the entire (pre-flare+flare) data set, weighted by the statistical errors, yields a χ² of 123.2 with 68 dof (${{P}_{{{\chi }^{2}}}}=6.6\times {{10}^{-5}}$). Restricting the analysis to the pre-flare data set only, the fit yields a χ² of 51.76 with 60 dof (${{P}_{{{\chi }^{2}}}}=0.77$), indicating again a flux increase detected by H.E.S.S. at the time of the flaring activity reported by Cortina (2012b).

Figure 2 shows the light curve during the flare together with the averaged integral fluxes above 300 GeV of both data sets. A fit with a constant to the H.E.S.S. light curve during the first night yields a χ² of 20.76 for 6 dof (${{P}_{{{\chi }^{2}}}}=2.0\times {{10}^{-3}}$), indicating intra-night variability. This is also supported by the use of a Bayesian block algorithm (Scargle 1998) that finds three blocks for the two nights at a 95% confidence level.

Zoom In Zoom Out Reset image size

The Fermi-LAT is detector converting γ-rays to e⁺e⁻ pairs (Atwood et al. 2009). The LAT is sensitive to γ-rays from 20 MeV to >300 GeV. In survey mode, in which the bulk of the observations are performed, each source is seen every 3 hr for approximately 30 minutes.

The Fermi-LAT data and software are available from the Fermi Science Support Center.⁴⁸ In this work, the ScienceTools V9R32P5 were used with the Pass 7 reprocessed data (Bregeon et al. 2013), specifically SOURCE class event (Ackermann et al. 2012a), with the associated P7REP_SOURCE_V15 instrument response functions. Events with energies from 300 MeV to 300 GeV were selected. Additional cuts on the zenith angle (<100°) and rocking angle (<52°) were applied as recommended by the LAT collaboration⁴⁹ to reduce the contamination from the Earth atmospheric secondary radiation.

The analysis of the LAT data was performed using the Enrico Python package (Sanchez & Deil 2013). The sky model was defined as a region of interest of 15° radius with PG 1553+113 in the center and additional point-like sources from the internal 4 yr source list. Only the sources within a 3° radius around PG 1553+113 and bright sources (integral flux greater than 5 × 10⁻⁷ ph cm⁻² s⁻¹) had their parameters free to vary during the likelihood minimization. The template files isotrop_4years_P7_V15_repro_v2_source.txt for the isotropic diffuse component, and template_4years_P7_v15_repro_v2.fits for the standard Galactic model, were included. A binned likelihood analysis (Mattox et al. 1996), implemented in the gtlike tool, was used to find the best-fit parameters.

As for the H.E.S.S. data analysis, two spectral models were used: a simple PWL and a LP. A likelihood ratio test was used to decide which model best describes the data. Table 3 gives the results for the two time periods considered in this work, and Figure 3 presents the γ-ray spectral energy distributions. The first one (pre-flare), before the H.E.S.S. exposures in 2012, includes more than 3.5 yr of data (from 2008 August 4 to 2012 March 1). The best fit model is found to be the LP (with a Test Statistic⁵⁰ of 11.3, ≈3.4σ). The second period (flare) is centered on the H.E.S.S. observation windows and lasts for 7 days. The best fit model is a power law, the flux being consistent with the one measured during the first 3.5 yr. Data points or light curves were computed within a restricted energy range or time range using a PWL model with the spectral index frozen to 1.70.

Zoom In Zoom Out Reset image size

To precisely probe the variability in HE γ-rays, 7-day time bins were used to compute the light curve of PG 1553+113 in an extended time window (from 2008 August 4 to 2012 October 30), to probe any possible delay of a HE flare with respect to the VHE one. While the flux of PG 1553+113 above 300 MeV is found to be variable in the whole period with a variability index of ${{F}_{\operatorname{var}}}=0.16\pm 0.04$ (Vaughan et al. 2003), there is no sign of any flaring activity around the 2012 H.E.S.S. observations. This result has been confirmed by using the Bayesian block algorithm, which finds no block around the H.E.S.S. exposures in 2012. Similar results were obtained when considering only photons with an energy greater than 1 GeV. No sign of enhancement of the HE flux associated to the VHE event reported here was found. This might be due to the lack of statistics at high energy in the LAT energy range.

The VHE data do not show any sign of variation of the spectral index (when comparing flare and pre-flare data sets with the same spectral model), and in HE no counterpart of this event can be found. The indication for intra-night variability is similar to other TeV HBLs (Mrk 421, Mrk 501 or PKS 2155-304) with, in this case, flux variations of a factor 3.

As noticed in previous works, PG 1553+113 presents a sharp break between the HE and VHE ranges (Abdo et al. 2010a) and the peak position of the γ-ray spectrum in the νf(ν) representation is located around 100 GeV. This is confirmed by the fact that the LP model better represents the pre-flare period in HE. Nonetheless, the precise location of this peak cannot be determined with the Fermi-LAT data only. Combining both energy ranges and fitting the HE and VHE data points with a power law with an exponential cutoff⁵¹ allows us to determine the νf(ν) peak position for both time periods. The functional form of the model is

$$\begin{eqnarray*}&&{{E}^{2}}\frac{dN}{dE}=N{{\left( \frac{E}{100\;{\rm GeV}} \right)}^{-{\Gamma }}}{\rm exp} (-E/{{E}_{{\rm c}}}).\end{eqnarray*}$$

For this purpose, Fermi-LAT and H.E.S.S. systematic uncertainties were taken into account in a similar way as in Abramowski et al. (2014) and added quadratically to the statistical errors. The Fermi-LAT systematic uncertainties were estimated by Ackermann et al. (2012a) to be 10% of the effective area at 100 MeV, 5% at 316 MeV and 15% at 1 TeV and above. For the VHE γ-ray range, they were taken into account by shifting the energy by 10%. This effect translates into a systematic uncertainty for a single point of $\sigma {{(f)}_{{\rm sys}}}=0.1\cdot \partial f/\partial E$ where f is the differential flux at energy E.

The results of this parameterization are given in Table 4. Using the pre-flare period, the peak position is found to be located at ${{{\rm log} }_{10}}({{E}_{{\rm max} }}/1\;{\rm GeV})=1.7\pm {{0.2}_{{\rm stat}}}\pm {{0.4}_{{\rm sys}}}$ with no evidence of variation during the flare and no spectral variation. This is consistent with the fact that no variability in HE γ-rays was found during the H.E.S.S. observations. This is also in agreement with the fact that HBLs are less variable in HE γ-rays than other BL Lac objects (Abdo et al. 2010b), while numerous flares have been reported in the TeV band.

Table 4.  Parametrization Results of the Two Time Periods (First Column) Obtained by Combining H.E.S.S. and Fermi-LAT

Period	N (E = 100 GeV)	Γ	${{{\rm log} }_{10}}({{E}_{{\rm c}}}/1\;{\rm GeV})$	${{{\rm log} }_{10}}({{E}_{{\rm max} }}/1\;{\rm GeV})$
	(10−11 ${\rm erg}\;{\rm c}{{{\rm m}}^{-2}}\;{{{\rm s}}^{-1}}$)
Pre-flare	9.6 ± 0.7stat ± 1.7sys	1.59 ± 0.02stat ± 0.03sys	2.03 ± 0.02stat ± 0.04sys	1.7 ± 0.2stat ± 0.4sys
Flare	13.0 ± 3.5stat ± 5.7sys	1.56 ± 0.08stat ± 0.11sys	2.16 ± 0.04stat ± 0.09sys	1.8 ± 0.7stat ± 1.3sys

Note. The second column gives the normalization at 100 GeV, while the third and the fourth present the spectral index and cut-off energy of the fitted power law with an exponential cut-off. The last column is the peak energy in a $\nu f(\nu )$ representation.

Download table as:  ASCIITypeset image

The EBL is a field of UV to far-infrared photons produced by the thermal emission from stars and reprocessed starlight by dust in galaxies (see Hauser & Dwek 2001, for a review) that interacts with very high energy γ-rays from sources at cosmological distances. As a consequence, a source at redshift z exhibits an observed spectrum ${{\phi }_{{\rm obs}}}(E)={{\phi }_{\operatorname{int}}}(E)\times {{e}^{-\tau (E,z)}}$ where ${{\phi }_{\operatorname{int}}}(E)$ is the intrinsic source spectrum and τ is the optical depth due to interaction with the EBL. Since the optical depth increases with increasing γ-ray energy, the integral flux is lowered and the spectral index is increased.⁵² In the following, the model of Franceschini et al. (2008) was used to compute the optical depth τ as a function of redshift and energy. In this section, the data taken by both instruments during the flare period are used, with the Fermi-LAT analysis restricted to the range 300 MeV < E < 80 GeV (see Table 3 for the results). In the modest redshift range of VHE emitters detected so far ($z\leqslant 0.6$), the EBL absorption is negligible below 80 GeV (${{\tau }_{\gamma \gamma }}\sim 0.1$ at 80 GeV for z = 0.6).

A measure of the EBL energy density was obtained by Ackermann et al. (2012b) and Abramowski et al. (2013b) based on the spectra of sources with a known z. In the case of PG 1553+113, for which the redshift is unknown, the effects of the EBL on the VHE spectrum might be used to derive constraints on its distance. Ideally, this would be done by comparing the observed spectrum with the intrinsic one but the latter is unknown. The Fermi-LAT spectrum, derived below 80 GeV, can be considered as a proxy for the intrinsic spectrum in the VHE regime, or at least, as a solid UL (assuming no hardening of the spectrum).

Table 5.  Calibrated 95% 1-Sided LL and UL (including systematic errors) on the Dispersion Parameter ${{\tau }_{n}}$ and Derived 95% one-sided Lower Limits on E_QG

	Limits on τn (s TeV−n)		Lower Limits on EQG (GeV)
n	${\rm L}{{{\rm L}}^{{\rm calib}+{\rm syst}}}$	ULcalib + syst	s = −1	s = +1
1	−838.9	576.4	2.83 × 1017	4.11 × 1017
2	−1570.5	1012.4	1.68 × 1010	2.10 × 1010

Download table as:  ASCIITypeset image

Following the method used by Abramowski et al. (2013a), it has been assumed that the intrinsic spectrum of the source in the H.E.S.S. energy range cannot be harder than the extrapolation of the Fermi-LAT measurement. From this, one can conclude that the optical depth cannot be greater than τ_max(E), given by:

$$\begin{eqnarray}&&{{\tau }_{{\rm max} }}(E)={\rm ln} \left[ \frac{{{\phi }_{\operatorname{int}}}}{(1-\alpha )({{\phi }_{{\rm obs}}}-1.64{\Delta }{{\phi }_{{\rm obs}}})} \right],\end{eqnarray} \tag{ 1 }$$

where ϕ_int is the extrapolation of the Fermi-LAT measurement toward the H.E.S.S. energy range. ${{\phi }_{{\rm obs}}}\pm {\Delta }{{\phi }_{{\rm obs}}}$ is the flux measured by H.E.S.S. The factor (1 − α) = 0.8 accounts for the systematic uncertainties of the H.E.S.S. measurement and the number 1.64 has been calculated to have a confidence level of 95% (Abramowski et al. 2013a). The comparison is made at the H.E.S.S. decorrelation energy where the flux is best measured.

Figure 4 shows the 95% UL on τ_max. The resulting UL on the redshift is z < 0.43. This method does not allow the statistical and systematic uncertainties of the Fermi-LAT measurement to be taken into account and does not take advantage of the spectral features of the absorbed spectrum (see Abramowski et al. 2013b).

Zoom In Zoom Out Reset image size

A Bayesian approach has been developed with the aim of taking all the uncertainties into account. It also uses the fact that EBL-absorbed spectra are not strictly power laws. The details of the model are presented in Appendix A and only the main assumptions and results are recalled here. Intrinsic curvature between the HE and VHE ranges that naturally arises due to either curvature of the emitting distribution of particles or emission effects (e.g., Klein–Nishina effects) is permitted by construction of the prior (Equation (A.1)). A spectral index softer than the Fermi-LAT measurement is allowed with a constant probability, in contrast with the previous calculation. It is assumed that the observed spectrum in VHE γ-rays cannot be harder than the Fermi-LAT measurement by using a prior that follows a Gaussian for indices harder than the Fermi-LAT one. The prior on the index is then:

$$\begin{eqnarray}&&P({\Gamma })\propto {{\mathcal{N}}_{G}}({\Gamma },{{{\Gamma }}_{{\rm Fermi}}},{{\sigma }_{{\Gamma }}})\end{eqnarray} \tag{ 2 }$$

if ${\Gamma }$ < ${\Gamma }$_Fermi and

$$\begin{eqnarray*}&&P({\Gamma })\propto 1\end{eqnarray*}$$

otherwise. ${\Gamma }$_Fermi is the index measured by Fermi-LAT and ${{\sigma }_{{\Gamma }}}$ is the uncertainty on this measurement that takes all the systematic and statistical uncertainties into account.

The most probable redshift found with this method is z = 0.49 ± 0.04, in good agreement with the independent measure of Danforth et al. (2010), who constrained the distance to be between 0.43 < z < 0.58. Figure 5 gives the posterior probability obtained with the Bayesian method compared with other measurements of z. Lower and upper limits at a confidence level of 95% can also be derived as 0.41 < z < 0.56. Note that this method allows the systematic uncertainties of both instruments (Fermi-LAT and H.E.S.S.) to be taken into account. The spectral index obtained when fitting the H.E.S.S. data with an EBL absorbed PWL using a redshift of 0.49 is compatible with the Fermi measurement below 80 GeV.

Zoom In Zoom Out Reset image size

As stated in Section 2.1, the H.E.S.S. data of the flare show an indication of intra-night variability, which is used here to test for a possible LIV. Some Quantum Gravity (QG) models predict a change of the speed of light at energies close to the Planck scale ($\sim {{10}^{19}}$ GeV). A review of such models can be found in Mattingly (2005) and Liberati (2013). An energy-dependent dispersion in vacuum is searched for in the data by testing a correlation between arrival times of the photons and their energies. For two photons with arrival times t₁ and t₂ and energies E₁ and E₂, the dispersion parameter of order n is defined as ${{\tau }_{n}}=\frac{{{t}_{2}}-{{t}_{1}}}{E_{2}^{n}-E_{1}^{n}}=\frac{{\Delta }t}{{\Delta }({{E}^{n}})}$. Here only the linear (n = 1) and quadratic (n = 2) dispersion parameters are calculated. Assuming no intrinsic spectral variability of the source, the dispersion τ_n can be related to the normalized distance of the source κ_n corrected for the expansion of the universe and an energy E_QG at which QG effects are expected to occur (Jacob & Piran 2008):

$$\begin{eqnarray}&&{{\tau }_{n}}=\frac{{\Delta }t}{{\Delta }({{E}^{n}})}\simeq {{s}_{\pm }}\frac{(1+n)}{E_{{\rm QG}}^{n}2{{H}_{0}}}{{\kappa }_{n}}\end{eqnarray} \tag{ 3 }$$

where H₀ is the Hubble constant and s_± = −1 (resp. +1) in the superluminal (resp. subluminal) case, in which the high-energy photons arrive before (resp. after) low-energy photons. The normalized distance κ_n is calculated from the redshift of the source z and the cosmological parameters Ω_m, ${{{\Omega }}^{{\Lambda }}}$ given in the introduction:

$$\begin{eqnarray}&&{{\kappa }_{n}}=\int _{0}^{z}\frac{{{(1+z^{\prime} )}^{n}}\;dz^{\prime} }{\sqrt{{{{\Omega }}_{m}}{{(1+z^{\prime} )}^{3}}+{{{\Omega }}_{{\Lambda }}}}}\end{eqnarray} \tag{ 4 }$$

Using the central value of z = 0.49 determined in Section 3.2, the distance κ_n for n = 1 and 2 is κ₁ = 0.541 and κ₂ = 0.677.

First, the dispersion measurement method will be described. It will then be applied to the H.E.S.S. flare dataset (Monte Carlo (MC) simulations and original dataset), in order to measure the dispersion and provide 95% 1-sided lower and upper limits on the dispersion parameter τ_n. These limits on τ_n will lead to lower limits (LLs) on E_QG using Equation (3).

3.3.1. Modified Maximum Likelihood Method

A maximum likelihood method, following Martinez & Errando (2009), was used to calculate the dispersion parameter τ_n. Albert et al. (2008) applied this method to a flare of Mkn 501, while Abramowski et al. (2011) applied it to a flare of PKS 2155-304. More recently, it was used by Vasileiou et al. (2013) to analyze Fermi data of four gamma-ray bursts. The data from Cherenkov telescopes are contaminated by π⁰ decay from proton showers, misidentified electrons, or heavy elements such as helium. In the case of PG 1553+113, and contrary to previous analyses, this background is not negligible: the signal-over-background ratio S/B is about 2, compared to 300 for the PKS 2155-304 flare event of 2006 July (Aharonian et al. 2007b). The background was included in the formulation of the probability density function (PDF) used in a likelihood maximization method. Given the times t_i and energies E_i of the gamma-like (ON) particles received by the detector, the unbinned likelihood function of the dispersion parameter τ_n is:

$$\begin{eqnarray}&&L({{\tau }_{n}})={{\mathop{\prod }\limits_{i=1}^{{{n}_{{\rm ON}}}}}^{{}}}P({{E}_{i}},{{t}_{i}}|{{\tau }_{n}}).\end{eqnarray} \tag{ 5 }$$

The PDF $P({{E}_{i}},{{t}_{i}}|{{\tau }_{n}})$ associated with each ON event is composed of two terms:

$$\begin{eqnarray}\begin{array}{rcl} P({{E}_{i}},{{t}_{i}}|{{\tau }_{n}}) & = & {{w}_{s}}\cdot {{P}_{{\rm Sig}}}({{E}_{i}},{{t}_{i}}|{{\tau }_{n}}) \\ {} & {} & +(1-{{w}_{s}})\cdot {{P}_{{\rm Bkg}}}({{E}_{i}},{{t}_{i}}) \\ \end{array}\end{eqnarray} \tag{ 6 }$$

with

$$\begin{eqnarray}\begin{array}{rcl} {{P}_{{\rm Sig}}}({{E}_{i}},{{t}_{i}}|{{\tau }_{n}}) & = & \frac{1}{N({{\tau }_{n}})}\;{{A}_{{\rm eff}}}({{E}_{i}},{{t}_{i}}) \\ {} & {} & \times {{{\Lambda }}_{{\rm Sig}}}({{E}_{i}})\;{{F}_{{\rm Sig}}}({{t}_{i}}-{{\tau }_{n}}\cdot E_{i}^{n}) \\ \end{array}\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&{{P}_{{\rm Bkg}}}({{E}_{i}},{{t}_{i}})=\frac{1}{N^{\prime} }\;{{A}_{{\rm eff}}}({{E}_{i}},{{t}_{i}})\;{{{\Lambda }}_{{\rm Bkg}}}({{E}_{i}})\;{{F}_{{\rm Bkg}}}\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&{{w}_{s}}=\frac{{{n}_{{\rm ON}}}-\alpha \;{{n}_{{\rm OFF}}}}{{{n}_{{\rm ON}}}}.\end{eqnarray} \tag{ 9 }$$

The PDF P_Sig includes the emission time distribution of the photons F_Sig determined from a parameterization of the observed light curve at low energies (discussed in the next section) and evaluated on $t-{{\tau }_{n}}\cdot {{E}^{n}}$ to take into account the delay due to a possible LIV effect, the measured signal spectrum Λ_Sig and the effective area A_eff. The PDF P_Bkg is composed of the uniform time distribution F_Bkg of the background events, the measured background spectrum Λ_Bkg and the effective area A_eff. No delay due to a possible LIV effect is expected in the background events of the ON data set. N(τ_n) and N' are the normalization factors of P_Sig and P_Bkg respectively, in the (E, t) range of the likelihood fit. The coefficient w_s corresponds to the relative weight of the signal events in the total ON data set, derived from the number of events in the ON region n_ON and the number of events in the OFF regions n_OFF weighted by the inverse number of OFF regions α. More details on the derivation of this function are given in Appendix B.1.

3.3.2. Specific Selection Cuts and Timing Model

The flare data set of the H.E.S.S. analysis (see Section 2.1) was used with additional cuts. To perform the dispersion studies, only uninterrupted data have been kept. Thus, the analysis was conducted on the first seven runs, taken during the night of April 26. Moreover, the cosmic ray flux increases substantially for the seventh run, due to a variation of the zenith angle during this night. This fact, along with its large statistical errors, leads us to discard this run from the analysis. The sixth run shows little to no variability and was therefore also removed from the LIV analysis. Since within the ON data set, the signal and the background spectra have different indices (${{{\Gamma }}_{{\rm Sig}}}=4.8$ for the signal and ${{{\Gamma }}_{{\rm Bkg}}}=2.5$ for the background), the ratio S/B is expected to decrease with increasing energy. An upper energy cut at E_max = 789 GeV was set, corresponding to the last bin with more than 3σ significance in the reconstructed photon spectrum (see the differential flux during the flare in Figure 1). A lower cut on the energy at E_min = 300 GeV was used in order to avoid large systematic effects arising from high uncertainties on the H.E.S.S. effective area at lower energies. The intrinsic light curve of the flare, needed in the formulation of the likelihood, can be obtained from a model of the timed emission or approximated from a subset of the data. To be as model-independent as possible, it was here derived from a fit of the measured light curve at low energies (with E < E_cut). The high-energy events ($E\gt {{E}_{{\rm cut}}}$) were processed in the calculation of the likelihood to search for potential dispersion. Here E_cut was set to E_cut = 400 GeV, which is approximately the median energy of the ON event sample. Other cuts on the energy did not introduce significant effects on the final results. The histogram and the fit (Figure 6) were obtained as follows: the main idea was to preserve the maximum detected variability in the PG 1553+113 flare, together with a significant response in each observed peak:

• 1.  
  The binning was chosen so that at least two adjacent bins of the distribution yield a minimum of 3σ excess with respect to the average value.
• 2.  
  Simple parameterizations have been tested on the whole data set (all energies): constant (χ²/dof = 25/12), single Gaussian (χ²/dof = 20/10) and double Gaussian (χ²/dof = 8.5/7) functions. The latter is preferred, since it improves the quality of the fit. This shape was chosen to fit the low energy subset of events. Choosing a single Gaussian parameterization would result in a decrease of the sensitivity to time-lag measurements by a factor of two.

Zoom In Zoom Out Reset image size

There is a gap of ∼2 minutes between each two consecutive runs. We did not consider the effect of these gaps as it is small with respect to the bin width of ∼10 minutes. More importantly, their occurrence is not correlated with the binning: one gap falls in the rising part of the light curve, one is at a maximum, two fall in the decreasing parts and none of the gaps is at the minimum.

Table B1 in Appendix B.2 shows the number of ON and OFF events for the different cuts applied to the data.

Table B1.  Selections Applied to the ON and OFF Data Sets

Selection	# of ON Events	Weighted # of OFF Events	S/B
Total sample	461 (100%)	144.3 (100%)	2.2
(1) = Time in 500–8500s	358 (77.7%)	95.8 (66.4%)	2.7
(1) and E in 0.3–0.789 TeV	154 (33.4%)	36.3 (25.1%)	3.2
(1) and E in 0.3–0.4 TeV (Template)	82 (17.8%)	14.2 (9.9%)	4.8
(1) and E in 0.4–0.789 TeV (LH fit)	72 (15.6%)	21.9 (15.2%)	2.3

Download table as:  ASCIITypeset image

3.3.3. Results: Limits on τ_n and E_QG

The maximum likelihood method was performed using high-energy events with E_i > E_cut. First, confidence intervals (CIs) corresponding to 95% confidence level (1-sided) were determined from the likelihood curve at the values of τ_n where the curve reaches 2.71, which corresponds to the 90% CL quantile of a χ² distribution. However, these CIs are derived from one realization only and do not take into account the "luckiness" factor of this measurement. To get statistically significant CIs ("calibrated CIs"), several sets were generated with MC simulations, with the same statistical significance, light curve model and spectrum as the original data set. No intrinsic dispersion was artificially added. Each simulated data set produces a LL and an UL on τ_n. The calibrated lower (upper) limit of the confidence interval is obtained from the mean of the distribution of the per-set individual lower (upper) limits. Both CIs (from the data only and from the simulated sets) are listed in Table B2. Sources of systematic errors include uncertainties on the light curve parameterization, the background contribution, the calculation of the effective area, the energy resolution, and the determination of the photon index (see Appendix B.4).

Table B2.  Linear (Top) and Quadratic (Bottom) Dispersion Parameters; from left to right: Best Estimate, LL and UL from Data (Cut on Likelihood Curve), LL and UL from MC Simulations (means of Per-set LL and UL Distributions), Calibrated LL and UL (Combination of Data and MC), Calibrated LL and UL Including Systematic Errors

n	$\tau _{n,{\rm best}}^{{\rm data}}$	${\rm LL}_{n}^{{\rm data}}$	${\rm UL}_{n}^{{\rm data}}$	$\tau _{n,{\rm best}}^{{\rm MC}}$	${\rm LL}_{n}^{{\rm MC}}$	${\rm UL}_{n}^{{\rm MC}}$	${\rm LL}_{n}^{{\rm calib}}$	${\rm UL}_{n}^{{\rm calib}}$	${\rm LL}_{n}^{{\rm calib}}$	${\rm UL}_{n}^{{\rm calib}}$
									With Systematics
1	−131.7	−806.7	554.7	99.1	−526.3	725.6	−757.1	494.8	−838.9	576.4
2	−287.5	−1449.9	853.6	217.2	−942.0	1395.0	−1446.7	890.3	−1570.5	1012.4

Note. Dispersion Parameters τ_n,best, LLs and ULs are in s TeV⁻ⁿ.

Download table as:  ASCIITypeset image

The resulting limits on the dispersion τ_n using the quadratic sum of the statistical errors from the simulations and the systematic errors determined from data and simulations were computed, leading to limits on the energy scale E_QG (Equation (3)). The 95% 1-sided LLs for the subluminal case (s = +1) are: E_QG,1 > 4.11 × 10¹⁷ GeV and ${{E}_{{\rm QG},2}}\gt 2.10\times {{10}^{10}}$ GeV for linear and quadratic LIV effects, respectively. For the superluminal case (s = –1) the limits are: ${{E}_{{\rm QG},1}}\gt 2.83\times {{10}^{17}}$ GeV and E_QG,2 > 1.68 × 10¹⁰ GeV for linear and quadratic LIV effects, respectively. Figure 7 shows a comparison of the different LLs on E_QG,1 and E_QG,2 for the subluminal case (s = +1) obtained with AGNs at different redshifts studied at VHE. All these limits, including the present results, have been obtained under the assumption that no intrinsic delays between photons of different energies occur at the source. For the linear/subluminal case, the most constraining limit on E_QG with transient astrophysical events has been obtained with GRB 090510: E_QG,1 > 6.3 × 10¹⁹ GeV (Vasileiou et al. 2013). The most constraining limits on E_QG with AGN so far have been obtained by Abramowski et al. (2011) with PKS 2155-304 data observed with H.E.S.S.: ${{E}_{{\rm QG},1}}\gt 2.1\times {{10}^{18}}$ GeV and E_QG,2 > 6.4 × 10¹⁰ GeV for linear and quadratic LIV effects, respectively (95% CL, 1-sided). Compared to the PKS 2155-304 limits, the limits on the linear dispersion for PG 1553+113 are one order of magnitude less constraining, but the limits on the quadratic dispersion are of the same order of magnitude since the source is located at a higher redshift. This highlights the interest in studying distant AGNs, in spite of the difficulties due to limited photon statistics.

Zoom In Zoom Out Reset image size

In previous LIV studies with AGN flares (Albert et al. 2008; Abramowski et al. 2011) the signal was clearly dominating over the background, whereas in the present study the signal-over-background ratio is about 2. The background has been included in the formulation of the PDF: in the most general case, for given numbers of signal and background events s and b in the observation region ("ON" region), for a given dispersion parameter τ_n, the unbinned likelihood is:

$$\begin{eqnarray}\begin{array}{rcl} L({{n}_{{\rm ON}}},{{n}_{{\rm OFF}}}|s,b,{{\tau }_{n}}) & = & {\rm Pois}({{n}_{{\rm ON}}}|s+b)\cdot {\rm Pois}\left( {{n}_{{\rm OFF}}}|\frac{b}{\alpha } \right) \\ {} & {} & \cdot {{\mathop{\prod }\limits_{i=1}^{{{n}_{{\rm ON}}}}}^{{}}}P({{E}_{i}},{{t}_{i}}|s,b,{{\tau }_{n}}) \\ \end{array}\end{eqnarray} \tag{ B.1 }$$

The PDF $P({{E}_{i}},{{t}_{i}}|s,b,{{\tau }_{n}})$ associated with each gamma-like particle characterized by its time t_i and energy E_i contains two terms (signal and background):

$$\begin{eqnarray}\begin{array}{rcl} P({{E}_{i}},{{t}_{i}}|s,b,{{\tau }_{n}}) & = & {{w}_{s}}\cdot {{P}_{{\rm Sig}}}({{E}_{i}},{{t}_{i}}|{{\tau }_{n}}) \\ {} & {} & +(1-{{w}_{s}})\cdot {{P}_{{\rm Bkg}}}({{E}_{i}},{{t}_{i}}) \\ \end{array}\end{eqnarray} \tag{ B.2 }$$

with

$$\begin{eqnarray}&&{{w}_{s}}=\frac{s}{s+b}.\end{eqnarray} \tag{ B.3 }$$

n_ON is the number of events detected in the source ON region included in the fit range $[{{E}_{{\rm cut}}};{{E}_{{\rm max} }}]\times [{{t}_{{\rm min} }};{{t}_{{\rm max} }}]$. n_OFF is the number of events in the OFF regions, in the same (E, t) range; α is the inverse number of OFF regions. ${\rm Pois}({{n}_{{\rm ON}}}|s+b)$ (${\rm Pois}({{n}_{{\rm OFF}}}|b/\alpha )$) is the Poisson distribution with index n_ON (n_OFF) and parameter s + b (b/α). The likelihood function can be simplified by fixing s and b from a comparison of ON and OFF sets: s = n_ON − αn_OFF and b = αn_OFF. In this case, the Poisson terms in Equation (B.2) are equal to 1. The probabilities P_Sig and P_Bkg are defined as:

$$\begin{eqnarray}&&{{P}_{{\rm Sig}}}({{E}_{i}},{{t}_{i}}|{{\tau }_{n}})=\frac{1}{N({{\tau }_{n}})}\cdot {{R}_{{\rm Sig}}}({{E}_{i}},{{t}_{i}}|{{\tau }_{n}})\end{eqnarray} \tag{ B.4 }$$

$$\begin{eqnarray}&&{{P}_{{\rm Bkg}}}({{E}_{i}},{{t}_{i}})=\frac{1}{N^{\prime} }\cdot {{R}_{{\rm Bkg}}}({{E}_{i}},{{t}_{i}})\end{eqnarray} \tag{ B.5 }$$

with

$$\begin{eqnarray}\begin{array}{rcl} {{R}_{{\rm Sig}}}(E,t|{{\tau }_{n}}) & = & \int _{{{E}_{{\rm true}}}=0}^{\infty }D(E,{{E}_{{\rm true}}}){{A}_{{\rm eff}}}({{E}_{{\rm true}}},t) \\ {} & {} & \times {{{\Lambda }}_{{\rm Sig}}}({{E}_{{\rm true}}}){{F}_{{\rm Sig}}}(t-{{\tau }_{n}}\cdot E_{{\rm true}}^{n})d{{E}_{{\rm true}}} \\ {} & {} & \quad \\ \end{array}\end{eqnarray} \tag{ B.6 }$$

$$\begin{eqnarray}\begin{array}{rcl} {{R}_{{\rm Bkg}}}(E,t) & = & \int _{{{E}_{{\rm true}}}=0}^{\infty }D(E,{{E}_{{\rm true}}}) \\ {} & {} & \times {{A}_{{\rm eff}}}({{E}_{{\rm true}}},t){{{\Lambda }}_{{\rm Bkg}}}({{E}_{{\rm true}}}){{F}_{{\rm Bkg}}}(t)d{{E}_{{\rm true}}}. \\ {} & {} & \quad \\ \end{array}\end{eqnarray} \tag{ B.7 }$$

${{P}_{{\rm Sig}}}({{E}_{i}},{{t}_{i}}|{{\tau }_{n}})$ is the probability that the event (E_i, t_i) is a photon emitted at the source and detected on Earth with a delay ${{\tau }_{n}}{{E}^{n}}$. It takes into account the emission (time distribution ${{F}_{{\rm Sig}}}(t)$ and energy spectrum ${{{\Lambda }}_{{\rm Sig}}}(E)$ at the source), the propagation (delay ${{\tau }_{n}}\cdot E_{i}^{n}$ due to possible LIV effect) and the detection of a photon by the detector (H.E.S.S. energy resolution $D(E,{{E}_{{\rm true}}})$ and effective area ${{A}_{{\rm eff}}}(E,t)$). ${{P}_{{\rm Bkg}}}({{E}_{i}},{{t}_{i}})$ is the probability that the event (E_i, t_i) is a background event; it is not expected to be variable with time, thus F_Bkg(t) is a uniform time distribution: ${{F}_{{\rm Bkg}}}(t)={{F}_{{\rm Bkg}}}$. The background energy distribution Λ_Bkg is measured from OFF regions. N(τ_n) (resp. N') is the normalization factor of the PDF P_Sig (resp. P_Bkg) in the range $[{{E}_{{\rm cut}}};{{E}_{{\rm max} }}]\times [{{t}_{{\rm min} }};{{t}_{{\rm max} }}]$ where the likelihood fit is performed.

Also, the energy resolution D(E, E_true) is assumed to be perfect in the range [E_cut; E_max].⁵³ This leads to simplified expressions of ${{P}_{{\rm Sig}}}({{E}_{i}},{{t}_{i}}|{{\tau }_{n}})$ and ${{P}_{{\rm Bkg}}}({{E}_{i}},{{t}_{i}})$:

$$\begin{eqnarray}\begin{array}{rcl} {{P}_{{\rm Sig}}}({{E}_{i}},{{t}_{i}}|{{\tau }_{n}}) & = & \frac{1}{N({{\tau }_{n}})}\cdot {{A}_{{\rm eff}}}({{E}_{i}},{{t}_{i}}) \\ {} & {} & \times {{{\Lambda }}_{{\rm Sig}}}({{E}_{i}}){{F}_{{\rm Sig}}}({{t}_{i}}-{{\tau }_{n}}\cdot E_{i}^{n}) \\ \end{array}\end{eqnarray} \tag{ B.8 }$$

$$\begin{eqnarray}&&{{P}_{{\rm Bkg}}}({{E}_{i}},{{t}_{i}})=\frac{1}{N^{\prime} }\cdot {{A}_{{\rm eff}}}({{E}_{i}},{{t}_{i}}){{{\Lambda }}_{{\rm Bkg}}}({{E}_{i}}){{F}_{{\rm Bkg}}}\end{eqnarray} \tag{ B.9 }$$

The best estimate of the dispersion parameter ${{\hat{\tau }}_{n}}$ is obtained by maximizing the likelihood L(τ_n).

Table B1 shows the effect of the selection cuts on the number of ON and OFF events. Other choices of E_min and E_cut did not introduce significant changes in the final results.

The method has been tested on MC simulated sets. Each set was composed of ${{n}_{{\rm ON}}}=72$ ON events, as in the real data sample:

• 1.  
  s = 50 signal events with times following the template light curve (Figure 6) shifted by a factor ${{\tau }_{n,{\rm inj}}}\cdot {{E}_{i}}$; energies follow a power law spectrum of photon index ${{{\Gamma }}_{{\rm Sig}}}=4.8$, degraded by the acceptance and convolved with the energy resolution.
• 2.  
  b = 22 background events with times following a uniform distribution and energies drawn from a power law spectrum of index ${{{\Gamma }}_{{\rm Bkg}}}=2.5$, degraded by the acceptance and convolvted with the energy resolution.

For a given injected dispersion, the maximum likelihood method is applied to each MC-simulated set. The initial light curve and energy spectrum were used as templates in the model instead of fitting them for each set.

Figure B1 shows the means of the reconstructed dispersion versus the real (injected) dispersion for n = 1; for a given injected dispersion, error bars correspond to the rms of the distribution of the best estimates $\widehat{{{\tau }_{1}}}$. The blue line shows the result of a linear fit. The slope roughly corresponds to the percentage of signal in the total ON data set. It is due to the loss of sensitivity resulting from the part of the data sets with no dispersion. A systematic shift is observed of about 100 s TeV⁻¹, well bellow 1σ value—the rms of the best estimate distribution is 361 s TeV⁻¹. The results in this paper have not been corrected for this bias.

Zoom In Zoom Out Reset image size

The coverage is not necessarily proper, i.e. the number of sets for which the injected dispersion value τ_inj lies between the set's LL and UL does not match the required 95% 1-sided confidence level. The common cut used on the likelihood curves to get the LLs/ULs has been iteratively adjusted to ensure a correct statistical coverage: using this new cut, 95% of the realizations provide CIs that include the injected dispersion ${{\tau }_{n,{\rm inj}}}$. The initial coverage was about 85% for a cut on $2{\rm ln} L$ of 2.71. The new common cut, found iteratively at 3.5, ensures the desired 90% 2-sided CL (approx. 95% one-sided CL). Figure B2 shows the distributions of the best estimates, the 95% 1-sided LLs and ULs for ${{\tau }_{1,{\rm inj}}}=0$ s TeV⁻¹ (linear case) and ${{\tau }_{2,{\rm inj}}}=0$ s TeV⁻² (quadratic case); the means of the lower and upper limit distributions, shown as a blue vertical line, are used to construct the "calibrated confidence interval."

Zoom In Zoom Out Reset image size

To get CIs from data, a maximum likelihood method is applied to the original data set and gives a best estimate $\tau _{{\rm best}}^{{\rm data}}$. The cut value determined from the simulations to ensure proper coverage is applied on the original data set to obtain LL^data and UL^data. The "calibrated" limits LL^calib and UL^calib, combining $\tau _{{\rm best}}^{{\rm data}}$ from data together with MC results, are taken as

$$\begin{eqnarray}\begin{array}{rcl} {\rm L}{{{\rm L}}^{{\rm calib}}} & = & \tau _{{\rm best}}^{{\rm data}}-|\tau _{{\rm best}}^{{\rm MC}}-{\rm L}{{{\rm L}}^{{\rm MC}}}| \\ {\rm U}{{{\rm L}}^{{\rm calib}}} & = & \tau _{{\rm best}}^{{\rm data}}+|\tau _{{\rm best}}^{{\rm MC}}-{\rm U}{{{\rm L}}^{{\rm MC}}}| \\ \end{array}\end{eqnarray} \tag{ B.10 }$$

with $\tau _{{\rm best}}^{{\rm MC}}$, ${\rm L}{{{\rm L}}^{{\rm MC}}}$ and ${\rm U}{{{\rm L}}^{{\rm MC}}}$ defined as the means of the per-set best-estimate distribution, LL distribution, and UL distribution respectively.

Table B2 lists the CIs determined in both ways, i.e., data-only and calibrated ones: ${\rm LL}_{n}^{{\rm data}}$ and ${\rm LL}_{n}^{{\rm calib}}$ (resp. ${\rm UL}_{n}^{{\rm data}}$ and ${\rm UL}_{n}^{{\rm calib}}$) are compatible within 10%. In this work, calibrated CIs have been used to derive the final LLs on E_QG. They are preferred over data-only CIs as they provide statistically well defined confidence levels. They also ensure coherent comparison with previous published results, e.g., with PKS 2155-304 by Abramowski et al. (2011) and GRB studies by Vasileiou et al. (2013).

Estimations of the systematic effects on the dispersion measurement were performed. It was found that the main systematic errors are due to the uncertainties on the light curve parameterization. Other sources of systematic errors include the contribution of the background, effect of the change of photon index, the energy resolution and the effective area determination of the detector. To study the following four contributions, new simulated data sets have been built, each one with different input parameters:

• 1.  
  background contribution: photons and background events have been reallocated within the ON data set in the fit range [E_cut; E_max], introducing a 1σ fluctuation in the number of signal event s in the ON data set;
• 2.  
  effective area: set to a constant, equal to 120000 m² for all energies and all times, which corresponds to a maximum shift of 10% (the actual effective area increases with energy);
• 3.  
  energy resolution: reconstructed energies have been replaced by the true energies; this corresponds to a shift of about 10% on the reconstructed energy values;
• 4.  
  photon index: changed by one standard deviation (±0.25).

For the determination of systematic errors arising from the light curve parameterization, the calibration of the CIs has been redone using successively the upper 1σ and the lower 1σ contours of the template, shown in Figure 6. The change in mean lower and ULs on the dispersion parameter τ_n gives an estimate of the systematic error associated to each contribution.⁵⁴ An additional systematic contribution comes from the shift arising from the method found with simulation (see Appendix B.3). Table B3 summarizes all studied systematic contributions. The overall estimated systematic error on τ_n is 330 s TeV⁻¹ for the linear case (n = 1) and 555 s TeV⁻² for the quadratic case (n = 2); they were included in the calculation of the limits on E_QG by adding the statistical and the systematic errors in quadrature.

Table B3.  Summary of all Studied Systematic Contributions. The Main Systematic Errors are due to the Uncertainties on the Light Curve Parametrization

	Estimated Error	τ1	τ2
	on Input Parameters	(s TeV−1)	(s TeV−2)
Background contribution	⋯	<45	<80
Acceptance factors	10%	<1	<1
Energy resolution	10%	<55	<85
Photon index	5%	<55	<50
Light curve parameterization	⋯	<300	<500
Systematic bias	⋯	∼100	∼200
Total: $\sqrt{{{\sum }_{i}}{\rm syst}_{i}^{2}}$		<330	<555

Download table as:  ASCIITypeset image