Magnetization of the intergalactic medium in the IllustrisTNG simulations: the importance of extended, outflow-driven bubbles

We study the effects of galaxy formation physics on the magnetization of the intergalactic medium (IGM) using the IllustrisTNG simulations. We demonstrate that large-scale regions affected by the outflows from galaxies and clusters contain magnetic fields that are several orders of magnitude stronger than in unaffected regions with the same electron density. Moreover, like magnetic fields amplified inside galaxies, these magnetic fields do not depend on the primordial seed, i.e. the adopted initial conditions for magnetic field strength. We study the volume filling fraction of these strong field regions and their occurrence in random lines of sight. As a first application, we use these results to put bounds on the photon-axion conversion from spectral distortion of the CMB. As photon-axion coupling grows with energy, stronger constraints could potentially be obtained using data on the propagation of gamma-ray photons through the IGM. Finally, we also briefly discuss potential applications of our results to the Faraday Rotation measurements.


INTRODUCTION
The Universe is magnetized on all scales -from planets and stars to galaxy clusters and beyond. Magnetic fields affect the propagation of charged particles and therefore play an important role in the physics of the Earth's atmosphere, the Sun and solar system, galaxies, and so on. However, the origin of magnetic fields in galaxies and, especially, galaxy clusters in the low-redshift Universe remains an important open problem (see e.g. Durrer & Neronov 2013 for a review).
Indeed, there are two stages of magnetogenesis -the first stage is the generation of weak seed magnetic fields and the second stage is their subsequent evolution during structure formation. Seed fields could have either primordial (produced before recombination) or astrophysical origin (produced after the formation of the first stars), see e.g. Subramanian (2016) for a review. Primordial seed fields fill aramburo@strw.leidenuniv.nl † kyrylo.bondarenko@cern.ch ‡ boyarsky@lorentz.leidenuniv.nl § dnelson@uni-heidelberg.de ¶ pillepich@mpia-hd.mpg.de anastasia.sokolenko@oeaw.ac.at the whole universe, although they are not necessarily constant, their correlation length depends on the production mechanism and the epoch when they were produced. An astrophysical seed magnetic field is generated by Biermann battery-type mechanisms when the curl of the electric field is created by non-parallel gradients of density and temperature, giving rise to a magnetic field via Faraday induction (Subramanian 2016). Such a mechanism could generate magnetic fields during star formation, reionization (Subramanian et al. 1994;Gnedin et al. 2000), or even during the later collapse of galaxies and halos in cosmological shocks (Kulsrud et al. 1997). For example, Vazza et al. (2017) simulate 25 different scenarios of initial magnetogenesis, both primordial and astrophysical.
In the second stage, magnetic fields evolve with the expansion of the Universe and structure formation. Outside structures, these fields dilute approximately as ∼ a −2 (Durrer & Neronov 2013), where a is the scale factor. In the regions where dense structures form, magnetic fields are adiabatically compressed and, moreover, strongly amplified (up to several µG Pakmor et al. 2014a;Rieder & Teyssier 2017) by different dynamo mechanisms driven by the baryonic physics of galaxy formation (Parker 1955;Ruzmaikin et al. 1988;Kulsrud 1999 Butsky et al. 2017;Steinwandel et al. 2019;Martin-Alvarez et al. 2018;Vazza et al. 2018). Amplification in filaments also occurs via shear flows (Birk et al. 1999;Dolag et al. 1999;Dolag et al. 2005). As a result, magnetic fields in galaxies and collapsed structures, amplified by many orders of magnitude by gravomagnetohydrodynamics (MHD) dynamos, "forget" the properties of the initial magnetic fields (see e.g. Pakmor et al. 2014a;Marinacci et al. 2015;Pillepich et al. 2018a, for cosmological simulations of galaxies). To the contrary, magnetic fields that are far from structures are much closer to the diluted initial fields and could therefore be used to infer information about the properties and the origin of the initial fields.
On the observational side, cosmic magnetic fields are relatively well studied in virialized objects like galaxies and galaxy clusters. A powerful method to detect these fields is Faraday Rotation Measure (RM) see e.g. Brentjens & de Bruyn (2005) and references therein. With the current generation of instruments, this method is efficient for magnetic fields with the strength of B O(1) nG (Durrer & Neronov 2013). Such magnetic fields exist mainly in the dense centers of collapsed structures -galaxies, galaxy groups, and clusters (e.g. Carilli & Taylor 2002;Laing et al. 2008;Beck 2015;van Weeren et al. 2019). However, these objects fill only a small fraction of the volume of the Universe.
Empirical constraints on the cosmic magnetic fields outside galaxies and clusters remain difficult. Attempts to measure magnetic fields in filaments by Faraday Rotation Measure with LO-FAR (O'Sullivan et al. 2020) place only an upper bound of a few nG, on Mpc scales, consistent with other works (Ravi et al. 2016;Vernstrom et al. 2019;Blasi et al. 1999;Pshirkov et al. 2016;Hackstein et al. 2016;Bray & Scaife 2018). In addition, a lower bound can be obtained with high-energy gamma-ray data (Neronov & Vovk 2010;Dermer et al. 2011;Tavecchio et al. 2011;Dolag et al. 2011;Taylor et al. 2011). This lower bound is as weak as B 10 −17 G on Mpc scales. In the future, the upper bound can be improved with next-generation radio telescopes such as the Square Kilometer Array (SKA) (Carilli & Rawlings 2004), while the lower bound may be improved by high-energy observatories including CTA (Acharya et al. 2018), which is expected to start obtaining data soon. For now, the observational uncertainty in the properties of the intergalactic magnetic fields (IGMF) outside galactic halos is large.
From the theoretical as well as observational perspective our knowledge on the magnetic fields outside galaxies and clusters is rather limited. At the same time, magnetic fields in the large volumes of the intergalactic medium (IGM) can play a profound role in many important problems in physics. For example, magnetic fields can strongly affect the propagation of light and the spectra of various astrophysical sources (see e.g. Brentjens & de Bruyn 2005;Neronov & Vovk 2010) as well as propagation of cosmic rays (Alves Batista et al. 2017). If the magnetic field along the lines of sight to certain classes of sources were better constrained, these effects (see examples below) can give us important insight into fundamental physics.
When light propagates towards us from remote sources, most of the intervening pathlength is typically not in virialized objects, as these occupy only a tiny fraction of the Universe by volume. Rather, photons propagate through the IGM -the space between dark matter halos occupied by less dense regions including cosmic voids, sheets, and filamentary structures. Both the free electron number density and magnetic field strength in the IGM evolve with time, and theoretical modeling of this evolution, through the epoch of reionization and down to the present day, remains a challenge. Therefore, to probe the effects of cosmic magnetic fields on light propagation it is not enough to model magnetic fields only within halos.
One regime where magnetic fields in the IGM play a crucial role is gamma-ray astronomy. The photons from high-energy gamma-ray sources create electron-positron pairs interacting with the extra-galactic background light. These charged particles can then emit secondary gamma-ray photons interacting with the CMB via the inverse Compton effect. The presence of a magnetic field results in a deviation of the charged particles and, therefore, in a change in the morphology of the signal (see e.g. Neronov & Semikoz 2009 for a more detailed discussion). 1 Magnetic fields can also affect light propagation in the presence of new, as of yet unobserved particles that are not included in the Standard Model of particle physics. A famous example is an axion or axion-like particle (ALP), initially introduced to explain why CP violation in QCD is so tiny (Weinberg 1978;Wilczek 1978). Axions have been theorized to play the role of dark matter (Preskill et al. 1983). Photons can be converted into ALPs, but only when they pass through magnetized regions of the Universe. In this case, the conversion probability depends sensitively on the strength of the magnetic field (Sikivie 1983;Raffelt & Stodolsky 1988).
In this paper, we use the IllustrisTNG suite of cosmological simulations (see Section 2 for details), as well as additional variation runs performed with different values of the initial magnetic fields and with and without feedback, to study the regions of the IGM that could be affected by galactic outflows. The starting point of our investigation is the idea that the strong magnetic fields generated deep within dark matter halos can affect and extend to much larger volumes, as baryonic outflows, caused by strong feedback processes, eject magnetized gas to regions extending far beyond halo scales (and not just beyond galactic scales, as studied e.g. in Marinacci et al. 2018;Steinwandel et al. 2020;Pakmor et al. 2020;Dubois & Teyssier 2010;Martin-Alvarez et al. 2020).
In this work, we concentrate our attention on the regions of the IGM affected by galactic outflows, their origin, and their impact on magnetic fields in the IGM. We discuss the properties and strength of the IGMF predicted in IllustrisTNG (Section 2), presenting our results in a form that may be used for different applications. We show that the magnetic fields affected by galactic outflows depend more on the MHD processes occurring within galaxies rather than on the primordial magnetic field seeds. In particular, we show that the predicted magnetic field strength in the regions affected by outflows is similar for runs spanning orders of magnitude different strength of the initial magnetic fields (Section 3). At the same time, magnetic fields in these regions are orders of magnitude larger than in other regions of the IGM with similar matter density. As a first example application, we apply our findings to constrain ALPs (Section 4) and summarize our results (Section 5).

Simulations
IllustrisTNG (TNG) is a suite of large-volume cosmological gravomagnetohydrodynamic simulations (Nelson et al. 2018;Springel et al. 2018;Pillepich et al. 2018b;Naiman et al. 2018; Marinacci . Each evolves initial conditions from z = 127 to the present time, following the evolution of gas, stars, and black holes (baryons), together with dark matter. The TNG simulations use the moving-mesh AREPO code (Springel 2010) to solve the coupled equations of self-gravity and ideal MHD (Pakmor et al. 2011;Pakmor & Springel 2013), and adopt cosmological parameters consistent with Planck 2015 (Planck Collaboration et al. 2016). The simulations include a comprehensive galaxy formation model incorporating astrophysical processes such as gas metal-line cooling and heating, star formation, stellar evolution, and heavy element enrichment, supermassive black hole growth, AGN feedback, and galactic winds launched by supernovae (Weinberger et al. 2017;Pillepich et al. 2018a). The TNG project currently spans three different volumes, TNG50, TNG100, and TNG300, each run with several different numerical resolutions. In this work, we mainly use the publicly available TNG100-1 simulation (Nelson et al. 2019a), the highest resolution run of TNG100, with a box side-length of L ∼ 100 cMpc (comoving Mpc). Containing 1820 3 dark matter particles and an equal number of gas cells, it has a mass resolution of mbar = 1.4 × 10 6 M , and mDM = 7.5 × 10 6 M , respectively. From now on, we refer to such simulation as TNG100.

Galaxy formation model
As we are particularly interested in the role of galactic-scale outflows in producing extended regions of high magnetic field strength, we describe the feedback physics briefly. A supermassive black hole (SMBH) is created in all dark matter halos which exceed a total mass of ∼ 7 × 10 10 M , by placing a SMBH at the potential minimum with an initial mass of ∼ 10 6 M . These black holes sub-sequently grow via binary mergers with each other during galaxy collisions, and via smooth gas accretion from the surrounding environment. Black hole accretion is calculated using the Bondi-Hoyle-Lyttleton assumption (Weinberger et al. 2017), which depends on the black hole mass, local gas density, and relative velocity between the black hole and its surroundings. The accretion rate onto SMBHs is limited to the Eddington rate. To model energetic feedback from SMBHs, a small fraction of the rest mass of accreted matter is available to be deposited back into the locally surrounding gas. This energy is injected in a dual-mode model, depending on this accretion rate: one mode is for the high-accretion state (above ∼ 10 percent of the Eddington rate), while the second, low-accretion state feedback mode, operates for accretion rates roughly below this value.
At high accretion rates, energy is deposited in a continuous manner, by thermally heating gas. At low accretion rates, kinetic energy is injected in a discrete rather than continuous fashion, such that feedback events occur once enough energy accumulates (see Weinberger et al. 2017 for additional details). In this mode each injection event is modeled as a high-velocity kinetic wind, which is oriented randomly for each event, producing a timeaverage isotropic energy injection. This sub-resolution model is based on theoretical expectations for low accretion rate black holes, i.e. below one percent of Eddington, which develop radiatively inefficient flows and thereby convert gravitational binding energy into a non-relativistic wind (Blandford & Begelman 1999;Yuan & Narayan 2014). In the TNG model, it is this population of low luminosity, slowly accreting SMBHs which drive the most powerful outflows (Nelson et al. 2019b). This occurs above a characteristic galaxy stellar mass (dark matter halo mass) of ∼ 10 10.5 M (∼ 10 12 M ), corresponding to the onset of quenching in the seed field values but with the same underlying TNG galaxy formation model, over 25 cMpc/h a-side volumes: 10 −12 cG (blue dashed), 10 −14 cG (red solid, the fiducial choice of the TNG simulation adopted throughout the paper), 10 −16 cG (green dotted), 10 −25 cG (magenta dot-dashed), where cG are comoving Gauss. While the magnitude |B| depends on the assumed initial field strength for very low field strengths, corresponding to underdense/void-like regions of the simulated universe (see Fig. 3), the magnetic field strength within collapsed structures at log (|B|/G) −10 is largely unchanged, due to the rapid amplification processes which effectively erase knowledge of the primordial seed fields explored in this work.

Mass weighted fraction
Volume weighted fraction  Fig. 2. We excluded values for two fractions for the largest seed magnetic field 10 −12 cG because they are significantly contaminated by the seed field.
galaxy population (Weinberger et al. 2018;Donnari et al. 2020). Lower mass galaxies can also produce strong outflows, via a model for supernovae-driven winds originating from SNII explosions associated with ongoing star formation (see Pillepich et al. 2018a for more details). In general, these outflows are slower and do not escape to distances as large as black hole-driven outflows (Nelson et al. 2019b).

Initial conditions and model variations
When specifying the initial conditions for the gas, in addition to small density and velocity perturbations required to realize the chosen cosmological constraints, TNG must also specify the initial conditions for the magnetic fields, which are given in terms of the strength and direction of the initial magnetic fields, which can be very different for the different production mechanisms, i.e. astrophysical and primordial (Garaldi et al. 2020). All TNG simulations to date have been run with the same configuration of the initial magnetic field, which is assumed to be a constant volume-filling field (which is an approximation for a primordial, i.e. inflational magnetogenesis, and which is the common practice in such simulations. In this work we mainly use the TNG100 box with the initial strength B0 = 10 −14 cG (comoving Gauss). At z = 127 this corresponds to a physical field strength of B0 = 0.16 nG.
In addition, a large number of 'TNG model variation' simulations have been run, each changing a single parameter value or model choice to assess its importance in the fiducial TNG galaxy formation model (Pillepich et al. 2018a). These simulate smaller 25 cMpc/h volumes, each realized at three resolution levels equivalent to those available for TNG100 itself. In this work, we primarily use four simulations which vary the initial magnetic field strength, adopting B0 = 10 −25 cG, 10 −16 cG, 10 −16 cG (the fiducial choice of the TNG100 flagship run) and 10 −12 cG. We also inspect the outcome of three other runs that, compared to the TNG fiducial model, include no SMBH feedback, no SMBH kinetic feedback, and no feedback of any type, respectively.

Analysis methodology
We analyze these simulations from 0 z 6 and measure the density of free electrons and magnetic field strength, avoiding z > 6 where these quantities become uncertain during the epoch of reionization (Heinrich et al. 2017;Aghanim et al. 2020). The electron number density is calculated from the helium and hydrogen number densities and their ionization states. We use the spectral synthesis code CLOUDY V17.01 (Ferland et al. 2017) to calculate the ion fractions of hydrogen and helium for gas exposed to a redshift-dependent UV background (Haardt & Madau 2012). We neglect the contribution of ionization states from elements heavier than He, as well as molecular gas phases.
The magnetic field and electron number density are smoothed onto a grid with (20 ckpc) 3 voxels. 2 One thousand lines of sight (LOS) are generated for each available snapshot with random orientation within the simulated volume, and we use these throughout as our fiducial set of sightlines. In Fig. 1 we show an example of the magnitude of the magnetic field, metallicity, and electron number density image of a given region together with a potential LOS. The regions of magnetic field enhancement extend beyond gravitationally collapsed structures (i.e. dark matter halos), c.f. also Fig. 7. They also extend, in some cases, to substantially larger distances, due to the ejection of magnetic fields in supernovae and black holedriven galactic outflows (Nelson et al. 2019b).

STATISTICAL PROPERTIES OF SIMULATED MAGNETIC FIELDS
We start by examining the three 25 cMpc/h variation boxes with different initial magnetic field seed conditions. In particular, we aim to determine to what extent the predicted IGMF depends on the value of the initial seed. In Fig. 2 (left panel) we show the mass-weighted distribution of the magnetic field for three different choices of the initial field. These global distributions across all gas cells contain two clear peaks -a low-B peak that has its center at the value of the initial field, and a strong-B peak, that are very similar for all four values of the initial conditions. It is important to keep in mind that the number of gas cells with a given value of the field does not represent the fraction of volume occupied by such a field. As the simulation is spatially adaptive and maintains a constant mass resolution, gas cells accumulate inside high-density regions. We therefore also show the volume-weighted distribution of magnetic field strength in Fig. 2 (right panel). Both cases show the same picture: that the high-B component of the distribution is weakly sensitive to the magnitude of the initial field. To characterize difference in strong-B peaks numerically we calculated mass and volume weighted fractions within this peaks for different seed fields, see Table 1. We see that while the seed field value changes by 13 orders of magnitude, all fractions change by less then a factor 1.3. This is due to the rapid amplification processes which effectively erase knowledge of the primordial seed fields (see also Pakmor et al. 2014b) explored in this work. Note that these distributions include not only galaxies and clusters but also larger volumes potentially affected by outflows, and the surrounding intergalactic medium, as we discuss in more details below.
To study the large-value component of the magnetic field more quantitatively we analyze the main simulation box TNG100 with a seed magnetic field 10 −14 cG. Let us consider its distribution versus electron number density at different redshifts, shown in Fig. 3 (see also Marinacci et al. 2018). At low redshift z 2 we see two distinct branches, corresponding to weak and strong magnetic fields, respectively. In both branches, the value of the magnetic field is correlated with the electron number density. Even in the regions with small values of electron number density ∼ 10 −8 cm −3 (onetenth of the average electron number density today), the magnetic field can be many orders of magnitude stronger than its average value (close to the value of the initial field). This strong-B (or "over-magnetized") branch in Fig. 3 corresponds to the primordial seed independent strong-B peak in Fig. 2. In Appendices A, B, C we present plots similar to Fig. 3 that illustrate that two branches also exist for z < 2 for different initial conditions and different box sizes, while for z > 2 there is generically only one branch in the B-ne plot, due to the time needed for exponential amplification via a small-scale turbulent dynamo. This amplification process is faster at higher numerical resolution, enabling the magnetic fields to reach their quasi-saturated values earlier, although the level of this final saturation is relatively unaffected by resolution (Pakmor et al. 2017). In such a relatively small volume, there are also no high-mass halos at early times (z 3), such that large outflowdriven bubbles have not yet formed. Either larger volumes, i.e. probing the environments of the largest overdensities, or higher resolution would therefore if anything enhance the importance of these structures.
Our goal is to understand the impact of the outflow-generated large-B component of the magnetic field distribution on the propagation of light through the Universe. Specifically, what is the probability for a photon to occupy a large-B region on its way towards an observer? Fig. 4 therefore shows the volume fraction of regions with large magnetic field values, excluding the regions where the electron concentration is larger than some given value (in this way we can, in particular, exclude the inner parts of the collapsed structures like galaxies). 3 At z = 0 in regions with ne < 100 ne the magnetic field is stronger than 10 −12 cG in ∼ 14% of the volume, while it is stronger than 10 −10 cG in ∼ 8% of the volume. We note that the value of the initial field is 10 −14 cG. Moving only to z = 0.5, we see that strong-B regions occupy half as much volume, indicating that strong-B regions are substantially enhanced by late time processes.
Alternatively, we can measure the fractional length, for a given line of sight, which intersects a strong magnetic field. Using our sample of 1000 random sightlines we show in Fig. 5 the distribution of length fractions having magnetic field strength B 10 −12 cG at four different redshifts. 4 At z = 0 more than half of all sightlines intersect such strong magnetic fields along more than 10% of their path length. The peak of the distribution shifts to smaller fractional path lengths towards higher redshift. At z 1, sightlines intersect magnetic fields of this strength only rarely.
For further applications, it is also interesting to consider longer lines of sight that do not fit into one ∼ 100 cMpc snapshot. Using the TNG simulation volume we constructed continuous lines of sight from z = 0 to z = 6 following the procedure described in Bondarenko et al. (2020). To construct the magnetic field along a continuous line of sight we take B(z) from additional random sightlines within the same snapshot, and assign it to any missing pathlength (between simulation snapshots) of the continuous LOS, rescaling as B ∝ (1 + z) 2 . An example of both electron number density and magnetic field strength along a single continuous sightline is shown in Fig. 6.

Connection between over-magnetized 'bubbles' and galaxy outflows
Visual inspection of Fig. 1   same values of ne, see Fig. 3 and discussion), metallicity, and the structures seen in the electron number density. This is consistent with the behavior seen in the TNG galaxy formation model where strong outflows can escape from galaxies and break out into the intergalactic medium, carrying mass, heavy elements, and magnetic energy density along the way (Nelson et al. 2019b). This behavior is also consistent with results for galactic and near-galactic magnetic fields discussed in Butsky et al. To explore the physical connection between large-scale overmagnetized 'bubbles' and galactic-scale feedback processes, Figure 7 shows the magnetic field strength in thin (20 kpc) slices of TNG100 at z = 0. In the left panel, we mark all dark matter halos with total mass > 10 11.5 M with red circles, where the marker size denotes 1.5 R200, where R200 is the virial radius of a halo. In the right panel, we instead mark all supermassive black holes which have injected a significant amount of energy in the lowaccretion state, E low > 10 58.5 erg, with stars. We see that most magnetized bubbles extending to Mpc scales are directly associated with massive halos and/or supermassive black holes near their center. This is, however, not always the case. Figure 8 zooms into a crowded region of the same slice, within which an over-magnetized bubble is evident. The three panels on the right show magnetic field strength, gas metallicity, and electron number density. There is no clear association between cosmic web filaments visible in ne and regions of strong |B|. We have specifically selected this region as having no supermassive black holes which satisfy our energy criterion. Their absence implies that either the SMBH sourcing this bubble is outside the ±2.5 Mpc vicinity of the slice or that this large bubble may be a collective effect of galactic-scale winds produced by ongoing supernovae explosions in lower mass, Figure 4. Volume fractions of the regions where the magnetic field is larger than B min (legend) and the electron number density is smaller than ne,max (x-axis) for z = 0 (left panel) and z = 0.5 (right panel). We here show TNG100 with a seed field of B 0 = 10 −14 cG using data along 200 random lines of sight with (20 ckpc) 3 voxels in the box (see Sec. 2.4). Comparing these two redshifts, a substantial growth of the volumes occupied by high magnetic field strengths is evident, implying that physical processes within the last 5 Gyr have had a strong impact. Results for other redshifts and simulation boxes are given in Appendices A, B, and C. Figure 5. Probability to find a fraction of length along the line of sight with magnetic field larger than 10 −12 comoving Gauss. We show results for TNG100 at redshifts z = 0, 0.2, 0.5, and 1, creating 1000 random lines of sight at each redshift. While there is negligible 'strong magnetic field path length' at high redshifts, by z = 0 roughly half of all sightlines intersect such strong magnetic fields along more than 10 percent of their length. Results for other redshifts and simulation boxes are given in Appendices B and C.   Fig. 1, from the TNG100 simulation. Within ±2.5 Mpc of this slice we indicate all halos with mass 10 11.5 M and above (red circles; left panel, radii corresponding to 1.5 R 200 ) and all black holes (stars; right panel) which have injected significant amount of energy back to the gas (with E low > 10 58.5 erg, see discussion in the text). Both are strongly correlated with the presence of large-scale magnetized bubbles.
star-forming galaxies, possibly in combination with the effect of a number of smaller black holes.
To better understand the origin of these bubbles, Fig. 9 shows the magnetic field, gas temperature, dark matter, and electron number densities in a region extending 12 Mpc × 14.1 Mpc × 15 Mpc that contains the same large over-magnetized bubble from the Fig. 8. In this volume, we see a clear filament of large-scale structure. We again mark the massive halos and energetic black holes. Particularly clear around the SMBHs in the center and upper right are signatures of a collimated, episodic outflow made up of successive gas shells, forming a butterfly-like morphology, which results from the SMBH kinetic wind feedback. Here the ejection of magnetic fields into these bubbles is directly caused by active SMBHs.
However, stellar feedback can also contribute, most notably through core-collapse supernovae (SN) explosions. These also produce galactic-scale outflows, albeit in lower mass galaxies. To isolate these two mechanisms, we turn to additional simulations of the 'TNG model variation' suite. In particular, Fig. 10 compares the fiducial model (upper left) to models with no SMBH kinetic wind feedback (lower left), no SMBH feedback of any kind (upper right), and no feedback whatsoever from either black holes or supernovae (lower right). In all cases, the identical 25 cMpc/h volume is shown. The regions present in the fiducial run, but missing when the kinetic winds are disabled, are due to the low-accretion feedback mode of the SMBHs. The overall extent of the central filamentary region is such an example. Furthermore, regions present in the runs without black hole feedback, but missing in the 'no feedback' case, are due to supernovae -some of the smaller bubbles towards the up- Figure 8. Zoom on a magnetized bubble region from TNG100 selected to have no supermassive black holes satisfying the energy threshold, E low > 10 58.5 erg, within 5 Mpc of the slice. Further, this region is not clearly associated with an overdensity in ne, implying that some large-scale magnetized bubbles can arise from the combined action of many lower mass galaxies, hosting less effective black holes and/or supernovae-driven outflows.
per right being prime examples. In general, bubbles inflated by SN rather than SMBHs tend to be smaller, have lower temperatures, and lower expansion velocities. Overall, we see that both SMBHs and SN contribute to extended regions of high magnetic energy, with the supernovae playing a sub-dominant role.

AXION CONVERSION PROBABILITY
In this section, we apply our findings to calculate the probability of photon-axion conversion during propagation through the IGM. This conversion occurs in the presence of an external magnetic field B due to the interaction term where ga is the ALP constant (with units of inverse energy), Bi are components of the magnetic field, and Ai are spatial components of Aµ in the gauge A0 = 0. For ALPs with energy Ea in an external magnetic field this interaction effectively gives mass-mixing between the ALP and the photon, where BT = |B| cos θ and θ is an angle between vector A and the magnetic field. This leads to oscillations between axions and photons (Raffelt & Stodolsky 1988). The strength of the mixing between axions and photons is proportional to the magnetic field as well as the energy of the axion. The conversion probability also depends on the axion mass ma and the effective photon mass in the medium. For soft enough photons propagating through IGM the effective photon mass is given by (see e.g. Mirizzi et al. 2009) where me is the electron mass, αEM is the fine-structure constant and ne is a free electron number density. It should be noted that for gamma-ray photons this description of the effect on the effective properties of photons may not be sufficient (see e.g. Dobrynina et al. 2015). However, in this paper we deal with less energetic photons and therefore we can use Eq.
In the case when the effective photon mass is equal to the ALP mass, ma = mA(ne), the conversion becomes resonant and the conversion probability significantly increases. The conversion probability for this case is (Mirizzi et al. 2009) where , Ea is the axion energy, BT is the component of the magnetic field orthogonal to the line of sight (direction of axion propagation), is the distance along the line of sight, and the derivative for the R factor is calculated at the point of the resonance. The case p 1 is called the adiabatic limit and the conversion probability is close to one, while the case p ≈ 1 is called the non-adiabatic limit and the conversion probability is given by Let us assume that the total conversion probability Ptot during axion propagation along the LOS is much smaller than one. In this case the total conversion probability is given by where Ea is the axion energy at z = 0, and the sum is taken for resonances along the line of sight: zi is the redshift and BT,i is the

Probability of axion conversion along random sightlines
As the mixing strength between axions and photons is proportional to the magnetic field strength squared, the probability of the axionphoton conversion is dominated by the contribution of the strong-B component of the IGMF distribution. As we have seen above, this part of the distribution is to a large extent universal -its dependence on the value of the initial seed field is negligible (at least within the IllustrisTNG framework and within the class of initial conditions used in this paper). Further study of the dependence on initial conditions from a wider class (e.g. those of Vazza et al. 2017) and with different models of baryonic physics will be considered in future work.
We proceed with our analysis based on the TNG simulations with initial magnetic field strength 10 −14 cG. We take into account only the contribution to conversion from the simulation pixels with the value of magnetic field B 10 −12 cG (see the black line in Figs. 3 and 6). As we discussed in Section 3, for such magnetic field values, the distribution of magnetic field only slightly depends on the properties of the primordial seed field for seed fields 10 −14 cG and below. As a result, this threshold on magnetic field strength (in comoving units) returns a somewhat conservative contribution to the axion conversion probability from the magnetic bubbles.
Using 500 simulated continuous lines of sight (as described in Section 2.4), we calculate the average conversion probability and its distribution between redshifts z = 0 and z = 6. The result as a function of resonant electron number density is shown in Fig. 11. We see that conversion probability is maximal near ne,res = 10 −5 cm −3 . The scatter grows for large and small values of ne,res because of the small amount of resonances along the line of sight. In the region of resonant electron number densities below 10 −8 cm −3 or above 10 −2 cm −3 the resonant condition occurs rarely (see Fig. 6). We therefore derive below the constraints on the axion-photon coupling for the axion mass range 4 · 10 −15 eV < ma < 4 · 10 −12 eV.

Constraints from CMB distortions
In this section, we consider the effect where a CMB photon converts into an axion. The probability of this conversion is proportional to the photon energy, and its occurrence induces distortions in the CMB spectrum. The strength of the effect depends on the axion coupling ga, so we can constrain the axion model from CMB spectrum measurements obtained by COBE/FIRAS (Fixsen et al. 1996), which determined the CMB spectrum in the frequency range from 68 to 637 GHz with a precision of up to ∆BE/BE ≈ 10 −4 , where BE is a measured spectral radiance and ∆BE is its uncertainty.
Resonant axion-photon conversion modifies the CMB spectrum as where B CMB E is the spectral radiance of the initial CMB spectrum and Ptot(E) is the average conversion probability. We see that conversion produces an energy-dependent modification of the Planck spectrum. We estimate the exclusion region of the COBE/FIRAS measurement by a simple condition that follows from (7), where TCMB = 2.7260(13) K (Fixsen 2009). The result for the exclusion is shown in Fig. 12, where we also add an estimation of the sensitivity of the future CMB distortion experiment PIXIE using the same condition (8) but taking an expected sensitivity for the PIXIE experiment of ∆BE/BE < 10 −7 (Chluba et al. 2019a,b). As we can see, the CMB-based constraints are not competitive with other existing constraints. The reason for this is clear -the probability of axion-photon conversion is proportional to energy (see Eq. (6)). Therefore, much stronger constraints can be obtained from sources of photons with higher energy, e.g. X-ray or, especially, gamma-rays, see e.g. Montanino et al. (2017); Reynolds et al. (2019), where however non-resonant conversion is discussed. X-ray and gamma-ray sources are not allsky, but individual point sources, and a study of the effect of the IGM on their spectra requires a different methodology that is beyond the scope of this paper. We expect however that our results can be used for such an analysis in the future.

SUMMARY, DISCUSSION AND CONCLUSIONS
In this paper, we have quantified the effects of galaxy evolution processes on the magnetization of the intergalactic medium (IGM).
To this end, we have used several simulations from the Illus-trisTNG (Nelson et al. 2018;Marinacci et al. 2018) suite which all include treatment of ideal magnetohydrodynamics coupled to a state-of-the-art model for galaxy formation physics and feedback. In these calculations we always assume a homogeneous initial seed field, which is then amplified throughout the process of collapse and structure formation. As demonstrated in previous analyses (see e.g. Pakmor et al. 2017), strong magnetic fields are produced inside galaxies due to small-scale MHD dynamos. In this paper, we have shown that such strongly amplified magnetic fields can be distributed to larger volumes due to galactic feedback, in particular feedback from supermassive black holes (see also Nelson et al. 2019b). These largescale bubbles produced by outflows from galaxies and clusters develop particularly at redshifts z < 2, and contain magnetic fields that are several orders of magnitude stronger than in the unaffected regions of the IGM with the same electron density.
As a result, similarly to the magnetic fields inside galaxies, these fields are to a certain extent invariant with respect to the assumed initial conditions for magnetic fields (see Fig. 2 and its discussion in the text). We show that these over-magnetized bubbles with |B| 10 −12 cG, enhanced metallicity, and with clear outflowing kinematic signatures, can be as large as tens of Mpc. Their existence, and extent, is directly related to feedback activity from supermassive black holes in the centers of massive galaxies. Supernovae explosions also produce similar albeit smaller magnetized bubbles around lower mass galaxies.
We study the volume filling fraction of these strong field regions and their distribution over random lines of sight. We find that a typical intergalactic line of sight at z = 0 extends for ∼ 10−15% of its length within these regions for seed magnetic fields lower than 10 −14 cG. This implies that strongly magnetized bubbles are important for the propagation of light and high-energy cosmic rays through the Universe.
We use these results from TNG to put bounds on the photonaxion conversion from spectral distortions of the CMB. The disappearance of CMB photons due to their resonant conversion into axions in the IGM introduces deviations from the black-body spectrum of the CMB. If the mass of the axion is in the range 4 · 10 −15 eV < ma < 4 · 10 −12 eV where many resonances happen along a typical LOS, this distortion can be above the limit of COBE/FIRAS. As the photon-axion coupling grows with energy, the bounds based on the CMB distortions are not competitive with other existing bounds obtained using much more energetic photons.
In addition to CMB spectrum distortion, the conversion of CMB photons into axions produces additional anisotropy in both the temperature and polarization of the CMB. These bounds could be stronger than those obtained using the spectrum averaged over the whole sky. Indeed, the effect of the axion-photon conversion in the IGM is dominated by the contributions of the over-magnetized bubbles discussed in the paper. On one hand, this makes this effect to some extent independent of the unknown properties of the initial seed fields and therefore justifies the use of simulations. On the other hand, this means that anisotropies are introduced on the scales that are small for the CMB, i.e. the signal peaks at relatively high-, where the power spectrum of CMB anisotropies is already suppressed. At the same time, the size of the bubbles is noticeably larger than that of dark matter halos themselves, so the resulting impact could be distinguished from e.g. the Sunyaev-Zeldovich effect. However, to perform such an analysis we cannot calculate the conversion probability only for infinitesimal sightlines, and must instead produce an anisotropy map for a part of the sky, which is a challenging task for future work.
Stronger constraints from photon-axion conversion in IGM could also leverage the energy dependence of the coupling constant in order to study the propagation of gamma-ray photons through the IGM (see e.g. Montanino et al. 2017 for a study of the non-resonant conversion of gamma-rays in the IGM using ENZO simulations). This involves different classes of sources and will also be studied elsewhere.
Finally, the over-magnetized bubbles may contribute to measurements of the Faraday rotation measure (RM). For example, O'Sullivan et al. (2020) studies 349 pairs of radio sources in order to separate the contribution to RM from the IGM versus from the contribution of the local environments of the source and the observer (Milky Way), aiming to put a bound on the value of the primordial seed fields. However, the contribution of the overmagnetized bubbles encountered on the way from the source to the observer may be significant and could hinder constraints on the contribution of the "true" intergalactic magnetic fields that exist in regions unaffected by outflows.
To illustrate this, we show in Fig 13 the value of the Faraday Rotation measure for 1000 LOS as a function of the distance to the source (blue shaded bands). To disentangle the contribution of the bubbles, we choose the lines of sight such that they do not have magnetic field B > 10 −12 cG at the beginning and at the end of the LOS. We also exclude the contribution of voxels with electron density ne > 0.01 cm −3 , to remove galaxies, and of voxels with magnetic fields B < 10 −12 cG, to include only the contribution of the bubbles. We also show the region that contains 90% of values of RM of the radio sources observed by LOFAR presented  68% 90% LOFAR 90% Figure 13. Distribution of the contribution to Faraday Rotation Measure from the magnetic fields in outflow-generated bubbles only (see text), for sources located at redshifts between z = 0 and 1 calculated using 1000 random lines-of-sight from the TNG100 simulation (see 2 for description We see that for the sources located at z 0.2 the contribution of outflow-driven bubbles to RM measurements can be significant. This contamination must therefore be taken into account when inferring intergalactic magnetic field values from RM measurements. On the other hand, a more detailed analysis could constrain the properties (i.e. size scales, abundance) of feedback-driven bubbles themselves, thus informing models of galaxy formation and feedback, particularly when radio source counterparts with measurable distances are identifiable.
We conclude by noting that the baryonic feedback models and physics employed in the IllustrisTNG simulations are necessarily simplified treatments. In the future, more sophisticated black hole feedback models, such as those modeling unresolved accretion disks, black hole spin, and relativistic jet production could produce different emergent outflows. Similarly, the TNG galaxy formation model neglects several physical processes, most notably cosmic rays, low-temperature cooling and chemistry below 10 4 K, and coupled radiation-hydrodynamical interactions, which could similarly impact the generation and propagation of the outflows driven by galaxies and their supermassive black holes. Figure A1. Distribution of the magnetic field magnitude and electron number density in the fiducial 25 cMpc/h simulation at redshift z = 0 for different values of seed field: 10 −12 cG, 10 −14 cG, 10 −16 cG, and 10 −25 cG. The red dashed line represents the average electron number density.
In Fig. A1 we present for these simulations the analogous result as in Fig. 3 for z = 0. The initial conditions of the magnetic field seed appear to affect mainly the lower branch while the upper branch occupies magnetic field values around 10 −9 G which is in concordance with our result presented in Fig. 2. In Fig. A2 we show the volume fractions for each variation as seen in Fig. 4. We see that the volume filling fraction weakly depends on the seed magnetic field value, and the results for the 25 cMpc/h volume are in good agreement with the fiducial TNG100 simulation adopted throughout.

APPENDIX B: EVOLUTION OF MAGNETIC FIELD STRENGTH WITH REDSHIFT
We present in Fig. B1 the time progression for the correlation of electron number density and the magnitude of the magnetic field for the TNG100 simulation. We can observe how the upper branch develops with time as the universe evolves and becomes clearly distinguishable for z 1. Similarly, Fig. B2 shows the evolution of volume filling fraction for large values of magnetic fields. Fig. B3 displays the evolution of the length fraction of the LOS that has a magnitude of a magnetic field larger than 10 −12 cG. We see that the distribution starts to broaden circa z ≈ 2.

APPENDIX C: TNG300: DEPENDENCE OF THE MAGNETIC FIELD ON BOX SIZE
To explore the impact of the finite simulation size, we repeat our analysis on the larger TNG300 simulation, which however is performed at somewhat lower mass resolution than the TNG100 box: namely, dark matter and baryonic mass resolution of 5.9 × 10 7 M and 1.1×10 7 M respectively. In Fig. C1 we present the magnitude of the magnetic field versus the electron number density. Comparing with the results for the same redshift shown in Fig.3, we observe the same bimodal distribution within the same value ranges for both the TNG100 and TNG300-1 simulations. Comparing Fig. C2 (left) and Fig. 4 we see that TNG300 has smaller volume filling fraction (for example, for B > 10 −12 the volume filling fraction in TNG100 is 0.14, while in TNG300 is 0.135). We believe that this is an effect of poorer resolution of the TNG300 simulation in comparison to the fiducial TNG100.

Volume fraction
B>B min and n e <n e,max , z=0, B 0 =10 -25 cG Figure A2. Volume fractions of the regions where the magnetic field is larger than B min (legend) and the electron number density is smaller than ne,max (x-axis) for z = 0 for different values of seed field: 10 −12 cG, 10 −14 cG, 10 −16 cG, and 10 −25 cG. To produce these figures we used all gas cells in each of the simulation volumes. Figure B1. Distribution of the magnetic field magnitude and electron number density in the TNG100 simulation using data along 200 random lines of sight in the box at redshifts z = 0.1, 0.4, 1, 2, 3, and 6. The white dashed line corresponds to the comoving magnetic field value 10 −12 cG that we use as smallest values of the outflow-generated magnetic field in this work. The red dashed line represents the average electron number density at a given redshift. The gray dashed line show a power law B ∝ n 2/3 e , that should work for the adiabatic evolution. The seed field is B 0 = 10 −14 cG. Figure B2. Volume fractions of the regions where magnetic field is larger than B min and electron number density is smaller than ne,max as a function of redshift (different panels), from z = 0.1 to z = 6 using data along 200 random lines of sight in the TNG100 box. The seed field is B 0 = 10 −14 cG. Figure B3. Probability to find a fractional length along the line of sight with magnetic field larger than 10 −12 comoving Gauss for the 100 Mpc box at redshifts from z = 0.1 to z = 6. At each figure 1000 random lines of sight were taken from the simulation. Figure C1. Distribution of the magnetic field strength and electron number density in the TNG300 simulation using data along 420 random lines of sight in the simulation box at z = 0, with the seed field of B 0 = 10 −14 cG. The dashed white line corresponds to the comoving magnetic field value 10 −12 cG that we use as the smallest value of outflow-generated magnetic fields in this work. The red dashed line represents the average electron number density at a given redshift. The gray dashed line show a power law B ∝ n 2/3 e for the adiabatic evolution. . Left panel: volume fractions of the regions where magnetic field is larger than B min and electron number density is smaller than ne,max for z = 0 in the 300 cMpc box. The seed field is B 0 = 10 −14 cG, as in TNG100. Right panel: Probability to find a given fractional length along the line of sight with magnetic field larger than 10 −12 comoving Gauss for the 300 Mpc box at redshift z = 0, based on 420 random lines of sight.