Constraining the Astrophysical Origin of Intergalactic Magnetic Fields

High-energy photons can produce electron–positron pairs upon interacting with the extragalactic background light. These pairs will in turn be deflected by the intergalactic magnetic field (IGMF), before possibly up-scattering photons of the cosmic microwave background, thereby initiating an electromagnetic cascade. The nonobservation of an excess of GeV photons and an extended halo around individual blazars due to this electromagnetic cascade can be used to constrain the properties of the IGMF. In this work, we use publicly available data of 1ES 0229+200 obtained with the Fermi Large Area Telescope and the High Energy Stereoscopic System to constrain cosmological MHD simulations of various magnetogenesis scenarios, and find that all models without a strong space-filling primordial component or overoptimistic dynamo amplifications can be excluded at the 95% confidence level. In fact, we find that the fraction of space filled by a strong IGMF has to be at least f ≳ 0.67, thus excluding most astrophysical production scenarios. Moreover, we set lower limits of B 0 > 5.1 × 10−15 G (B 0 > 1.0 × 10−14 G) for a space-filling primordial IGMF for a blazar activity time of Δt = 104 yr (Δt = 107 yr).


INTRODUCTION
The main hypothesis for the production of magnetic fields in astrophysical structures, from stars up to galaxy clusters, is that the magnetic fields are produced by amplifications of pre-existing magnetic fields by various dynamo effects or during gravitational collapse (Durrer & Neronov 2013).In order to explain the magnetic fields observed in galaxy clusters, one therefore usually presumes that the intergalactic magnetic field (IGMF) has a seed field that was present before galaxies were formed.This seed field is known as the primordial magnetic field, and will still be present in cosmic voids around its original strength, diluted only due to the expansion of the Universe and affected by (inverse) energy cascades.The production mechanism of the primordial field remains, however, still unknown.Thus, by constraining or detecting the IGMF, one can learn about unknown processes in the early Universe, e.g., deduce whether the IGMF was produced during inflation or during phase transitions (Durrer & Neronov 2013).Unfortunately, the primordial field will be weak, making it difficult to detect.Indeed, Faraday rotation measurements provide currently an upper limit on the order of 10 −9 G (Pshirkov et al. 2016).Alternatively, the IGMF could have an astrophysical origin (see below).
The IGMF can be probed by observing gamma-rays from blazars (Neronov & Semikoz 2009): High-energy photons can produce electron-positron pairs upon interacting with the extragalactic background light (EBL).The pairs will in turn be deflected by the IGMF, before they may up-scatter cosmic microwave background (CMB) photons via the inverse Compton scattering, thereby initiating an electromagnetic cascade.The deflection of the electrons and positrons leads to a time delayed arrival of cascade photons compared to the primary photons (Plaga 1995).In addition, the cascade photons will reach Earth slightly off-angle compared to the line-of-sight, leading to an extended halo around the primary point-source (Aharonian et al. 1994).The size of the halo and the time-delay depend on the strength of the IGMF: A strong IGMF implies a long time-delay and a wide cascade halo with low brightness.
The non-observation of an excess of GeV photons of individual blazars using Fermi LAT, first suggested by d' Avezac et al. (2007), has been used to set a lower limit on the IGMF of B ≳ 10 −15 -10 −16 G (e.g.Tavecchio et al. 2010;Neronov & Vovk 2010;Dolag et al. 2011;Podlesnyi et al. 2022).Likewise, using the nonobservation of cascade halos with imaging air Cherenkov telescopes have been used to rule out B ∼ 10 −16 using H.E.S.S. (Abramowski et al. 2014) and B ∼ 10 −14 using VERITAS (Archambault et al. 2017).The limits weaken, however, considerably for low blazar activity times (Dolag et al. 2011;Dermer et al. 2011) while becoming stronger for coherence lengths l B ≲ 0.1 Mpc.A recent and comprehensive review of existing limits on the IGMF is given by Alves Batista & Saveliev (2021).Currently, the most robust lower limit from blazar observations is B > 1.8 × 10 −14 G (B > 3.9 × 10 −14 G) for an activity time of 10 4 yr (10 7 yr), set by a combined analysis of Fermi LAT and H.E.S.S. blazar observations (Aharonian et al. 2023).
Even though current limits focus on a space-filling primordial field, it is likely that the IGMF has a complicated morphology induced by dynamo amplifications around galaxies and galaxy clusters, as well as outflows from magnetized galaxies.Indeed, an alternative production mechanism of the IGMF is that the magnetic fields are produced at late times by large-scale outflows from active magnetized galaxies (Durrer & Neronov 2013).This will allow for a miniscule magnetic field in the cosmic voids.Nevertheless, the lower limits from the blazar observations will in this case become tighter, and can be used to limit the filling fraction, f , of the IGMF (Dolag et al. 2011).
In this work, we make use of the recent Fermi LAT and H.E.S.S. observations of 1ES 0229+200 by Aharonian et al. (2023) to constrain the astrophysical origin of the IGMF.By comparing observations to Monte Carlo simulations performed using the photon propagation program ELMAG (Blytt et al. 2019;Kachelriess et al. 2012), we find a lower limit B 0 > 5.1 × 10 −15 G (B 0 > 1.0 × 10 −14 G) assuming a blazar activity time of ∆t = 10 4 yr (∆t = 10 7 yr), comparable to the results by Aharonian et al. (2023).The main emphasis in this work is, however, on the astrophysical origin of the IGMF.By using a top-hat distribution as a proxy for an astrophysical magnetic field, we find that the IGMF must have a filling fraction of at least f > 0.67, which excludes most astrophysical production scenarios.Furthermore, we perform the analysis on numerical models of magnetogenesis produced using cosmological MHD simulations (Vazza et al. 2017;Gheller & Vazza 2019).We find that none, even in the excessive galaxy formation feedback models, can magnetize the Universe in a way that is compatible with blazar data.Only models with a strong primordial magnetic field or over-optimistic dynamo amplifications cannot be excluded at the 95% confidence level (CL).

OBSERVATIONS OF 1ES 0229+200
The blazar 1ES 0229+200 is a high frequency peaked BL Lac object located at a redshift z = 0.14, with a hard spectrum extending to 10 TeV (Kaufmann et al. 2011).This makes it an ideal source to search for cascade photons and extended halos, and it has therefore been the primary source for setting bounds on the IGMF for a long time (e.g.Dolag et al. 2011;Abramowski et al. 2014).Recently, the observed Fermi LAT and H.E.S.S. spectra of 1ES 0229+200 were made publicly available (Aharonian et al. 2023) 1 , which we make use of in the present work.
The Fermi LAT observations comprise 11.5 yr of data leading to a detection of 1ES 0229+200 at 14.2 σ.The considered H.E.S.S. data consists of 144 hr (dead-time corrected) taken with the small 12 m diameter H.E.S.S. telescopes until 2016.The observations yield a source significance of 16.5 σ.Both Fermi LAT and H.E.S.S. spectral can be well described with a simple power law, ϕ(E) = N 0 (E/E 0 ) −Γ , with spectral index Γ = 1.77 ± 0.09 in the Fermi LAT energy band and Γ = 2.81 ± 0.11 in the H.E.S.S. energy band, respectively.Absorption of gamma rays on the EBL can account for the entire spectral softening.Further details on the analysis can be found in Aharonian et al. (2023).

IGMF MODELS
In this work, we make use of a selection of theoretical models of magnetogenesis produced using cosmological MHD simulations with the ENZO2 code (Bryan et al. 2014), known as the "Chronos++ suite" (Vazza et al. 2017(Vazza et al. , 2021) ) 3 .These explore numerous plausible scenarios for radiative gas cooling, star and AGN feedback and dynamo amplification, starting from the (non-radiative) evolution of large-scale magnetic fields with initial seed value of B 0 = 10 −9 G, B 0 = 10 −11 G or B 0 = 10 −18 G (comoving).All simulations start at redshift z = 38 and evolve at the constant comoving spatial resolution of 83.3 kpc per cell.More specific details of the tested astrophysical models are given in Appendix A. Portions of a line-of-sight magnetic field for five of the considered magnetogenesis models are visualized in Fig. 1: A spacefilling primordial field with seed value B 0 = 1 nG (P), a weak initial primordial field significantly enhanced by efficient dynamo amplification (DYN4), a model with extreme magnetization from astrophysical feedback (i.e. even beyond the likely feedback power of real low mass galaxies) (PCSF), and two with a weak initial primordial field and stellar feedback (CSFBH3 and CSF0).The models PCSF and DYN4 are, as we will see, in tension with current gamma-ray data, but cannot be excluded at the 95 % CL.Meanwhile, P cannot be excluded due to its high seed value, and CSFBH3 and CSF0 will be excluded due to their small filling fraction.
In addition to the MHD simulations, we consider some common simplified magnetic field models in order to draw some general conclusions about the allowed properties of the IGMF.For simplicity, we model the turbulent field as a simple domain-like magnetic field, where the field is split into homogeneous patches of size equal to the coherence length l B .The magnetic field strength in each patch has a constant strength B 0 , but a random direction.Since the cascade predictions remain independent of the coherence length for l B ≳ 0.1 Mpc, we fix for simplicity l B = 1 Mpc throughout this paper.As a proxy for an IGMF of a purely astrophysical origin, we consider a simple top-hat magnetic field, as shown in the bottom panel of Fig. 1.The magnetic field strength in the voids is B 0 = 0, while the magnetic field in the filaments is assumed to be sufficiently strong that the electrons and positrons in the electromagnetic cascade are completely isotropized.This model has then only two parameters: the width of filaments, W , and the filling fraction, f .In the particular example in Fig. 1, the filling fraction is f = W/D = 0.25, where D is the distance between the onset of two high magnetic field regions.

CASCADE PREDICTIONS
We simulate the electromagnetic cascade using the ELMAG 3.02 Monte Carlo code (Kachelriess et al. 2012;Blytt et al. 2019).For each magnetic field model, we inject photons with energies logarithmically distributed between 1 GeV and 100 TeV, and the results are stored in a multidimensional histogram depending on the primary energy of the initial photon, the secondary energy of the arriving photons at Earth, the time delay and the arrival angle.In this way, we can as a second step easily perform a fit to the experimental data by re-scaling the results by the source spectrum, as done in Aharonian et al. (2023).We consider the photons arriving within the 95% containment radius of the point-spread func- In contrast to the domain-like magnetic field, the effective filling fraction in the MHD simulations will depend on the chosen line-of-sight in the simulation volume.We therefore have to consider a sufficiently large amount of possible line-of-sights through the MHD simulation volume.To limit the necessary computational time, we have extracted 100 line-of-sight magnetic fields for each of the considered magnetic field models that we use in the analysis.Note that periodic boundary conditions are used since the size of the box in the MHD simulations, D = 83.3Mpc, is shorter than the distance to 1ES 0229+200, d ≈ 591 Mpc.Nevertheless, this will lead to a conservative estimation, since the cosmic variance of a larger simulation volume would lead to a smaller variation in the final cascade spectrum.Two examples of the cascade prediction are shown in Fig. 2 for visualization.

ANALYSIS AND RESULTS
In this work, we use the publicly available spectra from Aharonian et al. (2023) to constrain properties of the IGMF.As in Aharonian et al. (2023), we model the point source spectrum as where N 0 is the normalization, E 0 ≡ 1 TeV is the reference energy, Γ is the spectral index and E cut is the cut-off energy.The exponential cut-off at high energies will lead to conservative results in terms of the produced cascade emission, and will thus lead to conservative constraints on the IGMF.We construct a likelihood function by combining the likelihood profiles for each Fermi LAT energy bin4 with the χ 2 value of the H.E.S.S. data points.Using only the spectral information (and disregarding the spatial information of the extended halo) can lead, however, to a strong preference for a dominant cascade contribution.The reason is that, for a hard intrinsic spectral index and large E cut , the combination of the intrinsic emission and the cascade can perfectly describe the data (see also the discussion in Appendix C).However, this scenario would also predict a bright extended halo, which is not observed (Aharonian et al. 2023).Thus, to account partly for the non-observation of a cascade halo, we include an additional Gaussian prior on the spectral index, with a half-width5 σ = 0.05; the effect of the prior is detailed in Appendix C. The final likelihood function, i.e. the likelihood of observing the data D i with given model parameters, can be expressed as where B = (B 0 , f, . ..) refers to the parameters of the magnetic field model and θ = (N 0 , Γ, E cut ) contains the spectral parameters.Meanwhile, we fix the spectral index to the value obtained in the best fit case, Γ 0 = 1.64, which occurs for strong magnetic fields (i.e.no cascade halo).The likelihood fit and the various analyses are discussed in more depth in Appendix C. We now proceed to discuss the analyses performed.First, we set a lower limit on a space-filling cell-like IGMF with magnetic field strength B 0 .By virtue of Wilks' theorem, we can convert the likelihood function to a lower limit when 2(ln where the number α is the χ 2 -quantile.The subscript 'max' refers to the best fit magnetic field strength, which is achieved for a strong magnetic field (i.e.no detectable cascade photons).Since only one degree of freedom is considered, namely the magnetic field strength B 0 , we set a 95% CL upper limit with α = 2.71.We perform a grid search with B 0 = 10 −17 G, 10 −17.2 G, . . ., 10 −12.6 G, 10 −12.8 G and use cubic interpolation to estimate the lower limit.The final result is that B 0 > 5.1 × 10 −15 G is excluded at the 95% CL for an assumed blazar activity time of t max = 10 4 yr, while B 0 > 1.0 × 10 −14 G for t max = 10 7 yr.The results for other activity times are detailed in Appendix C.
Next, we set a lower limit to the filling fraction f of the top-hat IGMF model (see Fig. 1).Since B 0 in this case is assumed to be large and the cascade spectrum is independent of D ≳ 5 Mpc (see Appendix C.1), we are left with a single degree of freedom.As in the previous case, we perform a grid search with f = 0.00, 0.05, 0.10, . . ., 0.95, 1.00, and use cubic interpolation to estimate the lower limit.The result is that f < 0.67 is excluded with 95% confidence.
Finally, we compute the likelihood function for 100 line-of-sight magnetic fields for each of the models discussed in Sec. 3. In this case, however, the profile likelihood function cannot be directly applied to exclude the IGMF models since the exact number of degrees of freedom is unknown (we are not varying a single parameter like the magnetic field strength or the filling factor).The correct approach would be through Monte Carlo simulations, but it is quite computationally expensive and the instrumental response function (IRF) of H.E.S.S. is not publicly available.Instead, we note that the magnetic field distribution can be approximated by a top-hat distribution with a filling fraction, f , and mag- netic field strength, B, whose effects on the likelihood are furthermore correlated.The exact shape of the cascade spectrum will in general depend on the exact magnetic field distribution, in particular through the typical size between structures, D, and the correlation length of the magnetic field, l B .However, the cascade spectrum remains constant for D ≳ 5 Mpc (see Appendix C.1) and l B ≳ 100 kpc, both of which are the case in the considered IGMF models.Thus, we consider a magnetic field model to be excluded by at least 95% confidence if −2(ln L − ln L max ) > 4.61, corresponding to a χ 2 distribution with two degrees of freedom.The average and the 95% quantile of the likelihood estimates for each of the 100 line-of-sights magnetic fields are listed in Tab. 1.The estimated filling fraction of the models with |B| > 10 −15 G is also listed (see Appendix A.2).
We rule out a model if the 95% quantile of the delta log likelihood value is above 4.61.Of the considered MHD simulations without a strong primordial component, only DYN4 and PCSF survive our test.

DISCUSSION AND CONCLUSION
We have used the publicly available data of 1ES 0229+200 from Fermi LAT and H.E.S.S. (Aharonian et al. 2023) to set limits on the properties of the IGMF.By simulating the electromagnetic cascade using ELMAG, we find that a space-filling IGMF with B 0 > 5.1 × 10 −15 G (B 0 > 1.0 × 10 −14 G) for a blazar activity time of ∆t = 10 4 yr (∆t = 10 7 yr) is excluded by 95% CL.Note, however, that the source has only been observed for ∼ 20 years, which entails a lower limit on the activity time.In the conservative case ∆t = 20 yr, we obtain B 0 > 2.4 × 10 −16 G. Thanks to our prior on the spectral index, we are able to account (partially) for the spatial information of Fermi LAT and H.E.S.S.. Thus, our limits are among the most stringent limits on the IGMF, only slightly weaker than those obtained with full spatial information (Aharonian et al. 2023).Note that the fact that our results are perfectly compatible to what was found in the H.E.S.S. analysis by Aharonian et al. (2023) can be considered as an additional motivation for the prior.
The focus in this work has been on constraining an astrophysical origin of the IGMF.We have set new limits on the filling fraction f > 0.67 using a top-hat distribution of the magnetic field.Note that our result can be considered as conservative since we assume a strong magnetic field in the filaments: Decreasing B 0 would lead to a greater cascade spectrum, implying that a larger filling f is needed to avoid over-producing the observed spectrum.It is worth mentioning that this result is similar to the results obtained by Dolag et al. (2011) based on old upper limits of the spectrum of 1ES 0229+200 from Fermi LAT.However, the limits obtained in this work are statistically more robust since the spectral index, cut-off energy and normalization are included as nuisance parameters.The limit can in principle be used to get approximate bounds on the more realistic magnetic field models that we have considered in this work.Indeed, from the estimated filling fraction of the MHD models in Tab. 1, one can foresee that all purely astrophysical models we examined (CSFBH3, CSF0, CSF2 and CSF3) can be excluded with high significance, while DYN4 survives the likelihood test.Meanwhile, PCSF may survive the test due to cosmic variance, but will be in tension with the data.
We performed a likelihood analysis of a variety of simulated magnetogenesis scenarios.Due to the cosmic variance in the MHD simulations, the properties of the magnetic field (e.g., effective filling fraction) depend on the exact line of sight considered in the simulation volume.Therefore, we conducted the likelihood test on 100 random line-of-sights for each model, and consider the model excluded at the 95% CL if its 95% quantile has a log likelihood value −2∆ log L > 4.61, corresponding to a 95% exclusion of a χ 2 distribution with two degrees of freedom.Only the models DYN4 and PCSF cannot be excluded at the 95% CL.Both may likely be excluded in an improved analysis including more sources and complete spatial information from H.E.S.S., which currently is not publicly available.However, DYN4 and PCSF can be rejected as viable possibilities based on other complementary tests: The high efficiency dynamo amplification model DYN4 has been shown to yield a too large amount of Faraday rotation for polarized radio sources shining through the cosmic web, at variance with real LOFAR observations (e.g., O'Sullivan et al. 2020;Vazza et al. 2021).The PCSF model instead employs by construction a too large feedback efficiency for low mass halos, yielding inconsistent scaling relations for halos and star formation histories (Vazza et al. 2017) and it has been employed here just to gauge the maximum conceivable contribution from astrophysical seeding.
A few caveats of our results should be mentioned (see also the discussions in Aharonian et al. (2023)).Perhaps most importantly, the limits depend strongly on the assumed blazar activity time, reducing by about an order of magnitude for every second decade in activity time; the detailed dependence is given in Appendix C.Moreover, it has been argued that the electrons and positrons may predominantly lose their energy through plasma instabilities if the blazar activity time is ∆t ≳ O(10 2 ) yr (Broderick et al. 2012).If confirmed, it will invalidate the limits based on cascade emissions from TeV blazars.Safe to say, however, the effect is still debated (Rafighi et al. 2017;Alves Batista & Saveliev 2021).
In conclusion, with our work we have significantly narrowed down the possibility that the IGMF can be of pure astrophysical origin, because the necessary fraction of the space which must be filled by a strong magnetic field appears impossible to be obtained by all astrophysical scenarios for the magnetization of the cosmic web we tested.The upcoming Cherenkov Telescope Array (CTA) will improve the constraints further (Meyer et al. 2016;Abdalla et al. 2021).
The authors would like to thank M. Kachelrieß for fruitful discussions and for providing valuable comments on the manuscript.J.T. would like to express gratitude for the hospitality at the University of Agder (UiA).M.M. acknowledges the support from the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany's Excellence Strategy -EXC 2121 "Quantum Universe" -390833306 and from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program Grant agreement No. 948689 (Ax-ionDM).F.V. acknowledges financial support from the Cariplo "BREAKTHRU" funds Rif: 2022-2088 CUP J33C22004310003.In this work, we used the ENZO code (http://enzo-project.org), the product of a collaborative effort of scientists at many universities and national laboratories.Their commitment to open science has helped make this work possible.F.V. gratefully acknowledges the Swiss National Supercomputing Centre (CSCS), Switzerland, for the access to the Piz Daint, under the allotted project "s1096".• sub-grid dynamo: the small-scale dynamo amplification of magnetic fields in the simulation is computed at run-time via sub-grid modelling, based on the estimated rate of dissipation of solenoidal turbulence.These models attempt to bracket the possible residual amount of dynamo amplification which can be attained, on energy grounds, in the cosmic structures, but may be lost because of a lack of resolution and/or for the onset of small-scale plasma instabilities that cannot be captured from first principles in the simple hydro/MHD model.
• radiative cooling: standard equilibrium gas cooling implemented in ENZO, based on the chemical distribution of 7 main elements (HI, HII, HeI, HeII, HeIII, free electrons and "metals").
• star formation: star formation criterion (Kravtsov 2003) included in ENZO.Here, new star particles are formed on the fly if the local gas density exceeds a given threshold, n * , with a dynamical timescale, t * .The minimum star mass in the considered simulations is set to m * = 10 7 M ⊙ .Whenever (i) n ≥ n * , (ii) the gas is contracting (∇ • v < 0), (iii) the local cooling time is smaller than the dynamical timescale (t cool ≤ t * ), and (iv) the baryonic mass is larger than the minimum star mass (m b = ρ∆x 3 ≥ m * ), stars with a mass m * = m b ∆t/t * are set to form.This recipe was introduced to reproduce the observed Kennicutt's law (Kennicutt 1998).
• feedback from star formation: stars return thermal energy, gas and metals back into the gas phase.The feedback from the star forming process (e.g.supernova explosions) depends on the assumed fractions of energy/momentum/mass ejected for each formed star particle, E SN = ϵ SF m * c 2 .The feedback energy is assumed to be released entirely in thermal form (i.e.hot supernovae-driven winds), which typically turns into the formation of pressure-driven winds around the halos (v wind ∼ 10 − 10 2 km/s).For an overview of the performances and limitations of the recipes for star formation and feedback in ENZO we refer the reader to the recent survey of models by Skory et al. (2013) and Hlavacek-Larrondo et al. (2022).The limited spatial and mass resolution is expected to quench the formation of ≤ 10 8 M ⊙ galaxies, implying an important contributor to the cosmic star formation rate (e.g.Genel et al. 2014).However, the contribution to the chemical, thermal and magnetic enrichment of the intergalactic medium by galactic outflows is predicted to be dominated by ∼ 10 10 − 10 11 M ⊙ (e.g.Samui et al. 2018) galaxies, which are formed in the runs considered here.
• supermassive black holes (SMBHs): black hole (BH) seed particles are injected at z = 4 at the centre of all ≥ 10 9 M ⊙ halos identified in the simulation (Kim et al. 2011).For each seed particle, a fixed initial mass of M BH = 10 4 M ⊙ is assumed.Each BH particle accretes gas mass according to the Bondi-Hoyle rate, assuming all gas gets accreted at the fixed temperature of 3 × 10 5 K (because a precise estimate of gas temperature on the accretion radius is not feasible with the given resolution) and introducing a numerical boost to the measured accretion rate by α Bondi = 100 to overcome the effect of the limited gas density which can be resolved around the BH particles.
• supermassive black hole (SMBH) feedback: as an extra thermal energy output from each BH particle, they release thermal feedback on the surrounding gas, assuming an efficiency ϵ BH = ∆M c 2 ∆t/E BH between the accreted mass, ∆M , and the feedback energy, E BH .In the PCSF model variation, a "maximal" SMBH feedback scheme is used instead, in which every SMBH is forced to release the same fixed amount of feedback thermal energy, E BH = 5 × 10 59 erg, regardless of the mass accretion rate.This model is extreme, in the sense that it overpredicts the typical feedback energy in the lowest mass galaxies, up to ∼ 100-1000 factors in power; however, it is meant to bracket the maximum possible effect of astrophysical magnetization in the simulation volume.
• magnetic field seeding: For the injection of magnetic fields by stellar feedback, the generation of magnetic dipoles during each feedback episode is introduced in ENZO.The dipoles are randomly oriented along one of the coordinate axes of the simulation, and their total energy scales with the thermal feedback energy as E b,SN = ϵ SF,b •E SN , typically with ϵ SF,b ∼ 1-10%.Also for the release of magnetic energy from SMBH, magnetic fields from the BH region are injected (in the form of magnetic dipoles) assuming E b,AGN = ϵ b,AGN E BH , where values in the range ϵ b,AGN = 1-10% of the energy conversion into magnetic field dipoles at the jet base were explored.
While the fixed spatial resolution of the runs is not good enough to model properly the galaxy formation processes with sufficient detail (even massive galaxies are resolved within a few cells at most), it was shown by Vazza et al. (2017) how the star-formation recipes are specifically calibrated to reproduce properly the cosmic star formation history in the volume.Moreover, the combination of cooling and stellar/AGN feedback has been tuned so that the scaling relations of galaxy clusters and groups are well reproduced.The models also describe the innermost density, temperature, and entropy profiles of galaxy clusters fairly well in comparison to higher resolution simulations.In summary, previous work suggest that the feedback recipes implemented in the considered simulations can effectively mimic the large-scale effect and energetics of galaxy formation recipes at a higher resolution, hence they can be effectively used to model the large-scale effects of thermal and magnetic feedback on the surrounding distribution of filaments.

A.2. Properties of the magnetic fields
In a typical IGMF of astrophysical origin, the propagation of the electrons and positrons in an electromagnetic cascade can be split into two distinct regimes: When the electron Larmor radius, r L , is much smaller than the inverse Compton scattering length, λ IC , they move ballistically.This is the case in the voids, where the magnetic field is weak or potentially non-existent.On the other hand, when r L ≫ λ IC , they will be isotropized and contribute to the diffuse gamma-ray background.This is typically the case in the filaments.The typical transition between the two regimes is where the Larmor radius of the electrons is of the same order as the inverse Compton scattering length.That is, r L ∼ λ IC , which in turn leads to E e /TeV ∼ B/10 −15 G.In order to get a more intuitive understanding of the magnetic field models and their effect on the cascade spectrum, we compare them with a simple top-hat distribution of the line-of-sight magnetic field, wherein the filaments have a fixed size and mutual distance.Furthermore, it is assumed that B ∼ 10 −9 G in the filaments, so that r L ≪ λ IC .Thus, the model can be characterized by two numbers: the distance between filaments, D, and the filling fraction, f .We extract the numbers from the line-of-sights by defining filaments as regions with magnetic field strength larger than 10 −15 G. Furthermore, we demand that a void has to be larger than 333 kpc (i.e., four times the simulation grid).The distributions of D and f for the considered IGMF models are summarized in Tab. 3 and Fig. 3. Table 3. Summary of the average effective filling fraction, f , distance between magnetic field structures, D, and width of magnetic field structures, W (cf. Fig. 1) of the considered MHD models listed in Tab. 2.

B. DETAILS ON THE CASCADE SIMULATIONS WITH ELMAG
For each magnetic field model, we inject 10 5 -107 photons6 logarithmically randomly distributed between 1 GeV and 100 TeV.The results are stored as a four-dimensional array depending on the time-delay, ∆t, the primary energy of the injected photon, E i , the energy of the arriving photon, E f , and the arrival angle, θ.For the energies, 10 bins per decade are used.Meanwhile, for the time-delay, one bin per decade between ∆t = 1 and 10 7 years is used, in addition to ∆t < 1 and ∆t > 10 7 years.Finally, for the arrival direction, 500 bins are linearly distributed between θ = 0 • and 1.5 • .
Note that some changes to ELMAG 3.02 were needed for this work (see also Kachelriess & Tjemsland 2022, 2024): • Inclusion of a domain-like magnetic field • Inclusion of line-of-sight magnetic fields provided by a file; The field strength is set to the projection of the magnetic field along the line of sight • Export of results as multidimensional histograms • Inclusion of a top-hat magnetic field • Fine-tuning the ordinary differential equation solver, such that one can consider magnetic fields stronger than ∼ 10 −13 G C. DETAILS ON THE FIT We construct the log-likelihood estimate using the publicly available spectral information of the blazar 1ES 0229+200 from Fermi LAT and H.E.S.S. (Aharonian et al. 2023), and minimize it with a spectral reweighting (Ackermann et al. 2018) of the multidimensional histogram discussed in Appendix B using using the differential evolution algorithm in scipy.optimize.The resulting likelihood profile is plotted with dashed lines in the left panel of Fig. 4 for blazar activity times 10 4 yr (blue) and 107 yr (orange) for a domain-like magnetic field.The dependence of the lower limit on the activity time is shown in the right panel of Fig. 4. For comparison, we also plot the limits with a wider prior σ = 0.1 (open circles).The limits as a function of the activity times are well reproduced 7 as B 0 ∆t/yr, becoming constant at t 0 .The best fits8 , shown in Fig. 4, are {t 0 = 3.38 × 10 4 yr, B 0 = 5.46 × 10 −17 G} with the prior on the spectral index (σ = 0.05), and {t 0 = 7.08 × 10 4 yr, B 0 = 1.74 × 10 −18 G} without the prior.
We have tested that the results are consistent with those obtained using a Gaussian turbulent field, a rescaling of the pure primordial magnetogenesis scenario (P), and the 2.5D implementation in ELMAG.
In the likelihood profile (left panel in Fig. 4), it is apparent that there is an allowed region around B 0 ∼ 10 −15 G, which is due to a strong preference for a large cascade contribution in the fit, i.e. large cut-off energy and small spectral index.This is well visualized in the spectra plotted in Fig. 5.However, a significant cascade contribution has already been excluded due to the non-observation of the gamma-ray halo (Aharonian et al. 2023).To accommodate for this effect, we include a Gaussian prior on the spectral index with half-width σ = 0.05, a conservative prior given that the spectral index only vary by 0.002 for ∆t = 10 7 yr in Aharonian et al. (2023).In this case, the minimum vanishes, and the results are on par with those of Aharonian et al. (2023).The results are shown in Fig. 4 in solid lines.) of 1ES 0229+200: The source spectrum (brown dotted line) is attenuated by the EBL, so that an attenuated primary spectrum (red dashed line) is observed on Earth.In addition, the electromagnetic cascade may induce a spectrum of cascade photons (purple dotted dashed line) at a reduced energy.The source spectrum is adjusted such that the total spectrum (green solid line) minimizes the likelihood estimate of the measured data.The results are shown for B0 = 10 −13 G (left column) and B0 = 10 −15 G (left column), as well as without the Gaussian prior on the spectral index (top row) and with the Gaussian prior (bottom row).

Figure 1 .
Figure1.Visualization of the magnetic field strength of some of the magnetic field models considered: P, DYN4, PCSF, CSFBH3, CSF0 and a top-hat distribution.The models are described in more detail in the main text.

Figure 3 .
Figure3.Distributions of the degree of filling, f , distance between filaments, D, and size of filaments W of the 100 line-of-sights that are considered in the likelihood analysis of the MHD models.The different colours correspond to different models.

Figure 5 .
Figure5.Visualization of the contributions to the measured spectrum (orange error bars for Fermi LAT, blue for H.E.S.S.) of 1ES 0229+200: The source spectrum (brown dotted line) is attenuated by the EBL, so that an attenuated primary spectrum (red dashed line) is observed on Earth.In addition, the electromagnetic cascade may induce a spectrum of cascade photons (purple dotted dashed line) at a reduced energy.The source spectrum is adjusted such that the total spectrum (green solid line) minimizes the likelihood estimate of the measured data.The results are shown for B0 = 10 −13 G (left column) and B0 = 10 −15 G (left column), as well as without the Gaussian prior on the spectral index (top row) and with the Gaussian prior (bottom row).

Figure 6 .Figure 7 .
Figure6.Likelihood profile for the filling fraction for the top-hat distributed magnetic field.The result is shown both with (solid lines) and without (dashed lines) the prior on the spectral index, as well as for various values of the size of the magnetic field clusters, D. Values of f giving a likelihood value above the gray dotted line are excluded.

Table 1 .
Figure2.Example plot of the cascade prediction for the MHD models DYN4 (left) and PCSF (right).Each of the 100 lines correspond to a random line-of-sight from the simulation volume of the MHD simulation.The source spectrum is adjusted so that the final spectrum minimizes the maximum likelihood estimate of the experimental measurements by Fermi LAT (squares) and H.E.S.S (circles), as explained in Sec. 5. Table containing the estimated filling fraction f , mean log likelihood estimate and the 95% quantile of the log likelihood estimates for the 100 considered line-of-sights of the MHD simulations.