Very-high-energy gamma-rays from the Universe's middle age: detection of the z=0.940 blazar PKS 1441+25 with MAGIC

The flat-spectrum radio quasar PKS 1441+25 at a redshift of z = 0.940 is detected between 40 and 250 GeV with a significance of 25.5 {\sigma} using the MAGIC telescopes. Together with the gravitationally lensed blazar QSO B0218+357 (z = 0.944), PKS 1441+25 is the most distant very high energy (VHE) blazar detected to date. The observations were triggered by an outburst in 2015 April seen at GeV energies with the Large Area Telescope on board Fermi. Multi-wavelength observations suggest a subdivision of the high state into two distinct flux states. In the band covered by MAGIC, the variability time scale is estimated to be 6.4 +/- 1.9 days. Modeling the broadband spectral energy distribution with an external Compton model, the location of the emitting region is understood as originating in the jet outside the broad line region (BLR) during the period of high activity, while being partially within the BLR during the period of low (typical) activity. The observed VHE spectrum during the highest activity is used to probe the extragalactic background light at an unprecedented distance scale for ground-based gamma-ray astronomy.

is used to probe the extragalactic background light at an unprecedented distance scale for ground-based gamma-ray astronomy.
Key words: cosmic background radiationgalaxies: activegalaxies: jetsgamma rays: galaxiesquasars: individual (PKS 1441+25) is a known high-energy (HE; E 0.1 GeV 100 < < GeV) γ-ray flat spectrum radio quasar (FSRQ) (Abdo et al. 2010a;Nolan et al. 2012;Ackermann et al. 2013) located at z 0.9397 0.0003 stat =  82 . In 2015 January it became active from γ-rays to the near-infrared (Carrasco et al. 2015;Ojha 2015;Pacciani 2015;Pursimo & Ojha 2015). In April, the detection of the source with a hard spectral index with the Fermi Gamma-ray Space Large Area Telescope (LAT) together with increased multi-wavelength (MWL) emission triggered the MAGIC observations. They resulted in the first detection of this source at very high energy (VHE, E 100 > GeV; Mirzoyan 2015), later followed up by VERITAS (Mukherjee 2015). This detection makes PKS1441 +25 the 83 5th FSRQ with a firm classification detected at VHE, and the most distant known VHE source, along with QSO B0218+357 (z 0.944 0.002 =  ; Sitarek et al. 2015). In this letter, the MWL observations are discussed in the context of an external Compton model for four different states of activity, dubbed periods A (MJD 57125.0-57130.0), B (57130.0-57135.5), C (57135.5-57139.5) and D (57149.0-57156.0). Upper limits on the extragalactic background light (EBL) are obtained in the VHE band.

VHE γ-Ray Observations
MAGIC is a stereoscopic system consisting of two 17 m diameter Imaging Atmospheric Cherenkov Telescopes located at the Observatorio del Roque de los Muchachos, on the Canary Island of La Palma. The current sensitivity for lowzenith observations (zd 30 < ) above 220 GeV is 0.66 0.03%  of the Crab Nebulaʼs flux in 50 hr (Aleksić et al. 2015).
The MAGIC telescopes monitored PKS1441+25 from 2015 April 18 to 27 (MJD 57130-57139, for a total of 29.9 hr) and May 8-9 (MJD 57150-57151, for 1.8 hr), the observational gap being due to the full-moon break. The observations were performed in wobble mode with a 0 . 4  offset and four symmetric positions with respect to the camera center (Fomin et al. 1994). The data were collected in the zenith angle range of zd 3 38  < < . The analysis of the data is performed using the standard MAGIC analysis framework MARS (Zanin et al. 2013;Aleksić et al. 2015) and Monte Carlo (MC) simulations matching the night-sky background levels.
PKS1441+25 is detected with a significance of 25.5σ (γ-ray-like excess N 4010 160 ex =  ) during periods B+C. No significant emission was found in period D.
The differential VHE spectrum is measured from 40 to 250 GeV and 50-160 GeV in periods B and C respectively. A power-law (PWL) can describe both observed and EBL-corrected spectra using the model of Domínguez et al. (2011): where the normalization constant f 0 , the spectral index Γ and the goodness of the fit ( ndf 2 c and p-value) are: From a likelihood ratio test (LRT), a model with intrinsic curvature such as a log-parabola (LP) is preferred at 4.2 s during period B. It is defined as:

HE γ-Ray Observations
In nominal survey mode the LAT monitors the entire γ-ray sky every 3 hr in the energy range from 20 MeV to at least 300 GeV (Atwood et al. 2009). We select Pass 8 SOURCE class events collected from 2015 April 8 to May 23 (MJD 57120-57165) between 100 MeV to 500 GeV and within a 10°R egion of Interest (ROI) centered at the location of PKS1441 +25. In order to reduce contamination from the Earth Limb, a zenith angle cut of 90 <  is applied. The analysis is performed with the ScienceTools software package version v10r0p5 using the P8R2_SOURCE_V6 84 instrument response functions and the gll_iem_v06 and iso_P8R2_SOURCE_V6_v06 models 85 for the Galactic and isotropic diffuse emission, respectively.
The likelihood fit is performed using gtlike, including all 3FGL sources (Acero et al. 2015) within 20°from PKS1441 +25. The spectral indices and fluxes are left free for sources 82 From SDSS: http://skyserver.sdss.org/dr10/en/get/SpecById.ashx? id=6780257851631206400; see also Shaw et al. (2012). 83 http://tevcat.uchicago.edu 84 http://fermi.gsfc.nasa.gov/ssc/data/analysis/documentation/Cicerone/ Cicerone_LAT_IRFs/IRF_overview.html 85 http://fermi.gsfc.nasa.gov/ssc/data/access/lat/BackgroundModels.html within 10°, while sources from 10°to 20°have their parameters fixed to the catalog value. Both the flux and the spectral index of PKS1441+25 are left free for the light curve calculation, while the parameters for the rest of the sources in the ROI are fixed except the diffuse components. Five photons of energies 10-50 GeV were detected with a probability of association with PKS1441+25 larger than 99.6%, calculated with gtsrcprob. The spectrum of PKS1441+25 is well fit by a PWL (as in the 3FGL catalog) and no significant curvature was found. During the flare (period B+C), the spectral index is 1.75 0.06 G =  , harder than the 3FGL value 3FGL G = 2.13 0.07  .

Hard X-Ray Observations
NuSTAR (Nuclear Spectroscopic Telescope Array; Harrison et al. 2013) is a hard X-ray telescope operating in the energy range between 3 and 79 keV. PKS1441+25 was observed with NuSTAR on 2015 April25-26 (MJD 57137.1113) for a total (on-source) exposure of 40 ks. The data are processed using the standard nupipeline script (version 1.4.1) available in the NuSTARDAS software package (Perri et al. 2014). The source spectrum extends up to ;25keV, and can be described by a PWL with spectral index 2.30 0.08 . No significant variability is detected during the observation. The Swift-UVOT (Roming et al. 2005) flux in several bands was estimated using aperture photometry. De-reddening is performed using

Optical Observations
Optical R-band observations started on MJD 57130 and were performed using the 35 cm Celestron telescope attached to the KVA 86 60 cm telescope (La Palma, Canary Islands, Spain) and the 50 cm Hans-Haffner-Telescope (Hettstadt, Würzburg, Germany). 87 The data are analyzed using differential photometry and corrected for Galactic extinction (Schlafly & Finkbeiner 2011). The host galaxy contribution is negligible compared to the flux of the source during these observations. The optical emission shows a high degree of polarization, reaching a maximum of 37.7% on MJD57132 (Smith & Tutar Ozdarcan 2015).

Near-infrared Observations
NIR observations in the J, H, and K S bands started on MJD 57141 with CANICA, 88 a direct camera at the 2.1 m telescope of the Guillermo Haro Observatory located at Cananea, Mexico. The flux is estimated by means of differential photometry using the 2MASS catalog (Skrutskie et al. 2006).

Radio Observations
The observations of PKS1441+25 with the Metsähovi 13.7 m radio telescope started on MJD 57135. The measurements were made with a 1 GHz-band dual beam receiver centered at 37 GHz. A detailed description of the observation and analysis methods can be found in Teräsranta et al. (1998).

Multi-wavelength Flux Evolution
The MWL light curve is presented in Figure 1  . For X-rays, a halving flux time of 7.6 1.7 days  was found. The average flux in B is larger than in C by a factor of F F 1.80 0.27 B C =  in VHE. A similar pattern was found in X-rays (F F 1.58 0.17 B C =  ), optical (F F 1.23 0.02 B C =  ) and a hint in the HE (1.40 ± 0.29). No intra-night variability is detected.

Broadband Spectral Energy Distribution
The MWL spectral energy distributions (SEDs) shown in Figure 2 indicate a shift of both synchrotron and inverse-Compton (IC) peaks to higher energies during the 2015 observation campaign with respect to the archival data, accompanied by a significant variation of the X-ray and HE γ-ray spectral indices. This behavior resembles the less extreme outburst seen in PMNJ2345-1555 (Ghisellini et al. 2013), and can be explained by a change in the emitting region location: within the broad-line region (BLR) in the quiescent state to beyond the BLR during the outburst, where the external photon field is dominated by the optical-UV from the BLR or the IR thermal emission of a dusty torus, respectively (conventional framework for γ-loud FSRQ; e.g., Ghisellini & Tavecchio 2009;Tavecchio et al. 2011).
The consequences of this scenario are two-fold: (1) since the radiation energy density of the IR component is much lower than the one associated with the optical-UV photons from the BLR, the electron radiative cooling is less effective and the energy reachable by the acceleration process could be higher, accounting for the shift of the SED peak toward higher energies; (2) given the much lower energy of the external photons, absorption of γ-rays by pair production occurs only above several hundreds of GeV (e.g., Protheroe & Biermann 1997), enabling the detection of FSRQs at VHE. For an emission region well within the BLR, strong absorption features are expected for energies above tens of GeV (see e.g., Liu & Bai 2006;Donea & Protheroe 2003), which are not observed in the 2015 MWL SEDs.
According to this framework, the proposed SED model for the 2015 observations assumes that the emission region is located at a distance d R BLR > from the central compact object. Adopting the simple scaling proposed by Ghisellini & Tavecchio (2009) =´cm. The resulting emission is calculated using the code described in Maraschi & Tavecchio (2003). The emission region is assumed to be spherical (in the source frame) with radius R, in motion with bulk Lorentz factor Γ at angle v q with respect to the line of sight. It contains a tangled magnetic field with intensity B and a population of relativistic leptons described by a smoothed broken PWL energy distribution between Lorentz factors min g and max g , with a break at b g , slopes n 1 and n 2 and normalization K estimated at γ = 1. The external photon field (dominated by the IR torus emission) is assumed to follow a blackbody spectrum with Note. The other parameters are kept fixed (see the text). The IC peak frequency (in logarithmic scale) and the Compton Dominance (CD) are also reported. ; following Ghisellini & Tavecchio 2009) and temperature T, diluted within a region of radius R IR .
The model also includes γ-ray absorption within the IR radiation field of the torus. Assuming a temperature T 10 K . Given the large optical depth and the relatively broad spectrum of the target photons, the absorption is appreciable at few hundreds of GeV, i.e., 5% at 200 GeV and 50% at 300 GeV in the observer frame. Note also that an additional softening of the spectrum can be due to the fact that the emission in the VHE band is produced by scattering occurring in the Klein-Nishina (KN) regime (e.g., Blumenthal & Gould 1970;Zdziarski & Krolik 1993;Moderski et al. 2005).
To decrease the number of free parameters, we fix the bulk Lorentz and Doppler factor to 15 G = and 20 d = , close to the average obtained for a large sample of γ-loud FSRQ (Ghisellini & Tavecchio 2015). This implies a viewing angle of the jet 2 . 7 v q =  , and the aperture angle is fixed to 0.1 rad j q = (5 . 7  ). We assume that the emission region is located beyond but not very far from the BLR, d 5 10 17 =´cm, implying R 5 10 16 =´cm. The low-energy slope n 1 is fixed to the standard value of 2. The remaining parameters are chosen to reproduce the synchrotron and IC components (see Table 1). To account for the different flux states, an evolution in both the electron distribution and the magnetic field is required. For comparison, the archival data (representation of the quiescent state) were modeled considering the emitting region partially within the BLR (standard framework) at d 1.4 10 17 =´cm, so that the gg optical depth is small as indicated by the highest energy point of the 3FGL spectrum. The ratio between the peak luminosities (Compton Dominance, CD), are reported in Table 1. During the outburst, syn n lies more than an order of magnitude outside the FSRQ parameter space in the CD sequence proposed by Finke (2013), indicating a shift in the sequence during flares. The high degree of optical polarization suggests that the emission may come from a compressed region in the jet like an internal shock, which is also an ideal site for electron acceleration/injection.

EXTRAGALACTIC BACKGROUND LIGHT CONSTRAINTS
VHE γ-rays from distant blazars can interact with the optical-UV photons from the EBL via pair production (Gould & Schreder 1967;Stecker et al. 1992), resulting in an attenuation of the intrinsic VHE spectrum. The EBL imprint in the γ-ray spectra from distant blazars can be used to constrain the EBL density.
Measurements of the EBL absorption can be derived under some assumptions on the intrinsic spectrum of the source (see, e.g., Ackermann et al. 2012;Abramowski et al. 2013). With a redshift of z = 0.940 and a strong detection in the VHE band, PKS1441+25 allows us to probe EBL models at a distance never explored before in this energy regime with ground-based gamma-ray instruments. However, KN effects together with an expected intrinsic γ-ray absorption in the VHE band (see Section 3.2), can mimic the effect of EBL absorption, making it difficult to disentangle the two effects.
A LRT was used to compare a null hypothesis (no EBL absorption) with respect to the hypothesis of EBL absorption with a scaled opacity z E , ( ) a t as in Abdo et al. (2010b). Predicted opacities from Domínguez et al. (2011, D11 t ), Franceschini et al. (2008, F08 t ), Gilmore et al. (2012, G12 t ) and Scully et al. (2014, S14 t ) are considered, while α is a free scaling parameter. Different intrinsic spectral shapes were assumed: PWL dF dE E E 10 p Note.The normalization factor 10 p 0 is given in units of erg cm s [( ) )] where E is measured in GeV and E 100 GeV 0 = . The limits are reported in Table 2 and an example is given in Figure 3. A possible overall systematic error of 15%  in the absolute energy scale of the instrument is considered. Under the assumption that no curvature is present in the intrinsic VHE spectrum, the measured spectrum is compatible with the present generation of EBL models.
The 95% confidence level limit obtained in this work for Franceschini et al. (2008) is compatible with the one found in Ackermann et al. (2012) for z 0.5 1.6   , 1.3 0.4 a =  , which is obtained from observations with a wide range of redshift values while our UL is calculated for a precise redshift value.
The estimated scaling on the optical depth can be translated into EBL density constraints as shown in Domínguez et al. (2011) and Abramowski et al. (2013). The observed VHE spectrum allow us to constrain the EBL density between 0.

CONCLUSIONS
MAGIC has detected for the first time VHE emission from the z = 0.940 blazar PKS1441+25 during a MWL outburst in 2015 April. PKS1441+25 is, together with QSOB0218+357, the most distant VHE source detected so far. This allow us to study VHE blazars when the universe was only half of its current age.
The evolution of the MWL SED is studied in the framework of an external Compton emission model. The absence of intrinsic absorption features in the HE and the VHE regime constrains the localization of the emitting region to be just outside of the BLR during this period of high activity, while it is expected to be partially compatible with the BLR during the period of low activity. The SED evolution reveals changes in the electron distribution and the magnetic field.
For the first time, the VHE measurements are used to indirectly probe the EBL at redshifts out to z 1 with groundbased gamma-ray instruments. Although an internal cutoff cannot be excluded, the measured VHE spectrum is consistent with a steepening due to attenuation caused by the EBL. Employing state-of-the-art EBL models, upper limits to the EBL density are derived. The upper limits on the opacity calculated under the assumption of an intrinsic spectrum compatible with a PWL function for different EBL models result in z E , and z E , 3.41 S12 min ( ) t t < for EBL models from Domínguez et al. (2011), Franceschini et al. (2008, Gilmore et al. (2012) and maximum and minimum from Scully et al. (2014), respectively. The dotted and dashed lines show the best-fitting PWL, respectively. The gray shaded area accounts for the uncertainties derived by the use of different EBL models (Franceschini et al. 2008;Domínguez et al. 2011;Gilmore et al. 2012). (b) The probability of fit as a function of EBL relative opacity (Domínguez et al. 2011, D11). Only period Bwas considered (without upper limits). The best fit is marked with solid vertical lines and 95% confidence level upper limits with dashed vertical lines.