Measurement of the EBL spectral energy distribution using the VHE gamma-ray spectra of H.E.S.S. blazars

DOI: 10.1051/0004-6361/201731200

© E S O 2017

Astronomy

&

Astrophysics

Measurement of the EBL spectral energy distribution using the VHE y-ray spectra of H.E.S.S. blazars

H.E.S.S. Collaboration, H. Abdalla1, A. Abramowski2, F. Aharonian3,4’5, F. Ait Benkhali3, A. G. Akhperjanian6,5,t , T. Andersson10, E. O. Anguner21, M. Arakawa41, M. Arrieta15, P. Aubert23, M. Backes8, A. Balzer9, M. Barnard1, Y. Becherini10,

J. Becker Tjus11, D. Berge12, S. Bernhard13, K. Bernlohr3, R. Blackwell14, M. Bottcher1, C. Boisson15, J. Bolmont16, S. Bonnefoy35, P. Bordas3, J. Bregeon17, F. Brun25, P. Brun18* , M. Bryan9, M. Buchele34, T. Bulik19, M. Capasso27, J. Carr20, S. Casanova21,3, M. Cerruti16, N. Chakraborty3, R. C. G. Chaves17, A. Chen22, J. Chevalier23, M. Coffaro27, S. Colafrancesco22, G. Cologna24, B. Condon25, J. Conrad26, Y. Cui27, I. D. Davids1,8, J. Decock18, B. Degrange28, C. Deil3, J. Devin17, P. de Wilt14, L. Dirson2, A. Djannati-Atai29, W. Domainko3, A. Donath3, L. O’C. Drury4, K. Dutson31, J. Dyks32, T. Edwards3, K. Egberts33, P. Eger3, J.-P. Ernenwein20, S. Eschbach34, C. Farnier26,10, S. Fegan28, M. V. Fernandes2, A. Fiasson23, G. Fontaine28, A. Forster3,

S. Funk34, M. FuBling35, S. Gabici29, Y. A. Gallant17, T. Garrigoux1, G. Giavitto35, B. Giebels28, J. F. Glicenstein18, D. Gottschall27, A. Goyal36, M.-H. Grondin25, J. Hahn3, M. Haupt35, J. Hawkes14, G. Heinzelmann2, G. Henri30, G. Hermann3,

J. A. Hinton3, W. Hofmann3, C. Hoischen33, T. L. Holch7, M. Holler13, D. Horns2, A. Ivascenko1, H. Iwasaki41, A. Jacholkowska16, M. Jamrozy36, M. Janiak32, D. Jankowsky34, F. Jankowsky24, M. Jingo22, T. Jogler34, L. Jouvin29, I. Jung-Richardt34, M. A. Kastendieck2, K. Katarzyiiski37, M. Katsuragawa42, U. Katz34, D. Kerszberg16, D. Khangulyan41,

B. Khćlifi29, J. King3, S. Klepser35, D. Klochkov27, W. Kluzniak32, D. Kolitzus13, Nu. Komin22, K. Kosack18, S. Krakau11, M. Kraus34, P. P. Kruger1, H. Laffon25, G. Lamanna23, J. Lau14, J.-P. Lees23, J. Lefaucheur15, V. Lefranc18, A. Lemiere29,

M. Lemoine-Goumard25, J.-P. Lenain16, E. Leser33, T. Lohse7, M. Lorentz18 *, R. Liu3, R. López-Coto3, I. Lypova35, V. Marandon3, A. Marcowith17, C. Mariaud28, R. Marx3, G. Maurin23, N. Maxted14,43, M. Mayer7, P. J. Meintjes38, M. Meyer26,

A. M. W. Mitchell3, R. Moderski32, M. Mohamed24, L. Mohrmann34, K. Mora26, E. Moulin18, T. Murach35, S. Nakashima42, M. de Naurois28, F. Niederwanger13, J. Niemiec21, L. Oakes7, P. O’Brien31, H. Odaka42, S. Ohm35, M. Ostrowski36, I. Oya35, M. Padovani17, M. Panter3, R. D. Parsons3, N. W. Pekeur1, G. Pelletier30, C. Perennes16, P.-O. Petrucci30, B. Peyaud18, Q. Piel23,

S. Pita29, H. Poon3, D. Prokhorov10, H. Prokoph12, G. Puhlhofer27, M. Punch29,10, A. Quirrenbach24, S. Raab34, R. Rauth13, A. Reimer13, O. Reimer13, M. Renaud17, R. de los Reyes3, S. Richter1, F. Rieger3,39, C. Romoli4, G. Rowell14, B. Rudak32, C. B. Rulten15, V. Sahakian6,5, S. Saito41, D. Salek40, D. A. Sanchez23 *, A. Santangelo27, M. Sasaki34, R. Schlickeiser11, F. Schussler18, A. Schulz35, U. Schwanke7, S. Schwemmer24, M. Seglar-Arroyo18, M. Settimo16, A. S. Seyffert1, N. Shafi22, I. Shilon34, R. Simoni9, H. Sol15, F. Spanier1, G. Spengler26, F. Spies2, Ł. Stawarz36, R. Steenkamp8, C. Stegmann33,35, K. Stycz35,

I. Sushch1, T. Takahashi42, J.-P. Tavernet16, T. Tavernier29, A. M. Taylor4, R. Terrier29, L. Tibaldo3, D. Tiziani34, M. Tluczykont2, C. Trichard20, N. Tsuji41, R. Tuffs3, Y. Uchiyama41, D. J. van der Walt1, C. van Eldik34, C. van Rensburg1, B. van Soelen38, G. Vasileiadis17, J. Veh34, C. Venter1, A. Viana3, P. Vincent16, J. Vink9, F. Voisin14, H. J. Volk3, T. Vuillaume23, Z. Wadiasingh1,

S. J. Wagner24, P. Wagner7, R. M. Wagner26, R. White3, A. Wierzcholska21, P. Willmann34, A. Wornlein34, D. Wouters18, R. Yang3, D. Zaborov28, M. Zacharias1, R. Zanin3, A. A. Zdziarski32, A. Zech15, F. Zefi28, A. Ziegler34, and N. Zywucka36

(Affiliations can be fo u n d after the references) Received 19 M ay 2017 / Accepted 13 July 2017

ABSTRACT

Very high-energy y rays (VHE, E > 100 GeV) propagating over cosm ological distances can interact w ith the low-energy photons o f the extra- galactic background light (EBL) and produce electron-positron pairs. The transparency of the Universe to VHE y rays is then directly related to the spectral energy distribution (SED) of the EBL. The observation o f features in the VHE energy spectra o f extragalactic sources allows the EBL to be measured, w hich otherwise is very difficult. An EBL m odel-independent m easurem ent o f the EBL SED w ith the H.E.S.S. array of Cherenkov telescopes is presented. It was obtained by extracting the EBL absorption signal from the reanalysis o f high-quality spectra o f blazars.

From H.E.S.S. data alone the EBL signature is detected at a significance o f 9.5ix, and the intensity o f the EBL obtained in different spectral bands is presented together with the associated y-ray horizon.

Key words. gam m a rays: galaxies - BL Lacertae objects: general - cosm ic background radiation - infrared: diffuse background

1. Introduction

The extragalactic background light (EBL) is the second most intense background photon field in the Universe after the cos­

mic microwave background. This diffuse radiation stems from the integrated light emitted through thermal and non-thermal processes, and its reprocessing by the interstellar medium

0 Deceased.

* Corresponding authors: H.E.S.S. Collaboration, e-mail: c o n t a c t.h e s s @ h e s s - e x p e r i m e n t.e u

over cosmic history. It covers wavelengths ranging from the ultraviolet to far-infrared and submillimeter wavelengths. Di­

rect measurements of the EBL are very difficult because of fore­

ground contamination due to zodiacal light and diffuse Galactic light (Hauser et al. 1998). Lower limits have been derived from galaxy counts, and models have been developed to describe its spectral energy distribution (SED; see, e.g., Franceschini et al.

2008; Dominguez et al. 2011; Kneiske & Dole 2010; Finke et al.

2010; Gilmore et al. 2012). This SED is usually described with

two main components: an optical component due to starlight

em ission and an infrared com ponent due to the reprocessing o f starlight by dust. T he E B L SED contains unique inform ation about galaxy form ation and evolution; its study is therefore o f interest for cosm ology. This low -energy p h oton b ackground is responsible for th e lim ited horizon o f very high-energy (V H E , E > 100 G eV ) photons, since these y rays in teract w ith E B L photons through th e production o f electron-positron p airs, re ­ sulting in attenuated observed fluxes above the reaction th resh ­ old (N ikishov 1962; G ould & S chreder 1967) . W hile this affects the study o f extragalactic y -ra y sources, it also provides a w ay to p ro b e the E B L itself (for review s see, e.g., H auser & D w ek 2 0 0 1 ; D w ek & K rennrich 2 0 1 3 ; C ostam ante 2013) .

T he attenuation o f y rays on the E B L is an energy-dependent process w hich leads to a specific spectral signature. O bserva­

tions o f features in the V H E spectra o f extragalactic sources w ith C herenkov telescopes like H .E .S.S. can thus b e u sed to co n ­ strain the E B L u nder som e assum ptions on the intrinsic spec­

tra o f the co n sidered sources. Indeed, a m ajo r com plication in constraining the E B L w ith y rays com es from the in determ i­

nacy o f the intrinsic energy spectra o f sources, and consequently there is a possible degeneracy betw een intrinsic curvature and E B L attenuation. T he technique th at w as first applied to co n ­ strain the E B L w ith H .E.S.S. relied on the assum ption o f a th e ­ oretical lim it for th e hardness o f the intrinsic pow er-law spectra.

U pper lim its on the level o f E B L w ere obtained given the so ft­

ness o f th e observed spectra. U sing the tw o blazars H 2 3 5 6 -3 0 9 (z = 0.165) and 1ES 1 1 0 1 -2 3 2 (z = 0.186), H .E.S.S. show ed that th e U niverse was m o re transparent to y rays than expected at th at tim e from d irect E B L m easurem ents (A haronian et al.

2 0 0 6 a; M atsum oto e t al. 2005) . T he upper lim its on the E B L density o btained w ith H .E .S.S. turned out to b e close to the low er lim its derived from galaxy counts. A global reassessm ent o f E B L m odels follow ed. This H .E.S.S. study was follow ed by a m odel- d ependent determ ination o f the E B L (A bram ow ski et al. 2 013c) obtained by sim ultaneously fitting th e E B L optical depths and intrinsic spectra o f a sam ple o f extragalactic sources w ith a m axim um -likelihood m eth o d assum ing sm ooth concave in trin ­ sic shapes. T he shape o f the E B L SED w as frozen to the shape o f the m odel given in F ranceschini e t al. (2008) . O nly the n o r­

m alization to this m o d el was let free. T he overall test statistic led to an 8 .8 ^ detection o f E B L absorption w ith resp ect to no absorption, w ith a n orm alization factor relative to this m odel o f 1.27+018 (stat) ± 0.25(sys) in the 1.2 jum to 5.5 jum w avelength range. In th e follow ing, this study is referred to as H ESS2013.

O ther E B L constraints an d m easurem ents using y rays have been conducted as w ell, for exam ple w ith Ferm i l A t , M A G IC , and VERITAS (A ckerm ann et al. 2 0 1 2 ; A beysekara et al. 2 0 1 5 ; A hnen et al. 2016) .

T he new analysis p resen ted in th e p rese n t p aper follow s a different approach, focusing on the determ ination o f the shape o f th e E B L SED in addition to its overall norm alization. A n ex ­ tended sam ple o f blazars is used w ith resp ect to H E SS 2013, si­

m ultaneously fitting th eir V H E intrinsic spectra together w ith a generic attenuation. As the E B L is expected to leave a ty p i­

cal energy-dependent and redsh ift-d ep en d en t im print on th e o b ­ served spectra, the detection o f such a m odulation can b e used to translate the absorption pattern into spectrally reso lv ed E B L intensity levels. This analysis considers high-quality V H E spec­

tra and assum es featureless intrinsic spectra, allow ing for in trin ­ sic curvature. This approach aim s n o t only for a H .E.S.S. m e a­

surem ent o f the E B L SED in dependent o f any E B L m odel, but also fo r a generic characterization o f the U n iv erse’s tran s­

parency to V H E y rays w ith the few est possible priors. B e­

yond the interest in the E B L p er se, this is p articularly relevant

for th e study o f p o tential second-order processes in the p ro p ­ agation o f y rays over cosm ological distances. T hese include conversion into axion-like particles (e.g., S anchez-C onde e t al.

2 0 0 9 ; A bram ow ski e t al. 2 013a), L orentz invariance violation (e.g., S tecker & G lashow 2 0 0 1 ; Jacob & P iran 2008), o r cascade em ission in extragalactic m agnetic fields (e.g., A haronian e t al.

1994; Taylor e t a l. 2011) .

S pectral features are searched fo r in the reanalysis o f H .E.S.S. phase-I (four-telescope) data w ith a new m eth o d used to m easure th e E B L SED . U sing only H .E .S.S. data offers the possibility o f handling system atic uncertainties from different spectra in a hom ogeneous and w ell-controlled way. F urther­

m ore, pub lish ed spectral points are n o t usually released together w ith th eir covariance m atrix. U sing this additional inform a­

tion, these results are expected to b e m ore ro b u st than sim ilar studies only using publish ed spectra from different C herenkov telescopes (e.g., O r r e t a l . 2 0 1 1 ; M eyer e t a l. 2 0 1 2 ; S in h a e ta l.

2 0 1 4 ; B iteau & W illiam s 2015) .

This p aper is organized as follow s. T he b lazar sam ple and the data analysis are described in Sect. 2 . In Sect. 3 the need for an energy- and redsh ift-d ep en d en t m odulation in the energy spectra is dem onstrated. In Sect. 4 the E B L absorption process is p resented in detail, and in Sect. 5 the m ethod u sed to translate the m odulation seen in spectra in term s o f E B L is described. The results are p resented and discussed in Sect. 6 .

2. H.E.S.S. data analysis 2.1. D ata reduction

H .E.S.S. is an array o f five im aging atm ospheric C herenkov te le­

scopes located in th e K hom as H ighland, N am ib ia (23° 1 6 '1 8 "

S, 1 6 °3 0 '0 1 " E ), a t an elevation o f 1800 m above sea level. In this w ork, only data from th e four telescopes o f the first phase o f H .E .S.S. are used. This initial four-telescope array detects y rays above ~ 1 0 0 G eV w ith an energy resolution b etter than 15% (A haronian e t al. 2006b) . D ata reduction is p erform ed u s­

ing th e M o d el A n a ly sis technique (de N aurois & R olland 2009) in w hich reco rd ed air-show er im ages are com pared to tem plate im ages p re-calculated using a sem i-analytic m o d el and a log- likelihood optim ization technique. F o r a w ider energy coverage, the loose cuts o f the M o d el A n a ly sis are adopted, corresponding to a selection criterion on the im age charge o f a m inim um o f 40 photo-electrons. A cross-check analysis, p erform ed w ith the Im P A C T analysis (Parsons & H inton 2014) an d an independent calibration chain, yields com patible results.

2.2. B la za r s a m p le

T he data sets used are those o f blazars w ith know n red sh ift o b ­ served b y H .E.S.S. w ith a high significance (see Table 1) . They all belong to th e class o f high-frequency-peaked B L L ac objects.

T heir V H E em ission is therefore n o t expected to b e affected by the local b lazar environm ent. B lazars can som etim es show signs o f spectral variability correlated w ith their flux level, and this could bias the interpretation in term s o f th e E B L. To avoid this, d ata from sources w ith know n variability are d ivided into subsets w ithin logarithm ic flux bins, as in H E SS 2013. T hese subsets, la ­ beled by a num ber, are ordered b y increasing flux level. T he bulk o f the d ata sam ple is sim ilar to th at u sed in H E SS 2013; there are som e differences, w hich are m entio n ed below.

M rk 4 2 1 (z = 0.031, U lrich e t a l. 1975) is the first extra­

galactic source detected in the V H E dom ain (P unch et al. 1992) . This b rig h t and variable northern-sky blazar is observed by H .E.S.S. a t large zenith angles (>60°) (A haronian e t al. 2005) .

Table 1. Properties of the data sets used in this study, including the observation live time, the detection significance a , the source redshift z, and the energy range covered by the unfolded y-ray spectra.

D ata set L ive tim e a z E min " E max

(hours) (TeV)

M rk 4 2 1 (1) 4.9 89.6 0.031 1 . 41 - 14.9

M rk 421 (2) 3.8 122 0.031 1.2 2 - 15.9

M rk 421 (3) 2.9 123 0.031 1 .1 9 - 19.5

M rk 421 (4) 3.3 96.2 0.031 1.6 - 16.5

M rk 421 (5) 1.6 46.0 0.031 1 .5 - 15.2

M rk 5 0 1 1.8 66.7 0.034 1 .9 - 19.5

PKS 2 0 0 5 -4 8 9 (1) 71.2 28.8 0.071 0.29 - 1.6 PKS 2 0 0 5 -4 8 9 (2) 18.7 29.2 0.071 0.29 - 3 .0 PKS 2 1 5 5 - 3 0 4 (1) 7.4 94.8 0.116 0.24 - 4 .6 PKS 2 1 5 5 - 3 0 4 (2) 6.1 119 0.116 0 .2 4 - 1.98 PKS 2 1 5 5 - 3 0 4 (3) 5.5 187 0.116 0.24 - 3 .7 PKS 2 1 5 5 - 3 0 4 (4) 2 .6 135 0.116 0 .2 4 - 2.44 PKS 2 1 5 5 - 3 0 4 (5) 3.5 227 0.116 0.24 - 4 .6 PKS 2 1 5 5 - 3 0 4 (6 ) 1.3 172 0.116 0.29 - 4 .6 PKS 2 1 5 5 - 3 0 4 (7) 1.3 2 0 0 0.116 0.29 - 3 .6 PKS 2 1 5 5 - 3 0 4 (8 ) 25.4 111 0.116 0.19 - 3 .7

1ES 0 2 2 9+ 200 57.7 11.6 0.14 0 .4 - 2.8

H 2 3 5 6 -3 0 9 92.6 19.6 0.165 0 .1 9 - 1.98 1ES 1 1 0 1 -2 3 2 58.2 16.8 0.186 0 .1 9 - 1.98 1ES 0 3 4 7 -1 2 1 33.9 14.1 0.188 0.19 - 6 .9

1ES 0414+ 009 73.7 9.6 0.287 p VO 1 0.69

As a consequence, the energy threshold is high (> 1 TeV), due to the strong atm ospheric absorption o f C herenkov light. O n the other hand, the effective area is relatively large a t higher energies, resulting in a spectrum extending above 10 TeV. In addition to the 2004 observations on M rk 4 2 1 (labels 1 to 3), data taken during the 2010 high state (T luczykont 2010) are added (labels 4 and 5). M rk 5 0 1 (z = 0.034, M oles e t al. 1987) is the second extragalactic V H E source detected (Q uinn et al.

1996) and is also observed b y H .E .S.S. a t large zenith angles.

D ata taken during th e 2014 high state are used (C ologna et al.

2016) . T hese low -redshift spectra a t m ulti-TeV energies are key to probing the m id -IR region o f the E B L spectrum . F or PKS 2 0 0 5 -4 8 9 (z = 0.071, F alom o e t a l. 1987) , the d ata used here are identical to those used in H E SS 2013, and detailed in A cero (2010) and A bram ow ski et al. (2 0 1 1 ). PKS 2 1 5 5 -3 0 4 (z = 0.116, F alom o e t a l. 1993) is a very b rig h t b lazar ex ten ­ sively studied by H .E .S.S. As in H E SS 2013, the d ata o f the ex ­ ceptional Ju ly 2006 high state are u sed (labels 1 to 7), together w ith observations o f the 2008 low state (label 8). T hese very high-significance d ata sets yield excellent quality spectra th at are crucial for an unam biguous identification o f the E B L absorption pattern. F or 1ES 0 2 2 9 + 2 0 0 (z = 0.14, S chachter et al. 1993) , H 2 3 5 6 -3 0 9 (z = 0.165, Jones e t a l. 2009) , 1ES 1 1 0 1 -2 3 2 (z = 0.186, R em illard et al. 1989) , and 1ES 0 3 4 7 -1 2 1 (z = 0.188, W oo e t al. 2005), th e d ata are identical to those used in H E SS 2013. T he b lazar 1ES 0 4 1 4+ 009 (z = 0.287, H alpern et al.

1991) is also added to the sam ple. This d istan t source was o b ­ served by H .E .S.S. from 2005 to 2 009 (A bram ow ski e t al. 2012) .

2.3. S p e c tra l d ec o n vo lu tio n

P ublished spectra b y H .E .S.S. are usually obtained by m eans o f a forw ard-folding m eth o d (Piron et al. 2001) fo r w hich an assum ption on the spectral shape is required. T he results o f this procedure are spectral param eters an d their associated er­

rors. Spectral points can then b e constructed in different en ­ ergy bins (using the ratio o f the observed signal to the signal

p redicted b y the fitted shape in each bin), b u t these points are n o t a d irect resu lt o f the forw ard-folding procedure. The p resen t study follow s a different approach: the energy spectrum o f each d ata set is obtained using a B ayesian u nfolding technique b ased on A lb ert et al. (2007; an d already used by H .E .S.S. in A bram ow ski et al. 2013a) in o rder to directly obtain spectral points independently o f any a prio ri spectral shape, together w ith their correlations. This is a key asp ect o f this new analysis since these unfolded spectra allow the exploration o f spectral patterns.

T he energy threshold used in the spectral deconvolution is defined as the energy at w hich th e effective area reaches 15%

o f its m axim um value. This is a standard p rocedure used in H .E.S.S. (A bram ow ski e t al. 2 0 1 4 ; A bdalla et al. 20 1 7 ) . A fixed logarithm ic binning in energy is chosen for the deconvolution o f each spectrum , adapted to the energy resolution. A m inim um significance o f 2 a p er bin is required for a spectral p o in t to be defined. T he high-energy en d o f the ran g e indicated in Table 1 reflects the tail o f the significance distribution in energy. The consistency o f the unfolded spectral points have been verified and are in excellent agreem ent w ith the residual points o f the above-m entioned forw ard-folding procedure. T he 21 unfolded spectra o f the sam ple yield a total o f 247 spectral points.

3. Null hypothesis: fits without EBL

This determ ination o f the E B L is b ased on the observation o f fe a ­ tures in observed V H E spectra. A ssum ptions m ad e on intrinsic spectra are therefore essential. In this section, these assum ptions are presented, and it is show n th at the assum ption o f featureless spectra w ith n o E B L leads to a p o o r fit o f the data. This calls for additional degrees o f freedom in the d ata interpretation. T he best fit w ithout E B L w ill b e considered later as the null hypothesis. In the n ex t sections, it is show n th at w hen these additional degrees o f freedom reflect E B L levels in different bands, th e fit is sig­

nificantly im proved and interpreted as evidence o f a spectrally resolved E B L detection.

Fig. 1. Fit residuals of the featureless spectral shapes, as a function of energy. a) Residuals of the whole sample of spectra to the power-law fit.

b) Residuals of the whole sample to the log-parabola fit. c) Example residuals to the log-parabola fit for the subset PKS 2155-304 (5).

The simplest description of the energy spectrum of a non- thermal y-ray source like a blazar is the two-parameter power- law function

® p w l(e t ) - ®Q( E7/ E q ) a.

% P ( E y ) - ® 0 ( E y / E 0 ) " “ " ^ lo g (£ Y/E 0 ) ,

photon m om enta,

( 1)

Fitting each spectrum j of the blazar sample with a power-law yields overall £ j X /PWL - 1472.8, for 205 degrees of freedom.

The fit residuals have a large dispersion, as shown in Fig. 1a.

In addition, a visible modulation indicates the need for a more elaborate parameterization. Spectral curvature is introduced with the three-parameter log-parabola function

(3)

w here u - 1 - cos(0), and ethr(E r , z) - E is the threshold energy dictated b y kinem atics. T he cross section fo r p air p ro d u c­

tion (B reit & W heeler 1934; G ould & S chreder 1967) is defined as

(

2

(4)

with

0 the additional curvature parameter. The use of the log- parabola function improves the fit significantly (at the 3 0 ^ level) as 2 j

- 281.07, for 184 degrees of freedom. However, a modulation in the distribution of the fit residuals is still present, as illustrated in Fig. 1b. As an example, Fig. 1c shows a residual distribution isolated from Fig. 1b, for the log-parabola fit of the subset PKS 2155-304 (5).

The same observation can be made using other featureless spectral shapes, as long as no intrinsic irregularities are consid­

ered. It is then natural to try to interpret these modulations of the flux residuals as the effect of EBL attenuation. The properties of these energy-dependent modulations and their redshift depen­

dence are the central point of this study to measure the EBL SED.

In the following these modulations, which are not accounted for by featureless intrinsic shapes, are translated in terms of spec­

trally resolved EBL levels.

4. EBL optical depth

The extragalactic medium at a given redshift

is filled with EBL photons of proper number density « e b l(£. z) at proper energy e. The opacity of this medium for y rays of observed energy E

coming from a source at redshift zs is encoded in the optical depth t(E y, zs) (Gould & Schreder 1967; Stecker et al. 1992). It consists of an integration over z, e, and the angle 9 between the

where

is the velocity of the outgoing electron and positron in the center of a mass system, and

is the Thomson cross section.

A flat ACDM cosmology with Hubble constant H0, matter density parameter Om, and dark energy density parameter

is considered. The distance element in Eq. (3) then reads

(5)

where the values

70 km s Mpc , Om - 0.3, and Oa - 0.7 are assumed. The most influential cosmological parameter is the Hubble constant as

scales linearly with 1/H

. This generic choice of H

is in line with the latest

results obtained from cosmic microwave background data (Planck Collaboration XIII 2016) and with the results ob­

tained from more local constraints (Riess etal. 2016). The de­

pendence of the results on the precise choice for

is negligible with respect to the sensitivity of the method. For a detailed study of the influence of cosmological parameters on y-ray attenuation (see, e.g., Dominguez & Prada 2013).

The evolution of the EBL in Eq. (3) with redshift is ac­

counted for by decoupling the local

- 0) EBL SED and an evolution function,

d

d

d

—-— (

z)

d

— ---(

0)

/ ( ¾

z)

d

d

^ t

fP

f P

, ,

t ( E y, a ) -

J

X Jo duU2 ^ 77 [^(E y(1 + z), £, u ) ] ,

&yy (f>) - — (1 - P 2) (3 - ,^4 )ln | j __ ^ i - 2A 2 - f t 2) ,

dZ - c (H q(1 + z) V ^m(1 + z)3 + Qa) ,

w here e0 = e /(1 + z) is the E B L energy at z = 0. T he evo­

lution function f (e0, z) is extracted from the m o d el given in F ranceschini et al. (2008) using th e ratio o f the SED at a red- shift z to its value a t z = 0. T he influence on the results o f this m odel-dependent in gredient for the E B L evolution w ith redshift is w eak (see Sect. 6.2) .

T he observed energy spectrum $ obs( £ r ) o f an extragalac- tic source is the convolution o f its intrinsic spectrum ^ nt(E r ) w ith the E B L absorption effect $ obs(E r ) = ®int(Er ) x e -T(Er,zs).

E xtragalactic back g ro u n d light absorption then leaves a redshift and energy-dependent im print on the observed V H E spectra o f blazars.

Intrinsic spectral shapes are d escribed w ith log-parabolas.

This n aturally includes pow er law s in cases w here the fit p refers vanishing curvature param eters. This choice is intended to avoid attributing the entire origin o f spectral curvature to th e E B L. As only positive values o f a an d j3 are considered, another im plicit assum ption is th e non-convexity o f th e intrinsic spectra (e.g., D w ek & K rennrich 2 0 0 5 ; D w ek et al. 2005) . T hese sim ple co n ­ siderations ensure th at this E B L determ ination does n o t rely on specific assum ptions abo u t the underlying acceleration m e ch a­

nism b eh in d th e V H E y -ray em ission o f blazars. T he consider­

ation o f m ore com plex intrinsic shapes does n o t lead to an im ­ provem ent in individual fit qualities.

5. Method

5.1. P a ra m eteriza tio n o f th e EBL S E D a n d intrinsic b la za r sp e c tra

As show n in Sect. 3 , energy-dependent m odulations in the resid ­ uals o f spectral fits w ith featureless functions call for additional degrees o f freedom . T hese can b e interpreted in term s o f E B L ab ­ sorption. A determ ination o f the E B L is possible by confronting y-ray data w ith different E B L hypotheses, and in th e p rese n t ap ­ proach these hypotheses are independent o f E B L m odels. A p re ­ lim inary study testing E B L shapes as splines constructed upon a grid in energy density (as in M azin & R aue 2007) show ed in ­ deed th at the shape o f the E B L was accessible w ith H .E.S.S.

data (L orentz e t al. 2015) . In th e p rese n t study a different and m ore robust m eth o d is used: E B L rela ted degrees o f freedom are introduced as continuous levels o f E B L intensity in different bands over the ran g e o f interest for y -ray absorption. This ap ­ proach allow s a m o re accurate estim ation o f uncertainties and a m ore m eaningful statistical treatm ent o f data sets than does the use o f splines on a grid.

T he local E B L energy density is decom posed into connected energy bands w ith bounds [e;, e!+1] and content p ; > 0 ,

«0 d2 BL = Z -

w here - ( * , > = { 0 ^ i ^ ’2 <7 )

This p aram eterization o f the E B L SED is injected into the optical depth calculation (Eq. (3)). T he set o f E B L levels {p} is adjusted to fit th e absorption pattern in y -ra y data. This approach is sim ­ ilar to the one used in B iller et al. ( 1998) to derive u pper lim its on th e E B L SED.

T he local E B L SED is divided into fo u r energy (w avelength) bands w ith equal logarithm ic w idths. This sim ple choice is found to b e optim al in term s o f sensitivity. Increasing the n um ber o f subdivisions does n o t significantly im prove th e fit quality, but leads to an increase in the errors on the E B L levels. T he low - energy bo u n d eo,min (or equivalently th e high-w avelength bound

^o,max) o f the local E B L ran g e corresponds to the threshold for pair creation w ith the m o st energetic y rays o f the b lazar sam ple follow ing th e threshold relation previously m entioned:

e = hc = 2m2c4 (8)

e0,min = , = ^ / - 1 , \2 • (8)

^0,max EyP (1 + zs)

T he high-energy (low -w avelength) bo u n d o f the E B L ran g e is chosen b eyond th e peak o f the cross section for interaction w ith y rays in the low est energy spectral p o in t o f the sam ple and ad ­ ju ste d a p osteriori as th e energy at w hich th e sensitivity in this

band is seen to decrease.

5.2. J o in t fit

F or each individual d ata set, a jo in t fit o f th e E B L levels and intrinsic spectral param eters is perform ed. T he covariance m a ­ trix C T determ ined in the unfolding p rocedure is used in o rder to take into account th e correlations betw een spectral points in the X 2 m inim ization. T he m inim ized function is

X 2 = (O test - ^obs)TC T1(Otest - O obsX (9)

w here O obs is the vector o f observed spectral points and O test is the vector o f tested functions. T hese te st functions include both E B L param eters {p;} an d intrinsic spectral param eters ( 0 0, a,j3),

$ test(E r , zs) = $int(E r , ¢ 0, a ,P ) x e -T(EY,zs,{pi}), (10)

resulting in a seven-param eter fit.

T he four E B L levels are independent in the fit o f each in d i­

vidual spectrum and are com bined later, as d escribed in Sect. 6 .1 . This approach allow s a clear identification o f th e contribution o f each spectrum to th e overall results at different w avelengths.

6. Results

Fits o f all spectra are p erform ed follow ing the procedure d e­

scribed above. F or 1ES 0 4 1 4+ 009 the intrinsic spectrum is re ­ stricted to a pow er-law , due to the lim ited n um ber o f degrees o f freedom available.

A ll intrinsic spectral param eters are found to b e reasonable, in agreem ent w ith typical em ission m odels. It should b e noted th at this p aper focuses on the E B L m easurem ent, and discussion o f intrinsic spectral param eters w ill b e d etailed in a forthcom ing paper. Individual fit qualities obtained w ith Eq. ( 10) are show n in Table 2 . T he fits are n o t im proved w hen considering m o re co m ­ plex intrinsic functions such as ones w ith cutoffs. L ow fit qu ali­

ties can b e rela ted to sm all spectral irregularities that cannot be accounted fo r in the param eterization. T he unfolding covariance m atrix can also reduce the fit quality.

T he goodness-of-fit estim ator is £ j x 2 LP+EBL. Its value after the jo in t fit is 176.7. C onsidering the fo u r additional E B L d e­

grees o f freedom as com m on param eters and using W ilks’ th e o ­ rem , this can b e translated into an E B L detection significance o f 9 .5 ^ w ith resp ect to the log-parabola hypothesis w ithout E B L.

F igure 2 a displays th e accum ulated residuals from the fit w ith E B L. T he m odulation seen in Fig. 1b is reduced, show ing th at th e addition o f E B L related degrees o f freedom provide a better description o f the data. F igure 2b show s this effect fo r the subset PKS 2 1 5 5 - 3 0 4 (5), to b e com pared w ith F ig. 1c.

Table 2. Fit qualities of the different data sets after the jo in t fit.

D ata set _{X ,LP+EBL/}

ndf

M rk 421 (1) 5.17/5

M rk 421 (2) 17.3/6

M rk 421 (3) 8.83/6

M rk 421 (4) 16.37/5

M rk 421 (5) 6.19/5

M rk501 14.18/5

PKS 2 0 0 5 -4 8 9 (1) 13.3 /2 PKS 2 0 0 5 -4 8 9 (2) 3.17/5 PKS 2 1 5 5 -3 0 4 (1 ) 9.1/6 PKS 2 1 5 5 -3 0 4 (2) 6.71/4 PKS 2 1 5 5 -3 0 4 (3) 11.76/7 PKS 2 1 5 5 -3 0 4 (4) 10.3 /5 PKS 2 1 5 5 -3 0 4 (5) 3.23/7 PKS 2 1 5 5 -3 0 4 (6) 4.37/7 PKS 2 1 5 5 -3 0 4 (7) 12.29/6 PKS 2 1 5 5 -3 0 4 (8) 19.9 /7

1ES 0229+200 2.07/1

H 2356 -3 0 9 6.21/3

1ES 1101-232 1.9 /5

1ES 0347-121 1.4/3

1ES 0414+009 3.2/1

local EBL energy density

(in units of eV m 3 ) into spe­

cific intensity

(in units of nW m- 2 sr- 1 ) following the relation

e

=

4

^ 0

o .

Individual results per data set reflect the relative sensitivity range of the different spectra (see Fig. 3). For instance, subsets of PKS 2155-304 show optimal sensitivity in the 1.1-4.94 um band with precisely fitted EBL levels. Subsets of Mrk421 and Mrk 501 lead to precise measurements at larger wavelengths. Al­

ternatively, when the range covered by a given spectrum does not constrain the EBL in a given band, the corresponding uncer­

tainty on the fitted level is large. This behavior shows the ability of the method to probe the different ranges in EBL wavelength depending on the source properties (redshift, accuracy of spec­

tral points measurements, covered energy range). When no clear spectral modulation is identified, a softer or more curved intrin­

sic spectrum compatible with null levels of EBL are preferred.

This would not be the case if intrinsic curvature was forbidden.

The method thus prevents an overall spectral curvature from be­

ing interpreted as an EBL detection because its signature - if present - is imposed to be a more complex feature.

In the following, the results are presented by converting the

6.1. S p ec tra lly r e s o lv e d EBL le ve ls

Individual signals are combined to obtain a collective EBL mea­

surement. This approach takes advantage of the large sample of high-quality spectra from sources at various redshifts. Indeed, the EBL is expected to have a coherent effect that can be inter­

preted collectively.

The combined EBL level in a band is obtained as an error- weighted average over all data sets {

Fig. 2. a ) Residuals of the w hole sample of spectra for the log-parabola fit with EBL. b) Example residuals to the log-parabola fit w ith EBL for the subset PKS 2 1 5 5 -3 0 4 (5).

sample variance in each wavelength band by the corresponding reduced ^ 2:

(12)

(

11

The error on a combined EBL level takes into account both in­

dividual errors and the dispersion of the

individual measure­

ments around the averaged value. This is done by weighting the

Following this approach, the uncertainty on a combined EBL level is slightly corrected to yield conservative results.

Individual EBL levels should essentially be compatible with each other as the EBL is assumed to be a di

use isotropic back­

ground. Potential anisotropies of the EBL (Furniss etal. 2015;

Abdalla & Bottcher 2017) are estimated to be beyond the sensi­

tivity of this EBL measurement method. The dispersion of indi­

vidual values can reflect systematic uncertainties and potential limitations of the procedure. The latter can be due to the fit of patterns in spectra that might not be related to EBL absorption and not accounted for in the intrinsic spectrum parameterization.

However, to avoid introducing a bias, once the method is fixed, any kind of individual tuning of parameters is forbidden, and all fits are performed in one single blind procedure. Equation (12) ensures that significant deviations to the average value degrade the precision on the EBL measurement. Of the 84 points dis­

played in Fig. 3, only a few deviate from the average value. The thin red lines only represent the statistical uncertainties on the

_ 1 1 y [p ij - <Pi>)2 a{pP _ ) 2 j (n -

, \ 2 j Pi,j/ a '2pi,j {Pi) _ — — =2—

2 j ^ - 5

Fig. 3. Individual EBL levels and errors in the different wavelength bands obtained from the fit of each spectrum. The solid red lines represent the combined values and statistical errors.

Fig. 4. Measured EBL spectrum. Obtained levels are represented by the red points. Horizontal lines represent the bandwidth over which the inte­

grated EBL levels apply. The vertical solid lines represent the 1-ix (statistical) errors from Eq. (12). The dashed lines represent the systematic un­

certainties on fitted EBL levels conservatively estimated as explained in the text. Direct constraints on the EBL collected from Dwek & Krennrich (2013) and Biteau & Williams (2015) are shown; lower and upper limits are represented by green and brown arrows, respectively.

combined EBL levels. For the lowest wavelength band (0.25­

1.1 pm), Eq. (12) slightly reduces the size of the statistical un­

certainty on the combined level because of the underdispersion of individual levels due to their very large error bars. The large uncertainty on this combined level due to the poorly constrained individual measurements is properly taken into account with the consideration of systematic uncertainties in Sect. 6.2.

The combined EBL levels are summarized in Table 3 and shown in Fig. 4. Details on the estimation of systematic uncer­

tainties are given in Sect. 6.2. Comparisons with various EBL constraints and models are left for the discussion in Sect. 6.3.

6.2. S y s te m a tic u n c e rta in ties

In addition to the previously mentioned dispersion of indi­

vidual values, different sources of systematic uncertainties are investigated.

Table 3. Combined EBL levels (p, in nW m 2 sr 1) in the different wave­

length bands (A, in pm ).

A A

A

p

(sys) p

(sys)

0.52 0.25 1.11 6.42 4.02 (0) 8.82 (14.5) 2.33 1.11 4.94 8.67 7.63 (6.68) 9.71 (11.80) 10.44 4.94 22.07 3.10 2.62(1.16) 3.59(4.17) 46.6 22.07 98.60 4.17 3.71 (2.55) 4.63 (6.30)

Systematic uncertainties related to the EBL evolution hy­

pothesis are estimated considering the toy-model evolution pa­

rameterization used in other model-independent approaches

to determine the EBL (Meyer et al. 2012; Biteau & Williams

2015). It consists of a global rescaling of the photon density with

respect to the cosmological expansion (1 + z)3 ^ (1 + z)3-fevd,

where the value of f evoi is chosen in order to mimic the evolution

function o f E B L m odels. A dopting the typical value / evoi - 1.2 (R a u e & M a z in 2008) leads to differences o f less than 5% in E B L levels w ith resp ect to th e case w here th e evolution fu n c­

tion is extracted from F ranceschini et al. (2 008), w ithin statisti­

cal errors.

S ystem atic uncertainties related to th e energy scale in y-ray spectra m easurem ents cou ld originate from variations o f the C herenkov light yield, due to fluctuations in atm ospheric tran s­

parency n o t accounted for in the sim ulation, m ism atches b e ­ tw een real and sim ulated m irror reflectivities, etc. (H ahn et al.

2014) . A system atic energy shift o f ±15% is assum ed. It re p re ­ sents a conservative estim ate o f the absolute energy scale u ncer­

tainty w ith H .E.S.S. (A haronian e t al. 2006b) . This energy shift is applied a t th e spectrum level for all d ata sets an d th e w hole procedure is redone. T he E B L w avelength range is shifted a c ­ cordingly, so th at the p o ssibility is left for different E B L levels in different bands to induce identical patterns in spectra. T he o b ­ served variations on E B L levels are on the order o f 10%, sym ­ m etric w ith resp ect to the central value and global over the w ave­

length range, w ith sim ilar goodness o f fit in each case.

I f the E B L w avelength ran g e is n o t shifted according to the y -ray energy scale, a m ism atch betw een the spectrum en ­ ergy scale and the relative position o f the w avelength bands is introduced. This is equivalent to investigating the effects o f bin-shifting in the w avelength bands an d can lead to signifi­

cant changes in the m easured E B L levels, reflecting th e level o f degeneracy betw een intrinsic spectra an d fitted E B L levels.

T he com bined system atic effect o f shifts in the energy scale, changes in the w avelength ran g e and changes in fitted intrinsic spectra leads to significant variations from 10% to 70% o f the central E B L level. T hese uncertainties strongly depend on the band considered and are n o t sym m etric in intensity. In this way, conservative w avelength-dependent system atic uncertainties are obtained.

T he influence o f the w idth o f the w avelength bands and their nu m b er is also investigated. In addition to changes that can n a t­

urally arise from th e variations o f th e E B L SED over one band, integrated E B L levels can fluctuate due to th e different ab so rp ­ tion patterns available w hen using bands o f different size and num ber. T he n um ber o f bands is lim ited b y the degrees o f fre e­

dom available in the jo in t fit o f each spectrum . T he use o f m ore (and sm aller) bands is then only possible for spectra w ith suf­

ficient degrees o f freedom . U sing w avelength bands w ith larger and sm aller w idths, the changes in E B L levels are found to be d ependent on the w avelength ran g e considered, from negligible variations up to 40% variations.

A n alternative approach fitting sim ultaneously all spectra w ith com m on E B L param eters is also considered. W ith this global approach the com bination o f E B L levels obtained from individual spectra is n o t needed, b u t th e inform ation co n cern ­ ing the contribution o f each spectrum to the m easurem ent is lost.

C onsistent E B L levels are o btained from this global fit. T he level in the low est w avelength b an d (0.2 5 -1 .1 1 u m ) appears poorly constrained in the global fit and is co m patible w ith a null level o f E B L . This low est w avelength b an d is at the lim it o f th e sen ­ sitivity o f this study, as is already apparent from th e individual m easurem ents show n in Fig. 3 .

T he envelope o f largest variations corresponding to th e d if­

ferent kinds o f p otential system atic errors are rep resen ted in Fig. 4 as dashed lines. This b ehavior shows the relative in d e­

term inacy fo r th e 0 .2 5 -1 .1 1 u m and 4 .9 4 -2 2 .0 7 u m bands, but also the stronger signals in th e 0 .2 5 -1 .1 1 u m and in the 22.07­

98.6 u m bands w hich are clearly significant b eyond system atic uncertainties.

6.3. D isc u ssio n

S ensitivity to th e shape and norm alization o f the E B L SED is achieved using only V H E spectra obtained w ith H .E.S.S. A l­

though E B L levels w ere left free to cover a w ide ran g e o f p o ssi­

bilities and w ere thus n o t constrained in th e fits, the results do not conflict w ith the strict low er lim its from th e galaxy counts. This is an interesting p o in t as the tw o m ethods (y rays and galaxy counts) are com pletely independent, and the E B L levels in the presen t analysis w ere left free to vary betw een zero and an a rb i­

trary value. T he o btained results are consistent w ith state-of-the- art E B L m odels an d are in general agreem ent w ith o ther y-ray constraints, as show n in F ig. 5 . A gain, it should be noted that no prio r from the displayed m odels w as used in the H .E.S.S.

m easurem ent (apart from the evolution factor, w hich does not strongly influence th e z - 0 results). T he obtained results are com patible w ith the m o d el scaling o f H E SS 2013, although an extended d ata set is used and th at the treatm ent o f the data is very different. T he w avelength ran g e p ro b ed by H E SS 2013 was conservatively restricted to the central value o f the pair-creation cross section, neglecting its w idth. In the presen t study, the w ave­

length ran g e p robed is extended farther into the infrared because optical depth values m u st b e com puted over th e w hole k in e m at­

ically allow ed ran g e for pair-creation w ith th e m o st energetic y rays o f the sam ple, as described in Sect. 4 . T he obtained E B L levels close to low er lim its in th e optical ran g e are in line w ith Ferm i LAT results (A ckerm ann e t al. 2012) probing th e E B L at higher redshifts and a t low er w avelengths and also w ith the upper lim it obtained follow ing the d etection o f the h igh-redshift quasar PKS 1441+25 a t z - 0.94 by V ERITAS (A beysekara et al. 2015) and M A G IC (A h n e n e ta l. 2015) . T he results are also in g en ­ eral agreem ent w ith other constraints n o t represented in F ig. 5 obtained using y rays (A bram ow ski et al. 2 0 1 3 b ; A h n e n e ta l.

2016) o r w ith th e results o f em pirical approaches to the d e­

term ination o f the E B L SED (H elgason & K ashlinsky 2 0 1 2 ; Stecker e t al. 2016) .

W hile an im portant conclusion o f this w ork is to show that H .E.S.S. y -ray spectra alone contain enough inform ation to d e­

term ine the E B L shape and norm alization, the sensitivity o f this approach rem ains lim ited. T he coarse E B L binning achievable and the conservatively estim ated uncertainties on E B L levels show th at a fine spectroscopy o f the E B L SED (resolving fine substructures in the E B L spectrum , e.g., due to dust subcom po­

nents) is out o f reach using only the p rese n t V H E data. T he c o m ­ patibility betw een H .E .S.S. m easurem ents and the low er lim its from galaxy counts does n o t suggest a transparency anom aly o f the U niverse to V H E y rays (as hinted a t in H orns & M eyer 2 0 1 2) for the red sh ift ran g e considered.

T he results in term s o f E B L intensity can b e translated into a corresponding y -ray horizon. T he y -ray horizon for

t - 1 is a standard illustration o f the E B L absorption effects (Fazio & S tecker 1970) . It corresponds to th e typical attenuation length o f y rays a t a given observed energy. T he t - 1 energy- red sh ift horizon envelopes corresponding to the m easured E B L levels and th eir errors are rep resen ted in F ig. 6 . Is o -r curves o f selected m odels are show n for com parison. T hese results in term s o f y -ra y horizon o r optical depths are also com patible w ith the horizon derived from the SED o f state-of-the-art E B L m o d ­ els. H ere th e lim ited sensitivity o f the approach also appears, as the consideration o f system atic uncertainties significantly en ­ larges the w idth o f the horizon envelope. This shows the diffi­

culty in interpreting such results in term s o f transparency an o m a­

lies th at could b e due to second-order propagation effects.

Fig. 5. Combined EBL levels (red points) compared with various constraints and models. The dashed black line represents the local SED of the model given in Franceschini et al. (2008), which is used as the template for the HESS2013 model scaling (yellow area). Additional models repre­

sented are Dominguez et al. (2011; grey dot-dashed line), Finke et al. (2010; green dashed line), Gilmore et al. (2012; blue dot-dashed line), and the lower-limit model of Kneiske & Dole (2010; purple dot-dashed line). The model-independent upper limit using VHE and HE data Meyer et al.

(2012) is shown as a thin brown line. The model-independent measurement of Biteau & Williams (2015) restricted to the use of VHE data is represented by blue points.

Fig. 6. Obtained y-ray horizon for t = 1 in the redshift range covered by the blazar sample and comparison with selected models (Franceschini et al.

2008; Dominguez et al. 2011; Finke et al. 2010; Gilmore et al. 2012; Kneiske & Dole 2010).

The reduction of uncertainties would of course be possible using additional data and priors. For instance, using data from

LAT at lower energy, the degeneracy between intrinsic spectra and EBL absorption can be reduced. This assumes con­

tinuity in energy of the intrinsic spectrum behavior, and can also introduce additional systematics due to absolute flux level uncer­

tainties. Taking into account direct EBL measurements and strict lower limits from galaxy counts can de facto restrict the range of variations for EBL levels. Such strong priors were avoided here in order to address the question of EBL information contained only in VHE spectra obtained with H.E.S.S. The detailed study of the intrinsic spectra obtained with this EBL results is left for a forthcoming dedicated paper.

7. Summary and conclusion

A determination of the EBL SED with the H.E.S.S. array of Cherenkov telescopes is presented. This is achieved using a new method: coherent patterns in the high-quality unfolded spec­

tra of blazars observed by H.E.S.S. are translated into EBL intensity levels resolved in wavelength under the assumption that intrinsic spectra are described by smooth concave shapes.

The EBL signature is preferred at the 9 .5 a level compared to the null hypothesis. Combined EBL levels are compati­

ble with current constraints and models, and no indication of

an opacity anomaly is found. This robust result demonstrates

for the first time the capability of H.E.S.S. to measure the

E B L SED independently o f any existing E B L constraints and m odels.

Acknowledgements. The support of the Namibian authorities and of the Univer­

sity of Namibia in facilitating the construction and operation of H.E.S.S. is grate­

fully acknowledged, as is the support by the German Ministry for Education and Research (BMBF), the Max Planck Society, the German Research Foundation (DFG), the Alexander von Humboldt Foundation, the Deutsche Forschungsge- meinschaft, the French Ministry for Research, the CNRS-IN2P3 and the As- troparticle Interdisciplinary Programme of the CNRS, the UK Science and Tech­

nology Facilities Council (STFC), the IPNP of the Charles University, the Czech Science Foundation, the Polish National Science Centre, the South African De­

partment of Science and Technology and National Research Foundation, the Uni­

versity of Namibia, the National Commission on Research, Science & Technol­

ogy of Namibia (NCRST), the Innsbruck University, the Austrian Science Fund (FWF), and the Austrian Federal Ministry for Science, Research and Economy, the University of Adelaide and the Australian Research Council, the Japan So­

ciety for the Promotion of Science and by the University of Amsterdam. We appreciate the excellent work of the technical support staff in Berlin, Durham, Hamburg, Heidelberg, Palaiseau, Paris, Saclay, and in Namibia in the construc­

tion and operation of the equipment. This work benefited from services pro­

vided by the H.E.S.S. Virtual Organisation, supported by the national resource providers of the EGI Federation.

References

Abdalla, H., & Bottcher, M. 2017, ApJ, 835, 237

Abdalla, H., Abramowski, A., Aharonian, F., et al. 2017, A&A, 600, A89 Abeysekara, A. U., Archambault, S., Archer, A., et al. 2015, ApJ, 815, L22 Abramowski, A., Acero, F., Aharonian, F., et al. 2011, A&A, 533, A110 Abramowski, A., Acero, F., Aharonian, F., et al. 2012, A&A, 538, A103 Abramowski, A., Acero, F., Aharonian, F., et al. 2013a, Phys. Rev. D, 88,102003 Abramowski, A., Acero, F., Aharonian, F., et al. 2013b, A&A, 559, A136 Abramowski, A., Acero, F., Aharonian, F., et al. 2013c, A&A, 550, A4 Abramowski, A., Aharonian, F., Ait Benkhali, F., et al. 2014, A&A, 564, A9 Acero, F. 2010, A&A, 511, A52

Ackermann, M., Ajello, M., Allafort, A., et al. 2012, Science, 338, 1190 Aharonian, F. A., Coppi, P. S., & Voelk, H. J. 1994, ApJ, 423, L5 Aharonian, F., Akhperjanian, A. G., Aye, K.-M., et al. 2005, A&A, 437, 95 Aharonian, F., Akhperjanian, A. G., Bazer-Bachi, A. R., et al. 2006a, Nature,

440, 1018

Aharonian, F., Akhperjanian, A. G., Bazer-Bachi, A. R., et al. 2006b, A&A, 457, 899

Ahnen, M. L., Ansoldi, S., Antonelli, L. A., et al. 2015, ApJ, 815, L23 Ahnen, M. L., Ansoldi, S., Antonelli, L. A., et al. 2016, A&A, 590, A24 Albert, J., Aliu, E., Anderhub, H., et al. 2007, Nucl. Instr. Meth. A, 583, 494 Biller, S. D., Buckley, J., Burdett, A., et al. 1998, Phys. Rev. Lett., 80, 2992 Biteau, J., & Williams, D. A. 2015, ApJ, 812, 60

Breit, G., & Wheeler, J. A. 1934, Phys. Rev., 46, 1087

Cologna, G., Chakraborty, N., Mohamed, M., et al. 2016, PoS, ICRC 2015, 761 Costamante, L. 2013, Int. J. Mod. Phys. D, 22, 1330025

de Naurois, M., & Rolland, L. 2009, Astropart. Phys., 32, 231 Dominguez, A., & Prada, F. 2013, ApJ, 771, L34

Dominguez, A., Primack, J. R., Rosario, D. J., et al. 2011, MNRAS, 410, 2556 Dwek, E., & Krennrich, F. 2005, ApJ, 618, 657

Dwek, E., & Krennrich, F. 2013, Astropart. Phys., 43, 112 Dwek, E., Krennrich, F., & Arendt, R. G. 2005, ApJ, 634, 155

Falomo, R., Maraschi, L., Treves, A., & Tanzi, E. G. 1987, ApJ, 318, L39 Falomo, R., Pesce, J. E., & Treves, A. 1993, ApJ, 411, L63

Fazio, G. G., & Stecker, F. W. 1970, Nature, 226, 135

Finke, J. D., Razzaque, S., & Dermer, C. D. 2010, ApJ, 712, 238 Franceschini, A., Rodighiero, G., & Vaccari, M. 2008, A&A, 487, 837 Furniss, A., Sutter, P. M., Primack, J. R., & Dominguez, A. 2015, MNRAS, 446,

2267

Gilmore, R. C., Somerville, R. S., Primack, J. R., & Dominguez, A. 2012, MNRAS, 422, 3189

Gould, R. J., & Schreder, G. P. 1967, Phys. Rev., 155, 1404

Hahn, J., de los Reyes, R., Bernlohr, K., et al. 2014, Astropart. Phys., 54, 25 Halpern, J. P., Chen, V. S., Madejski, G. M., & Chanan, G. A. 1991, AJ, 101, 818 Hauser, M. G., & Dwek, E. 2001, ARA&A, 39, 249

Hauser, M. G., Arendt, R. G., Kelsall, T., et al. 1998, ApJ, 508, 25 Helgason, K., & Kashlinsky, A. 2012, ApJ, 758, L13

Horns, D., & Meyer, M. 2012, J. Cosmol. Astropart. Phys., 1202, 033 Jacob, U., & Piran, T. 2008, Phys. Rev. D, 78, 124010

Jones, D. H., Read, M. A., Saunders, W., et al. 2009, MNRAS, 399, 683 Kneiske, T. M., & Dole, H. 2010, A&A, 515, A19

Lorentz, M., Brun, P., & Sanchez, D. 2015, in Proc. 34th International Cosmic Ray Conference (ICRC 2015), 777

Matsumoto, T., Matsuura, S., Murakami, H., et al. 2005, ApJ, 626, 31 Mazin, D., & Raue, M. 2007, A&A, 471, 439

Meyer, M., Raue, M., Mazin, D., & Horns, D. 2012, AIP Conf. Proc., 1505, 602 Moles, M., Masegosa, J., & del Olmo, A. 1987, AJ, 94, 1143

Nikishov, A. I. 1962, Sov. Phys. JETP, 14, 393

Orr, M. R., Krennrich, F., & Dwek, E. 2011, ApJ, 733, 77 Parsons, R. D., & Hinton, J. A. 2014, Astropart. Phys., 56, 26 Piron, F., Djannati-Atai, A., Punch, M., et al. 2001, A&A, 374, 895 Planck Collaboration XIII. 2016, A&A, 594, A13

Punch, M., Akerlof, C. W., Cawley, M. F., et al. 1992, Nature, 358, 477 Quinn, J., Akerlof, C. W., Biller, S., et al. 1996, ApJ, 456, L83 Raue, M., & Mazin, D. 2008, Int. J. Mod. Phys. D, 17, 1515

Remillard, R. A., Tuohy, I. R., Brissenden, R. J. V., et al. 1989, ApJ, 345, 140 Riess, A. G., Macri, L. M., Hoffmann, S. L., et al. 2016, ApJ, 826, 56 Sanchez-Conde, M. A., Paneque, D., Bloom, E., Prada, F., & Dominguez, A.

2009, Phys. Rev. D, 79, 123511

Schachter, J. F., Stocke, J. T., Perlman, E., et al. 1993, ApJ, 412, 541

Sinha, A., Sahayanathan, S., Misra, R., Godambe, S., & Acharya, B. S. 2014, ApJ, 795, 91

Stecker, F. W., & Glashow, S. L. 2001, Astropart. Phys., 16, 97 Stecker, F. W., de Jager, O. C., & Salamon, M. H. 1992, ApJ, 390, L49 Stecker, F. W., Scully, S. T., & Malkan, M. A. 2016, ApJ, 827, 6 Taylor, A. M., Vovk, I., & Neronov, A. 2011, A&A, 529, A144 Tluczykont, M. 2010, PoS, TEXAS 2010, 197

Ulrich, M.-H., Kinman, T. D., Lynds, C. R., Rieke, G. H., & Ekers, R. D. 1975, ApJ, 198, 261

Woo, J.-H., Urry, C. M., van der Marel, R. P., Lira, P., & Maza, J. 2005, ApJ, 631, 762

1 Centre for Space Research, North-West University, Potchefstroom 2520, South Africa

2 Universitat Hamburg, Institut fur Experimentalphysik, Luruper Chaussee 149, 22761 Hamburg, Germany

3 Max-Planck-Institut fur Kernphysik, PO Box 103980, 69029 Heidelberg, Germany

4 Dublin Institute for Advanced Studies, 31 Fitzwilliam Place, Dublin 2, Ireland

5 National Academy of Sciences of the Republic of Armenia, Marshall Baghramian Avenue, 24, 0019 Yerevan, Republic of Armenia

6 Yerevan Physics Institute, 2 Alikhanian Brothers St., 375036 Yerevan, Armenia

7 Institut fur Physik, Humboldt-Universitat zu Berlin, Newtonstr. 15, 12489 Berlin, Germany

8 University of Namibia, Department of Physics, Private Bag 13301, Windhoek, Namibia

9 GRAPPA, Anton Pannekoek Institute for Astronomy, University of Amsterdam, Science Park 904, 1098 XH Amsterdam,

The Netherlands

10 Department of Physics and Electrical Engineering, Linnaeus University, 351 95 Vaxjo, Sweden

11 Institut fur Theoretische Physik, Lehrstuhl IV: Weltraum und Astro- physik, Ruhr-Universitat Bochum, 44780 Bochum, Germany 12 GRAPPA, Anton Pannekoek Institute for Astronomy and Institute of

High-Energy Physics, University of Amsterdam, Science Park 904, 1098 XH Amsterdam, The Netherlands

13 Institut fur Astro- und Teilchenphysik, Leopold-Franzens- Universitat Innsbruck, 6020 Innsbruck, Austria

14 School of Physical Sciences, University of Adelaide, Adelaide 5005, Australia

15 LUTH, Observatoire de Paris, PSL Research University, CNRS, Universitd Paris Diderot, 5 Place Jules Janssen, 92190 Meudon, France

16 Sorbonne Universitds, UPMC Universitd Paris 06, Universitd Paris Diderot, Sorbonne Paris, Citó, CNRS, Laboratoire de Physique Nucldaire et de Hautes Energies (LPNHE), 4 place Jussieu, 75252 Paris Cedex 5, France

17 Laboratoire Univers et Particules de Montpellier, Universitd Montpellier, CNRS/IN2P3, CC 72, Place Eugfene Bataillon, 34095 Montpellier Cedex 5, France

18 DSM/Irfu, CEA Saclay, 91191 Gif-Sur-Yvette Cedex, France 19 Astronomical Observatory, The University of Warsaw, Al.

Ujazdowskie 4, 00-478 Warsaw, Poland

20 Aix Marseille University CNRS/IN2P3, CPPM UMR 7346, 13288 Marseille, France

21 Instytut Fizyki Jadrowej PAN, ul. Radzikowskiego 152, 31-342 Kraków, Poland

22 School of Physics, University of the Witwatersrand, 1 Jan Smuts Avenue, Braamfontein, Johannesburg 2050, South Africa

23 Laboratoire d’Annecy-le-Vieux de Physique des Particules, Universitó Savoie Mont-Blanc, CNRS/IN2P3, 74941 Annecy-le- Vieux, France

24 Landessternwarte, Universitat Heidelberg, Konigstuhl, 69117 Heidelberg, Germany

25 Universitó Bordeaux, CNRS/IN2P3, Centre d’Etudes Nuclćaires de Bordeaux Gradignan, 33175 Gradignan, France

26 Oskar Klein Centre, Department of Physics, Stockholm University, Albanova University Center, 10691 Stockholm, Sweden

27 Institut fur Astronomie und Astrophysik, Universitat Tubingen, Sand 1, 72076 Tubingen, Germany

28 Laboratoire Leprince-Ringuet, Ecole Polytechnique, CNRS/IN2P3, 91128 Palaiseau, France

29 APC, AstroParticule et Cosmologie, Universitó Paris Diderot, CNRS/IN2P3, CEA/Irfu, Observatoire de Paris, Sorbonne Paris Citó, 10 rue Alice Domon et Lćonie Duquet, 75205 Paris Cedex 13, France

30 Univ. Grenoble Alpes, IPAG; CNRS, IPAG, 38000 Grenoble, France

31 Department of Physics and Astronomy, The University of Leicester, University Road, Leicester, LE1 7RH, UK

32 Nicolaus Copernicus Astronomical Center, Polish Academy of Sci­

ences, ul. Bartycka 18, 00-716 Warsaw, Poland

33 Institut fur Physik und Astronomie, Universitat Potsdam, Karl-Liebknecht-Strasse 24/25, 14476 Potsdam, Germany 34 Friedrich-Alexander-Universitat Erlangen-Nurnberg, Erlangen

Centre for Astroparticle Physics, Erwin-Rommel-Str. 1, 91058 Erlangen, Germany

35 DESY, 15738 Zeuthen, Germany

36 Obserwatorium Astronomiczne, Uniwersytet Jagiellonski, ul. Orla 171, 30-244 Kraków, Poland

37 Centre for Astronomy, Faculty of Physics, Astronomy and Informat­

ics, Nicolaus Copernicus University, Grudziadzka 5, 87-100 Torun, Poland

38 Department of Physics, University of the Free State, PO Box 339, Bloemfontein 9300, South Africa

39 ITA Universitat Heidelberg, 69120 Heidelberg, Germany 40 GRAPPA, Institute of High-Energy Physics, University of

Amsterdam, Science Park 904, 1098 XH Amsterdam, The Netherlands

41 Department of Physics, Rikkyo University, 3-34-1 Nishi-Ikebukuro, Toshima-ku, Tokyo 171-8501, Japan

42 Japan Aerpspace Exploration Agency (JAXA), Institute of Space and Astronautical Science (ISAS), 3-1-1 Yoshinodai, Chuo-ku, Sagamihara, Kanagawa 229-8510, Japan

43 Now at The School of Physics, The University of New South Wales, Sydney 2052, Australia