INTERGALACTIC MAGNETIC FIELDS AND GAMMA-RAY OBSERVATIONS OF EXTREME TeV BLAZARS

Timothy C. Arlen; Vladimir V. Vassilev; Thomas Weisgarber; Scott P. Wakely; S. Yusef Shafi

doi:10.1088/0004-637X/796/1/18

1. INTRODUCTION

Numerous observations have established the presence of magnetic fields in our own galaxy and in other galaxies and galaxy clusters on the order of a microgauss (see, e.g., Grasso & Rubinstein 2001; Widrow 2002; Carilli & Taylor 2002). Furthermore, there is increasing theoretical evidence based on numerical simulations of structure formation that nanogauss-order fields permeate filaments of the large-scale structure (see, e.g., Ryu et al. 2008). However, an unambiguous detection of the intergalactic magnetic field (IGMF), presumed to exist in cosmic voids, which represent a significant fraction of the volume of the universe, remains elusive. Such a field could be produced, for example, through astrophysical mechanisms such as bulk outflows of magnetized material from radio galaxies (Kronberg 1994; Kronberg et al. 2001), although it is unclear whether such processes could efficiently fill the entire volume of the voids (Zweibel 2006). Alternatively, processes such as the Biermann battery mechanism (Biermann 1950) operating during phase transitions in the early universe could produce the IGMF, provided that its correlation length is sufficiently large to overcome magnetic diffusion and survive to the present day (Grasso & Rubinstein 2001). "Primordial origin" hypotheses such as this one are particularly attractive because the IGMF could then play the role of the seed field necessary in magnetohydrodynamic models commonly invoked to explain the fields observed in galaxies and clusters (Widrow 2002; Kulsrud & Zweibel 2008). Consequently, the detection of the IGMF could provide important insights for solving outstanding problems of its origin and role in both the cosmology and astrophysics of structure formation.

The standard observational technique used to detect weak magnetic fields in galaxies, measuring the Faraday rotation of light from distant quasars, is inadequate for detecting the IGMF for two reasons. The first is that Faraday rotation measures the integrated magnetic field along the line of sight, and therefore the determination of the IGMF relies on the subtraction of the imperfectly measured Galactic magnetic field (Kronberg & Perry 1982; Blasi et al. 1999). The second, and perhaps more important one, is that a sufficiently weak IGMF strength will produce a Faraday rotation measure that is below the resolution limit of currently employed techniques. Existing Faraday rotation measurements place an upper limit of B <10⁻⁹ G on the strength of an IGMF with a correlation length greater than 1 Mpc, and this limit weakens as the correlation length decreases until the limits due to Zeeman splitting measurements of absorption lines in distant quasars become more constraining (as summarized in Neronov & Semikoz 2009). As a result of these two effects, until recently, only upper limits on the IGMF strength have been established.

A new measurement technique has emerged during the past few years that may become a more sensitive tool for the measurement of IGMF characteristics. This technique relies on observations of blazars, believed to be active galactic nuclei (AGNs) whose jets are oriented along the line of sight to Earth, in the gamma-ray energy range from 100 MeV to greater than 10 TeV, and are described by several authors (Neronov & Semikoz 2006; Eungwanichayapant & Aharonian 2009; Dolag et al. 2009; Elyiv et al. 2009). Briefly, teraelectronvolt-scale (TeV-scale) gamma-rays from the blazar interact with the extragalactic background light (EBL), producing an electron–positron pair. The electrons and positrons then undergo inverse Compton (IC) scattering on the cosmic microwave background (CMB) radiation, producing secondary gamma-rays of a lower energy than the primary. As a result, an electromagnetic cascade develops. Because the pairs' trajectories depend on the interactions with the magnetic field, the temporal and angular profiles of this cascade emission at the gigaelectronvolt scale carry information on the strength of the magnetic field in which the cascading occurs. If cascading develops in the voids, then the spectral, temporal, and angular characteristics of the secondary radiation will depend on the IGMF properties. A number of studies have characterized the spectral properties of the secondary photons as well as the temporal profile, commonly referred to as an "echo" (Plaga 1995; Ichiki et al. 2008; Murase et al. 2008), and the angular profile, or "halo," (Aharonian et al. 1994; Neronov & Semikoz 2009; Ahlers 2011). Comparison of the characteristics of secondary gamma-ray radiation with existing data can therefore become the methodology for studying properties of the IGMF.

Most published studies of the subject have focused on the spectral properties of the secondary radiation. For instance, utilizing simple geometric models for the cascade, Neronov & Vovk (2010) and Tavecchio et al. (2010) demonstrated that a lower limit on the IGMF strength can be placed by requiring that the secondary gamma-ray gigaelectronvolt-band emission does not exceed current measurements. In another study, Tavecchio et al. (2011) performed detailed modeling of the spectral energy distribution (SED) of four blazars, thereby reducing the dependence of the conclusions on assumptions about the properties of individual sources. Furthermore, Dermer et al. (2011) relaxed the assumption that the characteristic time to build up the secondary gamma radiation is less than the duration of the very high energy (VHE; energy ⩾100 GeV) activity of the source and derived a less-constraining lower limit on the IGMF. Expanding on these simplified geometric models, Huan et al. (2011) described a semianalytic model employing the energy-dependent distributions of electron and positron energies in the cascade and accounted for the effects of the source lifetime. Monte Carlo simulations were used by Dolag et al. (2011) and Taylor et al. (2011) to confirm and improve the results of simplified geometric semianalytic models. In particular, Taylor et al. (2011) employed a three-dimensional particle-tracking simulation to follow the cascade development and derived a lower bound on the strength of the IGMF that appeared to be consistent for the three blazars studied.

This paper reviews IGMF constraints derived from previous studies, particularly focusing on the conclusions of Taylor et al. (2011), Neronov & Vovk (2010), Tavecchio et al. (2011), Dolag et al. (2011), and Dermer et al. (2011). Because of the scientific importance of these results, we wish to systematically understand all of the ways in which these IGMF constraints may fail. Particular emphasis is given to the effects of major uncertainties, such as the spectral properties of the source, history of VHE activity, geometrical characteristics of the source, and uncertainty due to the poorly known EBL spectral energy density. The conclusions are derived based on the newly developed Monte Carlo code that is described in Section 2 together with models for the EBL, IGMF, and AGN source model. Section 3 describes the data utilized in this paper, and Section 4 provides a detailed analysis of individual sources and a comparison of the results of prior publications. Section 5, the discussion, concludes the paper with a brief review of alternative interpretations of the data.

2. NUMERICAL SIMULATIONS

VHE photons escaping a source such as an AGN jet interact with surrounding diffuse photon fields and generate electromagetic cascades. Cascading occurs in morphologically complex environments of photon and magnetic fields of the host galaxies, the large scale structure filaments, voids, etc. The amplitudes of the magnetic fields and the density of the photon fields in these structures sensitively affect the temporal, angular, and spectral evolution of the cascading, secondary photons.

For example, the highest energy of the escaping photon is determined by the spectral energy density of background photons of the host galaxy. A 30 TeV photon propagating through a Milky Way-like galaxy will have an optical depth of ∼1, based on rough estimates of the energy density of the Galaxy in the far infrared (∼100 μm). Photons with energies higher than this will either be absorbed by interactions with the galactic light or the CMB, on spatial scales less than the size of the galaxy (⩽1 Mpc). These photons will initiate cascades under the influence of galactic magnetic fields (∼10⁻⁵–10⁻⁶ G; see, e.g., Widrow 2002) which are strong enough to isotropize the secondary photons of the cascade (Aharonian et al. 1994).

Photons with energies low enough to escape the galaxy are expected to predominantly interact with the EBL. If this interaction occurs in the environments of either galaxy clusters with magnetic fields of order 10⁻⁶–10⁻⁷ G (see, e.g., Widrow 2002) or intervening large scale structure filaments, with magnetic fields ≳ 10⁻¹² G, then the secondary electrons will be isotropized, thereby dramatically attenuating the observable flux of the secondary photons produced by them. Effects of the electromagnetic cascading may become observable with present day instrumentation when the interaction with the EBL occurs in the voids with characteristic magnetic fields ≲ 10⁻¹² G.

In order to explore in detail the potentially observable spectral, angular, and temporal effects of cascading in the voids of the large-scale structure, we have developed a fully three-dimensional Monte Carlo code. It propagates individual particles of the cascade in a cosmologically expanding universe and accounts for all QED interactions with the EBL and CMB without simplifications. In the following, details regarding these numerical simulations are provided and some of its capabilities are illustrated with several characteristic examples.

2.1. Cosmology and EBL Model

Throughout this paper, a Friedmann–Robertson–Walker cosmology is used with critical matter, cosmological constant, and radiation densities given by Ω_M = 0.3, Ω_Λ = 1.0 − Ω_M, and Ω_R = 8.4 × 10⁻⁵ respectively. The radiation density is negligible for nearby sources of z <1. The standard EBL model chosen in the simulations at z = 0 is modeled with 6 energy density points (I_λ(240 μm) = 15.9 nW m⁻² sr⁻¹, I_λ(60) = 7.3, I_λ(12) = 2.0, I_λ(2.5) = 9.4, I_λ(0.38) = 3.2, I_λ(0.12) = 0.43), and a cubic spline is used to find the EBL energy density at intermediate wavelengths as adopted in Vassiliev (2000). Figure 1 shows a comparison of this EBL model (labeled EBL model 1) with one recently proposed, which is based on various empirical observational data as well as on EBL models existing in the literature (Domínguez et al. 2011). Two additional EBL models are displayed, one (labeled EBL model 2) with a substantially lower dust peak, and another (labeled EBL model 3) with a reduced stellar peak. These two models are explained in detail and used in Section 4.3.

**Figure 1.** Comparison between the EBL model of Domínguez et al. (2011) and the three EBL models used in the present work-a typical EBL model (EBL model 1), one chosen with a low dust peak (EBL model 2), and one chosen with a low stellar peak (EBL model 3).
Download figure:
Standard image High-resolution image

In this study, we investigate sources with redshift <0.2 and thus, we neglect evolution of the EBL in the comoving reference frame; only the effect of cosmological expansion on the EBL energy density and energy of EBL photons is taken into account. By comparison with Domínguez et al. (2011), we do note that evolutionary effects of the EBL may need to be taken into account for sources with z ≳ 0.3.

2.2. Magnetic Field Model

The observational effects of magnetic fields in the gamma-ray signal of extragalactic sources originate in the deflection of the charged particles from the trajectory of the primary photon and subsequent deflection of the direction of the secondary IC scattered photons. Different regimes of influence of the IGMF can be determined by exploring the interplay between the characteristic coherence lengths of the IGMF, the e^± IC cooling length, and the distance from the interaction point to the observer. For the production of secondary photons above 100 MeV, which is of relevance to this paper, the pair-produced e^± should have energies on the scale of ≳100 GeV. In the simulation, electrons are tracked down to energies of 75 GeV. A 1 TeV electron loses its energy to IC scattering on the CMB over a characteristic length of 0.4 Mpc, $\lambda _{{\rm IC}} \propto E_{e}^{-1}$ . For the coherence length of magnetic fields, λ_IGMF ≫ λ_IC, the IC scattering effectively happens in the environment of nearly constant B_IGMF. For λ_IGMF ≪ λ_IC, the deflection of charged particles occurs in the non-coherent regime and leads to smaller deflection. In this study, the coherence length of the IGMF is conservatively chosen to be 1 Mpc, which corresponds to coherent scattering for all photons of interest (with energies >100 MeV). It has been observed that the reversal field length of the magnetic fields in clusters of galaxies is on the scale of 10–100 kpc (Grasso & Rubinstein 2001) and reflects the spatial scales of the distribution of plasma. Thus, the magnetic fields in the voids with significantly larger characteristic plasma distribution scales, are likely to have coherence lengths much larger than this.

The IGMF is modeled in the code as a system of cubic cells with a size equal to the coherence length of the IGMF and magnetic field amplitudes, which are equal in value but randomly oriented in direction. To preserve cosmic variance, each cell is assigned an orientation when the first particle of the electromagnetic cascade propagates through it. If the cascade develops over a large number of these magnetic field cells, the observable effects of the IGMF are randomized. However, if the distance to the observer is comparable to the size of the magnetic field domain, the observational effects of the randomly chosen direction become significant. For this study, we analyze sources at distances greater than a few hundred megaparsecs (z > 0.1). The evolution of the IGMF is unknown, but is likely dominated by cosmological expansion for z < 1. Therefore, the size of the domains and the magnetic field value are evolved with the standard (1+z)⁻¹ and (1+z)² dependences. The values of the magnetic field reported in this paper, refer to the values at z = 0.

Figure 2 illustrates the mean time delay of secondary photons produced by a monoenergetic beam of 100 TeV primary photons at a redshift of z = 0.13. The photons arriving at the observer are integrated over an aperture radius of 10 fdg 0. For each of the energy bins (4 per decade), the mean delay time is computed for six magnetic fields including the zero field case. The time delay for the latter is due to QED scattering of the secondary particles of the pair production and IC scattering processes. For a 10 GeV photon, the time delay amounts to about a half hour and for 0.1 GeV photons, the time delay is about 10 hr. The saturation effect at low energies is due to the aperture cutoff. Figure 2(a) shows the mean time delay with no EBL photons as IC targets in the simulations. It is consistent with (Taylor et al. 2011, Figure 2), and follows the spectrum of T_delay ∝E^−5/2, derived analytically by making several simplifications, as explained in Neronov & Semikoz (2009). Figure 2(b) shows the result when the EBL field is included. The time delay for a non-zero field of secondary photons with energies above 10 GeV is significantly increased, because electrons can move farther from the position of pair production and still scatter higher energy EBL photons toward the observer increasing the average time delay in a given energy bin. The effect can be seen more clearly in Figure 3, which shows the distribution of arrival times of secondary photons in a single energy bin of 10.0–17.78 GeV with and without the target EBL photons for IC scattering.

**Figure 3.** Distribution of arrival times of secondary photons in a single energy bin of 10.0–17.78 GeV at B_IGMF = 10⁻¹⁷ G with and without the EBL photons included as a target for the IC scattering computation.
Download figure:
Standard image High-resolution image

2.3. Gamma Ray Source Model

The gamma-ray source model employed in the simulations is based on leading theoretical speculations about the nature of the TeV blazar source (see, e.g., Urry & Padovani 1995). In the reference frame of the blazar jet, the VHE photons are distributed isotropically with a power law spectrum of index α. Once this distribution is boosted into the reference frame of the host galaxy with a Doppler factor of Γ = (1 − β²)^−1/2, it can be parameterized as

$\epsilon \frac{dF}{d\epsilon }(\epsilon,\theta _{{\rm v}}) = F_o\delta ^{2}\left(\frac{\epsilon }{\epsilon _o\delta }\right)^{-\alpha +1}e^{-\epsilon /\epsilon _c},$

where epsilon is the photon energy in the reference frame of the host galaxy, δ = [Γ(1 −βcos θ_v)]⁻¹, θ_v is the viewing angle from the blazar jet axis to the line of sight of the observer, F_o is a flux normalization factor, and epsilon _o is an energy scale factor. To account for absorption of photons at the highest energies inside the host galaxy, an exponential cutoff at energy epsilon _c is introduced. This model is characterized by five physical parameters, and is sufficient to model the VHE part of the spectrum (≳100 GeV).

Since the energy range of interest for this study covers more than five decades of energy (from 0.1 GeV to >10 TeV) the HE part of the source spectrum (≲100 GeV) is allowed to obey a power law with a different spectral index γ

$\begin{equation} \epsilon \frac{dF}{d\epsilon } = F_0\delta ^{2} \left\lbrace \begin{array}{@{}ll} \left(\frac{\epsilon }{\epsilon _B\delta } \right)^{-\gamma +1} {\rm exp}\left(\frac{-\epsilon }{\epsilon _{c}} \right) & \frac{\epsilon }{\epsilon _B\delta } < 1 \\ \ms \ms \left(\frac{\epsilon }{\epsilon _B\delta } \right)^{-\alpha +1} {\rm exp}\left(\frac{-\epsilon }{\epsilon _{c}} \right) & \frac{\epsilon }{\epsilon _B\delta } > 1 \end{array}\right.. \end{equation} \tag{ 1 }$

The seventh parameter, epsilon _B, introduced in this equation, is the spectral break energy. This multi-parameter gamma-ray source spectrum is given at the redshift of the host galaxy and is necessary and sufficient to satisfy observational data of TeV blazars in both the HE and VHE regimes.

To illustrate the effects of different gamma-ray source model parameters, simulations are shown for z = 0.13, B = 10⁻¹⁵ G, epsilon _c = 30 TeV, α = γ = 1.5. Figure 4 shows the observed spectra of a source with the above parameters when it is viewed at different angles and different jet Doppler factors. All photons arriving within an aperture of 5° around the source are integrated. The black (solid) line shows the spectrum of the prompt photons reaching the observer, which is modified by absorption on the EBL. The prompt photon spectrum is normalized to the same level for different viewing angles and different Doppler factors. Figure 4(a) shows the spectrum of secondary photons when the source is viewed at θ_v = 0°, 2°, 5°, 10° and Figure 4(b) shows the spectrum of secondary photons for Γ = 5, 10, 30, 100. Viewing a source with the same prompt SED, but with increasing viewing angle or Doppler factor implies the power in the jet must increase.

The higher energy end of the secondary photon spectrum (≳20 GeV) is unchanged for different viewing angles and jet Doppler factors because photons at these energies are generated from primary photons leaving the source with nearly zero deflection from the angle to the observer and the flux of these photons is directly proportional to that of the prompt photons. The lower energy end of the spectrum (≲20 GeV) is generated by secondary electrons of lower energies, the trajectories of which are significantly deflected from that of the primary photon producing them. Additionally, the cooling length due to IC losses increases inversely proportionally to energy, allowing significantly larger deflection angles. The position of the peak of the SED is correlated to the value of the magnetic field, and its intensity is proportional to the overall power in the jet which increases with larger viewing angle or Doppler factors. When the characteristic angular size of the jet becomes larger than the viewing angle (1/Γ > θ_v), the SED is nearly independent of the Doppler factor.

Figure 5 displays a simulation of the photon arrival distribution from a source at z = 0.13, with gamma-ray source model parameters α = γ = 1.5, epsilon _c = 30 TeV, B_IGMF = 10⁻¹⁵ G, Γ = 30, at four different observing angles, θ_v = 0°, 2° 5°, and 10°. The main trend in these figures is that the overall luminosity of prompt and secondary emission rapidly declines as the observing angle increases, and the photon distribution around the source becomes increasingly axially non-symmetric, when θ_v ∼ 1/Γ. The luminosity of the secondary photons relative to the prompt emission rapidly declines with increasing viewing angle, thus, detecting non-axially symmetrical halos around AGN with existing instrumentation may prove challenging. Figure 5 is in qualitative agreement with previously reported findings in the study of these effects by Neronov et al. (2010).

2.4. Particle Propagation and Interaction

The development of cascades in intergalactic space initiated by VHE photons is modeled by full three-dimensional ray tracing and interaction of all particles: electrons and photons. The processes of pair production and IC scattering are treated with the use of full QED cross sections (Jelley 1966; Gould & Schréder 1967) implemented without simplification. For example, both the EBL and CMB are included in IC scattering as seed photon fields and scattering includes the Klein–Nishina regime. Both pair production and IC scattering are treated in the code through the sampling of marginal probability density functions. For example, the mean free path λ of electrons of energy E due to IC scattering is given by

$\begin{equation} \frac{1}{\lambda } = \frac{1}{2}\left(\frac{m_e^2}{2E}\right)^2\frac{1}{\beta }\int _0^\infty d\epsilon \frac{n_\epsilon (\epsilon)}{\epsilon ^2}\int _{\frac{2E\epsilon (1-\beta)}{m_e^2}}^{\frac{2E\epsilon (1+\beta)}{m_e^2}}dx\sigma _{\gamma e}(x)x, \end{equation} \tag{ 2 }$

where β is the speed of the electron, σ_γe(x) is the differential IC cross section, $x= (2E\epsilon /m_e^2)(1-\beta \cos \theta)$ where θ is the collision angle and n( epsilon ) is the spectral energy density of the isotropic seed photon field. The code samples the propagation length of the particle with the use of the mean free path, λ to determine the position of the interaction. It then generates the marginal probability density for an interaction of the electron with a given energy photon and samples the energy of the interacted photon. The process is continued by sampling the interaction angle of a given energy photon by generating a marginal probability density of the angular distribution of seed photons. The azimuthal angle is sampled randomly from a uniform distribution. At the completion of this process, the kinematics of the interaction is fully determined, allowing the computation of the energy-momentum vector of the outgoing photon and electron. The pair production process is treated similarly.

Propagation of both photons and electrons is simulated in a cosmologically expanding universe. Energy losses of photons are only due to cosmological expansion, while energy losses of electrons also include IC cooling. All photons with energies above 100 MeV are tracked until they reach the z = 0 surface, at which point the position and momentum 4-vectors are saved. All electrons are tracked until their energy decreases to less than 75 GeV, or until they reach the surface of z = 0. Particular care is given to the computation of time delays. The time delay of individual particles is computed relative to the arrival time of a putative photon propagating directly from the source to the current position of the particle. This procedure is adopted to maintain precision in numerical simulations down to minute timescales. The computation accuracy of time delays was verified using an arbitrary precision numerical integrator developed at Lawrence Berkeley National Laboratory.⁴ The full solution for the equations of motion of electrons in a cosmologically expanding universe and under the influence of constant magnetic field are used to track the position and momentum 4 vectors between IC interactions.

Previous models ranging from analytic to Monte Carlo codes have been reported in the literature to generate constraints on the IGMF (Neronov & Vovk 2010; Tavecchio et al. 2010; Dermer et al. 2011; Essey et al. 2011a; Dolag et al. 2011; Taylor et al. 2011; Huan et al. 2011). Each one has utilized various degrees of simplification through the use of solutions of one or two-dimensional kinetic equations or simplified analytical approximations of QED interactions or geometrical effects of cascade development. For example, the VHE secondary photons produced by the highest energy electrons are capable of generating second and higher orders of the EM cascades, which are typically neglected in analytic and semi-analytic codes, whereas the majority of Monte Carlo codes used in prior studies of this subject, have neglected a full three-dimensional simulation.

Figure 6 illustrates the importance of higher order cascading to correctly describe secondary photons with energies ≳ 200 GeV. This figure shows the results of the simulations of secondary emission produced for a source at z = 0.3, with B_IGMF = 10⁻¹⁶ G, θ_v = 0, and Γ = 10. The two panels display different intrinsic spectra for the source; panel (a) is obtained with α = 1.5 and epsilon _c = 30 TeV, and panel (b) corresponds to α = 1.5 and epsilon _c = 5 TeV. Both panels show the direct differential flux energy density (dFED) of the source together with the time-integrated secondary dFED, within 0 fdg 5 from the position of the source. The two lines of the secondary dFED shown on the figure are obtained with second order cascading included or not included in the simulation. It is evident that the secondary dFED ≳ 200 GeV is strongly affected by this assumption and if this simplification is made, the Monte Carlo simulations may over-predict the total dFED in the VHE energy range.

fdg — **Figure 6.** Comparison of spectra produced by both neglecting and including multiple generations in the cascade for a source at z = 0.3, with B_IGMF = 10⁻¹⁶ G, θ_v = 0, and Γ = 10. (integrated over a 05 aperture). The black (solid) line is the prompt, direct spectrum, the red (long dashed) line shows the secondary spectrum produced when multiple generations of the cascade are included, and the green (small dashed) line shows the secondary spectrum when only a single generation of cascading is allowed in the simulations. Two different source spectra were used (a) spectral index α = 1.5, cutoff energy, _c = 30 TeV and (b) α = 1.5, _c = 5 TeV.
Download figure:
Standard image High-resolution image

3. GAMMA-RAY DATA, INSTRUMENTATION, AND STRATEGY FOR DATA ANALYSIS

The gamma-ray data used in this study are compiled from the data sets of two types of instruments. In the HE regime, the data were obtained by the Fermi-Large Area Telescope (LAT) instrument and the data were processed utilizing publically available tools, version v9r23p1, with the update from 2011 November 6. For the VHE regime, data previously reported by the Imaging Atmospheric Cherenkov Telescopes (IACT) instruments VERITAS and HESS were used. In both the HE and VHE energy regimes, the software used for analysis interprets the gamma-ray flux as originating from a point source; no angular extension is assumed. Since the secondary photons of intergalactic cascades inherently represent an extended source, the methods of comparison of simulations and data are instrument-specific.

3.1. VHE Data Considerations

In most of the following discussion, the validity of the B_IGMF = 0 hypothesis (hereafter called H0) is examined. Under this assumption, the source of the gamma-ray photons has an angular extent due to the QED pair production and IC scattering processes. The angular size of this extension is much smaller than the gamma-ray point spread function (PSF) of both IACTs and the Fermi-LAT instrument, and therefore, the gamma-ray sources are point-source like. In the VHE regime, the point source assumption becomes invalid for B_IGMF ≳ 10⁻¹⁵ G because the PSF of IACTs and the angular extension of the source become comparable. This requirement, however, is dependent on the higher energy end of the spectral energy density (≳10 TeV) of the source and its history of activity.

To process simulated data in the VHE energy regime, the instrumental PSF of

$\begin{equation} \frac{dP}{d\theta ^2} = \frac{(1-\alpha)}{\theta _1^2} {\rm exp}\left(-\frac{\theta ^2}{\theta _1^2} \right) + \frac{\alpha }{\theta _2^2} {\rm exp}\left(-\frac{\theta ^2}{\theta _2^2} \right) \end{equation} \tag{ 3 }$

is used as suggested in Aharonian et al. (2006c). The typical values of parameters of the PSF are θ₁ = 0 fdg 06505, θ₂ = 0 fdg 1697, α = 0.505. These parameters are assumed to be energy independent in analyses of IACT data reported and therefore, we make a similar assumption for processing simulated data. For each simulated photon, Equation (3) is used to find the probability for the given photon to be reconstructed within the 68% containment radius from the position of the putative point source. The effective point source flux from an AGN is estimated as the total simulated flux within the 68% containment radius divided by 0.68. This method of flux evaluation for each energy bin accurately models the IACT data as reported in the literature.

The VHE data set used in this paper is summarized in Table 1. The data sets of the first three sources (RGB J0710+591, 1ES 1218+304, and 1ES 0229+200) are identical to those used in Taylor et al. (2011), except that an additional data set for 1ES 1218+304 obtained by the VERITAS Collaboration just before the launch of the Fermi satellite (on 2008 August 4 or MJD 54682) is considered in this study. As reported in Acciari et al. (2010b), the activity of the source is nearly identical during these non-overlapping periods, except for an elevated flux of the source peaking at the level of ∼20% Crab over a few nights of observations. The data set for 1ES 0229+200 was also obtained prior to the launch of the Fermi satellite. Based on the report from Perkins & VERITAS Collaboration (2010), the activity of the source as measured by VERITAS during the second year of the Fermi mission, appeared to resemble the reported SED by the HESS collaboration prior to the launch of Fermi satellite (Aharonian et al. 2007b). VERITAS has continued monitoring this source since the Fermi launch and tentatively detected flux variations on a yearly timescale (J. S. Perkins & VERITAS Collaboration 2012, private communication). Finally the data sets for four other extreme TeV blazars (1ES 0347−121, 1ES 1101−232, H 2356−309, and RGB J0152+017) were taken from the discovery publications by the HESS collaboration, and all of these sources were observed prior to the start of the Fermi mission.

Table 1. Data of RGB J0710+591 are Taken by VERITAS Acciari et al. (2010a)

IACT Data Summary
Source	z	IACT	Flux (10⁻¹² cm⁻² s⁻¹)	Index	MJD (Approx)	hr
RGB J0710+591	0.125	VERITAS	(>300 GeV) 3.9 ± 0.8	2.7 ± 0.3	54882–54892	22.1
1ES 1218+304	0.182	VERITAS	(>200 GeV) 18.4 ± 0.9	3.1 ± 0.3	54829–54944	27.2
1ES 1218+304	0.182	VERITAS	(>200 GeV) 12.2 ± 2.6	3.1 ± 0.1	54070–54220	17.4
1ES 0229+200	0.14	HESS	(>580 GeV) 0.94 ± 0.24	2.5 ± 0.2	53614–53649, 53967–54088	41.8
1ES 0347−121	0.188	HESS	(>250 GeV) 3.91 ± 1.1	3.1 ± 0.3	53948–54100	25
1ES 1101−232	0.186	HESS	(>200 GeV) 4.5 ± 1.2	2.9 ± 0.2	53111–53445	43
H 2356−309	0.165	HESS	(>200 GeV) 4.1 ± 0.5	3.1 ± 0.3	53157–53370	40
RBG J0152+017	0.08	HESS	(>300 GeV) 2.7 ± 1.0	2.9 ± 0.5	54403–54418	15

Notes. Source 1ES 1218+304 was observed by VERITAS and the results were published in two separate submissions (Acciari et al. 2010b, 2009). This source was also detected during a six day period around MJD 53750 by the MAGIC collaboration Albert et al. (2006), at a level similar to that of the reported flux by VERITAS about one year later. Observations of 1ES 0229+200, 1ES 0347−121, 1ES 1101−232, H 2356−309, and RGB J0152+017 were taken prior to the launch of the Fermi satellite, and are reported by the HESS collaboration (Aharonian et al. 2007a, 2007b, 2007c, 2006b, 2008).

Download table as: ASCII Typeset image

3.2. HE Data Considerations

For an AGN with significant power emitted at E ≳ few TeV, the angular extent of the source due to cascade emission may become comparable to the PSF of the Fermi-LAT at fields as low as B_IGMF ≳ 10⁻¹⁷ G. For B_IGMF magnitudes less than this, an AGN is effectively a point source for the LAT. To evaluate the effective point source flux of an AGN for the Fermi-LAT, a procedure similar to that of the IACT case was adopted. For every simulated photon of a given energy, the energy-dependent PSF of the Fermi-LAT was used to determine the probability of reconstruction of this photon within the 68% containment radius from the position of the putative point source. The flux evaluated within the 68% containment was again rescaled by 0.68⁻¹ to estimate the effective point source flux in each energy bin. The PSF used for this conversion was determined from a two year time-averaged sample of AGN.⁵

To compare simulated differential fluxes from an effective point source to the Fermi-LAT data, the standard analysis tools were applied but with a notable important distinction from previous studies, which are cited in Section 2.4. Since in the HE regime, the flux of gamma-ray photons can be dominated by either prompt or secondary emission, we first derive from simulations, the spectral index in each energy bin. This index is then used as a fixed parameter for the maximal likelihood evaluation⁶ of the flux in each energy bin in the Fermi data, within the 10° region of interest which also includes all nearby sources from the two year point source catalog and diffuse backgrounds. The HE point source fluxes or upper limits are derived using this procedure.

3.3. Data and Model Comparison

To compare the predictions of the Monte Carlo simulations to the data, we combine both the HE and VHE regimes to compute a χ²-like parameter, as further explained in this section. The simulated effective point-source flux is used as the model expectation value, and we assume that the statistical error of the vast amount of simulations is negligible compared to the observational error. The point-source fluxes derived with the use of the Fermi tools as explained above or obtained from IACT publications are used as data. The observational error obtained or reported is taken to be the primary source of discrepancy between the data and the model. A χ²-like parameter is used to estimate the goodness of fit of a given model and χ² statistics is used to convert it to a confidence level (or the probability of the model exclusion). For this conversion, an explicit assumption is being made, that the errors in both the HE and VHE regimes are dominated by statistics with a Gaussian distribution. Throughout the paper, our sole goal is to test the null hypothesis, H0, that B_IGMF = 0 is incompatible with the data and the simulations. The confidence level is the probability that the measurement cannot be obtained with the assumption of the given model. In what follows, we take the model to be incompatible if the confidence level exceeds 95%, corresponding to a 2σ deviation for a normal distribution. We do not claim any meaningful interpretation of a higher confidence level, due to the unknown behavior in the tails of the distribution of errors of each instrument.

The model spectral energy density is derived based on the full Monte Carlo simulations computed for monoenergetic primary photons with eight bins per decade, of equal width in logarithmic space. The simulated spectral energy density data are equally binned with eight bins per decade and each bin is centered on the energy of the primary photon monoenergetic line. The VHE data are reported in different publications at different energies and with different binning. We use simulated data to interpolate the flux value to the reported positions of the bins and their widths. In the three decades of the HE regime, four to six energy bins are generated depending on the given source luminosity and statistics. The simulated data are then used to interpolate the flux value to these energy bins. Moreover, we use the simulated data to find the spectral index for the power-law distribution of photons at each energy bin. This spectral index is taken to be fixed when we find the flux value and its error utilizing the Fermi tools.

To compare the HE and VHE data of each source with the simulations, the effective point source fluxes are generated for a set of gamma-ray source models with fixed values of four parameters—α, epsilon _c, Γ, and θ_v. For each model, α is chosen in the range (1.0, 2.5) with a step of 0.1, and epsilon _c is chosen in the range (600 GeV, 60 TeV) with eight bins per decade, equally spaced in logarithmic energy. Each model is characterized with default values of θ_v = 0, and Γ = 10, unless otherwise stated. Therefore, for each source, we test 240 individual source models.

The remaining three parameters of the seven parameter source model were determined as follows. Two of these parameters, the break energy epsilon _B and the spectral index γ are relevant only to the HE part of the spectrum. They were chosen by minimizing the $\chi _{\rm HE}^2$ value, by allowing epsilon _B to vary between ≳10 GeV to the lower edge of the lowest energy VHE data bin and γ between −5 and 5. This interval for the break energy is motivated by the fact that IACT instruments become insensitive in this energy regime and Fermi runs out of photon statistics, therefore allowing a possible knee feature in the spectrum to be undetectable. The spectral index γ in the HE regime may or may not be constrained by minimization of the $\chi _{{\rm HE}}^2$ value. It is evident that when secondary, cascade emission dominates in this energy regime, it is sufficient for γ to be larger than some value to keep prompt emission negligible. The flux normalization factor, F₀, is the only optimization parameter for a given model which is relevant to the $\chi _{{\rm VHE}}^2$ and which may or may not be relevant to the $\chi _{{\rm HE}}^2$ , depending on the relationship between the prompt and secondary emission of the source. We chose to determine F₀ by minimization of the VHE part of χ² by solving $\partial \chi _{{\rm VHE}}^2 / \partial F_0$ = 0 independently from the behavior of $\chi _{{\rm HE}}^2$ . Effectively, this means that we have made a stringent requirement of compatibility of the source model with the VHE data and have excluded some models which would be highly incompatible with the VHE measurements, but would allow statistical compatibility with an overall χ² = $\chi ^2_{{\rm HE}}$ + $\chi ^2_{{\rm VHE}}$ , just because of an increased number of data points and therefore number of degrees of freedom. We view this weighting procedure of HE and VHE parts of χ² in determination of F₀ as better physically motivated, since the highest energy data points in the VHE regime are of extreme importance for the production of the cascade emission, but from a statistical point of view, they are equal to any other point of N_HE or N_VHE measurements.

The χ² value for each model with four fixed parameters (α, epsilon _c, Γ, and θ_v) was converted into a confidence level, using χ² statistics with N_HE − 2 + N_VHE − 1 degrees of freedom, assuming that three parameters were optimized for each model (F₀, epsilon _B, γ). We make no attempt in our studies to evaluate the confidence intervals of the latter three parameters of each model. Our goal is exclusively to find a model or a set of source models which are compatible with H0.

As an illustration of a typical result of the data and model comparison, Figure 7(a) shows the χ² confidence level of the dFED, assuming four fixed and three free parameters for each model tested within the given range of α and epsilon _c parameters. The most favored model with value of α = 1.8 and epsilon _c = 3.16 TeV is incompatible with the data at the level of about 75%. If this choice of parameters α and epsilon _c is considered as an optimization process, in which case the number of free parameters in the model is 5, then this model is incompatible with the data at the level of 88%. Figure 7(b) shows the data points for the dFED and the best fit simulation result with these α, epsilon _c parameters. For this source's dFED, below 10 GeV, the spectrum is dominated by cascade emission, while above 10 GeV, it is dominated by prompt radiation.

4. IGMF CONSTRAINTS AND EFFECTS OF SYSTEMATIC UNCERTAINTIES

The primary goal of this paper is to provide precise numerical verification of the constraints on the IGMF as reported in Neronov & Vovk (2010); Tavecchio et al. (2010); Dolag et al. (2011); Dermer et al. (2011); Taylor et al. (2011), since most of these results were obtained with various simplifications in the analysis, and some were derived with semi-analytical approaches to qualitatively verify more detailed computations (Dermer et al. 2011). Particular focus is given to establishing the robustness of magnetic field limits when various systematics are taken into consideration.

Three sources (RGB J0710+591, 1ES 1218+304, and 1ES 0229+200) are used in this study which have been reported as having provided constraints on the IGMF, in the comprehensive study of Taylor et al. (2011 and references within). Hereafter, we refer to this paper as TVN11. An additional four hard-spectrum TeV blazars observed prior to the beginning of the Fermi mission are considered—1ES 0347−121, 1ES 1101−232, H 2356−309, and RGB J0152+017. The Fermi-LAT data for these sources are re-analyzed and are collected from the mission start time 2008 August 4 to 2012 February 14, and the updated P7SOURCE IRFs are used along with the Pass 7 data. The models for extragalactic and diffuse backgrounds were used together with the standard gamma-ray selection constraint of zenith angle <100° which eliminates earth limb gamma-rays.

Perhaps the main source of uncertainty in placing constraints on the IGMF stems from the unknown duty cycle of TeV blazars and particularly, the history of the highest energy TeV emission, as has been pointed out in Dermer et al. (2011). The sampling of the VHE activity of these sources reported by IACTs is limited to a few tens of hours dispersed over a period of a few weeks to a few years. In the regime of very low IGMF (B_IGMF < 10⁻²⁰ G), most of the secondary radiation from intergalactic cascades with energy >100 MeV, which originates from the primary VHE flux sampled by IACTs, would have reached the earth and would be detected by the Fermi-LAT (see Figure 2) within a few hours. This assumes that the HE flux from a given source sampled by the Fermi-LAT over the period of the mission (about four years) could be viewed as "contemporaneous" to IACT measurements, for the purposes of verification of H0. This explicitly assumes, however, that the duty cycle of a given source in the VHE regime = 1 over this same period.

4.1. Analysis of RGB J0710+591 Data

The VHE observations of RBG J0710+591 are summarized in Table 1. They include five energy data points reported by the VERITAS Collaboration in Acciari et al. (2010a), observed during the time period 2008 December–2009 March for a total of 22 hr. For the HE regime, six energy bins were used spanning from 200 MeV to 200 GeV. In each energy bin, if TS > 9, the flux point and 1σ error bars are displayed. Otherwise, an upper limit was computed. It is important to note that with this strategy and extended data set as compared to previously used in TVN11, the flux in the lowest energy bin now constitutes a flux point rather than an upper limit. This data point had been critical for rejecting H0.

Figure 7(a) shows the confidence level for rejecting H0, for a set of models characterized by the range of epsilon _c and α described in Section 3.3. It identifies the best fit model with values of α = 1.8 and epsilon _c = 3.16 TeV, which is incompatible with H0 at <75% confidence level assuming three free parameters (and <88% assuming five free parameters). The range of models in the vicinity of this point is not incompatible with H0. The best fit of the simulated dFED and observations is shown in Figure 7(b). It appears that the conclusion of TVN11 that H0 is ruled out at the 98.8% level is invalidated due to two factors. First, the Fermi-LAT data set underwent revision from the old Pass 6 version (P6) to the current Pass 7 (P7) and this allowed a detection to be made in the lowest energy bin, which is higher than the previously computed upper limit. The second factor may be due to the fitting algorithm applied in this present work which is different from that adopted in TVN11, in which the flux or upper limit determination in a given energy bin was fixed to the best fit index over the entire energy range.

Furthermore, the geometrical orientation of the jet with respect to the observer and the jet boost factor represents another source of uncertainty, and tuning these parameters can further improve the goodness of the χ² fit. For example, TVN11 assumes a viewing angle of 2° and an effective jet opening angle of 6°, corresponding to a boost factor of ∼10. As shown in Figure 4, however, lower boost factors or smaller viewing angles lead to lower total power of the jet at the highest energies, and therefore lead to reduced secondary flux.

4.2. Analysis of 1ES 1218+304 Data

The VHE observations of 1ES 1218+304 are summarized in Table 1, which includes two data sets. The first set, obtained during 2008 December–2009 May, is based on 27 hr of data and has nine energy data points reported by the VERITAS Collaboration in Acciari et al. (2010b). This data set was used in TVN11, and for consistency it is also used in this study. The Fermi-LAT data for this source were produced in the same way as for RGB J0710+591, with six energy bins spanning from 200 MeV to 200 GeV. This source has excellent statistics in each Fermi-LAT energy bin, with TS > 25. It is important to note that the extended exposure and updated Pass 7 data set used in this work shows no statistically significant difference compared to that of TVN11.

Figure 8(a) shows the confidence level for rejecting the H0 hypothesis on the α– epsilon _c plane. It suggests a best fit model with values of α = 1.8, epsilon _c = 3.16 TeV, a spectral break energy of epsilon _B = 10 GeV, and an index below the break energy of γ = 1.0, which is incompatible with H0 at the less than 65% confidence level, assuming three (and less than 80% assuming five) free model parameters. The fit of the simulated dFED and observations is shown in Figure 8(b). It is evident that the conclusion of TVN11, that H0 is ruled out with more than 99.99% probability, is purely due to the assumption that a single source spectral index holds for over five orders of magnitude in energy. Allowing a spectral break energy and an intrinsic spectral index below the break energy to vary as detailed in Section 2.3 makes it possible to interpret the 1ES 1218+304 data set as consistent with H0.

The amount of secondary radiation strongly depends on the power output of the TeV blazar at the highest energies. For this source, more so than the others, IACT observations demonstrate strong variability in the VHE spectrum. The VERITAS Collaboration reports that the 1ES 1218+304 data were sampled sparsely over a period of 115 days in late 2008–2009, and while the majority of the data are consistent with a steady baseline flux, the data set also includes a statistically significant flare which peaked at ∼20% Crab, and lasted a few nights. The flux at the peak of the flare was three to four times higher than the baseline flux and it significantly increases the average flux value observed over the entire period. Furthermore, evidence for variability of this source can be inferred from the VERITAS publication covering its two year activity which occurred prior to the Fermi mission (Acciari et al. 2009). The flux observed at that time constitutes about 60%–70% of that reported in the second data set (see Table 1). Overall, the IACT data to date suggest that the average observed VHE flux of this source could be lower than used in TVN11, yet the assumption of the higher average VHE flux is still compatible with H0.

4.3. Analysis of 1ES 0229+200 Data

The parameters of the data set of 1ES 0229+200 are given in Table 1 and include two sets of observations by the HESS collaboration during the period from 2005 September 1 to 2006 December 19 accumulating 41.8 hr exposure (Aharonian et al. 2007b). This data set provides the time-averaged dFED for 8 bins over the energy range spanning from 500 GeV to 16 TeV. This same data set was used in the previous study of TVN11. The Fermi-LAT dFED for 1ES 0229+200 utilizes four evenly spaced bins in log space in the range from 420 MeV to 300 GeV. Only the first energy bin in this data set provides a strong detection (TS >25), for all simulated models in which the secondary flux dominates the total flux. The TS for all other energy bins is typically found at >9 for the majority of simulated models but in some cases, only the upper limit can be established (TS <9). This indicates that the source detection in the HE regime is weak. Perhaps more so than for any other source, the dFED of 1ES 0229+200 does not resemble a power law in the HE energy regime, making the χ² fits relatively poor.

To evaluate the goodness of fit of the data to the Monte Carlo simulations, we assume that the data accumulated by the HESS collaboration over 2005–2006 is representative of the source activity during the first 3.5 yr of the Fermi-LAT data used in this work. The χ² fit obtained under this assumption and for H0 is shown in Figure 9(a). We confirm the result of TVN11 that this source does not have a viable source model that explains the combined HE–VHE data set and agree that H0 is ruled out at the 99.5% confidence level. The best fit (α = 1.3, epsilon _c = 1 TeV) model requires a dramatic spectral break just below 100 GeV and the dFED of 1ES 0229+200 below this energy is completely dominated by secondary flux as shown in Figure 9(b).

**Figure 9.** Source 1ES 0229+200 analyzed under the assumption of H0. (a) Confidence level, or the probability of exclusion of the gamma-ray source model (with fixed α and _c), where the remaining three parameters of the model (γ, _B, and F₀, see Equation (1)) are chosen so as to minimize χ². (b) Simulated dFED of the best-fit model (C.L. = 0.995) of α = 1.3, _c = 1 TeV, showing both the prompt and secondary cascade contributions to the total dFED, along with the HE and VHE data.
Download figure:
Standard image High-resolution image

The effects of two systematic uncertainties, viewing angle θ_v and Doppler factor Γ, were considered. The default assumptions in the analysis are θ_v = 0° and Γ = 10. Increasing the viewing angle, e.g., θ_v = 2° as used in TVN11, would further overproduce radiation in the HE regime (see Figure 4(a)). Increasing the Doppler factor, Γ, combined with the θ_v = 0° assumption would imply that the overall luminosity of the jet in the VHE band should be lower to fit the VHE observations since the jet is collimated into a smaller angle. Rescaling of the prompt radiation to fit the VHE data will however equally rescale secondary emission in the HE regime if the angular distribution of the prompt photons is significantly wider than the characteristic scattering angles acquired in the QED processes. Therefore, for reasonable values of Γ < 100, the change in the secondary radiation above 100 MeV for the best fit model was found to be negligible. The effect becomes of order 10% in the lower part of the HE spectral range only for Γ ∼10⁴–10⁵, and it cannot be used to reconcile the observational data of 1ES 0229+200 with H0.

Since the energy density of the EBL directly affects the propagation length of VHE photons, the uncertainty in the EBL model represents perhaps the most important source of systematic error. To investigate this, two additional EBL models were generated, EBL model 2 and 3, shown in Figure 1. EBL model 2 is characterized by a considerably lower energy density in the far infrared peak of the dust emission. This model is motivated by the recently resolved lower limit on the EBL which is based on the galaxy counts in the data obtained with the Spitzer (Béthermin et al. 2010; Dole et al. 2006), Herschel (Berta et al. 2010), and AKARI (Matsuura et al. 2011) satellites. This model corresponds to the lowest possible far infrared EBL energy density allowed within 2σ. It was found that even with such an extreme assumption about the far infrared EBL, the decrease of the secondary radiation in the model of 1ES 0229+200 was negligible. This conclusion is due to the fact that the source models providing the best χ² fits have HE cutoffs, epsilon _c, below ∼5 TeV, and are thus insensitive to EBL photon wavelengths ≳ 25 μm (kinematic threshold of pair production).

The EBL model 3 is based on the resolved EBL energy density of the starlight peak, which is derived from galaxy counts utilizing data from the Hubble Space Telescope (Madau & Pozzetti 2000) in the visible (from 0.36 to 2.2 μm), Spitzer in the near-IR (3.6–8 μm; Fazio et al. 2004), and ISO in the mid-IR (15–24 μm; Elbaz et al. 2002; Papovich et al. 2004). As compared to the default model 1, this EBL model has an energy density in the visible reduced by about 25% which is compatible with the galaxy counts results to within 1σ. The energy density in the near IR is reduced by about 50%. At 3.6 μm, model 3 has an energy density of 4.0 nW m⁻² sr⁻¹ which is within the 2σ limit (3.5 nW m⁻² sr⁻¹) from the galaxy counts result derived by Fazio et al. (2004). At 4.5 μm, model 3 has an energy density of 3.0 nW m⁻² sr⁻¹, which is within the 1σ limit from the Fazio et al. analysis. Given the uncertainty in the 5.8 and 8.0 μm galaxy counts by Fazio et al. at the bright fluxes, Franceschini et al. (2008) reanalyzed the Spitzer data to conclude that at 8.0 μm the energy density of the EBL is 1.92 nW m⁻² sr⁻¹ with a 2σ lower bound of 1.23 nW m⁻² sr⁻¹. The energy density of EBL model 3 is 1.4 nW m⁻² sr⁻¹ which is within the 2σ bound from the Franceschini et al. (2008) result. As has been pointed out by several authors (Mazin & Raue 2007; Kneiske & Dole 2010; Domínguez et al. 2011) the 5.8 μm result of Fazio et al. (2004) is likely to have been also contaminated by excessive contributions from bright local galaxies, but it has not yet been re-analyzed, unlike the 8 μm point. Nevertheless, all of these authors have recognized that that the Fazio et al. result at 5.8 μm should be corrected, and many have used the EBL energy density at this point, significantly lower than 3.6 nW m⁻² sr⁻¹ reported by Fazio et al. (2004). The 1σ lower bound used by Mazin & Raue (2007) is 2.4 nW m⁻² sr⁻¹ and it is 2.5 nW m⁻² sr⁻¹ in Domínguez et al. (2011). The 2.1 nW m⁻² sr⁻¹ assumed in model 3 is within the 2σ error bar from these later results. Thus, EBL model 3 is compatible with the galaxy counts results to within 2σ but it effectively does not allow any additional contribution to the EBL from unresolved or unknown sources.

The results of the simulations of intergalactic cascading for model 3 with H0 are shown in Figure 10(a). It was found that a number of 1ES 0229+200 source models are compatible with the combined VHE and HE data set, and one of the best fit examples characterized by α = 1.6, epsilon _c = 3.16 TeV is shown in Figure 10(b).

The strong sensitivity of the secondary photon flux to the EBL energy density in the near-IR combined with the conspicuous lack of a non-trivial EBL absorption feature in the VHE energy band (∼200 GeV–5 TeV) is due to the peculiar behavior of the EBL in this wavelength range. For λI_λ ∝ λ⁻¹, the optical depth is independent of the energy of a VHE photon (gray opacity). A small deviation from this proportionality results in a logarithmically slow dependence of the optical depth on the energy of the VHE photon, producing a power law, rather than exponential-like change in a blazar spectrum as discussed in Vassiliev (2000). The behavior of the SED in the mid-IR exactly satisfies this condition and explains the "invisibility" of EBL absorption effects in a blazar dFED. The effect however is strong and reflected in the change of the spectral index of this blazar. In fact, this feature was used to derive upper limits on the EBL energy density in the mid-IR, which are taken to be the values of model 1, by assuming that the spectral index of the source 1ES 1101−232 cannot be harder than 1.5 (e.g., Aharonian et al. 2006a). The 25%–50% lower mid-IR density of model 2 significantly softens the intrinsic spectrum of 1ES 0229+200 reducing the total energy available for the development of the intergalactic cascade, and therefore the flux of the secondary photons. Thus, there are a range of EBL models with mid-IR energy density bounded by the lower limits on the EBL to some SED slightly below that of model 1 which are compatible with the EBL lower limits and H0.

In a recent study of the dual constraints on the EBL and IGMF it was found that an EBL model similar to model 3, is still incompatible with H0 and would require a lower bound of B_IGMF = 6 × 10⁻¹⁸ G (Vovk et al. 2012). All source models analyzed in that paper had a single cutoff energy epsilon _C = 5 TeV and varying spectral index α in the range of 0–1.5 with a single power law dFED over the entire VHE and HE energy range. A similar assumption of a steady flux from 1ES 0229+200 over the lifetime of the Fermi-LAT was made. Figure 11 illustrates the confidence level for a wider range of source models analyzed in this work together with the range of models considered by Vovk et al. (2012). The figure confirms the incompatibility of H0 with the data given assumptions about the source model used in that paper. However, in an extended parameter space of the 1ES 0229+200 models, H0 can be reconciled with observations of this source.

**Figure 11.** Confidence level, or the probability of exclusion of the gamma-ray source model (with fixed α and _c), where the remaining three parameters of the model (γ, _B, and F₀, see Equation (1)) are chosen so as to minimize χ², for source 1ES 0229+200 obtained under the assumption of EBL model 3. The black line terminated with crosses represents the range of models analyzed in Vovk et al. (2012).
Download figure:
Standard image High-resolution image

Another avenue to make the HE–VHE observational data of 1ES 0229+200 compatible with H0 is to question whether or not the VHE flux of the source reported is representative of the average flux during the Fermi observations. The HESS measurements were accumulated in 2005 (6.8 hr)–2006 (35 hr) and are not strictly contemporaneous with the Fermi HE spectrum. The significance of the detection in 2005 reported is 2.7σ, while in 2006 it is 6.1σ with average photon fluxes above 580 GeV of 6.8 × 10⁻¹³ cm⁻² s⁻¹, and 10 × 10⁻¹³ cm⁻² s⁻¹, respectively. Due to the low flux of this source (1%–2% of the Crab nebula flux), and the small data set in the original 1ES 0229+200 discovery paper, statistically significant flux variability as observed in 2005 and 2006 was not detected. Although these observations are compatible with a constant flux, the variability hypothesis cannot be ruled out, based on the statistical and systematic errors reported. Furthermore, observations of this source in 2009, as reported by VERITAS (Perkins & VERITAS Collaboration 2010), were compatible with the average flux value of the HESS data set. However, the average flux obtained was dominated by a period of significantly higher "flaring" activity during a single dark run. In general, VHE observations above a few TeV (relevant for secondary photon production) require considerable integration time and so far, they are too sparse to claim that the HESS value of the flux is representative of the average flux during the Fermi mission. Further communication with the VERITAS Collaboration suggests that the flux level of 1ES 0229+200 has been steadily declining from 2009 to 2012 (J. S. Perkins 2012, private communication). To investigate the effect of a reduced duty cycle for 1ES 0229+200, the spectrum of this source was modified at the highest five energy points to half an order of magnitude of their reported values. It was found that the VHE–HE data set combined in this way does not rule out H0 at more than 95% confidence level. One of these compatible models, with α = 1.3 epsilon _c = 1.78 is shown in Figure 12. Therefore, the conclusion that the H0 is ruled out with high significance heavily rests on the assumption that the HESS measurements are representative of the average flux for E ≳ 2 TeV and a ∼10^−1/2 change of the highest energy part of the spectrum invalidates this conclusion.

**Figure 12.** Spectral model of 1ES 0229+200 of α = 1.3, _c = 1.78 TeV (C.L. ≲0.8), utilized for the study of the effect of the duty cycle of the TeV data. The five highest energy data points are scaled down by 10^−1/2 of their flux values to that reported in Aharonian et al. (2007b).
Download figure:
Standard image High-resolution image

4.4. Analysis of 1ES 0347−121, 1ES 1101−232, H 2356−309, and RGB J0152+017

Previously published IGMF studies (Neronov & Vovk 2010; Essey et al. 2011a; Tavecchio et al. 2011) have used additional extreme TeV blazars to derive constraints on the IGMF, four of which are considered in this section, 1ES 0347−121, 1ES 1101−232, H 2356−309, and RGB J0152+017. All of these sources were observed by the HESS collaboration prior to the beginning of the Fermi mission. Therefore, to investigate H0, it is necessary to assume that the VHE activity of these sources as characterized by the HESS collaboration typically during the period of 2004–2007 is representative of the VHE activity over the duration of the Fermi mission. The parameters of these data sets are summarized in Table 1. As before, the Fermi-LAT time-averaged dFED for all sources was derived from the beginning of the mission until 2012 February 14 and depending on the strength of the source, was computed for either four or six evenly spaced logarithmic bins over the energy range of 200 MeV–200 GeV.

The VHE data set of 1ES 0347−121 consists of a set of observations over the period of 2006 August–December, during which time 25.4 hr exposure was accumulated. A time-averaged dFED for seven bins over the energy range from 250 GeV–3.67 TeV was derived (Aharonian et al. 2007a). The flux found in the first and fourth (last) bins is weak (TS < 9) for all spectral indices tested, allowing only an upper limit to be established. The flux in the second and third bins is typically found with 9 < TS <25 for the simulated models where the secondary flux dominates the total flux. Figure 13(a) shows the confidence level of simulated models obtained for H0. The best fit models in the α– epsilon _C plane are found at epsilon _C values near 1 TeV, and the dFED for one of these models is illustrated in Figure 13(b). The relatively high confidence level of the 1ES 0347−121 simulated models is partially due to the poor fit of the highest energy bins of the VHE regime where the reported dFED tentatively exhibits a feature of increasing energy density. This trend in the dFED is not accounted for in the set of simulated models investigated. A similar spectral feature appears to be even more pronounced in the VHE data set of 1ES 1101−232, which perhaps may signal unmodeled physics process(es) or a systematic error in the data analyses.

**Figure 13.** Source 1ES 0347−121 analyzed under the assumption of H0. (a) Confidence level, or the probability of exclusion of the gamma-ray source model (with fixed α and _c), where the remaining three parameters of the model (γ, _B, and F₀, see Equation (1)) are chosen so as to minimize χ². (b) Simulated dFED of the best-fit model of α = 1.3, _c = 0.75 TeV, showing both the prompt and secondary cascade contributions to the total dFED, along with the HE and VHE data.
Download figure:
Standard image High-resolution image

The VHE data set of 1ES 1101−232 consists of three periods of observations spanning from 2004 April–2005 March, for a total of 43 hr. The time-averaged dFED is reported for 10 bins over the energy range 225 GeV–4 TeV. The Fermi-LAT dFED for 1ES 1101−232 is especially weak, and in fact, the first bin allows a determination of only an upper limit (TS < 9) for all simulated models tested. Figures 14(a) and (b) illustrate the compatibility of the 1ES 1101−232 VHE and HE data sets with H0, despite the fact that the previously described feature in the highest energy bins of the VHE regime is not well fitted by the models.

**Figure 14.** Source 1ES 1101−232 analyzed under the assumption of H0. (a) Confidence level, or the probability of exclusion of the gamma-ray source model (with fixed α and _c), where the remaining three parameters of the model (γ, _B, and F₀, see Equation (1)) are chosen so as to minimize χ². (b) Simulated dFED of the best-fit model of α = 1.2, _c = 0.75 TeV, showing both the prompt and secondary cascade contributions to the total dFED, along with the HE and VHE data.
Download figure:
Standard image High-resolution image

Observations of H 2356−309 were obtained over the period of 2004 June–September, for a total exposure of 40 hr. The time-averaged dFED was provided for eight bins over the energy range from 200 GeV–1.23 TeV (Aharonian et al. 2006b). The flux in the first energy bin is weakly detected for most simulated models (9 < TS < 25), but the second and third bins exhibit a strong detection (TS > 25). An upper limit is derived for the fourth bin in the HE dFED due to a weak signal present (TS < 9). As illustrated in Figure 15, this source has a very large set of models compatible with the H0. Most of these models, however, suggest that the flux in the two lowest energy bins of the HE regime is dominated by secondary radiation.

**Figure 15.** Source H 2356−309 analyzed under the assumption of H0. (a) Confidence level, or the probability of exclusion of the gamma-ray source model (with fixed α and _c), where the remaining three parameters of the model (γ, _B, and F₀, see Equation (1)) are chosen so as to minimize χ². (b) Simulated dFED of the best-fit model of α = 1.8, _c = 4.22 TeV, showing both the prompt and secondary cascade contributions to the total dFED, along with the HE and VHE data.
Download figure:
Standard image High-resolution image

RGB J0152+017 was observed over the period of 2007 October 30–November 14 for a total exposure of 14.7 hr. A time averaged dFED was reported for six bins over the energy range of 240 GeV–3.6 TeV (Aharonian et al. 2006b). The flux in all but the first energy bin of the Fermi data is strongly detected (TS > 25) for the majority of simulated models. Figures 16(a) and (b) illustrate that this source perhaps provides the weakest constraints on the IGMF with nearly the entire parameter space of simulated models compatible with the VHE and HE spectral data. In only some of these models, the HE part of the spectrum is dominated by the primary emission.

**Figure 16.** RGB J0152+017. (a) Confidence level, or the probability of exclusion of the gamma-ray source model (with fixed α and _c), where the remaining three parameters of the model (γ, _B, and F₀, see Equation (1)) are chosen so as to minimize χ². (b) Simulated dFED of the best-fit model of α = 1.9, _c = 3.16 TeV, showing both the prompt and secondary cascade contributions to the total dFED, along with the HE and VHE data.
Download figure:
Standard image High-resolution image

5. DISCUSSION

In this paper, we have investigated the HE–VHE energy spectrum of seven extreme TeV blazars for which radiation in the HE band may be dominated by the secondary photons produced through cascading. These sources are characterized by their observed hard spectra in the VHE band accompanied by redshifts of order ∼0.1 suggesting a very large energy output into pair production and subsequent cascading. For these sources, observations in the HE band can therefore limit the flux of secondary photons and establish a lower bound on the IGMF. This strategy has been utilized in several publications (Neronov & Vovk 2010; Tavecchio et al. 2010; Dermer et al. 2011; Essey et al. 2011a; Dolag et al. 2011; Taylor et al. 2011; Huan et al. 2011) to suggest B_IGMF ⩾ 10⁻¹⁷–10⁻¹⁸ G in the local cascading environment. In contrast to these studies, we systematically investigated effects of a wide range of uncertainties using detailed three-dimensional Monte Carlo simulations to conclude that H0 remains compatible with current observations.

Two of these blazars RGB J0710+591 and 1ES 1218+304, had VHE observations conducted during the Fermi mission and the measured VHE fluxes of these sources were assumed to be representative of the average VHE activity over the HE measurement period. For four other sources 1ES 0347−121, 1ES 1101−232, H 2356−309, and RGB J0152+017, the VHE component was measured prior to the Fermi mission and it was assumed to be representative of the average VHE activity during the Fermi observations. It was found that, for all of these sources, with Fermi-LAT Pass 7 data and a more general set of intrinsic source models of gamma-ray emission, H0 cannot be rejected at >95% confidence level (even with a typical EBL model; e.g., Domínguez et al. 2011, and models referenced within).

A single source, the extreme blazar 1ES 0229+200, did allow for a rejection of H0 at more than 99% confidence level, under the standard assumptions as discussed in Section 4.3. However, if the EBL model SED is decreased in the visible and near infrared band to within the uncertainty limits of measurements based on the galaxy counts, the B_IGMF lower limit cannot be maintained. H0 can also be reconciled with the data if the VHE measurements of its flux conducted over the period of 40 hours are not assumed to be representative of the past 3.25 yr of activity of this source. We must also point out that a thorough relative energy calibration of the Fermi-LAT and IACTs has not yet been reported and an error in the energy scale of VHE instruments may have significant impact on the conclusions regarding B_IGMF.

It should also be noted that a claim of the rejection of H0 based on the observational data of a single source should not be considered robust due to possible variations of the magnetic field in the local environment of the source. Figure 17 shows the number density of the secondary electrons from 1ES 0229+200 for a zero magnetic field as a function of distance from the source for electrons with energies >10 GeV, >100 GeV, and >1 TeV. It can be seen that a significant fraction of the secondary gamma-ray flux can be isotropized if cascading begins in the vicinity of several tens of megaparsecs of the source in the magnetic field of filaments which is expected to be orders of magnitude larger than the magnetic field in the voids. Recent work by Sutter et al. (2012) and Pan et al. (2012) on identifying and characterizing the Voids in the Sloan Digital Sky Survey Data release 7 (Abazajian et al. 2009), suggests that the volume occupied by them accounts for around 40%–60% of the total volume of the universe. It is therefore critical for studies of the B_IGMF to know the exact line of sight distribution of voids particularly in the few hundred megaparsec vicinity of the source.

**Figure 17.** Electron density as a function of distance from source 1ES 0229+200, assuming a normal EBL model (EBL model 1) and a spectral index in the VHE band (α) of 1.3 and cutoff energy (ε_c) 3.16 TeV. The black (solid) line traces the density of electrons with energies >10 GeV, the red (dashed) line represents the density of electrons with energies >100 GeV, and the blue (dot–dashed) line corresponds to energies >1 TeV.
Download figure:
Standard image High-resolution image

Another avenue to reduce the secondary gamma-ray flux and reconcile it with observations and H0 has recently been proposed by Broderick et al. (2012). The authors argue that beam-plasma instabilities could develop through interaction with the electrons and positrons in the cascade and electrons of the plasma in the voids. If such a collective interaction between two populations of electrons indeed exists, the energy dissipation rate of electron positron pairs into modes of the plasma waves in voids may become significantly larger than the energy loss through IC scattering.

Parameters of the plasma in the voids are largely unknown. To estimate an upper limit on the plasma electron number density in the voids, we assume that the mass of the baryonic matter in the voids does not exceed the difference between the baryonic mass identified through analysis of CMB fluctuations and the mass found in the galaxies and galaxy clusters. The cosmological density of the universe is 5.9 baryons m⁻³ of which 4.6% represents baryonic matter. Ninety percent of the 2.7 × 10⁻⁷ baryons cm⁻³ were identified in filaments and clusters of galaxies filling about half the volume of the universe. Therefore, the amount of baryonic matter in the voids cannot exceed the remaining ≈ 0.5 ×10%, making the density estimate of the electron plasma in the voids n_e ≈ 1.4 × 10⁻⁸ cm⁻³ or lower. The plasma frequency of these electrons, ω_{p, e} = 2πf_{p, e} = (4πn_ee²/m_e)^1/2, is given by f_{p, e} = 1Hz(n_e/10⁻⁸)^1/2.

To estimate the temperature of the electron plasma, we assume that the ionization of hydrogen in the voids occurs sometime during the end of the reionization epoch with redshift z between 6 and 10. At this point, the universe became mostly transparent to UV radiation, and UV photons from ionizing sources in galaxies could propagate to the voids. The photoionization cross section for hydrogen in the ground state by photons with energy ε > ε₀ = 13.6 eV, can be roughly approximated by σ_i(ε/ε₀)⁻³ where σ_i = 6.3 × 10⁻¹⁸ cm². The ionization is suppressed at higher photon energies, suggesting that the characteristic kinetic energy of the plasma electrons acquired by an electron in the ionization process is of order ε₀. These electrons undergo cosmological expansion which reduces the kinetic energy by a factor of (1+z), and they also transfer energy to CMB photons through the IC interaction. The heating of electrons by UV radiation accumulated and reprocessed in the EBL is negligible since the total energy density in the EBL (from the UV to far-IR) is an order of magnitude less than in the CMB. The average IC energy loss of non-relativistic electrons is given by

$\begin{equation} \frac{1}{E}\frac{dE}{dt} = \frac{8}{3}\sigma _{\mathrm{T}} c \frac{U_{{\rm cmb}}}{m_e c^2}, \end{equation} \tag{ 4 }$

where σ_T is the Thomson cross section and U_cmb is the CMB energy density which at the present epoch is U₀ = 0.26 eV cm⁻³. In the context of cosmological expansion, the rate of IC energy loss of non-relativistic particles is given by

$\begin{equation} \ln \left(\frac{E(z)}{E(0)} \right) = \frac{8}{3} \sigma _{\mathrm{T}} \frac{c}{H_0}\frac{U_0}{m_ec^2} \int _{0}^{z} \frac{(1+z^{\prime })^3 dz^{\prime }}{ \sqrt{\Omega _m(1 + z^{\prime })^3 + \Omega _{\Lambda } } }, \end{equation} \tag{ 5 }$

where H₀ is the Hubble constant. For z = 6, the right-hand side of this equation is equal to 1.1 and for z = 10, it is equal to 3.4. This, combined with the reduction of the kinetic energy due to cosmological expansion, suggests that the average kinetic energy of the plasma electrons is reduced by a factor of 5 × 10⁻² (z = 6) or 3 × 10⁻³ (z = 10), from ε₀, implying that the present day temperature of the electron plasma in the voids is a few thousand K.

Given the temperature and density of the electron plasma in the voids, the screening Debye radius is λ_D = (kT/8πn_ee²)^1/2 = $1.6 \;\times\; 10^6 \; {\rm cm} \, \sqrt{\vphantom{A^A}\smash{{{ (T/10^3 \; {\rm K})\, (10^{-8} \; {\rm cm}^{-3}/n_e) }}}}$ . We note that this screening radius is less than the typical distance between electrons in the beam, as estimated in our simulations shown in Figure 17. The bulk of the electrons in the beam is born near the threshold of pair production when a few hundred GeV primary photon interacts with the EBL. The typical Lorentz factor of the outgoing electron is ∼2 × 10⁵. These electrons dissipate ∼90% of their energy through IC scattering over the distance of ∼30 Mpc. For distances from the source larger than this, the typical separation of beam particles is $n_b^{-1/3}$ ≈ 10⁶–10⁷ cm. It appears that electron positron pairs of the beam are significantly screened by the electrons of the plasma in the voids, and therefore, collective interactions to generate a higher rate of energy transfer to plasma oscillations should be strongly suppressed by screening, which was not accounted for in Broderick et al. (2012). The plasma condition at these relatively small distances from the source (∼few megaparsecs where the density of the beam may be higher) is likely to be different from those assumed here for the voids, and the Debye radius estimated may not be applicable in this region.

One may argue that this screening effect can be alleviated if the density of electrons in the voids is considerably lower than the given upper bound estimate. Based on the results of simulations of density fluctuations quoted in Miniati & Elyiv (2013), the present-day density of the electrons in objects with spatial scales of several tens of megaparsecs might be as low as ≈4 × 10⁻¹⁰ cm⁻³. This would make the Debye radius, λ_D ≈ 10⁷ cm, which suggests that screening effects remain critically important for the rate of the energy losses into the excitation of plasma waves. To estimate another critical parameter in the problem, namely, the effect of the angular distribution of electrons in the beam, we consider a hydrodynamical approximation of the beam-plasma instability for the two populations of electrons with zero temperatures and one beam at relative velocity c. The dispersion relation for these waves is given by

$\begin{equation} 1 = \frac{\omega ^2_{p,e}}{\omega ^2} + \frac{\omega ^2_{p,b}}{ ({\bf k} \cdot {\bf c} - \omega)^2 } = \frac{1}{(x+iy)^2} + \frac{\kappa ^2}{(z-x-iy)^2}, \end{equation} \tag{ 6 }$

where ω_{p, b} is the plasma frequency of the electrons and positrons in the beam, and where we introduced the dimensionless variables x = Re(ω)/ω_{p, e}, y = Im(ω)/ω_{p, e}, z = k · c/ω_{p, e} and $\kappa = \omega _{p,b}/\omega _{p,e} = \sqrt{n_b/n_e} \ll 1$ . The real and imaginary part of this equation can be separated, and the unstable (y ≠ 0) solution should satisfy

$\begin{equation} z^2 = \frac{(1 + \varepsilon ^2)}{(1 - \varepsilon ^2)} \frac{(1 - \kappa \varepsilon) (\varepsilon - \kappa)}{\varepsilon }, \end{equation} \tag{ 7 }$

$\begin{equation} \frac{y^2}{x^2} = \varepsilon \frac{\kappa - \varepsilon ^3}{1 - \kappa \varepsilon }, \end{equation} \tag{ 8 }$

where $\varepsilon = \sqrt{z/x - 1}$ . The solution exists for κ < ε < κ^1/3 and for 0 < y/x < 3^1/2(κ/4)^2/3. Since κ ≪ 1, this would require that |z|/|x| ∼ 1 and that |z| ⩽ 1. The last condition can only be satisfied for plasma waves with wave vector precisely perpendicular to the beam velocity. The cosine of this angle must satisfy the following condition cos θ ⩽(λ_D/cω_{p, e})(kλ_D)⁻¹ = (k_BT/2m_ec²)^1/2(kλ_D)⁻¹ = 2.9 × 10⁻⁴(T/10³ K)^1/2(kλ_D)⁻¹, in which kλ_D ≫ 1 to enable collective interaction of plasma particles without significant screening. For the lowest estimate for n_e in the voids, we assume that kλ_D ∼ few, and therefore cos θ may need to be less than 10⁻⁴, but not significantly less than this.

According to our simulations, the angular distribution of electrons in the beam is shown in Figure 18. The bulk of the electrons is created with characteristic angles of order 1/γ ≈ 5 × 10⁻⁶. Although these electrons lose their energy through IC scattering, they retain the angular distribution of the parent particle, due to a small energy transfer to the CMB photons. Therefore, the angular distribution is relatively independent from the distance to the source, as illustrated in Figure 18. Under the condition of the lowest electron density in the voids, maintaining perpendicularity of the wave vector of the plasma waves with respect to the beam, appears to be possible, and the approximation of the beam as completely collimated seems applicable here. All particles of the beam in this reactive regime of interaction will be involved into the energy transfer to the plasma waves making energy losses potentially higher. If kλ_D becomes larger than about a few tens, then plasma waves will interact collectively only with a fraction of beam electrons in the kinetic regime of interaction, which has been pointed out in Miniati & Elyiv (2013). The energy losses in this regime are expected to be lower. Based on the parameters of the electron beam which we derived based on simulations, and reasonable assumptions about properties of electron plasma in the voids, we conclude that all effects of Debye screening and kinetic versus reactive descriptions of two beam instabilities are important. Perhaps the most detailed to date study of the relaxation of beam plasma instabilities in cosmic voids has been reported in Miniati & Elyiv (2013) with the conclusion that the relativistic pair beams of blazars remain stable on timescales much longer than the characteristic IC cooling time of electrons, and collective plasma-beam interaction effects in the voids are negligible. We find that the screening effects not accounted for in this work will only further validate their conclusions. However, the available parameter space in the regime of very low plasma density in the voids may enable the reactive interaction regime and therefore increase the energy loss rate.

An alternative approach to reconcile a lack of the secondary HE radiation in the observational data of extreme blazars is to suggest that a significant part of the VHE emission originates from proton-initiated cascades, rather than from direct photons, if CR protons were also accelerated by the same source (Essey et al. 2011b). The universe is nearly transparent to 10¹⁶–10¹⁹ eV protons on spatial scales of hundreds of megaparsecs and therefore, electromagnetic cascades can be initiated by these particles in a random location along the line of sight. Crossing small regions of intense magnetic fields such as those present in clusters and filaments represent difficulties for this mechanism since they rapidly destroy the correlations between the cosmic ray directions and the source. Moreover, it is likely a formidable task to devise a mechanism by which cosmic rays have been accelerated to ultra HEs in the source by intense magnetic fields satisfying the Hillas condition and then highly collimated along a particular direction, eventually entering regions of potentially extremely small magnetic fields of voids without being significantly disrupted by the intermediate magnetic fields. If this mechanism takes place in nature, it has two distinct observational characteristics, namely that the VHE radiation produced should show little evidence for variability and correlations with other wavelength bands such as X-ray, etc. and most importantly, should be accompanied by higher energy gamma-rays, with energies above 10s of TeV in the VHE spectra of extreme blazars with z ⩾ 0.1. No such observational evidence has been collected so far to necessitate considerations of this scenario, which also lacks an explanation of the origin of the collimation. This mechanism has been recently studied (Murase et al. 2012) in the context of 1ES 0229+200 with the conclusion that the detection of larger than 25 TeV photons would provide an indication of acceleration of ultra high energy cosmic rays in this source.

In summary, we have considered the observational data of seven extreme TeV blazars and performed detailed simulations of cascading in the IGMF of the cosmic voids. We find for all sources, except for 1ES 0229+200, no evidence to claim exclusion of H0 (B_IGMF = 0 hypothesis), due to existing multiple systematic uncertainties in source models. Furthermore, for the case of 1ES 0229+200, we find no definitive evidence that its data is in contradiction with the zero IGMF hypothesis, when astrophysical (EBL, local magnetic field environment) and systematic uncertainties (VHE duty cycle) are accounted for.

We thank the anonymous referee for their critical remarks and important suggestions, which improved the paper significantly. The simulations were performed on the UCLA hoffman2 cluster as well as the Joint Fermilab-KICP Supercomputing Cluster. V.V.V. gratefully acknowledges grants from UCLA and S.P.W. gratefully acknowledges grants from Fermilab, Kavli Institute for Cosmological Physics, and the University of Chicago. This research is supported by the U.S. National Science Foundation under grant Nos. (PHY-0969948, PHY-0422093, and PHY-0969529).

INTERGALACTIC MAGNETIC FIELDS AND GAMMA-RAY OBSERVATIONS OF EXTREME TeV BLAZARS

Article metrics

Permissions

Author e-mails

Author affiliations

ORCID iDs

Dates

ABSTRACT

1. INTRODUCTION