A THREE-YEAR MULTI-WAVELENGTH STUDY OF THE VERY-HIGH-ENERGY γ-RAY BLAZAR 1ES 0229+200

VERITAS is a ground-based imaging atmospheric Cherenkov telescope (IACT) array located at the Fred Lawrence Whipple Observatory in southern Arizona. The array consists of four 12 m diameter telescopes, each with a total mirror area of 110 m². Cameras, located 12 m in front of the dishes, contain 499 circular photomultiplier tubes (PMTs), yielding a field of view (FoV) of 35. Winston cones are installed in front of the PMTs to reduce the albedo and increase the light-collecting area of the camera by filling in the gaps between them. As an IACT, VERITAS detects the brief flashes of Cherenkov light produced by the particle shower induced when a γ-ray produces an electron–positron pair in the upper atmosphere. The reconstruction of the particle shower from the imaging of the Cherenkov flash gives the energy, time of detection and arrival direction of the initial photon. Overall, VERITAS can detect photons from 100 GeV up to 30 TeV at an energy resolution of 15% and angular resolution smaller than 01 at 1 TeV (Holder et al. 2008).

VERITAS has a three-level trigger to reduce the rate of background events from the night sky and local muons. Each shower that triggers the system is imaged by the array and stored to disk. These shower events are calibrated and cleaned using quality selections based on the number of triggered PMTs in each image and the position of the image in the camera. Then, the shape and orientation of the Cherenkov images are parameterized by their principal moments (Hillas 1985). These parameters are compared to Monte Carlo simulations of γ-ray-initiated air showers. Cuts based on the physical differences between γ-ray and hadronic showers, and optimized on data taken on the Crab Nebula, are used to reject a majority of the cosmic-ray events (which are vastly more numerous than the γ-ray showers). A bright source with a flux on the order of 10% of the Crab Nebula's flux can be detected at a significance of 5 standard deviations (σ) in 30 minutes, while a weaker source (1% of the Crab Nebula's flux) can be detected in ∼25 hr. More details on the VERITAS array, specifically the detection and analysis techniques, can be found in Acciari et al. (2008).

The VERITAS Collaboration has initiated a long-term science plan which includes the observation of relatively distant blazars with hard spectra. The goal of this strategy is to build up a database of spectral energy distributions (SEDs) from a variety of blazars whose emission can carry the signature of the EBL it traverses and to study the blazar population in greater detail. As part of this program, VERITAS observed 1ES 0229+200 for a total time of 54.3 hr from 2010 to 2012. These observations were taken over three seasons (27.9, 18.0 and 8.5 hr in 2009–2010, 2010–2011 and 2011–2012, respectively, after data quality selections for weather and other issues) and resulted in a strong detection of 11.7σ. For details on the observations see Table 1. VERITAS observed this source in a wobble configuration, where the telescopes are pointed 05 away from the source so that a simultaneous background sample can be taken along with the on-source observations (Fomin et al. 1994).

These observations resulted in 1917 on-source events and 15,704 off-source events. The off-source region is larger than the on-source region by a factor of eleven so the resulting excess is 489 events, corresponding to a γ-ray rate of (0.150 ± 0.014) photons per minute. This corresponds to an average integral flux above 300 GeV of (23.3 ± 2.8_stat ± 5.8_sys) × 10⁻⁹ photons m⁻² s⁻¹ or about 1.7% of the Crab Nebula's flux (as measured by Aharonian et al. 2006). On average, this is a similar flux to that seen by the H.E.S.S. Collaboration in 2005–2006 (1.8% of the Crab Nebula's flux) in 41.8 hr of observation (Aharonian et al. 2007c). A two-dimensional Gaussian fit to the VERITAS excess (VER J0232+202) is consistent with a point source located at R.A. = $02^{\rm h} 32^{\rm m} 48^{\rm s} \pm 2^{\rm s}_{\rm stat} \pm 6^{\rm s}_{\rm sys}$, decl. = $+20^{\circ }17^{\prime }22^{\prime \prime } \pm 23^{\prime \prime }_{\rm stat} \pm 1^{\prime }30^{\prime \prime }_{\rm sys}$. This is 98 away from the VLA position of the blazar (R.A. = 02^h32^m486, decl. = +20°17'17''; Schachter et al. 1993) and within the VERITAS point-spread function.

The spectrum shown in Figure 1 can be fitted with a simple power law, and the resulting normalization (at 1 TeV) and photon index are (5.54 ± 0.56_stat ± 1.10_sys) × 10⁻⁹ m⁻² TeV⁻¹ s⁻¹ and 2.59 ± 0.12_stat ± 0.26_sys respectively, with a χ² of 5.8 with 7 degrees of freedom (the spectral analysis procedure is described in Acciari et al. 2008). The systematic errors on the normalization and index based on observations of the Crab Nebula are 20% and 10%, respectively. These results are comparable with those seen by H.E.S.S., confirming the previously measured hardness. The data for the spectral points are given in Table 2. The spectral shapes were also derived individually for the first observing period (MJD 55118–55131), when the flux was high, and for the remaining low periods (MJD 55144–55951). These are shown as shaded regions in Figure 1. No significant change is observed in the photon index. The index is 2.53 ± 0.11_stat ± 0.25_sys during the high period and 2.64 ± 0.19_stat ± 0.26_sys during the low period while the normalization rose from (4.13 ± 0.65_stat ± 0.83_sys) × 10⁻⁹ m⁻² TeV⁻¹ s⁻¹ to (10.2 ± 1.0_stat ± 2.0_sys) × 10⁻⁹ m⁻² TeV⁻¹ s⁻¹.

Zoom In Zoom Out Reset image size

The full VHE light curve is shown in Figure 4. A fit of a constant flux to the VERITAS yearly points (shown as the solid red line in the figure, (22.9 ± 2.8) × 10⁻⁹ m⁻² s⁻¹) yields a χ² of 8.32 with 2 degrees of freedom (probability of 1.6%). The blue squares in Figure 4 show the data divided into individual observing periods, for those with statistical significance above 2σ (all points are given in Table 1). The normalized excess variance is 3.9  ±  1.7 for the yearly binned data and 0.038  ±  0.038 for the data binned by observing period. This corresponds to a fractional variability of 200% for the yearly lightcurve and 19% for the monthly lightcurve. A fit of a constant flux to the monthly data (including all observing periods, as listed in Table 1) results in a χ² value of 24.7 with 13 degrees of freedom or a probability of being constant of 2.5%. The evidence for variability in these data is not conclusive by itself, but, when considered in the context of a known variable source class and significant variability in the X-ray band, we consider it to be indicative of truly variable emission.

Zoom In Zoom Out Reset image size

The X-ray flux as measured by Swift/XRT is variable at a level of 9.1 standard deviations and the PCA data shown in Figure 3 show evidence for variability throughout the three seasons (fractional variability ∼20%), including a large flare preceding the 2011–2012 VERITAS observing season. A constant flux in the RXTE data is excluded at greater than 10 standard deviations in all three bands. The hard X-ray data from the BAT shown in Figure 2 display an interesting feature. Directly preceding the first season of VERITAS observations (where the highest VHE fluxes were measured), the BAT flux from 1ES 0229+200 reached a level not previously seen in the lifetime of the BAT instrument. The BAT flux then dropped to one of the lowest levels seen. The high flux was repeated at the end of the VERITAS observing season, where observations at VHE were not possible due to the full moon. If the same particle population is involved in both the high- and low-frequency emission, the X-ray variability seen in RXTE and Swift implies that variability should be seen at VHE.

Most previous studies using distant blazars to place a lower limit on the IGMF must assume that the measured VHE spectrum (exposure time on the order of tens of hours) is a good estimator of the time-averaged spectrum of the source over several years or more. Since the flux from 1ES 0229+200, and several other VHE blazars had previously been consistent with a constant-emission model, this assumption was made by some authors attempting to limit the strength of the IGMF (Arlen & Vassiliev 2012; Dermer et al. 2011; Huan et al. 2011; Neronov & Vovk 2010; Georganopoulos et al. 2010). However, for variable sources, the multi-year time-averaged differential flux is unknown and difficult to estimate with any reliability. Because of this inherent ambiguity, any lower limit on the IGMF derived using the measured VHE spectrum from variable sources is not robust. The observations presented here show that the constant-flux hypothesis may not be valid for 1ES 0229+200, as shown in Figure 4 and Table 1. Based on archival observations of 1ES 0229+200, Dermer et al. (2011) determined that B_IGMF ≳ 5 × 10⁻¹⁸G assuming a variability time scale of ∼3 yr. The observations presented here show that 1ES 0229+200 is variable on at least a yearly timescale. Since the spectral shape is not changing and the derived limit scales as the square root of the time scale, we can assume that this reduces the limit by a factor of $\sqrt{1/3}$ to B_IGMF ≳ 3 × 10⁻¹⁸ G. The detection of variability also modifies the conclusions of Arlen & Vassiliev (2012; who assumed that the original H.E.S.S. measurement was characteristic of the average flux). They ruled out a zero IGMF hypothesis (H0) at 99% confidence based upon the 1ES 0229+200 spectrum but stated that if variability is detected in this blazar the H0 hypothesis is not ruled out at more than 95% confidence (1ES 0229+200 is the only blazar in their sample that could have rejected the H0 hypothesis).

Several authors have developed quantitative models of the plasma physics of pair cascades in the IGMF. Broderick et al. (2012) and Schlickeiser et al. (2012) place doubts on whether measurements of blazar spectra can be used to place constraints on the IGMF even if the source is non-variable since plasma instability losses dominate over IC loses. However Broderick et al. (2012) concedes that, during a flare, IC losses might dominate since instabilities have not promptly set in. Miniati & Elyiv (2013) contradict this by stating that the relaxation time of a plasma beam is much longer than the IC cooling time so that the beam can be stable and allow for secondary γ-ray emission. In conclusion, even if 1ES 0229+200 is non-variable, it is debatable that a meaningful constraint on the magnitude of the IGMF could be derived depending on the exact physics of the pair cascade plasma beam.

The broadband SED of 1ES 0229+200 is shown in Figure 6 and the synchrotron peak is detailed in Figure 5. The Fermi-LAT data, as well as the best-fit LAT bow-tie are reproduced from Vovk et al. (2012), where they reported a detection at the level of 7 standard deviations in almost three years of observations. Vovk et al. (2012) fitted the high energy spectrum from 1 to 300 GeV using a power law with a spectral index of Γ = 1.36 ± 0.25 and normalization at 20 GeV of N₀ = (1.4 ± 0.5) × 10⁻¹⁵ MeV cm⁻² s⁻¹. Note that the Swift/BAT, Fermi-LAT, and VERITAS data are long-term average spectra (70 months for the BAT, three years for the LAT, and VERITAS), while the Swift/XRT and UVOT data are short-term averages taken during the initial VERITAS observing season in 2009–2010. The UVOT spectral points plotted in Figure 6 have been corrected for host-galaxy contribution, using the correction factors from Kaufmann et al. (2011). The availability of Swift, Fermi, and VERITAS data allows both the low-energy and high-energy peaks to be constrained. The Swift/XRT spectrum is especially hard, indicating that the synchrotron peak is located above the XRT energy band (see Figure 5), but the additional information from the BAT suggests that the peak is located between the two bands at E ≃ 10 keV.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The SED is modeled using the one-zone synchrotron-self-Compton (SSC) code of Katarzyński et al. (2001), taking into account EBL attenuation based upon the calculations of Franceschini et al. (2008)³⁷ The SSC parameter space is constrained by the algorithm described in Cerruti et al. (2013), which can be seen as a numerical extension of the constraints defined by Tavecchio et al. (1998), using in addition the information from Fermi-LAT and VERITAS. The basic idea is to define a system of equations linking SSC parameters and physical observables, and to solve this system in order to determine the set of SSC parameters which correctly describes the SED. The system of equations is obtained numerically, simulating a grid of SSC models, determining for each of them the expected values of the observables, and then performing a fit to find the best parameterization of each observable as a function of the model parameters. Given the uncertainty in the physical observables, the system of equations is solved iteratively, spanning each observable in the range ±1σ.

The SSC model by Katarzyński et al. (2001) assumes a spherical emission volume of radius R moving toward the observer with Doppler factor δ, and filled with a tangled, homogeneous magnetic field B and a nonthermal population of electrons and/or positrons N_e(γ_e). The particle distribution is parameterized by a single power-law function (with index α and normalization K), between minimum and maximum Lorentz factors γ_min and γ_max, respectively. Note that the modeling presented here ignores any possible contribution to the SED from photons reprocessed by the IGMF.

Thus, the SSC model has seven free parameters: δ, B and R, for the emitting region, and α, γ_{min; max} and K, for the particle distribution, where K is defined as the electron number density at γ_max. In order to determine the set of solutions which correctly describe the SED, we used seven physical observables: the synchrotron peak frequency and flux, the X-ray spectral index, the Fermi and VERITAS fluxes at their respective decorrelation energies, and the Fermi and VERITAS spectral indices.

Following Cerruti et al. (2013), we then simulated SSC models within the following parameter space: δ ∈ [40, 100], B ∈ [0.0005, 0.01]G, R ∈ [10¹⁵, 10¹⁷] cm, γ_min ∈ [10⁴, 10⁵], γ_max ∈ [10⁶, 2 × 10⁷], and K ∈ [10⁻¹², 10⁻⁹] cm⁻³. The value of α is fully constrained by the value of the Swift/XRT spectral index: to take into account its uncertainty we computed three different sets of solutions, for α = 2.18, 2.24, and 2.30. For each simulated SSC model we determined the values of the observables, producing a grid containing, for each combination of the six parameters, the corresponding values of the six observables. It is important to note that, when computing the simulated spectral index observed in the VERITAS energy range, we excluded the upper limit at E ≃ 16 TeV. In fact, given the strong EBL absorption at E ⩾ 10 TeV, the simulated SSC models significantly under-estimate the VERITAS power-law fit at E ≃ 16 TeV, even though they are fully consistent with the flux upper limit.

The grid is then fitted, determining the system of equations linking parameters and observables, which is then solved for the specific case of 1ES 0229+200. We iteratively solved the system spanning the range (in logarithm, except for the Fermi and VERITAS spectral indices): ν_sync-peak ∈ [18.45, 18.55] Hz ([11.66, 14.67] keV), νF_{ν; sync-peak} ∈ [ − 10.99, −10.93] erg cm⁻² s⁻¹, Γ_Fermi ∈ [1.08, 1.64], Γ_VERITAS ∈ [2.30, 2.88], νF_{ν; Fermi} ∈ [ − 12.25, −11.91] erg cm⁻² s⁻¹ and νF_{ν; VERITAS} ∈ [ − 12.15, −11.96] erg cm⁻² s⁻¹. For the γ-ray observables, the uncertainty includes systematic errors, summed in quadrature with the statistical errors. The system of equations includes an inequality relating the variability time-scale (fixed to 33 days, corresponding to the doubling time-scale measured in the Swift data) to the size and the Doppler factor of the emitting region, and we reject solutions with δ > 100 which would represent a strong violation of the constraints determined from radio observations of relativistic jets (see, for example, Lister et al. 2013) or from the unification model of AGN (Henri & Saugé 2006). In order to reproduce the break observed by UVOT, we introduced in the algorithm a new observable, the break frequency ν_break, defined as the intersection of the two power-law functions fitted between 10¹⁰ and 10¹² Hz, and 10¹⁶ and 10¹⁷ Hz. We then computed a new equation linking ν_break to the six free parameters. Once the set of solutions is obtained, we then select only those solutions which are characterized by log₁₀ν_break ∈ [14.4, 14.6] Hz.

To verify the accuracy of our result we produce an SSC model for each solution and compute the χ² with respect to the observational data (XRT, LAT and VERITAS only) in order to determine the solution that minimizes the χ² and to check that all the solutions are included in a 1σ confidence interval.

The best-fit parameter values are given in Table 4, alongside two previous modeling efforts by Tavecchio et al. (2009) and Kaufmann et al. (2011). Figure 7 details the possible values for B and δ. The main difference with respect to these previous models comes from the value of the synchrotron peak frequency, which, in our SED, is located between the XRT and the BAT energy bands, at ν_sync-peak ≃ 3 × 10¹⁸ Hz (12.4 keV), an order of magnitude less than that reported by Kaufmann et al. (2011; 3.5 × 10¹⁹ Hz), but more in line with what Tavecchio et al. (2009) found (9.1 × 10¹⁸ Hz). This difference arises mainly from the higher statistics in the 70 month BAT spectrum compared with the 58 month spectrum previously used. Another difference compared to the previous modeling attempts of 1ES 0229+200 in an SSC scenario is that neither Tavecchio et al. (2009) nor Kaufmann et al. (2011) had the information from the Fermi-LAT detection.

Zoom In Zoom Out Reset image size

The first important result of the current modeling is that the minimum value of the Doppler factor required to fit the SED of 1ES 0229+200 is δ ⩾ 53. The main observational constraint on this parameter is the hard VHE spectral index, and the solutions characterized by the lowest values of δ are the ones with the softest VHE emission. This value is higher than the ones commonly assumed in SSC modeling of HSPs (see, for example, Abdo et al. 2011a, 2011b; Abramowski et al. 2012) but in agreement with the one adopted by Tavecchio et al. (2009; δ = 50), while Kaufmann et al. (2011) adopted δ = 40. It should be noted that δ = 40 is the smallest Doppler factor that could reproduce the spectra of 1ES 0229+200 in Kaufmann et al. (2011) and larger values could also have been used which would have been more in line with this work. Our solutions show some important differences with respect to the model fit performed by Kaufmann et al. (2011); the size of the emitting region (located between 5 × 10¹⁵cm and 3 × 10¹⁶cm or 1.62–9.72 mpc) is three orders of magnitude lower than that previously derived. As a consequence, the magnetic field assumed by Kaufmann et al. (2011) is several orders of magnitude lower than the one assumed here and in Tavecchio et al. (2009). An emitting-region size of the order of 10¹⁵–10¹⁶cm is similar to the ones inferred for the VHE HSPs 1ES 1218+304 (Weidinger & Spanier 2010) and 1ES 1101−232 (Aharonian et al. 2007a).

In our modeling, we found that a parameterization of the electron distribution by a single power-law function provides a good description of the SED. However, a break in the spectrum is expected in the presence of synchrotron cooling, but it is possible that the break energy is above the value of γ_max (or that the break is coincident with γ_max, i.e., that the particle distribution extends above γ_max with an index α₂ = 2.3  +  1 = 3.3). Following the study presented in Cerruti et al. (2013), we compared our values of γ_max to the expectations from synchrotron cooling, and we found that γ_max would be consistent with a synchrotron break only if the injected particles are escaping from the emitting region with a speed v comprised between c/200 and c/50. If instead the particles are escaping faster (v ≃ c), the synchrotron break is expected to be at a Lorentz factor higher than γ_max, in agreement with our modeling. Another interesting aspect is the energy budget of the emitting region. For all the solutions, we compute the values of the magnetic energy density u_B = B²/8π (in CGS units) and the particle energy density u_e = mc²∫dγ_eγ_eN_e(γ_e). The equipartition factor u_e/u_B is between 2 × 10⁴ and 10⁵, implying an emitting region significantly out of equipartition (see Figure 7).

An additional point is that the value of γ_min, constrained between 2.5 × 10⁴ and 4.5 × 10⁴, is unusually high compared to standard SSC modeling of blazars. The fact that the modeling of hard-VHE-spectrum HSPs requires such a high value of γ_min has been previously noted by Katarzyński et al. (2006), who claimed that γ_min ⩾ 10⁴ can be a characteristic of this kind of source. Katarzyński et al. proposed two alternatives to explain a high γ_min: either the injected particle population is characterized by a low-energy cut off and no cooling mechanism is efficient enough to populate the low-energy part of the spectrum, or there is an equilibrium between the cooling and the reheating of particles due to stochastic particle acceleration.

Finally, as stated previously, if the same population of electrons is responsible for both the X-ray and VHE emission then variability at X-ray energies should imply VHE variability. In the SSC scenario a correlation between X-rays and VHE photons is naturally expected if the scattering occurs in the Thomson regime. The onset of Klein–Nishina effects can be computed as a function of the Doppler factor of the emitting region (see Tavecchio et al. 1998). For the modeling presented here, assuming δ = 53, we obtain that 100 GeV photons are produced by soft X-rays in the Thomson regime, while hard X-rays (above 5 keV) are already Comptonized in the Klein–Nishina regime. On the other hand, for 10 TeV photons the scattering of X-rays is entirely (above 0.05 keV) happening in the Klein–Nishina regime. Therefore, we do expect a correlation between X-rays and γ-rays if the variability in X-rays is characterized by a simple variation in the overall normalization of the synchrotron component; on the other hand, if the peak of the synchrotron component shifts toward higher energies, this would not affect the measured VHE flux, nor the spectral index.