The 2012 flare of PG 1553+113 seen with H.E.S.S. and Fermi-LAT
A. Abramowski, F. Aharonian, F. Ait Benkhali, A.G. Akhperjanian, E.O. Angüner, M. Backes, S. Balenderan, A. Balzer, A. Barnacka, Y. Becherini, J. Becker Tjus, D. Berge, S. Bernhard, K. Bernlöhr, E. Birsin, J. Biteau, M. Böttcher, C. Boisson, J. Bolmont, P. Bordas, J. Bregeon, F. Brun, P. Brun, M. Bryan, T. Bulik, S. Carrigan, S. Casanova, P.M. Chadwick, N. Chakraborty, R. Chalme-Calvet, R.C.G. Chaves, M. Chrétien, S. Colafrancesco, G. Cologna, J. Conrad, C. Couturier, Y. Cui, M. Dalton, I.D. Davids, B. Degrange, C. Deil, P. deWilt, A. Djannati-Ataï, W. Domainko, A. Donath, L.O'C. Drury, G. Dubus, K. Dutson, J. Dyks, M. Dyrda, T. Edwards, K. Egberts, P. Eger, P. Espigat, C. Farnier, S. Fegan, F. Feinstein, M.V. Fernandes, D. Fernandez, A. Fiasson, G. Fontaine, A. Förster, M. Füßling, S. Gabici, M. Gajdus, Y.A. Gallant, T. Garrigoux, G. Giavitto, B. Giebels, J.F. Glicenstein, D. Gottschall, M.-H. Grondin, M. Grudzińska, D. Hadsch, S. Häffner, J. Hahn, J. Harris, G. Heinzelmann, G. Henri, G. Hermann, O. Hervet, A. Hillert, J.A. Hinton, W. Hofmann, P. Hofverberg, M. Holler, D. Horns, A. Ivascenko, A. Jacholkowska, C. Jahn, M. Jamrozy, M. Janiak, F. Jankowsky, I. Jung, M.A. Kastendieck, K. Katarzyński, U. Katz, S. Kaufmann, B. Khélifi, M. Kieffer, et al. (127 additional authors not shown)
aa r X i v : . [ a s t r o - ph . H E ] J a n To be submitted to Astrophysical Journal
The 2012 flare of PG 1553+113 seen with H.E.S.S. and
Fermi -LAT: Constraints on the source redshift and Lorentzinvariance violation
H.E.S.S. Collaboration, A. Abramowski , F. Aharonian , , , F. Ait Benkhali ,A.G. Akhperjanian , , E.O. Ang¨uner , M. Backes , S. Balenderan , A. Balzer ,A. Barnacka , , Y. Becherini , J. Becker Tjus , D. Berge , S. Bernhard ,K. Bernl¨ohr , , E. Birsin , J. Biteau , , M. B¨ottcher , C. Boisson , J. Bolmont ,P. Bordas , J. Bregeon , F. Brun , ∗ , P. Brun , M. Bryan , T. Bulik , S. Carrigan ,S. Casanova , , P.M. Chadwick , N. Chakraborty , R. Chalme-Calvet , R.C.G. Chaves , M. Chr´etien , S. Colafrancesco , G. Cologna , J. Conrad , , C. Couturier , ∗ ,Y. Cui , M. Dalton , , I.D. Davids , , B. Degrange , C. Deil , P. deWilt ,A. Djannati-Ata¨ı , W. Domainko , A. Donath , L.O’C. Drury , G. Dubus , K. Dutson , J. Dyks , M. Dyrda , T. Edwards , K. Egberts , P. Eger , P. Espigat ,C. Farnier , S. Fegan , F. Feinstein , M.V. Fernandes , D. Fernandez , A. Fiasson , G. Fontaine , A. F¨orster , M. F¨ußling , S. Gabici , M. Gajdus , Y.A. Gallant ,T. Garrigoux , G. Giavitto , B. Giebels , J.F. Glicenstein , D. Gottschall ,M.-H. Grondin , , M. Grudzi´nska , D. Hadsch , S. H¨affner , J. Hahn , J. Harris ,G. Heinzelmann , G. Henri , G. Hermann , O. Hervet , A. Hillert , J.A. Hinton ,W. Hofmann , P. Hofverberg , M. Holler , D. Horns , A. Ivascenko , A. Jacholkowska , C. Jahn , M. Jamrozy , M. Janiak , F. Jankowsky , I. Jung , M.A. Kastendieck , K. Katarzy´nski , U. Katz , S. Kaufmann , B. Kh´elifi , M. Kieffer , S. Klepser ,D. Klochkov , W. Klu´zniak , D. Kolitzus , Nu. Komin , K. Kosack , S. Krakau ,F. Krayzel , P.P. Kr¨uger , H. Laffon , G. Lamanna , J. Lefaucheur , ∗ , V. Lefranc , A. Lemi`ere , M. Lemoine-Goumard , J.-P. Lenain , ∗ , T. Lohse , A. Lopatin ,C.-C. Lu , V. Marandon , A. Marcowith , R. Marx , G. Maurin , N. Maxted ,M. Mayer , T.J.L. McComb , J. M´ehault , , P.J. Meintjes , U. Menzler , M. Meyer , A.M.W. Mitchell , R. Moderski , M. Mohamed , K. Mor˚a , E. Moulin ,T. Murach , M. de Naurois , J. Niemiec , S.J. Nolan , L. Oakes , H. Odaka , S. Ohm , B. Opitz , M. Ostrowski , I. Oya , M. Panter , R.D. Parsons , M. Paz Arribas ,N.W. Pekeur , G. Pelletier , J. Perez , P.-O. Petrucci , B. Peyaud , S. Pita ,H. Poon , G. P¨uhlhofer , M. Punch , A. Quirrenbach , S. Raab , I. Reichardt ,A. Reimer , O. Reimer , M. Renaud , R. de los Reyes , F. Rieger , L. Rob ,C. Romoli , S. Rosier-Lees , G. Rowell , B. Rudak , C.B. Rulten , V. Sahakian , ,D. Salek , D.A. Sanchez , ∗ , A. Santangelo , R. Schlickeiser , F. Sch¨ussler ,A. Schulz , U. Schwanke , S. Schwarzburg , S. Schwemmer , H. Sol , F. Spanier , 2 –G. Spengler , F. Spies , L. Stawarz , R. Steenkamp , C. Stegmann , , F. Stinzing ,K. Stycz , I. Sushch , , J.-P. Tavernet , T. Tavernier , A.M. Taylor , R. Terrier ,M. Tluczykont , C. Trichard , K. Valerius , C. van Eldik , B. van Soelen ,G. Vasileiadis , J. Veh , C. Venter , A. Viana , P. Vincent , J. Vink , H.J. V¨olk ,F. Volpe , M. Vorster , T. Vuillaume , P. Wagner , R.M. Wagner , M. Ward ,M. Weidinger , Q. Weitzel , R. White , A. Wierzcholska , P. Willmann ,A. W¨ornlein , D. Wouters , R. Yang , V. Zabalza , , D. Zaborov , M. Zacharias ,A.A. Zdziarski , A. Zech , H.-S. Zechlin * Corresponding authors:D.A. Sanchez, [email protected], F. Brun, [email protected], C. Couturier,[email protected], J. Lefaucheur, [email protected], J.-P. Lenain, [email protected] Universit¨at Hamburg, Institut f¨ur Experimentalphysik, Luruper Chaussee 149, D 22761 Hamburg, Ger-many Max-Planck-Institut f¨ur Kernphysik, P.O. Box 103980, D 69029 Heidelberg, Germany Dublin Institute for Advanced Studies, 31 Fitzwilliam Place, Dublin 2, Ireland National Academy of Sciences of the Republic of Armenia, Marshall Baghramian Avenue, 24, 0019Yerevan, Republic of Armenia Yerevan Physics Institute, 2 Alikhanian Brothers St., 375036 Yerevan, Armenia Institut f¨ur Physik, Humboldt-Universit¨at zu Berlin, Newtonstr. 15, D 12489 Berlin, Germany University of Namibia, Department of Physics, Private Bag 13301, Windhoek, Namibia University of Durham, Department of Physics, South Road, Durham DH1 3LE, U.K. GRAPPA, Anton Pannekoek Institute for Astronomy, University of Amsterdam, Science Park 904, 1098XH Amsterdam, The Netherlands Obserwatorium Astronomiczne, Uniwersytet Jagiello´nski, ul. Orla 171, 30-244 Krak´ow, Poland now at Harvard-Smithsonian Center for Astrophysics, 60 Garden St, MS-20, Cambridge, MA 02138,USA Department of Physics and Electrical Engineering, Linnaeus University, 351 95 V¨axj¨o, Sweden Institut f¨ur Theoretische Physik, Lehrstuhl IV: Weltraum und Astrophysik, Ruhr-Universit¨at Bochum,D 44780 Bochum, Germany GRAPPA, Anton Pannekoek Institute for Astronomy and Institute of High-Energy Physics, Universityof Amsterdam, Science Park 904, 1098 XH Amsterdam, The Netherlands Institut f¨ur Astro- und Teilchenphysik, Leopold-Franzens-Universit¨at Innsbruck, A-6020 Innsbruck, Aus-tria Laboratoire Leprince-Ringuet, Ecole Polytechnique, CNRS/IN2P3, F-91128 Palaiseau, France now at Santa Cruz Institute for Particle Physics, Department of Physics, University of California atSanta Cruz, Santa Cruz, CA 95064, USA Centre for Space Research, North-West University, Potchefstroom 2520, South Africa LUTH, Observatoire de Paris, CNRS, Universit´e Paris Diderot, 5 Place Jules Janssen, 92190 Meudon,France LPNHE, Universit´e Pierre et Marie Curie Paris 6, Universit´e Denis Diderot Paris 7, CNRS/IN2P3, 4Place Jussieu, F-75252, Paris Cedex 5, France Institut f¨ur Astronomie und Astrophysik, Universit¨at T¨ubingen, Sand 1, D 72076 T¨ubingen, Germany Laboratoire Univers et Particules de Montpellier, Universit´e Montpellier 2, CNRS/IN2P3, CC 72, PlaceEug`ene Bataillon, F-34095 Montpellier Cedex 5, France DSM/Irfu, CEA Saclay, F-91191 Gif-Sur-Yvette Cedex, France Astronomical Observatory, The University of Warsaw, Al. Ujazdowskie 4, 00-478 Warsaw, Poland Instytut Fizyki J¸adrowej PAN, ul. Radzikowskiego 152, 31-342 Krak´ow, Poland School of Physics, University of the Witwatersrand, 1 Jan Smuts Avenue, Braamfontein, Johannesburg,2050 South Africa Landessternwarte, Universit¨at Heidelberg, K¨onigstuhl, D 69117 Heidelberg, Germany Oskar Klein Centre, Department of Physics, Stockholm University, Albanova University Center, SE-10691 Stockholm, Sweden Wallenberg Academy Fellow, Universit´e Bordeaux 1, CNRS/IN2P3, Centre d’´Etudes Nucl´eaires de Bordeaux Gradignan, 33175Gradignan, France Funded by contract ERC-StG-259391 from the European Community, School of Chemistry & Physics, University of Adelaide, Adelaide 5005, Australia APC, AstroParticule et Cosmologie, Universit´e Paris Diderot, CNRS/IN2P3, CEA/Irfu, Observatoirede Paris, Sorbonne Paris Cit´e, 10, rue Alice Domon et L´eonie Duquet, 75205 Paris Cedex 13, France Univ. Grenoble Alpes, IPAG, F-38000 Grenoble, FranceCNRS, IPAG, F-38000 Grenoble, France Department of Physics and Astronomy, The University of Leicester, University Road, Leicester, LE17RH, United Kingdom Nicolaus Copernicus Astronomical Center, ul. Bartycka 18, 00-716 Warsaw, Poland Institut f¨ur Physik und Astronomie, Universit¨at Potsdam, Karl-Liebknecht-Strasse 24/25, D 14476 Pots-dam, Germany Laboratoire d’Annecy-le-Vieux de Physique des Particules, Universit´e de Savoie, CNRS/IN2P3, F-74941Annecy-le-Vieux, France DESY, D-15738 Zeuthen, Germany Universit¨at Erlangen-N¨urnberg, Physikalisches Institut, Erwin-Rommel-Str. 1, D 91058 Erlangen, Ger-many Centre for Astronomy, Faculty of Physics, Astronomy and Informatics, Nicolaus Copernicus University,Grudziadzka 5, 87-100 Torun, Poland Department of Physics, University of the Free State, PO Box 339, Bloemfontein 9300, South Africa Charles University, Faculty of Mathematics and Physics, Institute of Particle and Nuclear Physics, V
ABSTRACT
Very high energy (VHE,
E >
100 GeV) γ -ray flaring activity of the high-frequency peaked BL Lac object PG 1553+113 has been detected by the H.E.S.S. telescopes. The flux of the source increased by a factor of 3 during the nightsof 2012 April 26 and 27 with respect to the archival measurements with hintof intra-night variability. No counterpart of this event has been detected inthe
Fermi -LAT data. This pattern is consistent with VHE γ ray flaring beingcaused by the injection of ultrarelativistic particles, emitting γ rays at the highestenergies. The dataset offers a unique opportunity to constrain the redshift ofthis source at z = 0 . ± .
04 using a novel method based on Bayesian statistics.The indication of intra-night variability is used to introduce a novel method toprobe for a possible Lorentz Invariance Violation (LIV), and to set limits on theenergy scale at which Quantum Gravity (QG) effects causing LIV may arise.For the subluminal case, the derived limits are E QG , > . × GeV andE QG , > . × GeV for linear and quadratic LIV effects, respectively.
Subject headings:
Galaxies: active – BL Lacertae objects: Individual: PG 1553+113– Gamma rays: observations – Quantum Gravity – Lorentz invariance breaking
1. Introduction
Blazars are active galactic nuclei (AGN) with their jets closely aligned with the lineof sight to the Earth (Urry & Padovani 1995). Among their particularities is flux variabil-ity at all wavelengths on various time scales, from years down to (in some cases) minutes(Gaidos et al. 1996; Aharonian et al. 2007a). Flaring activity of blazars is of great interest forprobing the source-intrinsic physics of relativistic jets, relativistic particle acceleration andgeneration of high-energy radiation, as well as for conducting fundamental physics tests. Onthe one hand, exploring possible spectral variability between flaring and stationary stateshelps to understand the electromagnetic emission mechanisms at play in the jet. On theother hand, measuring the possible correlation between photon energies and arrival timesallows one to test for possible Lorentz invariance violation (LIV) leading to photon-energy-dependent variations in the speed of light in vacuum.
Holeˇsoviˇck´ach 2, 180 00 Prague 8, Czech Republic GRAPPA, Institute of High-Energy Physics, University of Amsterdam, Science Park 904, 1098 XHAmsterdam, The Netherlands F /F ) > − .
5, Osterman et al. 2006), which placesit among the most extreme HBLs (Rector et al. 2003). The object was observed in X-rays by multiple instruments in different flux states. Its 2–10 keV energy flux ranges from0 . × − erg cm − s − (Osterman et al. 2006) to 3 . × − erg cm − s − (Reimer et al.2008) but no fast variability (in the sub-hour time scale) has been detected so far.PG 1553+113 was discovered at very high energies (VHE, E >
100 GeV) by H.E.S.S.(Aharonian et al. 2006a, 2008) with a photon index of Γ = 4 . ± .
6. At high energies(HE, 100 MeV < E <
300 GeV) the source has been detected by the
Fermi -LAT (Abdo et al.2009b, 2010a) with a very hard photon index of Γ = 1 . ± .
03, making this object theone with the largest HE – VHE spectral break (∆Γ ≈ . Fermi -LAT was found by Abdo et al. (2009b, 2010a) on daily or weekly time scales, butusing an extended data set of 17 months, Aleksi´c et al. (2012) reported variability above1 GeV with flux variations of a factor of ∼ γ rays with only modest flux variations (from 4 to 11 % of the CrabNebula flux). In addition to the high X-ray variability, this behavior can be interpreted asevidence for Klein-Nishina effects (Abdo et al. 2010a) in the framework of a synchrotron self-Compton model. The source underwent VHE γ -ray flares in 2012 March (Cortina 2012a)and April (Cortina 2012b), detected by the MAGIC telescopes. During the March flare,the source was at a flux level of about 15% of that of the Crab Nebula, while in April itreached ≈ γ -ray flares, also a brightening in X-ray, UV and op-tical wavelengths has been noticed by the MAGIC collaboration. A detailed study of theMAGIC telescopes and multi-wavelength data is in press (Aleksi´c et al. 2014). The latterevent triggered the H.E.S.S. observations reported in this work.Despite several attempts to measure it, the redshift of PG 1553+113 still suffers from un-certainties. Different attempts, including optical spectroscopy (Treves et al. 2007; Aharonian et al.2008) or comparisons of the HE and VHE spectra of PG 1553+113 (Prandini et al. 2009;Sanchez et al. 2013), were made. Based on the assumption that the EBL-corrected VHEspectral index is equal to the
Fermi -LAT one, Prandini et al. (2009) derived an upper limitof z < .
67. Comparing PG 1553+113 statistically with other known VHE emitters and tak-ing into account a possible intrinsic γ -ray spectral break through a simple emission model,Sanchez et al. (2013) constrained the redshift to be below 0.64. The best estimate to-datewas obtained by Danforth et al. (2010) who found the redshift to be between 0.43 and 0.58 7 –using far-ultraviolet spectroscopy.This paper concentrates on the HE and VHE emission of PG 1553+113 and is dividedas follow: Sections 2.1 and 2.2 present the H.E.S.S. and
Fermi -LAT analyses. The discus-sion, in section 3, includes the determination of the redshift using a novel method and theconstraints derived on LIV using a modified likelihood formulation. Throughout this papera ΛCDM cosmology with H = 70 . ± . − Mpc − , Ω m = 0 . ± .
03, Ω Λ = 0 . ± .
2. Data analysis2.1. H.E.S.S. observations and analysis
H.E.S.S. is an array of five imaging atmospheric Cherenkov telescopes located in theKhomas highland in Namibia (23 ◦ ′ ′′ S, 16 ◦ ′ ′′ E), at an altitude of 1800 m above sealevel (Hinton & the H.E.S.S. Collaboration 2004). The fifth
H.E.S.S. telescope was addedto the system in 2012 July and is not used in this work, reporting only on observations priorto that time.PG 1553+113 was observed with
H.E.S.S. in 2005 and 2006 (Aharonian et al. 2008). Novariability was found in these observations, which will be referred to as the “pre-flare” dataset in the following. New observations were carried out in 2012 April after flaring activityat VHE was reported by the MAGIC collaboration (“flare” data set, Cortina 2012b).The pre-flare data set is composed of 26 . ∼
28 minutes each were taken during the nightsof 2012 April 26 and 27, corresponding to 3 . . ◦ Model analysis (de Naurois & Rolland 2009) with
Loosecuts . This method–based on the comparison of detected shower images with a pre-calculatedmodel–achieves a better rejection of hadronic air showers and a better sensitivity at lowerenergies than analysis methods based on Hillas parameters. The chosen cuts, best suitedfor sources with steep spectra such as PG 1553+113 , require a minimum image charge of40 photoelectrons, which provides an energy threshold of ∼
217 GeV for the pre-flare and PG 1553+113 has one of the steepest spectra measured at VHE. ∼
240 GeV for the flare data set . All the results presented in this paper were cross-checkedwith 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, α J2000 = 15 h m . s , δ J2000 = 11 ◦ ′ . ′′ (Green et al. 1986), with a maximum squaredangular distance of 0 . , are used for the analysis. In order to estimate the backgroundin this region, the reflected background method (Berge et al. 2007) is used to define the OFFregions. The excess of γ rays in the ON region is statistically highly significant (Li & Ma1983): 21 . σ for the pre-flare period and 22 . σ for the flare. Statistics are summarized inTable 1.The differential energy spectrum of the VHE γ -ray emission has been derived usinga forward-folding method (Piron et al. 2001). For the observations prior to 2012 April, apower law (PWL) model fitted to the data gives a χ of 51 . χ probability of P χ = 0 . , with a χ of 37 . P χ =0 . . σ using the log-likelihoodratio test. Note that systematic uncertainties, presented in Table 2, have been evaluated byAharonian et al. (2006b) for the PWL model and using the jack-knife method for the LPmodel. The jack-knife method consist in removing one run and redoing the analysis. Thisprocess is repeated for all runs.For the flare data set, the log-parabola model does not significantly improve the fit andthe simple PWL model describes the data well, with a χ of 33 . P χ = 0 . ∼ χ of 123 . P χ = 6 . × − ). Restricting the The difference of energy threshold between the two data set is due to the changing observation conditions,e.g., zenith angle and optical efficiency. The log-parabola is defined by dN/dE = Φ ( E/E ) − a − b log( E/E ) . r ). The excess and the corresponding significance aregiven, as well as the energy threshold and the mean zenith angle of the source during theobservations. The last column presents the probability of the flux to be constant within theobservations (see text).Data set ON OFF r Excess Significance E th [GeV] Zenith angle P cst χ Pre-Flare 2205 13033 0.100 901.7 21.5 217 34 ◦ ◦ . × − True energy (TeV)1 ) - T e V - s - F l u x ( c m -13 -12 -11 -10 Pre-Flare data set ob s n σ ) / m ode l - n ob s ( n -4-3-2-101234 Reconstructed energy (TeV)0.3 0.4 0.5 0.6 0.7 1 2 True energy (TeV)1 ) - T e V - s - F l u x ( c m -13 -12 -11 -10 -9 Flare data set ob s n σ ) / m ode l - n ob s ( n -4-3-2-101234 Reconstructed energy (TeV)0.3 0.4 0.5 0.6 0.7 1 2
Fig. 1.— Differential fluxes of PG 1553+113 during the pre-flare (left) and flare (right)periods. Error contours indicate the 68 % uncertainty on the spectrum. Uncertainties onthe spectral points (in black) are given at 1 σ level, and upper limits are computed at the99 % confidence level. The gray squares were obtained by the cross-check analysis chain andare presented to visualize the match between both analyses. The gray error contour on theleft panel is the best-fit power law model. The lower panels show the residuals of the fit,i.e. the difference between the measured ( n obs ) and expected numbers of photons ( n model ),divided by the statistical error on the measured number of photons ( σ n obs ). 10 – Time [MJD]56044 56044.2 56044.4 56044.6 56044.8 56045 ] - s - c m - ( E > G e V ) [ φ Fig. 2.—
H.E.S.S. light curve of PG 1553+113 during the 2 nights of the flare period. Thecontinuous line is the measured flux during the flare period while the dashed one correspondsto the pre-flare period (see Table 2 for the flux values). Gray areas are the 1 σ errors.analysis to the pre-flare data set only, the fit yields a χ of 51.76 with 60 d.o.f. ( P χ = 0 . H.E.S.S. at the time of the flaring activityreported by Cortina (2012b).Figure 2 shows the light curve during the flare together with the averaged integral fluxesabove 300 GeV of both data sets. A fit with a constant to the
H.E.S.S. light curve duringthe first night yields a χ of 20 .
76 for 6 d.o.f. ( P χ = 2 . × − ), indicating intra-nightvariability. This is also supported by the use of a Bayesian block algorithm (Scargle 1998)that finds three blocks for the 2 nights at a 95% confidence level. 11 – Fermi -LAT analysis
The
Fermi
Large Area Telescope (LAT) is detector converting γ ray to e + e − pairs(Atwood et al. 2009). The LAT is sensitive to γ rays from 20 MeV to >
300 GeV. In surveymode, in which the bulk of the observations are performed, each source is seen every 3 hoursfor approximately 30 minutes.The
Fermi -LAT data and software are available from the
Fermi
Science Support Cen-ter . 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 theassociated
P7REP SOURCE V15 instrument response functions (IRFs). Events with energiesfrom 300 MeV to 300 GeV were selected. Additional cuts on the zenith angle ( < ◦ ) androcking angle ( < ◦ ) were applied as recommended by the LAT collaboration to reducethe contamination from the Earth atmospheric secondary radiation.The analysis of the LAT data was performed using the Enrico
Python package (Sanchez & Deil2013). The sky model was defined as a region of interest (ROI) of 15 ◦ radius with PG 1553+113in the center and additional point-like sources from the internal 4-years source list. Only thesources within a 3 ◦ radius around PG 1553+113 and bright sources (integral flux greater that5 × − 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 diffusecomponent, 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 aLP. A likelihood ratio test was used to decide which model best describes the data. Table 3gives the results for the two time periods considered in this work, and Figure 3 presentsthe γ -ray SEDs. The first one (pre-flare), before the H.E.S.S. exposures in 2012, includesmore than 3.5 years of data (from 2008 August 4 to 2012 March 1). The best fit model isfound to be the LP (with a TS of 11.3, ≈ . σ ). The second period (flare) is centered onthe H.E.S.S. observations windows and lasts for seven days. The best fit model is a powerlaw, the flux being consistent with the one measured during the first 3.5 years. Data pointsor light curves were computed within a restricted energy range or time range using a PWL http://fermi.gsfc.nasa.gov/ssc/data/analysis/ http://fermi.gsfc.nasa.gov/ssc/data/analysis/documentation/Cicerone/index.html Here the TS is 2 times the difference between the log-likelihood of the fit with a LP minus the log-likelihood with a PL.
12 –model with the spectral index frozen to 1.70.To precisely probe the variability in HE γ rays, seven-day time bins were used to computethe light curve of PG 1553+113 in an extended time window (from 2008 August 4 to 2012October 30), to probe any possible delay of a HE flare with respect to the VHE one. Whilethe flux of PG 1553+113 above 300 MeV is found to be variable in the whole period with avariability index of F var = 0 . ± .
04 (Vaughan et al. 2003), there is no sign of any flaringactivity around the 2012
H.E.S.S. observations. This result has been confirmed by usingthe 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 than1 GeV. No sign of enhancement of the HE flux associated to the VHE event reported herewas found. This might be due to the lack of statistic at high energy in the LAT energyrange.
E [TeV] -4 -3 -2 -1
10 1 ] - s - d N / d E [ e r g c m E -12 -11 -10 Flare data setPre-Flare data set
Fig. 3.— Spectral energy distribution of PG 1553+113 in γ rays as measured by the Fermi -LAT and H.E.S.S. Red (blue) points and butterflies have been obtained during the flare(pre-flare) period. The
Fermi and H.E.S.S. data for the pre-flare are not contemporaneous.H.E.S.S. data were taken in 2005-2006 while the
Fermi data were taken between 2008 and2012. 13 –Table 2. Summary of the fitted spectral parameters for the pre-flare and the flare datasets and the corresponding integral flux I calculated above 300 GeV. The last column givesthe decorrelation energy.Data Set (Model) Spectral Parameters I (E >
300 GeV) E dec [10 − ph cm − s − ] [GeV]Pre-Flare (PWL) Γ = 4 . ± . stat ± . sys . ± . stat ± . sys a = 5 . ± . stat ± . sys . ± . stat ± . sys · · · b = 4 . ± . stat ± . sys Flare (PWL) Γ = 4 . ± . stat ± . sys . ± . stat ± . sys Fermi -LAT data analysis for the pre-flare and flare periods. Forthe latter, the analysis has been performed in two energy ranges (see 3.2). The firstcolumns give the time and energy windows and the third the corresponding test statistic(TS) value. The model parameters and the flux above 300 MeV are given in the lastcolumns. The systematic uncertainties were computed using the IRFs bracketing method(Abdo et al. 2009a).MJD range Energy range TS Spectral Parameters I (E > − [ph cm − s − ]54682-55987 0.3-300 7793.7 a = 1 . ± . stat ± . sys . ± . stat ± . sys b = 3 . ± . stat ± . sys . ± . stat ± . sys . ± . stat ± . sys . ± . stat ± . sys . ± . stat ± . sys
14 –
3. Discussion of the results3.1. Variability in γ -rays The VHE data do not show any sign of variation of the spectral index (when comparingflare and pre-flare data sets with the same spectral model), and in HE no counterpart of thisevent 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 andVHE 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 log-parabolamodel better represents the pre-flare period in HE. Nonetheless, the precise location of thispeak cannot be determined with the Fermi -LAT data only. Combining both energy rangesand fitting the HE and VHE data points with a power law with an exponential cutoff allowsus to determine the νf ( ν ) peak position for both time periods. The functional form of themodel is E dNdE = N (cid:18) E
100 GeV (cid:19) − Γ exp( − E/E c ) . For this purpose,
Fermi -LAT and
H.E.S.S. systematic uncertainties were taken intoaccount in a similar way as in Abramowski et al. (2014) and added quadratically to the sta-tistical 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 andabove. 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 σ (f) sys = 0 . · ∂ f /∂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 log ( E max / . ± . stat ± . sys withno evidence of variation during the flare and no spectral variation. This is consistent withthe fact that no variability in HE γ rays was found during the H.E.S.S. observations. Thisis also in agreement with the fact that HBLs are less variable in HE γ rays than other BLLac objects (Abdo et al. 2010b), while numerous flares have been reported in the TeV band. A fit with a LP model has been attempted, but the power law with an exponential cutoff leads to abetter description of the data.
Table 4. Parametrization results of the two time periods (first column) obtained by combining
H.E.S.S. and
Fermi -LAT. The second column gives the normalization at 100 GeV, while the third and the fourth present thespectral index and cut-off energy of the fitted power law with an exponential cut-off. The last column is the peakenergy in a νf ( ν ) representation.Period N ( E =100 GeV) Γ log ( E c / ( E max / − [erg cm − s − ]Pre-Flare 9 . ± . stat ± . sys . ± . stat ± . sys . ± . stat ± . sys . ± . stat ± . sys Flare 13 . ± . stat ± . sys . ± . stat ± . sys . ± . stat ± . sys . ± . stat ± . sys
16 –
The extragalactic background light (EBL) is a field of UV to far infrared photons pro-duced by the thermal emission from stars and reprocessed starlight by dust in galaxies (seeHauser & Dwek 2001, for a review) that interacts with very high energy γ rays from sourcesat cosmological distances. As a consequence, a source at redshift z exhibits an observedspectrum φ obs ( E ) = φ int ( E ) × e − τ ( E,z ) where φ int ( E ) is the intrinsic source spectrum and τ is the optical depth due to interaction with the EBL. Since the optical depth increaseswith 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 opticaldepth τ as a function of redshift and energy. In this section, the data taken by both instru-ments during the flare period are used, with the Fermi -LAT analysis restricted to the range300 MeV < E <
80 GeV (see Table 3 for the results). In the modest redshift range of VHEemitters detected so far ( z ≤ . τ γγ ∼ . z = 0 . z . In the caseof PG 1553+113, for which the redshift is unknown, the effects of the EBL on the VHEspectrum might be used to derive constraints on its distance. Ideally, this would be doneby 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 intrinsicspectrum in the VHE regime, or at least, as a solid upper limit (assuming no hardening ofthe spectrum).Following the method used by Abramowski et al. (2013a), it has been assumed that theintrinsic spectrum of the source in the H.E.S.S. energy range cannot be harder than theextrapolation of the
Fermi-LAT measurement. From this, one can conclude that the opticaldepth cannot be greater than τ max ( E ), given by: τ max ( E ) = ln (cid:20) φ int (1 − α )( φ obs − . φ obs ) (cid:21) , (1)where φ int is the extrapolation of the Fermi-LAT measurement towards the H.E.S.S. energyrange. φ obs ± ∆ φ obs is the measured flux by H.E.S.S. The factor (1 − α ) = 0 . For sake of simplicity it is assumed here that the best-fit model is a power law, an assumption which istrue for most of the cases due to limited statistics in the VHE range
17 –calculated to have a confidence level of 95% (Abramowski et al. 2013a). The comparison ismade at the
H.E.S.S. decorrelation energy where the flux is best measured.
E [TeV] -1
10 1 10 ( E ) m a x τ O p t i c a l D ep t h M a x UL (95% CL)z= 0.43 (Fr08)
Fig. 4.— Values of τ max as a function of the photon energy. The black line is the 95% ULobtained with the H.E.S.S. data and the red line is the optical depth computed with themodel of Franceschini et al. (2008) for a redshift of 0.43. The blue line is the decorrelationenergy fo the H.E.S.S. analyse. The gray lines are the value of optical depth for differentredshift.Figure 4 shows the 95 % UL on τ max . The resulting upper limit on the redshift is z < .
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 spectralfeatures of the absorbed spectrum (see Abramowski et al. 2013b).A Bayesian approach has been developed with the aim of taking all the uncertaintiesinto 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 andresults are recalled here. Intrinsic curvature between the HE and VHE ranges that naturallyarises 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 (Eq. A1): A spectral index 18 –softer than the
Fermi -LAT measurement is allowed with a constant probability, in contrastwith the previous calculation. It is assumed that the observed spectrum in VHE γ rayscannot be harder than the Fermi -LAT measurement by using a prior that follows a Gaussianfor indices harder than the
Fermi -LAT one. The prior on the index is then: P (Γ) ∝ N G (Γ , Γ Fermi , σ Γ ) (2)if Γ < Γ Fermi and P (Γ) ∝ Fermi is the index measured by
Fermi -LAT and σ Γ is the uncertainty on thismeasurement that takes all the systematic and statistical uncertainties into account. redshift z0.1 0.2 0.3 0.4 0.5 0.6 0.7 P r obab ili t y D en s i t y ( A . U . ) Posterior Probability (This work)Danforth (2010)Sanchez et al (2013)
Fig. 5.— Posterior probability density as a function of redshift (red). The blue area rep-resents the redshift range estimated by Danforth et al. (2010) while the green dashed lineindicates the limit of Sanchez et al. (2013).The most probable redshift found with this method is z = 0 . ± .
04, in good agreementwith the independent measure of Danforth et al. (2010), who constrained the distance to bebetween 0 . < z < .
58. Figure 5 gives the posterior probability obtained with the Bayesian 19 –method compared with other measurements of z . Lower and upper limits at a confidencelevel of 95 % can also be derived as 0 . < z < .
56. Note that this method allows thesystematic uncertainties of both instruments (
Fermi -LAT and
H.E.S.S. ) to be taken intoaccount. The spectral index obtained when fitting the H.E.S.S. data with an EBL absorbedPWL using a redshift of 0.49 is compatible with the
Fermi measurement below 80 GeV.
As stated in section 2.1, the
H.E.S.S. data of the flare show a indication of intra-nightvariability, which is used here to test for a possible Lorentz Invariance Violation (LIV). SomeQuantum Gravity (QG) models predict a change of the speed of light at energies close to thePlanck scale ( ∼ GeV). A review of such models can be found in Mattingly (2005) andLiberati (2013). An energy-dependent dispersion in vacuum is searched for in the data bytesting a correlation between arrival times of the photons and their energies. For two pho-tons with arrival times t and t and energies E and E , the dispersion parameter of order n is defined as τ n = t − t E n − E n = ∆ t ∆( E n ) . Here only the linear (n = 1) and quadratic (n = 2) dis-persion 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 forthe expansion of the Universe and an energy E QG at which Quantum Gravity effects areexpected to occur (Jacob & Piran 2008): τ n = ∆ t ∆( E n ) ≃ s ± (1 + n )E n QG κ n (3)where H is the Hubble constant and s ± = − κ n is calculated from the redshift of the source z and the cosmo-logical parameters Ω m , Ω Λ given in the introduction: κ n = Z z (1 + z ′ ) n dz ′ p Ω m (1 + z ′ ) + Ω Λ (4)Using the central value of z = 0 .
49 determined in section 3.2, the distance κ n f or n = 1 and2 is κ = 0 .
541 and κ = 0 . H.E.S.S. flare dataset (MC simulations and original dataset), in order to measure thedispersion and provide 95 % 1-sided lower and upper limits on the dispersion parameter τ n .These limits on τ n will lead to lower limits on E QG using equation 3. 20 – A maximum likelihood method, following Martinez & Errando (2009), was used to cal-culate the dispersion parameter τ n . Albert et al. (2008) applied this method to a flare ofMkn 501, while Abramowski et al. (2011) applied it to a flare of PKS 2155-304. More re-cently, it was used by Vasileiou et al. (2013) to analyse Fermi data of four gamma-ray bursts.The data from Cherenkov telescopes is contaminated by π decay from proton showers,misidentified electrons, or heavy elements such as helium. In the case of PG 1553+113, andcontrary to previous analyses, this background is not negligible: the signal-over-backgroundratio S/B is about 2, compared to 300 for the PKS 2155–304 flare event of July 2006(Aharonian et al. 2007b). The background was included in the formulation of the proba-bility 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 unbinnedlikelihood, function of the dispersion parameter τ n is: L ( τ n ) = n ON Y i =1 P ( E i , t i | τ n ) . (5)The PDF P ( E i , t i | τ n ) associated with each ON event is composed of two terms: P ( E i , t i | τ n ) = w s · P Sig ( E i , t i | τ n ) + (1 − w s ) · P Bkg ( E i , t i ) (6)with P Sig ( E i , t i | τ n ) = 1 N ( τ n ) A eff ( E i , t i ) Λ Sig ( E i ) F Sig ( t i − τ n · E ni ) (7) P Bkg ( E i , t i ) = 1 N ′ A eff ( E i , t i ) Λ Bkg ( E i ) F Bkg (8) w s = n ON − α n OFF n ON . (9)The PDF P Sig includes the emission time distribution of the photons F Sig determined froma parametrization of the observed light curve at low energies (discussed in the next section)and evaluated on t − τ n · E n to take into account the delay due to a possible LIV effect, themeasured signal spectrum Λ Sig and the effective area A eff . The PDF P Bkg is composed of theuniform 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 thebackground 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 correspondsto the relative weight of the signal events in the total ON data set, derived from the numberof events in the ON region n ON and the number of events in the OFF regions n OFF weightedby the inverse number of OFF regions α . More details on the derivation of this function aregiven in Appendix B.1. 21 – The flare data set of the H.E.S.S. analysis (see section 2.1) was used with additionalcuts. To perform the dispersion studies, only uninterrupted data have been kept. Thus, theanalysis was conducted on the first 7 runs, taken during the night of April 26th. Moreover,the cosmic ray flux increases substantially for the 7th run, due to a variation of the zenithangle during this night. This fact, along with its large statistical errors, leads us to discardthis run from the analysis. The 6th run shows little to no variability and was thereforealso removed from the LIV analysis. Since within the ON data set, the signal and thebackground spectra have different indices (Γ
Sig = 4 . Bkg = 2 . E max = 789 GeV was set, corresponding to the last bin with more than 3 σ significancein the reconstructed photon spectrum (see the differential flux during the flare in Fig. 1).A lower cut on the energy at E min = 300 GeV was used in order to avoid large systematiceffects arising from high uncertainties on the H.E.S.S. effective area at lower energies. Theintrinsic light curve of the flare, needed in the formulation of the likelihood, can be obtainedfrom a model of the timed emission or approximated from a subset of the data. To be asmodel-independent as possible, it was here derived from a fit of the measured light curveat low energies (with E < E cut ). The high-energy events (
E > E cut ) were processed in thecalculation 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 onthe energy did not introduce significant effects on the final results. The histogram and thefit (Fig. 6) were obtained as follows: the main idea was to preserve the maximum detectedvariability in the PG 1553+113 flare, together with a significant response in each observedpeak: • The binning was chosen so that at least two adjacent bins of the distribution yield aminimum of 3 σ excess with respect to the average value. • Simple parameterization have been tested on the whole data set (all energies): constant( χ /d.o.f=25/12), single Gaussian ( χ /d.o.f=20/10) and double Gaussian ( χ /d.o.f=8.5/7)functions. The latter is preferred, since it improves the quality of the fit. This shape waschosen to fit the low energy subset of events. Choosing a single Gaussian parametriza-tion would result in a decrease of the sensitivity to time-lag measurements by a factorof two.There is a gap of ∼ ∼
10 min. More importantly,their occurrence is not correlated with the binning: one gap falls in the rising part of the 22 – histONOFFLightcurveAcc
Time (s)0 2000 4000 6000 8000 10000 E xc e ss histONOFFLightcurveAcc Fig. 6.— Time distribution of the excess ON − αOF F in the first 6 runs (70971-70976), withenergies between 300 GeV and 400 GeV. T = 0 corresponds to the time of the first detectedevent in run 70971. The vertical bars correspond to 1 σ statistical errors; the horizontal barscorrespond to the bin width in time. The best fit, in red, was used as the template lightcurve in the maximum likelihood method; the ± σ error envelope is shown in green. 23 –light curve, one is at a maximum, two fall in the decreasing parts and none of the gaps is atthe minimum.Table 6 in Appendix B.2 shows the number of ON and OFF events for the different cutsapplied to the data. τ 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 de-termined from the likelihood curve at the values of τ n where the curve reaches 2.71, whichcorresponds to the 90% C.L. quantile of a χ distribution. However, these CIs are derivedfrom one realization only and do not take into account the “luckiness” factor of this mea-surement. To get statistically significant CIs (“calibrated CIs”), several sets were generatedwith Monte Carlo simulations, with the same statistical significance, light curve model andspectrum as the original data set. No intrinsic dispersion was artificially added. Each simu-lated data set produces a lower limit and an upper limit on τ n . The calibrated lower (upper)limit of the confidence interval is obtained from the mean of the distribution of the per-setindividual lower (upper) limits. Both confidence intervals (from the data only and from thesimulated sets) are listed in Table 7. Sources of systematic errors include uncertainties on thelight curve parameterization, the background contribution, the calculation of the effectivearea, the energy resolution, and the determination of the photon index (see Appendix B.4).The resulting limits on the dispersion τ n using the quadratic sum of the statistical errorsfrom the simulations and the systematic errors determined from data and simulations werecomputed, leading to limits on the energy scale E QG (Eq. 3). The 95 % 1-sided lower limitsfor the subluminal case (s = +1) are: E QG , > . × GeV and E QG , > . × GeVTable 5. Calibrated 95% 1-sided LL and UL (including systematic errors) on thedispersion parameter τ n and derived 95% 1-sided lower limits on E QG .Limits on τ n (s TeV − n ) Lower limits on E QG (GeV)n LL calib+syst U L calib+syst s= −
24 –for linear and quadratic LIV effects, respectively. For the superluminal case (s = –1) thelimits are: E QG , > . × GeV and E QG , > . × GeV for linear and quadraticLIV effects, respectively. Fig. 7 shows a comparison of the different lower limits on E QG , and E QG , for the subluminal case (s = +1) obtained with AGN at different redshifts studiedat very high energies. All these limits, including the present results, have been obtainedunder the assumption that no intrinsic delays between photons of different energies occurat the source. For the linear/subluminal case, the most constraining limit on E QG withtransient astrophysical events has been obtained with GRB 090510: E QG , > . × GeV(Vasileiou et al. 2013). The most constraining limits on E QG with AGN so far have beenobtained by Abramowski et al. (2011) with PKS 2155-304 data observed with H.E.S.S.:E QG , > . × GeV and E QG , > . × GeV for linear and quadratic LIV effects,respectively (95% CL, 1-sided). Compared to the PKS 2155-304 limits, the limits on thelinear dispersion for PG 1553+113 are one order of magnitude less constraining, but thelimits on the quadratic dispersion are of the same order of magnitude since the source islocated at a higher redshift. This highlights the interest in studying distant AGN, in spiteof the difficulties due to limited photon statistics.
4. Conclusions
A VHE γ -ray flaring event of PG 1553+113 has been detected with the H.E.S.S. telescopes, with a flux increasing of a factor of 3. No variability of the spectral index hasbeen found in the data set, but indication of intra-night flux variability is reported in thiswork. In HE γ rays, no counterpart of this event can be identified, which may be interpretedas the sign of injection of high energy particles emitting predominantly in VHE γ rays. Suchparticles might not be numerous enough to have a significant impact on the HE flux duringeither their acceleration or cooling phases.The data were used to constrain the redshift of the source using a new approach based onthe absorption properties of the EBL imprinted in the spectrum of a distant source. Takinginto account all the instrumental systematic uncertainties, the redshift of PG 1553+113 isdetermined as being z = 0 . ± . H.E.S.S. data of a flare of PG 1553+113.This analysis relies on the indication of the intra-night variability of the flare at VHE. Nosignificant dispersion was measured, and limits on the E QG scale were derived, in a region of 25 –Fig. 7.— Lower limits on E QG , from linear dispersion (left) and on E QG , from quadraticdispersion (right) for the subluminal case (s = +1) obtained with AGN as a function ofredshift. The limits are given in terms of E Planck . The constraints from Mkn 421 have beenobtained by Biller et al. (1999), from Mkn 501 by Albert et al. (2008), and from PKS 2155-304 by Abramowski et al. (2011). 26 –redshift unexplored until now. Limits on the energy scale at which QG effects causing LIVmay arise, derived in this work, are E QG , > . × GeV and E QG , > . × GeV forthe subluminal case. Compared with previous limits obtained with the PKS 2155-304 flareof 2006 July, the limits for PG 1553+113 for a linear dispersion are one order of magnitudeless constraining while limits for a quadratic dispersion are of the same order of magnitude.With the new telescope placed at the center of the
H.E.S.S. array that provides an energythreshold of several tens of GeV, a better picture of the variability patterns of AGN flaresshould be obtained. The future Cherenkov Telescope Array (CTA) will increase the numberof flare detections (Sol et al. 2013) with better sensitivity, allowing for the extraction of evenmore constraining limits on the LIV effects.
Acknowledgements
The support of the Namibian authorities and of the University of Namibia in facilitatingthe construction and operation of H.E.S.S. is gratefully acknowledged, as is the support bythe German Ministry for Education and Research (BMBF), the Max Planck Society, theFrench Ministry for Research, the CNRS-IN2P3 and the Astroparticle Interdisciplinary Pro-gramme of the CNRS, the U.K. Particle Physics and Astronomy Research Council (PPARC),the IPNP of the Charles University, the South African Department of Science and Technol-ogy and National Research Foundation, and by the University of Namibia. We appreciatethe excellent work of the technical support staff in Berlin, Durham, Hamburg, Heidelberg,Palaiseau, Paris, Saclay, and in Namibia in the construction and operation of the equipment.The
Fermi
LAT Collaboration acknowledges generous ongoing support from a numberof agencies and institutes that have supported both the development and the operation ofthe LAT as well as scientific data analysis. These include the National Aeronautics andSpace Administration and the Department of Energy in the United States, the Commis-sariat l’Energie Atomique and the Centre National de la Recherche Scientifique / InstitutNational de Physique Nuclaire et de Physique des Particules in France, the Agenzia SpazialeItaliana and the Istituto Nazionale di Fisica Nucleare in Italy, the Ministry of Education,Culture, Sports, Science and Technology (MEXT), High Energy Accelerator Research Or-ganization (KEK) and Japan Aerospace Exploration Agency (JAXA) in Japan, and the K.A. Wallenberg Foundation, the Swedish Research Council and the Swedish National SpaceBoard in Sweden.Additional support for science analysis during the operations phase is gratefully acknowl-edged from the Istituto Nazionale di Astrofisica in Italy and the Centre National d’EtudesSpatiales in France. 27 –DS work is supported by the LABEX grant enigmass. The authors want to thanks F.Krauss for her useful comments. 28 –
A. Bayesian model used to constrain the redshift
A Bayesian approach has been used to compute the redshift value of PG 1553+113in Section 3.2. The advantage of such a model is that systematic uncertainties, whichare important in Cherenkov astronomy, can easily be included in the calculation. In thefollowing, the notation Θ for the model parameters and Y for the data set is adopted.All normalization constants are dropped in the development of the model, and the finalprobability is normalized at the end.Bayes’ Theorem, based on the conditional probability rule, allows us to write the poste-rior probability P (Θ | Y ) for the model parameters Θ as the product of the likelihood P ( Y | Θ)and the prior probability P (Θ): P (Θ | Y ) ∝ P (Θ) P ( Y | Θ) . The likelihood is the quantity that is maximized during determination of the best-fitspectrum (Piron et al. 2001). It is at this step that the
H.E.S.S. data, taken during theflare, were actually used. The spectrum model here is a simple power law corrected for theEBL absorption: φ = N × ( E/E ) − Γ × e − τ ( E,z ) . The model parameters are then N , Γ and z .The prior is the most difficult and most interesting part of the model. To derive it, N and Γ are assumed to be independent from each other and independent of the redshift.In contrast, the prior on the redshift might depend on N and Γ. Then, the prior can besimplified using the conditional probability rule: P (Θ) = P ( z | N, Γ) P ( N ) P (Γ)As much as possible, weak assumptions should be made to write a robust prior thenoften flat (i.e. P ∝ const) are used. Priors should also be based on a physical meaningand not contradict the physical and observed properties of the objects. For the purpose ofthis model, the prior on N is assumed to be flat and the prior on the spectral index is atruncated Gaussian P (Γ) ∝ N G (Γ , Γ Fermi , σ Γ ) if Γ < Γ Fermi and P (Γ) = ∝ const otherwise.The values of Γ Fermi and σ Γ are obtained by analyzing LAT data below 80 GeV (see section3 and Table 3). Here, it is assumed that the intrinsic spectrum in the VHE range cannotbe harder than the Fermi -LAT measurement. σ Γ takes into account the statistical andsystematic uncertainties on the Fermi -LAT measurement and also the systematic uncertainty 29 –on the
H.E.S.S. spectrum ( σ = 0 .
20, see Aharonian et al. 2006b) added quadratically and σ Γ = 0 .
33 for a mean value of Γ
Fermi = 1 . z is much more difficult to determine. A flat prior has no physical motiva-tions since the probability to detect sources at TeV energy decreases with the redshift. Thenumber of sources detected at TeV energy is not sufficient to use the corresponding redshiftdistribution as a prior.A prior which takes into account the EBL, can be derived assuming a population ofsources with a constant spatial density. In the small space element 4 πz dz , the numberof such sources scales ∝ z . For any given luminosity, their flux (which scales with theprobability to detect them) is scaled by z − exp ( − τ ( z )). Lacking a proper knowledge ofthe intrinsic luminosity function of VHE γ -ray blazars, a reasonable assumption on thedetection probability of a blazar at any redshift is a scaling proportional to the flux fora given luminosity, i.e., ∝ z − exp ( − τ ( z )). Putting everything together, the prior on theredshift reads P ( z | N, Γ) = P ( z ) ∝ exp ( − τ ( z ))Finally, the prior we use for our analysis is: P (Θ) ∝ exp ( − τ ( z )) N G (Γ , . , .
33) (A1)if Γ < .
72 and P (Θ) ∝ exp ( − τ ( z ))otherwise. Putting all the components of the model together and marginalizing over thenuisance parameters N and Γ, the probability on the redshift can be computed numerically.The obtained mean value is z = 0 . ± .
04. At a confidence level of 95 %, the redshift isbetween 0 . < z < . B. Development of the LIV methodB.1. Modified maximum likelihood method
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 thesignal-over-background ratio is about 2. The background has been included in the formula-tion of the probability density function (PDF): in the most general case, for given numbersof signal and background events s and b in the observation region (“ON” region), for a givendispersion parameter τ n , the unbinned likelihood is: L ( n ON , n OFF | s, b, τ n ) = Pois( n ON | s + b ) · Pois (cid:18) n OFF | bα (cid:19) · n ON Y i =1 P ( E i , t i | s, b, τ n ) (B1)The PDF P ( E i , t i | s, b, τ n ) associated with each gamma-like particle characterized by its time t i and energy E i contains two terms (signal and background): P ( E i , t i | s, b, τ n ) = w s · P Sig ( E i , t i | τ n ) + (1 − w s ) · P Bkg ( E i , t i ) (B2)with w s = ss + b . (B3) n ON is the number of events detected in the source ON region included in the fit range[ E cut ; E max ] × [ t min ; t 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. Pois( n ON | s + b ) (Pois( n OFF | b/α ))is the Poisson distribution with index n ON ( n OFF ) and parameter s + b ( b/α ). The likeli-hood 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 Eq. B2 are equal to 1.The probabilities P Sig and P Bkg are defined as: P Sig ( E i , t i | τ n ) = 1 N ( τ n ) · R Sig ( E i , t i | τ n ) (B4) P Bkg ( E i , t i ) = 1 N ′ · R Bkg ( E i , t i ) (B5)with R Sig ( E, t | τ n ) = Z ∞ E true =0 D ( E, E true ) A eff ( E true , t ) Λ Sig ( E true ) F Sig ( t − τ n · E n true ) dE true (B6) R Bkg ( E, t ) = Z ∞ E true =0 D ( E, E true ) A eff ( E true , t ) Λ Bkg ( E true ) F Bkg ( t ) dE true . (B7) 31 – P Sig ( E i , t i | τ n ) is the probability that the event ( E i , t i ) is a photon emitted at the source anddetected on Earth with a delay τ n E n . It takes into account the emission (time distribution F Sig ( t ) and energy spectrum Λ Sig ( E ) at the source), the propagation (delay τ n · E ni due topossible LIV effect) and the detection of a photon by the detector (H.E.S.S. energy resolution D ( E, E true ) and effective area A eff ( E, t )). P 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 uniformtime distribution: F Bkg ( t ) = F Bkg . The background energy distribution Λ
Bkg is measured fromOFF regions. N ( τ n ) (resp. N ′ ) is the normalization factor of the PDF P Sig (resp. P Bkg ) inthe range [ E cut ; E max ] × [ t min ; t 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 Sig ( E i , t i | τ n ) and P Bkg ( E i , t i ): P Sig ( E i , t i | τ n ) = 1 N ( τ n ) · A eff ( E i , t i )Λ Sig ( E i ) F Sig ( t i − τ n · E ni ) (B8) P Bkg ( E i , t i ) = 1 N ′ · A eff ( E i , t i )Λ Bkg ( E i ) F Bkg (B9)The best estimate of the dispersion parameter b τ n is obtained by maximizing the likelihood L ( τ n ). B.2. Selection cuts
Table 6 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 actual energy resolution is of the order of 10 % in this range.
Table 6. Selections applied to the ON and OFF data setsSelection 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 32 – B.3. Test of the method, confidence intervals
The method has been tested on Monte Carlo (MC) simulated sets. Each set was com-posed of n ON = 72 ON events, as in the real data sample: • s = 50 signal events with times following the template light curve (Fig. 6) shifted bya factor τ n, inj · E i ; energies follow a power law spectrum of photon index Γ Sig = 4 . • b = 22 background events with times following a uniform distribution and energiesdrawn from a power law spectrum of index Γ Bkg = 2 .
5, degraded by the acceptanceand convoluted by the energy resolution.For a given injected dispersion, the maximum likelihood method is applied to each MC-simulatedset. The initial light curve and energy spectrum were used as templates in the model insteadof fitting them for each set.Figure 8 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 thedistribution of the best estimates ˆ τ . The blue line shows the result of a linear fit. The sloperoughly corresponds to the percentage of signal in the total ON data set. It is due to theloss of sensitivity resulting from the part of the data sets with no dispersion. A systematicshift is observed of about 100 s TeV − , well bellow 1 σ value – the RMS of the best estimatedistribution is of 361 s TeV − . The results in this paper have not been corrected for thisbias.The coverage is not necessarily proper, i.e. the number of sets for which the injecteddispersion value τ inj lies between the set’s lower limit (LL) and upper limit (UL) does notmatch the required 95 % 1-sided confidence level. The common cut used on the likelihoodcurves to get the LLs/ULs has been iteratively adjusted to ensure a correct statistical cov-erage: using this new cut, 95 % of the realizations provide CIs that include the injecteddispersion τ n, inj . The initial coverage was about 85 % for a cut on 2 ln L of 2.71. The newcommon cut, found iteratively at 3.5, ensures the desired 90 % 2-sided CL (approx. 95 %1-sided CL). Figure 9 shows the distributions of the best estimates, the 95% 1-sided LLs andULs for τ , inj = 0 s TeV − (linear case) and τ , inj = 0 s TeV − (quadratic case); the means ofthe lower and upper limit distributions, shown as a blue vertical line, are used to constructthe “calibrated confidence interval”.To get CIs from data, a maximum likelihood method is applied to the original dataset and gives a best estimate τ databest . The cut value determined from the simulations to 33 – ) -1 (s TeV τ Injected -2000 -1500 -1000 -500 0 500 1000 1500 2000 ) - ( s T e V τ R e c o n s t r u c t e d -2000-1500-1000-5000500100015002000 Fig. 8.— Means of the reconstructed dispersion versus the real (injected dispersion) for thelinear case n = 1; for a given injected dispersion, errors bars correspond to the means of thedistribution of the upper and lower limits (90 % 2-sided ≃
95 % 1-sided). The blue line is alinear fit to the points. The red line shows the ideally obtained curve τ recontructed = τ injected obtained in the case S/B = ∞ .ensure proper coverage is applied on the original data set to obtain LL data and U L data . The“calibrated” limits LL calib and U L calib , combining τ databest from data together with MC results,are taken as LL calib = τ databest − | τ MCbest − LL MC | (B10) U L calib = τ databest + | τ MCbest − U L MC | with τ MCbest , LL MC and UL MC defined as the mean of the per-set best-estimate distribution,LL distribution, and UL distribution respectively.Table 7 lists the CIs determined in both ways, i.e., data-only and calibrated ones: LL data n and LL calib n (resp. U L data n and U L calib n ) are compatible within 10 %. In this work, calibratedCIs have been used to derive the final lower limits on E QG . They are preferred over data-onlyCIs as they provides statistically well defined confidence levels. They also ensure coherentcomparison with previous published results, e.g. with PKS 2155–304 by Abramowski et al.(2011) and GRB studies by Vasileiou et al. (2013). 34 – B.4. Estimation of the systematics
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 lightcurve parametrization. Other sources of systematic errors include the contribution of thebackground, effect of the change of photon index, the energy resolution and the effective areadetermination of the detector. To study the following four contributions, new simulated datasets have been built, each one with different input parameters: • background contribution: photons and background events have been reallocated withinthe ON data set in the fit range [ E cut ; E max ], introducing a 1 σ fluctuation in the numberof signal event s in the ON data set; • 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 withenergy); • energy resolution: reconstructed energies have been replaced by the true energies; thiscorresponds to a shift of about 10 % on the reconstructed energy values; • photon index: changed by one standard deviation ( ± . σ andthe lower 1 σ contours of the template, shown in Fig. 6. The change in mean lower and upperlimits on the dispersion parameter τ n gives an estimate of the systematic error associatedto each contribution . An additional systematic contribution comes from the shift arisingfrom the method found with simulation (see Appendix B.3). Table 8 summarizes all studiedsystematic contributions. The overall estimated systematic error on τ n is 330 s TeV − forthe linear case (n = 1) and 555 s TeV − for the quadratic case (n = 2); they were includedin the calculation of the limits on E QG by adding the statistical and the systematic errors inquadrature. In particular the errors on the peak positions constitute the most important part of the uncertaintyon the template light curve contributing to the likelihood fit – see previous works, e.g. Abramowski et al.(2011). Therefore, the covariance matrix of the fit of the template was studied in detail ; the peak positionswere varied by values of ± σ extracted from the covariance matrix. This study led to an increase in overallsystematics of the order of 20% for τ and 40% for τ , and a decrease of maximum 7% and 2% of limits onE QG , and E QG , respectively.
35 –Table 7. Linear (top) and quadratic (bottom) dispersion parameter; from left to right:best estimate, LL and UL from data (cut on likelihood curve), LL and UL from MCsimulations (means of per-set LL and UL distributions), calibrated LL and UL(combination of data and MC), calibrated LL and UL including systematic errors.Dispersion parameters τ n,best , LLs and ULs are in s TeV − n .n τ data n,best LL data n U L data n τ MC n,best LL MC n U L MC n LL calib n U L calib n LL calib n U L calib n with systematics1 -131.7 -806.7 554.7 99.1 -526.3 725.6 -757.1 494.8 -838.9 576.42 -287.5 -1449.9 853.6 217.2 -942.0 1395.0 -1446.7 890.3 -1570.5 1012.4Table 8. Summary of all studied systematic contributions. The main systematic errors aredue to the uncertainties on the light curve parametrization.Estimated error τ τ on input parameters (s TeV − ) (s TeV − )Background contribution < < < < < < < < < < ∼ ∼ pP i syst i < <
555 36 – hTauEntries 1000Mean 99.1RMS 361.1Underflow 0Overflow 0 / ndf 2 χ ± ± ± (s/TeV) τ -1000 0 1000 N u m be r o f s e t s hTauEntries 1000Mean 99.1RMS 361.1Underflow 0Overflow 0 / ndf 2 χ ± ± ± χ ± ± -524.3 Sigma 7.7 ± (s/TeV) τ -2000 -1000 0 1000 N u m be r o f s e t s hTauLL90Entries 1000Mean -526.3RMS 362.9Underflow 0Overflow 0 / ndf 2 χ ± ± -524.3 Sigma 7.7 ± χ ± ±
731 Sigma 8.4 ± (s/TeV) τ -1000 0 1000 2000 N u m be r o f s e t s hTauUL90Entries 1000Mean 725.6RMS 367.5Underflow 0Overflow 0 / ndf 2 χ ± ±
731 Sigma 8.4 ± χ ± ± ± ) (s/TeV τ -2000 0 2000 N u m be r o f s e t s hTauEntries 1000Mean 217.2RMS 676.1Underflow 0Overflow 0 / ndf 2 χ ± ± ±
664 hTauLL90Entries 1000Mean -942RMS 673.4Underflow 0Overflow 0 / ndf 2 χ ± ± -939.1 Sigma 16.1 ± ) (s/TeV τ -2000 0 2000 N u m be r o f s e t s hTauLL90Entries 1000Mean -942RMS 673.4Underflow 0Overflow 0 / ndf 2 χ ± ± -939.1 Sigma 16.1 ± χ ± ± ± ) (s/TeV τ N u m be r o f s e t s hTauUL90Entries 1000Mean 1395RMS 704.6Underflow 0Overflow 1 / ndf 2 χ ± ± ± Fig. 9.— Distributions of the best estimates, the 95% 1-sided lower and upper limits fromsimulations in case of no injected dispersion ( τ n, inj = 0 s TeV − n ), for n = 1 (top) and n = 2(bottom); dispersion values are in s TeV − n . The blue vertical line on the LL (resp. UL)distribution shows LL MC (resp. U L MC ), defined as the mean of the distribution. 37 – REFERENCES
Abdo, A. A., et al. 2009a, ApJ, 707, 1310—. 2009b, ApJS, 183, 46—. 2010a, ApJ, 708, 1310—. 2010b, ApJ, 722, 520Abramowski, A., et al. 2011, Astroparticle Physics, 34, 738Abramowski, A., et al. 2013a, A&A, 552, A118—. 2013b, A&A, 550, A4—. 2014, Ap Librae paper Submitted to A&AAckermann, M., et al. 2012a, ApJS, 203, 4—. 2012b, Science, 338, 1190Aharonian, F., et al. 2006a, A&A, 448, L19—. 2006b, A&A, 457, 899—. 2007a, ApJ, 664, L71—. 2007b, ApJ, 664, L71—. 2008, A&A, 477, 481Ajello, M., et al. 2014, ApJ, 780, 73Albert, J., et al. 2008, Physics Letters B, 668, 253Aleksi´c, J., et al. 2012, ApJ, 748, 46—. 2014, ArXiv e-printsAtwood, W. B., et al. 2009, ApJ, 697, 1071Becherini, Y., Djannati-Ata¨ı, A., Marandon, V., Punch, M., & Pita, S. 2011, AstroparticlePhysics, 34, 858Berge, D., Funk, S., & Hinton, J. 2007, A&A, 466, 1219 38 –Biller, S. D., et al. 1999, Physical Review Letters, 83, 2108Bregeon, J., Charles, E., & M. Wood for the Fermi-LAT collaboration. 2013, ArXiv 1304.5456Cortina, J. 2012a, The Astronomer’s Telegram, 3977, 1—. 2012b, The Astronomer’s Telegram, 4069, 1Danforth, C. W., Keeney, B. A., Stocke, J. T., Shull, J. M., & Yao, Y. 2010, The Astrophys-ical Journal, 720, 976de Naurois, M., & Rolland, L. 2009, Astroparticle Physics, 32, 231Franceschini, A., Rodighiero, G., & Vaccari, M. 2008, A&A, 487, 837Gaidos, J. A., et al. 1996, Nature, 383, 319Giommi, P., Ansari, S. G., & Micol, A. 1995, A&AS, 109, 267Green, R. F., Schmidt, M., & Liebert, J. 1986, ApJS, 61, 305Hauser, M. G., & Dwek, E. 2001, ARA&A, 39, 249Hinton, J. A., & the H.E.S.S. Collaboration. 2004, New A Rev., 48, 331Jacob, U., & Piran, T. 2008, Journal of Cosmology and Astroparticle Physics, 2008, 031Komatsu, E., et al. 2011, The Astrophysical Journal Supplement Series, 192, 18Li, T.-P., & Ma, Y.-Q. 1983, ApJ, 272, 317Liberati, S. 2013, Classical and Quantum Gravity, 30, 133001Martinez, M., & Errando, M. 2009, Astroparticle Physics, 31, 226Mattingly, D. 2005, Living Reviews in Relativity, 8Mattox, J. R., et al. 1996, ApJ, 461, 396Osterman, M. A., et al. 2006, AJ, 132, 873Piron, F., et al. 2001, A&A, 374, 895Prandini, E., Dorner, D., Mankuzhiyil, N., Mariotti, M., Mazin, D., & for the MAGICCollaboration. 2009, ArXiv 0907.0157Rector, T. A., Gabuzda, D. C., & Stocke, J. T. 2003, AJ, 125, 1060 39 –Reimer, A., Costamante, L., Madejski, G., Reimer, O., & Dorner, D. 2008, ApJ, 682, 775Sanchez, D. A., & Deil, C. 2013, ArXiv 1307.4534Sanchez, D. A., Fegan, S., & Giebels, B. 2013, A&A, 554, A75Scargle, J. D. 1998, ApJ, 504, 405Sol, H., et al. 2013, Astroparticle Physics, 43, 215Treves, A., Falomo, R., & Uslenghi, M. 2007, A&A, 473, L17Urry, C. M., & Padovani, P. 1995, PASP, 107, 803Vasileiou, V., et al. 2013, Physical Review D, 87, 122001Vaughan, S., Edelson, R., Warwick, R. S., & Uttley, P. 2003, MNRAS, 345, 1271