# Constraints on the intergalactic magnetic field using Fermi-LAT and H.E.S.S. blazar observations

F Aharonian, J Aschersleben, M Backes, V. Barbosa Martins, R Batzofin, Y Becherini, D Berge, B Bi, M Bouyahiaoui, M Breuhaus, et al.

► To cite this version:

F Aharonian, J Aschersleben, M Backes, V. Barbosa Martins, R Batzofin, et al. Constraints on the intergalactic magnetic field using Fermi-LAT and H.E.S.S. blazar observations. The Astrophysical journal letters, 2023, 950 (2), pp.L16. 10.3847/2041-8213/acd777. hal-04137311

HAL Id: hal-04137311

https://hal.science/hal-04137311

Submitted on 23 Jun 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
# Constraints on the Intergalactic Magnetic Field Using Fermi-LAT and H.E.S.S. Blazar Observations

F. Aharonian<sup>1,2</sup>①, J. Aschersleben<sup>3</sup>, M. Backes<sup>4,5</sup>①, V. Barbosa Martins<sup>6</sup>①, R. Batzofin<sup>7</sup>①, Y. Becherini<sup>8,9</sup>①, D. Berge<sup>6,10</sup>①, B. Bi<sup>11</sup>, M. Bouyahiaoui<sup>2</sup>, M. Breuhaus<sup>2</sup>①, R. Brose<sup>1</sup>①, F. Brun<sup>12</sup>①, B. Bruno<sup>13</sup>, T. Bulik<sup>14</sup>①, C. Burger-Scheidlin<sup>1</sup>, T. Bylund<sup>9</sup>①, S. Caroff<sup>15</sup>①, S. Casanova<sup>16</sup>①, J. Celic<sup>13</sup>, M. Cerruti<sup>8</sup>①, T. Chand<sup>5</sup>, S. Chandra<sup>5</sup>, A. Chen<sup>17</sup>①, J. Chibueze<sup>5</sup>①, O. Chibueze<sup>5</sup>, G. Cotter<sup>18</sup>①, M. de Bony<sup>15</sup>, K. Egberts<sup>7</sup>, J.-P. Ernenwein<sup>19</sup>, G. Fichet de Clairfontaine<sup>20</sup>①, M. Filipovic<sup>21</sup>①, G. Fontaine<sup>22</sup>①, M. Füssling<sup>6</sup>, S. Funk<sup>13</sup>①, S. Gabici<sup>8</sup>, S. Ghafourizadeh<sup>23</sup>, G. Giavitto<sup>6</sup>①, D. Glawion<sup>13</sup>①, J. F. Glicenstein<sup>12</sup>①, P. Goswami<sup>5</sup>, M.-H. Grondin<sup>24</sup>, L. Haerer<sup>2</sup>, T. L. Holch<sup>6</sup>①, M. Holler<sup>25</sup>, D. Horns<sup>26</sup>, M. Jamrozy<sup>27</sup>①, F. Jankowsky<sup>23</sup>, V. Joshi<sup>13</sup>①, I. Jung-Richardt<sup>13</sup>, E. Kasai<sup>4</sup>, K. Katarzyński<sup>28</sup>, R. Khatoon<sup>5</sup>, B. Khélifi<sup>8</sup>①, W. Kluźniak<sup>29</sup>, Nu. Komin<sup>17</sup>①, K. Kosack<sup>12</sup>①, D. Kostunin<sup>6</sup>①, R. G. Lang<sup>13</sup>, S. Le Stum<sup>19</sup>, F. Leitl<sup>13</sup>, A. Lemière<sup>8</sup>, J.-P. Lenain<sup>30</sup>①, F. Leuschner<sup>11</sup>①, T. Lohse<sup>10</sup>, A. Luashvili<sup>20</sup>①, I. Lypova<sup>23</sup>, J. Mackey<sup>1</sup>①, D. Malyshev<sup>11</sup>①, D. Malyshev<sup>13</sup>①, V. Marandon<sup>2</sup>①, P. Marchegiani<sup>17</sup>①, A. Marcowith<sup>31</sup>, G. Martí-Devesa<sup>25</sup>①, R. Marx<sup>23</sup>①, M. Meyer<sup>26,39</sup>①, A. Mitchell<sup>13</sup>①, R. Moderski<sup>29</sup>, L. Mohrmann<sup>2</sup>①, A. Montanari<sup>12</sup>①, E. Moulin<sup>12</sup>①, J. Muller<sup>22</sup>①, T. Murach<sup>6</sup>①, K. Nakashima<sup>13</sup>, J. Niemiec<sup>16</sup>①, S. Ohm<sup>6</sup>①, L. Olivera-Nieto<sup>2</sup>①, E. de Ona Wilhelmi<sup>6</sup>, S. Panny<sup>25</sup>①, M. Panter<sup>2</sup>, R. D. Parsons<sup>10</sup>①, G. Peron<sup>8</sup>, D. A. Prokhorov<sup>32</sup>, H. Prokoph<sup>6</sup>, G. Pühlhofer<sup>11</sup>①, M. Punch<sup>8</sup>①, A. Quirrenbach<sup>23</sup>, P. Reichherzer<sup>12</sup>①, A. Reimer<sup>25</sup>①, O. Reimer<sup>25</sup>, B. Reville<sup>2</sup>①, F. Rieger<sup>2</sup>, G. Rowell<sup>33</sup>①, B. Rudak<sup>29</sup>①, E. Ruiz-Velasco<sup>2</sup>①, V. Sahakian<sup>34</sup>①, D. A. Sanchez<sup>15</sup>, M. Sasaki<sup>13</sup>①, F. Schüssler<sup>12</sup>①, H. M. Schutte<sup>5</sup>①, U. Schwanke<sup>10</sup>, J. N. S. Shapopi<sup>4</sup>①, H. Sol<sup>20</sup>, S. Spencer<sup>13</sup>①, S. Steinmassl<sup>2</sup>①, H. Suzuki<sup>35</sup>, T. Takahashi<sup>36</sup>, T. Tanaka<sup>35</sup>①, A. M. Taylor<sup>6</sup>①, R. Terrier<sup>8</sup>①, C. Thorpe-Morgan<sup>11</sup>, M. Tsirou<sup>6</sup>, N. Tsuji<sup>37</sup>①, Y. Uchiyama<sup>38</sup>, C. van Eldik<sup>13</sup>①, J. Veh<sup>13</sup>①, C. Venter<sup>5</sup>, S. J. Wagner<sup>23</sup>①, R. White<sup>2</sup>, A. Wierzcholska<sup>16</sup>①, Yu Wun Wong<sup>13</sup>, M. Zacharias<sup>5,23</sup>①, D. Zargaryan<sup>1</sup>①, A. A. Zdziarski<sup>29</sup>①, S. Zouari<sup>8</sup>①, N. Żywucka<sup>5</sup>①

H.E.S.S. Collaboration,

and

M. Meyer<sup>26,39</sup>①

Fermi-LAT Collaboration

<footnote n="1"><sup>1</sup> Dublin Institute for Advanced Studies, 31 Fitzwilliam Place, Dublin 2, Ireland</footnote>

<footnote n="2"><sup>2</sup> Max-Planck-Institut für Kernphysik, P.O. Box 103980, D-69029 Heidelberg, Germany</footnote>

<footnote n="3"><sup>3</sup> Kapteyn Astronomical Institute, University of Groningen, Landleven 12, 9747 AD Groningen, The Netherlands</footnote>

<footnote n="4"><sup>4</sup> University of Namibia, Department of Physics, Private Bag 13301, Windhoek 10005, Namibia</footnote>

<footnote n="5"><sup>5</sup> Centre for Space Research, North-West University, Potchefstroom 2520, South Africa</footnote>

<footnote n="6"><sup>6</sup> DESY, D-15738 Zeuthen, Germany</footnote>

<footnote n="7"><sup>7</sup> Institut für Physik und Astronomie, Universität Potsdam, Karl-Liebknecht-Strasse 24/25, D-14476 Potsdam, Germany</footnote>

<footnote n="8"><sup>8</sup> Université de Paris, CNRS, Astroparticule et Cosmologie, F-75013 Paris, France</footnote>

<footnote n="9"><sup>9</sup> Department of Physics and Electrical Engineering, Linnaeus University, SE-351 95 Växjö, Sweden; contact.hess@hess-experiment.eu</footnote>

<footnote n="10"><sup>10</sup> Institut für Physik, Humboldt-Universität zu Berlin, Newtonstr. 15, D-12489 Berlin, Germany</footnote>

<footnote n="11"><sup>11</sup> Institut für Astronomie und Astrophysik, Universität Tübingen, Sand 1, D-72076 Tübingen, Germany</footnote>

<footnote n="12"><sup>12</sup> IRFU, CEA, Université Paris-Saclay, F-91191 Gif-sur-Yvette, France</footnote>

<footnote n="13"><sup>13</sup> Friedrich-Alexander-Universität Erlangen-Nürnberg, Erlangen Centre for Astroparticle Physics, Erwin-Rommel-Str. 1, D-91058 Erlangen, Germany
<sup>14</sup> Astronomical Observatory, The University of Warsaw, Al. Ujazdowskie 4, 00-478 Warsaw, Poland</footnote>

<footnote n="15"><sup>15</sup> Université Savoie Mont Blanc, CNRS, Laboratoire d'Annecy de Physique des Particules—IN2P3, F-74000 Annecy, France</footnote>

<footnote n="16"><sup>16</sup> Instytut Fizyki Jądrowej PAN, ul. Radzikowskiego 152, 31-342 Kraków, Poland</footnote>

<footnote n="17"><sup>17</sup> School of Physics, University of the Witwatersrand, 1 Jan Smuts Avenue, Braamfontein, Johannesburg, 2050 South Africa</footnote>

<footnote n="18"><sup>18</sup> University of Oxford, Department of Physics, Denys Wilkinson Building, Keble Road, Oxford OX1 3RH, UK</footnote>

<footnote n="19"><sup>19</sup> Aix Marseille Université, CNRS/IN2P3, CPPM, Marseille, France</footnote>

<footnote n="20"><sup>20</sup> Laboratoire Univers et Théories, Observatoire de Paris, Université PSL, CNRS, Université de Paris, F-92190 Meudon, France</footnote>

<footnote n="21"><sup>21</sup> School of Science, Western Sydney University, Locked Bag 1797, Penrith South DC, NSW 2751, Australia</footnote>

<footnote n="22"><sup>22</sup> Laboratoire Leprince-Ringuet, École Polytechnique, CNRS, Institut Polytechnique de Paris, F-91128 Palaiseau, France</footnote>

<footnote n="23"><sup>23</sup> Landessternwarte, Universität Heidelberg, Königstuhl, D-69117 Heidelberg, Germany</footnote>

<footnote n="24"><sup>24</sup> Université Bordeaux, CNRS, LP2I Bordeaux, UMR 5797, F-33170 Gradignan, France</footnote>

<footnote n="25"><sup>25</sup> Institut für Astro- und Teilchenphysik, Leopold-Franzens-Universität Innsbruck, A-6020 Innsbruck, Austria</footnote>

<footnote n="26"><sup>26</sup> Universität Hamburg, Institut für Experimentalphysik, Luruper Chaussee 149, D-22761 Hamburg, Germany</footnote>

<footnote n="27"><sup>27</sup> Obserwatorium Astronomiczne, Uniwersytet Jagielloński, ul. Orla 171, 30-244 Kraków, Poland</footnote>

<footnote n="28"><sup>28</sup> Institute of Astronomy, Faculty of Physics, Astronomy and Informatics, Nicolaus Copernicus University, Grudziadzka 5, 87-100 Torun, Poland</footnote>

<footnote n="29"><sup>29</sup> Nicolaus Copernicus Astronomical Center, Polish Academy of Sciences, ul. Bartycka 18, 00-716 Warsaw, Poland</footnote>

<footnote n="30"><sup>30</sup> Sorbonne Université, Université Paris Diderot, Sorbonne Paris Cité, CNRS/IN2P3, Laboratoire de Physique Nucléaire et de Hautes Energies, LPNHE, 4 Place Jussieu, F-75252 Paris, France</footnote>

<footnote n="31"><sup>31</sup> Laboratoire Univers et Particules de Montpellier, Université Montpellier, CNRS/IN2P3, CC 72, Place Eugène Bataillon, F-34095 Montpellier Cedex 5, France</footnote>

<footnote n="32"><sup>32</sup> GRAPPA, Anton Pannekoek Institute for Astronomy, University of Amsterdam, Science Park 904, 1098 XH Amsterdam, The Netherlands</footnote>

<footnote n="33"><sup>33</sup> School of Physical Sciences, University of Adelaide, Adelaide 5005, Australia</footnote>

<footnote n="34"><sup>34</sup> Yerevan Physics Institute, 2 Alikhanian Brothers St., 375036 Yerevan, Armenia</footnote>

<footnote n="35"><sup>35</sup> Department of Physics, Konan University, 8-9-1 Okamoto, Higashinada, Kobe, Hyogo 658-8501, Japan</footnote>

<footnote n="36"><sup>36</sup> Kavli Institute for the Physics and Mathematics of the Universe (WPI), The University of Tokyo Institutes for Advanced Study (UTIAS), The University of Tokyo, 5-1-5 Kashiwa-no-Ha, Kashiwa, Chiba 277-8583, Japan</footnote>
<sup>37</sup>RIKEN, 2-1 Hirosawa, Wako, Saitama 351-0198, Japan

<sup>38</sup>Department of Physics, Rikkyo University, 3-34-1 Nishi-Ikebukuro, Toshima-ku, Tokyo 171-8501, Japan
*Received 2023 February 27; revised 2023 May 8; accepted 2023 May 22; published 2023 June 21*

**Abstract**

Magnetic fields in galaxies and galaxy clusters are believed to be the result of the amplification of intergalactic seed fields during the formation of large-scale structures in the universe. However, the origin, strength, and morphology of this intergalactic magnetic field (IGMF) remain unknown. Lower limits on (or indirect detection of) the IGMF can be obtained from observations of high-energy gamma rays from distant blazars. Gamma rays interact with the extragalactic background light to produce electron–positron pairs, which can subsequently initiate electromagnetic cascades. The gamma-ray signature of the cascade depends on the IGMF since it deflects the pairs. Here we report on a new search for this cascade emission using a combined data set from the Fermi Large Area Telescope and the High Energy Stereoscopic System. Using state-of-the-art Monte Carlo predictions for the cascade signal, our results place a lower limit on the IGMF of $B > 7.1 \times 10^{-16}$ G for a coherence length of 1 Mpc even when blazar duty cycles as short as 10 yr are assumed. This improves on previous lower limits by a factor of 2. For longer duty cycles of $10^4$ (10<sup>7</sup>) yr, IGMF strengths below $1.8 \times 10^{-14}$ G (3.9 × 10<sup>-14</sup> G) are excluded, which rules out specific models for IGMF generation in the early universe.

*Unified Astronomy Thesaurus concepts: Active galaxies (17); High energy astrophysics (739); Extragalactic magnetic fields (507)*

**1. Introduction**

Astrophysical magnetic fields in galaxies and galaxy clusters are ubiquitously observed and believed to have been seeded by a turbulent intergalactic magnetic field (IGMF), which should still be present at its seed field value and coherence length in cosmic voids (see, e.g., Durrer & Neronov 2013; Alves Batista & Saveliev 2021, for reviews). Such a seed field should have a small field strength in the voids, making it extremely difficult to observe directly. Currently, only upper limits of the order of $10^{-9}$ G are known from Faraday rotation measurements (Pshirkov et al. 2016).

Detecting or constraining the IGMF would provide insight into many processes in the early universe; it might be possible to determine whether the IGMF was created during inflation or during phase transitions and to probe contributions of beyond-the-standard-model physics. Alternatively, the IGMF could be found to have an astrophysical origin.

Observations of blazars (active galactic nuclei (AGNs) with relativistic bulk outflows closely aligned with the line of sight (LOS)) at gamma-ray energies offer a unique alternative probe of the IGMF, with different systematics from other established methods (e.g., Durrer & Neronov 2013) During propagation, gamma rays interact with photons of background radiation fields, most notably the extragalactic background light (EBL), and produce electron–positron pairs. These pairs can inverse Compton (IC) scatter cosmic microwave background (CMB) photons up to gamma-ray energies, thereby initiating a cascade as these upscattered photons again pair-produce. The deflection of the pairs in the IGMF delays the arrival of the gamma rays produced in the cascades (compared to the primary gamma rays; Plaga 1995), which should appear as an extended halo around the blazar (Aharonian et al. 1994) and enhance its apparent brightness at lower gamma-ray energies (Neronov & Semikoz 2009).

The higher the value of the IGMF, the stronger the deflection of the pairs and, in turn, the larger the pair halos (with lower surface brightness), the lower the GeV excess of cascade photons, and the longer the time delays. The nonobservation of an excess of GeV photons in individual blazar observations with the Fermi Large Area Telescope (LAT) has thus been used to place lower limits on the IGMF of the order of $B \gtrsim 10^{-15}$ G (e.g., Neronov & Vovk 2010). Similarly, searches for a halo with imaging air Cherenkov telescopes (IACTs) have ruled out an IGMF strength between $\sim 10^{-16.5}$ and $\sim 10^{-15.5}$ G (H.E.S.S. Collaboration et al. 2014), as well as between $\sim 10^{-14.5}$ and $\sim 10^{-13.5}$ G (Archambault et al. 2017). The gamma-ray constraints are considerably relaxed when short activity times of blazars are considered (Dermer et al. 2011).<sup>40</sup> However, even if blazars have only been active for the past 10 yr, a joint analysis of Fermi-LAT and archival IACT spectra from several sources is still able to exclude $B \lesssim 3 \times 10^{-16}$ G for a coherence length of $l_B \gtrsim 10$ kpc (Ackermann et al. 2018). A comprehensive overview of current constraints is provided by Alves Batista & Saveliev (2021).

Here we perform an updated search for the cascade emission of hard-spectrum blazars. In contrast to previous cascade searches, we use a combination of Fermi-LAT data and observations carried out with the High Energy Stereoscopic System (H.E.S.S.) based on the full Poisson likelihood for gamma-ray event statistics and calculate pair halo templates with state-of-the-art Monte Carlo simulations.

**2. Observations**

**2.1. Source Selection**

Among the blazars observed with both H.E.S.S. and the LAT, we aim to select those for which the cascade emission is expected to be maximal. Extreme high-frequency-peaked BL Lac–type objects (EHBLs) with an observed synchrotron peak $\gtrsim 1$ keV are the most promising for such searches. Their hard spectra are measured up to TeV energies without finding evidence for a

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

<footnote n="39"><sup>39</sup>Now at CP3-Origins, University of Southern Denmark, Campusvej 55, 5230 Odense M, Denmark.</footnote>

<footnote>Original content from this work may be used under the terms of the [Creative Commons Attribution 4.0 licence](https://creativecommons.org/licenses/by/4.0/). Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.</footnote>

<footnote n="40"><sup>40</sup>We define the activity time (duty cycle) as the length of the time during which the blazar has been emitting gamma rays. Throughout this work, we assume that a blazar has emitted a flux that is constant in time, which can be described with the average spectrum observed today.</footnote>
<table><caption>**Table 1**
Selected Sources and Observation Summaries</caption><thead><tr><th rowspan="2">Source</th><th rowspan="2">z</th><th colspan="3">H.E.S.S.</th><th colspan="2">Fermi-LAT</th></tr><tr><th>Live Time (hr)</th><th>Significance (σ)</th><th>Γ<sub>int</sub></th><th>√TS</th><th>Γ<sub>int</sub></th></tr></thead><tbody><tr><td>1ES 0229+200</td><td>0.140</td><td>144.1</td><td>16.5</td><td>1.76 ± 0.12</td><td>14.2</td><td>1.73 ± 0.12</td></tr><tr><td>1ES 0347-121</td><td>0.188</td><td>59.2</td><td>16.1</td><td>2.12 ± 0.15</td><td>18.3</td><td>1.81 ± 0.11</td></tr><tr><td>PKS 0548-322</td><td>0.069</td><td>53.9</td><td>10.2</td><td>1.92 ± 0.12</td><td>11.7</td><td>2.06 ± 0.17</td></tr><tr><td>1ES 1101-232</td><td>0.186</td><td>71.9</td><td>18.7</td><td>1.66 ± 0.09</td><td>19.5</td><td>1.65 ± 0.09</td></tr><tr><td>H2356-309</td><td>0.165</td><td>150.5</td><td>23.4</td><td>2.10 ± 0.09</td><td>21.7</td><td>1.81 ± 0.09</td></tr></tbody></table>

**Note.** Second column: spectroscopic redshift, *z*, taken from (in order of increasing R.A.) Fosbury & Disney (1976), Remillard et al. (1989), Woo et al. (2005), Jones et al. (2009), and Peña-Herazo et al. (2021). Third column: total observation time taking data quality cuts and dead time into account. Fourth and sixth columns: detection significance in units of σ. Fifth and seventh columns: power-law indices for EBL-corrected spectra.

source-intrinsic spectral cutoff. Therefore, a lot of energy is potentially available for reprocessing into the cascade. To select EHBL sources, we start from the fourth catalog of AGNs observed with the LAT (4LAC; Lott et al. 2020). We select blazars that (i) follow a power-law spectrum in energy $E$, $dN/dE = N_0(E/E_0)^{-\Gamma}$, with spectral index $\Gamma < 2$ within statistical uncertainties; (ii) have a known spectroscopic redshift $z$; and (iii) have a synchrotron peak $\nu_{\text{peak}} > 10^{17}$ Hz. Additionally, we demand that the flux have a chance probability of being variable at less than 99% in the LAT energy band on monthly timescales. This choice excludes sources that exhibit strong variability. Otherwise, the past temporal flux evolution would have to be taken into account to predict the time-delayed cascade emission. This selection results in five sources that have been previously detected with H.E.S.S. (Aharonian et al. 2006b, 2007a, 2007b, 2010). These are listed in Table 1. Except for PKS 0548-322, all of the selected sources are also included in the analysis of Ackermann et al. (2018). In addition, Ackermann et al. (2018) considered the H.E.S.S.-detected blazar 1ES 0414+009; however, this source fails our cuts on the synchrotron peak frequency and the probability for variability.

## 2.2. Fermi-LAT Observations

The Fermi-LAT is a pair-conversion telescope that detects gamma rays from 20 MeV to beyond 300 GeV (Atwood et al. 2009). For our five sources, we analyze 11.5 yr of LAT data between 1 GeV and 3 TeV extracted from a 6° × 6° region of interest (ROI) centered on the nominal source position. The data selection closely follows Ackermann et al. (2018). Using events and instrumental response functions of the P8R3_SOURCE_V3 class, we find the optimal model for the ROI that includes all point sources listed in the second data release of the fourth LAT point-source catalog (4FGL-DR2; Ballet et al. 2020), as well as templates for the diffuse Galactic and isotropic backgrounds (see Appendix A for further details). All of our blazars are modeled with power laws in our initial ROI model. After an initial optimization of the model, we change the spectral model to a power law with index $\Gamma_{\text{int}}$ that includes absorption on the EBL, given by $\exp(-\tau_{\gamma\gamma})$, where $\tau_{\gamma\gamma}$ is the optical depth. Throughout this article, the EBL absorption is taken from the model of Domínguez et al. (2011). The obtained spectral indices and source significances, expressed as the square root of the test statistic (TS),$^{41}$ are listed in Table 1. In Figure 1, we show a sky map with TS values for our optimized ROI model for

<figure><img src="image_3.png" /><figcaption>**Figure 1.** A TS map for the LAT analysis of the ROI around 1ES 0229+200. Point sources within the ROI are marked with plus signs. The 68% containment radii at 1 GeV for the LAT point-spread function for the events with the worst and best angular reconstruction are shown as a dashed and solid circle, respectively. The thin white contour indicates a region where $\sqrt{\text{TS}} = 3$.</figcaption></figure>

1ES 0229+200. In each pixel, the significance of a putative source with a power-law index equal to 2 is calculated. A potential halo could reveal itself as a TS excess around the central source. This is, however, not observed for any of the studied sources (see also Appendix A).

## 2.3. H.E.S.S. Observations

H.E.S.S. is an array of five IACTs located in the Khomas highlands in Namibia. The array is composed of four telescopes with an effective mirror diameter of 12 m, called CT1-CT4, and one with an effective mirror diameter of 28 m, called CT5. New cameras were installed between 2015 and 2016 for CT1-CT4 and in 2019 for CT5. H.E.S.S. can detect very high energy gamma rays between ≈50 GeV and ≈100 TeV. All sources considered here have been detected with H.E.S.S. For all EHBLs except PKS 0548-322, additional data have been taken since the original publications during different phases of the experiment (with and without CT5, before and after the camera exchange campaigns). In order to have a data set with a hardware setup as homogeneous as possible, we restrict ourselves to runs$^{42}$ with

<footnote>$^{41}$ The test statistic is defined as $TS = -2(\ln L_1 - \ln L_0)$, where $L_1$ ($L_0$) is the likelihood of the model including (excluding) a source.</footnote>

<footnote>$^{42}$ One run is a consecutive observation of a source with a duration of approximately 28 minutes.</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Simulated halo photon flux per solid angle, $dN/d\Omega$, integrated above 1 GeV for the best-fit spectrum (see Section 3.2) for 1ES 0229+200. An activity time of $t_{max} = 10^4$ yr and a field strength of $B = 3.2 \times 10^{-16}$ G are assumed. The same 68% containment radius as in Figure 1 is shown for comparison, along with the On region for the H.E.S.S. observations (dotted-dashed circle). The hot spots correspond to individual simulated particles.</figcaption></figure>

good data quality in which three or four of the smaller telescopes participated. For the same reason, we only use data taken with the cameras originally mounted on the 12 m telescopes. The energies of the recorded gamma-ray-like events are reconstructed using the ImPACT method and standard selection cuts (Parsons & Hinton 2014). The spectral analysis is performed using gammapy version 0.18.2 (Deil et al. 2017). Source counts are extracted from an “On” region of radius of 0.07 centered on the source location, which corresponds to the 68% confidence radius of the point-spread function. The observed spectra are modeled with power laws with EBL attenuation; models with spectral curvature are not preferred. The resulting spectral indices are provided in Table 1, together with observation live times and detection significances. Further details on the analysis are provided in Appendix B.

# 3. Search for the Pair Halo and Constraints on the IGMF

## 3.1. Cascade Predictions

We generate templates for the pair halo for each redshift using the CRPropa Monte Carlo simulation package version 3.1.6$^{43}$ (Alves Batista et al. 2016a). In the simulation, we include all relevant particle interactions and consider cascade photons arriving with a time delay shorter than a maximum blazar activity time $t_{max}$. Using a spectral reweighting (Ackermann et al. 2018) and parallel transport (Alves Batista et al. 2016b), the simulations are converted into energy-dependent sky maps for arbitrary input spectra and angles $\theta_{obs}$ between the jet axis and the LOS (see Appendix C for further details on the simulation).$^{44}$ These serve as input templates for our data analysis. For simplicity, $\theta_{obs} = 0^\circ$ is assumed throughout. For $0 < \theta_{obs} < \theta_{jet}/2$, the templates will be asymmetric. However, in the considered energy range, the asymmetry is small, and our results do not depend strongly on

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Observed LAT and H.E.S.S. spectra of 1ES 0229+200 (blue and red markers and lines, respectively), together with the total best-fit spectrum (orange solid band, including uncertainties) for a point-source model of the primary gamma-ray flux (dashed orange line) and the halo flux (dashed-dotted orange line). Spectral points are shown as upper limits if the detection significance is below $2\sigma$. The same values for $t_{max}$ and $B$ are assumed as in Figure 2.</figcaption></figure>

$\theta_{obs}$ and the orientation of the halo. The IGMF is modeled as cell-like; each cell has a side length of $l_B = 1$ Mpc, and the B-field orientation changes randomly from one cell to the next. The templates are generated for each source and for seven values of the field strength, $B = 10^{-16}, 10^{-15.5}, \dots, 10^{-13}$ G. For higher values of $B$, the pairs are quickly isotropized and the cascade emission would appear as an additional component to the isotropic gamma-ray background in the LAT energy band. An example of the simulated energy-integrated sky map is shown in Figure 2, and an output spectrum is shown in Figure 3. Spectra for all other sources and different choices of the IGMF are presented in Figure 6 in Appendix D.

## 3.2. Combined LAT and H.E.S.S. Analysis

We test for the presence of a halo for a given field strength using a likelihood ratio test. Different fixed values of $t_{max} = 10, 10^4$, and $10^7$ yr are assumed, which correspond to a conservative, intermediate, and optimistic scenario, respectively, in terms of the amount of cascade emission. Our model is the sum of the point-source flux $\phi_{src}(E, \theta_i)$ for central source $i$ and spectral parameters $\theta_i$ and the diffuse halo flux $\phi_{halo}(E, p, \theta_i, B)$, which also depends on the sky position $p$ and the IGMF. In contrast to Section 2, the point-source spectra are now modeled with a power law with a cutoff at energy $E_{cut}$,

$$ \phi_{src} = N_0(E/E_0)^{-\Gamma_{int}} \exp(-E/E_{cut} - \tau_{\gamma\gamma}), \quad (1) $$

i.e., $\theta = (N_0, \Gamma_{int}, E_{cut}, E_0)$. With this choice, we obtain plausible constraints on the IGMF since the high-energy cutoff reduces the number of gamma rays that can initiate the cascade. Assuming instead a simple power law as the intrinsic spectrum would bear the risk of overestimating the flux of primary gamma rays and thus of the cascade.

The likelihood of observing the data $D_i$ given the model parameters is the product of the individual Poisson likelihoods for LAT and H.E.S.S. data, $\mathcal{L}_i \equiv \mathcal{L}(B, \theta_i|D_i) = \mathcal{L}(B, \theta_i|D_{LAT,i}) \mathcal{L}(B, \theta_i|D_{H.E.S.S.,i})$. The log-likelihood ratio test

<footnote>$^{43}$ https://github.com/CRPropa/CRPropa3</footnote>

<footnote>$^{44}$ A Python package performing these calculations is publicly available at https://github.com/me-manu/simCRPropa.</footnote>
can then be written as

$$ \lambda(B) = -2 \sum_{i} \ln \left( \frac{\mathcal{L}(B, \hat{\theta}_{i} | \mathcal{D}_{i})}{\mathcal{L}(\hat{B}, \hat{\theta}_{i} | \mathcal{D}_{i})} \right), \quad (2) $$

where $\hat{\theta}_i$ denotes the spectral parameters maximizing $\mathcal{L}_i$ for a fixed field strength, whereas $\hat{\theta}_i$ and $\hat{B}$ are the parameters that maximize $\mathcal{L}_i$ unconditionally. We first extract $\ln \mathcal{L}(B, \theta_i | \mathcal{D}_{\text{LAT},i})$ over a grid of IGMF values and spectral parameters by adding the reweighted cascade template to the ROI model (see Section 2.2) and reoptimizing the spectral parameters of the background sources (see Appendix A.2 for further details). In a second step, we determine the spectral parameters that maximize the total $\mathcal{L}_i$ for a given value of B using the MINUIT fit routine in gammapy. For the H.E.S.S. data, the model includes the point source and the pair halo intensity within the On region. In each step of the parameter optimization, the cascade template is reweighted to account for the changing spectrum, and the extracted value of $\ln \mathcal{L}(B, \theta_i | \mathcal{D}_{\text{LAT},i})$ is added to $\ln \mathcal{L}(B, \theta_i | \mathcal{D}_{\text{H.E.S.S.,i}})$. In this way, for a given set of IGMF parameters, the best-fit spectrum of each source for both the LAT and H.E.S.S. energy ranges is determined. An example of the fit for 1ES 0229+200 is shown in Figure 3. Figures for other sources are shown in Appendix D.1, where we also report the best-fit parameters.<fnref n="45" />

In general, high values of B provide a better fit to the data than low IGMF values, implying halos with large extensions and a low surface brightness. For low values of B, the fit prefers low values of $E_{\text{cut}}$ to minimize the cascade contribution. For instance, for $t_{\text{max}} \ge 10^4$ yr and $B = 10^{-16}$ G the cutoff energy is below $E_{\text{cut}} \lesssim 3$ TeV for all sources except PKS 0548–322. However, even the fits assuming high IGMF values are not significantly preferred over the fit without cascade emission. The largest value for the log-likelihood ratio between the fit without and with a halo is found to be TS ~ 0.9 for $B = 10^{-13}$ G and $t_{\text{max}} = 10^4$ yr. Thus, we do not find any evidence for the presence of a halo and proceed to set lower limits on the IGMF strength. As the field strength constitutes one additional degree of freedom, we set 95% confidence limits when $\lambda > 2.71$. We find that we can rule out an IGMF with a field strength $B < 3.9 \times 10^{-14}$ G for $t_{\text{max}} = 10^7$ yr. For a more conservative choice of the blazar activity time, we exclude $B < 1.8 \times 10^{-14}$ G for $t_{\text{max}} = 10^4$ yr and $B < 7.1 \times 10^{-16}$ G for $t_{\text{max}} = 10$ yr. The strongest constraints are provided by 1ES 0229+200 and 1ES 1101–232 (see Appendix D.2 for the full likelihood profiles). These limits, together with previous constraints and theoretically motivated parameters, are shown in Figure 4.

## 4. Discussion

The analysis presented here improves on previous analyses of a combination of IACT and LAT data in two ways. First, we use the CRPropa code to generate predictions for the spectral and spatial shape of the halo. As the code traces the full trajectories of the particles, we are not limited by the assumption of, e.g., small deflections of the pairs in the IGMF, which limited the validity range of the results of

Ackermann et al. (2018). Second, we combine LAT and H.E.S. S. data on the likelihood level and self-consistently include the halo contribution. As a result, we are able to improve on the constraints found by Ackermann et al. (2018) by a factor of two for blazar activity times as short as 10 yr. To our knowledge, these are the tightest constraints on the IGMF from gamma-ray observations to date (compare Figure 4).

The assumed activity time of the sources poses the largest systematic uncertainty in our study. Changing $t_{\text{max}}$ from 10 to $10^4$ yr strengthens the constraint by almost a factor of 30 (activity times of $10^7–10^8$ yr could be realistic; Parma et al. 2002). For $t_{\text{max}} = 10^7$, we cannot rule out as high IGMF values as Archambault et al. (2017). The likely reason is that we do not include the spatial halo extension in the H.E.S.S. energy range, as we only consider the cascade photons that arrive within the On region. We are instead able to probe low IGMF values thanks to the combination of LAT and H.E.S.S. data. We have deliberately selected sources that do not experience strong variability in the LAT data. However, at least H2356 –309 and 1ES 0229+200 show flux variations by factors ~2–3 above 100 GeV on monthly and yearly timescales, respectively (H.E.S.S. Collaboration et al. 2010; Acciari et al. 2023). Our reweighting framework developed in Appendix C could in principle be extended to study how long-term variability affects the predictions for the cascade. This is left for future study. Additional systematic uncertainties with a smaller effect on the results compared to the change in activity time, such as the choice of the EBL model, the choice of $\theta_{\text{jet}}$, and the different absolute energy scales of the LAT and H.E.S. S. (which can be treated self-consistently in our likelihood framework), are discussed in Appendix D.3.

Under the assumptions that (a) the sources have been active for 10 yr $\lesssim t_{\text{max}} \lesssim 10^4$ yr and (b) our constraints scale as $\ell_B^{-1/2}$ in the case when $\ell_B$ is shorter than the mean free path of the pairs for IC scattering (Neronov & Semikoz 2009; see, however, Caprini & Gabici 2015), our results challenge particular models for the IGMF generation. For example, we can rule out models where the seed fields are created through axion-driven magnetogenesis (which predicts a present-day IGMF with $B \sim 10^{-13}$ G and $\ell_B \sim 20$ pc; see Miniati et al. 2018) and the low end of allowed B-field values generated during or shortly after inflation with values $B \gtrsim 10^{-15}$ G on megaparsec scales (Kobayashi 2014).

It could be possible that the pairs predominantly lose their energy by heating the intergalactic medium through plasma instabilities, which would suppress the cascade (Broderick et al. 2012). Whether such instabilities develop is, however, still debated (see, e.g., Alves Batista & Saveliev 2021). We can introduce an additional normalization for our halo templates, $0 \le N_{\text{halo}} \le 1$, to mimic the possibility that the pairs lose part of their energy through this process ($N_{\text{halo}} = 1$ corresponds to the case discussed so far). The fit tends to select $N_{\text{halo}} \sim 0$, as the halo is not present in the data, i.e., no energy is dissipated through IC scattering. In this case, no limits on the IGMF are possible. We note, however, that our constraints for $t_{\text{max}} = 10$ yr are not affected, as the energy loss due to beam instabilities would only dominate over IC cooling after hundreds of years of blazar activity (Broderick et al. 2012).

In the future, observations with the upcoming Cherenkov Telescope Array (CTA) should be able to detect an IGMF as strong as $10^{-13}$ G for sufficiently long blazar activity times (see Figure 4) thanks to its improved angular resolution, point-

<footnote n="45">⁴⁵ We provide the Fermi-LAT and H.E.S.S. spectra, best-fit values, as well as the cascade simulation output electronically at 10.5281/zenodo.8014311.</footnote>
<figure><img src="image_6.png" /><figcaption>**Figure 4.** The IGMF parameter space probed with gamma-ray instruments. Lower limits at 95% confidence level on the IGMF derived from the combined Fermi-LAT and H.E.S.S. analysis for different assumed blazar activity times are shown as orange filled and hatched regions. A scaling of the limits with the coherence length as $\ell_B^{-1/2}$ is assumed. Exclusion regions from H.E.S.S., MAGIC, and VERITAS are shown as gray filled regions (Aleksić et al. 2010; H.E.S.S. Collaboration et al. 2014; Archambault et al. 2017). The CTA could exclude field strengths below the dashed gray line (Abdalla et al. 2021). The dark-gray line shows combined Fermi-LAT and IACT lower limits for $t_{max} = 10$ yr (Ackermann et al. 2018). Light-gray shaded regions are disfavored by theory. Blue lines and regions indicate theoretically favored regions for the generation of the IGMF during either electro-weak phase transitions (EWPT) or QCD phase transitions (QCDPT), by Galactic winds or induced by cosmic-ray (CR) streaming (see Durrer & Neronov 2013, and references therein).</figcaption></figure>

source sensitivity, and covered energy range (Meyer et al. 2016; Abdalla et al. 2021).

## Acknowledgments

The support of the Namibian authorities and of the University of Namibia in facilitating the construction and operation of H.E.S.S. is gratefully acknowledged, as is the support by the German Ministry for Education and Research (BMBF), the Max Planck Society, the German Research Foundation (DFG), the Helmholtz Association, the Alexander von Humboldt Foundation, the French Ministry of Higher Education, Research and Innovation, the Centre National de la Recherche Scientifique (CNRS/IN2P3 and CNRS/INSU), the Commissariat à l'énergie atomique et aux énergies alternatives (CEA), the U.K. Science and Technology Facilities Council (STFC), the Irish Research Council (IRC) and the Science Foundation Ireland (SFI), the Knut and Alice Wallenberg Foundation, the Polish Ministry of Education and Science, agreement No. 2021/WK/06, the South African Department of Science and Technology and National Research Foundation, the University of Namibia, the National Commission on Research, Science & Technology of Namibia (NCRST), the Austrian Federal Ministry of Education, Science and Research and the Austrian Science Fund (FWF), the Australian Research Council (ARC), the Japan Society for the Promotion of Science, the University of Amsterdam, and the Science Committee of Armenia grant 21AG-1C085. We appreciate the excellent work of the technical support staff in Berlin, Zeuthen, Heidelberg, Palaiseau, Paris, Saclay, and Tübingen, as

well as in Namibia in the construction and operation of the equipment. This work benefited from services provided by the H.E.S.S. Virtual Organisation, supported by the national resource providers of the EGI Federation.

The Fermi-LAT Collaboration acknowledges support for LAT development, operation, and data analysis from NASA and DOE (United States), CEA/Irfu and IN2P3/CNRS (France), ASI and INFN (Italy), MEXT, KEK, and JAXA (Japan), and the K. A. Wallenberg Foundation, the Swedish Research Council, and the National Space Board (Sweden). Science analysis support in the operations phase from INAF (Italy) and CNES (France) is also gratefully acknowledged. This work performed in part under DOE contract DE-AC02-76SF00515.

**Facilities:** Fermi-LAT, H.E.S.S.

**Software:** astropy (Astropy Collaboration et al. 2013), numpy (Harris et al. 2020), scipy (Virtanen et al. 2020), CRPropa3 (Alves Batista et al. 2016a), fermipy (Wood et al. 2017), fermitools (Fermi Science Support Development Team 2019), gammapy (Deil et al. 2017).

# Appendix A
## Details on the Fermi-LAT Data Analysis

### A.1. Data Selection and Initial Model Optimization

In this appendix, we provide further details on the Fermi-LAT data selection and analysis choices. Events of the P8R3_SOUR-CE_V3 class with arrival times between 2008 August 4 and 2020 January 4 are extracted from a 6° × 6° sky region in the energy
range from 1 GeV to 3 TeV that arrived with a zenith angle
≤100° in order to minimize contamination of our sample with
gamma rays from the Earth limb. The standard data filter
DATA QUAL > 0 && LAT_CONFIG == 1 is applied. Following
Ackermann et al. (2018), the events are then binned using a
spatial binning of 0.025 pixel⁻¹ and a logarithmic binning in
energy with 8 bins per decade. Furthermore, we split the data
sample into two subsamples using the quality of the angular
reconstruction (the so-called PSF event classes; see Atwood
et al. 2013). The first subsample contains all events belonging to
the PSF event classes 0–2 (event type 28), and the second
sample contains the events with the best reconstructed arrival
directions, PSF event class 3 (event type 32). A combined
likelihood analysis is performed on these two subsamples using
fermipy<fnref n="46" /> (Wood et al. 2017) version 0.19.0 and the
fermitools version 1.2.23⁴⁷ with the MINUIT optimizer.

For the initial optimization of the ROI model we include the
point sources listed in the 4FGL-DR2 up to an angular distance
of 10° from the ROI center. Templates for the galactic diffuse
(gll_iem_v07) and isotropic diffuse (iso_P8R3_SOUR-
CE_V3_v1) emission are added to the model as well.<fnref n="48" /> All
spectral parameters of sources up to a radius of 6° are left free
to vary, and only spectral normalizations are free parameters
for sources at larger distances. After an initial optimization
using the fermipy optimize routine, we fix all spectral
source parameters (spectral shape parameters, i.e., all para-
meters but the spectral normalization) for sources that are
detected with a TS < 5 (TS < 9) or that have a predicted
number of photons less than 10⁻². The flux normalization
parameter of sources with TS > 50 that was originally frozen is
also left free to vary during the fit. After a first fit to the data, a
TS map is computed. Putative point sources (modeled with a
power law) are added consecutively to the ROI to hot spots
with TS > 25. This procedure could tentatively mask any
additional unaccounted flux from the cascade halo. Only one
additional source had to be added in the vicinity of 1ES 1101
−232 at a distance of larger than 1°, which should not have an
effect on the halo search. From this it appears unlikely that our
analysis procedure could mask an additional halo contribution.
The TS maps for 1ES 0229+200 are shown in Figure 1 and for
all other sources in Figure 5. The above analysis choices result
in generally flat spectral and spatial residuals.

The best-fit spectral parameters for a point-source analysis
with and without EBL absorption are provided in Table 2. For
the fit with EBL, we only provide the spectral index, Γint,
which, for all sources, is compatible with the fit without
absorption (with index Γ) within 1σ uncertainties. The reason is
the limited number of observed gamma rays by the LAT in the
energy range where the optical depth τγγ ≳ 1. The obtained
point-source spectra for all sources are shown as blue squares
and solid lines in Figure 6 in Appendix D.1. We note that the
LAT spectrum of 1ES 0347-121 exhibits evidence for a
detection (√TS ≳ 4) of the source above 1 TeV. This result
is due to two photons detected above 1 TeV in the vicinity of
the source.

A.2. Extracting Likelihood Values in the Presence of the Halo

In the presence of a halo, it is nontrivial to extract the
likelihood values, $\mathcal{L}(B, \theta_i | \mathcal{D}_{\text{LAT},i})$, using the built-in function-
ality of fermipy. The reason is that the cascade flux
depends on the intrinsic spectral parameters of the blazar.
Therefore, we loop over the spectral parameters of the central
blazar and fix them at each step. Then, we add the
corresponding reweighted halo template and extract
the likelihood values. More specifically, we start from the
optimized ROI model of the LAT data for each source $i$ as
described in Section 2.2. We loop through values of the
source-intrinsic spectral index, $\Gamma_{\text{int}} = \Gamma_{\text{int},0} + \Delta\Gamma_{\text{int}}$, where
$\Gamma_{\text{int},0}$ denotes the best-fit value from LAT data (see Table 1)
and $\Delta\Gamma_{\text{int}}$ ranges from −0.75 to 0.75 with a step size of
0.075. We further sample the exponential cutoff energy of the
spectrum, $E_{\text{cut}} = 1, 3, \dots, 13$ TeV (which itself is not
constrained from LAT data) and the spectral normalization
$N_0$. The choice is motivated by the energy range of the
H.E.S.S. spectra. For each combination of the spectral
parameters, we reweight the cascade template, add it to the
ROI, and extract the multidimensional profile likelihood cube
for the model with the halo emission to describe the data. The
normalizations of the other (diffuse) sources in the ROI are
optimized in each step of the loop.

<footnote n="46">⁴⁶ https://fermipy.readthedocs.io/en/latest/</footnote>

<footnote n="47">⁴⁷ https://fermi.gsfc.nasa.gov/ssc/data/analysis/</footnote>

<footnote n="48">⁴⁸ See https://fermi.gsfc.nasa.gov/ssc/data/access/lat/BackgroundModels.html. For the isotropic diffuse emission, we calculate an average of the templates for the PSF classes 0–2.</footnote>
<figure><img src="image_7.png" /><figcaption>**Figure 5.** TS maps as in Figure 1 for all other considered sources. Again, no significant excess is observed close to any considered source. The additional source in the ROI of 1ES 1101-232 is marked with a lowercase j and its coordinates.</figcaption></figure>

<table><caption>**Table 2**

Best-fit Parameters for the LAT Spectra of the Considered EHBL Sources without a Halo Component</caption><thead><tr><th>Source Name</th><th>N<sub>0</sub> (10<sup>-14</sup> MeV<sup>-1</sup> cm<sup>-2</sup> s<sup>-1</sup>)</th><th>E<sub>0</sub> (GeV)</th><th>Γ</th><th>Γ<sub>int</sub></th></tr></thead><tbody><tr><td>1ES 0229+200</td><td>1.02 ± 0.14</td><td>5.80</td><td>1.77 ± 0.09</td><td>1.73 ± 0.12</td></tr><tr><td>1ES 0347-121</td><td>1.90 ± 0.21</td><td>4.82</td><td>1.81 ± 0.08</td><td>1.81 ± 0.11</td></tr><tr><td>PKS 0548-322</td><td>2.17 ± 0.32</td><td>3.29</td><td>1.99 ± 0.12</td><td>2.06 ± 0.17</td></tr><tr><td>1ES 1101-232</td><td>2.09 ± 0.23</td><td>4.84</td><td>1.68 ± 0.07</td><td>1.65 ± 0.09</td></tr><tr><td>H2356-309</td><td>3.49 ± 0.35</td><td>3.70</td><td>1.78 ± 0.07</td><td>1.81 ± 0.09</td></tr></tbody></table>

## Appendix B
### Details on the H.E.S.S. Data Analysis

In this appendix, details on the analysis choices for the H.E. S.S. observations are given in addition to the main selection criteria described in Section 2.3. For each source, the runs are stacked onto one data set per source for runs that pass our selection given in the main text in combination with an offset between telescope pointing and source direction of ≲2° and a zenith angle of Z ≲ 60°. The background in the On region is estimated from multiple reflected “Off” regions of the same size that have a minimum separation to the source location of

0°2 (Aharonian et al. 2006a). The energy threshold $E_{thr}$ of this stacked data set is computed from the condition that the energy bias of the energy dispersion matrix is below 10%. The median zenith angles $\langle Z \rangle$ and values of $E_{thr}$ are provided in Table 3. We choose an energy binning of the recorded events of 24 bins per decade. Using the MINUIT back end in gammapy, we fit a simple power law without (with) EBL absorption to the observations. The best-fit parameters, $N_0$ ($N_{0,\text{int}}$) and $\Gamma$ ($\Gamma_{\text{int}}$), are also provided in Table 3. The EBL absorption effect is clearly visible as $\Gamma_{\text{int}} < \Gamma$ for all sources. The energy normalization $E_0$ is determined from the fit without EBL absorption such that the off-diagonal terms of the covariance matrix are
<table><caption>**Table 3**
Mean Zenith Angle (Z), Energy Threshold $E_{thr}$, and Best-fit Parameters for the H.E.S.S. Spectra of the Considered EHBL Sources without a Halo Component</caption><thead><tr><th>Source Name</th><th>⟨Z⟩ (deg)</th><th>E<sub>thr</sub> (TeV)</th><th>N<sub>0</sub> (10<sup>-12</sup> TeV<sup>-1</sup> cm<sup>-2</sup> s<sup>-1</sup>)</th><th>N<sub>0</sub><sub>int</sub> (10<sup>-12</sup> TeV<sup>-1</sup> cm<sup>-2</sup> s<sup>-1</sup>)</th><th>E<sub>0</sub> (TeV)</th><th>Γ</th><th>Γ<sub>int</sub></th></tr></thead><tbody><tr><td>1ES 0229+200</td><td>44.4</td><td>0.46</td><td>0.32 ± 0.02</td><td>1.64 ± 0.11</td><td>1.14</td><td>2.81 ± 0.11</td><td>1.76 ± 0.12</td></tr><tr><td>1ES 0347-121</td><td>16.1</td><td>0.22</td><td>4.70 ± 0.35</td><td>11.94 ± 0.09</td><td>0.43</td><td>3.40 ± 0.14</td><td>2.12 ± 0.15</td></tr><tr><td>PKS 0548-322</td><td>11.0</td><td>0.22</td><td>0.44 ± 0.05</td><td>0.85 ± 0.09</td><td>0.89</td><td>2.44 ± 0.11</td><td>1.92 ± 0.12</td></tr><tr><td>1ES 1101-232</td><td>21.3</td><td>0.16</td><td>2.75 ± 0.18</td><td>10.13 ± 0.64</td><td>0.54</td><td>2.93 ± 0.09</td><td>1.66 ± 0.09</td></tr><tr><td>H2356-309</td><td>13.5</td><td>0.22</td><td>3.60 ± 0.18</td><td>9.14 ± 0.46</td><td>0.46</td><td>3.21 ± 0.09</td><td>2.10 ± 0.09</td></tr></tbody></table>

minimal. The resulting point-source spectra are shown in Figure 6 in Appendix D.1 as a red bow tie and red points. For the spectral points, we use 9 bins between $E_{thr}$ and $10^{1.25}$ TeV. A cross-check analysis with an independent analysis chain (de Naurois & Rolland 2009) yields consistent results.

## Appendix C
### Details on the Cascade Simulations with CRPropa

The CRPropa simulations are carried out in the large sphere observer framework, with the source placed at the origin of a sphere with radius equal to the comoving distance D to the source. We use a standard $\Lambda$CDM cosmology with $H_0 = 67.3$ km s$^{-1}$ Mpc$^{-1}$ and a matter density of $\Omega_m = 0.315$ (these are default CRPropa values). From the source, gamma rays are injected into a cone with opening angle $\theta_{jet} = 6^\circ$ (which is in the range of typically observed angles; e.g., Jorstad et al. 2017), with its central axis pointing along the negative x-axis in the simulation coordinate system. We include energy losses through pair production, IC scattering, synchrotron emission of electrons in the IGMF, and propagation of particles in an expanding universe. Within the propagation module of CRPropa, we set the tolerance parameter, which determines the numerical precision, to $10^{-9}$ and set the minimum step length of the simulation such that time delays $\tau$ down to 1 yr can be resolved. In order to reduce the computational time, we stop tracing particles if their energy drops below 0.1 GeV or if their total propagating distance becomes larger than 4 Gpc. We inject monochromatic gamma rays with energies $\epsilon$ equal to the central energies of bins $\Delta\epsilon$ between 1 GeV and $\approx 31.6$ TeV using 4 bins per decade. The value of the maximum energy is limited by the available computational resources for the simulations. As shown in Section 3.2, the intrinsic blazar spectra will tend to have lower cutoff energies so that this choice does not affect our results. The time delay $\tau$ of the cascade gamma rays is calculated through $c\tau \approx (d-D)$, where c is the speed of light. To mitigate simulation noise introduced through a finite number of injected gamma rays, the sky maps are adaptively smoothed using the ASMOOTH algorithm (Ebeling et al. 2006). From the CRPropa simulation, we obtain a multidimensional histogram of photon counts $N$. We can convert this histogram into the number of counts arriving per injected energy $\epsilon$, observed energy $E$, solid angle $\Omega$, time delay interval $d\tau$, and injected particle in the interval $\Delta\epsilon$ as

$$ \frac{dN}{d\epsilon dEd\tau d\Omega} = \frac{1}{N_{\text{inj}}(\Delta\epsilon) \Delta\epsilon \Delta E \Delta\tau \Delta\Omega}, \quad (C1) $$

where $\Delta X$ is the chosen bin width of the CRPropa histogram for quantity X. To speed up computations, the number of injected particles $N_{inj}$ decreases for larger injected energies and thus depends on the injected energy. For a single simulation, $N_{inj} = 10$ for $\epsilon < 1$ TeV, $N_{inj} = 5$ for $1 \le \epsilon < 10$ TeV, and

$N_{inj} = 3$ for $\epsilon \ge 10$ TeV. For magnetic fields $\log_{10}(B/G) \le -15$, 1000 independent simulations are performed that are combined into a single spectrum. For higher B fields we run more simulations as more particles are deflected at large angles and discarded (see Section C.1), which increases the noise of the simulations. We increase the number of simulations by 200 for each tested IGMF, resulting in 1800 simulations for $\log_{10}(B/G) = -13$.

To obtain the cascade flux that arrives within some maximum time delay $\tau_{max}$ for an arbitrary injected (time-averaged) photon spectrum $dN/d\epsilon$ (e.g., in units $s^{-1}$ cm$^{-2}$ eV$^{-1}$), we need to convolve the initial cascade spectrum in Equation (C1) with the injected spectrum and integrate over the delay times,

$$ \begin{aligned} \frac{dN}{dEd\Omega} &= \int_{0}^{\infty} d\epsilon \int_{0}^{\infty} d\epsilon' \frac{dN}{d\epsilon'} \int_{0}^{\tau_{\max}} d\tau \frac{dN}{d\epsilon' dEd\tau d\Omega} \\ &\approx \sum_{i} \sum_{j} \Delta\epsilon_i \Delta\tau_j w_i \left( \frac{dN}{d\epsilon dEd\tau d\Omega} \right)_{ij} \end{aligned} \quad (C2) $$

In the last step, we have discretized the two integrals and introduced the spectral weights $w_i$ for the ith energy bin of the injected spectrum $\Delta\epsilon_i$,

$$ w_i = \int_{\Delta\epsilon_i} \frac{dN}{d\epsilon} d\epsilon. \quad (C3) $$

### C.1. Applied Rotations in the Large Sphere Observer Framework

Each event is recorded when it hits the sphere at some position $x$ in the global simulation coordinate system $S_{sim}$ with Cartesian unit vectors $\hat{x}_{sim}$, $\hat{y}_{sim}$, and $\hat{z}_{sim}$. For each event the final (initial) unit momentum vector $p$ ($p_0$), final (initial) energy $E(\epsilon)$, and propagation length $d$ are recorded. The sphere has a radius equal to the comoving distance $D$ to the source, and the source is located at the center of the sphere with Cartesian coordinates $x_0 = (0, 0, 0)$. For the simulation, we assume a source that emits particles uniformly into a cone with aperture $\theta_{jet}$ with the jet axis $j$ pointing in the $(-1, 0, 0)$ direction.

We adopt the large sphere observer framework and follow the same method as in Alves Batista et al. (2016b), placing an observer at each hit position $x$ on the surface of the sphere and defining a local coordinate system $S_{obs}$. In $S_{obs}$, the unit vector of the $z_{obs}$-axis points along the connection between the hit location and the source, i.e., the LOS, defined by the unit vector $\xi = \xi_x \hat{x}_{sim} + \xi_y \hat{y}_{sim} + \xi_z \hat{z}_{sim} = (x - x_0)/|x - x_0|$. With this definition, we see that $\hat{z}_{obs} = -\xi$. The $\hat{x}_{obs}$ and $\hat{y}_{obs}$ unit vectors are then defined in such a way that $S_{obs}$ is a right-handed
coordinate system. Introducing the definitions,

$$ \rho = \sqrt{\xi_x^2 + \xi_y^2} \qquad (C4) $$

$$ r = \sqrt{\xi_x^2 + \xi_y^2 + \xi_z^2} \qquad (C5) $$

$$ \cos \phi = \frac{\xi_x}{\rho} \qquad (C6) $$

$$ \sin \phi = \frac{\xi_y}{\rho} \qquad (C7) $$

$$ \cos \vartheta = \frac{\xi_z}{r} \qquad (C8) $$

$$ \sin \vartheta = \frac{\rho}{r}, \qquad (C9) $$

the unit vectors $\hat{x}_{\text{obs}}$ and $\hat{y}_{\text{obs}}$ in $S_{\text{sim}}$ are found to be

$$ \hat{x}_{\text{obs}} = (-\sin \phi, \cos \phi, 0)^T, \qquad (C10) $$

$$ \hat{y}_{\text{obs}} = (\cos \phi \cos \vartheta, \sin \phi \cos \vartheta, -\sin \vartheta)^T. \qquad (C11) $$

We proceed by projecting the momentum vectors $\mathbf{p}$ and $\mathbf{p}_0$ into the new basis,

$$ \mathbf{p}_{\text{obs}} = \langle \mathbf{p}, \hat{x}_{\text{obs}} \rangle \hat{x}_{\text{obs}} + \langle \mathbf{p}, \hat{y}_{\text{obs}} \rangle \hat{y}_{\text{obs}} + \langle \mathbf{p}, \hat{z}_{\text{obs}} \rangle \hat{z}_{\text{obs}}, \qquad (C12) $$

$$ \mathbf{p}_{0,\text{obs}} = \langle \mathbf{p}_0, \hat{x}_{\text{obs}} \rangle \hat{x}_{\text{obs}} + \langle \mathbf{p}_0, \hat{y}_{\text{obs}} \rangle \hat{y}_{\text{obs}} + \langle \mathbf{p}_0, \hat{z}_{\text{obs}} \rangle \hat{z}_{\text{obs}}. \qquad (C13) $$

We can now consider arbitrary orientations of the jet axis with respect to the LOS by defining the jet vector in $S_{\text{obs}}$ to be

$$ \mathbf{j}_{\text{obs}} = \cos \phi_{\text{jet}} \sin \vartheta_{\text{jet}} \hat{x}_{\text{obs}} + \sin \phi_{\text{jet}} \sin \vartheta_{\text{jet}} \hat{y}_{\text{obs}} + \cos \vartheta_{\text{jet}} \hat{z}_{\text{obs}}, \qquad (C14) $$

with some values for $\phi_{\text{jet}}$ and $\vartheta_{\text{jet}}$. The $\vartheta_{\text{jet}}$ angle is related to the angle $\theta_{\text{obs}}$ between the LOS and the jet axis through $\theta_{\text{obs}} = \pi/2 - \vartheta_{\text{jet}}$ (here we choose $\phi_{\text{jet}} = \pi$ and $\vartheta_{\text{jet}} = \pi/2$).

Only photons will be included in the cascade flux for which the angle $\alpha_{\text{obs}}$ between the initial photon momentum $\mathbf{p}_{0,\text{obs}}$ and the jet axis $\mathbf{j}_{\text{obs}}$ is smaller than the assumed half-opening angle of the jet, $\theta_{\text{jet}}/2$,

$$ \cos \alpha_{\text{obs}} = \langle -\mathbf{j}_{\text{obs}}, \mathbf{p}_{0,\text{obs}} \rangle \geq \cos(\theta_{\text{jet}}/2). \qquad (C15) $$

In this way, we can use the same simulation for arbitrary angles between the jet axis and the LOS.

## Appendix D
### Details on Fits with Halo Components

In this appendix, we provide additional details on the combined fits of LAT and H.E.S.S. data in the presence of a halo.

#### D.1. Best-fit Parameters

The best-fit parameters for the spectrum in Equation (1) for all sources for all tested values of the IGMF and $t_{\text{max}}$ are reported in Table 4. For simplicity, $E_0$ is fixed to 1 TeV. Example spectra for all sources for three different values of the IGMF and $t_{\text{max}} = 10^7$ yr are shown in Figure 6 (same as

Figure 3 in the main text). From Figure 6 and the best-fit values of $\kappa$ in the last column of Table 4, it becomes clear that exponential cutoffs at lower energies are preferred for low magnetic fields. This again shows the absence of a halo, as a lower-energy cutoff leads to fewer cascade photons. Figure 6 also illustrates that an increasing magnetic field leads to a suppression of the cascade at lower energies. Under the assumption that the coherence length of the magnetic field is much larger than the IC cooling length, the deflection of the pairs roughly scales as $E_e^{-2}B$, where $E_e$ is the electron energy (Neronov & Semikoz 2009). As $B$ increases, the low-energy pairs in particular will exhibit a large deflection and the upscattered cascade emission will not reach the observer. As a consequence, the peak of the cascade emission shifts toward higher energies as can be seen in each row of Figure 6.

#### D.2. Likelihood Profiles

In Figure 7, we show the likelihood ratios $\lambda(B)$ combined for all sources as defined in Equation (2) for our choice of different blazar activity times (black solid lines). In addition, we show the likelihood ratios $-2\Delta \ln L_i$ for the individual sources as colored lines. It is evident that shorter activity times lead to weaker constraints on the IGMF, as the likelihood curves are consistently shifted toward lower values of $B$.

#### D.3. Additional Systematic Uncertainties

In the following, the effects of systematic uncertainties in addition to the unknown blazar activity time are discussed.

Increasing the jet opening angle $\theta_{\text{jet}}$ has a similar effect to increasing $t_{\text{max}}$. Larger jet opening angles will lead to larger halos and larger cascade fluxes as fewer gamma rays are rejected in our parallel transport prescription (see Section C.1). This makes sense, as photons emitted along the edge of the jet could be deflected toward the observer. However, given the observed spread of $\theta_{\text{jet}}$ (e.g., Jorstad et al. 2017), the uncertainty is subdominant compared to the uncertain activity time. Similarly, the choice of the EBL model introduces a negligible uncertainty, as it typically affects the number of generated pairs of the order of ~5% in the energy and redshift range under consideration (Meyer et al. 2016).

Another potentially important systematic uncertainty is the difference in the absolute energy scales of H.E.S.S. and the LAT. For the LAT, the error on the energy scale is of the order of a few percent (Ackermann et al. 2012), whereas the reconstructed energies of gamma rays observed by IACTs have an uncertainty of typically 10% (Hofmann et al. 2000). We can allow for a relative shift in energy between the two instruments by adding a scaling parameter $s$ in our fitting procedure, such that the model fluxes fitted to H.E.S.S. data are changed to $\phi(E) \rightarrow \phi(E(1+s))$. We constrain $s$ to lie within $\pm 0.2$. Repeating the fits for all sources, we find that for low magnetic fields $s$ tends to its lowest allowed value. This minimizes the cascade emission, as the H.E.S.S. data are effectively shifted to lower fluxes. With this systematic term included, the limits on the magnetic field are weakened and $B < 5.3 \times 10^{-16}$ G is ruled out for $t_{\text{max}} = 10$ yr.
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Best-fit spectra assuming a point-source model in addition to the halo flux templates for a combined LAT and H.E.S.S. analysis. An activity time of $t_{max} = 10^7$ yr is assumed for all sources. The columns show different field strengths $B \in [10^{-16}, 10^{-15}, 10^{-14}]$ G, and each source is shown in one row. A coherence length of 1 Mpc is assumed throughout, together with the EBL model of Domínguez et al. (2011). Blue and red markers and lines show the point-source spectra measured with the LAT and H.E.S.S., respectively.</figcaption></figure>
<table><caption>**Table 4**

Best-fit Parameters for the Combined LAT and H.E.S.S. Data in the Presence of a Halo</caption><thead><tr><th>Source<br>Name</th><th>IGMF Strength<br>(log<sub>10</sub>(B/G))</th><th>N<sub>0,int</sub><br>(10<sup>-12</sup> TeV<sup>-1</sup> cm<sup>-2</sup> s<sup>-1</sup>)</th><th>Γ<sub>int</sub></th><th>E<sub>cut</sub><br>(TeV)</th></tr></thead><tbody><tr><td colspan="5"><strong>t<sub>max</sub> = 10 yr</strong></td></tr><tr><td>1ES 0229+200</td><td>-16.0</td><td>2.114 ± 0.223</td><td>1.632 ± 0.004</td><td>8.103 ± 3.022</td></tr><tr><td>1ES 0229+200</td><td>-15.5</td><td>2.160 ± 0.066</td><td>1.633 ± 0.003</td><td>8.865 ± 2.414</td></tr><tr><td>1ES 0229+200</td><td>-15.0</td><td>2.236 ± 6.016</td><td>1.639 ± 0.297</td><td>8.970 ± 10.878</td></tr><tr><td>1ES 0229+200</td><td>-14.5</td><td>2.260 ± 0.073</td><td>1.635 ± 0.007</td><td>10.031 ± 3.545</td></tr><tr><td>1ES 0229+200</td><td>-14.0</td><td>2.261 ± 0.085</td><td>1.639 ± 0.009</td><td>10.996 ± 0.500</td></tr><tr><td>1ES 0229+200</td><td>-13.5</td><td>2.300 ± 0.105</td><td>1.638 ± 0.011</td><td>11.000 ± 0.080</td></tr><tr><td>1ES 0229+200</td><td>-13.0</td><td>2.308 ± 0.092</td><td>1.638 ± 0.010</td><td>10.973 ± 1.487</td></tr><tr><td>1ES 0347-121</td><td>-16.0</td><td>2.945 ± 0.063</td><td>1.633 ± 0.005</td><td>3.970 ± 1.244</td></tr><tr><td>1ES 0347-121</td><td>-15.5</td><td>2.902 ± 0.119</td><td>1.633 ± 0.007</td><td>5.346 ± 2.173</td></tr><tr><td>1ES 0347-121</td><td>-15.0</td><td>2.907 ± 0.104</td><td>1.655 ± 0.007</td><td>5.901 ± 2.799</td></tr><tr><td>1ES 0347-121</td><td>-14.5</td><td>3.020 ± 0.093</td><td>1.633 ± 0.006</td><td>5.134 ± 2.421</td></tr><tr><td>1ES 0347-121</td><td>-14.0</td><td>3.137 ± 0.087</td><td>1.631 ± 0.005</td><td>4.019 ± 1.721</td></tr><tr><td>1ES 0347-121</td><td>-13.5</td><td>3.244 ± 0.066</td><td>1.633 ± 0.004</td><td>3.610 ± 1.410</td></tr><tr><td>1ES 0347-121</td><td>-13.0</td><td>2.797 ± 0.148</td><td>1.659 ± 0.012</td><td>11.000 ± 0.014</td></tr><tr><td>1ES 1101-232</td><td>-16.0</td><td>3.360 ± 0.111</td><td>1.632 ± 0.003</td><td>5.948 ± 1.609</td></tr><tr><td>1ES 1101-232</td><td>-15.5</td><td>3.109 ± 0.025</td><td>1.627 ± 0.002</td><td>10.906 ± 1.529</td></tr><tr><td>1ES 1101-232</td><td>-15.0</td><td>3.438 ± 0.056</td><td>1.633 ± 0.002</td><td>9.652 ± 3.680</td></tr><tr><td>1ES 1101-232</td><td>-14.5</td><td>3.625 ± 0.174</td><td>1.633 ± 0.008</td><td>9.517 ± 3.018</td></tr><tr><td>1ES 1101-232</td><td>-14.0</td><td>3.607 ± 0.184</td><td>1.627 ± 0.012</td><td>10.998 ± 0.315</td></tr><tr><td>1ES 1101-232</td><td>-13.5</td><td>3.706 ± 3.540</td><td>1.633 ± 0.010</td><td>10.777 ± 88.543</td></tr><tr><td>1ES 1101-232</td><td>-13.0</td><td>3.738 ± 0.540</td><td>1.632 ± 0.012</td><td>10.198 ± 14.941</td></tr><tr><td>H2356-309</td><td>-16.0</td><td>2.754 ± 0.075</td><td>1.717 ± 0.005</td><td>2.832 ± 0.573</td></tr><tr><td>H2356-309</td><td>-15.5</td><td>2.754 ± 0.690</td><td>1.726 ± 0.046</td><td>2.924 ± 0.895</td></tr><tr><td>H2356-309</td><td>-15.0</td><td>2.780 ± 0.591</td><td>1.727 ± 0.038</td><td>2.907 ± 1.988</td></tr><tr><td>H2356-309</td><td>-14.5</td><td>2.890 ± 0.024</td><td>1.710 ± 0.002</td><td>2.637 ± 0.503</td></tr><tr><td>H2356-309</td><td>-14.0</td><td>2.997 ± 0.049</td><td>1.707 ± 0.005</td><td>2.383 ± 0.437</td></tr><tr><td>H2356-309</td><td>-13.5</td><td>2.997 ± 0.046</td><td>1.697 ± 0.003</td><td>2.383 ± 0.436</td></tr><tr><td>H2356-309</td><td>-13.0</td><td>3.001 ± 0.058</td><td>1.696 ± 0.004</td><td>2.401 ± 0.450</td></tr><tr><td>PKS 0548-322</td><td>-16.0</td><td>7.520 ± 0.071</td><td>1.774 ± 0.008</td><td>7.621 ± 0.077</td></tr><tr><td>PKS 0548-322</td><td>-15.5</td><td>7.158 ± 0.161</td><td>1.803 ± 0.004</td><td>10.203 ± 4.209</td></tr><tr><td>PKS 0548-322</td><td>-15.0</td><td>7.449 ± 0.310</td><td>1.794 ± 0.008</td><td>9.973 ± 5.004</td></tr><tr><td>PKS 0548-322</td><td>-14.5</td><td>7.449 ± 0.294</td><td>1.789 ± 0.007</td><td>10.906 ± 3.623</td></tr><tr><td>PKS 0548-322</td><td>-14.0</td><td>7.853 ± 0.313</td><td>1.780 ± 0.008</td><td>9.150 ± 5.184</td></tr><tr><td>PKS 0548-322</td><td>-13.5</td><td>7.853 ± 0.306</td><td>1.780 ± 0.008</td><td>9.150 ± 5.177</td></tr><tr><td>PKS 0548-322</td><td>-13.0</td><td>7.853 ± 0.307</td><td>1.780 ± 0.009</td><td>9.150 ± 5.211</td></tr><tr><td colspan="5"><strong>t<sub>max</sub> = 10<sup>4</sup> yr</strong></td></tr><tr><td>1ES 0229+200</td><td>-16.0</td><td>2.257 ± 0.061</td><td>1.633 ± 0.001</td><td>3.171 ± 0.456</td></tr><tr><td>1ES 0229+200</td><td>-15.5</td><td>2.274 ± 0.038</td><td>1.632 ± 0.002</td><td>4.106 ± 0.690</td></tr><tr><td>1ES 0229+200</td><td>-15.0</td><td>2.142 ± 0.046</td><td>1.633 ± 0.001</td><td>6.321 ± 1.301</td></tr><tr><td>1ES 0229+200</td><td>-14.5</td><td>2.146 ± 0.067</td><td>1.632 ± 0.004</td><td>8.070 ± 2.014</td></tr><tr><td>1ES 0229+200</td><td>-14.0</td><td>2.206 ± 0.083</td><td>1.633 ± 0.005</td><td>8.673 ± 2.458</td></tr><tr><td>1ES 0229+200</td><td>-13.5</td><td>2.217 ± 0.344</td><td>1.646 ± 0.013</td><td>9.921 ± 9.287</td></tr><tr><td>1ES 0229+200</td><td>-13.0</td><td>2.227 ± 0.106</td><td>1.644 ± 0.011</td><td>10.979 ± 1.248</td></tr><tr><td>1ES 0347-121</td><td>-16.0</td><td>4.398 ± 0.052</td><td>1.520 ± 0.003</td><td>1.000 ± 0.000</td></tr><tr><td>1ES 0347-121</td><td>-15.5</td><td>3.618 ± 0.041</td><td>1.626 ± 0.002</td><td>1.269 ± 0.210</td></tr><tr><td>1ES 0347-121</td><td>-15.0</td><td>2.997 ± 0.065</td><td>1.633 ± 0.002</td><td>2.976 ± 0.319</td></tr><tr><td>1ES 0347-121</td><td>-14.5</td><td>2.948 ± 0.009</td><td>1.637 ± 0.005</td><td>4.149 ± 0.013</td></tr><tr><td>1ES 0347-121</td><td>-14.0</td><td>2.965 ± 0.056</td><td>1.636 ± 0.004</td><td>5.046 ± 1.309</td></tr><tr><td>1ES 0347-121</td><td>-13.5</td><td>2.964 ± 0.042</td><td>1.644 ± 0.003</td><td>5.222 ± 1.402</td></tr><tr><td>1ES 0347-121</td><td>-13.0</td><td>2.960 ± 0.085</td><td>1.656 ± 0.006</td><td>5.022 ± 2.266</td></tr><tr><td>1ES 1101-232</td><td>-16.0</td><td>3.534 ± 0.087</td><td>1.633 ± 0.002</td><td>1.747 ± 0.258</td></tr><tr><td>1ES 1101-232</td><td>-15.5</td><td>3.366 ± 0.022</td><td>1.633 ± 0.001</td><td>2.897 ± 0.519</td></tr><tr><td>1ES 1101-232</td><td>-15.0</td><td>3.364 ± 0.059</td><td>1.632 ± 0.003</td><td>4.254 ± 0.951</td></tr><tr><td>1ES 1101-232</td><td>-14.5</td><td>3.378 ± 0.069</td><td>1.631 ± 0.004</td><td>6.024 ± 1.563</td></tr><tr><td>1ES 1101-232</td><td>-14.0</td><td>3.197 ± 0.024</td><td>1.625 ± 0.001</td><td>10.368 ± 2.810</td></tr><tr><td>1ES 1101-232</td><td>-13.5</td><td>3.365 ± 0.093</td><td>1.632 ± 0.005</td><td>10.964 ± 1.380</td></tr><tr><td>1ES 1101-232</td><td>-13.0</td><td>3.538 ± 0.084</td><td>1.631 ± 0.005</td><td>10.991 ± 0.715</td></tr></tbody></table>
<table><caption>**Table 4**
(Continued)</caption><thead><tr><th>Source<br>Name</th><th>IGMF Strength<br>(log<sub>10</sub>(B/G))</th><th>N<sub>0,int</sub><br>(10<sup>-12</sup> TeV<sup>-1</sup> cm<sup>-2</sup> s<sup>-1</sup>)</th><th>Γ<sub>int</sub></th><th>E<sub>cut</sub><br>(TeV)</th></tr></thead><tbody><tr><td>H2356-309</td><td>-16.0</td><td>3.555 ± 3.139</td><td>1.633 ± 0.002</td><td>1.255 ± 2.016</td></tr><tr><td>H2356-309</td><td>-15.5</td><td>3.696 ± 0.036</td><td>1.633 ± 0.002</td><td>1.328 ± 0.148</td></tr><tr><td>H2356-309</td><td>-15.0</td><td>3.494 ± 0.043</td><td>1.634 ± 0.002</td><td>1.640 ± 0.214</td></tr><tr><td>H2356-309</td><td>-14.5</td><td>2.789 ± 0.004</td><td>1.725 ± 0.006</td><td>2.704 ± 0.004</td></tr><tr><td>H2356-309</td><td>-14.0</td><td>2.795 ± 7.856</td><td>1.726 ± 0.159</td><td>2.770 ± 3.188</td></tr><tr><td>H2356-309</td><td>-13.5</td><td>2.841 ± 0.034</td><td>1.723 ± 0.002</td><td>2.703 ± 0.524</td></tr><tr><td>H2356-309</td><td>-13.0</td><td>2.900 ± 0.037</td><td>1.703 ± 0.003</td><td>2.638 ± 0.510</td></tr><tr><td>PKS 0548-322</td><td>-16.0</td><td>7.793 ± 0.080</td><td>1.745 ± 0.001</td><td>5.245 ± 1.428</td></tr><tr><td>PKS 0548-322</td><td>-15.5</td><td>7.877 ± 0.129</td><td>1.745 ± 0.002</td><td>5.870 ± 1.744</td></tr><tr><td>PKS 0548-322</td><td>-15.0</td><td>7.680 ± 0.250</td><td>1.768 ± 0.008</td><td>6.844 ± 2.431</td></tr><tr><td>PKS 0548-322</td><td>-14.5</td><td>7.610 ± 0.237</td><td>1.768 ± 0.006</td><td>8.011 ± 3.147</td></tr><tr><td>PKS 0548-322</td><td>-14.0</td><td>7.539 ± 0.284</td><td>1.783 ± 0.007</td><td>9.001 ± 4.146</td></tr><tr><td>PKS 0548-322</td><td>-13.5</td><td>7.543 ± 0.301</td><td>1.782 ± 0.009</td><td>10.057 ± 5.377</td></tr><tr><td>PKS 0548-322</td><td>-13.0</td><td>7.554 ± 0.407</td><td>1.786 ± 0.011</td><td>10.222 ± 5.494</td></tr><tr><td>t<sub>max</sub> = 10<sup>7</sup> yr</td><td></td><td></td><td></td><td></td></tr><tr><td>1ES 0229+200</td><td>-16.0</td><td>2.217 ± 0.042</td><td>1.633 ± 0.001</td><td>2.880 ± 0.408</td></tr><tr><td>1ES 0229+200</td><td>-15.5</td><td>2.166 ± 0.030</td><td>1.633 ± 0.001</td><td>3.667 ± 0.552</td></tr><tr><td>1ES 0229+200</td><td>-15.0</td><td>2.144 ± 0.044</td><td>1.633 ± 0.001</td><td>5.209 ± 0.959</td></tr><tr><td>1ES 0229+200</td><td>-14.5</td><td>2.129 ± 0.165</td><td>1.633 ± 0.003</td><td>7.099 ± 2.412</td></tr><tr><td>1ES 0229+200</td><td>-14.0</td><td>2.201 ± 0.104</td><td>1.632 ± 0.005</td><td>7.814 ± 2.125</td></tr><tr><td>1ES 0229+200</td><td>-13.5</td><td>2.225 ± 0.124</td><td>1.635 ± 0.011</td><td>9.109 ± 3.505</td></tr><tr><td>1ES 0229+200</td><td>-13.0</td><td>2.225 ± 0.080</td><td>1.635 ± 0.008</td><td>10.584 ± 3.262</td></tr><tr><td>1ES 0347-121</td><td>-16.0</td><td>4.133 ± 0.418</td><td>1.520 ± 0.001</td><td>1.000 ± 0.000</td></tr><tr><td>1ES 0347-121</td><td>-15.5</td><td>4.462 ± 0.060</td><td>1.520 ± 0.002</td><td>1.000 ± 0.005</td></tr><tr><td>1ES 0347-121</td><td>-15.0</td><td>3.257 ± 0.055</td><td>1.632 ± 0.003</td><td>1.638 ± 0.324</td></tr><tr><td>1ES 0347-121</td><td>-14.5</td><td>2.960 ± 0.317</td><td>1.632 ± 0.004</td><td>3.020 ± 0.482</td></tr><tr><td>1ES 0347-121</td><td>-14.0</td><td>2.824 ± 0.117</td><td>1.632 ± 0.002</td><td>5.152 ± 2.231</td></tr><tr><td>1ES 0347-121</td><td>-13.5</td><td>2.655 ± 1.179</td><td>1.672 ± 0.088</td><td>8.368 ± 6.833</td></tr><tr><td>1ES 0347-121</td><td>-13.0</td><td>2.656 ± 0.241</td><td>1.672 ± 0.023</td><td>10.998 ± 0.100</td></tr><tr><td>1ES 1101-232</td><td>-16.0</td><td>5.335 ± 0.435</td><td>1.520 ± 0.002</td><td>1.001 ± 0.012</td></tr><tr><td>1ES 1101-232</td><td>-15.5</td><td>5.183 ± 0.081</td><td>1.520 ± 0.002</td><td>1.167 ± 0.139</td></tr><tr><td>1ES 1101-232</td><td>-15.0</td><td>4.893 ± 0.042</td><td>1.520 ± 0.001</td><td>1.766 ± 0.262</td></tr><tr><td>1ES 1101-232</td><td>-14.5</td><td>4.651 ± 0.149</td><td>1.553 ± 0.010</td><td>2.357 ± 0.224</td></tr><tr><td>1ES 1101-232</td><td>-14.0</td><td>3.810 ± 0.044</td><td>1.599 ± 0.002</td><td>4.753 ± 1.205</td></tr><tr><td>1ES 1101-232</td><td>-13.5</td><td>3.843 ± 0.080</td><td>1.600 ± 0.004</td><td>5.675 ± 1.756</td></tr><tr><td>1ES 1101-232</td><td>-13.0</td><td>3.858 ± 0.070</td><td>1.604 ± 0.004</td><td>6.623 ± 2.599</td></tr><tr><td>H2356-309</td><td>-16.0</td><td>3.709 ± 0.025</td><td>1.632 ± 0.001</td><td>1.085 ± 0.100</td></tr><tr><td>H2356-309</td><td>-15.5</td><td>3.411 ± 0.096</td><td>1.632 ± 0.003</td><td>1.421 ± 0.176</td></tr><tr><td>H2356-309</td><td>-15.0</td><td>3.486 ± 0.042</td><td>1.632 ± 0.002</td><td>1.527 ± 0.186</td></tr><tr><td>H2356-309</td><td>-14.5</td><td>3.331 ± 0.046</td><td>1.663 ± 0.003</td><td>1.739 ± 0.239</td></tr><tr><td>H2356-309</td><td>-14.0</td><td>2.875 ± 0.107</td><td>1.706 ± 0.007</td><td>2.571 ± 0.526</td></tr><tr><td>H2356-309</td><td>-13.5</td><td>2.877 ± 0.042</td><td>1.703 ± 0.003</td><td>2.642 ± 0.501</td></tr><tr><td>H2356-309</td><td>-13.0</td><td>2.878 ± 0.040</td><td>1.709 ± 0.003</td><td>2.646 ± 0.510</td></tr><tr><td>PKS 0548-322</td><td>-16.0</td><td>7.790 ± 0.097</td><td>1.745 ± 0.001</td><td>5.039 ± 1.201</td></tr><tr><td>PKS 0548-322</td><td>-15.5</td><td>7.781 ± 0.084</td><td>1.745 ± 0.001</td><td>5.821 ± 1.684</td></tr><tr><td>PKS 0548-322</td><td>-15.0</td><td>8.004 ± 0.249</td><td>1.745 ± 0.004</td><td>5.895 ± 1.845</td></tr><tr><td>PKS 0548-322</td><td>-14.5</td><td>7.836 ± 0.969</td><td>1.765 ± 0.023</td><td>7.142 ± 3.571</td></tr><tr><td>PKS 0548-322</td><td>-14.0</td><td>7.483 ± 0.336</td><td>1.769 ± 0.009</td><td>9.153 ± 4.287</td></tr><tr><td>PKS 0548-322</td><td>-13.5</td><td>7.482 ± 0.297</td><td>1.783 ± 0.008</td><td>10.372 ± 5.814</td></tr><tr><td>PKS 0548-322</td><td>-13.0</td><td>7.489 ± 0.412</td><td>1.792 ± 0.011</td><td>10.988 ± 2.337</td></tr></tbody></table>
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Likelihood profiles as a function of the IGMF strength for (a) $t_{max} = 10$ yr, (b) $t_{max} = 10^4$ yr, and (c) $t_{max} = 10^7$ yr. Magnetic field values resulting in differences in the likelihood $-2\Delta \ln L > 2.71$ (indicated by the gray dotted line) are ruled out at the 95% confidence level, which is shown by the gray shaded region and arrow. The likelihood curves are smoothed using a quadratic spline interpolation.</figcaption></figure>
ORCID iDs

F. Aharonian https://orcid.org/0000-0003-1157-3915
M. Backes https://orcid.org/0000-0002-9326-6400
V. Barbosa Martins https://orcid.org/0000-0002-5085-8828
R. Batzofin https://orcid.org/0000-0002-5797-3386
Y. Becherini https://orcid.org/0000-0002-2115-2930
D. Berge https://orcid.org/0000-0002-2918-1824
M. Breuhaus https://orcid.org/0000-0003-0268-5122
R. Brose https://orcid.org/0000-0002-8312-6930
F. Brun https://orcid.org/0000-0003-0770-9007
T. Bulik https://orcid.org/0000-0003-2045-4803
T. Bylund https://orcid.org/0000-0003-2946-1313
S. Caroff https://orcid.org/0000-0002-1103-130X
S. Casanova https://orcid.org/0000-0002-6144-9122
M. Cerruti https://orcid.org/0000-0001-7891-699X
A. Chen https://orcid.org/0000-0001-6425-5692
J. Chibuee https://orcid.org/0000-0002-9875-7436
G. Cotter https://orcid.org/0000-0002-9975-1829
G. Fichet de Clairfontaine https://orcid.org/0000-0003-1143-3883
M. Filipovic https://orcid.org/0000-0002-4990-9288
G. Fontaine https://orcid.org/0000-0002-6443-5025
S. Funk https://orcid.org/0000-0002-2012-0080
G. Giavitto https://orcid.org/0000-0002-7629-6499
D. Glawion https://orcid.org/0000-0003-4865-7696
J. F. Glicenstein https://orcid.org/0000-0003-2581-1742
T. L. Holch https://orcid.org/0000-0001-5161-1168
M. Jamrozy https://orcid.org/0000-0002-0870-7778
V. Joshi https://orcid.org/0000-0003-4467-3621
B. Khélifi https://orcid.org/0000-0001-6876-5577
Nu. Komin https://orcid.org/0000-0003-3280-0582
K. Kosack https://orcid.org/0000-0001-8424-3621
D. Kostunin https://orcid.org/0000-0002-0487-0076
J.-P. Lenain https://orcid.org/0000-0001-7284-9220
F. Leuschner https://orcid.org/0000-0001-9037-0272
A. Luashvili https://orcid.org/0000-0003-4384-1638
J. Mackey https://orcid.org/0000-0002-5449-6131
D. Malyshev https://orcid.org/0000-0001-9689-2194
D. Malyshev https://orcid.org/0000-0002-9102-4854
V. Marandon https://orcid.org/0000-0001-9077-4058
P. Marchegiani https://orcid.org/0000-0001-7487-8287
G. Martí-Devesa https://orcid.org/0000-0003-0766-6473
R. Marx https://orcid.org/0000-0002-6557-4924
M. Meyer https://orcid.org/0000-0002-0738-7581
A. Mitchell https://orcid.org/0000-0003-3631-5648
L. Mohrmann https://orcid.org/0000-0002-9667-8654
A. Montanari https://orcid.org/0000-0002-3620-0173
E. Moulin https://orcid.org/0000-0003-4007-0145
J. Muller https://orcid.org/0000-0003-0004-4110
T. Murach https://orcid.org/0000-0003-1128-5008
J. Niemiec https://orcid.org/0000-0001-6036-8569
S. Ohm https://orcid.org/0000-0002-3474-2243
L. Olivera-Nieto https://orcid.org/0000-0002-9105-0518
S. Panny https://orcid.org/0000-0001-5770-3805
R. D. Parsons https://orcid.org/0000-0003-3457-9308
G. Pühlhofer https://orcid.org/0000-0003-4632-4644
M. Punch https://orcid.org/0000-0002-4710-2165
P. Reichherzer https://orcid.org/0000-0003-4513-8241
A. Reimer https://orcid.org/0000-0001-8604-7077
B. Reville https://orcid.org/0000-0002-3778-1432
G. Rowell https://orcid.org/0000-0002-9516-1581
B. Rudak https://orcid.org/0000-0003-0452-3805

E. Ruiz-Velasco https://orcid.org/0000-0001-6939-7825
V. Sahakian https://orcid.org/0000-0003-1198-0043
M. Sasaki https://orcid.org/0000-0001-5302-1866
F. Schüssler https://orcid.org/0000-0003-1500-6571
H. M. Schutte https://orcid.org/0000-0002-1769-5617
J. N. S. Shapopi https://orcid.org/0000-0002-7130-9270
S. Spencer https://orcid.org/0000-0001-5516-1205
S. Steinmassl https://orcid.org/0000-0002-2865-8563
T. Tanaka https://orcid.org/0000-0002-4383-0368
A. M. Taylor https://orcid.org/0000-0001-9473-4758
R. Terrier https://orcid.org/0000-0002-8219-4667
N. Tsuji https://orcid.org/0000-0001-7209-9204
C. van Eldik https://orcid.org/0000-0001-9669-645X
J. Veh https://orcid.org/0000-0003-4736-2167
S. J. Wagner https://orcid.org/0000-0002-7474-6062
A. Wierzcholska https://orcid.org/0000-0003-4472-7204
M. Zacharias https://orcid.org/0000-0001-5801-3945
D. Zargaryan https://orcid.org/0000-0002-2876-6433
A. A. Zdziarski https://orcid.org/0000-0002-0333-2452
S. Zouari https://orcid.org/0000-0002-5333-2004
N. Żywucka https://orcid.org/0000-0003-2644-6441
M. Meyer https://orcid.org/0000-0002-0738-7581

References

Abdalla, H., Abe, H., Acero, F., et al. 2021, JCAP, 2021, 048
Acciari, V. A., Agudo, I., Aniello, T., et al. 2023, A&A, 670, A145
Ackermann, M., Ajello, M., Allafort, A., et al. 2012, ApH, 35, 346
Ackermann, M., Ajello, M., Baldini, L., et al. 2018, ApJS, 237, 32
Aharonian, F., Akhperjanian, A. G., Anton, G., et al. 2010, A&A, 521, A69
Aharonian, F., Akhperjanian, A. G., Barres de Almeida, U., et al. 2007a, A&A, 473, L25
Aharonian, F., Akhperjanian, A. G., Bazer-Bachi, A. R., et al. 2006a, A&A, 457, 899
Aharonian, F., Akhperjanian, A. G., Bazer-Bachi, A. R., et al. 2006b, Natur, 440, 1018
Aharonian, F., Akhperjanian, A. G., Bazer-Bachi, A. R., et al. 2007b, A&A, 470, 475
Aharonian, F. A., Coppi, P. S., & Voelk, H. J. 1994, ApJL, 423, L5
Aleksić, J., Antonelli, L. A., Antoranz, P., et al. 2010, A&A, 524, A77
Alves Batista, R., Dundovic, A., Erdmann, M., et al. 2016a, JCAP, 2016, 038
Alves Batista, R., & Saveliev, A. 2021, Univ, 7, 223
Alves Batista, R., Saveliev, A., Sigl, G., & Vachaspati, T. 2016b, PhRvD, 94, 083005
Archambault, S., Archer, A., Benbow, W., et al. 2017, ApJ, 835, 288
Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33
Atwood, W., Albert, A., Baldini, L., et al. 2013, arXiv:1303.3514
Atwood, W. B., Abdo, A. A., Ackermann, M., et al. 2009, ApJ, 697, 1071
Ballet, J., Burnett, T. H., Digel, S. W., & Lott, B. 2020, arXiv:2005.11208
Broderick, A. E., Chang, P., & Pfrommer, C. 2012, ApJ, 752, 22
Caprini, C., & Gabici, S. 2015, PhRvD, 91, 123514
de Neurois, M., & Rolland, L. 2009, APh, 32, 231
Deil, C., Zanin, R., Lefaucheur, J., et al. 2017, in 35th ICRC, Gammapy - A prototype for the CTA science tools, 301 (Trieste: SISSA), 766
Dermer, C. D., Cavadini, M., Razzaque, S., et al. 2011, ApJL, 733, L21
Domínguez, A., Primack, J. R., Rosario, D. J., et al. 2011, MNRAS, 410, 2556
Durrer, R., & Neronov, A. 2013, A&ARv, 21, 62
Ebeling, H., White, D. A., & Rangarajan, F. V. N. 2006, MNRAS, 368, 65
Fermi Science Support Development Team 2019, Fermitools: Fermi Science Tools, Astrophysics Source Code Library, ascl:1905.011
Fosbury, R. A. E., & Disney, M. J. 1976, ApJL, 207, L75
Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Natur, 585, 357
H. E. S. S. Collaboration, Abramowski, A., Acero, F., et al. 2010, A&A, 516, A56
H.E.S.S. Collaboration, Abramowski, A., Aharonian, F., et al. 2014, A&A, 562, A145
Hofmann, W., Lampeitl, H., Konopelko, A., & Krawczynski, H. 2000, APh, 12, 207
Jones, D. H., Read, M. A., Saunders, W., et al. 2009, MNRAS, 399, 683
Jorstad, S. G., Marscher, A. P., Morozova, D. A., et al. 2017, ApJ, 846, 98
Kobayashi, T. 2014, JCAP, 2014, 040

Lott, B., Gasparrini, D., & Ciprini, S. 2020, arXiv:2010.08406

Meyer, M., Conrad, J., & Dickinson, H. 2016, ApJ, 827, 147

Miniati, F., Gregori, G., Reville, B., & Sarkar, S. 2018, PhRvL, 121, 021301

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

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

Parma, P., Murgia, M., de Ruiter, H. R., & Fanti, R. 2002, NewAR, 46, 313

Parsons, R. D., & Hinton, J. A. 2014, APh, 56, 26

Peña-Herazo, H. A., Massaro, F., Gu, M., et al. 2021, AJ, 161, 196

Plaga, R. 1995, Natur, 374, 430

Pshirkov, M. S., Tinyakov, P. G., & Urban, F. R. 2016, PhRvL, 116, 191302

Remillard, R. A., Tuohy, I. R., Brissenden, R. J. V., et al. 1989, ApJ, 345, 140

Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, NatMe, 17, 261

Woo, J.-H., Urry, C. M., van der Marel, R. P., Lira, P., & Maza, J. 2005, ApJ, 631, 762

Wood, M., Caputo, R., Charles, E., et al. 2017, in 35th ICRC, Fermipy: An open-source Python package for analysis of Fermi-LAT Data, 301 (Trieste: SISSA), 824