# Intergalactic Magnetic Fields and Gamma Ray Observations of Extreme TeV Blazars

###### Abstract

The intergalactic magnetic field (IGMF) in cosmic voids can be indirectly probed through its effect on electromagnetic cascades initiated by a source of TeV gamma-rays, such as active galactic nuclei (AGN). AGN that are sufficiently luminous at TeV energies, “extreme TeV blazars” can produce detectable levels of secondary radiation from inverse Compton (IC) scattering of the electrons in the cascade, provided that the IGMF is not too large. We review recent work in the literature which utilizes this idea to derive constraints on the IGMF for three TeV-detected blazars-1ES 0229+200, 1ES 1218+304, and RGB J0710+591, and we also investigate four other hard-spectrum TeV blazars in the same framework. Through a recently developed detailed 3D particle-tracking Monte Carlo code, incorporating all major effects of QED and cosmological expansion, we research effects of major uncertainties such as the spectral properties of the source, uncertainty in the UV - far IR extragalactic background light (EBL), undersampled Very High Energy (VHE; energy GeV) coverage, past history of gamma-ray emission, source vs. observer geometry, and jet AGN Doppler factor. The implications of these effects on the recently reported lower limits of the IGMF are thoroughly examined to conclude that presently available data are compatible with a zero IGMF hypothesis.

###### Subject headings:

cosmology: large-scale structure of universe, cosmology: cosmic background radiation, gamma rays: general, magnetic fields, methods: numerical, radiation mechanisms: non-thermal## 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 micro Gauss (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 nano Gauss 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 Gauss 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 which may become a more sensitive tool for the measurement of IGMF characteristics. This technique relies on observations of blazars, believed to be AGN whose jet is oriented along the line of sight to Earth, in the gamma-ray energy range from 100 MeV to greater than 10 TeV and is described by several authors (Neronov & Semikoz (2006); Eungwanichayapant & Aharonian (2009); Dolag et al. (2009); Elyiv et al. (2009)). Briefly, TeV-scale gamma-rays from the blazar interact with the EBL, producing an electron-positron pair. The electrons and positrons then undergo 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 GeV scale carry information on the strength of the magnetic field in which the cascading occurs. If cascading develops in the voids, then spectral, temporal, and angular characteristics of the secondary radiation will depend on 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 GeV-band emission does not exceed current measurements. In another study, Tavecchio et al. (2011) performed detailed modeling of the spectral energy distribution 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 semi-analytic 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 semi-analytic 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 which appeared to be consistent for 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 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 which 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, while section 4 provides 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 progagating through a Milky Way-like galaxy will have an optical depth of , based on rough estimates of the energy density of the Galaxy in the far infrared ( ). 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 ( Mpc). These photons will initiate cascades under the influence of galactic magnetic fields (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 G (see, e.g. Widrow (2002)) or intervening large scale structure filaments, with magnetic fields G, then the secondary electrons will be isotropized, thereby dramatically attenuating the observable flux of the secondary photons produced by them. Effects of the electromagentic cascading may become observable with present day instrumentation when the interaction with the EBL occurs in the voids with characteristic magnetic fields 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 3-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 , , and respectively. The radiation density is negligible for nearby sources of z . The standard EBL model chosen in the simulations at z = 0 is modeled with 6 energy density points (m) = nW m sr, , , , , ), 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.

In this study, we investigate sources with redshift 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 .

### 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, . For the coherence length of magnetic fields, , the IC scattering effectively happens in the environment of nearly constant B. For , 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 Mpc (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.0. For each of the energy bins (4 per decade), the mean delay time is computed for 6 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 hours. The saturation effect at low energies is due to the aperture cutoff. Figure 2a 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 , derived analytically by making several simplifications, as explained in Neronov & Semikoz (2009). Figure 2b 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 towards 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.

### 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 , it can be parameterized as

where is the photon energy in the reference frame of the host galaxy, = (1 - )], is the viewing angle from the blazar jet axis to the line of sight of the observer, is a flux normalization factor, and is an energy scale factor. To account for absorption of photons at the highest energies inside the host galaxy, an exponential cutoff at energy is introduced. This model is characterized by 5 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 5 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

(1) |

The seventh parameter, , 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, = 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 4a shows the spectrum of secondary photons when the source is viewed at = 0, 2, 5, 10 and Figure 4b shows the spectrum of secondary photons for = 5, 10, 30, 100. Viewing a source with the same prompt spectral energy density (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 (), 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, = 30 TeV, B = 10 gauss, = 30, at four different observing angles, = 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 the photon distribution around the source becomes increasingly axially non-symmetric, when . 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 intiated by VHE photons is modeled by full 3-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 due to IC scattering is given by

(2) |

where is the speed of the electron, is the differential IC cross section, where is the collision angle and 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 arbritrary precision numerical integrator developed at Lawrence
Berkeley National
Laboratory^{2}^{2}2http://crd-legacy.lbl.gov/$∼$dhbailey/mpdist/. 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 3-dimensional simulation.

Figure 6 illustrates the importance of higher order cascading to correctly describe secondary photons with energies GeV. This figure shows the results of the simulations of secondary emission produced for a source at z=, with B G, , and . The two panels display different intrinsic spectra for the source; panel a) is obtained with and TeV, and panel b) corresponds to and TeV. Both panels show the direct differential flux energy density (dFED) of the source together with the time-integrated secondary dFED, within 0.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 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.

## 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-LAT instrument and the data were processed utilizing publically available tools, version v9r23p1, with the update from November 6, 2011. 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 hypothesis (hereafter called ) 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 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

(3) |

is used as suggested in Aharonian et al. (2006c). The typical values of parameters of the PSF are = 0.06505, = 0.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 August 4, 2008 or MJD 54682) is considered in this study. As reported in Acciari et al. (2010a), 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., 2007c). VERITAS has continued monitoring this source since the Fermi launch and tentatively detected flux variations on a yearly time scale (private communication, J.S. Perkins and VERITAS collaboration). 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.

IACT data summary | ||||||
---|---|---|---|---|---|---|

Source | z | IACT | Flux [10cm s] | index | MJD (approx) | Hrs |

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 |

### 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
10 G. For 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 2 year time-averaged sample of
AGN^{3}^{3}3http://fermi.gsfc.nasa.gov/ssc/data/analysis/documentation/

Pass7_usage.html.

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^{4}^{4}4http://fermi.gsfc.nasa.gov/ssc/data/analysis/
of the flux in each energy bin in the Fermi data, within the
10 region of interest (ROI) which also includes all nearby
sources from the 2 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, , that 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 8 bins per decade, of equal width in logarithmic space. The simulated spectral energy density data are equally binned with 8 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 3 decades of the HE regime, 4-6 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-, , , and . For each model, is chosen in the range (1.0, 2.5) with a step of 0.1, and is chosen in the range (600 GeV, 60 TeV) with 8 bins per decade, equally spaced in logarithmic energy. Each model is characterized with default values of = 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 and the spectral index are relevant only to the HE part of the spectrum. They were chosen by minimizing the value, by allowing 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 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 and which may or may not be relevant to the , depending on the relationship between the prompt and secondary emission of the source. We chose to determine by minimization of the VHE part of by solving = 0 independently from the behavior of . 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 = + , 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 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 or measurements.

The value for each model with four fixed parameters (, , , and ) was converted into a confidence level, using statistics with degrees of freedom, assuming that three parameters were optimized for each model (, , ). 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 .

As an illustration of a typical result of the data and model comparison, Figure 7a shows the confidence level of the dFED, assuming four fixed and three free parameters for each model tested within the given range of and parameters. The most favored model with value of = 1.8 and = 3.16 TeV is incompatible with the data at the level of about 75 %. If this choice of parameters and 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 7b shows the data points for the dFED and the best fit simulation result with these , 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 4 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 August 4, 2008 to February 14, 2012, 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 ( 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 Fig. 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 4 years) could be viewed as “contemporaneous” to IACT measurements, for the purposes of verification of . 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 5 energy data points reported by the VERITAS collaboration in Acciari et al. (2010b), observed during the time period December 2008 - March 2009 for a total of 22 hours. For the HE regime, 6 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 .

Figure 7a shows the confidence level for rejecting , for a set of models characterized by the range of and described in section 3.3. It identifies the best fit model with values of = 1.8 and = 3.16 TeV, which is incompatible with 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 . The best fit of the simulated dFED and observations is shown in Figure 7b. It appears that the conclusion of TVN11 that is ruled out at the 98.8 % level is invalidated due to two factors. First, the Fermi-LAT dataset 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 2 data sets. The first set, obtained during December 2008 - May 2009, is based on 27 hours of data and has 9 energy data points reported by the VERITAS collaboration in Acciari et al. (2010a). 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 6 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 8a shows the confidence level for rejecting the hypothesis on the - plane. It suggests a best fit model with values of = 1.8, = 3.16 TeV, a spectral break energy of = 10 GeV, and an index below the break energy of = 1.0, which is incompatible with 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 8b. It is evident that the conclusion of TVN11, that 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 .

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 3-4 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 2 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 souce could be lower than used in TVN11, yet the assumption of the higher average VHE flux is still compatible with .

### 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 September 1, 2005 to December 19, 2006, accumulating 41.8 hours exposure (Aharonian et al. (2007c)). 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 ), for all simulated models in which the secondary flux dominates the total flux. The TS for all other energy bins is typically found at for the majority of simulated models but in some cases, only the upper limit can be established (TS ). 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 years of the Fermi-LAT data used in this work. The fit obtained under this assumption and for is shown in Figure 9a. 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 is ruled out at the 99.5 % confidence level. The best fit (, 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 9b.

The effects of two systematic uncertainties, viewing angle and Doppler factor , were considered. The default assumptions in the analysis are and . Increasing the viewing angle, e.g. as used in TVN11, would further overproduce radiation in the HE regime (see Figure 4a). Increasing the Doppler factor, , combined with the 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 , 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 - , and it cannot be used to reconcile the observational data of 1ES 0229+200 with .

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 & 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 . 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 high energy cutoffs, , below TeV, and are thus insensitive to EBL photon wavelengths 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 HST (Madau & Pozzetti (2000)) in the visible (from 0.36 - 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 . 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 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 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 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 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 error bar from these later results. Thus, EBL model 3 is compatible with the galaxy counts results to within 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 are shown in Figure 10a. 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 , TeV is shown in Figure 10b.

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 , 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 .

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 and would require a lower bound of = G (Vovk et al. (2012)). All source models analyzed in that paper had a single cutoff energy = 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 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, can be reconciled with observations of this source.

Another avenue to make the HE-VHE observational data of 1ES 0229+200 compatible with 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 hrs) - 2006 (35 hrs) 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 cm s, and 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 - 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 5 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 at more than 95% confidence level. One of these compatible models, with = 1.3 = 1.78 is shown in Figure 12. Therefore, the conclusion that the is ruled out with high significance heavily rests on the assumption that the HESS measurements are representative of the average flux for E TeV and a 10 change of the highest energy part of the spectrum invalidates this conclusion.

### 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 , 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 February 14, 2012, and depending on the strength of the source, was computed for either 4 or 6 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 August - December 2006, during which time 25.4 hours exposure was acculumlated. A time-averaged dFED for 7 bins over the energy range from 250 GeV - 3.67 TeV was derived (Aharonian et al. (2007b)). 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 TS for the simulated models where the secondary flux dominates the total flux. Figure 13a shows the confidence level of simulated models obtained for . The best fit models in the - plane are found at values near 1 TeV, and the dFED for one of these models is illustrated in Figure 13b. 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.

The VHE data set of 1ES 1101-232 consists of 3 periods of observations spanning from April 2004 - March 2005, for a total of 43 hours. 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 14a and b illustrate the compatibility of the 1ES 1101-232 VHE and HE data sets with , despite the fact that the previously described feature in the highest energy bins of the VHE regime is not well fitted by the models.

Observations of H 2356-309 were obtained over the period of June - September 2004, for a total exposure of 40 hours. 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 . Most of these models, however, suggest that the flux in the two lowest energy bins of the HE regime is dominated by secondary radiation.

RGB J0152+017 was observed over the period of October 30 - November 14 2007 for a total exposure of 14.7 hours. 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. Figure 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.

## 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 G in the local cascading environment. In contrast to these studies, we systematically investigated effects of a wide range of uncertainties using detailed 3D Monte Carlo simulations to conclude that 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, cannot be rejected at 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 at more than 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 lower limit cannot be maintained. can also be reconciled with the data if the VHE measurements of its flux conducted over the period of forty hours are not assumed to be representative of the past 3.25 years 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 .

It should also be noted that a claim of the rejection of 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 GeV, GeV, and 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 Mpc 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 to know the exact line of sight distribution of voids particularly in the few hundred Mpc vicinity of the source.

Another avenue to reduce the secondary gamma-ray flux and reconcile it with observations and 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 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 10 %, making the density estimate of the electron plasma in the voids cm or lower. The plasma frequency of these electrons, , is given by .

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 – 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 eV, can be roughly approximated by where = 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 inverse Compton (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

(4) |

where is the Thomson cross section and is the CMB energy density which at the present epoch is eV/cm. In the context of cosmological expansion, the rate of IC energy loss of non-relativistic particles is given by

(5) |

where is the Hubble constant. For z = 6, the RHS 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 (z = 6) or (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 = . 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 electon is . These electrons dissipate % of their energy through IC scattering over the distance of Mpc. For distances from the source larger than this, the typical separation of beam particles is 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 Mpc 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 (2012), the present-day density of the electrons in objects with spatial scales of several tens of Mpc might be as low as cm. This would make the Debye radius, 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 paramter 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

(6) |

where is the plasma frequency of the electrons and positrons in the beam, and where we introduced the dimensionless variables , , and . The real and imaginary part of this equation can be separated, and the unstable (y 0) solution should satisfy

(7) |

(8) |

where . The solution exists for and for . Since , this would require that and that . 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 = , in which to enable collective interaction of plasma particles without significant screening. For the lowest estimate for in the voids, we assume that few, and therefore may need to be less than , but not significatnly less than this.

According to our simulations, the angular distribution of electrons in the beam is shown in Fig. 18. The bulk of the electrons is created with characteristic angles of order . 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 Fig. 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 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 (2012). 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 vs. 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 (2012) 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 - eV protons on spatial scales of hundreds of Mpc 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 high energies 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 7 extreme TeV blazars and performed detailed simulations of cascading in the intergalactic magnetic field of the cosmic voids. We find for all sources, except for 1ES 0229+200, no evidence to claim exclusion of ( 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 grants no. (PHY-0969948, PHY-0422093, and PHY-0969529).

## References

- Abazajian et al. (2009) Abazajian, K. N., et al. 2009, ApJS, 182, 543
- Acciari et al. (2009) Acciari, V. A., et al. 2009, ApJ, 695, 1370
- Acciari et al. (2010a) —. 2010a, ApJ, 709, L163
- Acciari et al. (2010b) —. 2010b, ApJ, 715, L49
- Aharonian et al. (2006a) Aharonian, F., et al. 2006a, Nature, 440, 1018
- Aharonian et al. (2006b) —. 2006b, A&A, 455, 461
- Aharonian et al. (2006c) —. 2006c, A&A, 457, 899
- Aharonian et al. (2007a) —. 2007a, A&A, 470, 475
- Aharonian et al. (2007b) —. 2007b, A&A, 473, L25
- Aharonian et al. (2007c) —. 2007c, A&A, 475, L9
- Aharonian et al. (2008) —. 2008, A&A, 481, L103
- Aharonian et al. (1994) Aharonian, F. A., Coppi, P. S., & Voelk, H. J. 1994, ApJ, 423, L5
- Ahlers (2011) Ahlers, M. 2011, Phys. Rev. D, 84, 063006
- Albert et al. (2006) Albert, J., et al. 2006, ApJ, 642, L119
- Berta et al. (2010) Berta, S., et al. 2010, A&A, 518, L30
- Béthermin et al. (2010) Béthermin, M., Dole, H., Beelen, A., & Aussel, H. 2010, A&A, 512, A78
- Biermann (1950) Biermann, L. 1950, Zeitschrift Naturforschung Teil A, 5, 65
- Blasi et al. (1999) Blasi, P., Burles, S., & Olinto, A. V. 1999, ApJ, 514, L79
- Broderick et al. (2012) Broderick, A. E., Chang, P., & Pfrommer, C. 2012, ApJ, 752, 22
- Carilli & Taylor (2002) Carilli, C. L., & Taylor, G. B. 2002, ARA&A, 40, 319
- Dermer et al. (2011) Dermer, C. D., Cavadini, M., Razzaque, S., Finke, J. D., Chiang, J., & Lott, B. 2011, ApJ, 733, L21
- Dolag et al. (2009) Dolag, K., Kachelrieß, M., Ostapchenko, S., & Tomàs, R. 2009, ApJ, 703, 1078
- Dolag et al. (2011) Dolag, K., Kachelriess, M., Ostapchenko, S., & Tomàs, R. 2011, ApJ, 727, L4+
- Dole et al. (2006) Dole, H., et al. 2006, A&A, 451, 417
- Domínguez et al. (2011) Domínguez, A., et al. 2011, MNRAS, 410, 2556
- Elbaz et al. (2002) Elbaz, D., Cesarsky, C. J., Chanial, P., Aussel, H., Franceschini, A., Fadda, D., & Chary, R. R. 2002, A&A, 384, 848
- Elyiv et al. (2009) Elyiv, A., Neronov, A., & Semikoz, D. V. 2009, Phys. Rev. D, 80, 023010
- Essey et al. (2011a) Essey, W., Ando, S., & Kusenko, A. 2011a, Astroparticle Physics, 35, 135
- Essey et al. (2011b) Essey, W., Kalashev, O., Kusenko, A., & Beacom, J. F. 2011b, ApJ, 731, 51
- Eungwanichayapant & Aharonian (2009) Eungwanichayapant, A., & Aharonian, F. 2009, International Journal of Modern Physics D, 18, 911
- Fazio et al. (2004) Fazio, G. G., et al. 2004, ApJS, 154, 10
- Franceschini et al. (2008) Franceschini, A., Rodighiero, G., & Vaccari, M. 2008, A&A, 487, 837
- Gould & Schréder (1967) Gould, R. J., & Schréder, G. P. 1967, Phys. Rev., 155, 1404
- Grasso & Rubinstein (2001) Grasso, D., & Rubinstein, H. R. 2001, Phys. Rep., 348, 163
- Huan et al. (2011) Huan, H., Weisgarber, T., Arlen, T., & Wakely, S. P. 2011, ApJ, 735, L28+
- Ichiki et al. (2008) Ichiki, K., Inoue, S., & Takahashi, K. 2008, ApJ, 682, 127
- Jelley (1966) Jelley, J. V. 1966, Phys. Rev. Lett., 16, 479
- Kneiske & Dole (2010) Kneiske, T. M., & Dole, H. 2010, A&A, 515, A19
- Kronberg (1994) Kronberg, P. P. 1994, Nature, 370, 179
- Kronberg et al. (2001) Kronberg, P. P., Dufton, Q. W., Li, H., & Colgate, S. A. 2001, ApJ, 560, 178
- Kronberg & Perry (1982) Kronberg, P. P., & Perry, J. J. 1982, ApJ, 263, 518
- Kulsrud & Zweibel (2008) Kulsrud, R. M., & Zweibel, E. G. 2008, Reports on Progress in Physics, 71, 046901
- Madau & Pozzetti (2000) Madau, P., & Pozzetti, L. 2000, MNRAS, 312, L9
- Matsuura et al. (2011) Matsuura, S., et al. 2011, ApJ, 737, 2
- Mazin & Raue (2007) Mazin, D., & Raue, M. 2007, A&A, 471, 439
- Miniati & Elyiv (2012) Miniati, F., & Elyiv, A. 2012, ArXiv e-prints
- Murase et al. (2012) Murase, K., Dermer, C. D., Takami, H., & Migliori, G. 2012, ApJ, 749, 63
- Murase et al. (2008) Murase, K., Takahashi, K., Inoue, S., Ichiki, K., & Nagataki, S. 2008, ApJ, 686, L67
- Neronov et al. (2010) Neronov, A., Semikoz, D., Kachelriess, M., Ostapchenko, S., & Elyiv, A. 2010, ApJ, 719, L130
- Neronov & Semikoz (2006) Neronov, A., & Semikoz, D. V. 2006, ArXiv Astrophysics e-prints
- Neronov & Semikoz (2009) —. 2009, Phys. Rev. D, 80, 123012
- Neronov & Vovk (2010) Neronov, A., & Vovk, I. 2010, Science, 328, 73
- Pan et al. (2012) Pan, D. C., Vogeley, M. S., Hoyle, F., Choi, Y.-Y., & Park, C. 2012, MNRAS, 421, 926
- Papovich et al. (2004) Papovich, C., et al. 2004, ApJS, 154, 70
- Perkins & VERITAS Collaboration (2010) Perkins, J. S., & VERITAS Collaboration. 2010, in Bulletin of the American Astronomical Society, Vol. 42, AAS/High Energy Astrophysics Division #11, 708
- Plaga (1995) Plaga, R. 1995, Nature, 374, 430
- Ryu et al. (2008) Ryu, D., Kang, H., Cho, J., & Das, S. 2008, Science, 320, 909
- Sutter et al. (2012) Sutter, P. M., Lavaux, G., Wandelt, B. D., & Weinberg, D. H. 2012, ArXiv e-prints
- Tavecchio et al. (2011) Tavecchio, F., Ghisellini, G., Bonnoli, G., & Foschini, L. 2011, MNRAS, 414, 3566
- Tavecchio et al. (2010) Tavecchio, F., Ghisellini, G., Foschini, L., Bonnoli, G., Ghirlanda, G., & Coppi, P. 2010, MNRAS, 406, L70
- Taylor et al. (2011) Taylor, A. M., Vovk, I., & Neronov, A. 2011, A&A, 529, A144
- Urry & Padovani (1995) Urry, C. M., & Padovani, P. 1995, PASP, 107, 803
- Vassiliev (2000) Vassiliev, V. V. 2000, Astroparticle Physics, 12, 217
- Vovk et al. (2012) Vovk, I., Taylor, A. M., Semikoz, D., & Neronov, A. 2012, ApJ, 747, L14
- Widrow (2002) Widrow, L. M. 2002, Reviews of Modern Physics, 74, 775
- Zweibel (2006) Zweibel, E. G. 2006, Astronomische Nachrichten, 327, 505