GAMMA-RAYS FROM THE QUASAR PKS 1441+25: STORY OF AN ESCAPE

The extragalactic gamma-ray sky is dominated by the emission of blazars, active galactic nuclei whose jets are pointed within a few degrees of Earth. About sixty blazars have been detected at very high energy (VHE; E > 100 GeV),⁴⁵ with only four belonging to the class of flat-spectrum radio quasars (FSRQs): 3C 279 (z = 0.536, MAGIC Collaboration et al. 2008), PKS 1510–089 (z = 0.361, H.E.S.S. Collaboration et al. 2013a), PKS 1222+216 (z = 0.432, Aleksić et al. 2011), and S3 0218+35 (z = 0.944, Sitarek et al. 2015).

FSRQs are believed to host radiatively efficient disks that enrich the environment of the supermassive black hole with ultraviolet-to-optical photons. This photon field, the reprocessed emission from the clouds of the broad line region (BLR), and the infrared radiation from the "dusty torus" can all interact with gamma-rays through pair production, preventing the escape of VHE radiation from the base of the jet (Donea & Protheroe 2003).

VHE gamma-rays that do escape will face pair production on the extragalactic background light (EBL), which encompasses the ultraviolet-to-infrared emission of all stars and galaxies since the epoch of reionization (z ≲ 10).⁴⁶ Direct measurements of the EBL are prone to contamination from the bright local environment, while strict lower limits are derived from galaxy surveys, measuring the light emitted by known populations of sources (Madau & Pozzetti 2000).

VHE detections of high-redshift FSRQs constrain both the EBL and the gamma-ray emission regions in blazars. This letter reports the detection of VHE gamma-rays from the FSRQ PKS 1441+25 (z = 0.939, Shaw et al. 2012) by VERITAS. This source was observed in 2015 April and May following the VHE discovery by MAGIC (Mirzoyan 2015), triggered by multiwavelength activity (Pacciani 2015) and a spectral hardening at high energies (HE, 100 MeV < E < 100 GeV; D. Thompson 2015, private communication on behalf of the Fermi-LAT team).

PKS 1441+25 was detected from 2015 April 21 (MJD 57133) to April 28 (MJD 57140) with VERITAS, an array of four imaging atmospheric Cherenkov telescopes located in southern Arizona (Holder 2011). VERITAS imaged gamma-ray induced showers from the source above 80 GeV, enabling the detection of PKS 1441+25 (VER J1443+250) at a position consistent with its radio location and at a significance of 7.7 standard deviations (σ) during the 15.0 hr exposure (2710 ON-source events, 13780 OFF-source events, OFF normalization of 1/6). Using a standard analysis with cuts optimized for low-energy showers (Archambault et al. 2014, and references therein), we measure an average flux of Φ(>80 GeV) = (5.0 ± 0.7) × 10⁻¹¹ cm⁻² s⁻¹ with a photon index Γ_VHE = 5.3 ± 0.5 up to 200 GeV,⁴⁷ corresponding to an intrinsic index of 3.4 ± 0.5 after correction for the EBL (Gilmore et al. 2012, "fixed"). The day-by-day lightcurve is compatible with constant emission in that period (χ²/ndf = 7.4/6), and fractional variability F_var < 110% at the 95% confidence level (Vaughan et al. 2003).⁴⁸ Subsequent observations in May (MJD 57155–57166, 3.8 hr exposure) showed no significant excess (660 ON-source events, 3770 OFF-source events, OFF normalization of 1/6), resulting in an upper limit of Φ(>80 GeV) < 4.3 × 10⁻¹¹ cm⁻² s⁻¹ at the 99% confidence level. These results have been cross checked with an independent calibration and analysis. Monte-Carlo simulations indicate systematic uncertainties on the VHE energy scale and photon index of 20% and 0.2, respectively. The systematic uncertainty on the flux of this source is estimated to be 60%, including the energy-scale uncertainty discussed in Archambault et al. (2014).

The LAT pair-conversion telescope onboard the Fermi satellite has surveyed the whole sky in the HE band since 2008 August (Atwood et al. 2009). We analyzed the LAT data using the public science tools v10r0p5 (Pass-8) leaving free the parameters of sources from the 3FGL (Acero et al. 2015) within a region of interest of 10° radius and fixing them for sources 10°–20° away. We reconstruct the spectrum of PKS 1441+25 between 100 MeV and 100 GeV in four-week (MJD 54705–57169, top panel in Figure 1) and two-week (MJD 57001–57169, middle panel) bins assuming a power-law model with a free normalization and photon index (purple points), as well as in one-day bins (pink points) fixing the photon index to its best-fit average value in MJD 57001–57169, Γ_HE = 1.97 ± 0.02, slightly harder than in the 3FGL, 2.13 ± 0.07. The source was in a high state during MJD 57001–57169, with an integrated 100 MeV–100 GeV flux that is one to two orders of magnitude above the 3FGL value, (1.3 ± 0.1) × 10⁻⁸ cm⁻² s⁻¹. During the period contemporaneous with the VERITAS detection (MJD 57133–57140), the source shows a flux of (34 ± 4) × 10⁻⁸ cm⁻² s⁻¹ and a hard index of 1.75 ± 0.08. Although a power-law model is used for robustness in the lightcurve determination, the spectrum shows a hint of curvature, with a log parabola preferred over a power law by 3.2σ (see Figure 2). The curvature is resilient to changes in the analysis and the temporal window, and fits in smaller energy ranges confirm the hint.

X-ray observations with NuSTAR and Swift were triggered following the VHE detection. NuSTAR, a hard-X-ray instrument sensitive to 3–79 keV photons (Harrison et al. 2013), observed the source on MJD 57137 for an exposure of 38.2 ks. The data were reduced using the NuSTARDAS software v1.3.1. Swift-XRT (Gehrels et al. 2004) observed PKS 1441+25 between 0.3 and 10 keV in 2010 June (MJD 55359), in 2015 January (MJD 57028 and 57050), in 2015 April (MJD 57127–57138), and in 2015 May (MJD 57155 and 57160). Data taken in photon-counting mode were calibrated and cleaned with xrtpipeline using CALDB 20140120 v.014. ON-source and background events were selected within regions of 20-pixel (∼46 arcsec) and 40-pixel radius, respectively. The XRT and NuSTAR spectral analyses were performed with XPSEC v12.8.2, requiring at least 20 counts per bin. The NuSTAR spectrum is matched by a power law with a photon index of 2.30 ± 0.10 and an integrated 3–30 keV flux of (1.25 ± 0.09) × 10⁻¹² erg cm⁻² s⁻¹. No intranight variability is detected. Swift-XRT did not significantly detect the source in 2010 June, but the 2015 observations reveal a power-law spectrum with no detectable spectral variability and an average 0.3–10 keV photon index of 2.35 ± 0.24, using an absorbed model with a hydrogen column density of 3.19  ×  10²⁰ cm⁻² (Kalberla et al. 2005). Significant flux variations are detected in the period contemporaneous with VERITAS observations (χ²/ndf = 25.9/3, F_var = 22.6 ± 0.9%), with a flux-halving time of 13.9 ± 1.4 days based on an exponential fit to the data in MJD 57127–57155 (χ²/ndf = 8.6/6).

Simultaneously with the XRT observations, Swift-UVOT (Roming et al. 2005) took photometric snapshots of PKS 1441+25 in six optical-to-ultraviolet filters. Flux densities were extracted using uvotmaghist and circular ON-source and background regions of 5 and 15 arcsec radius, respectively. PKS 1441+25 has been observed in the V-band since 2012 January within the All-Sky Automated Survey for Supernovae (ASAS-SN; Shappee et al. 2014), using the quadruple 14-cm "Brutus" telescope in Hawaii. The fluxes from both experiments in Figure 2 are dereddenned using E(B − V) = 0.043, consistent with the column density used for the XRT analysis (Jenkins & Savage 1974). The 0.68-m Catalina Schmidt Telescope (AZ) has also performed long-term unfiltered optical observations of PKS 1441+25 since 2005 within the Catalina Real-time Transient Survey (CRTS, Drake et al. 2009). Observed magnitudes were converted into the V-band using the empirical method described in Drake et al. (2013). The SPOL spectropolarimeter (Schmidt et al. 1992) has monitored the linear optical polarization of PKS 1441+25 in 5000–7000 Å, with observations at the 1.54-m Kuiper Telescope, at the 6.5-m MMT, and at the Steward Observatory 2.3-m Bok Telescope (AZ). The source shows a high degree of polarization, with values ranging from 37.7 ± 0.1% to 36.2 ± 0.1% between MJD 57133 and MJD 57140.

The OVRO 40-m telescope (Richards et al. 2011) has monitored PKS 1441+25 at 15 GHz since late 2009. A 15 GHz VLBA image obtained by the MOJAVE program (Lister et al. 2009) on 2014 March 30 (MJD 56381) shows a compact core and a bright, linearly polarized jet feature located 1.2 mas downstream, at position angle −68°. Both features have relatively high fractional polarization (∼10%), and electric vectors aligned with the jet direction, at an angle of 102° similar to that measured by SPOL, indicating a well-ordered transverse magnetic field. The fractional polarization level of the core feature is among the highest seen in the MOJAVE program (Lister et al. 2011).

The 2008–2015 observations of PKS 1441+25 shown in Figure 1 reveal a brightening of the source in the radio, optical, and HE bands starting around MJD 56900. A simple Pearson test (see caveats in Max-Moerbeck et al. 2014a) applied to the radio and HE long-term lightcurves shows a correlation coefficient r = 0.75 ± 0.02, differing from zero by 5.4σ based on the r-distribution of shuffled lightcurve points. Similarly, the analysis of the optical and HE lightcurves yields r = 0.89 ± 0.02, differing from zero by 4.8σ. The discrete correlation functions display broad, zero-centered peaks with widths of ∼100 days, indicating no significant time lags beyond this time scale. During the period marked by gray dashed lines in Figure 1, observations on daily timescales from optical wavelengths to X-rays reveal fractional flux variations smaller than 25%, compatible with the upper limits set by Fermi-LAT and VERITAS (30% and 110% at the 95% confidence level, respectively). Such flux variations are small with respect to the four orders of magnitude spanned in νF_ν, enabling the construction of a quasi-contemporaneous spectral energy distribution in Section 3.

Zoom In Zoom Out Reset image size

The spectral energy distribution, with the X-ray-to-VHE data averaged over the active phase in 2015 April (MJD 57133–57140), is shown in Figure 2. The optical-to-X-ray spectrum is well described by a power law with photon index Γ = 2.29 ± 0.01 from 2 to 30 keV, including a 10% intrinsic scatter in the fit procedure that accounts for the small-amplitude optical-to-UV variability. This spectrum suggests a single synchrotron component peaking below 2 eV ∼ 5 × 10¹⁴ Hz, created by an electron population of index p = 2Γ − 1 ∼ 3.58 ± 0.02. As expected in FSRQs (Fossati et al. 1998), the emission of PKS 1441+25 is dominated by the gamma-ray component, well-described by a single component peaking at ${3.3}_{-1.1}^{+1.8}\;{\rm{GeV}}.$

Zoom In Zoom Out Reset image size

The detection of gamma-rays up to 200 GeV, about 400 GeV in the galaxy's frame, suggests that the emitting region is located beyond the BLR, or else pair production would suppress any VHE flux even for a flat BLR geometry (Tavecchio & Ghisellini 2012). The elevated radio state, correlated with the optical and HE brightening, also suggests synchrotron emission outside of the BLR where synchrotron self-absorption is smaller. The hypothesis of large-scale emission is strengthened by the week-long duration of the optical-to-gamma-ray flare. This behavior contrasts with other observations of bright FSRQs, displaying different flux variations at different wavelengths (e.g., Abdo et al. 2010), more in line with multi-component scenarios. The flare of PKS 1441+25 appears to be one of the few events whose detailed temporal and spectral multiwavelength features are consistent with the emission of a single component beyond the BLR.

The BLR size can be derived using the estimated black-hole mass, ${M}_{{\rm{BH}}}={10}^{7.83\pm 0.13}{M}_{\odot }$ (Shaw et al. 2012), assuming ${r}_{{\rm{BLR}}}\simeq {10}^{17}\;{\rm{cm}}\times \sqrt{{L}_{{\rm{disk}}}/{10}^{45}\;{\rm{erg}}\ {{\rm{s}}}^{-1}}$ (Kaspi et al. 2007) and an accretion disk luminosity that is a fraction η = 10% of the Eddington luminosity. Alternatively, L_disk can be estimated from the BLR luminosity as L_disk ≃ 10 × L_BLR, with L_BLR = 10^44.3 erg s⁻¹ (Xiong & Zhang 2014). Both estimates yield r_BLR ≃ 0.03 pc, setting a lower limit on the distance between the black hole and the emitting region of r ≳ 5000 Schwarzschild radii.

We show in Figure 2 that a stationary model can reproduce the data, using the numerical code of Cerruti et al. (2013) and the EBL model of Gilmore et al. (2012). Non-stationary Klein–Nishina effects on electron cooling could be important for this source (Moderski et al. 2005). Nonetheless, the straight power law observed from optical to X-rays suggests that the upscattering electrons see photon energies that are low enough to stay out of the Klein–Nishina regime. This again indicates that the emission is outside of the BLR. We parametrize the electron–positron population by a fixed broken power law with indices p₁ = 2 and p₂ = 3.58 between γ_min = 1 and γ_max = 7 × 10⁵, with a break at γ_break = 1.2 × 10⁴ (jet frame). The particle energy density is on the order of the magnetic energy density, u_e/u_B = 1.5, with a total luminosity of 6 × 10⁴⁵ erg s⁻¹, and an infrared photon energy density of 1.4 × 10⁻⁵ erg cm⁻³ at T = 10³ K (galaxy frame), characteristic of thermal radiation from the torus. We broke the degeneracies among the parameters of the model by requiring that the minimum variability timescale in the observer frame, t_var = (1 + z)/δ × R/c, be comparable to the flux-halving timescale of about two weeks observed in X rays, where we have the best statistics to probe variability. The data are then modeled with an emitting region of radius R = 4 × 10¹⁷ cm ∼ 0.1 pc and Doppler factor δ ∼ 18. The magnetic field, B ∼ 80 mG, is tangled within the emitting frame, but compressed transversely to the motion within the observer's frame, which would explain the high optical-polarization degree, PD. Following Sasada et al. (2014), the ratio of PD over the maximum theoretical polarization Π_S = (p + 1)/(p + 7/3) constrains the angle θ at which the emitting region is viewed. We find that the Doppler factor and the geometry of the system are well reproduced by a jet Lorentz factor Γ_jet ∼ 12 and θ ∼ 26, which is within the jet opening angle, θ_jet = 1/Γ_jet ∼ 48.

The location of the emission can be roughly estimated assuming that the whole cross-section of the jet contributes to the radiation (Tavecchio et al. 2010), as $r\sim R/\mathrm{tan}{\theta }_{{\rm{jet}}}\;\sim 1.5\;{\rm{pc}},$ or 200,000 Schwarzschild radii. The region is not expected to be much more compact than in this model, as no fast large-amplitude variability is seen from optical to X-ray wavelengths, despite the statistics being sufficient to detect doubling timescales as short as days. The region could still be further away from the black hole if its size was only fraction of the jet cross-section.

The model parameters are similar to those obtained by MAGIC Collaboration et al. (2008), Böttcher et al. (2008), Tavecchio et al. (2011), and Barnacka et al. (2014) for other FSRQs, but a remarkably high break energy in the electron spectrum is needed to explain the optical-to-X-ray synchrotron emission. The break in the electron distribution is consistent with radiative cooling, but is pushed to higher energies due to the magnetic field being two to ten times lower than that inferred for other VHE FSRQs. The magnetic field is also in equipartition with the particle population, minimizing the energy budget required to produce the synchrotron emission. Finally, the jet Lorentz factor is two to four times smaller than that required for other VHE FSRQs, highlighting again the reasonable energetics of this scenario.

The redshift of PKS 1441+25, z = 0.939, and its detection up to 200 GeV provide an exceptional opportunity to study the EBL.

Gamma-rays interact with EBL photons through pair production, yielding an observed spectrum that is softer than the intrinsic spectrum. Imposing a maximum intrinsic VHE hardness can then constrain the EBL intensity (Aharonian et al. 2006). In a scenario where the HE and VHE energy photons originate from the same component (but see Stern & Poutanen 2014, for emission within the BLR), the unattenuated HE observations set an upper limit on the intrinsic hardness. Based on Section 3, we neglect any additional emission component and fit, as in Biteau & Williams (2015), an absorbed power law with free EBL normalization, α, to the VERITAS spectrum. We minimize

$$\begin{eqnarray}{\chi }^{2} & = & \displaystyle \sum _{i=1..n}\displaystyle \frac{{({\phi }_{i}-{\phi }_{0}\times {e}^{-{\rm{\Gamma }}\mathrm{log}({E}_{i}/{E}_{0})-\alpha \tau ({E}_{i})})}^{2}}{{\sigma }_{{\phi }_{i}}^{2}}\\ & & +{\rm{\Theta }}({{\rm{\Gamma }}}_{{\rm{LAT}}}-{\rm{\Gamma }})\times \displaystyle \frac{{({{\rm{\Gamma }}}_{{\rm{LAT}}}-{\rm{\Gamma }})}^{2}}{{\sigma }_{{\rm{\Gamma }}}^{2}},\end{eqnarray} \tag{ 1 }$$

where α, ϕ₀, and Γ are free parameters, E₀ is fixed to 120 GeV, ${\{{\phi }_{i}\}}_{i=1..n}$ are the fluxes measured by VERITAS at energies ${\{{E}_{i}\}}_{i=1..n},$ and Θ is the Heaviside function. The Fermi-LAT log-parabolic spectrum in Figure 2 shows a photon index Γ_LAT = 2.76 ± 0.43 at 30 GeV, where the absorption by the EBL is smaller than 5% (Franceschini et al. 2008; Domínguez et al. 2011; Gilmore et al. 2012). We account for the systematic uncertainty on the VERITAS photon index, 0.20, by imposing a maximum hardness with uncertainty σ_Γ = 0.43 ⊕ 0.20 = 0.47, where ⊕ indicates a quadratic sum. We finally marginalize the equivalent likelihood, exp(−χ²/2), over the VERITAS energy scale. The logarithm of the latter is assumed to be Gaussian, with zero mean and width 0.2, corresponding to a 20% systematic uncertainty. Equation (1) allows for an intrinsic VHE spectrum that is softer, but not harder, than the HE spectrum. This yields an EBL normalization that is consistent with α = 0 and constrained to α < 1.5 at the 95% confidence level for the model of Gilmore et al. (2012). This result is almost independent of the choice of model (Franceschini et al. 2008; Domínguez et al. 2011).

Considering both the peak and full width at half maximum (FWHM) of the cross section integrated along the line of sight (as in Biteau & Williams 2015, with various evolutions tested), the VERITAS observations constrain the near-ultraviolet to near-infrared EBL. Our constraint shown in Figure 3 is compatible and competitive below 1 μm with other state-of-the-art gamma-ray measurements from Ackermann et al. (2012), H.E.S.S. Collaboration et al. (2013b), and Biteau & Williams (2015). Although α is compatible with zero, there is no significant tension with local constraints (see Biteau & Williams 2015, brown and orange arrows in Figure 3), since the differences are only 1.5σ and 1.7σ in the peak and FWHM regions, respectively.

Zoom In Zoom Out Reset image size

The low energy threshold of VERITAS enabled the detection above 80 GeV of one of the most distant VHE gamma-ray sources (z = 0.939), in a redshift range previously accessible only to space-borne gamma-ray observatories. We obtain stringent constraints on the EBL intensity below 1 μm and conclude that galaxy surveys have resolved most, if not all, of the sources of the EBL in this region. This provides an excellent baseline for studies above 1 μm where the redshifted ultraviolet emission of primordial stars could be detected (Dwek et al. 2005; Biteau & Williams 2015).

The VHE detection of the highly polarized source PKS 1441+25 is contemporaneous with a period of hard HE emission and of enhanced flux at all wavelengths. The correlation between the radio, optical, and HE lightcurves, unusual for this class of sources (Max-Moerbeck et al. 2014b), together with slow multiwavelength variability, suggest that the multi-band flare was produced by a single region located ∼10⁴–10⁵ Schwarzschild radii away from the black hole, which is consistent with the VHE-gamma-ray escape condition.

PKS 1441+25 is by far the dimmest HE emitter of all VHE-detected FSRQs listed in the 3FGL catalog. While HE activity remains a prime trigger of VHE observations, searches for new VHE-emitting quasars could also factor in radio-to-optical brightening and synchrotron-dominated X-ray emission, as reported for PKS 1441+25. These criteria will be of particular interest if applied to distant FSRQs, possibly opening a new observational window on the jets of blazars and on the transformation of the universe's light content with cosmic time.

This research is supported by grants from the U.S. Department of Energy Office of Science, the U.S. National Science Foundation and the Smithsonian Institution, and by NSERC in Canada, with additional support from NASA Swift GI grant NNX15AR38G. We acknowledge the excellent work of the technical support staff at the Fred Lawrence Whipple Observatory and at the collaborating institutions in the construction and operation of the instrument. The VERITAS Collaboration is grateful to Trevor Weekes for his seminal contributions and leadership in the field of VHE gamma-ray astrophysics, which made this study possible.

ASAS-SN thanks LCOGT, NSF, Mt. Cuba Astronomical Foundation, OSU/CCAPP and MAS/Chile for their support.

The observations at Steward Observatory are funded through NASA Fermi GI grant NNX12AO93G.

CRTS is supported by the NSF grants AST-1313422 and AST-1413600.

The OVRO 40-m monitoring program is supported in part by NASA grants NNX08AW31G and NNX11A043G, and NSF grants AST-0808050 and AST-1109911.

The National Radio Astronomy Observatory is a facility of NSF operated under cooperative agreement by Associated Universities, Inc.

This research has made use of data from the MOJAVE database that is maintained by the MOJAVE team (Lister et al. 2009).