Properties of the Intergalactic Magnetic Field Constrained by Gamma-Ray Observations of Gamma-Ray Bursts

The source, in our case a GRB, emits TeV range ($0.1\lesssim E={E}_{\mathrm{TeV}}$ TeV ≲ 100 TeV) photons that interact with EBL photons on the pair-production distance scale ${\lambda }_{\gamma \gamma }(E)$. This can be calculated for a specific EBL model that provides an optical depth $\tau (E,z)$ at photon energy E and redshift z, noting that ${\lambda }_{\gamma \gamma }(E)\cong D/\tau (E,z)$, where D is the distance to the source.

This interaction yields an electron–positron pair, each with Lorentz factor ${\gamma }_{e}={10}^{6}{\gamma }_{6}$ and energy ${E}_{e}={m}_{e}{c}^{2}{\gamma }_{e}$, that are on average half of the energy of the TeV range photon, so that ${\gamma }_{6}\approx {E}_{\mathrm{TeV}}$. These pairs interact with the CMB photons and lose energy through inverse Compton (IC) scattering. The electron has a mean free path ${\lambda }_{e}\ =\ {({n}_{\mathrm{CMB}}{\sigma }_{T})}^{-1}=0.40\,\mathrm{kpc}\,{(1+z)}^{-3}$ between scatterings, where ${n}_{\mathrm{CMB}}\approx 409\,{(1+z)}^{3}\,{\mathrm{cm}}^{-3}$ is the CMB photon number density for temperature ${T}_{\mathrm{CMB}}\,=2.725\,(1+z)\,{\rm{K}}$ (Fixsen 2009). The electrons lose their energy through repeated scatterings on a cooling length scale of

$$\begin{eqnarray}&&{\lambda }_{T}({\gamma }_{e})=\displaystyle \frac{{m}_{e}{c}^{2}}{\tfrac{4}{3}{\sigma }_{T}{u}_{\mathrm{CMB}}{\gamma }_{e}}=0.72\,{(1+z)}^{-4}\,{\gamma }_{6}^{-1}\,\mathrm{Mpc},\end{eqnarray} \tag{ 1 }$$

where ${u}_{\mathrm{CMB}}\approx 4.2\times {10}^{-13}{(1+z)}^{4}\,\mathrm{erg}\,{\mathrm{cm}}^{-3}$ is the energy density of the CMB. Electrons with ${\gamma }_{6}\approx 1$ scatter on the order of ${\lambda }_{T}/{\lambda }_{e}\sim 1000$ times before losing a substantial part of their energy. The upscattered photon's energy is ${E}_{\mathrm{IC}}\approx 4/3\,(2.7\,{k}_{B}\,{T}_{\mathrm{CMB}}){\gamma }_{e}^{2}\cong 0.8\,{\gamma }_{6}^{2}\,\mathrm{GeV}$. The scattering is in the Thomson regime provided $4{\gamma }_{e}(E/{m}_{e}{c}^{2})\ll 1$, which implies that ${\gamma }_{e}\ll 2\times {10}^{8}$. Klein–Nishina effects therefore become important only for $\gtrsim 100$ TeV photons and, though we use the full Klein–Nishina kernel in our calculations, can be neglected in our study since the cutoff energy of the GRB emission is always assumed to be $\lesssim 30$ TeV.

The process is sketched in Figure 1. The angle between the direction of the VHE photon and the upscattered CMB photon is

$$\begin{eqnarray}&&{\theta }_{B}={\lambda }_{T}/{R}_{L}=1.3\times {10}^{-5}\,{(1+z)}^{-4}\,{B}_{-20}\,{\gamma }_{6}^{-2},\end{eqnarray} \tag{ 2 }$$

where ${R}_{L}=\gamma {m}_{e}{c}^{2}/{q}_{e}B=55{\gamma }_{6}\,{B}_{-20}^{-1}\,\mathrm{Gpc}$ is the Larmor radius in a uniform field of strength ${10}^{-20}{B}_{-20}$ G. This equation is valid when the IGMF is coherent on scales larger than the IC-cooling length ${\lambda }_{T}$, given by Equation (1). In the opposite case, the deflection angle is modified by a factor of $\sqrt{{R}_{\mathrm{coh}}/{\lambda }_{T}}$, reflecting the random walk of the electron as it cools.

Zoom In Zoom Out Reset image size

Following the notation of Dermer et al. (2011), the time delay of the echo photons arriving at the detector compared to the arrival time of photons that are emitted in the direction to the observer arrive can be calculated from the differences in the path length, giving

$$\begin{eqnarray}&&c{\rm{\Delta }}t=\lambda +x-D\approx \displaystyle \frac{1}{2}{\lambda }_{\gamma \gamma }{\theta }_{B}^{2}\left(1-\displaystyle \frac{{\lambda }_{\gamma \gamma }}{D}\right),\end{eqnarray} \tag{ 3 }$$

where $\lambda ={\lambda }_{\gamma \gamma }+{\lambda }_{{\rm{T}}}\cong {\lambda }_{\gamma \gamma }$, and we have expressed x through the sine theorem: $x=D\sin {\theta }_{1}/\sin {\theta }_{B}$ (see Figure 1).

Because the time delay can vary over many orders of magnitude depending on the values of $B,{R}_{\mathrm{coh}},E$, and z, it is instructive to give order-of-magnitude estimates for the time delay in specific cases. For $E=0.5\,\mathrm{TeV}$ and z = 0.34 (the redshift of GRB 130427A), we have ${\lambda }_{\gamma \gamma }\approx 700\,\mathrm{Mpc}$ and ${\rm{\Delta }}t\cong 0.3{B}_{20}^{-2}$ year. In the case of an $E=10\,\mathrm{TeV}$ photon, the mean free path for pair production with photons of the EBL is ${\lambda }_{\gamma \gamma }\approx 110\,\mathrm{Mpc}$ and the time delay is ${\rm{\Delta }}t\approx 9{B}_{-20}^{2}$ s.

We generate a distribution of photons in the specified VHE energy range assumed to be described by a power-law distribution with high-energy cutoff. To mimic the EBL absorption and to select the population for pair creation, we retain photons for which a randomly generated number (uniformly distributed between 0 and 1) exceeds $1-{e}^{-\tau (E)}$. The distance at which each photon pair produces is drawn from a random exponential distribution with ${\lambda }_{\gamma \gamma }(E)$ as the average.

The TeV photons produce pairs by interacting with the EBL. For analytical calculations (e.g., Dermer et al. 2011) and some previous numerical treatments (e.g., Takahashi et al. 2011; Fitoussi et al. 2017), it is customary to assume that the pairs equally share the energy of the parent TeV photon. Here, we use the appropriate distribution function for the pairs from Equation (B1) in Zdziarski (1988). The distribution in energies of an electron with energy ${\gamma }_{e}$, given a VHE photon with energy $\epsilon {m}_{e}{c}^{2}$, is given by

$$\begin{eqnarray}P({\gamma }_{e},\epsilon ) & = & {\displaystyle \int }_{\tfrac{\epsilon }{4{\gamma }_{e}{\gamma }_{e}^{{\prime} }}}^{\infty }d{\epsilon }_{\mathrm{EBL}}n({\epsilon }_{\mathrm{EBL}})\displaystyle \frac{3{\sigma }_{{\rm{T}}}c}{4{\epsilon }^{2}\gamma }\\ & & \times \left(r-(2+r)\displaystyle \frac{\epsilon }{4{\epsilon }_{\mathrm{EBL}}{\gamma }_{e}{\gamma }_{e}^{{\prime} }}\,+\,\right.2{\left(\displaystyle \frac{\epsilon }{4{\epsilon }_{\mathrm{EBL}}{\gamma }_{e}{\gamma }_{e}^{{\prime} }}\right)}^{2}\\ & & +\left.\displaystyle \frac{\epsilon }{2{\epsilon }_{\mathrm{EBL}}{\gamma }_{e}{\gamma }_{e}^{{\prime} }}\mathrm{ln}\displaystyle \frac{4{\epsilon }_{\mathrm{EBL}}{\gamma }_{e}{\gamma }_{e}^{{\prime} }}{\epsilon }\right),\end{eqnarray} \tag{ 4 }$$

where ${\gamma }_{e}^{{\prime} }=\epsilon -{\gamma }_{e}$, $r=({\gamma }_{e}/{\gamma }_{e}^{{\prime} }+{\gamma }_{e}^{{\prime} }/{\gamma }_{e})/2$, and ${\epsilon }_{\mathrm{EBL}}$ and ${n}_{\mathrm{EBL}}$ are the EBL photon energy and number density, respectively, provided in Finke et al. (2010). For different VHE photon energies, the distribution is plotted in Figure 2. The average value is indeed at ${\gamma }_{e}=\epsilon /2$. At energies $\gtrsim \,\mathrm{TeV}$, however, the distribution of secondary pairs starts to become increasingly unequal. For each VHE photon with energy ${E}_{\mathrm{TeV}}=\epsilon \,{{\rm{m}}}_{e}\,{{\rm{c}}}^{2}$, we draw from this distribution and retain ${\gamma }_{e}$ and ${\gamma }_{e}^{{\prime} }$. The result of taking the appropriate pair-production cross-section versus the simple approach of assigning half of the TeV photons' energy to pairs will be an increased number of the highest-energy pairs, and consequently the GeV spectrum will reach to higher energies, e.g., if we imagine a monochromatic beam of photons, with 1 TeV energy in the simple scenario the electrons will all have Lorentz factors of ${\gamma }_{e}\approx {10}^{6}$, and the scattered photons will have an energy of ≈0.8 GeV. Using the exact cross-section, we can have a significant fraction of pairs at γ_e ≈ (0.85/0.5)10⁶ (see Figure 2, 0.85 corresponds to the second peak of the red curve showing 1 TeV photons). Thus, the scattered photons with the highest energy will reach ${(0.85/0.5)}^{2}\,\times \,0.8$ GeV = 2.3 GeV. This effect is more pronounced above 1 TeV: a single electron is more likely to carry a high fraction of the VHE photons with increasing energy.

Zoom In Zoom Out Reset image size

The pairs travel an average distance ${\lambda }_{T}$ before losing their energy by scattering CMB photons. In our simulation, we again draw from an exponential distribution with ${\lambda }_{T}({\gamma }_{e})$ as the mean. Only one generation of secondaries is followed. The Compton-scattered radiation may again be susceptible to pair attenuation by the CMB for the highest-energy photons. For the 30 TeV maximum energy of GRB photons assumed in this study, the second generation emission in most of the cases can be neglected. The most energetic echo photons come from the more energetic original photons, which interact with the EBL close to the source, namely, at the same redshift. Where appropriate, we apply EBL absorption to the echo spectrum. We use the model described in Finke et al. (2010) for the EBL, which is similar to currently favored models (Biteau & Williams 2015; Khaire & Srianand 2015; Stecker et al. 2016).

The differential distribution of pairs created from the interaction of TeV photons with the EBL is denoted by ${{dN}}_{0}({\gamma }_{e})/d{\gamma }_{e}$. The echo radiation spectrum is calculated by integrating the electron IC power over the electron distribution using the expression

$$\begin{eqnarray}&&\displaystyle \frac{{d}^{2}{N}_{\mathrm{echo}}}{{{dE}}_{\mathrm{IC}}{dt}}=\int d{\gamma }_{e}\displaystyle \frac{{dN}({\gamma }_{e})}{d{\gamma }_{e}}\displaystyle \frac{{d}^{2}{N}_{{\rm{E}},\mathrm{IC}}}{{dEdt}}.\end{eqnarray} \tag{ 5 }$$

The single-electron power when scattering CMB photons is

$$\begin{eqnarray}&&\displaystyle \frac{{d}^{2}{N}_{{\rm{E}},\mathrm{IC}}}{{dEdt}}=\displaystyle \frac{3{\sigma }_{T}c}{4{\gamma }_{e}^{2}}\int \displaystyle \frac{{{dE}}_{\mathrm{CMB}}}{{E}_{\mathrm{CMB}}}{n}_{{\rm{E}},\mathrm{CMB}}f(x),\end{eqnarray} \tag{ 6 }$$

where $f(x)=2x\mathrm{ln}x+x+1-2{x}^{2}$ and $x=\epsilon {m}_{e}{c}^{2}/4{\gamma }_{e}^{2}{E}_{\mathrm{CMB}}$ (Blumenthal & Gould 1970). The time-integrated IC photon flux can be obtained by multiplying Equation (6) with the local IC-cooling time of the pairs:

$$\begin{eqnarray}&&{\rm{\Delta }}t{{\prime} }_{\mathrm{IC}}({\gamma }_{e})=\displaystyle \frac{{m}_{e}{c}^{2}}{\tfrac{4}{3}{\gamma }_{e}{\sigma }_{T}c\,{u}_{\mathrm{CMB}}}=7.3\times {10}^{13}\,{\gamma }_{6}^{-1}\,{\rm{s}}\end{eqnarray} \tag{ 7 }$$

(Fan et al. 2004).

In the analytic treatment of this problem, we need to link the distribution of the electrons (${{dN}}_{0}/d{\gamma }_{e}$) to the pairs that contribute within the observing window (${dN}/d{\gamma }_{e}$; e.g., Dai & Lu 2002; Dai et al. 2002; Razzaque et al. 2004). To accomplish this, one has to consider the maximum timescale ${\rm{\Delta }}{t}_{\mathrm{obs}}(\gamma )$ of the angular, magnetic-deflection, IC-cooling, and GRB timescales as a function of ${\gamma }_{e}$ (Razzaque et al. 2004), giving ${dN}/d{\gamma }_{e}=({\rm{\Delta }}{t}_{\mathrm{IC}}^{{\prime} }/{\rm{\Delta }}{t}_{\mathrm{obs}}(\gamma )){{dN}}_{0}/{\gamma }_{e}$ (Dai et al. 2002). In a numerical treatment of this process (Takahashi et al. 2011), the integration is performed between the locations of pair production and energy loss, provided that the deflection angle is sufficient for radiating emission into the observer's direction.

By contrast, in the MC approach used here we start from the distribution of pairs generated using Equation (4) and shown in Figure 2. We select those individual electrons whose associated delay time matches the observational criteria, and then we calculate their contribution to the echo radiation.

We assume the source emits photons with a power-law spectrum up to ${E}_{\mathrm{cut}}$, which we typically choose to be 3 or 30 TeV in our examples. Photons from this spectrum are absorbed by the EBL. We calculate the isotropic equivalent energy absorbed from the difference between the unabsorbed and the absorbed fluxes (as shown in the numerical calculations below).

$$\begin{eqnarray}&&{{ \mathcal E }}_{\mathrm{abs}}={\int }_{10\,\mathrm{GeV}}^{{E}_{\mathrm{cut}}}E\displaystyle \frac{{dN}(E)}{{dE}}(1-{e}^{-\tau (E,z)}){dE}.\end{eqnarray} \tag{ 8 }$$

The lower limit on the integration is set because the universe is transparent to 10 GeV photons at all redshifts since the epoch of galaxy formation (z ∼ 10). The number of pairs involved in our simulations is ${ \mathcal O }({10}^{5})$, but varies for different calculations. We use the total absorbed energy ${{ \mathcal E }}_{\mathrm{abs}}$ to scale our simulation to the actual differential distribution of pairs. The simulated differential pair distribution is related to the real distribution through a scale factor C from the expression

$$\begin{eqnarray}&&C{\int }_{1}^{\infty }d{\gamma }_{e}{({{dN}}_{{\gamma }_{e}}/d{\gamma }_{e})}_{\mathrm{sim}}{\gamma }_{e}{m}_{e}{c}^{2}=C\displaystyle \sum _{i=1}^{{N}_{\mathrm{pairs}}}{\gamma }_{i}{m}_{e}{c}^{2}={{ \mathcal E }}_{\mathrm{abs}}.\end{eqnarray} \tag{ 9 }$$

Next, we apply the observational criteria specifying start and stop times of observations. We choose the number of TeV photons, the cut energy ${E}_{\mathrm{cut}}$, $z,{B}_{\mathrm{IGMF}}$, and ${R}_{\mathrm{coh}}$. We calculate ${\lambda }_{\gamma \gamma }$ and ${\lambda }_{T}$. Using these values, we determine ${\theta }_{\mathrm{dfl}}$ and the time delay ${\rm{\Delta }}t$. We calculate the individual electrons' contribution to the spectrum and sum. Although we follow individual electrons, the resulting spectrum will be smooth because we convolve the electron distribution with the Compton kernel.

The observing angle θ of IC photons with energy ${E}_{\mathrm{GeV}}$ GeV that reach the observer can be calculated from the relation $\sin \theta =({\lambda }_{\gamma \gamma }/D)\sin {\theta }_{\mathrm{dfl}}$. In order to be observed by Fermi, this angle has to be less than the energy-dependent PSF of the LAT. For this, we use the expression ${\theta }_{\mathrm{PSF}}\approx 1\buildrel{\circ}\over{.} 68{E}_{\mathrm{GeV}}^{-0.77}\,+0\buildrel{\circ}\over{.} 2{e}^{-10{E}_{\mathrm{GeV}}^{-1}}$, which is an analytical approximation for the 95% containment PSF that is sufficiently accurate for our purposes (see Ackermann et al. 2013 for a detailed discussion of the Fermi PSF). This constraint is important for large IGMF values $\gtrsim {10}^{-17}$ G. Neglecting this effect can introduce inaccuracies into the low-energy part of the echo radiation spectrum.

To assess if our code behaves correctly, we calculate the time evolution of the echo flux for a cutoff of 3 TeV. Figure 3 shows results of our calculations for $B={10}^{-20}$ G. Here and throughout, we fix the coherence length of the magnetic field to ${R}_{\mathrm{coh}}=1$ Mpc, while recognizing that this is a poorly known quantity. Derivations with different ${R}_{\mathrm{coh}}$ can be calculated; however, since our data consist of upper limits, we focus this study on the IGMF for an assumed value of ${R}_{\mathrm{coh}}$ to demonstrate the method.

Zoom In Zoom Out Reset image size

The spectra at different times have a low-energy power-law segment followed by a break and a rapidly dropping flux (Figure 3, top panel). The echo spectrum at low energies behaves roughly as $\nu {F}_{\nu }\propto {E}^{1.5}$, corresponding to a cooling spectrum of a narrow electron distribution injected at high energies. The evolution of the break with time is a consequence of the more rapid reprocessing of higher-energy photons, which make higher-energy pairs that are less deflected by the IGMF, so that the higher-energy Compton-scattered photons have a shorter path length to the observer.

Figure 3 (bottom panel) shows the passing of the spectral peak in different energy bands. At lower energies the power-law slope is shallower and the break occurs later. The peak of the $\nu {F}_{\nu }$ spectrum decreases ∝t⁻¹. This behavior broadly follows from the fact that the total echo fluence (F × t) is determined by the amount of VHE flux absorbed by the EBL and is constant.

GRB 130427A had the largest γ-ray fluence of any GRB yet observed, with high-energy emission detectable up to ∼1 day after the burst trigger with Fermi-LAT. VERITAS followed up and derived an upper limit for the flux at 100 GeV at 0.82 days after the trigger (Aliu et al. 2014). Here, we use the VERITAS upper limit to place constraints on the IGMF. We use the emission episode from 11.5 to 33 s after the burst trigger, where the spectrum can be described as a power law with photon index ${\rm{\Gamma }}=-1.66\pm 0.13$, to define the γ-ray spectrum from this GRB. Since there was no sign of a spectral cutoff, we assume the spectrum extended up to 3 and 30 TeV to give estimates for the echo radiation flux.

In Figure 4, we calculate the echo radiation using our simulation and compare it with the VERITAS measurements. VERITAS provides a range of upper limits based on the assumed spectral shape. By comparing the VERITAS limits with the simulated spectra, we find that either ${B}_{\mathrm{IGMF}}\gtrsim {10}^{-17}$ G or ${B}_{\mathrm{IGMF}}\lesssim 3\times {10}^{-19}$ G for the ${E}_{\mathrm{cut}}=30\,\mathrm{TeV}$ case (otherwise the calculated echo flux would violate the VERITAS upper limits). The ${E}_{\mathrm{cut}}=3\,\mathrm{TeV}$ case does not give any constraints. Based on current estimates for CTA sensitivity, more stringent constraints are expected for observations similar to that of VERITAS (see Figure 4). HAWC is more likely to observe the direct prompt radiation than the echo flux.

Zoom In Zoom Out Reset image size

Another intriguing method to constrain the IGMF is provided by Fermi-LAT observations of GRB 130427A. The $E=32\,\mathrm{GeV}$ photon observed at ${\rm{\Delta }}t=34.4$ ks after the GRB trigger is difficult to interpret in the framework of synchrotron radiation from the forward shock (Ackermann et al. 2014). Nonetheless, it might be associated with an IC component, even though there is no evidence for the spectral and temporal break expected for a transition from synchrotron to Compton emission.

If we assume this photon originates from the echo radiation, there is a straightforward way to estimate the IGMF using Equation (3). A simple calculation yields ${B}_{\mathrm{IGMF}}=2.2\,\times {10}^{-19}$ ${({\rm{\Delta }}t/34.4\mathrm{ks})}^{1/2}{({\lambda }_{\gamma \gamma }/146\mathrm{Mpc})}^{-1/2}$ ${((1+z)/1.34)}^{4}(E/32\,\mathrm{GeV})\,{\rm{G}}$ with ${\lambda }_{\gamma \gamma }\approx 150\,\mathrm{Mpc}$ and $E\sim 6.3\,\mathrm{TeV}$ as the energy of the primary photon. Equation (3) can be solved by Monte Carlo methods to gauge the error on this IGMF value, and we get ${\mathrm{log}}_{10}{B}_{\mathrm{IGMF}}=-18.8\pm 0.3$. This represents the average of the simulated values and it differs somewhat from the estimated value, partly due to the asymmetry of the distribution of simulated values.

The delayed echo radiation potentially lasts for a long time compared to the time during which high-energy GRB afterglow radiation is expected. Since Fermi-LAT is an all-sky monitor, we explore the possibility that the echo radiation can be detected with long-exposure observations starting from one day after the GRB to several years, on the order of the lifetime of Fermi-LAT. We therefore calculated the flux of the echo radiation from GRB 130427A for observations starting one day after the GRB, where the direct GeV range afterglow radiation already went undetected. This observation spans 1000 days (see Figure 5).

Zoom In Zoom Out Reset image size

To compare the simulated echo spectrum with observations, we use Fermi-LAT data, specifically PASS 8 photons with P8R2_SOURCE_V6 filters. The diffuse galactic foreground was accounted for using gll_iem_v6, and the isotropic diffuse radiation by iso_P8R2_SOURCE_V6. We used the fermipy⁴ package and its routines to derive Fermi-LAT upper limits.

Fermi-LAT upper limits for the long-term (1000 days) observations in the direction of GRB 130427A are on the order of 10⁻¹³ erg cm⁻² s⁻¹ (see Figure 5 for the exact values). By comparing the upper limits to the MC simulations, we can rule out the cases where ${10}^{-18}\lesssim {B}_{\mathrm{IGMF}}({\rm{G}})\lesssim {10}^{-20}$. Stronger magnetic fields (${{\rm{B}}}_{\mathrm{IGMF}}\gtrsim {10}^{-17}$ G) do not violate the upper limits.

The power-law spectral components with photon index harder than −2 will have a cutoff energy. This cutoff energy of the delayed hard component is generally outside of the Fermi-LAT range. To gauge the effects of different cutoff energies, here we use Fermi-LAT observations starting from one day after the GRB trigger, where the afterglow has faded, ending at 1001 days after (similarly to the previous subsection). In Figure 6, we make a calculation of the echo radiation using the parameters of GRB 130427A, with varying cutoff energies, assuming a magnetic field strength of 10⁻¹⁸ G. This is near the lower limit of the IGMF inferred from observations of blazars (Dermer et al. 2011). We overplot the expected echo radiation with different cutoff energies and show that long-term LAT observations provide meaningful constraints assuming this value of B. In particular for $B={10}^{-18}$ G, ${E}_{\mathrm{cut}}\gtrsim 3\,\mathrm{TeV}$ of the predicted echo flux violates the LAT upper limits.

Zoom In Zoom Out Reset image size