<figure><img src="image_1.png" /></figure>

# Constraints on the Intergalactic Magnetic Field from Bow Ties in the Gamma-Ray Sky

Paul Tiede¹,²,³ ID, Avery E. Broderick¹,²,³ ID, Mohamad Shalaby¹,³,⁴,⁵ ID, Christoph Pfrommer⁴ ID, Ewald Puchwein⁴,⁶ ID,
Philip Chang⁷ ID, and Astrid Lamberts⁸ ID

¹ Department of Physics and Astronomy, University of Waterloo, 200 University Avenue West, Waterloo, ON, N2L 3G1, Canada; ptiede@perimeterinstitute.ca

² Waterloo Centre for Astrophysics, University of Waterloo, Waterloo, ON N2L 3G1 Canada; abroderick@perimeterinstitute.ca

³ Perimeter Institute for Theoretical Physics, 31 Caroline Street North, Waterloo, ON, N2L 2Y5, Canada; mshalaby@perimeterinstitute.ca

⁴ Leibniz-Institut für Astrophysik Potsdam (AIP), An der Sternwarte 16, D-14482 Potsdam, Germany; cpfrommer@aip.de

⁵ Department of Physics, Faculty of Science, Cairo University, Giza 12613, Egypt

⁶ Institute of Astronomy and Kavli Institute for Cosmology, University of Cambridge, Madingley Road, Cambridge, CB3 0HA, UK; epuchwein@aip.de

⁷ Department of Physics, University of Wisconsin-Milwaukee, 3135 N. Maryland Avenue, Milwaukee, WI 53211, USA; chang65@uwm.edu

⁸ Université Côte d'Azur, Observatoire de la Côte d'Azur, CNRS, Laboratoire Lagrange, Laboratoire ARTEMIS, France; astrid.lamberts@oca.eu

Received 2020 January 7; revised 2020 February 5; accepted 2020 February 5; published 2020 April 3

## Abstract

Pair creation on the cosmic infrared background and subsequent inverse-Compton scattering on the cosmic microwave background potentially reprocesses the TeV emission of blazars into faint GeV halos with structures sensitive to intergalactic magnetic fields (IGMF). Previous work has shown that these halos are then highly anisotropic and extended. If the coherence length of the IGMF is greater than the inverse-Compton cooling length of the pairs, then the orientation of the gamma-ray halo will be correlated with the direction of the magnetic field which is unknown and expected to change for each source. In order to constructively add each source we then use angular power spectra which are insensitive to the jet orientation. By looking at known GeV blazars detected by Fermi, we exclude the existence of an IGMF with coherence lengths >100 Mpc at greater than 3.9σ with current-day strengths in the range 10⁻¹⁶ to 10⁻¹⁵ G, and at 2σ from 10⁻¹⁷ to 10⁻¹⁴ G. This provides a direct measurement of the nonexistence of gamma-ray halos, providing an important check on previous results.

Unified Astronomy Thesaurus concepts: BL Lacertae objects (158); Gamma-rays (637); High energy astrophysics (739)

## 1. Introduction

Astronomical magnetic fields appear on scales ranging from the terrestrial to those of galaxy clusters, the latter spanning many megaparsecs. Within all of these systems magnetic fields are believed to play a variety of important roles. These include providing additional pressure components (Beck 2016), mediating the magnetohydrodynamic turbulence (Goldreich & Sridhar 1995), providing and limiting the effective viscosity that transports angular momentum in accretion disks (Balbus & Hawley 1991) and mediates collapse (Crutcher 2012; Hennebelle & Falgarone 2012), contributing to the acceleration of high-energy particles, e.g., cosmic rays (Blandford & Eichler 1987), and dictating their subsequent propagation (Hanasz & Lesch 2003; Pakmor et al. 2016) as well as diffusion of thermal energy and momentum in weakly collisional and collisionless systems (Schekochihin et al. 2009).

Observed astrophysical magnetic fields are likely amplified from weak primordial seed fields. These seed fields may originate from Biermann battery processes (Kulsrud & Zweibel 2008), electroweak or quantum chromo-dynamics phase transitions in the early universe, or during cosmic inflation. Biermann battery processes can naturally generate magnetic fields when temperature and density gradients are not parallel, resulting in electric fields that have a curl and providing a source for the magnetic field. Typically, these processes have a correlation length that reflects the length scales associated with the temperature and density gradients (Kulsrud & Zweibel 2008; Widrow et al. 2012), and generate small-scale (~10 kpc) fields (B ~ 10⁻²⁰ G), which can be amplified by the action of a dynamo to µG fields observed in galaxies today (Schober et al. 2013; Beck 2016; Pakmor et al. 2017).

By their nature, fields generated in primordial processes may be present in the intergalactic medium as the intergalactic magnetic fields (IGMF). While the general nature of these fields in still unknown, causal constraints limit the coherent length scales of primordial fields generated from post-inflation mechanisms, e.g., phase transitions, to small comoving scales, typically <10 kpc (Widrow et al. 2012). Small-scale fields may seed large-scale fields through the inverse cascade up to ~0.1 Mpc, i.e., the sonic scale of 10⁴ K gas over the age of the universe (Ryu et al. 2012).

To generate large-scale (≥1 Mpc) volume-filling fields, seed fields must have been generated prior to or during inflation (Biermann 1950; Turner & Widrow 1988) or from second-order interactions between the photon and electron fluid prior to recombination (Ichiki et al. 2006). In the former case, additional physics of superadiabatic magnetic amplification (Tsagas & Kandus 2005) is needed to keep these fields relevant today. This leaves the strength of these inflationary fields unconstrained (Widrow et al. 2012). In the latter case, the expected comoving field strength is B ~ 10⁻²⁴ G.

Probing these primordial fields is difficult, but some constraints can be set on a large-scale IGMF. The best direct upper limits on the strength of a putative large-scale IGMF to date come from the impact of such fields on the cosmic microwave background (CMB), roughly 10⁻⁹ G, followed closely by limits arising from Faraday rotation measurements, roughly 10⁻⁹ to 10⁻⁸ G depending weakly on the coherence length of the field λB (Neronov & Vovk 2010). Indirect limits on large-scale IGMFs can be obtained from the rotation measures of galaxy clusters, constraining the IGMF on scales exceeding 1 Mpc to be less than 10⁻¹² G (Dolag et al. 2005; Neronov & Vovk 2010). In contrast,
indirect lower limits can be placed by the observed absence of inverse-Compton gamma-ray cascades.

Very high-energy gamma-rays (VHEGRs), with energies exceeding $E > 100$ GeV, lie above the pair-creation threshold with the infrared background (Gould & Schréder 1967; Salamon & Stecker 1998). As a result, the mean free path of VHEGRs in the nearby universe is $10^2(1+z)^{-4.5}(E/6 \text{ TeV})^{-1}$ Mpc (Broderick et al. 2012). That is, while VHEGRs can traverse intergalactic distances, they are absorbed on cosmological distances, producing highly relativistic beams of electron-positron pairs (i.e., Lorentz factors of order $10^6$). This redshift-dependent absorption of high-energy gamma-rays has been observed (Ackermann et al. 2012a) and subsequently used to probe the infrared background (Domínguez et al. 2011, 2013; Gilmore et al. 2012).

In the absence of any other process, e.g., collective beam-plasma instabilities (Broderick et al. 2012), these beams will inverse-Compton cool by upscattering CMB photons to GeV energies. That is, this chain of events would effectively reprocess the TeV emission of active galactic nucleus (AGN) into GeV emission. Thus, bright, nearby, extragalactic TeV sources should exhibit a GeV bump that is correlated with their TeV spectral energy distribution (SED). Stringent limits have been placed by Fermi on the existence of such a GeV bump already, effectively excluding such a feature in several TeV sources (Neronov & Semikoz 2009; Neronov & Vovk 2010).

Presuming that inverse-Compton cooling dominates the evolution of the beams, these nondetections provide indirect, circumstantial evidence for the existence of a strong IGMF, with $B > 10^{-15}$ to $10^{-17}$ G, depending on the particular extragalactic background light model employed (Neronov & Vovk 2010; Taylor et al. 2011; Vovk et al. 2012). This assumes that the duty cycle of the gamma-ray blazars is comparable to that observed for the radio mode in AGN, roughly $10^6$ yr; shorter duty cycles could result in GeV echos that lag behind the initial VHEGR emission by long times, resulting in a weaker limit: $B > 10^{-17}$ G (Dermer et al. 2011; Taylor et al. 2011). Such an IGMF would deflect the resulting pairs, whose inverse-Compton emission would then be highly beamed away from us, naturally removing the GeV bump. However, at the same time, this would open the door to the direct detection of the IGMF through the observation of GeV halos surrounding intrinsically gamma-ray bright objects, for which the inverse-Compton emission is beamed into our line of sight.

A number of searches for inverse-Compton halos using the publicly available Fermi observations have been performed, yielding contradictory claims but no convincing evidence for their presence (Ando & Kusenko 2010; Neronov et al. 2011; Ackermann et al. 2013; Chen et al. 2015). These have focused on the large radial extent of the gamma-ray halos, typically subtending many degrees. However, this is generally degenerate with the point-spread function (PSF) of Fermi's Large Area Telescope (LAT), which until recently was highly uncertain (Neronov et al. 2011; Ackermann et al. 2012b). As result, even when detections are claimed they remain plagued by large systematic uncertainties.

Previous work in Broderick et al. (2018) stacked, rotated Fermi images of Fanaroff-Riley class I and II objects from the VLA FIRST survey (Helfand et al. 2015), to limit the strength of the IGMF to $<10^{-15}$ G for large correlation length magnetic fields. A key difference in this work is that we use Fermi blazars and their full spectral information in the GeV-TeV

range (see Tiede et al. 2016) to attempt to measure gamma-ray halos. This is different from the work by Broderick et al. (2018) who statistically sampled the blazar SEDs using the distribution of Fermi and VHEGR spectra. Due to this, the work in this paper provides an additional and independent constraint on the IGMF strength.

The layout of this paper is as follows. In Section 2, we review our method first described in two companion publications (Broderick et al. 2016; Tiede et al. 2016), where we used angular power spectra to coherently stack angular structure from gamma-ray halos. In Section 3 we describe a model for the gamma-ray sky background, which adds power at low multipoles. Section 4 details the main results followed by a discussion of the impact of various potential systematics. Finally, Section 5 summarizes the conclusions of this work.

## 2. Method

In two companion publications (Broderick et al. 2016; Tiede et al. 2016) we have demonstrated that these inverse-Compton halos are expected to be highly anisotropic (see also Neronov et al. 2010; Long & Vachaspati 2015). This originates from different physical origins for a small-scale, tangled IGMF ($\lambda_B \ll 3$ Mpc) and a large-scale, uniform IGMF ($\lambda_B > 100$ Mpc). For the latter, of interest here, this is due to extreme beaming of the inverse-Compton emission coupled with the geometric conditions upon the gyration of the particles. This typically produces a bilateral, narrowly beamed halo, the two lobes being due to the opposite gyration directions of the electrons and positrons (Broderick et al. 2016).

In Tiede et al. (2016) we demonstrated that existing Fermi observations should permit the detection of a large-scale IGMF with strengths comparable to those implied by the absence of GeV bumps in nearby TeV sources, i.e., $10^{-16}$ to $10^{-15}$ G. This makes use of stacked one-dimensional angular power spectra that amplify the impact of anisotropy<fnref n="9" />:

$$P_m = \frac{1}{N_{\text{src}}} \sum_j \frac{1}{N^2} \left| \sum_j e^{im\theta_j} \right|^2, \quad (1)$$

where $\theta_j$ is the polar angle of the jth photon defined relative to the position of the central source and an arbitrary direction, and $N$ the number of photons for the source. For this we make use of the Pass 8R2_V6 ULTRACLEANVETO event class and limit ourselves to gamma-ray energies between 1 and 100 GeV to avoid large variations in the instrument response. We further apply an energy-dependent and conversion-location-dependent mask (i.e., Front and Back), equal to the 68% inclusion region of the Pass 8R2_V6 PSF, thus excising the region most dominated by the direct emission from the source.

The existence of a bilateral, jet-like feature appears primarily as an excess of quadrupolar power; with a sufficient number of photons it would produce a sawtooth-like structure in the power spectrum in which even multipoles exhibit an excess. Based upon the source SEDs and Monte Carlo modeling of the Fermi source population we have identified sets of Fermi objects that are optimized for the detection of the anisotropic halo feature. These necessarily depend upon the assumed strength of the IGMF, with larger IGMFs producing observable

<footnote n="9">⁹ Note that Broderick et al. (2018) used the VLA survey so that the jet direction is known, so stacking angular histograms coherently was possible.</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Aitoff projection of the 18 source (red) and 507 high-latitude background (gray) field locations shown in equatorial coordinates. The region with $|b| \le 10^\circ$ is shown in blue and the location of the Galactic center is shown by the asterisk. For comparison, the location of 1ES 0229+200, the primary source employed by Taylor et al. (2011), is shown by the bright red star; the bright red circle indicates the projected 100 Mpc region about 1ES 0229+200. All field locations are circles when viewed directly, the apparent asymmetries are due solely to projection.</figcaption></figure>

halos in more sources and therefore permitting a larger sample that may be profitably stacked. These optimized source lists may be found in Table 1 of Tiede et al. (2016), with the number of sources ranging from 4 to 18 for current-day strengths of $B_0 = 10^{-17}$ to $10^{-14}$ G (note that the source-frame field strength is given by $B = B_0(1+z)^2$, though due to the small redshifts of the gamma-ray sources employed the two differ by at most a factor of 2). For each, an uncontaminated region about each source has been identified via visual inspection. Note that as a result these fields do typically contain many dim sources identified in the Fermi-LAT Third Source Catalog (3FGL Acero et al. 2015).

Critically, the anisotropy signal is unlikely to be confused with a number of anticipated systematic contaminants. In principle, these may arise from unknown physics at the TeV source, the LAT response, or unmodeled structures in the gamma-ray sky. In practice, chief among the potential contaminants is the degree-scale substructures in Galactic contribution to the diffuse gamma-ray background. To assess the degree to which these contribute to the stacked angular power spectra we now construct a model of the gamma-ray background.

### 3. Modeling the Background Contributions to the Power Spectra

In this section we compute the contributions of the photons from background and contaminating 3FGL sources to the stacked angular power spectra. To assess, quantify, and guide the modeling of said contribution, we first identified 507 fields with Galactic latitudes more than 20° from the Galactic plane. The locations of these are shown in Figure 1. Also shown for comparison are the source fields (red). Note that while the background fields appear to nearly cover the sky, the bright point sources in the Fermi sky have effectively been excised according to the criteria used in Tiede et al. (2016) to determine the size of the source fields. That is, all of the bright sources fall into the gaps between the source and background fields. Note that these fields are not devoid of 3FGL sources—sources with fluences less than roughly 40 photons do appear in these background fields and in the source fields. The stacked power

spectrum of these 507 fields comprise the background power spectrum that we sought to successfully model.

The model we constructed consists of four components:

1. 3FGL sources.

2. Completion of the 3FGL below the survey sensitivity limit.

3. Large-scale gradients (>9°).

4. Small-scale gradients (≈2°).

We will discuss the construction of each individually.

#### 3.1. Inspecting the Fermi Sky

The construction of all four components is facilitated by the direct inspection of the Fermi sky, obtained from 239557417 s (the official start of the Fermi data release) to 487012885 s in mission elapsed time. A density map of the 1–100 GeV photons in the Front- and Back-converted Pass 8R2_V6 ULTRACLEANVETO events class is smoothed with a Gaussian kernel. The size of this kernel is informed by typical structure of the point sources, set by the Pass 8R2_V6 PSF, below which no discernible structure appears. Thus, as our finest resolution map we choose a kernel full-width, half-max (FWHM) of 1°2, which is shown in the top row of Figure 2. The multitude of point sources that comprise the 3FGL are immediately evident.

#### 3.2. Component 1: 3FGL Sources

Because the period from which the 3FGL was constructed is roughly half that over which these maps are generated, both exposure differences and source variability will produce significant differences in the anticipated number of photons from each 3FGL object. Therefore, we begin by refitting the individual 3FGL source photon fluences.

We do this by comparing the observed gamma-ray maps smoothed on 1°2 and 4°7 scales to the expected response from the individual point sources. That is, we begin by computing the observed photon density at each of the 3FGL source locations, given by $n_j^{1.2}$ and $n_j^{4.7}$, respectively, for the $j$th source. We then simulate the corresponding anticipated fractional response of each individual source at each of the 3FGL positions, i.e., the number density of photons located at the $j$th source location associated with $N_i$ photons from the $i$th source are $f_{ji}^{1.2}$ $N_i$ and $f_{ji}^{4.7}$ $N_i$ for the two smoothing scales, respectively. In doing this we employ the Pass 8R2_V6 ULTRACLEANVETO PSFs and assume the spectral index listed in the 3FGL.

In addition to the point sources in the 3FGL there exists a diffuse component associated with gamma-ray emission from Galactic leptons and hadrons. Unlike point sources, this contribution should be similar at both smoothing scales. Thus, we have

$$n_j^s = \bar{n}_j + \sum_i f_{ji}^s N_i \qquad (2)$$

for some fixed set of diffuse background photon densities, $\bar{n}_j$, at both smoothing scales. This may be concurrently solved for the $\bar{n}_j$ and $N_i$, yielding the set of linear equations for the latter,

$$\sum_i (f_{ji}^{1.2} - f_{ji}^{4.7}) N_i = (n_j^{1.2} - n_j^{4.7}). \qquad (3)$$

We solve this for the $N_i$, and subsequently verify that the implied photon density associated with the diffuse background
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Aitoff projected maps in Galactic coordinates of the gamma-ray sky at 1–100 GeV for events converted in the front (left) and back (right) detectors. Top: the flux density smoothed with a Gaussian filter with an FWHM of 1°2. Upper middle: the flux density smoothed with a Gaussian filter with an FWHM of 1°2 after subtracting the 3FGL sources with their refitted fluences. Lower middle: the difference of 3FGL-subtracted maps smoothed over 1°2 and 9°4, where the latter is a measure of the large-scale structure in the map, producing an image of the small-scale structure after removal of the 3FGL. Bottom: the fractional difference of the 3FGL-subtracted maps smoothed over 1°2 and 9°4, where the latter is taken as the normalization. In all plots lines of $|b| = 20°$ are shown in white, and the region for which $|b| \le 10°$ has been masked out to remove the Galactic plane.</figcaption></figure>

is positive. While this method for estimating the 3FGL source fluences does account for the differing exposures and source variability, it does not account for the possibility of new sources. It is also poorly suited to correcting the fluences of extended sources. We do not anticipate that either of these limits will present significant obstacles.

With the individual source fluences we then construct a realization of the 3FGL, smooth it over the appropriate scale, and subtract it from the observed photon density. This is shown for a smooth scale of 1°2 for Front- and Back-converted Pass 8R2_V6 ULTRACLEANVETO events in the second row of Figure 2. Apart from a handful of extended Galactic sources,
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Flux distribution of 3FGL sources with Galactic latitudes above 20° (black). For comparison the flux distribution of 3FGL sources without any Galactic latitude cut (blue points) and an extreme latitude cut ($|b| > 60°$, red points) are also shown; these are offset for clarity. In green we show the single power-law model we assume for the intrinsic flux distribution. The purple dotted line shows the assumed detection efficiency, $w(F_{35})$, and the dark-green dashed line shows the associated expected flux distribution, both evaluated at $b = 90°$.</figcaption></figure>

and up to Poisson fluctuations in the number of sources for the Back-converted map, point sources have been reasonably well removed. In practice, because the background fields specifically exclude bright sources, the details of the 3FGL removal are relatively unimportant. Simply subtracting the 3FGL sources assuming the 3FGL fluxes and fixed exposure produces qualitatively identical results.

### 3.3. Component 2: Completion of the 3FGL

The second component, the unresolved completion of the 3FGL, can only be inserted in a statistical sense. Ultimately, a complement of sources with the appropriate fluxes are selected and randomly placed within the background/source fields in a manner consistent with the 9°4-smoothed background. The source fluxes are chosen from the power-law distribution shown in Figure 3:

$$N^{-1} \frac{dN}{d \log_{10} F_{35}} = 0.45 \left( \frac{F_{35}}{10^{-9} \text{ ph s}^{-1} \text{ cm}^{-2}} \right)^{-1.1} \quad (4)$$

This is supplemented with a detection efficiency, which we assume has a log-normal cutoff, i.e.,

$$w(F_{35}) = \begin{cases} \exp[-12.88(\log_{10} F_{35}/F_{\text{th}})^2] & F_{35} < F_{\text{th}} \\ 1 & \text{otherwise.} \end{cases} \quad (5)$$

While for the $|b| > 20°$ population the threshold flux, $F_{\text{th}}$, is well described by a single number, in detail it does depend on Galactic latitude. This is clear in Figure 4, in which there is a trend in the lower envelope of 3FGL fluxes with Galactic latitude, scaling roughly as $|b|^{-0.33}$. Therefore, we set

$$F_{\text{th}} = 4.27 \times 10^{-10} |b/90^{\circ}|^{-0.33} \text{ ph cm}^{-2} \text{ s}^{-1}. \quad (6)$$

Note that, despite the fact that the assumed flux distribution produces a formally infinite flux (and is thus unphysical), the condition that photons are observed in integer numbers, i.e., at least one, limits the number of photons observed in practice.

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Flux vs. Galactic latitude for the 3FGL catalog. Blue points indicate the estimated position of the lowest-flux quartile in the indicated latitude bin; error estimates are obtained via bootstrap and incorporate the underlying flux uncertainties. The threshold for detection evolves with latitude as $b^{-0.33}$, shown for an appropriate normalization by the red line.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 5.** Stacked angular power spectra of the background fields from the data (black), and for realizations drawn from Gaussian smoothed background maps with various smoothing scales. As the smoothing scale increases power is suppressed at progressively smaller $m$. Power at $m \le 2$ is associated with structures on roughly a few degree scales.</figcaption></figure>

That is, for a flux $F_{35}$ the probability that at least a single photon will be observed is

$$p_1(F_{35}) = 1 - p_0(F_{35}) = 1 - e^{-F_{35}\epsilon} \approx F_{35}\epsilon, \quad (7)$$

where $\epsilon$ is the exposure. Thus, at small fluxes, the probability of observing at least one photon is linearly suppressed. Therefore, any flux distribution flatter than $F_{35}^{-2}$ produces a convergent number of photons.

In practice, the contribution of the unresolved completion of the 3FGL is negligible.

### 3.4. Component 3 and 4: Gradient Fields

The third component consists of randomly drawing the remaining background photons with a spatial distribution that follows the 9°4-smoothed map. This modifies the dipole component significantly, though again makes a negligible difference at $m > 1$ as shown in Figure 5.

Finally, the fourth component consists of an attempt to model the degree-scale fluctuations in the diffuse Galactic background. Stacked angular power spectra generated from
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Stacked angular power spectra for background fields above a given Galactic latitude from the background fields from the data (red), and for simulated realizations that include the 3FGL, an approximation of its completion, and the diffuse structure on scales above 9°4 (blue). The comparison becomes increasingly good as the Galactic latitude increases, strongly suggesting that the disparity is a result of structure in the Galactic foreground.</figcaption></figure>

realizations of the background fields using smoothed maps of the gamma-ray sky are shown in Figure 5. Unsurprisingly, larger smoothing scales correspond to less power at large *m*. Specifically, at smoothing scales of 1°2 the power above *m* = 5 is noticeably suppressed. Beyond smoothing scales of 2°4 this suppression extends to *m* ≥ 3; by 9°4 only the dipole power remains. This suggests that any missing power at *m* = 2 is most sensitive to a few degree-scale structures in the diffuse gamma-ray background. Therefore, to clearly identify these we generate difference maps between high (1°2) and low (9°4) smoothing scales. These are shown in the third row of Figure 2.

Generally, these exhibit a latitude dependent amplitude, falling rapidly with distance from the Galactic plane. This supports their association with Galactic structures. A similar conclusion holds regarding their impact on the stacked angular power spectra. Figure 6 shows a direct comparison between the observed and simulated stacked power spectra, including only the first three components, from the background fields assuming different latitude cuts. When restricted to fields with $|b| > 60°$ the two are very similar, as anticipated by the large reduction of the Galactic diffuse emission at high latitudes.

In contrast to their absolute amplitudes, the fractional amplitudes in the photon densities are roughly fixed, however (shown in the bottom row of Figure 2). The PDF of the fractional density fluctuations is well described by the normal distributions, shown in Figure 7. The mean is offset significantly from zero.

<figure><img src="image_8.png" /><figcaption>**Figure 7.** Fractional flux probability distribution function for fractional fluctuations in the 3FGL-subtracted maps of Front-converted (top) and Back-converted (bottom) events, smoothed over 1°2 and 9°4. A Gaussian fit is shown in red for each, with mean -0.035 and -0.033, standard deviation 0.114 and 0.115, for the Front- and Back-converted event distributions, respectively. In both cases the relative norm of the Gaussian is 0.94. For comparison, the blue dotted line is the PDF of fractional fluctuations obtained without subtracting the 3FGL sources, giving an indication of the impact of the unresolved completion of the 3FGL. That is, the negative mean and substantial tails at large deviations are likely to be associated with the incomplete source subtraction.</figcaption></figure>

This is a natural result from imperfect point source subtraction, likely arising from the failure to subtract unidentified sources. This is supported by the PDF of fractional density fluctuations constructed without any point source subtraction, shown in the dotted blue lines in Figure 7, which also shows a noted shift to negative fractional flux fluctuations.

The impact of the degree-scale structures in the Galactic diffuse emission is to increase the power at low *m*. We quantify this by generating two stacked angular power spectra for the background fields from over-resolved realizations of the gamma-ray sky that include only the initial three components: 3FGL, its completion, and large-scale gradients (Figure 8). These differ in the maps used to generate the last component, the large-scale gradients; one using the map smoothed on 9°4 and the other smoothed on 1°2. The latter map is statistically indistinguishable from the stacked angular power spectrum obtained from the data directly. The difference between the two is entirely due to the otherwise unresolved degree-scale structures of interest.

### 3.5. Constructing the Background Power Spectra

We model the angular power due to these small-scale background fluctuations statistically as a multipole-dependent, log-normal additional component with means and standard deviations chosen to reproduce the difference in power and difference in the variances of the low- and high-resolution angular power spectra. That is, we assume that the power in the high- and low-resolution angular power spectra, $P_m^h$ and $P_m^l$,
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Comparison of the stacked angular power spectra for simulated background fields using the smoothed, 3FGL-subtracted sky maps smoothed on 1°2 and 9°4 to reconstruct the degree-scale Galactic contribution. To prevent the Poisson noise from dominating the estimate, we randomly draw 100 times more photons from the flux maps than present in the background fields. We model this component as a multipole-dependent log-normal additional component with means and standard deviations chosen to reproduce the difference in power and difference in the variances of the low- and high-resolution angular power spectra. The black triangles show the angular power spectra of the Galactic background of the 3FGL source subtracted *Fermi* sky.</figcaption></figure>

respectively, are related via

$$P_m^h = P_m^l + P_m^{ss} \qquad (8)$$

where the small-scale structure contribution $P_m^{ss}$ is a log-normal random variable, i.e., $P_m^{ss} = \delta_m e^{\sigma_m X}$ for a zero-mean, unit-variance Gaussian random variable $X$. The values of $\delta_m$ and $\sigma_m$ are fully defined by the means and variances of the $P_m^h$ and $P_m^l$; it is straightforward to show that choosing at each $m$

$$\delta_m = \sqrt{\frac{\Delta_m^4}{\Delta_m^2 + \Sigma_m^2}}, \quad \sigma_m^2 = \ln\left(\frac{\Sigma_m^2 + \Delta_m^2}{\Delta_m^2}\right), \qquad (9)$$

where $\Delta_m$ is the difference in the mean powers and $\Sigma_m^2$ is the difference in the variances at high and low resolution at a given $m$, uniquely reproduces the means and variance of $P_m^h$ and $P_m^l$. This is done independently for angular power spectra constructed using Front-converted events, using Back-converted events, and using both Front- and Back-converted events, resulting in slightly different choices for $\delta_m$ and $\sigma_m$. For the last case, the resulting $\delta_m$ are shown in Figure 8 for the red and blue points. This is compared to the power spectra for the Galactic background of the 3FGL source subtracted *Fermi* sky (black points).

In summary, we have modeled the background power with a combination of the contaminating 3FGL sources, a random realization of the unresolved completion of the 3FGL source population, the large-scale (>9°) gradients, and random realizations of the small-scale gradients. At very high Galactic latitudes ($|b| > 50°$) this last component becomes insignificant, strongly implicating the Galactic contribution as its source. The resulting background power spectra are shown in Figure 9, subdivided into their origins.

# 4. Results

## 4.1. Comparing the Power Spectra to the Background

The stacked angular power spectrum for the optimized source set for $B_0 = 10^{-15}$ G is shown in Figure 10. As anticipated

<figure><img src="image_10.png" /><figcaption>**Figure 9.** Comparison of the observed (red) and reconstructed (black) stacked angular power spectra for the background fields. The formal $1\sigma$ errors for all 507 background fields are shown. The reconstructed power spectra are constructed using simulated realizations that include the 3FGL, an approximation of its completion, the diffuse structure on scales above 9°4, and the degree-scale background fluctuations.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 10.** Comparison of the stacked angular power spectra for the source (optimized for $10^{-15}$ G) and 507 background fields. The latter have been linearly shifted to adjust for the differing Poisson noise level and the error bars indicate the formal $1\sigma$ errors for 15 fields assuming that these are Gaussian (for the subsequent simulated backgrounds, this is not found to be the case).</figcaption></figure>

by the power spectrum from the background fields, there is excess power at low $m$. However, this is entirely consistent with that expected for 15 background fields, adjusted via a constant shift for the different Poisson noise level (red points in Figure 10). Thus, there is no evidence for any excess power at low multipoles.

This is further evident in Figure 11, which shows the probability distribution of the simulated and observed power about the simulated median. Not only is each value consistent with being drawn from the background, the distributions about the mean is itself similar. Conversely, the ability to reproduce the statistics of the observed stacked angular power provides additional confidence in our ability to adequately model the background model. This remains the case regardless of which optimized source sample is used.

Additionally, we expect different low $m$ to be correlated when pair halos exist (Tiede et al. 2016). While we have made no attempt to quantify the discrepancies, we found that the cross-correlations are consistent with the expected null distribution, and qualitatively look discrepant in the model that accounts for inverse-Compton halos (see the Appendix).
<figure><img src="image_12.png" /><figcaption>**Figure 11.** Probability distribution of the observed power about the mean of the expected power for the $B_0 = 10^{-15}$ G optimized source sample for $m < 10$ (brown) and $m < 100$ (orange). For reference, the distribution observed in the simulated realizations is also shown by the solid purple line. The three are statistically indistinguishable.</figcaption></figure>

## 4.2. Quantifying the Lack of Gamma-Ray Halos

To quantitatively assess the limits upon the IGMF strength that may be placed by the nondetection of excess structure, we employ our background model to generate multiple realizations of optimized *Fermi* source lists, with and without inverse-Compton halos present. From these we identify the one-sided confidence level at which a given value of the power may be excluded, i.e., the confidence level (CL), in each direction. Because the stacked powers at different multipoles are strongly correlated in the presence of inverse-Compton halos (see the Appendix), we focus on the quadrupole exclusively when reporting constraints. In principle, we can improve upon these by including information at multiple *m*; in practice, this serves to obscure the nature of the test for marginal benefit. Note that because the optimized source lists vary with $B_0$, this procedure must be repeated for each IGMF strength.

In the top panel of Figure 13 we show the 2σ (95.5%) and 3σ (99.8%) one-sided confidence levels for each multipole with (red) and without (purple) an inverse-Compton halo with $B_0 = 10^{-15}$ G. Even without modeling and including the background in the power spectra, we can exclude a large-scale IGMF with strength $B_0 = 10^{-15}$ G by 2σ. When the background is included we can exclude our IGMF model at 3.9σ. Figure 12 shows this result does not depend on which detector, front or back, is used, suggesting this is not from telescope systematics.

Moving to $B_0 = 10^{-18}$ G, $10^{-17}$ G, $10^{-16}$ G, $10^{-14}$ G, and $10^{-13}$ G, we generally find similar results. For $B_0 = 10^{-18}$ G and $10^{-13}$ G we used the optimized source list described in Broderick et al. (2016) for $B_0 = 10^{-17}$ G and $10^{-14}$ G, respectively (for their power spectra see Figure 13). These may be excluded at 1.9σ, 2.4σ, 4.0σ, 2.0σ, and 1.6σ, respectively. That is, for present-day field strengths between $10^{-16}$ and $10^{-15}$ G a large-scale IGMF can be excluded at more than 3.9σ. At greater than 2σ present-day field strengths slightly more than $10^{-18}$ to $10^{-14}$ G a large-scale IGMF can be excluded.

<figure><img src="image_13.png" /><figcaption>**Figure 12.** Power spectra for gamma-ray halos with $B_0 = 10^{-15}$ G and the null case split into front (top) and back (bottom) detector contributions. The solid lines show the one-sided 95% confidence intervals, while the dotted lines show the 99% confidence intervals. No significant difference is observed between front and back detectors.</figcaption></figure>

These regions are shown in Figure 14 in comparison to the limits in Ackermann et al. (2018) and Broderick et al. (2018). Beyond those limits the inverse-Compton halos either become sufficiently small that they are confused with the source (small $B_0$) or extended that they lie beyond the regions about the sources considered (large $B_0$).

## 5. Conclusions

Combined with the SED lower limits, and making modest assumptions regarding the source duty cycles, we can exclude at
<figure><img src="image_14.png" /><figcaption>**Figure 13.** Power spectra for halo models with the magnetic field strengths considered, i.e., $B_0 = 10^{-13}$ G, $10^{-14}$ G, $10^{-15}$ G, $10^{-16}$ G, $10^{-17}$ G, $10^{-18}$ G. In all instances the power spectra are consistent with the null hypothesis that no gamma-ray halos exist. The null power spectra are different in each panel since they use different source samples. For the source list see Tiede et al. (2016).</figcaption></figure>

more than $3.9\sigma$ the existence of an IGMF of $10^{-16} < B_0 < 10^{-15}$ G, and at $2\sigma$ field strengths of $10^{-17} < B_0 < 10^{-14}$, with $\lambda_B > 100$ Mpc. The constraints presented in this work complement the results of Broderick et al. (2018) and improve the limit of large-scale magnetic fields to be $<10^{-17}$ G at $2\sigma$. However, in this paper we used known gamma-ray bright sources, i.e., FSRQ and BL-Lac objects, from the *Fermi* 3FGL source list. As a result, while the constraints on the IGMF in this paper are statistically weaker than the Broderick et al. (2018), they require less assumptions about the nature of AGN and their high-energy emission. Namely, the high-energy spectrum for each source is measured and we do not require AGN unification. Therefore, these results form an independent test of the IGMF when compared to Broderick et al. (2018) and is more direct. Furthermore, we have demonstrated that even if the orientation of the gamma-ray halo is unknown, i.e., the magnetic field

orientation for each source is unknown, we can still constructively stack sources. This work strongly constrains models of magnetogenesis that invoke inflationary processes in the early universe, which generate fields on large scales.

It is nevertheless worth pointing out that all of the above limits can be avoided if inverse-Compton cascades are preempted by other physical processes. Chief among these are large-scale beam-plasma instabilities in the intergalactic medium. The linear growth rates for these far exceed the inverse-Compton cooling times, a conclusion that appears to be robust even when higher order perturbative effects are included (Broderick et al. 2012; Schlickeiser et al. 2012; Chang et al. 2014, 2016). This may have important implications for the gamma-ray sky (Broderick et al. 2014a, 2014b), the thermal history of the intergalactic medium, and cosmic structure formation in low-density environments (Chang et al. 2012;
<figure><img src="image_15.png" /><figcaption>**Figure 14.** Significance with which various large-scale ($\lambda_B > 100$ Mpc) IGMF strengths may now be excluded. The blue bars show the limits presented here. The two red regions show the limits from Ackermann et al. (2018). The dark red assumes a duty cycle of 10 yr while, the lighter red is for $10^7$ yr. Above $10^{-15}$ G Broderick et al. (2018) was able to exclude the existence of large-scale magnetic fields.</figcaption></figure>

Pfrommer et al. 2012; Puchwein et al. 2012; Lamberts et al. 2015). However, it does remain debated whether these instabilities saturate at sufficiently large amplitudes to efficiently extract the energy of the pairs in the nonlinear regime (Miniati & Elyiv 2013; Chang et al. 2014, 2016; Sironi & Giannios 2014; Kempf et al. 2016; Shalaby et al. 2017; Vafin et al. 2018, 2019).

A.E.B., P.T., and M.S. received financial support from the Perimeter Institute for Theoretical Physics and the Natural Sciences and Engineering Research Council of Canada through

a Discovery Grant and Alexander Graham Bell scholarship (PT). Research at Perimeter Institute is supported by the Government of Canada through Industry Canada and by the Province of Ontario through the Ministry of Research and Innovation. C.P., M.S., and E.P. acknowledge support by the European Research Council under ERC-CoG grant CRAGS-MAN-646955. P.C. is supported by the NASA ATP program through NASA grant NNH17ZDA001N-ATP. A.L. acknowledges support by the Programme National des Hautes Energies (France). E.P. acknowledges support by the Kavli Foundation.

Software: *astroy* (Price-Whelan et al. 2018), *matplo-tlib* (Hunter 2007), *scipy* (Jones et al. 2001).

## Appendix
Power Correlations

Generally the powers at different multipoles are highly correlated when inverse-Compton halos are present (Tiede et al. 2016). As a result, we have not made an attempt to leverage the observed broad inconsistency with the halo model expectation at low *m*. Nevertheless, in Figure 15 we show the joint probability distributions of the angular power in various low *m* multipoles in comparison to the observed values. From the comparison when inverse-Compton halos are present (Figure 15, left) it is clear that the observed powers lie well outside the joint distribution for all instances shown. That is, the observations are even more highly inconsistent than implied by the quadrupolar comparison. In contrast, when inverse-Compton halos are absent, i.e., our null case (Figure 15, right), the discrepancy disappears.

<figure><img src="image_16.png" /><figcaption>**Figure 15.** Correlations in stacked angular power spectrum of the inverse-Compton halo (left) and Null (right) models (color map), in comparison to the observed values for a large-scale $10^{-15}$ G IGMF (white star). In all panels the red dashed line shows the one-to-one relationship.</figcaption></figure>
ORCID iDs

Paul Tiede https://orcid.org/0000-0003-3826-5648

Avery E. Broderick https://orcid.org/0000-0002-3351-760X

Mohamad Shalaby https://orcid.org/0000-0001-9625-5929

Christoph Pfrommer https://orcid.org/0000-0002-7275-3998

Ewald Puchwein https://orcid.org/0000-0001-8778-7587

Philip Chang https://orcid.org/0000-0002-2137-2837

Astrid Lamberts https://orcid.org/0000-0001-8740-0127

References

Acero, F., Ackermann, M., Ajello, M., et al. 2015, ApJS, 218, 23

Ackermann, M., Ajello, M., Albert, A., et al. 2012b, ApJS, 203, 4

Ackermann, M., Ajello, M., Allafort, A., et al. 2012a, Sci, 338, 1190

Ackermann, M., Ajello, M., Allafort, A., et al. 2013, ApJ, 765, 54

Ackermann, M., Ajello, M., Baldini, L., et al. 2018, ApJS, 237, 32

Ando, S., & Kusenko, A. 2010, ApJL, 722, L39

Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214

Beck, R. 2016, A&ARv, 24, 4

Biermann, L. 1950, ZNatA, 5, 65

Blandford, R., & Eichler, D. 1987, PhR, 154, 1

Broderick, A. E., Chang, P., & Pfrommer, C. 2012, ApJ, 752, 22

Broderick, A. E., Pfrommer, C., Puchwein, E., & Chang, P. 2014a, ApJ, 790, 137

Broderick, A. E., Pfrommer, C., Puchwein, E., Chang, P., & Smith, K. M. 2014b, ApJ, 796, 12

Broderick, A. E., Tiede, P., Chang, P., et al. 2018, ApJ, 868, 87

Broderick, A. E., Tiede, P., Shalaby, M., et al. 2016, ApJ, 832, 109

Chang, P., Broderick, A. E., & Pfrommer, C. 2012, ApJ, 752, 23

Chang, P., Broderick, A. E., Pfrommer, C., et al. 2014, ApJ, 797, 110

Chang, P., Broderick, A. E., Pfrommer, C., et al. 2016, ApJ, 833, 118

Chen, W., Buckley, J. H., & Ferrer, F. 2015, PhRvL, 115, 211103

Crutcher, R. M. 2012, ARA&A, 50, 29

Dermer, C. D., Cavadini, M., Razzaque, S., et al. 2011, ApJL, 733, L21

Dolag, K., Grasso, D., Springel, V., & Tkachev, I. 2005, JCAP, 1, 009

Domínguez, A., Finke, J. D., Prada, F., et al. 2013, ApJ, 770, 77

Domínguez, A., Primack, J. R., Rosario, D. J., et al. 2011, MNRAS, 410, 2556

Gilmore, R. C., Somerville, R. S., Primack, J. R., & Domínguez, A. 2012, MNRAS, 422, 3189

Goldreich, P., & Sridhar, S. 1995, ApJ, 438, 763

Gould, R. J., & Schréder, G. P. 1967, PhRv, 155, 1408

Hanasz, M., & Lesch, H. 2003, A&A, 412, 331

Helfand, D. J., White, R. L., & Becker, R. H. 2015, ApJ, 801, 26

Hennebelle, P., & Falgarone, E. 2012, A&ARv, 20, 55

Hunter, J. D. 2007, CSE, 9, 90

Ichiki, K., Takahashi, K., Ohno, H., Hanayama, H., & Sugiyama, N. 2006, Sci, 311, 827

Jones, E., Oliphant, T., Peterson, P., et al. 2001, SciPy: Open source scientific tools for Python, http://www.scipy.org/

Kempf, A., Kilian, P., & Spanier, F. 2016, A&A, 585, A132

Kulsrud, R. M., & Zweibel, E. G. 2008, RPPh, 71, 046901

Lamberts, A., Chang, P., Pfrommer, C., et al. 2015, ApJ, 811, 19

Long, A. J., & Vachaspati, T. 2015, JCAP, 9, 065

Miniati, F., & Elyiv, A. 2013, ApJ, 770, 54

Neronov, A., Semikoz, D., Kachelriess, M., Ostapchenko, S., & Elyiv, A. 2010, ApJL, 719, L130

Neronov, A., & Semikoz, D. V. 2009, PhRvD, 80, 123012

Neronov, A., Semikoz, D. V., Tinyakov, P. G., & Tkachev, I. I. 2011, A&A, 526, 90

Neronov, A., & Vovk, I. 2010, Sci, 328, 73

Pakmor, R., Gomez, F. A., Grand, R. J. J., et al. 2017, MNRAS, 469, 3185

Pakmor, R., Pfrommer, C., Simpson, C. M., & Springel, V. 2016, ApJL, 824, L30

Pfrommer, C., Chang, P., & Broderick, A. E. 2012, ApJ, 752, 24

Price-Whelan, A. M., Sipőcz, B. M., Günther, H. M., et al. 2018, AJ, 156, 123

Puchwein, E., Pfrommer, C., Springel, V., Broderick, A. E., & Chang, P. 2012, MNRAS, 423, 149

Ryu, D., Schleicher, D. R. G., Treumann, R. A., Tsagas, C. G., & Widrow, L. M. 2012, SSRv, 166, 1

Salamon, M. H., & Stecker, F. W. 1998, ApJ, 493, 547

Schekochihin, A. A., Cowley, S. C., Dorland, W., et al. 2009, ApJS, 182, 310

Schlickeiser, R., Ibscher, D., & Supsar, M. 2012, ApJ, 758, 102

Schober, J., Schleicher, D. R. G., & Klessen, R. S. 2013, A&A, 560, A87

Shalaby, M., Broderick, A. E., Chang, P., et al. 2017, ApJ, 848, 81

Sironi, L., & Giannios, D. 2014, ApJ, 787, 49

Taylor, A. M., Vovk, I., & Neronov, A. 2011, A&A, 529, A144

Tiede, P., Broderick, A. E., Shalaby, M., et al. 2016, ApJ, 850, 157

Tsagas, C. G., & Kandus, A. 2005, PhRvD, 71, 123506

Turner, M. S., & Widrow, L. M. 1988, PhRvD, 37, 2743

Vafin, S., Deka, P. J., Pohl, M., & Bohdan, A. 2019, ApJ, 873, 10

Vafin, S., Rafighi, I., Pohl, M., & Niemiec, J. 2018, ApJ, 857, 43

Vovk, I., Taylor, A. M., Semikoz, D., & Neronov, A. 2012, ApJL, 747, L14

Widrow, L. M., Ryu, D., Schleicher, D. R. G., et al. 2012, SSRv, 166, 37