Discovery of Very High Energy Gamma Rays from 1ES 1440+122

Active galactic nuclei (AGN) are observed to emit electromagnetic radiation from radio waves up to very high energy (VHE; $E>100$ GeV) $\gamma$-rays. These objects, which make up only a small fraction of the total number of observed galaxies, are very luminous, extremely compact, and can exhibit large luminosity variability. Although AGN differ widely in their observed characteristics, a unified picture has emerged in which AGN are powered by accretion onto a super-massive black hole (10⁷–10${}^{9}~$M_⊙). Near the black hole is a hot accretion disc surrounded by a thick torus of gas and dust. In some AGN (the radio-loud population, $\sim$15%), a highly-relativistic outflow of energetic particles form a highly collimated jet generating non-thermal emission. Blazars are thought to be the case where the jet is aligned with our line of sight (Urry & Padovani, 1995).

Blazar spectral energy distributions (SEDs) are dominated by non-thermal radiation. This emission has a flat radio spectrum, radio and optical polarisation, and is often highly variable. 1ES 1440+122 belongs to the BL Lacertae (BL Lac) subclass of blazars. BL Lacs do not have broad emission lines present, unlike flat-spectrum radio quasars. Blazar SEDs are characterised by two broad peaks, with a significant fraction of the power often being emitted in the $\gamma$-ray band. The low-energy peak in the SED is well understood as synchrotron emission from relativistic electrons. However, there are competing models to explain the high-energy peak emission as dominated by either leptonic or hadronic interactions (Blandford & Levinson, 1995; Bloom & Marscher, 1996; Mannheim, 1998; Pohl & Schlickeiser, 2000; Aharonian et al., 2000). BL Lacs have been further classified depending on the position of their lower-energy peak. Padovani & Giommi (1995) originally proposed two classes. Class definitions were extended to include an intermediate case by Nieppola et al. (2006) and Abdo et al. (2010), though there is not agreement on where to place the boundaries between classes. Based on a parabolic fit in log-log space to archival data, Nieppola et al. (2006) determined the location of synchrotron peak in 1ES 1440+122 to be at $\nu_{{\rm peak}}=10^{16.4}$ Hz, lying near the border between their intermediate-frequency-peaked and high-frequency-peaked BL Lac (IBL and HBL, respectively) class definitions. According to the classification scheme in Abdo et al. (2010), this synchrotron-peak frequency sets the classification of the source as a high synchrotron peaked (HSP) BL Lac.

1ES 1440+122 was initially classified as an AGN in X-rays in the Einstein Slew Survey (Elvis et al., 1992). It is surrounded by $\sim$20 galaxies within approximately 200 kpc (Heidt et al., 1999), which suggests that it may belong to a small cluster of galaxies. Indeed it is likely that the blazar is interacting with an elliptical galaxy with a projected separation of $\sim$4 kpc (Sbarufatti et al., 2006). The host galaxy has been resolved in several imaging studies, and high-resolution HST imaging (Scarpa et al., 1999) reveals a very close companion ($\sim 0.3^{\prime\prime}$) now known to be a foreground star (Giovannini et al., 2004). The optical spectrum of 1ES 1440+122 is well measured, and a redshift of $z$=0.163 is obtained from the identification of three spectral lines (Sbarufatti et al., 2006).

1ES 1440+122 was identified as a likely VHE emitter on the basis of its SED (Costamante & Ghisellini, 2002). It was observed by the H.E.S.S. array of imaging atmospheric Cherenkov telescopes (IACTs) for 11.2 hr between 2004 and 2009 resulting in a 99.9% confidence level integral flux upper limit of 1.66 $\times$ 10^-12 cm^-2 s^-1 (1.0% of the Crab Nebula’s flux) above 290 GeV (Abramowski et al., 2014). 1ES 1440+122 belongs to the 3FGL catalogue (3FGL J1442.8+1200) with spectral index, $\Gamma$ = 1.80 $\pm$ 0.12 (Acero et al., 2015). The Fermi Large Area Telescope (LAT) spectrum shows no sign for a cutoff, and its extrapolation to the VHE band predicts a $\sim$2% Crab Nebula flux between 200 GeV and 1 TeV, including extragalactic background light (EBL) absorption effects (Franceschini et al., 2008). The 3FGL lists the source as being non-variable (Acero et al., 2015).

VHE $\gamma$-ray emission from 1ES 1440+122 was discovered by VERITAS (Ong et al., 2010) during the 2008–2010 observing seasons. After the confirmation of the initial excess from this source, Swift X-ray Telescope (XRT) target of opportunity observations were triggered in order to provide a detailed X-ray spectrum, which is crucial for constraining models of emission. In combination with Fermi-LAT data, all data were used to construct an SED spanning 10 decades in energy.

VERITAS consists of four 12-m diameter IACTs located in southern Arizona (Holder et al., 2006). The array detects $\gamma$-ray emission from astrophysical objects in the energy range from $\sim$85 GeV to $\sim$30 TeV. VERITAS has an energy resolution of $\sim$15% and angular resolution (68% containment) of $\sim$0.1^∘ per event. The current sensitivity of the array allows for a 1% Crab Nebula flux source detection in 25 hours (5$\sigma$ detection), while a 10% Crab Nebula flux is detected in 0.5 hours. Note that for sources with a softer spectrum than the Crab, the observing time required for detection will be longer. The field of view of the VERITAS telescope has a diameter of 3.5^∘. More information on VERITAS and the IACT technique can be found in Holder et al. (2008).

VERITAS observed 1ES 1440+122 for about 78 hours from 2008 May to 2010 June, which includes two observing seasons, 2008–2009 and 2009–2010. The earlier data were taken with the original VERITAS telescope configuration, while the later data were taken after one of the telescopes was relocated in order to increase the array sensitivity (Perkins et al., 2009). These observations cover the zenith angle range from 19^∘ to 38^∘. All the observations were performed in a mode where the source is offset by 0.5^∘ from the centre of the field of view. This offset allows for simultaneous background estimation with good precision while maintaining a high signal efficiency. This is known as “wobble” mode (Fomin et al., 1994). Observations affected by poor weather or hardware problems were removed, and the remaining 53 hours of data were processed with two independent analysis packages (Daniel et al., 2007) yielding consistent results.

Images of the showers were first calibrated in gain and timing at the pixel level using nightly calibration data using an artificial light source. Following calibration, images from each telescope were parameterized using fits to two-dimensional Gaussian distributions. This technique is similar to the frequently-used moment analysis (Hillas, 1985). Tests performed using $\gamma$-ray simulations have shown that the use of the Gaussian fit leads to several improvements with respect to the moment analysis. Truncated images, where some part of the shower is not contained in the field of view of the camera, are reconstructed with better angular and energy resolution. This also leads to an increase in the rate of events passing selection at high energy. No image cleaning (thresholding) was used, but a constant offset represents the night-sky background. As a result, background rejection was improved at low energy. More details of this image-fitting technique are given in Christiansen et al. (2012).

From the results of the image fits, parameters were calculated and used for event reconstruction and selection. The event selection criteria were optimised beforehand using observations of the Crab Nebula with the excess counts scaled to match a 1% Crab Nebula flux source. The source region was defined by a 0.1^∘ radius circle centred on the source coordinates, and all the $\gamma$-ray-like events within this region were considered the ON counts. The reflected-region model (Berge et al., 2007) was used for background subtraction, where the background was estimated from eleven identically-sized regions reflected from the source region around the camera centre. The events found in these regions were considered the OFF counts. The Li & Ma Formula 17 (Li & Ma, 1983) was used to calculate the significance at the source location.

An excess of 166 events was observed in VERITAS data from the direction of 1ES 1440+122 (954 ON events, 8673 OFF events with an off-source normalisation ratio of $\alpha=0.0909$). The excess corresponds to a statistical significance of 5.5$\sigma$. The significance map in the vicinity of 1ES 1440+122 is shown in Figure 1 while the distribution of events with respect to the source location is shown in Figure 2. The VERITAS point spread function is 6^′ for 68% containment radius at these zenith angles making the distribution of events consistent with a point-like source. Fitting a symmetric 2-dimensional Gaussian to the uncorrelated excess counts map results in a best-fit centroid at R.A. = 14^h 43^m 15^s and Dec. = +12^∘ 00^′ 11^′′. The new TeV source is catalogued as VER J1443+120. The statistical uncertainty in the position of 1^′ and systematic uncertainty of 25^′′ make the $\gamma$-ray emission consistent with the position of 1ES 1440+122, which is reported in the SDSS¹¹1http://www.sdss.org/dr5/ as R.A. = 14^h 42^m 48.3^s and Dec. = +12^∘ 00^′ 40^′′. We fitted the integrated flux light curve $>$200 GeV in monthly bins with a constant-flux hypothesis, resulting in no statistically significant evidence of variability ($\chi^{2}$/dof =15.2/9 or a 9% chance of being generated given a constant flux hypothesis).

The spectrum of 1ES 1440+122 is described well by a power law of the form $dN/dE=I_{0}\,(E/0.5~\mathrm{TeV})^{-\Gamma}$, with $I_{0}=(1.47\pm 0.62_{\mathrm{stat}})\times 10^{-12}$ cm^-2 s^-1 TeV^-1 and $\Gamma=3.1\pm 0.4_{\mathrm{stat}}$, resulting in $\chi^{2}$/dof = 2.4/2.

The integral photon flux above 200 GeV is $\Phi_{E>200~\mathrm{GeV}}=(2.8\pm 0.7)\times 10^{-12}$ cm^-2 s^-1. Using the parameterization by Aharonian et al. (2006), this is equivalent to 1.2% of the flux from the Crab Nebula above 200 GeV. We estimate the systematic errors on the flux normalisation constant and the photon index to be $\Delta I_{0}$/I₀ = 30% and $\Delta\Gamma=0.2$. The VERITAS spectral points are shown in Figure 3. For a strict comparison with the H.E.S.S integral flux upper limit of 1.66 $\times$ 10^-12 cm^-2 s^-1 above 290 GeV (Abramowski et al., 2014), we recalculate our integral flux above the same energy threshold to be $(1.5\pm 0.7)$ $\times$ 10^-12 cm^-2 s^-1, which shows they are consistent with one another.

The Fermi-LAT is a pair-conversion $\gamma$-ray detector sensitive to photons in the energy range from about 30 MeV to 300 GeV (Atwood et al., 2009). The Fermi-LAT data for this analysis were taken between 2008 August 04 and 2010 August 02. We used the likelihood tools distributed with the standard Science Tools v9r32p5 package available from the Fermi Science Support Center²²2http://fermi.gsfc.nasa.gov/ssc. Events were required to have zenith angle $<$105^∘ in order to limit contamination by the Earth albedo effect. Only events in the energy range of 300 MeV to 100 GeV and within a circular region of 12^∘ radius centred on the source were selected. The background was modelled with a galactic diffuse emission model³³3gll_iem_v02.fit and an isotropic component. Catalogued sources within 17^∘ of the target location were included in the model. The fluxes were determined using the instrument response functions P7REP_SOURCE_V15 (a check was performed using updated response functions released after modeling was complete but these did not yield significant changes). The systematic uncertainty on the flux is approximately 5% at 560 MeV and under 10% at 10 GeV and above (Ackermann et al., 2012a).

A point source was detected at the position of 1ES 1440+122 with a significance of more than 9 standard deviations (TS$=$89.6). The highest-energy photon detected by the Fermi-LAT with a high probability of association with the source (99.3%) has an energy of 62 GeV. The time-averaged Fermi-LAT spectrum computed assuming a constant power law includes five points and is shown in Figure 4. The spectrum is fit by a power law that can be described as $dN/dE=N_{0}(1-\Gamma)E^{-\Gamma}/(E_{\mathrm{max}}^{1-\Gamma}-E_{\mathrm{min}}^{1-\Gamma})$, where $\Gamma=1.52\pm 0.16_{\mathrm{stat}}$, $N_{0}=(5.39\pm 1.18_{\mathrm{stat}})\times 10^{-10}$ photons $\textrm{cm}^{-2}$ $\textrm{s}^{-1}$, $E_{\mathrm{min}}=$1 GeV and $E_{\mathrm{max}}=$100 GeV. The values reported in the 3FGL are: $\Gamma=1.80\pm 0.12$ and $N_{0}=(5.59\pm 0.86)\times 10^{-10}$ photons $\textrm{cm}^{-2}$ $\textrm{s}^{-1}$ with the same min and max energy (Acero et al., 2015). There is some tension between the spectral index we found and the 3FGL value. The period of our observations matches the 2FGL catalogue more closely, and our spectral index is in agreement with the 2FGL value of $\Gamma=1.41\pm 0.18$ (Nolan et al., 2012).

Swift-XRT (Gehrels et al., 2004) observations of 1ES 1440+122 were performed on 2008 June 12 and 2010 March 9. The Swift-XRT data were analysed with HEAsoft⁴⁴4http://heasarc.nasa.gov/lheasoft 6.13. Both observations were completed in photon-counting mode, showing count rates of 0.40 $\pm$ 0.01 and 0.38 $\pm$0.02 counts s^-1. With these low count rates, photon pile-up is negligible and systematic uncertainties on the flux are negligible compared to the statistical errors.

XSPEC version 12.8.0 was used for the XRT spectral analysis. The data were combined and grouped into bins with a minimum of $20$ counts per bin, enabling the use of $\chi^{2}$ spectral-model fitting. The time-averaged 0.3-10 keV data were fitted with an absorbed power-law model (tbabs(po) in XSPEC) where the hydrogen column density $N_{H}$ was fixed at 1.58 $\times$ 10²⁰ cm^-2, taken from the LAB survey of Galactic HI (Kalberla et al., 2005). The data were reasonably well fit ($\chi^{2}/$dof $=61.9/45$) by an absorbed power law with normalisation at 1 keV of (3.2 $\pm$ 0.1) $\times$ 10^-3 s^-1 keV^-1 and photon index of 1.95 $\pm$ 0.04. In order to represent the intrinsic X-ray emission, the absorption-corrected spectrum was used for SED modelling. The fluxes during both periods are consistent with one another, so an average was used.

The Ultraviolet/Optical Telescope (UVOT; Roming et al. (2005)) onboard Swift observed 1ES 1440+122 simultaneously with the Swift-XRT time periods. Source photons in each of six filters (V, B, U, UVW2, UVM2, and UVW1) were extracted from a circular aperture of radius 5.0^′′ centred on the source. The background was estimated in a 30^′′ radius circular region located away from the blazar. The fluxes were computed using the uvotsource⁵⁵5HEASOFT v6.13, Swift_Rel4.0(Bld29)_14Dec2012 with calibrations from Breeveld et al. (2011). tool. Corrections for interstellar absorption were made using the extinction curve of Fitzpatrick (1999) and assuming an $E(B-V)$ value of 0.0239$\pm$0.0006, determined from the IPAC Extragalactic Database (Schlafly & Finkbeiner, 2011). The average fluxes from the two observations were used for the SED modelling.

The quasi-simultaneous SED of 1ES 1440+122 with data from VERITAS, Fermi-LAT, and Swift-XRT is shown in Figure 4 and parameters of the three models, discussed below, are listed in Table 1. To these data, we added archival low-frequency radio and optical data from the NED⁶⁶6http://nedwww.ipac.caltech.edu as well as a Swift-BAT point⁷⁷7http://tools.asdc.asi.it. The vertical bar between the two B-band points illustrates the amount of historical variability at that frequency. The IR-through-optical emission of 1ES 1440+122 is clearly dominated by thermal emission from the host galaxy, modelled as a blackbody spectrum here. The UV data are under-represented in all models that we investigated. This could be additional contamination from the host galaxy (whose emission has been parametrized by a simple blackbody though it is known to be more complex) or extended jets unrelated to the VHE emission. The latter scenario has been considered for PKS 2155-304, which does not show correlated variability between the two bands (Abramowski et al., 2012). Optical polarisation data may be useful for disentangling contamination from other parts of the jet (de Almeida et al., 2014).

We produce models of the SED of 1ES 1440+122 with both leptonic and hadronic jet models. The VHE emission is corrected for EBL absorption according to Franceschini et al. (2008), which is in agreement with the most recent constraints from gamma-ray observations (Abramowski et al., 2013; Ackermann et al., 2012b; Biteau & Williams, 2015). In leptonic models for blazar emission, a population of relativistic electrons is responsible for both the lower-frequency component of the SED (via synchrotron emission) as well as the higher-frequency emission (via Compton scattering). Potential soft photon fields that can serve as targets for Compton scattering are either the synchrotron photons (SSC = synchrotron self-Compton), or radiation fields produced externally to the jet (EC = External Compton). For our models, we use a steady-state scenario in the fast-cooling regime based on the time-dependent blazar jet radiation transfer code of Böttcher & Chiang (2002), as described in detail in Böttcher et al. (2013). In this model, the emission originates from a spherical region of radius $R$, moving along the jet with a Lorentz factor $\Gamma$, corresponding to a jet speed $\beta_{\Gamma}c$. The jet is oriented at an angle $\theta_{\rm obs}$ with respect to the line of sight, resulting in Doppler boosting characterised by the Doppler factor $D=(\Gamma[1-\beta_{\Gamma}\cos\theta_{\rm obs}])^{-1}$.

Non-thermal electrons are injected and accelerated into a power-law distribution at a rate $Q(\gamma)=Q_{0}\gamma^{-q}$ between a low- and high-energy cutoff, $\gamma_{1,2}$. A value of $q=3.0$ was chosen for all models, though it is not well constrained by the observations and can take a wide range of values depending on the obliquity and shock velocity (Summerlin & Baring, 2012). An equilibrium between this particle injection, radiative cooling and particle escape is established self-consistently with the radiation mechanisms considered. Particle escape was parameterized through an escape parameter $\eta$ such that the escape time scale $t_{\rm esc}=\eta\,R/c$. The resulting particle distribution will correspond to a power $L_{e}$ in electrons streaming along the jet (see Acciari et al., 2009). The synchrotron emission is evaluated assuming the presence of a tangled magnetic field $B$, corresponding to a power in Poynting flux, $L_{B}$. For each model calculation, our code evaluates the equipartition parameter $\epsilon_{Be}\equiv L_{B}/L_{e}$. Because of a lack of observational constraints and in order to reduce the number of free parameters, we choose the observing angle as the critical angle, for which $\Gamma=D$, i.e., $\cos\theta_{\rm obs}=\beta_{\Gamma}$.

In a pure SSC model, only synchrotron photons play the role as targets for Compton scattering. The SSC model satisfactorily produces the non-thermal SED with plausible parameters. The size of the emission region used in the model implies a minimum variability time scale allowed by the model given by $t_{\rm var}^{\rm min}=R(1+z)/(cD)=1.7$ hours. The required magnetic field energy density is a factor of approximately 35 below equipartition with the non-thermal electron distribution.

For a model including an external radiation field as target for Compton scattering, we have improved the model presented in Acciari et al. (2009) by allowing for isotropic (in the rest frame of the AGN) radiation fields with arbitrary spectra. Guided by recent results of EC modelling of SEDs of other VERITAS-detected IBLs, such as W Comae (Acciari et al., 2008, 2009) and 3C66A (Abdo et al., 2011), we consider a thermal infrared radiation field, possibly originating in a dusty torus around the central engine, as an appropriate choice for an external radiation field. The energy density $u_{\rm ext}$ of this photon field and temperature $T_{\rm ext}$ of the dusty torus used in the model are poorly constrained but consistent with expectations. This model also satisfactorily represents the SED and allows for the choice of parameters very close to equipartition between the magnetic-field and non-thermal electron energy densities. The minimum variability time scale is 4.2 hours.

In addition to the purely leptonic models described above, we also consider a lepto-hadronic model, in which ultrarelativistic protons contribute significantly to the high-energy emission through proton-synchrotron radiation and $p\gamma$ pion production. The spectra of $\pi^{0}$ decay photons as well as the final decay products of charged pions are evaluated using the templates of Kelner & Aharonian (2008), accounting for secondary cascades as described in Böttcher (2010). In our model, in addition to the SSC model outlined above, we assume a power-law distribution of relativistic protons, $n(\gamma)\propto\gamma^{-p}$ between a low- and high-energy cutoff, $\gamma_{\rm p}^{\rm min,max}$, normalised to a total kinetic luminosity $L_{p}$ of the proton population propagating along the jet. We then evaluate the energy partition fractions $\epsilon_{Bp}\equiv L_{B}/L_{p}$ and $\epsilon_{ep}\equiv L_{e}/L_{p}$. As for the other models, the result is shown in Figure 4, and the model parameters are listed in Table 1. This model also adequately produces the non-thermal SED. It requires a strongly-magnetically dominated jet with $\epsilon_{Be}=1.4\times 10^{6}$ and $\epsilon_{Bp}=24$. Under these conditions, $\gamma$-ray emission is dominated by proton-synchrotron radiation. The minimum variability time scale from the size of the emission region is just 4.2 hours. However, the radiative cooling time of ultrarelativistic protons is on the order of several days, excluding variability on shorter time scales (not yet seen) under this model.

1ES 1440+122 was detected by VERITAS at a significance level of 5.5$\sigma$ during the 2008–2010 observing seasons. In this paper, we described VHE observations of 1ES 1440+122 along with the quasi-simultaneous observations with Swift in optical, UV, and X-rays and Fermi-LAT in high-energy $\gamma$-rays. The observed non-thermal SED of 1ES 1440+122 is consistent with purely leptonic (SSC and EC) models as well as with a hadronic origin. A leptonic model with an external infrared radiation field as target for Compton scattering allowed for parameters close to equipartition between the relativistic electron population and the magnetic field. The other models did not allow for partition fractions near unity. The model parameters are comparable to those obtained from other studies of VHE BL Lacs.

Our model SEDs all result in synchrotron peak frequencies contained in the Swift-XRT band, very close to $\nu_{\rm{synch}}=3\times 10^{17}$ Hz, classifying the source as an HBL according to the scheme of Nieppola et al. (2006) or an HSP according to Abdo et al. (2010). The Compton-peak frequencies are all close to $\nu_{\rm{Compton}}=3\times 10^{25}$ Hz. The SED is dominated by the lower synchrotron peak (i.e. low Compton dominance), which is generally observed in other VHE HBLs (Fossati et al., 1998; Ghisellini et al., 1998).

Our results show that the SED alone does not allow us to confidently distinguish between different models for the high-energy emission from 1ES 1440+122. Future observations, such as probing for intraday variability of VHE $\gamma$-rays, may aid in distinguishing leptonic and hadronic models. The radiative cooling time scales of ultrarelativistic protons are of the order of several days, whereas those of ultrarelativistic electrons are typically of the order of hours or less. Therefore rapid (intraday) VHE $\gamma$-ray variability would be an indication of leptonic processes dominating the $\gamma$-ray output.

As a VHE source with a relatively hard spectrum for its redshift, 1ES 1440+122 could be useful in studies of the EBL. Several studies have used similar blazar spectra to examine lower limits on the intergalactic magnetic fields (Neronov & Vovk, 2010; Huan et al., 2011; Dermer et al., 2011; Arlen & Vassiliev, 2012). These studies look for emission in the Fermi-LAT band that might have cascaded down from the VERITAS band. However, there is an ongoing debate in the literature regarding the validity of these limits given the possibility of plasma instability energy losses dominating over the inverse Compton losses, resulting in less lower-energy cascade emission (Broderick et al., 2012; Schlickeiser et al., 2012; Miniati & Elyiv, 2013).

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, by NSERC in Canada, by Science Foundation Ireland (SFI 10/RFP/AST2748) and by STFC in the U.K. 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. M. Böttcher acknowledges support by the South African Department of Science and Technology through the National Research Foundation under NRF SARChI Chair grant no. 64789. 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.