Search for Magnetically Broadened Cascade Emission from Blazars with VERITAS

The intergalactic magnetic field (IGMF) is a postulated weak magnetic field permeating the voids between cosmological filaments. It provides a plausible seed field for the strong magnetic fields observed in galaxies and galaxy clusters, and is thus relevant for developing a complete picture of large-scale structure formation (see Durrer & Neronov (2013) for a review). Charged cosmic rays will be deflected by the IGMF, complicating attempts to search for correlations between observed ultra-high-energy cosmic rays and potential extragalactic sources, such as active galactic nuclei (AGNs) (see, e.g., Sigl et al. 2004).

A number of mechanisms have been discussed for the generation of the IGMF (see Durrer & Neronov (2013) for a review). The field could have been generated in the early universe during the epoch of inflation, the electroweak phase transition, or during recombination (Grasso & Rubinstein 2001). A non-primordial IGMF could be generated via injection of magnetized plasma into the intergalactic medium by galactic winds (Bertone et al. 2006). Many of these scenarios result in specific predictions for the IGMF strength and correlation length (the maximum length over which the magnetic field can be treated as coherent) that distinguish them from alternative models, although degeneracies exist for some combinations of strength and correlation length. Precise constraints on the parameters of the IGMF are needed to constrain models of the field's generation.

Though the strength and correlation length of the IGMF are constrained observationally, a broad swathe of values for these quantities remains both theoretically and observationally allowed. Upper limits on the strength of the IGMF have been set with three methods: Zeeman splitting measurements of spectral lines, and Faraday rotation measurements of distant quasars and of the cosmic microwave background (CMB) (Blasi et al. 1999; Heiles & Troland 2004; Ade et al. 2015a, 2016). For correlation lengths of 1 Mpc, the upper limits are on the order of 10⁻⁹ G. Lower limits on the IGMF strength have been set based on studies of the gamma-ray emission from distant AGNs, described in more detail below. Assuming correlation lengths of at least 1 Mpc, field strengths below about 10⁻¹⁹ G have been excluded (Finke et al. 2015), in addition to a small range of field strengths around 10⁻¹⁵ G (Abramowski et al. 2014).

Observations of AGNs provide a means to probe the IGMF strength across much of the allowed range. Over 50 AGNs are detected in the very-high-energy (VHE; >100 GeV) gamma-ray range²⁵ , most of which are high-frequency-peaked BL Lacertae objects (HBLs). The high-energy photons emitted by these sources can be used as probes of the IGMF via their interactions with the extragalactic background light (EBL; see, e.g., Dwek & Krennrich 2013) en route to the observer. As multi-TeV photons travel to the observer, they interact with EBL photons and produce electron–positron pairs. The trajectory of the electrons and positrons will be bent by the IGMF, extending the path length of the cascade emission with respect to the unscattered primary emission. The mean free path for the electrons and positrons is on the order of 10 kpc for the energies and redshifts considered in the study described below (Eungwanichayapant & Aharonian 2009). The electrons/positrons eventually up-scatter low-energy CMB photons to GeV energies via inverse-Compton scattering. The GeV photons can again pair-produce on the EBL, leading to an electromagnetic cascade. Due to the deflection of the electrons and positrons before inverse-Compton scattering, the cascade emission will be angularly broadened and time delayed (Aharonian et al. 1994; Plaga 1995). The time delay varies from hours to years, depending on the energy of the cascade photons and the IGMF strength. Furthermore, the cascade emission will have a lower average energy than the primary emission.

The energy distribution and angular and temporal properties of the cascade emission provide observable signatures, but measurements are challenging, given the limited sensitive energy range of existing gamma-ray instruments. Several efforts have been made to constrain the IGMF based on the shape of spectral energy distributions of AGNs in the GeV to TeV energy range (Taylor et al. 2011; Arlen et al. 2014; Finke et al. 2015), using spectral measurements by Fermi-LAT below ∼100 GeV and imaging atmospheric Cherenkov telescope (IACT) arrays above ∼100 GeV. The interpretation of results is complicated by the different sensitivities of the instruments used in the measurements and the use of data from non-contemporaneous observations of variable sources.

The currently operating IACT arrays can be used independently to search for an IGMF-induced angular broadening of cascade emission. However, to produce cascade emission above the energy thresholds of these instruments, the initial photons must have multi-TeV energies (see Equation (27) of Neronov & Semikoz (2009) for an approximate relation between the initial and cascade photon energies). The best candidates for searches for cascade emission with IACTs are thus extreme-HBLs, whose spectral energy distributions exhibit a high-frequency peak at ∼1 TeV and hard spectral indices (Bonnoli et al. 2015). For several of these sources, no evidence of a spectral break/high-energy cutoff is observed in the intrinsic energy spectrum. It is possible that the primary emission follows an unbroken power-law distribution to several tens of TeV for these sources, as will be discussed in Section 3.

The magnitude of the angular broadening and the time delay varies with the IGMF strength and can be divided into three regimes. For 10⁻¹² G ≲ B < 10⁻⁷ G, the electron–positron pairs will be isotropized in the vicinity of the blazar, forming a pair halo that manifests itself as a broader spatial emission than expected for a point source (Eungwanichayapant & Aharonian 2009). As the mean free path before production of the electron–positron pairs does not depend on the magnetic field strength, neither will the size of the predicted pair halo. However, as a stronger field will isotropize higher energy electron–positron pairs, in the event of a pair halo detection the energy distribution of the cascade emission will provide information about the field strength.

For a field strength of 10⁻¹⁶ G < B ≲ 10⁻¹² G, the bulk of the electron–positron pairs do not isotropize, but as in the pair halo regime, the cascade produces an angularly broadened emission component in addition to the unscattered emission (Elyiv et al. 2009). The magnitude of the broadening is proportional to the field strength.

For B < 10⁻¹⁶ G, the predicted angular broadening is too small to be resolved by currently operating IACTs, and in the VHE range the cascade can only be detected via the observation of a time delayed component following a source flare (Plaga 1995; Dermer et al. 2011) or by the observation of angularly broadened emission in the GeV band. The study here focuses on the search for angular broadening rather than time delays, and is thus insensitive to field strengths of B < 10⁻¹⁶ G.

In addition to the dependence on the intrinsic source spectrum, the projected sensitivity to angularly broadened cascade emission depends on the source redshift, including the evolution of the EBL with redshift. At redshifts z ≳ 0.2, the EBL intensity is high and results in a short mean free path for both gamma-rays and electron–positron pairs, producing cascade emission that is not easily distinguished from the primary emission in a spatial analysis. In contrast, for nearby (z ≲ 0.1) sources, the cascade emission is too broad to be easily distinguishable from the isotropic cosmic-ray background (Eungwanichayapant & Aharonian 2009). Additionally, the distance between the source and the first pair production interaction typically exceeds the distance between the source and Earth for sources that are ∼100 Mpc away, resulting in a small predicted cascade fraction.

Previous searches for angularly broadened emission around blazars have been performed by MAGIC using Mrk 501 and Mrk 421 (Aleksic et al. 2010), Fermi-LAT using a large blazar sample (Ackermann et al. 2013), VERITAS using Mrk 421 (Fernández Alonso et al. 2013), and H.E.S.S. using 1ES 1101-232, 1ES 0229+200, and PKS 2155-304 (Abramowski et al. 2014). All searches resulted in non-detections, with H.E.S.S. excluding IGMF strengths of (0.3–3) × 10⁻¹⁵ G at the 99% confidence level (CL). Chen et al. (2015) claim a pair halo detection based on Fermi-LAT blazar observations, corresponding to an IGMF strength of 10⁻¹⁷–10⁻¹⁵ G. Although not confirmed, the suggested range partially overlaps the expected sensitivity range for IACT searches for angularly broadened emission.

It has been questioned whether cascade emission can be expected to reach the observer, as it is possible that the energy of the cascade could be entirely dissipated by collective behavior of the charged particles in the cascade (Broderick et al. 2012). In this scenario, energy losses due to plasma beam instabilities would dominate over the cooling of the electron–positron pairs by inverse-Compton scattering. This proposal has been argued against elsewhere in the literature, and thus the impact of plasma instabilities remains an open question (Schlickeiser et al. 2012; Miniati & Elyiv 2013). Plasma instabilities are not considered in the cascade simulations used in this study.

VERITAS is an array of four 12 m IACTs (Holder et al. 2006) located at the Fred Lawrence Whipple Observatory in southern Arizona, USA (+31° 40'3021, −110° 57' 777, 1268 m above sea level). The four telescopes are of Davies–Cotton design (Davies & Cotton 1957). The VERITAS cameras are each instrumented with 499 photomultiplier tubes (PMTs). Dead space between the PMTs is minimized with the use of light cones. VERITAS detects Cherenkov emission induced by particle showers in the atmosphere and is sensitive to gamma-rays with energies from ∼85 GeV up to greater than 30 TeV. One of the telescopes was moved in 2009 to create a more symmetric array, improving the instrument sensitivity. A major camera upgrade was completed in 2012, which decreased the lower bound on the energy range from ∼100 GeV (Kieda et al. 2013). The instrument has a field of view of 35, an energy resolution of 15%–25%, and an angular resolution (given as the 68% containment radius) of <01 at 1 TeV (Park et al. 2015).

Data were collected in wobble pointing mode, with the camera center offset by 05 from the source position, allowing for the simultaneous collection of data from the source and background regions within the same field of view (Fomin et al. 1994).

The data used in this analysis were collected between 2009 and 2012, after the array upgrade but before the camera upgrade. Obtaining the best possible angular resolution and the lowest possible energy threshold is critical for this measurement. Consequently, the data sample was restricted to observations made with zenith angles <30° and all four telescopes operating. However, characterization of the instrument performance is equally critical, thus a mature data set was used for this work. A study of observations taken after the camera upgrade will be considered in a future publication.

The sources used in this analysis were selected for optimal sensitivity to magnetically broadened emission. The amount of cascade emission falling within the sensitive energy range of VERITAS is greatest for hard-spectrum sources with intrinsic emission to ∼10 TeV. The intrinsic emission in the VHE band is attenuated by interaction with the EBL, resulting in an observed spectral index that is softer than the intrinsic value. Attenuation in the high-energy (HE; <100 GeV) range is expected to be minimal. Thus, spectral indices measured by Fermi-LAT are used as proxies for the intrinsic source indices, with values taken from the 3LAC catalog (Ackermann et al. 2015). Sources with indices harder than ∼2 were selected.

To minimize statistical uncertainties, only sources with a detection significance of >7.5σ significance (calculated using Equation (17) of Li & Ma 1983) were selected. Although sources in this optimal range of redshifts were preferentially selected, sources at higher and lower redshifts were also included in the analysis. In the case of a detection of angularly broadened emission, this would allow tests of the redshift dependence of the broadening. An angular broadening should not be detectable for the most distant and closest sources, and the detection of a spatial extension would suggest an underestimated gamma-ray point-spread function (PSF).

A test was performed for an energy cutoff/spectral break below the highest-energy spectral point measured by VERITAS. The observed VHE spectra were corrected for EBL absorption with the fiducial model of Gilmore et al. (2012). This model was selected based on its consistency with observational constraints on the EBL intensity (Biteau & Williams 2015). Each EBL-deabsorbed spectrum was fit with a power law, and the spectral index compared to the Fermi-LAT index. If the two indices disagreed by more than 1σ, this was considered an indication of a spectral break or exponential cutoff within the energy range covered by the two spectra (note that this results in a conservative prediction of the cascade emission fraction). For PG 1553+113, the source redshift is not well constrained, making it difficult to estimate the impact of the EBL absorption on the intrinsic VHE spectrum and consequently to check for the presence of a spectral break. Consequently, only limits on the flux due to extended emission are extracted for PG 1553+113, which does not require an assumption about the energy cutoff.

During strong flares, it is expected that the primary emission will be significantly brighter than the cascade emission, decreasing the overall sensitivity to the cascade component. Furthermore, rapid spectral variability can occur during flares, increasing the uncertainty on the predicted fraction of cascade emission. Thus, periods of high source activity were removed from the data sets of the highly variable sources Mrk 421 and Mrk 501 (observations with integral fluxes below 9 × 10⁻¹⁰ cm⁻² s⁻¹ and 2 × 10⁻¹⁰ cm⁻² s⁻¹ above 200 GeV, respectively, were retained). However, the high flux data were used to verify the point source simulation procedure described in Section 5. The remaining sources did not show significant flux variability within the selected data sets. The impact of spectral variability is further addressed in Section 7.1.

The final source list is shown in Table 1, with the assumed intrinsic spectral index, indication for a spectral cutoff or break, redshift, and detection significance in the VERITAS data sample used in this study. For the assumptions about the intrinsic spectrum and EBL intensity taken, the predicted cascade fraction is below the VERITAS sensitivity for all of the sources other than 1ES 1218+304. However, these assumptions are subject to large uncertainties which strongly affect the predicted cascade fraction, as will be shown in Section 7.1. Consequently, although only 1ES 1218+304 will be used to set limits on the cascade fraction and IGMF strength, the other sources are included in the analysis, in light of the possibility that the cascade fraction could be underpredicted for the nominal set of assumptions. Limits on the flux due to extended emission will be derived for all the sources.

All data were processed with the standard VERITAS calibration and shower parameterization pipelines (Daniel et al. 2007). To achieve the best possible angular resolution, events with shower images recorded in only two telescopes were discarded, while events with three or four images were retained. Gamma–hadron separation was achieved using selection on Hillas parameters (Hillas 1985; Weekes et al. 1989), with information from multiple images combined into the mean reduced scaled width and length as described by Daniel et al. (2007). For the sources Mrk 421, Mrk 501, 1ES 1218+304, and PG 1553+113, gamma-ray events were selected with box cuts on the individual gamma–hadron separation parameters. The event selection was optimized for a gamma-ray source with a soft spectral index (Γ > 3.5), thus maximizing the sensitivity to low-energy cascade emission.

The remaining sources—H 1426+428, 1ES 0229+200, and VER J0521+211—were not detected with high significance using the standard gamma–hadron separation described above. For these sources a boosted decision tree (BDT) analysis was used, resulting in substantial improvements in the source detection significances (Krause et al. 2016). The BDT analysis incorporates the same gamma–hadron separation variables used in the standard analysis into a single discriminator.

The detection significances quoted here differ from previously published VERITAS results, due to differences in the data samples and event selection. The results on VERITAS and multiwavelength observations on H 1426+428 will be presented in a separate paper (in preparation).

Both the data analysis and the simulations described below are restricted to an energy range from 160 GeV to 1 TeV. The lower limit is motivated by the energy threshold after analysis cuts (defined as the energy above which the energy bias falls below 10%). The upper limit is imposed to minimize the systematic uncertainties associated with high-energy reconstruction, while retaining sensitivity to the low-energy cascade emission.

The distribution of the parameter θ², defined as the squared angular distance between an air shower's reconstructed arrival direction and the target's estimated location, characterizes the angular profile of a source. A straightforward test for the presence of angularly broadened emission is to compare the measured θ² distribution of a source against that of a known point source. For a point source, the width of the θ² distribution is due only to the instrument's PSF.

Each point source was modeled with the standard VERITAS Monte Carlo simulation pipeline, which uses the CORSIKA program (Heck et al. 1998) to model the interaction of gamma-rays in the atmosphere and the GrISU package²⁶ to model the detector response. The VERITAS gamma-ray PSF depends on the reconstructed gamma-ray energy, the zenith and azimuthal angles of observation (as the PSF is impacted by the Earth's geomagnetic field), and the night sky background level. Consequently, the Monte Carlo simulations were weighted event-by-event to match their data counterparts' energy and azimuthal angular distributions. The range and mean value of the simulated night sky background level were matched to observations. The simulated sources were generated assuming a zenith angle of observations of 20°, which only approximates the zenith-angle distributions of the data. To correct for the simplification, a function PSF(Ze) (where Ze is the zenith angle of observations) was derived from a large sample of Crab Nebula observations. With this function, the width w_obs corresponding to the observed zenith-angle distribution was calculated. The difference (${w}_{\mathrm{obs}}\mbox{--}{w}_{{20}^{^\circ }}$) was used to correct the fitted width of the simulated sources. The magnitude of the zenith-angle corrections ranges from 00009 to 00023.

The uncertainty in the telescope pointing of 25'' (Griffiths 2015) translates into a broadening of the θ² distribution by 00005, which was determined by shifting the assumed source position in a large Crab Nebula data sample from its nominal value within the pointing uncertainty.

The uncertainty in the energy scale of VERITAS introduces a further systematic uncertainty on the θ² distribution. An average energy-scale uncertainty of 20% was assumed. The impact of the energy-scale uncertainty on the θ² distributions varies depending on the source spectrum, and thus was calculated for each source separately. The energy distributions used to weight the simulation were shifted up and down by 20%, and the change in width of the resulting θ² distributions taken as the uncertainty. The resulting uncertainties are listed in Table 2. The energy-scale uncertainty dominates the uncertainty in the simulations, exceeding the statistical errors and pointing uncertainty by a few to several factors.

Table 2.  Fitted Width of the θ² Distribution for Data and Simulations. In the Final Column, $s=({w}_{\mathrm{data}}-{w}_{\mathrm{sim}})/\sqrt{({\sigma }_{\mathrm{data}}^{2}+{\sigma }_{\mathrm{sim}}^{2})}$

Source Name	${w}_{\mathrm{data}}\pm {\sigma }_{\mathrm{stat}}^{\mathrm{data}}$	${w}_{\mathrm{sim}}\pm {\sigma }_{\mathrm{stat}}^{\mathrm{sim}}\pm {\sigma }_{\mathrm{pointing}}\pm {\sigma }_{\mathrm{energy}\mathrm{scale}}$	s
Mrk 421	00496 ± 00003	00484 ± 00002 ± 00005 ± 00017	0.7
Mrk 501	00495 ± 00004	00481 ± 00003 ± 00005 ± 00011	1.1
VER J0521+211	00477 ± 00019	00451 ± 00002 ± 00005 ± 00024	0.8
H 1426+428	00447 ± 00047	00547 ± 00003 ± 00005 ± 00027	−1.8
1ES 0229+200	00395 ± 00040	00461 ± 00003 ± 00005 ± 00016	−1.5
1ES 1218+304	00512 ± 00012	00507 ± 00003 ± 00005 ± 00012	0.3
PG 1553+113	00497 ± 00011	00521 ± 00002 ± 00005 ± 00022	−1.0

Download table as:  ASCIITypeset image

It is worth noting that the systematic uncertainties are small in comparison with the width of the θ² distributions, on the order of several percent, indicating that the angular resolution of VERITAS is the limiting factor in the sensitivity to angularly broadened emission.

Predictions for the angular profiles of the magnetically broadened cascade emission were calculated via a dedicated Monte Carlo simulation. This simulation is based on the code presented in Weisgarber (2012) but includes substantial modifications improving both the speed and accuracy of the calculations. The code tracks the three-dimensional trajectories of electrons, positrons, and gamma-rays in a ΛCDM spacetime, and it employs the full relativistic cross sections for the pair-production and inverse-Compton interactions with the redshift-dependent CMB and EBL populations. Complete kinematic modeling for both interactions enables accurate predictions for arbitrarily small IGMF strengths, while an assumption that the electron–positron pairs do not isotropize limits the code's range of validity to IGMF strengths below ∼10⁻¹² G. 13 IGMF strengths logarithmically spaced in the range B = 10⁻¹⁶–10⁻¹³ G are considered in the following.

For each particle, the simulation also calculates the amount of time accumulated between the injection of the initial gamma-ray and the particle's instantaneous position. This time is then compared to the amount of time that would have been accumulated by a non-interacting gamma-ray propagating directly from the source to the same position. The difference between these two times is tracked by the simulation at all points along every particle's trajectory.

To model the IGMF, a simple redshift scaling ${\boldsymbol{B}}{(z)={B}_{0}(1+z)}^{2}\hat{{\boldsymbol{b}}}$ is assumed, with B₀ fixed to a constant value throughout all space and the unit vector $\hat{{\boldsymbol{b}}}$ encapsulating the direction of the field. A cubic grid with sides of length 1 Mpc in comoving coordinates represents the correlation length of the IGMF. This value of the correlation length is experimentally allowed for a broad range of IGMF strengths (Taylor et al. 2011), and has been used in similar studies (e.g., Abramowski et al. 2014). Within a cube, the code selects a fixed random direction for $\hat{{\boldsymbol{b}}}$, independent of the directions in the neighboring cubes.

For each source redshift, the code samples gamma-ray energies from a distribution uniform in logarithmic energy between 0.15 and 500 TeV, injects the gamma-rays at the redshift of the source, and tracks the resulting cascades. Particles are recorded once their comoving distance from the source is equal to the comoving distance between the source and Earth. The cascades are thinned to improve the statistical independence of the results. Weights are applied to the recorded events to account for the cascade thinning and to allow predictions to be obtained for arbitrary intrinsic spectra. The simulated EBL-corrected spectra were compared to measured spectra for several test cases, and found to be in good agreement. The predictions are valid for spectra that cut off at observed energies much lower than $500/(1+z)$ TeV. The bulk Lorentz factor and the viewing angle of the blazar jet were set as 10° and 0°, respectively.

The agreement between the measured and simulated θ² distributions was first assessed based on the derived residual distributions and a χ² probability test. Figure 1 shows the θ² distributions for Mrk 501 and 1ES 1218+304 and their simulated counterparts. The distributions are plotted for θ² = 0.00–0.10 deg² in order to show detail, but the χ² probability test is applied to the residual distributions between θ² = 0.00 deg² and θ² = 0.24 deg². The upper boundary on the θ² range was selected based on the predicted width of the θ² distribution when including a cascade component. The distributions for the remaining sources are given in Figure 5 of Appendix A. Note that the systematic uncertainties and zenith-angle corrections are not accounted for in Figure 1. Even ignoring these effects, the χ² probability test indicates good agreement between data and simulations. The p-values from the χ² probability tests are shown in Table 1. Only for PG 1553+113 does the p-value indicate a mismatch between data and simulation, at the 3σ level, a tension that is resolved when the systematic uncertainties and zenith correction are considered.

Zoom In Zoom Out Reset image size

In order to include the systematic uncertainties and the zenith-angle correction in the comparison of the simulation and data, fits of the θ² distributions were performed. The simulated θ² distributions are well-described by a polynomial multiplied with a hyberbolic secant of width w (Zitzer et al. 2013):

$$\begin{eqnarray}&&P(\theta )=({c}_{0}+{c}_{1}{\theta }^{2}+{c}_{2}{\theta }^{4})\,sech (\theta /w).\end{eqnarray} \tag{ 1 }$$

In order to facilitate comparison, the data distributions were fit with the same function, but with the allowed range of the parameters c_n restricted to the 1σ uncertainty band on the fitted values from simulation, ${c}_{n}^{\mathrm{sim}}\pm {\sigma }_{{c}_{n}^{\mathrm{sim}}}$. Note that allowing a broader range for the parameters c_n (i.e., the 2σ uncertainty band) did not change the fitted values of w.

Figure 2 shows the fitted θ² distributions for data and simulation for 1ES 1218+304. The figures for the other sources are given in Figures 6–11 of Appendix B. The best-fit width parameters w_data and w_sim are compared in Table 2. Note that for w_sim, the zenith-angle correction has been applied. The statistical and systematic uncertainties on the simulated widths are added in quadrature to produce the total uncertainty σ_sim. There is no significant discrepancy between data and simulations, nor are the data distributions systematically broader or narrower than the simulated distributions. The agreement is quantified in the figure of merit $s=({w}_{\mathrm{data}}-{w}_{\mathrm{sim}})/\sqrt{({\sigma }_{\mathrm{data}}^{2}+{\sigma }_{\mathrm{sim}}^{2})}$, shown in the final column of Table 2.

Zoom In Zoom Out Reset image size

7.1. Limits on the IGMF Strength

The projected sensitivity to broadening of the source angular distribution due to a cascade emission component hinges heavily on the intrinsic spectrum of the source. Based on the cascade simulations, the predicted cascade fraction (f_c; ratio of cascade emission to total emission) must be ≳10% to produce an angular broadening that exceeds the statistical and systematic uncertainties on the widths in the data sets studied here. Evidence for an intrinsic cutoff below several TeV leads to a predicted f_c of less than 1% for all sources but the extreme-HBLs 1ES 0229+200 and 1ES 1218+304. The source 1ES 0229+200, although not showing evidence of an intrinsic cutoff, has the softest spectral index of the sources studied, at 2.025 ± 0.150 in the HE range (Ackermann et al. 2015), which also results in a low predicted value of f_c. For 1ES 1218+304, the predicted f_c is 10%–25% for the range of magnetic fields considered, as shown in Figure 4. Consequently, of all the sources studied, only 1ES 1218+304 is used to place constraints on the IGMF strength. While a stacked analysis of all sources at similar redshifts was feasible, the combined limit would be entirely dominated by the contributions of 1ES 1218+304. Hence, a stacked analysis was not attempted.

However, several uncertainties on the predicted cascade emission remain, and their impact must be examined when deriving a limit on the IGMF strength.

• 1.   
  Intrinsic cutoff; a cutoff at energies above the highest energy VERITAS spectral point cannot be excluded. Limits on the IGMF strength were derived assuming an exponential cutoff in the intrinsic spectrum at several energies: E_C = 5, 10, and 20 TeV.
• 2.  
  Spectral variability; the assumed value of the intrinsic spectral index of 1ES 1218+304, Γ = 1.660 (Ackermann et al. 2015) is measured from the full Fermi-LAT data set. However, this does not account for any spectral variability that occurred either within this data set or over the lifetime of the blazar (this is relevant as the cascade emission can, for a high IGMF strength, experience a time delay longer than the time for which VERITAS has been operating). The dependence of the IGMF limits on Γ was tested by assuming Γ = 1.460 and Γ = 1.860, while fixing the cutoff energy to 10 TeV.
• 3.  
  EBL model; the development of the cascade depends on the photon density predicted by the input EBL model. Limits on the IGMF strength were nominally derived assuming the fiducial model of Gilmore et al. (2012). To estimate the sensitivity of the IGMF constraints to the EBL model, limits were also derived with the model of Franceschini et al. (2008), while keeping Γ = 1.660 and E_C = 10 TeV fixed. These models were selected for their consistency with the EBL measurement of Biteau & Williams (2015).

For each assumed IGMF strength and set of model assumptions, a function describing the total (cascade and primary) emission was produced for 100 values of f_c between 0 and 1. The θ² function for the cascade emission was derived by convolving the simulated cascade emission's θ² distribution with the PSF measured from the simulated point source for 1ES 1218+304. The total emission for a given f_c is described by

$$\begin{eqnarray}\mathrm{total}\,\mathrm{emission}({\theta }^{2}) & = & [(1-{f}_{c})\times \mathrm{primary}\,\mathrm{emission}\\ & & +{f}_{c}\times \mathrm{cascade}\,\mathrm{emission}],\end{eqnarray} \tag{ 2 }$$

where the function describing the primary emission is again taken from the simulated point source for 1ES 1218+304. For each value of f_c, an angular distribution for the total emission was simulated with ∼1 million events (matching the number of events for the simulated point source). The distributions were then fit with Equation (1).

The left panel of Figure 3 shows the widths extracted from the fits versus f_c for an IGMF strength of B = 10⁻¹³ G, assuming the Gilmore 2012 fiducial model, Γ = 1.660 and E_C = 10 TeV. The uncertainty bands are given by the uncertainties on the simulated width shown in Table 2, added in quadrature with the statistical uncertainty ${\sigma }_{\mathrm{stat}}^{\mathrm{data}}$. The measured width w_data for 1ES 1218+304 is shown by the black vertical line. The red dashed horizontal line shows the upper limit on the cascade fraction; values of f_c above this line are excluded at the 95% CL.

Zoom In Zoom Out Reset image size

For each set of model assumptions, the 95% CL upper limit on f_c as a function of IGMF strength was compared against the predicted cascade fraction obtained from the cascade simulations. The results are shown in Figure 4. IGMF strengths for which the upper limit falls below the predicted f_c are excluded at the 95% CL. The exclusion ranges on the IGMF strength for each set of model assumptions are summarized in Table 3.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 3.  The 95% Confidence Level Exclusion Ranges on the IGMF Strength for each Set of Model Assumptions

Γ	EC (TeV)	EBL Model	IGMF Excluded (G)
1.660	10	Gilmore2012 (fid)	5.5 × 10−15–7.4 × 10−14
1.460	10	Gilmore2012 (fid)	4.5 × 10−15–1.0 × 10−13
1.860	10	Gilmore2012 (fid)	non-constraining
1.660	5	Gilmore2012 (fid)	non-constraining
1.660	20	Gilmore2012 (fid)	5.4 × 10−15–1.0 × 10−13
1.660	10	Francheschini2008	9.1 × 10−15–5.6 × 10−14

Download table as:  ASCIITypeset image

7.2. Limits on the Flux from Extended Emission

Upper limits on the integrated flux between 160 GeV and 1 TeV from angularly broadened emission are set for all sources. The bulk of the primary emission is expected to fall in the range θ² = 0.00–0.01 deg², thus excess counts due to angularly broadened emission were calculated from the difference $\int {\theta }_{\mathrm{data}}^{2}-\int {\theta }_{\mathrm{sim}}^{2}$ within the integration range θ² = 0.01–0.24 deg². The integration range was chosen to match the ranges used in similar calculations performed by Abramowski et al. (2014) and Aleksic et al. (2010). Upper limits on the number of gamma-ray events due to angularly broadened emission are calculated using the frequentist method of Rolke (Rolke et al. 2005), and translated into an upper limit on the rate by dividing by the deadtime-corrected exposure time.

Zoom In Zoom Out Reset image size

Translating the upper limit on the rate into an upper limit on the integrated flux requires an assumption about the spectral index of the angularly broadened emission. A spectral index (Γ+2)/2 was assumed, accounting for a slight softening of the cascade emission compared to the primary emission as inverse-Compton scattering proceeds in the Thompson limit. The resulting 95% CL upper limits on the integrated flux due to angularly broadened emission for an energy range between 160 GeV and 1 TeV are shown in Table 4.

Zoom In Zoom Out Reset image size

Table 4.  Limits on the Integrated Flux from Angularly Broadened Emission in an Energy Range between 160 GeV and 1 TeV, Assuming a Spectral Index of (Γ+2)/2 for the Cascade Emission

Source Name	95% CL (10−12 cm−2 s−1)
Mrk 421	1.4
Mrk 501	4.2
VER J0521+211	0.6
H 1426+428	2.1
1ES 0229+200	1.2
1ES 1218+304	0.9
PG 1553+113	3.2

Download table as:  ASCIITypeset image

A search for source extension due to cascade emission broadened by the IGMF was performed with VERITAS observations of seven blazars. No indication of angularly broadened emission was observed. Limits were set on the fraction of the total emission due to cascade emission (f_c) for the blazar with the largest predicted cascade fraction, 1ES 1218+304. IGMF strengths between 10⁻¹⁶ and 10⁻¹³ G, and an IGMF coherence length of 1 Mpc were assumed. Exclusion regions on the IGMF strength were determined under different sets of assumptions about the source intrinsic spectrum and the EBL intensity. For a nominal set of assumptions (spectral index Γ = 1.660 and cutoff energy E_C = 10 TeV for 1ES 1218+304, EBL model of Gilmore et al. 2012), an IGMF strength of 5.5 × 10⁻¹⁵ G–7.4 × 10⁻¹⁴ G can be excluded at the 95% CL. This shows a similar sensitivity to measurements from other instruments, as well as complementarity to previous results. Namely, H.E.S.S. ruled out an IGMF strength in the range (0.3–3) × 10⁻¹⁵ G at the 99% CL, using observations of PKS 2155-304 and slightly different model assumptions but otherwise similar methodology (Abramowski et al. 2014). Taken together, the H.E.S.S. and VERITAS constraints rule out an IGMF strength falling in much of the range between 10⁻¹⁶ and 10⁻¹³ G. The VERITAS exclusion region, however, does not rule out the IGMF range suggested by the claimed detection of Chen et al. (2015).

Zoom In Zoom Out Reset image size

Varying the assumptions on the intrinsic spectrum of the source, namely the spectral index and the high-energy cutoff, substantially alters the extracted limits on f_c. Softening the assumed intrinsic spectral index by 0.2 or decreasing the energy of the exponential cutoff from 10 to 5 TeV resulted in limits on f_c falling above the predicted cascade fraction. In these cases, the IGMF strength is not constrained. The assumed shape of the cutoff impacts the constraints as well: assuming a super exponential cutoff power law exp($-{(E/{E}_{C})}^{\gamma }$) will produce stronger constraints for 0 < γ < 1 (softer cutoff) and weaker constraints for γ > 1 (sharper cutoff) than for the assumed exponential cutoff power law (γ = 1). Finally, the EBL model assumed when simulating the cascade process affects the limits. It was observed that using the model of Gilmore et al. (2012) in the cascade simulations produced a broader IGMF exclusion region than using the model of Franceschini et al. (2008).

Zoom In Zoom Out Reset image size

As the cascade emission is time delayed by years for the IGMF strengths considered here, the flux variability of the source over its lifetime will impact the limits on the IGMF strength, as the predicted f_c is derived assuming the currently observed flux of primary emission. If the source exhibited a lower (higher) flux at the time that the cascade emission reaching the observer today was produced, f_c will be lower (higher) than expected based on the current flux of 1ES 1218+308. Based on the ratio of the observed upper limits on f_c and the predicted values, the differential flux of 1ES 1218+308 at 1 TeV would have to be ∼70% of its current value on average over its lifetime to invalidate the entire nominal IGMF exclusion range.

Zoom In Zoom Out Reset image size

It should be noted that while the assumed intrinsic source spectrum and EBL model affect the IGMF limits, the results are not expected to be sensitive to the choice of Doppler factor or viewing angle of the jet, nor to the choice of correlation length. It was demonstrated in Arlen et al. (2014) that the cascade spectrum above 100 GeV is unaffected by variation of the bulk Lorentz factor between 5 and 100, or variation of the jet viewing angle between 0° and 10°. The IGMF limits are insensitive to the correlation length provided that the correlation length exceeds the inverse-Compton cooling length (Neronov & Semikoz 2009). The cooling length and primary gamma-ray energy scale inversely. The primary gamma-rays must have energies above 1 TeV to produce cascade emission with energies above the VERITAS energy threshold. Thus, the majority of the primary gamma-rays will have cooling lengths of less than a few hundred kpc, much less than the correlation length of 1 Mpc used in the cascade simulation code.

For the cutoff energy of 10 TeV assumed for the intrinsic spectrum of 1ES 1218+308, the first pair production interaction occurs >10 Mpc from the source. Consequently, this study probes the magnetic field strength in areas distant from the source, sampling cosmic voids, rather than matter-rich regions.

Limits were set on the integrated flux due to angularly broadened emission for all sources, resulting in 95% CL upper limits of (0.6–4.2) × 10⁻¹² cm⁻² s⁻¹ for an energy range between 160 GeV and 1 TeV. A spectral index of (Γ+2)/2 was assumed for the angularly broadened emission.

The IGMF constraints presented here are limited by a number of factors, however the dominant limitation on the sensitivity to angularly broadened emission is the instrument PSF. The Cherenkov Telescope Array (CTA) is projected to begin taking data in several years with a sensitivity greater than currently operating instruments, and in particular, with a substantially better PSF compared to VERITAS (Acharya et al. 2013; Meyer et al. 2016), and consequently an improved ability to probe weaker IGMF strengths and smaller cascade fractions.

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. E.P. acknowledges the support of a Marie Curie Intra-European Fellowship within the 7th European Community Framework Programme, and thanks Andrew Taylor for useful discussion. 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.

Figure 5 shows the comparison of the angular profiles for the sources Mrk 421, VER J0521+211, H 1426+428, 1ES 0229+200, and PG 1553+113 and their simulated counterparts.

Fits to the θ² distributions of Mrk 421, Mrk 501, VER J0521+211, H 1426+428, 1ES 0229+200 and PG 1553+113 are shown in Figures 6–11 for the sources and for their simulated counterparts.