VERY HIGH ENERGY γ-RAYS FROM THE UNIVERSE'S MIDDLE AGE: DETECTION OF THE z = 0.940 BLAZAR PKS 1441+25 WITH MAGIC

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 low-zenith observations (${zd}\lt 30^\circ$) above 220 GeV is $0.66\pm 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\buildrel{\circ}\over{.} 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 $3^\circ \lt {zd}\lt 38^\circ$.

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}_{\mathrm{ex}}=4010\pm 160$) 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):

$$\begin{eqnarray}&&\displaystyle \frac{{dF}}{{dE}}={f}_{0}{\left(\displaystyle \frac{E}{100\;\mathrm{GeV}}\right)}^{-{\rm{\Gamma }}},\end{eqnarray} \tag{ 1 }$$

where the normalization constant f₀, the spectral index Γ and the goodness of the fit (${\chi }^{2}/{ndf}$ and p-value) are:

• 1.  
  Period B:

  • (i)  
    Observed: ${f}_{0}=(1.14\pm {0.06}_{\mathrm{stat}}\pm {0.20}_{\mathrm{sys}})\;\times \;$ ${10}^{-9}{\mathrm{cm}}^{-2}\;{{\rm{s}}}^{-1}\;{\mathrm{TeV}}^{-1}$ and ${\rm{\Gamma }}=4.62\pm {0.11}_{\mathrm{stat}}\pm {0.18}_{\mathrm{sys}}$ (${\chi }^{2}/{ndf}=22.9/6$, $P=8.4\times {10}^{-4}$)
  • (ii)  
    EBL-corrected: ${f}_{0}=(2.7\pm {0.1}_{\mathrm{stat}}\pm {0.5}_{\mathrm{sys}})\;\times$ ${10}^{-9}{\mathrm{cm}}^{-2}\;{{\rm{s}}}^{-1}\;{\mathrm{TeV}}^{-1}$ and ${\rm{\Gamma }}=3.18\pm {0.15}_{\mathrm{stat}}\pm {0.18}_{\mathrm{sys}}$ (${\chi }^{2}/{ndf}=5.6/6$, P = 0.47)
• 2.  
  Period C:

  • (i)  
    Observed: ${f}_{0}=(0.82\pm {0.09}_{\mathrm{stat}}\pm {0.13}_{\mathrm{sys}})\;\times \;$ ${10}^{-9}{\mathrm{cm}}^{-2}\;{{\rm{s}}}^{-1}\;{\mathrm{TeV}}^{-1}$ and ${\rm{\Gamma }}=3.7\pm {0.4}_{\mathrm{stat}}\pm {0.2}_{\mathrm{sys}}$ (${\chi }^{2}/{ndf}=2.7/3$, P = 0.44)
  • (ii)  
    EBL-corrected: ${f}_{0}=(1.7\pm {0.2}_{\mathrm{stat}}\pm {0.3}_{\mathrm{sys}})\;\times \;$ ${10}^{-9}{\mathrm{cm}}^{-2}\;{{\rm{s}}}^{-1}\;{\mathrm{TeV}}^{-1}$ and ${\rm{\Gamma }}=2.5\pm {0.4}_{\mathrm{stat}}\pm {0.2}_{\mathrm{sys}}$ (${\chi }^{2}/{ndf}=4.3/3$, P = 0.23).

From a likelihood ratio test (LRT), a model with intrinsic curvature such as a log-parabola (LP) is preferred at $4.2\;\sigma$ during period B. It is defined as:

$$\begin{eqnarray}&&\displaystyle \frac{{dF}}{{dE}}={f}_{0}{\left(\displaystyle \frac{E}{100\;\mathrm{GeV}}\right)}^{-{{\rm{\Gamma }}}_{\mathrm{LP}}-b{\mathrm{log}}_{10}\frac{E}{100\;\mathrm{GeV}}},\end{eqnarray} \tag{ 2 }$$

where ${f}_{0}=(1.39\pm {0.09}_{\mathrm{stat}}\pm {0.24}_{\mathrm{sys}})\times {10}^{-9}\;{\mathrm{cm}}^{-2}\;{{\rm{s}}}^{-1}{\mathrm{TeV}}^{-1}$, ${{\rm{\Gamma }}}_{\mathrm{LP}}=4.69\pm {0.16}_{\mathrm{stat}}$ and $b=3.2\pm {1.0}_{\mathrm{stat}}$ (${\chi }^{2}/{ndf}=5.2/5$, P = 0.39). A full description of the MAGIC systematic uncertainties can be found in Aleksić et al. (2015) and references therein.

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° Region 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 $\lt 90^\circ$ is applied. The analysis is performed with the ScienceTools software package version v10r0p5 using the P8R2_SOURCE_V6 ⁸⁴ instrument response functions and the gll_iem_v06 and iso_P8R2_SOURCE_V6_v06 models⁸⁵ 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 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 ${\rm{\Gamma }}=1.75\pm 0.06$, harder than the 3FGL value ${{\rm{\Gamma }}}_{3\mathrm{FGL}}=2.13\pm 0.07$.

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 ${\rm{\Gamma }}=2.30\pm 0.08$ (${\chi }^{2}/{ndf}=10.4/7$). No significant variability is detected during the observation.

A Swift (Gehrels 2004) target of opportunity started on 2015 April 15. Swift-XRT (Burrows et al. 2005) observed the source in photon-counting mode. Standard filtering and analysis of the data were employed. The source exhibited a soft X-ray photon index (from 1.94 ± 0.16 to 2.55 ± 0.24) and is described by an absorbed PWL model, with the Galactic absorption fixed to ${N}_{H}=3.2\;\times \;$ ${10}^{20}\;{\mathrm{cm}}^{-2}$ (Kalberla et al. 2005) during April–May 2015. For comparison, the observations on 2010 June 12 (MJD 55359) can be fit with a PWL with spectral index 1.2 ± 0.3.

The Swift-UVOT (Roming et al. 2005) flux in several bands was estimated using aperture photometry. De-reddening is performed using $E(B-V)=0.033$ (Schlafly & Finkbeiner 2011) and ${A}_{V}/E(B-V)=3.1$ (Schultz & Wiemer 1975).

Optical R-band observations started on MJD 57130 and were performed using the 35 cm Celestron telescope attached to the KVA⁸⁶ 60 cm telescope (La Palma, Canary Islands, Spain) and the 50 cm Hans-Haffner-Telescope (Hettstadt, Würzburg, Germany).⁸⁷ 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).

NIR observations in the J, H, and K_S bands started on MJD 57141 with CANICA,⁸⁸ 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).

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).

The MWL light curve is presented in Figure 1. In the VHE band, the no-variability hypothesis can be discarded as it results in a ${\chi }^{2}/{ndf}=52.5/11$ (${P}_{\mathrm{const}}^{B-D}=2.2\times {10}^{-7}$) for $B+C+D$. A constant fit is also unlikely for the flare in April (B+C) with a ${\chi }^{2}/{ndf}=26.0/9$ (${P}_{\mathrm{const}}^{B+C}=2.1\times {10}^{-3}$). We gauge the characteristic variability timescale by heuristically fitting the VHE light curve with a Gaussian function. The fit provides a standard deviation $\sigma =5.5\pm 1.6\;\mathrm{days}$ (halving flux time of $6.4\pm 1.9\;\mathrm{days}$) and a peak flux of $(8.8\pm 0.6)\times {10}^{-11}\;{\mathrm{cm}}^{-2}\;{{\rm{s}}}^{-1}$ (${\chi }^{2}/{ndf}=5.3/9$, ${P}_{\mathrm{Gaussian}}^{B-D}=0.80$). For X-rays, a halving flux time of $7.6\pm 1.7\;\mathrm{days}$ was found.

Zoom In Zoom Out Reset image size

The average flux in B is larger than in C by a factor of ${F}_{B}/{F}_{C}=1.80\pm 0.27$ in VHE. A similar pattern was found in X-rays (${F}_{B}/{F}_{C}=1.58\pm 0.17$), optical (${F}_{B}/{F}_{C}=1.23\pm 0.02$) and a hint in the HE (1.40 ± 0.29). No intra-night variability is detected.

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).

Zoom In Zoom Out Reset image size

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\gt {R}_{{\rm{BLR}}}$ from the central compact object. Adopting the simple scaling proposed by Ghisellini & Tavecchio (2009), ${R}_{{\rm{BLR}}}$ depends only on the disk luminosity, ${R}_{{\rm{BLR}}}={10}^{17}{({L}_{\mathrm{disk}}/{10}^{45})}^{1/2}$ cm. The latter can be inferred from the luminosity of the optical broad lines, ${L}_{{\rm{disk}}}\simeq 2.0\times {10}^{45}$ erg s⁻¹ (Ghisellini & Tavecchio 2015), resulting in ${R}_{{\rm{BLR}}}=1.4\times {10}^{17}\ \mathrm{cm}$. In the same way, the size of the dusty torus can be inferred from a similar scaling law, ${R}_{{\rm{IR}}}=3.5\times {10}^{18}$ 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 ${\theta }_{{\rm{v}}}$ 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 ${\gamma }_{{\rm{min}}}$ and ${\gamma }_{{\rm{max}}}$, with a break at ${\gamma }_{b}$, slopes n₁ and n₂ and normalization K estimated at γ = 1. The external photon field (dominated by the IR torus emission) is assumed to follow a blackbody spectrum with luminosity ${L}_{\mathrm{IR}}=\xi {L}_{{\rm{disk}}}$ ($\xi =0.6$; 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\approx {10}^{3}\ {\rm{K}}$ for the IR torus, the maximum absorption is reached at $E={({m}_{{\rm{e}}}{c}^{2})}^{2}/2.8\;{kT}\simeq 1\ \mathrm{TeV}$ in the source frame, with an optical depth ${\tau }_{\gamma \gamma }\approx ({\sigma }_{{\rm{T}}}/5)({U}_{\mathrm{IR}}/h{\nu }_{\mathrm{IR}}){R}_{{\rm{IR}}}\;$ $\approx 250$. 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\;\mathrm{GeV}$ and $50\%$ at $300\;\mathrm{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 ${\rm{\Gamma }}=15$ and $\delta =20$, close to the average obtained for a large sample of γ-loud FSRQ (Ghisellini & Tavecchio 2015). This implies a viewing angle of the jet ${\theta }_{{\rm{v}}}=2\buildrel{\circ}\over{.} 7$, and the aperture angle is fixed to ${\theta }_{{\rm{j}}}=0.1\;{\rm{rad}}$ ($5\buildrel{\circ}\over{.} 7$).

We assume that the emission region is located beyond but not very far from the BLR, $d=5\times {10}^{17}$ cm, implying $R=5\times {10}^{16}$ cm. The low-energy slope n₁ 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\times {10}^{17}$ cm, so that the $\gamma \gamma$ 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, ${\nu }_{\mathrm{syn}}$ 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.

Table 1.  Input Parameters for the Emission Models Shown in Figure 2

Period	MJD	${\gamma }_{{\rm{min}}}$	${\gamma }_{{\rm{b}}}({10}^{4})$	${\gamma }_{{\rm{max}}}({10}^{6})$	n2	B (G)	K (${10}^{3}\;{\mathrm{cm}}^{-3}$)	${\nu }_{\mathrm{IC}}[{\rm{H}}z]$	CD
A	57125.0–57130.0	80	1.0	1.0	3.55	0.15	2.80	24.2	24
B	57130.0–57135.5	80	1.0	1.0	3.70	0.15	4.00	24.1	25
C	57135.5–57139.5	50	0.8	1.0	3.75	0.17	3.35	24.0	21
D	57149.0–57156.0	50	0.5	0.2	3.90	0.23	2.00	23.6	13
Archival	...	20	10−2	3 × 10−2	3.05	0.35	70	22.4	7

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.

Download table as:  ASCIITypeset image