Fermi-LAT high $z$ AGN and the Extragalactic Background Light

TL;DR

This paper investigates how the extragalactic background light affects the spectra of high-redshift gamma-ray sources observed by Fermi-LAT, using spectral fitting and Monte Carlo simulations to test EBL models and predict future observations with CTA.

Contribution

It introduces a method to estimate the intrinsic spectra of high-z AGN and tests EBL models against VHE observations, enhancing understanding of the EBL's impact on gamma-ray astronomy.

Findings

VHE photons from distant sources are consistent with current EBL models.

Spectral fitting below 10 GeV can estimate unabsorbed source spectra.

Upcoming CTA observations will improve EBL constraints.

Abstract

Observations of distant gamma-ray sources are hindered by the presence of the extragalactic background light (EBL). In order to understand the physical processes that result in the observed spectrum of sources, it is imperative that a good understanding of the EBL is included. In this work, an investigation into the imprint of the EBL on the observed spectra of high redshift Fermi-LAT AGN is presented. By fitting the spectrum below $\sim$10 GeV, an estimation of the un-absorbed intrinsic source spectrum is obtained; by applying this spectrum to data up to 300 GeV, it is then possible to derive a scaling factor for different EBL models. A second approach uses 5 sources (PKS 0426-380, 4C +55.17, Ton 116, PG 1246+586 and RBS 1432) which were found to exhibit very high energy emission ($E_{\gamma}>100$ GeV). Through Monte Carlo simulations it is shown that the observation of VHE photons, despite the large distances of these objects, is consistent with current EBL models. Many of these sources would be observable with the upcoming ground based observatory the Cherenkov Telescope Array (CTA), leading to a better understanding of the EBL.

Tables

Table 1. Table 1: Summary table of source data used in this study.

Source	Type	RA	Dec	$z$	$z_{err}$	$z_{ref}$
PKS 0426-380	BLL	67.18	37.93	1.111	…	Heidt et al. (2004)
4C +55.17	FSRQ	149.41	55.38	0.901	0.00019	SDSS Collaboration (2016)
Ton 116	BLL	190.80	36.48	1.066	0.00150	Abazajian et al. (2009)
PG 1246+586	BLL	192.08	58.34	0.847	0.00169	Abazajian et al. (2009)
RBS 1432	BLL	221.50	36.35	1.565	+0.275-0.125	Richards et al. (2009)
TXS 1452+516	BLL	223.61	51.41	1.522	0.00183	Alam et al. (2015)
GB6 J0043+3426	FSRQ	10.95	34.44	0.966	…	Shaw et al. (2012)
B0218+357	FSRQ	35.27	35.94	0.960	…	Lawrence (1996)
PKS 0235+164	BLL	39.66	16.62	0.940	…	Mao (2011)
MG2 J043337+2905	BLL	68.41	29.10	0.970	…	Massaro et al. (2009)
PKS 0454-234	FSRQ	74.26	-23.41	1.003	…	Stickel et al. (1989)
PKS 0537-441	BLL	84.70	-44.09	0.892	…	Peterson et al. (1976)
TXS 0628-240	BLL	97.75	-24.11	1.6	+0.10-0.05	Rau et al. (2012)
OJ 014	BLL	122.86	1.78	1.148	…	Sbarufatti et al. (2005)
PKS B1424-418	FSRQ	216.98	-42.11	1.522	0.002	White et al. (1988)
B2 2114+33	BLL	319.06	33.66	1.596	…	Shaw et al. (2013)

Table 2. Table 2: Intrinsic model parameters for each source derived from the binned likelihood analysis below E m ​ a ​ x subscript 𝐸 𝑚 𝑎 𝑥 E_{max} . The parameters correspond to those found in equation 2 . PL = power law model and LGP = log-parabola model.

	$N_0$ $\times {10}^{-12}$	$\Gamma$	($a,b$)	$E_{max}$
Source	ph cm${}^{-2}$s${}^{-1}$MeV-1	(PL)	(LGP)	(GeV)
PKS 0426-380	54.09	-1.99	1.99	9.00
			0.06
4C +55.17	10.64	-1.90	1.92	10.58
			0.06
Ton 116	1.57	-1.62	1.64	9.42
			-0.02
PG 1246+586	4.27	-1.81	1.82	11.20
			-0.02
RBS 1432	1.21	-1.69	1.64	6.96
			0.07
TXS 1452+516	3.3	-1.98	1.99	6.89
			0.02
GB6 J0043+3426	2.62	-1.90	1.91	9.94
			-0.03
B0218+357	10.25	-2.25	2.29	10.13
			0.04
PKS 0235+164	14.80	-2.08	2.12	10.17
			0.07
MG2 J043337+2905	3.26	-2.00	1.99	9.91
			0.03
PKS 0454-234	27.09	-2.12	2.19	9.65
			0.08
PKS 0537-441	23.95	-2.03	2.07	10.66
			0.05
TXS 0628-240	3.23	-1.73	1.72	8.12
			0.02
OJ 014	3.84	-2.00	2.01	8.75
			0.04
PKS B1424-418	47.58	-2.09	2.15	6.75
			0.05
B2 2114+33	2.42	-1.64	1.57	6.56
			0.10

Table 3. Table 3: Derived EBL correction factor for each of the 5 models used in this work. The ‘Best Fit’ column is based on whichever of the power law or log parabola models is a better fit for each source, while the last column uses only log parabola spectral models for every source, regardless of whether this is the best fit model.

	EBL correction Factor
Model	Best Fit	LogParabola
(Finke et al., 2010)	1.31$\pm$0.08	1.24$\pm$0.08
(Kneiske & Dole, 2010)	1.31$\pm$0.07	1.24$\pm$0.07
(Gilmore et al., 2012)	0.95$\pm$0.05	0.90$\pm$0.05
(Domínguez et al., 2011)	1.85$\pm$0.11	1.75$\pm$0.11
(Franceschini et al., 2008)	1.85$\pm$0.11	1.71$\pm$0.11

Table 4. Table 4: The highest energy photons and their associated probability of originating from the proposed source.

source	E1	P(E1)	$\sigma$(E1)	E2	P(E2)	$\sigma$(E2)
PKS 0426-380	122.0	0.99983	3.76	115.8	0.99991	3.92
4C +55.17	150.6	0.99990	3.89	94.1	0.99840	3.16
Ton 116	148.8	0.96454	2.10	132.2	0.99941	3.44
PG 1246+586	251.8	0.99983	3.76	198.1	0.99826	3.13
RBS 1432	140.4	0.99798	3.09	111.0	0.99932	3.40

Equations

$${\frac{dN}{dE}}_{obs}=e^{-\tau (E,z,n)}{\frac{dN}{dE}}_{int}$$

$$\frac{dN}{dE}=N_0{\left(\frac{E}{E_0}\right)}^{\Gamma },\frac{dN}{dE}=N_0{\left(\frac{E}{E_0}\right)}^{-(a+b\text{log}(E/E_0))}$$

$${\frac{dN}{dE}}_{obs}=e^{-\alpha \cdot \tau (E,z,n)}\cdot {\frac{dN}{dE}}_{int}.$$

$$\text{TS}=2\text{log}[{\mathcal{L}}_{\text{max}}(LogParabola)/{\mathcal{L}}_{\text{max}}(PowerLaw)].$$

$$\text{TS}=2\text{log}[\mathcal{L}(\alpha )/\mathcal{L}(\alpha =0)].$$

$$\bar{\alpha }=\int \text{pdf}(\alpha )\cdot \alpha d\alpha$$

$$Var=\int (\alpha -\bar{\alpha }{)}^2\cdot \text{pdf}(\alpha )d\alpha .$$

Full text

††pagerange: Fermi-LAT high $z$ AGN and the Extragalactic Background Light–A††pubyear: 2017

1 Introduction

A large fraction of the sources that have so far been detected in the very high energy (VHE, E${}_{\gamma}>100$ GeV) regime are blazars (Wakely & Horan, 2008). These active galactic nuclei (AGN) have relativistic jets orientated in such a way that the emission is beamed towards Earth, providing a large boost to the observed flux. Despite this, our ability to detect blazars at energies greater than some 10’s to 100’s of GeV is hindered by the existence of the Extragalactic Background Light (EBL), which attenuates gamma-rays by way of pair production with lower energy photons (Nikishov (1962), Gould & Schréder (1967b), Gould & Schréder (1967a) and Dwek & Krennrich (2013) for a recent review).

The EBL is thought to consist of UV-FIR (far infrared) radiation from stars (with UV-optical representing the main attenuation band for gamma-rays at energies below 100 GeV), either produced directly from their surfaces or via reprocessing by dust within their host galaxies, with an additional component coming from optically bright AGN (Abdo et al., 2010). The evolution of the EBL density is of considerable interest as it probes models of galaxy and star formation/evolution. However, direct measurements are difficult due to the presence of foreground zodiacal and Galactic light (Hauser & Dwek, 2001). Additionally, measurements of the EBL in the local universe ($z=0$) provide little information about the evolution of the EBL through past epochs. This evolution is especially important when considering the attenuation of gamma-rays from distant sources and its understanding therefore represents one of the major science goals of the Fermi-LAT space based instrument, which is able to observe gamma-rays from 100 MeV to above 300 GeV (Atwood et al., 2009).

For AGN at a redshift of z=1, the attenuation from the EBL quickly becomes significant above 10 GeV, therefore the Fermi-LAT instrument is well suited to observe, and place constraints on, the effect of the EBL. By modelling the un-absorbed part of a given AGN spectrum, it is possible to obtain an indication of the intrinsic source spectrum. Combining this with different EBL models, the amount of EBL absorption, and therefore density, can be inferred by evaluating the magnitude of the softening of the spectrum above $\sim$10 GeV. This method has been adopted to derive EBL limits by several authors (Abdo et al. (2010), Ackermann et al. (2012) H.E.S.S. Collaboration (2013), Biteau & Williams (2015), MAGIC Collaboration et al. (2016), Ahnen et al. (2016) and Mazin et al. (2016)).

An alternative method which can be used to set limits on the EBL is to measure the maximum observed photon energy from distant gamma-ray emitters, such as AGN or gamma-ray bursts (GRBs). The observation of high energy photons from large redshifts, where the universe is thought to be optically thick to gamma-rays, can challenge current EBL models. Currently, the most distant AGN detected in the VHE band are the gravitationally-lensed blazar B0 0218+357 ($z=0.94$) (MAGIC Collaboration et al., 2016) and the FSRQ PKS 1441+25 ($z=0.939$) (Ahnen et al., 2015), both detected with the MAGIC telescope, with the latter also detected by VERITAS (Benbow, 2015). Additionally, the second Fermi-LAT catalogue of hard sources (50 GeV - 2 TeV) contains a large sample of 271 AGN with a redshift range out to 2.1 (Ackermann et al., 2016). The two GRBs, GRB 090902B ($z=1.822$) (Abdo et al., 2009) and GRB 08916C ($z=4.35$) (Abdo et al., 2009), have also aided in providing limits to the EBL density.

In an attempt to increase the sample of VHE AGN, previous work (Armstrong et al. (2015), Armstrong et al. (2016)) has focused on developing clustering methods to identify extragalactic sources with $E_{\gamma}>100$ GeV in the Fermi-LAT data set. The chosen algorithm, dbscan, was first proposed in Ester et al. (1996) and was designed to efficiently detect clusters of arbitrary shape in noisy data sets. The algorithm works by evaluating the number of events around a detected photon within the data set; if these exceed a given number then they are added to a cluster and in turn evaluated, accumulating further events to the cluster that satisfy the same condition111For the dbscan sources presented in this work that were found in Armstrong et al. (2016), the search radius around a given point was 0.4 degrees, the PSF at 100 GeV, and the number of required events started at 2 and scaled upward with the integrated Galactic background emission. The dbscan algorithm makes no assumption regarding the underlying source spectrum. The ability to identify clusters within sparse noisy data sets made dbscan an ideal choice for analysing Fermi-LAT VHE data. So far this method has proven to be successful, discovering emission with $E_{\gamma}>100$ GeV from 45 sources. The sources detected are all either AGN or unassociated sources. One of the aims of this work is to draw attention to six high redshift AGN, PKS 0426-380 ($z=1.11$), 4C +55.17 ($z=0.899$), Ton 116 ($z=1.088$), PG 1246+586 ($z=0.857$), RBS 1432 ($z=1.508$) and TXS 1452+516222The Blazar TXS 1452+516 ($z=1.522$) was found in a separate work [in prep], using the same method but for $50~{}\textrm{GeV}>E_{\gamma}>2$ TeV., that were found using this method and could have implications for EBL limits.

This paper is organised as follows, in Section 2 a brief overview of the EBL and EBL models used in this work will be given. In Section 3 a description of each source of interest, with a focus on the reliability of their redshift determination, will be presented. Section 4 will describe the data selection process and Sections 5 and 6 will discuss the spectral and ‘highest energy photon’ analyses respectively. Finally Section 7 will conclude and summarise the results presented.

2 EBL Models

Since direct measurements of the EBL are hindered due to the presence of foreground emission, and because this reveals little about the evolution on cosmological scales, it is necessary to use complex models to estimate the density of EBL photons and therefore the absorption of gamma-rays. There is currently a wide range of models available which adopt different methods to determine the evolution of the EBL density as a function of redshift.

In this work, five models have been considered: Gilmore et al. (2012) (GIL12) which uses a semi-analytic model of the star formation rate, initial mass function and dust extinction, Finke et al. (2010) (FIN10) and Kneiske & Dole (2010) (K&D10) which use forward evolution models based on observations of Spitzer, ISO, Hubble, COBE, BLAST and GALEX from which the cosmic star formation rate is inferred, and finally Franceschini et al. (2008) (FRA08) and Domínguez et al. (2011) (DOM11) which use backward evolution models to model the redshift evolution of the luminosity function of galaxies based on number counts. For an overview of the different model types, see Dwek & Krennrich (2013).

The absorption of distant gamma-rays depends on the optical depth of the EBL, expressed as $\tau(E,z,n)$, which is dependent on the gamma-ray energy ($E$), the distance ($z$) and the density ($n$) of EBL photons, the last being defined by the choice of model. The total is based on an integral along the line of sight to the target source. For an individual source, the spectrum is therefore attenuated as a function of its energy such that

[TABLE]

where $dN/dE_{obs}$ is the observed spectrum and $dN/dE_{int}$ is the intrinsic/unabsorbed spectrum. For the models used in this work, tabulated data of $\tau$ as a function of energy, available for a range of redshift values, were obtained from online resources333http://www.phy.ohiou.edu/~finke/EBL/ for Finke et al. (2010), http://www.ias.u-psud.fr/irgalaxies/cib.php for Kneiske & Dole (2010), http://physics.ucsc.edu/~joel/EBLdata-Gilmore2012 for Gilmore et al. (2012) and http://www.astro.unipd.it/background/ for Franceschini et al. (2008). Data from Dominguez were obtained from https://github.com/me-manu/eblstud/blob/master/ebl/ebl_model_files/tau_dominguez10.dat due to the inactivity of the original link in Domínguez et al. (2011).. By interpolating these data sets, a function for determining $\tau$ based on the energy and redshift is obtained.

3 Source Information

An overview of each source analysed in this work is given here, with a focus on determining the most reliable redshift from literature. The sources are divided into two categories: 1) DBSCAN VHE sources, which were identified in Armstrong et al. (2016) to have significant emission in the 100 GeV - 3 TeV energy range and 2) a sample of 10 2FHL objects with the highest quoted redshifts within the catalogue which also coincided with the redshift range of the first sample. Two sources where however excluded from the selection, MG4 J000800+4712 which is quoted to have a redshift of 2.1 but is found with z=0.28 in all other literature, and PKS 0823-223 due to its proximity to the Galactic plane. The second category of sources was included in order to obtain a more reliable limit for the EBL when considering the fit to the intrinsic and absorbed spectra (see Section 5).

3.1 DBSCAN VHE Sources

1. PKS 0426-380:

Classed as a BL Lac in the 3FGL and situated in the southern hemisphere at RA=67.18*∘, Dec=-37.93∘*. The 2FHL redshift is quoted as $z=1.111$ as found in Heidt et al. (2004), where the existence of a closer source at $z=0.559$ implies micro-lensing may be present. This source has previously been noted for producing VHE photons in Tanaka et al. (2013) and Neronov et al. (2015), which both use Pass 7 data. The new Pass 8 data have reclassified the energy of these photons. 2. 4C +55.17:

Defined as a FSRQ, and located in the northern hemisphere at RA=149.41*∘, Dec=55.38∘*, in the 2FHL the redshift is quoted as 0.899, most likely obtained from Hewett & Wild (2010). However, more recent measurements from the 13th SDSS data release (SDSS Collaboration, 2016) finds a redshift of 0.901$\pm$0.00019 which will be adopted in this work. This source was considered in Neronov et al. (2015) where a redshift of 0.8955 was adopted. 3. Ton 116:

A BL Lac located in the northern hemisphere at RA=190.80*∘, Dec=36.46∘* with a redshift of 0.0 in the 2FHL indicating that no redshift was available at the time or that it was determined to be unreliable. However, measurements from SDSS indicate a redshift of 1.066$\pm$0.00150 (data release 7, (Abazajian et al., 2009)) or 1.182$\pm$0.00132 (data release 13 (SDSS Collaboration, 2016)); both have data flags associated with the spectrum. This source was also considered in Neronov et al. (2015) in which a redshift of 1.065 was used. Our work adopts the SDSS data release 7 value, which is the same within errors. 4. PG 1246+586

A BL Lac located in the northern hemisphere at RA=192.08*∘, Dec=58.34∘* and with an unknown or uncertain redshift in the 2FHL. Measurements from SDSS (data release 7 (Abazajian et al., 2009)) show a redshift of 0.847$\pm$0.00169. This source was considered in Neronov et al. (2015) who applied the same redshift as used here. 5. RBS 1432:

A northern hemisphere BL Lac at RA=221.50*∘, Dec=36.35∘* for which the 2FHL again defines the redshift as unknown or uncertain. The redshift in the literature is generally accepted as 1.565 which originates from Richards et al. (2009), although the redshift is quoted as being between 1.440 and 1.840. 6. TXS 1452+516:

This source was found in a comparative study of DBSCAN and the 2FHL [paper in prep], using the same energy range and observational period. As it was not found to be significant in VHE range (E${}_{\gamma}>$100 GeV) range, as with the rest of this sample, we do not claim in to be a VHE source. However as it is not in the 2FHL (but is in the 3FGL) it is included in this list. It is a BL Lac located in the northern hemisphere at RA=223.61*∘, Dec=51.41∘* and has a redshift of 1.522$\pm$0.00183, which originates from the 12th SDSS data release (Alam et al., 2015).

3.2 High-$z$ 2FHL Sample

1. GB6 J0043+3426:

A northern hemisphere FSRQ at RA=10.95*∘, Dec=34.44∘* with a redshift of 0.966 based on observations at the W. M. Keck Observatory and presented in Shaw et al. (2012). 2. B0218+357:

A northern hemisphere FSRQ at RA=35.27*∘, Dec=35.94∘* which is thought to be a gravitationally lensed object with two images of the FSRQ at a redshift of 0.960 (Lawrence, 1996) and a lensing source at 0.685 (Browne et al., 1993). This FSRQ has been observed with MAGIC (MAGIC Collaboration et al., 2016) and was used to determine an EBL correction factor. 3. PKS 0235+164:

A northern hemisphere BL Lac at RA=39.66*∘, Dec=16.62∘* which has an associated redshift of 0.940 as stated in Mao (2011). 4. MG2 J043337+2905:

A northern hemisphere BL Lac at RA=68.41*∘, Dec=29.10∘* with a redshift of 0.970 found in Massaro et al. (2009). 5. PKS 0454-234:

A southern hemisphere FSRQ at RA=74.26*∘, Dec=-23.41∘*. The redshift is given as 1.003 based on absorption lines in Stickel et al. (1989). 6. PKS 0537-441:

A southern hemisphere BL Lac at RA=84.71*∘, Dec=-44.09∘* with a redshift of 0.892 (Peterson et al., 1976). 7. TXS 0628-240:

A southern hemisphere BL Lac at RA=97.75*∘, Dec=-24.11∘* which has an associated photometric redshift of 1.6${}^{+0.10}_{-0.05}$ (Rau et al., 2012) or limits of 1.239 to 1.91 based on spectroscopic measurements (Shaw et al., 2013). In this work we use the photometric redshift. 8. OJ 014:

A southern hemisphere BL Lac at RA=122.86*∘, Dec=1.78∘* with a redshift of 1.148 from Sbarufatti et al. (2005). 9. PKS B1424-418:

A southern hemisphere FSRQ at RA=216.98*∘, Dec=-42.11∘* with redshift 1.522$\pm$0.002 from White et al. (1988). 10. B2 2114+33:

A southern hemisphere BL Lac located at RA=319.06*∘, Dec=33.66∘*. A redshift is given as 1.596 based on significant broad emission feature identified with C${}_{\textrm{IV}}$ consistent with a weak bump in the far blue at Ly$\alpha$. It is suggested that if the emission feature is false that this would then at least represent a lower limit (Shaw et al., 2013).

A summary of all these sources and their basic data can be found in Table 1.

4 Data Selection

Each source presented in this work was evaluated using the Pass 8 processed data (Atwood et al., 2013), which provides several improvements over the previous Pass 7 reprocessed data set (Bregeon et al., 2013)444Pass 7 was previously used to evaluate some of the sources presented in this work (Tanaka et al. (2013) and Neronov et al. (2015)).. These include better energy and angular resolution, increased effective area, an extended observable energy range (around 10 MeV < $E_{\gamma}$ <3 TeV) and better background characterisation, resulting in improved point source sensitivity. As a result of this re-evaluation of the data, the energies of many events have been updated.

To investigate each source, 8 years of data ranging from August 2008 to September 2016 (MET 239557417 - 494930839) and with energies from 100 MeV to 3 TeV were downloaded from the Fermi data servers555http://fermi.gsfc.nasa.gov/ssc/data/. All source class events were retained for both front and back converting photons. Additionally, the recommended filter expression ‘DATA_QUAL$>$0 && LAT_CONFIG==1’ was applied in order to remove any sub-optimal data affected by spacecraft events. Finally a zenith cut of 90*∘* was implemented in order to remove any $\gamma$-rays induced by cosmic-ray interactions in the Earth’s atmosphere.

5 Spectral Analysis

To investigate the effect of the EBL absorption on the spectrum of each source presented in Section 3, a binned likelihood analysis was performed. An initial fit was obtained for each source below a threshold energy $E_{max}$, here defined to be the point where there is less than 0.1% of emission attenuated by the EBL, in order to represent the intrinsic spectrum. For a source at redshift $\sim$1 this corresponds to an energy of $\sim$10 GeV. The values for each source can be found in Table 2.

From each source data set comprising a 10*∘* region of interest (ROI) around the source, an initial model is constructed from the 3FGL which includes sources out to a further 10*∘* fixed to the 3FGL values. The normalisations and spectral indices of the sources within the ROI were allowed to vary. Also included were the extragalactic diffuse emission model (iso_P8R2_SOURCE_V6_v06.txt) with a free normalisation, and the Galactic diffuse template (gll_iem_v06.fits), which was multiplied by a power law in energy and the normalisation of which was free to vary. The sources of interest were modelled with both a power law and a log parabola in the form

[TABLE]

where $N_{0}$ is the normalisation, $\Gamma$ is the spectral index of the power law, $a$ is the log parabola index, $b$ is the curvature and $E_{0}$ is a scaling factor. In the analysis of each source, all parameters excluding $E_{0}$ were left free to vary and a full binned likelihood analysis was performed for each ROI, returning the best fit model. The best fit intrinsic model parameters can be found in Table 2.

In order to investigate the level of EBL absorption present, the initial fit was scaled with $e^{-\alpha\cdot\tau(E,z,n)}$, as in Equation 1, where the scaling factor $\alpha$ has now been introduced (here a value of $\alpha=1$ would return the original EBL model absorption level). By taking the scaled spectral model, and fitting this to to data between 100 MeV and 300 GeV, which includes both the absorbed and unabsorbed sections of the spectrum, it is possible to scan through values of $\alpha$ and identify the optimal level of EBL absorption (as was performed in Ackermann et al. (2012) and H.E.S.S. Collaboration (2013); a more detailed description of the process is given in Appendix A).

For each EBL scaling factor, a likelihood was generated for the attenuated power law and log parabola models. The latter was chosen if the Test Statistic (TS) between the two models, with the maximum likelihood value of $\alpha$ applied, was greater than 16 (the power law is the null hypothesis, i.e. TS=2log$[\mathcal{L}_{\textrm{max}}(LogParabola)/\mathcal{L}_{\textrm{max}}(PowerLaw)]$)666This is the value used by the Fermi-LAT collaboration in the 3FGL to chose a given spectral model over a power law Acero et al. (2015).. In Figures 1 and 2 the delta logLikelihood can be seen for each energy bin of the full spectral energy distribution along with the intrinsic spectrum (black dotted line) which is calculated based on events with energies below that indicated by the vertical dotted line. The black dashed line shows the absorbed spectrum calculated using the chosen EBL model (in this case GIL12) and the solid black line the resulting the best maximum likelihood fit when GIL12 is modified by a scaling factor $\alpha$. The individual EBL scaling factors for each source are given in the Figures.

From the likelihood distribution for each source as a function of the EBL correction factor, the TS was calculated with respect to the null hypothesis that there is no absorption from the EBL (i.e. TS=2log$[\mathcal{L}(\alpha)/\mathcal{L}(\alpha=0)]$). These distributions were then summed in order to obtain a combined TS distribution, allowing an overall EBL correction to be calculated. The result of this process applied to the GIL12 model can be seen in Figure 3 from which an optimal correction factor of $\alpha_{\textrm{GIL12}}=0.95\pm 0.05$ was derived (1 standard deviation uncertainty). The EBL normalisation for the other models gives $\alpha_{\textrm{K\&D10}}=1.31\pm 0.07$, $\alpha_{\textrm{FIN10}}=1.31\pm 0.08$, $\alpha_{\textrm{DOM11}}=1.85\pm 0.11$ and $\alpha_{\textrm{FRA08}}=1.85\pm 0.11$ for the redshift range $0.897<z<1.596$. It is clear that the GIL12 model is consistent with the results, while it is suggested that the other models underestimate the EBL density.

As noted above, the quoted errors are statistical only. One of the main sources of systematic error in this analysis comes from the choice of spectral model, where it is clear that an incorrect choice will strongly affect the derived EBL scaling factor. This was addressed by Ackermann et al. (2012) who concluded that the assumption that spectral cutoffs existed in the EBL absorption energy range was acceptable. However, Ackermann et al. (2012) assumed each source to be best represented by a log parabola. We have therefore repeated our analysis using only log parabola models, which provides a lower limit for the EBL correction factor. This leads to the following derived scaling factors: $\alpha_{\textrm{GIL12}}=0.90\pm 0.05$, $\alpha_{\textrm{K\&D10}}=1.24\pm 0.07$, $\alpha_{\textrm{FIN10}}=1.24\pm 0.08$, $\alpha_{\textrm{FRA08}}=1.71\pm 0.11$ and $\alpha_{\textrm{DOM11}}=1.75\pm 0.11$. The derived EBL scaling factors are summarised in Table 3.

In each panel of Figures 1 and 2, the expected sensitivity for 50 hours of observation of the two arrays of the future ground based gamma-ray observatory, CTA, are shown (publicly available at https://www.cta-observatory.org/science/cta-performance/, see Bernlöhr et al. (2013) for description of method). Since CTA will possess a much improved sensitivity compared to current telescopes, it will be well placed to develop further our understanding of the EBL by observing AGN in the range $0<z<1$ and at higher energies than is possible with Fermi. This work shows that several of our sources should be detectable with 50 hours of observation, while others may also be detectable with deeper observations. Further observations with CTA would aid in reducing systematics in the choice of best fit spectral model.

6 VHE Photons

An alternative method for probing the EBL is to identify the highest energy photons associated with distant AGN. Here, detection of a photon from a high-$z$ source at energies where the density of the EBL is thought to attenuate the signal almost entirely can bring EBL models into question. In this section, we investigate whether this is the case for the sources presented in the DBSCAN VHE sample which have already exhibited significant emission above 100 GeV (which therefore excludes TXS 1452+516), despite their large redshifts.

For this analysis, it is particularly important to ensure that the VHE photons are genuinely associated with the AGN. To verify this, a subset of data was selected for each source covering energies above 100 GeV only for a 1*∘* region around the source position (much larger than the PSF at these energies). The Fermi science tool gtsrcprob was applied to this dataset in order to calculate the probability that each photon is associated with a source in or around the ROI. This is based on the likelihood of each photon and convolves a source model with the instrument response functions. For the source model, the best fit derived in the previous section without EBL absorption was used. Additionally, in order to account for diffuse components such as the Galactic diffuse model, the Fermi science tool gtdiffrsp was applied in order to add a diffuse response to the input data. The energies of the highest energy photons and the probabilities of their being associated with the source in question can be seen in Table 4. It is worth noting that, if the attenuation from the EBL is included, these probabilities will reduce. For example, the highest energy photon probability associated with PKS 0426-380 drops from 0.99983 to 0.99899 when the GIL12 model is included. For out purpouse we assume this to be a negligible difference.

As an initial consideration, the maximum energies found in Table 4 and the redshift of the proposed sources were compared to objects in the rest of the 2FHL (where the 2FHL catalogue contains the maximum energy detected). This can be seen in Figure 4 in which different levels of optical depth, starting at $\tau=1$ (where the universe becomes optically thick to gamma-rays) have also been shown. What is immediately apparent is that the majority of 2FHL sources lie within, or close to, the optically thin ($\tau<1$) region, while the sources found in the DBSCAN VHE source sample are pushing out to larger values of $\tau$.

Of the VHE sources, the largest derived optical depth corresponds to $\tau_{GIL12}=4.39$ for RBS 1432, indicating that a photon was detected despite 98.8% of the flux being attenuated. In order to evaluate the probability of observing the highest energy photons from each source, Monte Carlo simulations were performed as in Abdo et al. (2010). Here, spectral parameters were drawn from a distribution based on the best fit intrinsic model and corresponding errors from section 5, including the spectral EBL absorption from the GIL12 model. Then using the Fermi science tool gtobsim, eight year observations between 10 GeV and 1 TeV were simulated, taking into account the instrument response functions and spacecraft pointing history. Using 1000 simulations for each source, the number of observations of photons of a given energy was determined. An example of this for PG 1246+586 can be seen in Figure 5.

From these simulations, it is evident that, in the 8 years of Fermi-LAT observations, the highest energy photons detected are within expectation, with all but the 251 GeV photon from PG 1246+586 at 2.2$\sigma$ found at an energy less than one standard deviation from the simulated mean value. This result confirms that found in Section 5, that the observations presented in this work strengthen our confidence in the EBL models.

Given that it has been shown that photons with energies greater than 100 GeV are expected from these sources, this also provides further evidence that they should be observable with ground based gamma-ray observatories which, with their much greater sensitivity, could provide stricter limits on EBL models. An important caveat here is that the analysis carried out in this work makes use of the near continuous Fermi-LAT observations over 8 years. The assumption is therefore that these sources are non-variable over that time scale, or more importantly that the spectrum remains unchanged. The published variability index from the 3FGL suggests that this assumption is reasonable for all but PKS 0426-380 and RBS 1432, with the former exhibiting strong evidence for variability. Additionally, many of the sources from Section3.2 also show evidence of variability. It is expected that this may bias the results presented in this and the previous section. In a following paper, the effect of variability on the estimate of the EBL will be studied in detail.

7 Discussion and Conclusion

Fermi-LAT’s excellent energy coverage enables it to be used to observe attenuation effects from the EBL. In this work, analysis of 16 high redshift sources has been presented, in which the EBL absorption is apparent. By modelling the energy spectrum from 100 MeV to around $\sim$10 GeV (redshift dependent), over which the effect of the EBL is negligible, it is possible to obtain an idea of the intrinsic spectrum of the source of interest. Using this, the full data set from 100 MeV to 300 GeV and the absorption from a selection of EBL models, scaled by a correction factor $\alpha$, a measurement of the EBL scale was derived. By using a combination of power law and log parabola models, a best fit to the combined TS of each individual source revealed scaling factors of $\alpha_{\textrm{GIL12}}=0.95\pm 0.05$, $\alpha_{\textrm{K\&D10}}=1.31\pm 0.07$, $\alpha_{\textrm{FIN10}}=1.31\pm 0.08$, $\alpha_{\textrm{FRA08}}=1.85\pm 0.11$ and $\alpha_{\textrm{DOM11}}=1.85\pm 0.11$ for the redshift range $0.897<z<1.596$. All models apart from GIL12 predict an EBL density less than suggested by this study. As the choice of spectral model introduces an unknown uncertainty, we also derived a conservative EBL scaling factor based solely on the log parabola model, From this it was found that $\alpha_{\textrm{GIL12}}=0.90\pm 0.05$ $\alpha_{\textrm{K\&D10}}=1.24\pm 0.07$, $\alpha_{\textrm{FIN10}}=1.24\pm 0.08$, $\alpha_{\textrm{FRA08}}=1.71\pm 0.11$ and $\alpha_{\textrm{DOM11}}=1.75\pm 0.11$.

As mentioned previously, the source B0 218+357 has been observed using the ground-based MAGIC telescopes (MAGIC Collaboration et al., 2016). They used a combination of MAGIC and Fermi-LAT data for a period of activity between the 11th of July 2014 and the 6th of August 2014, and modelling B0 218+357 with a power law, to obtain EBL correction factors for the following models777In the analysis of B0 218+357, a larger range of models were investigated. We present only their results for the models also considered in this work.: $\alpha_{\textrm{FIN10}}=0.91\pm 0.32_{stat}\pm 0.19_{sys}$, $\alpha_{\textrm{FRA08}}=1.19\pm 0.42_{stat}\pm 0.25_{sys}$, $\alpha_{\textrm{DOM11}}=1.19\pm 0.42_{stat}\pm 0.25_{sys}$ and $\alpha_{\textrm{{GIL12}}}=0.99\pm 0.34_{stat}~{}^{+0.15sys}_{-0.18sys}$.

In addition to B0 218+357, the MAGIC collaboration used 11.8 hours of observations from a 2014 flare of 1ES 1011+496 ($z=0.212$) (Ahnen et al., 2016) to obtain a scaling factor for the EBL models $\alpha_{\textrm{DOM11}}=1.07~{}(-0.20,+0.24)_{stat+sys}$ and $\alpha_{\textrm{FRA08}}=1.14~{}(-0.14,+0.09)_{stat}$. This was followed by a study using 8 high frequency peaked BL Lac and 4 FSRQs within the redshift range of $0.031<z<0.944$ (which also includes the observations of the two aforementioned sources) where a scaling factor of $\alpha_{\textrm{DOM11}}=0.99~{}(-0.56,+0.15)_{stat+sys}$ was derived (Mazin et al., 2016)

Correction factors for the EBL have also been derived for less distant AGN by the HESS collaboration (H.E.S.S. Collaboration, 2013). Here data for 7 bright Blazars (Mrk 421 [z=0.031], PKS 2005-489 [z=0.071], PKS 2155-304 [z=0.116], 1ES 0229+200 [z=0.14], H 2356-309 [z=0.165], 1ES 1101-232 [z=0.186], and 1ES 0347-121 [z=0.188]), over different periods of activity, were modelled in a similar way to that presented here. By fitting a range of different spectral models to the observed data, scanning the EBL correction parameter space, a value was obtained that maximises the likelihood. In the paper, only the Franceschini et al. (2008) model was considered, for which a correction factor of $\alpha_{\textrm{FRA08}}=1.27_{-0.15_{stat}}^{+0.18_{stat}}\pm 0.24_{sys}$ was obtained.

The most encompassing study using ground based observations was performed in Biteau & Williams (2015). Here spectra from 38 sources in the redshift range $0.019<z<0.64$ were evaluated, using data from HESS, MAGIC, VERITAS, Whipple, ARGO-YBJ, TACTIC, HEGRA, Tibet and CAT. Using these data sets, along with local constraints taken from Dwek & Krennrich (2013), the spectra were fitted with power law, log parabola and exponential cutoff powerlaw models, resulting in a derived EBL spectrum. Also presented were correction factors for current EBL models: $\alpha_{\textrm{GIL12}}=1.13\pm 0.07$ $\alpha_{\textrm{K\&D10}}=1.52\pm 0.14$, $\alpha_{\textrm{FIN10}}=1.48\pm 0.07$, $\alpha_{\textrm{FRA08}}=1.05\pm 0.07$ and $\alpha_{\textrm{DOM11}}=1.16\pm 0.05$.

One of the main scientific objectives of the Fermi-LAT instrument was to measure the effect of the EBL attenuation of high energy gamma-rays. In Ackermann et al. (2012), data from the first 46 months of observation were used to model the EBL attenuation of 150 BL Lac type blazars. The method to determine the scaling of the EBL was the same as presented in this work, first modelling the intrinsic blazar spectrum and then fitting the higher energies, up to 500 GeV. In doing this, scaling factors for a large range of EBL models were determined. For the models used in this work, the corresponding weighted averages were found to be $\alpha_{\textrm{GIL12}}=0.67\pm 0.14$ $\alpha_{\textrm{K\&D10}}=0.90\pm 0.19$, $\alpha_{\textrm{FIN10}}=0.86\pm 0.23$, $\alpha_{\textrm{FRA08}}=1.02\pm 0.23$ and $\alpha_{\textrm{DOM11}}=1.02\pm 0.23$. Also presented was the binned results for the Franceschini et al. (2008) model, giving $\alpha_{\textrm{FRA08}}=1.18^{+0.94}_{-0.81}$ ($z<0.2$), $\alpha_{\textrm{FRA08}}=0.82^{+0.41}_{-0.30}$ ($0.2<z<0.5$) and $\alpha_{\textrm{FRA08}}=1.29^{+0.43}_{-0.36}$ ($0.5<z<1.6$).

The results from MAGIC Collaboration et al. (2016), Ahnen et al. (2016), Mazin et al. (2016), H.E.S.S. Collaboration (2013), Ackermann et al. (2012) and Biteau & Williams (2015) are shown in comparison to the results derived in this paper in Figure 6. It can be seen that, for all but the GIL12, our results imply a larger density of the EBL at high redshifts. For all models we find a larger estimate than found by the Fermi-LAT in Ackermann et al. (2012); however they are consistent with the results found by the MAGIC collaboration (MAGIC Collaboration et al., 2016). When considering the binned Fermi-LAT results for the FRA08 model, the last bin which encompasses the redshift range under investigation in this paper, is consistent with our results for log parabola models. When considering the EBL model types, it is suggested that the backward evolution models DOM11 and FRA08, which are based on the luminosity function of galaxies based on number counts, underestimate the EBL at large redshifts. This is not an unexpected result as these models match observed lower limits which are based on incomplete surveys. In a following paper we will attempt to investigate these trends in more detail, taking into account variability and source class.

The second part of this work focused on investigating the highest energy photons from the VHE sources found with DBSCAN. By considering the often used Fazio-Stecker relation, it seemed that these sources were pushing out to large optical depths. Upon further investigation using Monte Carlo simulations, it was shown that the observed photon energies were not unexpected (less than 1 standard deviation from the mean energy or $2.2\sigma$ for PG 1246+586), illustrating the pitfalls of relying solely on the Fazio-Stecker relation.

Given the expected performance of the future ground based observatory, CTA, which will be able to place even stronger constraints on the EBL, it was shown that several of the sources presented in the work should be observable by CTA (See Figure 1 and 2). Indeed, the observation of VHE emission from PKS 0426-380, 4C +55.17, Ton 116, PG 1246+586 and RBS 1432 is expected, providing further evidence that these sources should be observable with CTA.

Acknowledgments

TPA would like to acknowledge the support from the UK Science and Technology Facilities Council grant ST/K501979/1. AMB would like to acknowledge the financial support of Durham University. This work has made use of publicly available Fermi-LAT data from the High Energy Astrophysics Science Archive Research Center (HEASARC), provided by NASAs Goddard Space Flight Center. Finally, we would like to thank the reviewer for their insightful comments which helped improve the quality of this work.

Appendix A Detailed Source Analysis

To provide further clarity of the method by which the EBL scaling factor was determined, a more detailed description is given here. For simplicity, the source Ton 116 (z=1.066) is used as an example throughout.

For each source, an initial fit was obtained from 100 MeV to an energy where the EBL absorption is less than 0.1%. For Ton 116, at a redshift of $z=1.066$, this corresponds to an energy of 9.24 GeV using the GIL12 model. The source of interest in the model file is then replaced with a file function based on the intrinsic spectral model convolved with the EBL absorption including the scaling factor $\alpha$:

[TABLE]

A series of model files were then created using the scaled spectral model for $\alpha$ values in the range 0 to 3 at 0.01 intervals. The sources within the ROI and the background models are fixed to the initial fit results and a second fit is applied, returning the log likelihood for the entire ROI. In the left panel of Figure 7 the intrinsic power law and log parabola spectral models are compared to both the data and the absorbed spectral models ($\alpha=1$ for the GIL12 model), in the right panel the resulting log likelihood for each $\alpha$ value is shown. The standard method for choosing between models, as specified by the Fermi collaboration, is to calculate the Test Statistic (TS) that the proposed model is preferred over the power law model, where the TS is defined as

[TABLE]

For Ton 116, the TS value calculated is 1.24 and therefore the default power law model is retained (a TS$>$16 is required to justify the choice of an alternative model). Following this, the likelihood distribution is converted into a TS distribution, where this is now the TS of a given $\alpha$ value over the case where there is no EBL absorption ($\alpha=0$), i.e.

[TABLE]

Due to the asymmetrical nature of the likelihood/TS distribution, the mean and variance were calculated as in D’Agostini (2004), where the pdf is derived from the TS distribution by calculating $e^{TS(\alpha)}$ which is normalised based on the total area. From this the expected value is derived as

[TABLE]

and the variance as

[TABLE]

For Ton 116, this corresponds to $\alpha=0.90\pm 0.16$. Lastly, in order to determine an overall EBL scaling factor, the TS distributions for each source are summed to obtain an combined TS where the expected value and variance are calculated as before (See Figure 3).