© The Authors 2025

1. Introduction

Astrophysical magnetic fields can be found across a huge range of different scales, reaching from stars at small scales to the intracluster medium (ICM) in galaxy clusters (approximately megaparsecs). These magnetic fields can be inferred from Faraday rotation measures (RMs, Carilli & Taylor 2002; Osinga et al. 2022) and radio synchrotron observations (van Weeren et al. 2019). The discovery of radio megahalos even provides evidence for volume-filling magnetic fields out to the virial boundaries of galaxy clusters (Cuciti et al. 2022; Botteon et al. 2022). These observations have revealed that galaxy clusters are filled with magnetic fields with strengths of a few microgauss and coherence lengths at the order of a few kiloparsecs (Vogt & Enßlin 2005; Kuchar & Enßlin 2011). However, the generation and amplification mechanisms of magnetic fields within galaxy clusters remain uncertain.

Understanding cluster magnetic fields is crucial for understanding the diverse and rich radio phenomena in galaxy clusters that are linked to the transport and cooling of cosmic rays along magnetic field lines (Ruszkowski & Pfrommer 2023). Yet, the exact origins of these phenomena such as radio halos, relics, and diffuse radio emission near cluster centers, which reveal interactions between active galactic nuclei (AGNs) feedback and the magnetized ICM, remain uncertain (Brunetti & Jones 2014; van Weeren et al. 2019; Wittor et al. 2023). Recent radio observations conducted by the LOFAR Two-metre Sky Survey have unveiled that galaxy clusters at redshifts around z ∼ 1 have magnetic field strengths similar to low redshift clusters (Di Gennaro et al. 2021, 2023). This implies an early phase of magnetic field creation and amplification within galaxy clusters. Observations of high-redshift protoclusters in the redshift range 2.2 < z < 4.25 have unveiled interactions between the brightest cluster galaxy (BCG) and the magnetized ICM (Anderson et al. 2022; Cordun et al. 2023; Chapman et al. 2024), which gives hints that they might play a role in the high redshift magnetization of the protocluster.

There is a consensus that magnetic field growth in galaxy clusters starts with a small seed field that becomes amplified. The seed field origin can be attributed to magnetic field generation in the early Universe, such as through the Biermann battery mechanism that turns baroclinic motions in ionized gas into magnetic fields (Biermann 1950; Kulsrud et al. 1997; Kunz et al. 2022), the Weibel instability (Weibel 1959; Huntington et al. 2015; Zhou et al. 2024), or cosmological phase transitions in the early Universe (Widrow 2002; Durrer & Neronov 2013). The most widely discussed amplification scenario involves the turbulent small-scale dynamo, which converts turbulent kinetic energy into magnetic energy (Kazantsev 1968; Kraichman 1966; Rincon 2019; Kunz et al. 2022). Dynamo theory makes use of the fact that a turbulent magnetic field flux-frozen into the thermal plasma can evolve below the viscous scales, where the turbulence in the velocity field is dissipated, and that it grows rapidly at the resistive scales where timescales are very short (see, e.g., Kunz et al. 2022). Turbulence within the ICM arises mostly from mergers and interactions with galaxies (Kulsrud et al. 1997; Norman & Bryan 1999; Subramanian et al. 2006; Creasey et al. 2013; Wittor et al. 2017; Vazza et al. 2017; Li et al. 2020; Bacchini et al. 2023). It can also result from plasma instabilities such as the magneto-thermal instability (Balbus 2000; Perrone et al. 2024a,b), pressure anisotropies (Rappaz & Schober 2024), and magnetic pumping (Ley et al. 2024), which may all be able to self-regulate magnetic field growth.

According to the theory of hierarchical structure formation, galaxy clusters form through frequent mergers with galaxies and smaller clusters over cosmic time (Kravtsov & Borgani 2012). In this paper, we demonstrate how the magnetic evolution of galaxy clusters is linked to that of the galaxies themselves. During the collapse of the protocluster, the primordial magnetic field is adiabatically compressed and turbulence is generated, which initiates dynamo processes as seen in cosmological simulations of galaxy clusters (e.g., Vazza et al. 2014; Domínguez-Fernández et al. 2019; Steinwandel et al. 2024). In galaxies, gravitational infall and stellar feedback processes drive turbulence, which further contributes to dynamo processes (e.g., Schober et al. 2013; Käpylä et al. 2018a; Pfrommer et al. 2022). After a seed magnetic field has been amplified through a small-scale dynamo in galaxies, it can be transported and expanded into intergalactic space by powerful galactic winds or jets or ram pressure stripping (Rees & Setti 1968; Rees 1987; Rephaeli 1988; Furlanetto & Loeb 2001; Fujita 2004; Hou et al. 2014; Dolag et al. 2016; Martín-Navarro et al. 2019; Serra et al. 2023). These mechanisms impact galaxies within the cluster but also those in pre-merging substructures, a phenomenon known as “pre-processing” (Fujita 2004; Serra et al. 2023) in the context of metal enrichment. Observations and simulations support this, showing that many galaxies in cluster environments lose their self-bound gas reservoirs (Butcher & Oemler 1984; Hou et al. 2014; Dolag et al. 2016; Martín-Navarro et al. 2019; Hasan et al. 2023). Cluster mergers additionally inject turbulence into the ICM, which enables further dynamo activity (e.g., Vazza et al. 2014; Domínguez-Fernández et al. 2019; Steinwandel et al. 2024).

Numerical simulations serve as valuable tools for investigating the evolution of magnetic fields in galaxy clusters (see, e.g., the review by Donnert et al. 2018). Previous studies have explored the creation and amplification of magnetic fields in isolated galaxy clusters (without radiative physics, Roh et al. 2019), galaxy clusters in cosmological simulations (without radiative physics, Dolag et al. 1999; Miniati & Beresnyak 2015; Vazza et al. 2018; Steinwandel et al. 2022, 2024), and cosmological simulations with passive magnetic fields (Adduci Faria et al. 2020) or active magnetic fields while including radiative physics and subgrid models for galaxy formation (Marinacci et al. 2015, 2018; Nelson et al. 2024). Most simulations successfully reproduce the magnetic field evolution to match observed values of a few microgauss within the age of the Universe.

Such simulations necessitate the coverage of a wide range of scales, spanning from small events such as supernovae (SNe) to large events such as cluster mergers. In ideal magneto-hydrodynamical (MHD) simulations, the growth rates depend on the cell sizes, as these determine the shortest eddy turnover scales (see e.g., Pakmor et al. 2017; Whittingham et al. 2021; Pfrommer et al. 2022). Therefore, simulations require high resolution to reproduce the strong magnetic fields observed at high redshifts. While Steinwandel et al. (2022) found that the magnetic field only recently became saturated, the higher resolved simulations of Steinwandel et al. (2024) observed a saturated magnetic field already at z ≈ 4.

However, increasing the resolution beyond the viscous scales theoretically necessitates the inclusion of kinetic plasma effects (e.g., Berlok et al. 2020; St-Onge et al. 2020; Steinwandel et al. 2024). A related question is whether dynamo theory is applicable in the ICM, given that Kazantsev theory is a fluid theory. Due to the weakly collisional nature of the ICM (see, e.g., Schekochihin & Cowley 2006; Kunz et al. 2011, 2012 and references therein), the fluid approximation is formally invalid, and thus any MHD description of the ICM would also be invalid. Additionally, microscale kinetic instabilities, specifically the firehose and mirror instabilities on small scales (see, e.g., Schekochihin & Cowley 2006 and Majeski et al. 2024), can increase the effective collision rate and reduce the effective viscosity. Cluster observations indicate that the turbulent dissipation scale is smaller than expected from Spitzer viscosity (Zhuravleva et al. 2019; Heinrich et al. 2024), which is theorized to be due to these small-scale instabilities. However, in high-beta plasmas such as the ICM (Carilli & Taylor 2002), it has been argued that the pressure anisotropy will self-regulate by reorganizing turbulent motions (Squire et al. 2019; Majeski et al. 2024) and that this so-called magneto-immutability (Squire et al. 2019) can reduce the viscosity without triggering the small scale instabilities (Majeski et al. 2024). Simulations that include kinetic plasma instabilities via MHD with anisotropic viscous stress (Braginskii 1965) show only minor deviations from the saturated state of the normal MHD dynamo (St-Onge & Kunz 2018; Squire et al. 2019; St-Onge et al. 2020), and hybrid-kinetic particle-in-cell simulations have shown a similar magnetic field growth to the MHD dynamo model (Achikanath Chirakkara et al. 2024).

However, these studies only considered magnetic dynamo action in weakly collisional conditions that are found in the ICM at the current epoch and did not self-consistently account for the magnetic dynamo in the plasma that evolves into the ICM throughout cosmic evolution. MHD simulations that include radiative physics can naturally resolve higher growth rates due to higher densities and smaller cell sizes in the interstellar medium (ISM). In the ISM, the fluid description is valid again due to the higher particle collision rates. It has been shown that simulations of magnetic fields in galaxies reveal rapid early growth (Gressel et al. 2008; Schleicher et al. 2010; Pakmor et al. 2014, 2017). This fast growth is linked to SNe from early generations of stars and gravitational collapse due to structure formation processes (Schober et al. 2013; Käpylä et al. 2018a; Pfrommer et al. 2022). Cosmological galaxy simulations have also shown that magnetized galactic outflows can occur and magnetize their surroundings (Pakmor et al. 2020; van de Voort et al. 2021). This provides a method to magnetize the ICM in radiative simulations, as magnetic fields can grow rapidly in the ISM, and the magnetized gas is then transported outward via feedback mechanisms.

Marinacci et al. (2015) demonstrated that incorporating radiative physics results in faster magnetic field growth in the ICM compared to non-radiative simulations, although they did not investigate the exact sources of this accelerated growth. This leaves the contributions of galaxies to the magnetic field evolution in galaxy clusters unclear. In this paper, we expand on the work of Marinacci et al. (2015) by comparing radiative and non-radiative simulations of galaxy clusters. We employ higher resolution and updated galaxy formation physics. Furthermore, we extend the dynamo analysis to shed light on the mechanisms of magnetic field growth and transport.

In this study, we conduct two sets of each of the four cosmological galaxy cluster simulations. In one set, the radiative physics are enabled, and in the other set, the radiative physics are disabled. We show that the magnetic field grows faster in the radiative simulations due to the enhanced growth rate inside galaxies. Furthermore, we show how the magnetic field in the clusters grows with subsequent mergers of galaxies. The paper is organized as follows. In Sect. 2, we describe our simulation setup. We give properties of the simulated clusters, and we introduce the cosmological code and model as well as the galaxy formation model. In Sect. 3, we show how the magnetic field evolves with time in the two simulation sets. In Sect. 4, we investigate the origin of the magnetic field growth, including a detailed analysis on the dynamo activity in our simulations. In Sect. 5, we present mock high redshift Faraday RM observations of the protocluster. Finally, in Sect. 6, we discuss our findings, summarize our most important results, and give an outlook on testing early magnetic field amplification in protoclusters, motivated by upcoming observational surveys.

2. Methods

We analyze the magnetic field evolution in four different galaxy clusters which have been simulated with two distinct physical models using the zoom-in technique. One type of simulation uses MHD and gravity (‘non-radiative’ simulations), while the other one additionally includes subgrid models for galaxy formation and feedback physics (‘radiative’ simulations). These simulations are part of the PICO-Clusters simulation suite, which will be described in detail in a forthcoming paper (Berlok et al, in prep.). The PICO-Clusters simulations employ the moving-mesh code AREPO (Springel 2010; Pakmor et al. 2016b; Weinberger et al. 2020) for evolving the MHD equations (Pakmor & Springel 2013) in comoving form and uses an updated version called AREPO-2 that includes the SUBFIND-HBT algorithm of GADGET4 (Springel et al. 2021, 2022). To simulate radiative physics, the galaxy formation model used in the IllustrisTNG simulations (Weinberger et al. 2017b; Pillepich et al. 2018) is employed. Brief descriptions of AREPO and its ability to correctly evolve magnetic fields (Sect. 2.1), IllustrisTNG (Sect. 2.2) and PICO-Clusters (Sect. 2.3) are provided below. We summarize the details most important to this paper in order to keep it self-contained. An overview of the physical characteristics of the four selected clusters at redshift zero is shown in Table 1, including the virial radius, virial mass, and the number of subhalos at z = 0. Here the virial radius (R₂₀₀) is the radius that contains a mean density of 200 times the critical density, and the virial mass (M₂₀₀) is the mass contained in a sphere of this virial radius. We also show the number of high resolution dark matter particles in each zoom12 simulation (see Sect. 2.3), each of which have a mass of 5.9 × 10⁷ M_⊙. This makes the mass resolution in our simulations comparable¹ to the TNG300 cosmological volume simulation (Springel et al. 2018) and the TNG-Cluster zoom-in simulations (Nelson et al. 2024). The halo numbers are assigned by sorting the M₂₀₀ masses of all clusters with M₂₀₀ > 10¹⁵ M_⊙ at z = 0 in the parent simulation, starting with the most massive cluster, Halo0. As a result, “Halo3” represents the fourth most massive halo in our simulation.

We chose to use a sample of four clusters to test for generality of the suggested mechanism of the cluster magnetic dynamo and limit the impact of selection effects. We then perform a more in-depth analysis of one of the clusters (Halo3) in order to elucidate the detailed physics of the mechanism of magnetic field growth. Halo3 and 4 are two of the most massive clusters in the PICO-Clusters sample and were selected for re-simulation because of their massive, low redshift, galaxy cluster mergers. Halo194 and 260 were randomly selected among the PICO-Clusters that have M₂₀₀ ≈ 10¹⁵ M_⊙ and were included in the analysis in this paper in order to ensure that our results are applicable to somewhat smaller clusters and not just to very massive clusters undergoing major mergers.

In our simulation analysis, we will distinguish between ICM and ISM by designating star-forming cells as ISM and non-star-forming cells within the cluster virial radius as ICM. When we refer to the “center” of the galaxy cluster, this is to be specifically understood as a spherical region of cluster-centric physical radius 20 kpc, which corresponds to the virial radius of the central galaxy in Halo3 at z = 9.5. We make this separation to distinguish between actions occurring in the BCG at high redshift and activity in the entire cluster², both of which contribute to the magnetic field enrichment. While the simulations output quantities in comoving units (Sect. 2.1), we show proper physical quantities throughout this paper if not otherwise mentioned. For velocity calculations in an expanding universe, the total (physical) velocity is given by the sum of the peculiar velocity and the Hubble flow:

$$\begin{array}{c}\hfill \mathit{v}=\dot{\mathit{r}}=a\frac{d\mathit{x}}{\mathit{dt}}+\mathit{x}\frac{\mathit{da}}{\mathit{dt}}={\mathit{v}}_{\mathrm{pec}}+H(a)\mathit{r},\end{array}$$(1)

where x is the comoving position, r = ax is the physical position, and H(a) = ȧ/a the Hubble function. While the Hubble flow contributes to the total velocity, it does not contribute to the vorticity,

$$\begin{array}{cc}\hfill \mathit{\omega }& =\mathbf{\nabla }\mathbf{\times }\mathit{v}=\mathbf{\nabla }\mathbf{\times }{\mathit{v}}_{\mathrm{pec}},\hfill \end{array}$$(2)

because it is irrotational, i.e. ∇ × r = 0.

2.1. AREPO and magnetic fields

AREPO discretizes MHD quantities on a moving unstructured mesh that is defined by the Voronoi tessellation of a set of discrete points that can be moved arbitrarily. If those mesh-generating points follow the motion of the gas, the code inherits the advantages of Lagrangian methods. The Voronoi mesh can adaptively refine and de-refine. The standard refinement criteria in AREPO ensures that cell masses remain within a factor of 2 of some target mass, but additional criteria can be used to enhance the resolution in areas of interest relevant to the particular simulation in question. Dark matter, stars, and black holes are included as collisionless particles. Gravitational interactions are computed using a TreePM method (Springel 2005) to separate long- and short-range gravitational forces. Here the long-range forces are determined using a Fourier transform on a mesh, while the short-range forces are solved using an oct-tree algorithm (Barnes & Hut 1986).

The comoving ideal MHD equations are solved in AREPO on the moving Voronoi mesh (Pakmor et al. 2011; Pakmor & Springel 2013). The equations are solved using an HLLD Riemann solver (Miyoshi & Kusano 2005) and the code is able to reproduce comoving analytic wave solutions for both Alfvén and magnetosonic waves (Berlok 2022). To mitigate divergence errors in the magnetic field, the Powell 8-wave divergence control scheme (Powell et al. 1999) is utilized, which subtracts the divergence error from the initial state, enhancing overall accuracy and ensuring stability. Simulations that use this MHD scheme can reproduce analytic growth rates and magnetic curvature statistics expected for a fluctuating dynamo (Pfrommer et al. 2022) as well as different observables, for example, magnetic field saturation strengths and radial profiles in disk galaxies (Pakmor et al. 2014, 2017), Faraday RMs in spiral galaxies and the Milky Way (Pakmor et al. 2018; Reissl et al. 2023), and RM values of fast radio burst in normal star forming galaxies (Mannings et al. 2023). Furthermore, the MHD setup is able to reproduce observed magnetic fields in dwarf galaxies, galaxies and groups (Pakmor et al. 2024), and in the circumgalactic medium (Pakmor et al. 2020) via dynamo activity. The divergence error due to the Powell cleaning method remains low in dynamic galaxy mergers (Whittingham et al. 2021, 2023), which shows that magnetic field amplification is not an artifact of the code.

We initialized a homogeneous magnetic seed field of comoving field strength 10⁻¹⁴ G at z = 127. Our chosen initial field strength is strong enough to saturate at z = 2 in galaxies at our numerical resolution (Pakmor et al. 2014), yet still weak enough to not suppress star or galaxy formation (Marinacci & Vogelsberger 2016). However, within a broad range, variations in seed strength show little to no impact on the later evolution as the initial memory is quickly erased by dynamo activity and the similar saturation strength (Pakmor & Springel 2013; Pakmor et al. 2014; Marinacci et al. 2015). This property in combination with our finite numerical resolution justifies our choice of seed field strength, which is larger than expected from cosmological Biermann battery processes (Kulsrud & Zweibel 2008).

The efficiency of the fluctuating dynamo, as well as the morphology of the resulting magnetic field, depend on the smallest scales, particularly the viscous and resistive scales. While the ideal MHD equations assume a perfectly conducting fluid (neglecting resistivity and viscosity), numerical effects stemming from finite volume and time discretization introduce numerical viscosity and resistivity. The numerical resolution defines the numerical viscous and resistive scales, equivalent to the smallest eddy scales. The utilized (quasi-)Lagrangian moving mesh with an approximately constant gas mass per resolution element yields higher resolution in dense regions (with cell sizes of ∼100 pc in the center) and degrades the resolution in low-density regions (reaching resolutions of ∼10 − 100 kpc in the outskirts). This enables us to resolve resolution-dependent Reynolds numbers of ∼100 − 1000 in the ICM.

2.2. The IllustrisTNG galaxy formation model

We simulated radiative physics using the IllustrisTNG galaxy formation model (Weinberger et al. 2017b; Pillepich et al. 2018), which encompasses star formation and feedback (Vogelsberger et al. 2013; Pillepich et al. 2018), supermassive black hole seeding, growth, and resulting AGN feedback (Springel et al. 2005a,b; Vogelsberger et al. 2013; Weinberger et al. 2017b), as well as gas cooling and heating. Simulations using the IllustrisTNG model have been shown to reproduce observed stellar-halo mass functions of groups and galaxy clusters (Pillepich et al. 2018), and clustering properties, such as two-point correlation functions between low-redshift galaxies that match measurements from SDSS data (Springel et al. 2018). Magnetic field strengths of galaxy clusters, simulated using this model, are consistent with observations (Marinacci et al. 2018). In this paper, we expand on the work of Marinacci et al. (2015), who used the Illustris galaxy formation model (Vogelsberger et al. 2014) and additionally account for the magnetic field evolution using ideal MHD.

To model star formation and the ISM, we employed an effective model (Springel & Hernquist 2003; Vogelsberger et al. 2013; Pillepich et al. 2018). The formation of star particles occurs stochastically when the density exceeds a threshold of n_H > 0.106 cm⁻³ adhering to the Kennicutt-Schmidt relation that is based on observations (Schmidt 1959; Kennicutt 1989), and using a Chabrier initial mass function (Chabrier 2003). Mass loss and metal enrichment are calculated in each time step, considering the initial mass function, as outlined in Pillepich et al. (2018). Stellar feedback employs a kinetic wind model, wherein wind particles are decoupled from hydrodynamics, interacting only gravitationally. They recouple hydrodynamically upon encountering gas cells that fulfill a density threshold (n_H < 0.0053 cm⁻³, which corresponds to 5 percent of the star formation threshold) or, alternatively, after a defined travel time 0.025 × 1/H(z), where 1/H(z) is the Hubble time at the simulation redshift. At very high redshift, recoupling is primarily governed by the time criterion due to shorter Hubble times and higher critical densities. Upon recoupling, wind particles inject mass, metals, thermal energy, and momentum to gas cells. Wind particles form stochastically from star forming cells and are launched isotropically, with velocity correlated to local dark matter velocity dispersion and adjusted for redshift.

We include gas cooling and heating following the method outlined in Vogelsberger et al. (2013). A treatment for primordial cooling and heating through ionization equations based on cooling, recombination, and ionization rates for primordial gas was employed (Cen 1992; Katz et al. 1996) alongside Compton cooling off the cosmic microwave background (Ikeuchi & Ostriker 1986). Metal line cooling uses CLOUDY cooling tables (Wiersma et al. 2009). A spatially uniform, redshift-dependent ionizing UV background is included for gas heating (Faucher-Giguère et al. 2009) with corrections for self-shielding in dense gas (Rahmati et al. 2013).

The AGN treatment in IllustrisTNG is based on the AGN model developed by Springel et al. (2005b) and further developed by Vogelsberger et al. (2013) and Weinberger et al. (2017b). Each halo that reaches a friends-of-friends halo mass of 5 × 10¹⁰ M_⊙ is assigned a black hole at the position of the densest cell in the halo (Weinberger et al. 2017b). The black holes in the simulation accrete gas from their surrounding environment. A fraction of the accreted gas contributes the black hole mass. A fraction of the rest mass energy of the accreted gas is expelled as energetic feedback, some of which couples to the surrounding gas. The feedback model employed in this study is a two-state model, where the feedback mode depends on the accretion rate. This differentiation is based on the mass accretion rate in units of the Eddington rate (Ṁ_BH/Ṁ_Edd) that fulfills a threshold value that scales with black hole mass, but should not exceed 0.1. Below this threshold, the black hole operates in a kinetic AGN feedback mode, generating black hole-driven winds (Weinberger et al. 2017b). Otherwise, the black hole operates in quasar mode, releasing thermal energy isotropically to neighboring gas cells (Springel et al. 2005a,b). The quasar mode is predominantly active at higher redshifts.

2.3. PICO-Clusters and halo selection

Plasmas In COsmological Clusters (PICO-Clusters; Berlok et al., in prep.) is a new simulation suite of cosmological zoom-in simulations of galaxy clusters using the IllustrisTNG galaxy formation model with the moving mesh code AREPO-2 (an updated version of AREPO that evolved from the MilleniumTNG project (Pakmor et al. 2023); for details, see below). It consists of a cosmological volume simulation of a periodic box with side length ≈1.5 Gpc that contains 272 galaxy clusters with masses above 10¹⁵ M_⊙ and re-simulations of 25 clusters using the zoom-in technique. The long-term goal of PICO-Clusters is to understand the role of including various (plasma-)physical models in simulations, such as anisotropic viscosity (Berlok et al. 2020), anisotropic conduction (Talbot et al. 2025), cosmic ray protons (Pakmor et al. 2016a; Pfrommer et al. 2017) and electrons (Winner et al. 2019), and AGN feedback via a low-density momentum-driven jet (Weinberger et al. 2017a; Ehlert et al. 2023). The base-line physical model of PICO-Clusters is the IllustrisTNG galaxy formation model (see Sect. 2.2). In the present paper, we start this upcoming systematic investigation of the role of the physical model employed by comparing non-radiative simulations with the base-line IllustrisTNG model with a focus on magnetic fields.

All PICO-Cluster simulations are initialized at z = 127 and have simulation properties stored in snapshots. The four clusters analyzed in this paper have a high output frequency, i.e., around 260 full snapshots per simulation. This allows for a detailed analysis of the time evolution. In addition, the AREPO-2 code has better memory management and load balancing, and includes the SUBFIND-HBT algorithm (Springel et al. 2021, 2022), which allows group and subfind properties to be stored even more frequently than the full snapshots and merger trees to be easily calculated. We used the on-the-fly merger trees for tracking the centers of the galaxy clusters (Springel et al. 2022).

PICO-Clusters uses a Planck-2018 cosmology (Planck Collaboration 2020) with the following cosmological parameters Ω_Λ = 0.684, Ω_m0 = 0.316, Ω_b0 = 0.049 and the Hubble parameter h = 0.673. We used these parameters to create initial conditions for the parent simulation. We generated a power spectrum using the CAMB code (Lewis et al. 2000; Lewis & Bridle 2002), and N-GenIC (Springel 2015) was then used to create initial conditions based on Gaussian random fields, adhering to the power spectrum prescription (as detailed in Springel et al. 2005c; Angulo et al. 2012). The parent simulation had 1024³ dark matter particles, each with a mass of 1.0 × 10¹¹ M_⊙. This cosmological box was evolved to z = 0. The most massive halos at z = 0 were identified using a friends-of-friends algorithm (FoF).

The initial conditions for the re-simulations were created using a new code that will be described in detail in Puchwein et al., in prep. In brief, we select all dark matter particles in the parent simulation object out to 5R₂₀₀ and trace them back to the starting redshift, z = 127, of the parent simulation. For the zoom-in initial conditions, we then refine this region to the desired resolution using a tree structure. The rest of the volume is initially set to a coarse base resolution sampling the whole parent simulation volume. We then introduce boundary layers, in which we refine the particle distribution to intermediate resolution levels based on the distance to the nearest edge of the high-resolution region. This results in a particle load with a high-resolution region and boundary layers that reflect the irregular shape of the Lagrangian volume of the selected object, allowing us to concentrate the numerical effort in the re-simulations in the regions where it is needed, while avoiding contamination of the objects with low-resolution particles out to several virial radii. Perturbations are added based on the Zel’dovich approximation using the same phases and mode amplitudes as in the parent simulation, but adding additional small scale power as appropriate. For each resolution level of the particle distribution we add power up to the corresponding Nyquist frequency.

As a consequence, the clusters analyzed in this paper have no low resolution particles inside 3R₂₀₀ at z = 0. We emphasize that contamination by heavy dark matter particles in the high-resolution region would not only have the unfortunate property of artificially heating the high-resolution dark matter particles and the ICM via numerical dynamical friction but also would signal missing structure in the Lagrangian region in the initial conditions and an incomplete high-resolution formation history of the simulated cluster³. Each dark matter particle has a mass of 5.9 × 10⁷ M_⊙ inside a high-resolution region (for reference, this corresponds to the TNG300 resolution). The results presented in this paper are based on the zoom12 simulations, where the zoom-in regions have a mass resolution that is 12³ times higher than that of the parent simulation. The target gas mass is 7 × 10⁶ M_⊙. The high resolution dark matter, star and black hole particles have comoving softening lengths of 4.8 kpc above redshift 1 and a physical softening of 2.4 kpc at lower redshift. The gas softening depends on the cell size with a minimum comoving value of 0.6 kpc.

The use of the IllustrisTNG galaxy formation model makes the 25 base-line PICO-Clusters simulations roughly comparable to the TNG300 and TNG-Cluster simulations (Pillepich et al. 2018; Nelson et al. 2024). The differences between those simulations and our PICO-Cluster project are: 1) our initial conditions code, which ensures a contamination-free high-resolution region out to 2.5 R₂₀₀, and out to 3R₂₀₀ for the four clusters studied in this work, 2) the simulation code, where we use AREPO-2 while TNG300 and TNG-Cluster use AREPO, and 3) the scientific focus as we will be investigating plasma effects in the ICM in PICO-Clusters. We also simulate clusters with a better mass resolution so that we can address numerical convergence of ICM turbulence and magnetic field strength (see Berlok et al. in prep.).

3. Magnetic field evolution with time

In this section, we will compare magnetic field growth in the radiative and non-radiative simulations and explore environmental effects for the differences such as interactions between the ICM and galaxies.

3.1. Projections of magnetic field and vorticity

We present in Fig. 1 thin projections of Halo3 at z = 3 (upper panels) and at z = 0 (lower panels) in both simulation setups. The virial radius (R₂₀₀ = 380 kpc at z = 3 and R₂₀₀ = 3 Mpc at z = 0) is indicated with a white circle. From left to right, we subsequently zoom into the same projection. In each panel we show the magnetic field strength (left) and vorticity (right). At z = 3, the magnetic field is significantly larger in the radiative simulation whereas it looks very similar between the two simulations setups at z = 0.

Fig. 1. Thin projections through the center of Halo3 displaying the root mean square magnetic field strength and vorticity, respectively, at z = 3 in the non-radiative (first row) and radiative (second row) simulations. The bottom two rows show the same but at z = 0. The virial radius R200 is indicated with a white circle. From left to right, the panels have sidelengths L = 4R200,  2R200,  0.5R200 and a thickness of L/10. At z = 3, the magnetic field in the radiative simulation is nearly saturated. The BCG appears as a bright spot with a strong magnetic field. An object infalling from the left has also developed a strong magnetic field that will contribute to the cluster’s final magnetic field. By contrast, at z = 3 the magnetic field in the non-radiative simulation remains weak. The infalling object is not visible, and there is no BCG. This underscores the magnetic field growth in galaxies and its transport to the ICM via merging halos.

At z = 3 (upper panel of Fig. 1), a strong magnetic field with 30 μG has developed in the center of the radiative simulation. The magnetic field fills the volume in the cluster, with field strengths decreasing towards the outskirts to ∼10⁻² μG in the radiative simulation. The location of the BCG in the center shows a bright spot of magnetic field and reflects the very efficient amplification in the central galaxy. A merging cluster is visible on the left side. The merging cluster has grown its own strong magnetic field that will contribute to the main halo’s magnetic field after merging. Both findings highlight the physical mechanisms at work to amplify the magnetic field: the growth in galaxies and enrichment of the ICM via subsequent mergers. By comparison, in the non-radiative simulation, the magnetic field remains weak and only grows to ∼10⁻² μG in the center, where the densities are highest. The merging cluster has not grown a magnetic field that is above the plotting threshold of 10⁻² μG.

Turbulence and magnetic field strength in the ICM are intricately linked through potential dynamo interactions. Turbulent magnetic amplification depends on resolution, which is set to maintain a target gas mass per cell within a factor of 2. Since cluster density profiles follow a single/double beta profile (see Appendix A.1), the highest resolution is achieved in the central region and within galaxies (i.e., in the ISM). Resolution decreases outward following a power-law scaling. We show projections of the root mean square vorticity, which itself is defined as $\omega =\mathrm{\nabla }\times \mathbf{v}$, where the eddy turnover time is defined as $t_{\mathrm{eddy}}=1/|\omega |$. We compute the curl of the velocity from the velocity derivative matrix ${\partial }_i{v}_j$ and use the slope-limited velocity derivatives to construct this matrix. The projections look comparable in both simulations at z = 3. Remarkably, central filaments with the highest vorticity also coincide with regions displaying the strongest magnetic fields in the radiative simulation. Eddy turnover rates in the central regions are approximately 1 Myrs⁻¹ in both simulations, indicating a similar vorticity profile within and beyond the virial radius in the two simulations. This turbulence is primarily attributed to frequent mergers and accretion shocks in both simulations, rather than originating from radiative physics, which would naturally be only visible in the radiative simulation. Various studies come to the same conclusions (Norman & Bryan 1999; Vazza et al. 2011, 2017). Galaxies in the radiative simulation are visible as a few bright spots, reflecting the enhanced level of turbulence injected by stellar and AGN feedback as well as the small cell sizes.

In the lower panels, we compare the magnetic field strength and vorticity at z = 0. Both look very similar in the two simulations. In both runs, the highest magnetic field values – approximately 30 μG – are observed in the central regions coinciding with the highest density regions. Magnetic field strength exponentially decreases towards the outskirts in both simulations (we show radial exponential profiles in Appendix A.1). Beyond the virial radius, the magnetic field is slightly more extended in the radiative simulation. The highest vorticity values occur at the center, where density is greatest, resulting in smaller cell sizes and elevated vorticity. Ram pressure stripped galaxies are visible as bright spots with tails that are moving towards the cluster center. These phenomena are also observed in low redshift galaxies within galaxy clusters and groups (e.g., Poggianti et al. 2017) as well as in cosmological simulations (e.g., Rohr et al. 2023) and dedicated wind-tunnel simulations (e.g., Tonnesen & Bryan 2012; Sparre et al. 2024a,b) using similar galaxy formation physics.

The substantial vorticity levels in the ICM, and the strong spatial correlation between vorticity and magnetic field, is a first hint at the presence of a fluctuating dynamo mechanism in the ICM. Furthermore, this dynamo effectively amplifies magnetic fields to comparable strengths within the virial radius in both simulations at z = 0. Subsequent sections will quantitatively analyze magnetic field evolution and identify the sources of magnetization.

3.2. Radial magnetic field evolution in all four halos

In order to 1) demonstrate that the magnetic field grows faster in the radiative simulations of all four halos, 2) connect the growth in the ISM and the ICM, and 3) distinguish between different phases of magnetic field growth, we analyze the time evolution of the magnetic field strength. Figure 2 shows the time evolution of the volume-averaged root mean square magnetic field strength in the radiative simulation (upper panels), and in the non-radiative simulation (lower panels). In the left panels, we show as an example the evolution of the magnetic field in the ICM for Halo3 in more detail, including different radial ranges within the virial radius, indicated with the blue color scheme. The magnetic field in the ISM (i.e., by our definition, star-forming cells) is shown in red.

Fig. 2. Evolution of magnetic field strength (derived from volume-weighted magnetic energy densities) in radiative (upper panels) and non-radiative (lower panels) simulations within R200. The left panels focus on Halo3; the right panels show the magnetic field strength across four halos. The magnetic field amplification is faster in radiative simulations. Its evolution in the entire ICM closely mirrors that of the ISM across all 4 halos To show this in more detail, we mark 4 phases of magnetic field evolution in the radiative simulation of Halo3: I) the collapse of the first galaxies, II) the amplification in the BCG, III) the amplification in the host cluster by mergers that contribute their magnetically pre-enriched gas to the host halo, IV) low redshift mergers that constantly inject turbulence, which drives a fluctuating small-scale dynamo.

In the right panels, we show the magnetic field in the ICM averaged across the virial radius for all four halos with solid lines. Dashed lines show the volume-averaged root mean square magnetic field strength in the ISM. The initial decay of the magnetic field strength at high redshift (shown in gray) reflects adiabatic expansion in an isotropically expanding Universe while conserving magnetic flux (B ∝ (1 + z)²). We measure this in the entire high-resolution region. At around z = 20 the protocluster starts to collapse and adiabatically compresses the field. The black dotted lines indicate the expected power-law fits to the adiabatic expansion before the protocluster collapse.

Substantial differences between both simulation setups emerge during the high-redshift interval 20 ≳ z > 3. We highlight four key stages in the magnetic field evolution within the radiative simulation of Halo3 (upper left panel).

1. 20 ≳ z > 8.5: This phase marks the collapse of the first galaxies, which adiabatically compresses the initial seed magnetic field. Structure formation shocks introduce turbulence, which amplifies the magnetic field via a small-scale dynamo. This works very efficiently in the central, cool regions with short eddy turnover times.

2. 8.5 > z > 5.5: At these high redshifts, the protocluster only consists of a few galaxies – the progenitors that eventually merge to form the BCG. Supersonic cooling flows penetrate deep into the central regions and into star-forming galaxies, and inject turbulence. This leads to an exponential magnetic field growth in the cold ISM driven by the gravitational dynamo (we demonstrate this in Sect. 4). Stellar feedback produces additional turbulence, mainly within R ≲ 0.5R₂₀₀, and efficiently mixes the magnetically enriched gas in the galaxies with the ICM (shown in Sect. 4.3.1). As a result, the radiative simulation displays strong magnetic growth in the ISM, followed by the inner radial bin, and less pronounced growth in the outer radial bin.

3. 5.5 > z > 3: This phase is marked by the first series of major mergers. Subsequent magnetic field enrichment results from mergers, contributing pre-enriched gas and turbulence to the ICM, leading to exponential growth (as shown in Sect. 4.4).

4. 3 > z: Left alone, the magnetic field would relax into a minimum energy configuration (Taylor 1974). The accretion of primordial (less) magnetized gas from the intergalactic medium would furthermore cause the magnetic field to decrease as a result of adiabatic cooling and mixing. Instead, merging structures at lower redshifts continuously inject turbulence and drive a small-scale dynamo (see Sect. 4.4). This keeps the magnetic field strength relatively constant in both the ICM and ISM.

The evolution in the non-radiative simulation (lower left panel) differs significantly. The magnetic field strength grows at a reduced rate at high redshifts. In phase I, the central regions are less dense (we show this in Fig. A.1) due to the lack of cooling, which reduces the efficiency of the dynamo. During phase II there are no central galaxies to develop and mix a strong magnetic field. Magnetic field growth only continues at an increased rate after the third phase, marked by major merger events, but still at a relatively low growth rate, compared to the radiative simulation. This results from the lower initial magnetic fields available for a dynamo, as there are no merging galaxies to provide their magnetic fields as an effective seed field. Notably, magnetic field growth accelerates only after z ≈ 3, and saturates at z ∼ 1.

In the right panels of Fig. 2, we examine the magnetic field evolution in all four clusters. Magnetic field growth is more rapid in the radiative simulations, saturating at around z ≈ 4 in the higher mass clusters, Halo3 and Halo4, and at around z ≈ 3 in the lower mass clusters, Halo194 and Halo260. The saturation field strengths are 1 μG volume-averaged within R₂₀₀ in all four clusters. Notably, magnetic field evolution in the ICM of the radiative simulations closely mirrors that in the ISM, indicating a strong connection between the two. The ISM field saturates at around 20 μG in all halos. The ISM field in all halos saturates before the ICM. The ISM field in the most massive halo, Halo3 saturates slightly earlier in comparison to the other halos. In contrast, the non-radiative simulation achieves magnetic field saturation later, at around z ≈ 2 in the higher mass halos and around z ≈ 1 in the lower mass halos. The saturation field strengths are 1 μG in the higher mass clusters and a factor of two lower in the smaller clusters.

In Fig. 3, we address the important topic of numerical resolution on the evolution of the magnetic field strength. We compare two different resolution levels, ‘zoom8’ (8³ = 512 times better mass resolution compared to the parent simulation), and ‘zoom12’ (12³ = 1728 times better mass resolution compared to the parent simulation). As the numerical Reynolds number scales with the inverse cell size, this would give an eight times higher Reynolds number in the zoom8 simulation, compared to the parent simulation and a twelve times higher Reynolds number in the zoom12 simulation. The growth rate for the dynamo, in turn scales as Γ ∝ Re^1/2 for incompressible turbulence, which would result in a 2.8 (3.5) times higher growth rate in the zoom8 (zoom12) simulations, respectively. We find no difference between the magnetic field strength at different levels of resolution for adiabatic compression of the magnetic field due to protocluster collapse at z ∼ 20. In our radiative simulations, the saturated field is well converged while magnetic field growth is only somewhat slowed down in the lower resolution simulation during the dynamo growth phase. By contrast, in the non-radiative simulation the smaller numerical Reynolds number in combination with the finite time that a sizeable amount of turbulence is injected during merger events causes significantly slower dynamo growth. While the saturated magnetic field is just converged in the two massive clusters at z = 0, it fails to converge in the two smaller clusters. In the following sections, we investigate the different growth phases in the radiative simulation in more detail, using Halo3 as our main example. To quantify our results, we show fit growth rates (see equation 8) for all 16 simulations in Table 2. We determine the time intervals for the fit visually. A more detailed study of the exponential magnetic growth in Halo3, and a comparison of analytic and fit growth rates will be presented in Sect. 4.3.3.

Fig. 3. Resolution study of the evolution of the magnetic field strength averaged within R200 for the radiative (left) and non-radiative (right) simulations. We present the evolution for all four halos (shown in different colors) and compare two resolution levels: “zoom8” and “zoom12” (see text). The radiative simulation shows no significant differences between the two zoom levels for all four halos, reflecting the fact that the magnetic field is amplified predominantly in the ISM, which is well resolved in both resolutions. In the non-radiative simulation, significant differences emerge between the two zoom levels, with zoom8 simulations showing smaller growth rates and lower final field strengths for the two smaller clusters.

The rate of magnetic field growth in a small-scale dynamo depends on the amount of available kinetic (turbulent) energy. In Fig. 4, we show the gas-relaxation parameter 2E_kin, gas/|E_pot| for all four halos in the radiative simulation. We compute the gas kinetic energy E_kin, gas in the cluster rest system within R₂₀₀. We compute the gravitational potential energy by summing over all particles within the cluster (gas, stars, and dark matter), treating the cluster as an isolated system and assuming Newtonian gravity with the potential normalized to zero at infinity. For this, we used a gravitational tree code that includes only particles within R₂₀₀:

Fig. 4. Evolution of the gas-relaxation parameter 2Ekin, gas/|Epot|, which quantifies the dynamical state of the ICM for all four halos in the radiative (left) and non-radiative (right) simulations. In the radiative simulations, the high-redshift phase is dominated by frequent mergers and strong radiative gas cooling, which results in accretion and feedback events. In consequence, the ICM finds itself in an unrelaxed state in the forming protocluster. The non-radiative simulations display lower values of 2Ekin, gas/|Epot| at high redshift due to the lack of feedback and cooling. Mergers at late times can be seen as upwards fluctuations in this quantity in both simulation runs.

$$\begin{array}{c}\hfill |E_{\mathrm{pot}}|=-{\displaystyle \sum _{\begin{array}{c}i,j=1\\ i>j\end{array}}^N}\frac{G{m}_i{m}_j}{|{\mathit{r}}_i-{\mathit{r}}_j|}.\end{array}$$(3)

The total relaxation parameter 2E_kin/|E_pot| is commonly defined as twice the kinetic energy of all components (dark matter, gas, stars and black holes) divided by their total potential energy. It serves as an indicator of the clusters dynamical state, where 2E_kin/|E_pot|> 1 indicates that the cluster is dynamically unrelaxed as it is the case after a merger event, 2E_kin/|E_pot| = 1 reflects a virialized cluster (with a negligible effective pressure due to infalling matter at the cluster boundary, The & White 1986), and 2E_kin/|E_pot|< 1 indicates a system dominated by potential energy, such as the initial stages of collapsing structures⁴. Here, we show specifically only the gas-relaxation parameter 2E_kin, gas/|E_pot| to emphasize the differences between the ICM kinetic energies⁵ of the radiative/non-radiative simulations of the same halo because the potential energies are approximately the same in both types of simulations. Note that for a relaxed massive cluster one would expect that 2(E_kin, gas + E_therm, gas + E_B)/|E_pot|≈Ω_b0/Ω_m0 ≈ 0.15, where E_therm, gas is the thermal energy of the gas, and E_pot is the total potential energy including all matter species, as above. Figure 4 demonstrates that the gas-relaxation parameter, 2E_kin/|E_pot|, is at most about half of that and smaller than 0.08 for all redshifts and simulations. This is because, upon passing through the cluster accretion shock, a significant fraction of the gas kinetic energy is dissipated and converted into thermal energy, lowering E_kin, gas. Kinetic energy is then further reduced by dynamo action. The comparison in Fig. 4 shows that the kinetic energy in the radiative simulations is significantly increased at high redshift in comparison to the non-radiative simulations. This is a result of cooling-induced infall and kinetic energy injection as a result of stellar and AGN feedback.

3.3. Magnetic pre-enrichment in galaxies

As demonstrated, rapid magnetic field growth takes place in galaxies. Two stages are to be distinguished with regards to this magnetic pre-enrichment: We explore the influence of pre-enrichment in the BCG (phases I and II) and in merging substructures (phase III and IV). This provides an interesting view of the magnetic field growth in galaxy clusters, showing that magnetic fields are intrinsically linked to its structure formation processes.

We show this in Fig. 5, where we examine the time evolution of the spherically averaged radial profiles of the mass-weighted vorticity (left), the mass-weighted metallicities (middle) (indicative of stellar activity), and the volume-averaged root mean square magnetic field strength (right) of Halo3. In the upper panel, we focus on merger activity, showing the region 3R₂₀₀(z) across the redshift range 20 > z > 0.1. Dashed lines denote the virial radius (R₂₀₀) and the central cluster region (0.25 R₂₀₀), both in physical units. The lower panels focus on the redshift range 11 > z > 4 to emphasize the early enrichment by metals and magnetic fields in the central galaxy and out to a fixed physical radius of 150 kpc, which corresponds to twice the virial radius of the central galaxy at z = 6.6.

Fig. 5. Redshift evolution of the physical radius corresponding to 3R200 (≈9 Mpc) surrounding the cluster in the radiative simulation of Halo3. From left to right, we color code the spherically averaged fields of vorticity, metallicity, and magnetic field strength. The upper row shows the redshift range 20 > z > 0.1, highlighting the late time evolution dominated by mergers, while the bottom row shows the redshift range 11 > z > 4, focusing on the evolution in the central galaxy at higher redshifts. The dashed lines indicate the virial radius (R200) and the central cluster region (0.25R200). We keep the maximum radius in the lower panels at a fixed physical value of 150 kpc (see text). The lower panels show how vorticity and metallicity are enhanced in the center already at high redshifts, z > 9. The central magnetic field amplification sets in shortly afterwards at z ≈ 7. Turbulence from gravitational infall of cool gas amplifies the magnetic field and drives a small-scale dynamo. Stellar feedback additionally injects turbulence, distributes metals, and mixes magnetically enriched gas with less enriched gas. The upper panel highlights how merging substructures contribute their magnetically enriched gas to the host cluster. These merging structures also exhibit higher metallicity, indicating the connection between strong magnetic fields and the metal-rich ISM.

First, we examine magnetic field growth in the central galaxy at higher redshifts in the lower panels. Initially, at z > 8, magnetic field strength, metallicity, and vorticity have a low amplitude. At approximately z ≈ 9, metallicity and vorticity notably increase to Z = 1 Z_⊙⁶, and ω = 10⁻¹ Myr⁻¹ in the central region, where the BCG is located. A strong magnetic field B ≈ 1 μG emerges slightly later at around z ≈ 7 in the central region. This shows dynamo activity with a slight time lag following the onset of turbulence. Supersonic inflows of accreted gas generate turbulence in the cold, dense filaments, amplifying the magnetic field. Stellar feedback further interacts with the gas, distributing metals and creating additional turbulence, which mixes magnetized and less magnetized gas at larger radii. This leads to inside-out growth of the magnetic field, with initial amplification in the cluster center, followed by gradual enrichment of the outer regions.

The important role of merger activity is evident in the upper panels of Fig. 5. Mergers appear as regions with increased magnetic field, metallicity, and vorticity beyond the virial radius, subsequently entering the cluster. Merging substructures exhibit high magnetic field and metallicity, emphasizing the role of stellar activity in redistributing the gas, and highlighting the correlation between metals and magnetic fields in intergalactic space. Mergers contribute their pre-enriched gas from within galaxies to the ICM of the main cluster⁷. In Appendix A, we present radial profiles showing how gas pre-processing in merging substructures enriches the cluster outskirts in both magnetic field strength and metallicity.

Presumably, also ram pressure stripping plays an important role in removing the gas from infalling galaxies. This happens already at large radii, beyond the virial radius, as ∼90% of all galaxies arrive at the virial radius without any self bound gas at lower redshifts. Although the majority of the galaxies fulfills the Gunn et al. (1972) criterion for ram pressure stripping to occur, this picture is probably too simplistic and can only give an upper limit because of insufficient numerical resolution to fully resolve this process (Sparre et al. 2024b) and because of the pressurized effective equation of state, which responds differently to ram pressure in comparison to a multi-phase ISM.

Pre-processing of metals and quenching of galaxies has been studied in observations and simulations of clusters and their vicinities, highlighting the role of ram pressure stripping and galactic winds (Butcher & Oemler 1984; Fujita 2004; Tonnesen & Bryan 2012; Boselli et al. 2014; Ebeling et al. 2014; Bahé & McCarthy 2015). This process is observed not only near the cluster center but also at large radii, and prior to merging (Fujita 2004; Bahé et al. 2013; Haines et al. 2015). From Fig. 5 it is evident that the pre-processing also plays an important role for the magnetic field enrichment in galaxy clusters.

4. Magnetic field growth in the ISM and in the ICM

In this section, we distinguish between magnetic field growth in the ISM and ICM. At high redshifts, we show that a dynamo amplifies the magnetic field in the BCG and other early galaxies through compressible turbulence. Galactic winds transport and mix this magnetically pre-enriched gas into the ICM. In the ICM, a small-scale fluctuating dynamo driven by vortical turbulence further amplifies the magnetic field in the pre-enriched gas at low redshifts.

4.1. Overview of fluctuating dynamos

First, we provide an overview of the magnetic small-scale dynamo and the emergent scaling laws. A fluctuating small-scale dynamo operates as magnetic field lines stretch, twist, fold, and merge to amplify weak seed magnetic fields (Brandenburg & Subramanian 2005; Donnert et al. 2018; Rincon 2019). Such a dynamo operates in initially weakly magnetized, turbulent media like the ICM (Wittor et al. 2017; Vazza et al. 2017; Donnert et al. 2018; Wittor & Gaspari 2020) and exhibits exponential growth.

One sign of dynamo activity is the emergence of specific scaling laws that link magnetic power to wavenumber. In the kinematic phase of the dynamo, assuming an incompressible velocity field, Kazantsev (1968) demonstrated that magnetic power scales as P(k)∝k^3/2 above the resistive scales. When the magnetic field becomes strong and back-reacts on the small scales, (see, e.g., Vainshtein & Cattaneo 1992 and Zweibel & Yamada 2009), larger eddies shift the power spectrum peak to larger scales, following Kolmogorov scaling for incompressible turbulence (Kolmogorov 1941), P(k)∝k^−5/3, on scales below the power spectrum peak. By doing so, also the largest magnetic correlation length of the dynamo amplified field shifts to larger scales. Theories and simulations show that the Kazantsev slope shifts to larger scales in this so-called non-linear regime (Schober et al. 2015).

Studies have demonstrated that dynamo activity is also possible in compressible turbulence (Kazantsev et al. 1985; Brandenburg & Subramanian 2005; Schober et al. 2012b, 2015). Supersonically induced turbulence on large scales decays to incompressible turbulence below the sonic scale (Fouxon & Mond 2020; Federrath et al. 2021; Beattie et al. 2023). This can be understood such that compressible velocity fields caused by, for example, oblique shocks produce a non-negligible component of vortical turbulence (see, e.g., Käpylä et al. 2018a,b, who investigated turbulence caused by clustered SN bubbles). This implies that the growth still takes place in Kolmogorov-type turbulence (with the addition of shock compression of magnetic fields), such that Kazantsev’s theory is still applicable. According to simulations (see Brandenburg & Subramanian 2005 and references therein), a dynamo can develop, exhibiting Kazantsev scaling on large scales and Burgers scaling for compressible turbulence (Burgers 1969), P(k)∝k⁻², on smaller scales. However, the saturation strengths and growth rates are found to be lower, in comparison to Kolmogorov turbulence (Brandenburg & Subramanian 2005; Schober et al. 2012b; Federrath 2016). We will use this information to compare our results to the scaling laws described above in order to identify dynamo activity and in order to differentiate between dynamos induced by vortical or by compressive turbulence in Sect. 4.2.

The final saturation levels of the fluctuating small-scale dynamo are still not fully understood. Saturation is caused by the Lorentz force backreacting on the flow. This can be understood such that with increasing magnetic field strength, the backreaction of the Lorentz force onto the gas increases, causing an increasing velocity drift that lowers the effective magnetic Reynolds number (Subramanian 1999). When it falls below the critical Reynolds number (Re_c ≈ 60 − 120) and the dynamo activity comes to a halt (Subramanian 1999; Seta et al. 2020; Achikanath Chirakkara et al. 2024). It has been furthermore shown that the final saturation values of the magnetic-to-turbulent energy density depend on the magnetic Prandtl number and on the form of turbulence, where higher Prandtl numbers, as well as turbulence in vortical form yield higher saturation values, and lower Prandtl numbers and turbulence in compressive form yield lower saturation values (Schober et al. 2015). We are working in the Pr_m ≈ 1 regime, as both our viscosity and magnetic diffusion derive from numerical dissipation effects at approximately the scale of our numerical Voronoi cells. The effects that stem from the different forms of turbulence (solenoidal or compressive), will be discussed in the upcoming sections.

There is ongoing debate about whether the dynamo always saturates at the same fraction of available kinetic energy (see, e.g., Miniati & Beresnyak 2016). This idea stems from the fact that most kinetic and turbulent energy in the ICM arises from accretion and merger events. If the gravitational potential energy from a merger is converted into bulk and turbulent kinetic energy, the available kinetic energy should scale with the cluster mass. However, magnetic energy depends on density profiles and dynamo activity. Since density profiles in clusters are believed to be self-similar (McDonald et al. 2017), magnetic fields may be comparable due to adiabatic gas compression. Dynamo activity, on the other hand, is driven by the amount of turbulent kinetic energy, meaning the ratio of magnetic to kinetic energy should reflect the dynamo’s efficiency and be roughly universal across clusters (Miniati & Beresnyak 2015). We will investigate the efficiency of the dynamo amplification in our simulations in Sect. 4.4.3.

4.2. Power spectra for different regimes of magnetic field growth

To differentiate between dynamos driven by compressible or vortical turbulence in the high/low redshift regime, we analyze power spectra. We have computed power spectra for Halo3 in order to connect its magnetic field evolution with the theory and literature outlined above. In Fig. 6, we show the magnetic and kinetic power spectra across different redshifts (8.5 > z > 0) for Halo3. The power spectrum is computed from a zero-padded fast Fourier transform of the components of $\sqrt{\rho }\mathbf{v}$ and $\mathbf{B}/\sqrt{4\pi }$. Calculations are performed in a sphere with a physical radius of 20 kpc (which corresponds to the virial radius of the central galaxy at z = 9.5) in the left panels, illustrating magnetic growth in the BCG, and a radius of 0.5R₂₀₀ in the right panels, excluding the central 20 kpc. We show the radiative simulation in the upper panels and the non-radiative simulation in the lower panels. Dotted lines represent various scaling behaviors: Burgers’ model (∝k⁻²) for compressible turbulence, Kolmogorov’s model (∝k^−5/3) for incompressible turbulence, the Kazantsev model at large scales (∝k^3/2) indicating dynamo activity, and a Lorentzian profile (that scales with ∝k² on large scales) as an indicator of an exponential magnetic field profile in real space.

Fig. 6. Kinetic (dashed) and magnetic (solid) power spectra of Halo3. The left panels show the power spectra in the central regions within a fixed physical radius of R200(z = 9) = 20 kpc for four different redshifts (z = 4.5, 5.5, 7, 8.5) for both the radiative (upper) and the non-radiative (lower) simulations aimed at tracing the high-redshift dynamo in the BCG. The right panels show the ICM power spectra within 0.5R200(z), excluding the central 20 kpc, at six different redshifts (z = 0, 1.5, 2.5, 3.5, 4.5, 6.5) focusing on dynamo activity in the ICM within the radiative and non-radiative simulations. Dotted lines illustrate different scaling behaviors (see text). The left panels show that a small-scale dynamo (with the Kazantsev scaling ∝k3/2) driven by compressive Burgers turbulence (∝k−2) is only active in the central regions of the radiative simulation. By contrast, the right panels show the existence of a small-scale dynamo driven by incompressible turbulence (∝k−5/3) throughout the entire ICM within 0.5R200 in both simulations, albeit at a faster rate in the radiative simulation. The contribution of a Lorentzian profile, which we show for z = 4.5 in both simulations and additionally z = 0 in the radiative simulation (it scales as ∝k2 on scales larger than the scale height), at lower redshifts is minimal because the sharp cutoff at the characteristic exponential scale is located at rather large scales.

First, we analyzed magnetic power spectra in the central cluster regions at high redshifts (8.5 > z > 4.5) in Fig. 6 (left panels). Generally, magnetic energy in the radiative simulation exceeds that of the non-radiative simulation at all redshifts. At z = 8.5, the magnetic energy remains low in both simulations. At z = 6.5, the magnetic energy begins to grow on small scales in the radiative simulation, with the Kazantsev slope emerging on larger scales. The first peak of the magnetic power spectrum occurs at a characteristic scale of 2 kpc in the radiative simulation and shifts to larger scales of 10 kpc by z = 4.5. At lower redshifts, the wavenumber scaling becomes somewhat steeper than the Kazantsev slope. On smaller scales, the spectrum follows the Burgers slope for compressible turbulence in both simulations. The compressible turbulence in the radiative simulation is due to gravitational, supersonic accretion of cold gas. In the non-radiative simulation the compressive turbulence is also due to cosmological accretion and mergers that move supersonically towards the center, but without cooling, which barely drives dynamo activity, as the eddy turnover times are larger. The turbulent energy is injected at larger scales (ℒ ≈ 10 kpc) as seen in the kinetic power spectra.

The right panels of Fig. 6, showing the entire ICM within 20 kpc < r < 0.5R₂₀₀ over a broader redshift range 6.5 > z > 0, demonstrate that magnetic energy increases over time, with its peak shifting to larger scales as expected in both simulations. In the non-radiative simulation, the magnetic power grows more continuously, whereas in the radiative simulation, it grows rapidly between z = 8.5 and z = 6.5. While the power spectrum follows the Kolmogorov scaling ∝k^−5/3 on small scales, it is ab initio unclear whether the cluster magnetic field profile or the Kazantsev scaling of the fluctuating dynamo sets the large-scale behavior of the power spectrum (cf. Fig. 7 of Pfrommer et al. 2022).

We show in Appendix A that the radial profiles of the magnetic field strength in the radiative and non-radiative simulations follow (double) exponential profiles, which implies (double) Lorentzian profiles in Fourier space. The associated one-dimensional power spectrum for a single-exponential profile is given by (equation (A.4))

$$\begin{array}{c}\hfill P(k)=4\pi k^2{P}_{3\mathrm{D}}(k)\propto \frac{{\left(B_0{k}_{0}k\right)}^2}{{(k^2+k_0^2)}^4},\end{array}$$(4)

where k₀⁻¹ is the characteristic scale length of the exponential profile. This one-dimensional power spectrum scales ∝k² for small wave numbers and ∝k⁻⁶ for large wave numbers. This implies that the contribution of the steeply rising Lorentzian form factor to the power spectra that we measure in our simulations on large scales could be flattened to the Kazantsev scaling ∝k^3/2 below the characteristic scale of the exponential profile because of the sharp drop-off of the squared Lorentzian towards small scales. Most importantly, we expect the Kolmogorov spectrum to dominate for incompressible turbulence on scales smaller than the magnetic coherence length, defined as

$$\begin{array}{c}\hfill {\lambda }_B=\frac{\int 2\pi k^{-1}P(k)\mathrm{d}k}{\int P(k)\mathrm{d}k}.\end{array}$$(5)

In Fig. 6 we show the Fourier-space double Lorentzian profile of the magnetic field in the radiative simulation at z = 4.5, a single Lorentzian profile in the non-radiative simulation at z = 4.5, and a single Lorentzian profile in the radiative simulation at z = 0. The corresponding (double) exponential fits are shown in Appendix A.

Figure 6 shows that over time the magnetic coherence scale grows in both simulations. On scales smaller than the exponential scale length (but still in the kinematic limit, i.e., on scales larger than the magnetic coherence scale), the logarithmic power spectrum slope is consistent with the Kazantsev spectrum, suggesting the emergence of an active small-scale dynamo in the ICM for redshifts z < 3.5 in the radiative simulation and for all measured redshifts in the non-radiative simulation. By z = 0, the exponential scale length of the magnetic profiles is two orders of magnitude larger than the magnetic coherence length, so that we can see a clear Kazantsev scaling on scales smaller than the magnetic coherence length in both simulations.

In the radiative simulations for 6.5 ≳ z ≳ 3.5, the magnetic field strength profile follows a double exponential structure, which has an inner and outer characteristic scale. As shown in Fig. A.2, for z ≳ 3.5, the inner scale of the magnetic profile is comparable to the magnetic coherence length λ_B as defined in equation (5). However, we argue that the magnetic coherence scale λ_B, the characteristic scale of the exponential 2π/k₁, and the outer scale of the turbulence ∼2π/k_peak lie close together, so that the peak in the power spectrum is not only a feature of the inner magnetic profile. If it were, we would expect a discontinuity or bump-like features in the measured power spectrum, arising from the superposition of two distinct components – the Lorentzian profile and the turbulent spectrum – each with its own characteristic peak. Instead, the observed power spectra show that the Lorentzian profiles smoothly connect to the Kolmogorov scaling expected for incompressible turbulence, without introducing any noticeable discontinuities.

Notably, the injection scale of turbulence measured in the analyzed sphere is approximately the virial radius corresponding to the given time (ℒ ≈ R₂₀₀). Most of the turbulent energy comes from mergers that excite supersonic velocities beyond the virial radius and transitioning to subsonic turbulence when entering the virial radius, minor contributions are stemming from stellar and AGN feedback.

4.3. Magnetic field growth in the ISM

At high redshifts, z > 5.5, most of the ICM gas is star-forming, with magnetic field growth primarily occurring in the ISM. In this section, we demonstrate that the turbulence required for the dynamo process in the ISM is driven by supersonic cooling flows, while stellar and AGN feedback generate magnetized outflows that effectively enrich the remaining cluster volume with magnetic fields.

4.3.1. Gravitationally driven dynamo and mixing induced by stellar and AGN feedback

Figure 7 presents thin projections from the radiative simulation of Halo3, showing gas speed, thermal sound speed, Mach numbers, star formation rate surface density, vorticity, the ratio of free-fall time to the eddy-turnover time, particle number density, magnetic field strength, and metallicity at z = 6.2. This analyzed time lies in phase II, during which the magnetic field is mainly amplified in the BCG and then partially mixed with the surrounding ICM. The free-fall time t_ff is calculated in N = 50 radial annuli with spacing Δr, following the expression $t_{\mathrm{ff}}(r)=\sqrt{3\pi /(32\mathrm{G}{\rho }_{\text{mean}}(r<(R_{200}-N\mathrm{\Delta }r)))}$ which assumes that the gas collapses from R₂₀₀ − NΔr to r = 0 starting with an initial velocity of 0 and assuming a mean density that takes into account dark matter, gas and stars. This approach does not account for gas that has already fallen into the cluster and has an initial, non-zero velocity, which leads to an overestimation of t_ff at smaller radii, as the enclosed mass inside R₂₀₀ − NΔr decreases significantly. From this panel, we observe that the ratio remains close to unity, particularly in the outskirts of the cluster. This indicates that turbulence in these regions is primarily driven by gravitational infall. In contrast, regions near the cluster center, where the magnetic field is stronger, exhibit ratios around ∼10, which is due to our $v(t=0)=0$ assumption. In these central regions, magnetic field amplification is more efficient due to higher gas densities and correspondingly shorter eddy turnover times. This effect is consistent with the findings of Irshad et al. (2025), who show that magnetic growth is accelerated in collapsing systems owing to enhanced eddy turnover rates. It is furthermore evident from the velocity panel that there is a gas inflow towards the center. The gas flows along dense, star-forming filaments with a high velocity, which creates strong shocks in these cold filaments. Shocks encountering density inhomogeneities produce a high amount of vorticity. Due to the high densities, the resolved eddy turnover rates are very high, ranging from ω = 0.1 Myr⁻¹ − 1 Myr⁻¹ in the dense regions, as evident from the middle panel, which enables efficient dynamo amplification.

Fig. 7. Thin projections through the center of Halo3 in the radiative simulation with sidelength 2.2R200 at z = 6.2 and thickness 0.5R200. From top left to bottom right, the nine panels display the gas velocity in the cluster rest frame, sound speed, Mach number, star formation rate surface density, vorticity, the ratio of the free-fall time to the eddy-turnover time, particle number density, magnetic field strength, and metallicity. The velocity field is shown together with its absolute value in the first panel. Regions with cold, supersonic gaseous inflows create vorticity and display strong magnetic fields. Star formation in the high density regions creates stellar feedback driven winds that interact with the inflowing gas and produce vorticity. This vorticity efficiently mixes highly magnetized gas with weaker magnetized gas, distributing metals across the cluster.

Regions with a high star formation rate launch stellar feedback driven galactic winds which impact the inner cluster region (r ≲ 0.5R₂₀₀). These produce turbulence and can interact with the infalling gas, as can be seen in the upper left corner of the velocity panel, close to the virial radius. Additionally, the gravitational potential of the forming cluster accretes gas that forms a shock close to the virial radius upon encountering the dense ICM. The inhomogeneous density distribution causes curved virial shocks, which generates vortical turbulence according to Crocco’s theorem. The turbulence stemming from both sources can efficiently mix stronger magnetized gas from lower radii with less magnetized gas at higher radii due to the shallow entropy profile in the cluster outskirts (see the middle right panel and the shallow entropy profiles at larger radii in Appendix A.1). Metals are widely distributed throughout the entire ICM which is a testament to the effective turbulent mixing due to galactic winds, as shown in the lower right panel.

We investigate the main energy injection sources into the gas to identify drivers of galactic winds that mix and distribute high magnetic field gas. In Fig. 8, we show the energy injection rates of different feedback channels in the radiative simulation of Halo3: stellar feedback, thermal AGN feedback, kinetic AGN feedback, and compare this to the energy loss rate due to cooling. In the left panel of Fig. 8 (representative of the BCG), we display the energy injection rates that are launched from the central 20 kpc, while the right panel shows the energy injection rates that are launched from the remaining part of the cluster (excluding the center) within the virial radius. Our calculations follow the methodology outlined in Weinberger et al. (2018) (see their section 2.3 for AGN feedback and their section 2.4 for stellar feedback, the models used to simulate AGN and stellar feedback are described in Sect. 2.2).

Fig. 8. Energy injection rates into gas cells of Halo3 in the radiative simulation are shown. The panels differentiate between feedback launched from the BCG center (r < 20 kpc; left panel) and from the rest of the cluster, excluding the center (right panel). We present different energy injection rates that originate from stellar feedback (yellow) as well as thermal (purple) and kinetic (orange) AGN feedback. The injected thermal AGN feedback (dashed purple) is numerically reduced in star-forming cells where an effective equation of state is enforced, implying a rapid energy loss and a smaller effective thermal AGN feedback energy rate (solid purple). Radiative cooling losses (blue) are displayed for the ICM in the right panel. At high redshifts, both in the BCG and the entire cluster, the primary energy injection rate is due to stellar feedback. Subsequently, AGN feedback also becomes a significant contributor.

As observed in the left panel, stellar feedback serves as the primary energy injection source launched from the center until approximately z ≈ 6. Subsequently, there is a brief phase (5 > z > 3) where also thermal AGN feedback becomes a main contributor while kinetic AGN feedback starts to dominate the energy injection after z = 3. Examining the right panel of Fig. 8, we discern two distinct phases of feedback injection in the entire ICM: at z > 5, stellar feedback is the primary energy source, while at z < 5, contributions from AGN, both thermal and kinetic, dominate the energy injection. The feedback acts in such a way that it removes gas from galaxies and mixes it into the ICM. Apart from the feedback channels listed above, cluster galaxies, being located in high-density regions, experience gas loss due to ram pressure stripping during and after phase III (see Sect. 3.3).

4.3.2. Adiabatically expanding galactic magnetic fields into the ICM

We demonstrate that as galaxies fall into the cluster potential, processes that transport gas from galaxies into the ICM, such as magnetized outflows and ram pressure stripping, can significantly contribute to the magnetization of clusters. We do this by demonstrating that the simulated magnetic field strength in the ICM can be largely explained by adiabatically expanding the magnetic field from galaxies. We first calculate the volume of the ISM, V_ISM, and the mean magnetic energy inside of it, ε_ISM. We then imagined that the ISM volume is adiabatically and isotropically expanded to fill the ICM volume while conserving magnetic flux. Because B ∝ ρ^2/3 during this process, the magnetic energy scales as ε_B ∝ B² ∝ ρ^4/3 ∝ V^−4/3, and we have ε_BV^4/3 = const. We can use this relation to make a prediction for the mean magnetic energy density in the ICM if it has been purely generated as a result of adiabatic expansion of the ISM:

$$\begin{array}{c}\hfill {\epsilon }_{\text{ICM},\text{expanding}\text{ISM}}={\epsilon }_{\mathrm{ISM}}{\left(\frac{V_{\mathrm{ISM}}}{V_{\mathrm{ICM}}}\right)}^{4/3}.\end{array}$$(6)

We can also calculate the magnetic energy as a result of compressing/expanding the pre-collapsed protocluster magnetic field and density to its current density, which yields

$$\begin{array}{c}\hfill {\epsilon }_{\text{ICM},\text{collapsing}\text{cluster}}={\epsilon }_{B,0}{\left(\frac{\rho }{{\rho }_0}\right)}^{4/3},\end{array}$$(7)

where ε_B, 0 and ρ₀ are the mean magnetic energy density and density at z = 15, respectively.

Figure 9 illustrates the evolution of magnetic field energy density between z = 15 and 4.5 for Halo3. Solid lines represent the measured mean energy density with a volume (as explained below) while dashed lines show the expectations as a result of adiabatically changes of the gas from the initial magnetic field strength and gas density after the protocluster collapse (equation (7)). The radius of the measurement volume expands as the halo’s virial radius. Once the energy density within this radius exceeds the expectations from adiabatic compression for the first time by a factor of 10², we set the virial radius at this time as a fixed radius for further measurements, following the method of Pakmor et al. (2024). This radius is R₂₀₀ = 74 kpc at z ≈ 6.5. The left panel displays magnetic energy densities for the ISM (red), while the right panel shows the magnetic energy density for the ICM (blue). Additionally, we represent the adiabatically expanded magnetic field from the ISM in the right panel with dotted red lines (equation (6)).

Fig. 9. Left panel: Magnetic energy density in the ISM (red). We measured the energy density in an expanding sphere with a radius equal to the virial radius at the time. When the magnetic energy density inside this sphere grew for the first time by a factor of 102 compared to the adiabatically compressed field right after the protocluster collapse, we fixed the radius at this redshift. Right panel: Magnetic energy density in the ICM (blue) compared to the magnetic energy density expected from expanding the ISM magnetic field purely adiabatically to the volume of the ICM (red dotted). This estimate even overshoots the total magnetic energy in the ICM. The lower measured ICM magnetic field strengths, compared to adiabatically expanded ISM fields, may result from depletion during star formation, retention in the ISM, or unaccounted mixing losses.

From the left panel of Fig. 9, it is evident that a dynamo is operational, causing the magnetic field to surpass the expectations from adiabatic compression by multiple orders of magnitude and amplifying the magnetic field exponentially. In the right panel of Fig. 9, the energy density in the ICM similarly exceeds that achievable through adiabatic compression alone by several orders of magnitude. However, expanding the ISM magnetic field into the ICM overshoots the actual values in the ICM⁸. This discrepancy can arise because this estimate does neither account for depletion of the magnetic field during star formation events nor for the magnetic field that is retained in the ISM. Note that this argument assumes full mixing of the adiabatically expanded magnetic field with the weaker magnetized ICM as well as no magnetic dissipation during this process, which will also lead to estimates being higher than the true value.

Clearly, all these effects can decrease the magnetic field strengths in the ICM below the measured values, thereby requiring a small-scale dynamo to amplify this already pre-enriched magnetic field (see Sect. 4.4). In any case, the pre-enrichment scenario obviates the need for high dynamo growth rates in the ICM to explain a strong magnetic field at high redshift. Interestingly, the adiabatic evolution of the initial magnetic field shows little evolution as adiabatic dilution due to cosmic expansion is approximately balanced by adiabatic compression due to gravitational collapse.

4.3.3. Growth rates in the ISM

Figure 10 shows the evolution of the magnetic energy density in the ISM of all four halos (left panel). The contributions due to adiabatic compression of the gas have been linearly subtracted in order to isolate the amplification via dynamo activity. A fluctuating small-scale dynamo shows a characteristic exponential growth:

Fig. 10. Left: Time evolution of the magnetic energy density with fits of the exponential growth rates for the ISM in all four halos (dotted lines). The fits started when the magnetic energy density for the first time grew to a value at or above 4 × 10−7 erg cm−3 and ended when the magnetic energy density for the first time reached a value of at least 40 erg cm−3. We subtracted the contributions from adiabatic compression following the relation given in equation (7), where we fixed B0 and ρ0 to be the measured values after the cluster collapse, and εB(t), ρ(t) were measured in each snapshot. Right: Fit growth rates (circles) along with the expected growth rates (crosses) that include error bars as described in the text. Our analytic and measured growth rates agree within a factor of two, with the expected growth rates slightly larger than the measurements.

$$\begin{array}{c}\hfill {\epsilon }_B(t)={\epsilon }_{B_0}{\mathrm{e}}^{\mathrm{\Gamma }t},\end{array}$$(8)

where ε_B0 is the initial magnetic field strength and Γ is the growth rate. We fit magnetic growth rates in the ISM of all four halos (shown as dotted lines) and compare these fit growth rates against analytical growth rates. The fit to the curves is a least square fit minimizing (Δlog ε_B)² and it is performed linearly in time. We restrict the fit to the interval between the time when the renormalized magnetic energy density first reaches a value of 4 × 10⁻⁷ erg cm⁻³ and the time when it first exceeds a value of 40 erg cm⁻³. The fits for Halo3 and Halo260 approximate the overall growth of magnetic energy density in the ISM very well. In contrast, the fits for Halo4 and Halo194 are less accurate due to spikier and irregular growth patterns, likely resulting from feedback events that eject magnetized gas from the ISM or lock it up in stars. Additionally, the ISM likely contains multiple superimposed dynamos (Pfrommer et al. 2022), making it difficult to fit a single growth rate.

In order to compare our measured growth rates against theoretical models, we adopted the equation of Schober et al. (2012c):⁹

$$\begin{array}{c}\hfill \mathrm{\Gamma }=\frac{163-304\theta }{60}\frac{\mathcal{V}}{\mathcal{L}}{\mathrm{Re}}^{(1-\theta )/(1+\theta )},\end{array}$$(9)

where Re is the Reynolds number, 𝒱 is the injection velocity of the turbulence, and ℒ is the injection scale of the turbulence. The parameter θ is characterized by the type of turbulence, with θ = 1/3 for Kolmogorov and θ = 1/2 for Burgers turbulence. Because we encounter compressible turbulence in the ISM, we used θ = 1/2.

We used an injection scale of ℒ = 10 kpc, as observed in the kinetic power spectra. The injection velocity 𝒱 is calculated as the mass-weighted root mean square velocity of all cells in the center of mass frame. This approach, is based on the assumption that turbulence originates mainly from gravitational infall that cascades towards the center, but also from various other sources such as star-forming regions and AGN feedback. We assume that these different sources of turbulence are well mixed in the center and that the mean velocity is somewhat connected to the injection velocity. The numerical Reynolds number can be estimated via (Pfrommer et al. 2022)

$$\begin{array}{c}\hfill {\mathrm{Re}}_{\mathrm{num}}=\frac{\mathcal{L}\mathcal{V}}{{\nu }_{\mathrm{num}}}\sim \frac{3\mathcal{L}\mathcal{V}}{d_{\mathrm{cell}}{v}_{\mathrm{th}}}.\end{array}$$(10)

Here, ν_num is the numerical viscosity, and d_cell is the diameter of the cell which we calculate from the cell volume, assuming the cell is spherical. The velocity $v_{\mathrm{th}}=\sqrt{2k_{\mathrm{B}}T/[\overline{m}(\gamma -1)]}$ is the thermal velocity of the cell. For v_th, we adopt the mass weighted thermal velocity of all star forming cells that lie within the measured radius. We take the median cell diameter and set the 25th percentile and 75th percentile cell diameters as upper and lower uncertainty regions.

In the right panel of Fig. 10, we compare the analytical growth rates making use of equations (9) and (10) with the fit growth rates obtained from the left panel. We find that the fit and calculated growth rates closely align for all halos. Generally, the fit growth rates are lower than the analytical ones, possibly due to losses as a result of the locking up of the magnetic field during star formation. For Halo3 and Halo260, where our fits approximated the growth very well, the fit and analytic growth rates are close to each other. For Halo4 and Halo194, where our fits did not perform as well, the fit growth rates are lower (by factors of around two) than the actual measured ones.

4.4. Magnetic field growth in the ICM

At lower redshifts (z < 3.5), the cluster has significantly grown, reducing the impact of individual galaxies on magnetic pre-enrichment. In this section, we provide further evidence for dynamo activity within the ICM.

4.4.1. Dynamo amplified magnetic correlation lengths

The magnetic field strength in the ICM remains constant during phase IV, as seen in Fig. 2. However, the accretion of less magnetized gas, as well as the decay of the magnetic field, would lower the field strength during this phase. The fact that it stays constant instead, shows that there are additional dynamo processes that continuously grow the field. The power spectra in Fig. 11 are calculated within a sphere of 0.5R₂₀₀ and illustrate the growth of magnetic energy over the last 10 Gyr in both radiative and non-radiative simulations. As predicted by dynamo theory, the peak of magnetic energy shifts to larger scales. The similarity in the evolution of the magnetic correlation length (equation (5)) in both simulations indicates that kinetic energy from cosmic accretion and merger events that are present in both cases, is the source. A major merger event at t = 10.7 Gyr (z = 0.3) moves the peak magnetic energy to smaller scales, where smaller-scale turbulence effectively amplifies the fields. Subsequently, the peak shifts back to larger scales while the total magnetic energy remains nearly constant. At z = 0, the largest dynamo amplified correlation length inside 0.5R₂₀₀ is ≈200 kpc. The kinetic power spectra follow a Kolmogorov scaling on larger scales (10⁻² kpc⁻¹  < k < 10⁻¹ kpc⁻¹) but steepen at smaller scales. At even smaller scales, we observed a flattening that could indicate an additional injection of small-scale kinetic energy, for example, as a result of feedback. This phenomenon can be studied in a future investigation. The magnetic power continues to grow together with the cluster, as (less magnetized) mergers bring new mass and kinetic energy, which acts to keep the magnetic field strength constant in the growing cluster volume.

Fig. 11. Kinetic (dashed) and magnetic (solid) power spectra of Halo3 presented within a radius of 0.5R200 at six different cosmic times, t, that are equally spaced (3.6 Gyr < t < 13.7 Gyr with approximately 2 Gyr separation in between; this corresponds to 1.5 > z > 0) for both the radiative (left) and the non-radiative (right) simulations. These times correspond to our phase IV (defined in Fig. 2), where the small-scale dynamo operates in the ICM. This fluctuating dynamo amplifies the field ejected from galaxies and grows it in strength and correlation scale via turbulence injected by cluster merger events. This process also occurs in the non-radiative simulations, albeit at a much reduced level because in this model there is no galactic pre-processing of magnetic fields. Over time, the peak magnetic energy shifts to larger scales.

4.4.2. Phase space diagrams

To further investigate the magnetic field growth through the kinematic and saturated phases of the dynamo, we analyze the phase space diagrams of Halo3 in the upper two rows of Fig. 12. These diagrams show mass-weighted histograms in the two-dimensional plane spanned by magnetic field strength and density for four distinct redshifts.

Fig. 12. Mass weighted histograms in the magnetic field versus gas density plane (upper two rows) in the ICM. Adiabatic compression of a flux-frozen magnetic field implies B ∝ ρ2/3, and a self-gravitating system with a saturated dynamo scales as B ∝ ρ1/2 (see black dashed lines). The radiative simulation starts to deviate early and first at high densities (corresponding to small scales) from the adiabatic relation, hinting toward dynamo activity at high redshifts. A rapid growth and early saturation is visible. By contrast, in the non-radiative simulation, magnetic amplification is significantly delayed and saturation is only reached at z = 0. The lower row shows the magnetic field vs. metallicity histogram in the ICM. The contours of constant density are indicated by brown-orange lines. These panels illustrate how the magnetic field is first amplified in high metallicity gas (in the center) and then, subsequently, in lower metallicity gas. The correlation between both quantities highlights the contributions of pre-enriched gas (in magnetic fields and metals) from merging galaxies in enriching the outskirts with both magnetic fields and metals.

Adiabatic compression in an isotropically collapsing halo and magnetic flux freezing yields the scaling B ∝ ρ^2/3. Observations of nearby clusters revealed a scaling of B ∝ ρ^0.47 (Bonafede et al. 2010; Govoni et al. 2017). Indeed, in a saturated small-scale dynamo, the specific magnetic energy ${\mathcal{E}}_B\propto v_{\text{A}}^2=B^2/(4\pi \rho )$ remains constant provided the turbulent energy does not change, which implies B ∝ ρ^1/2 (Xu & Lazarian 2020). Thus, we adopt this scaling for the saturated dynamo stage.

At z = 6.6, the high-density region (with the shortest dynamical time) in the radiative simulation first departs from the adiabatic relationship, hinting at dynamo activity or mixing of the ICM with magnetically enriched ISM. By z = 3.5, the cluster experiences major merger series, which causes the entire distribution in the radiative panel to deviate from the adiabatic expectation, especially at higher densities. In the non-radiative simulation, the denser region begins to gradually diverge from adiabatic behavior. At z = 1.5, the distribution in the radiative simulation enters the saturated regime, while the non-radiative simulation exhibits two power-law segments. These may represent the influence of two superposed dynamos, possibly triggered by the two major merger events. By z = 0, both distributions closely adhere to the scaling B ∝ ρ^1/2 characteristic of a saturated dynamo. In the non-radiative simulation, there are some regions that are not yet saturated, even though most of the ICM is in the saturated regime. Steinwandel et al. (2022) who simulate a non-radiative cosmological cluster show in their phase space diagrams (their Figure 6) that at z = 0 some parts of the cluster adhere to the B ∝ ρ^1/2 scaling, while other parts still follow the adiabatic relation, which is explained by their lower resolution. These observations summarize our findings: 1) magnetic fields grow initially via adiabatic compression, 2) there is a small-scale fluctuating dynamo that starts acting in the densest regions with the fastest eddy turnover rates and which is clearly much faster in the radiative simulation, and 3) eventually both simulations enter the saturated dynamo regime.

In the lower panel of Fig. 12, 2D histograms of magnetic field strength and metallicity are presented for the entire ICM (excluding star forming regions). Density contours are indicated with brown-orange lines. In our model, metals are transported from the ISM to the ICM in form of stellar feedback driven outflows. While in reality the outflows also carry magnetic fields, this does not happen in our subgrid model, therefore our results represent lower limits to the ICM magnetization at high redshifts. At z = 6.6, adiabatic compression of the proto-galaxy has increased the magnetic field to ∼10⁻⁴ μG while the magnetic field in the ICM starts to be amplified by a dynamo in the central regions that have been polluted with metals. By z = 3.5, the cluster has experienced many merger events with galaxies that contribute their high metallicity, strong magnetic field gas to the ICM yielding a positive correlation between both quantities. Additionally, the subsonic dynamo in the ICM starts to operate and magnetically amplifies the lower-metallicity gas, resulting in saturated magnetic fields by z = 0.

4.4.3. Saturation strengths

Another piece of evidence for a small-scale dynamo in the saturated regime is the constant ratio between magnetic and kinetic energy. There is an ongoing debate surrounding the existence of a universal saturation ratio between magnetic field energy density and kinetic energy density, both of which are expected to scale with cluster mass (see Sect. 4.1).

In Fig. 13, we present magnetic-to-kinetic energy density ratios for various radial bins (indicated by blue shadings) of Halo3. The kinetic energy is the sum of kinetic energies of each cell in a radial bin in the center of mass system. The ratio stays constant with time for all radial bins, indicating the presence of a dynamo. The ICM in the innermost radial bin (0 < r < 0.25R₂₀₀) saturates at 10% of kinetic energy, and the outermost radial bin (R₂₀₀ < r < 1.5R₂₀₀) at 1%. The reason for the decreasing magnetic-to-kinetic energy density ratio towards the outer radial bins can be twofold. 1) There is a progressively larger fraction of bulk kinetic energy in comparison to the turbulent energy density at larger radii. 2) The saturated field strength of a small-scale dynamo scales as B ∝ ρ^1/2 (see Fig. 12), implying lower magnetic field strengths at larger radii (see Fig. A.1). A merger causes the magnetic-to-kinetic energy ratio to drop at all radii at z ∼ 0.3 due to an injection of kinetic energy and gas with a lower magnetization into the host halo. We intend to improve the turbulent kinetic energy identification in future work.

Fig. 13. Ratios of the magnetic-to-kinetic energy density in the ICM calculated by averaging within spherical shells for the radiative simulation of Halo3. The progressively lighter colors indicate radial shells that move to larger radii. The magnetic energy saturates at around 10% in the innermost bin. The constant ratios at low redshifts indicate dynamo activity, while the lower values in the outer radial bins are likely due to a higher proportion of bulk kinetic in comparison to turbulent energy, as well as weaker magnetic fields in lower-density gas (see text).

4.5. A magnetic dynamo on magneto-hydrodynamics scales

The question arises whether the use of the MHD approximation is justified, as it assumes that the scales on which the small-scale dynamo operates in our simulation are collisional. This condition can be refined by requiring that the magnetic coherence length peak exceeds the particle mean free path during magnetic field amplification. Here, we demonstrate that in our proposed two-stage model of magnetic field evolution – high redshift growth in the ISM and low redshift growth in the ICM – amplification occurs entirely on collisional scales.

To start, we need to assume a generation process of magnetic fields. There are several processes that produce macroscopic fields on kpc scales or larger, including cosmological phase transitions in the early universe (Widrow 2002; Durrer & Neronov 2013) or the Biermann battery mechanism in ionization fronts or at shocks of strong, supernova-driven outflows (Attia et al. 2021). Turbulent motions then redistribute and amplify this seed field via dynamo activity.

In our proposed picture, field growth advanced through the following three stages:

1. In Section 4.3, we demonstrated that at high redshifts (between z = 8.5 and 5.5), the magnetic field is primarily amplified in the ISM, which is well approximated by the MHD approximation (see, e.g., Schekochihin et al. 2009). That is, the mean free path for typical ISM values is 1 AU in the cold neutral medium (Ryden & Pogge 2021) and ≈0.05 AU in the warm ionized medium (Schekochihin et al. 2009), which are much shorter than our numerical resolution. The dynamo exponentially amplifies the magnetic field on the scale of the fastest resolved eddies, which corresponds to our cell sizes. In a fluctuating small-scale dynamo, the scale of peak magnetic energy shifts to larger scales after the smaller scales have saturated.

2. The magnetic coherence length is further amplified by galactic winds and ram pressure stripping, which adiabatically expand the field and transport it into the ICM, as shown in Fig. 9.

3. By the time exponential amplification of the magnetic field develops in the ICM (at around z ≈ 5.5), the magnetic coherence length has already substantially grown to λ_B ≈ 5 kpc, as can be inferred from the power spectra in Fig. 6, so that the dynamo in the ICM operates on these scales. The ICM particle mean free path at these high redshifts is smaller in comparison to nearby clusters due to the higher densities in the early Universe. The particle mean free path in the ICM, assuming it to be fully ionized, can be calculated as (Braginskii 1965)

   $$\begin{array}{c}\hfill {\lambda }_{\mathrm{mfp}}=\frac{3\sqrt{2}}{4\sqrt{\pi }}\frac{k_{\mathrm{B}}^2{T}_{\mathrm{i}}^2}{ln{\mathrm{\Lambda }}_{\mathrm{i}}{n}_{\mathrm{i}}{e}^4}.\end{array}$$(11)

   Here, T_i, n_i, lnΛ_i are the ion temperature, number density, and Coulomb logarithm, respectively. The resulting particle mean free path in the ICM is shown in Fig. 14 for z < 5.5 (the blue curve shows the median and the 25^th and 75^th percentiles). Clearly, it is significantly smaller than the magnetic coherence length (shown in orange), which implies that the ICM dynamo remains in the fully collisional regime until z = 0. This justifies the use of the MHD approximation for modeling a dynamo in the ICM.

   Fig. 14. Comparison of the particle mean free path in the fully ionized ICM (equation (11)) and the magnetic coherence length (equation (5), which is close to the peak of the magnetic power spectrum from Fig. 6). This shows that by the time exponential amplification in the ICM begins (at z = 5.5), the magnetic coherence length has already significantly increased as a result of the dynamo activity in the ISM and the successive transport of the magnetic field into the ICM by means of galactic winds or ram pressure stripping, which further expanded the coherence length. This demonstrates that the dynamo in the ICM operates in the fully collisional regime throughout this process, rendering our MHD approximation valid.

This presents an alternative view to the plasma dynamo (St-Onge & Kunz 2018; Squire et al. 2019; St-Onge et al. 2020; Zhou et al. 2022, 2024), where the magnetic seed field in the ICM is generated by plasma instabilities. In this scenario, large-scale turbulence in a weakly collisional plasma induces pressure anisotropies, triggering the Weibel instability, which creates magnetic seed fields. Kinetic instabilities, such as the firehose and mirror instabilities, transport these fields to larger scales. These instabilities also provide effective viscosity at Larmor radius scales, where the small-scale dynamo can then further amplify the magnetic field. Future research is needed to demonstrate a fully self-consistent computation of the plasma dynamo on astrophysically relevant scales. Additionally, further work is required to identify observational consequences that may differ from the proposed scenario.

5. Faraday rotation maps

We simulate Faraday RMs at high redshifts. Our main finding is that magnetic fields are intrinsically linked to galaxy clusters as they form, originating in early galaxies and growing through subsequent mergers. Understanding these processes at high redshifts and in proto-BCGs is crucial and will be an active research topic of future observations.

Numerous radio observations of high redshift galaxies and protoclusters show that their magnetic field strength are comparable to that observed in nearby clusters (Di Gennaro et al. 2021; Anderson et al. 2022; Cordun et al. 2023; Chapman et al. 2024; Di Gennaro et al. 2023). However, precise field strength estimates from radio brightness observations are challenging due to uncertainties associated with the equipartition assumption in galaxy clusters, i.e., whether the ratio of cosmic ray proton-to-magnetic energy density is indeed unity and whether the radio-emitting electrons constitute 1% of the cosmic ray protons, as observed in the Milky Way (Ruszkowski & Pfrommer 2023). Faraday rotation describes the effect of rotation of the plane of polarization in a magnetized plasma due to its birefringent property. In such a plasma, the rotation angle is given by ϕ = ϕ₀ + RM λ², where ϕ and ϕ₀ are the measured and intrinsic polarization angle, respectively, and λ is the wavelength of radiation. The RM is given by (Xu & Han 2022)

$$\begin{array}{c}\hfill \mathrm{RM}=\frac{e^3}{2\pi m_{\mathrm{e}}^2{c}^4}{\displaystyle {\int }_0^L}\frac{n_{\mathrm{e}}\mathit{B}\mathbf{·}\mathrm{d}\mathit{l}}{{(z+1)}^2}.\end{array}$$(12)

Here, e represents the electron charge, m_e the electron mass, c the speed of light, n_e the number density of free electrons, B the magnetic field, and L the physical distance from the radio source to the observer. The (1 + z)⁻² factor arises because of the wavelength stretching in an expanding universe. These quantities are a direct output of our simulation, hence providing the possibility to create more accurate mock observations in comparison with, for example, radio maps that would require information about cosmic rays that are not simulated. Knowing the electron density along the line of sight (e.g., through thermal X-ray observations) enables us to infer the magnetic field strength, although this is however complicated by the non-monotonicity of the RM along the line of sight for a turbulent magnetic field and an inhomogeneous electron distribution (Hutschenreuter et al. 2024).

For our mock observations, we choose to mimic the Square Kilometre Array (SKA) due to its high angular resolution, sensitivity, and large sky coverage. We selected SKA-MID1 Band 2, which offers an angular resolution of ∼0.03 − 4″ (Swart et al. 2022) in a frequency range of 950 − 1740 MHz with a frequency resolution of 1 MHz and with a sensitivity of 4 μJ/beam (Heald et al. 2020). Although low-frequency observations are generally preferred for RM studies due to their higher frequency resolution and the need for narrow frequency sampling, the higher frequency range of SKA-MID1 Band 2 offers distinct advantages: it provides better angular resolution, and many sources that depolarize at low frequencies remain polarized at higher frequencies (Heald et al. 2020). Since we do not model cosmic ray electron injection and transport, which is required to realistically model synchrotron emission, we cannot predict the emerging polarized radio flux. Instead, we base our telescope specifications on the study by Heald et al. (2020), which proposes a polarization survey, including investigations of the early Universe and galaxy clusters. We generate RM maps using a filter that smooths the image by performing a convolution with a Gaussian kernel with a standard deviation equivalent to a 0.4″ resolution. We compute the angular size of the galaxy cluster using the angular distance between the cluster redshift and the observer at z = 0 using the cosmology that we used to create the simulation introduced in Section 2.3. For our adopted ISM model (Hernquist & Springel 2003), we calculate the free electron number density following Pakmor et al. (2018), where for star-forming regions, the fraction of the cold phase contributing to the density is subtracted. Additionally, we assume that the line-of-sight integral of RM is completely dominated by the virial region of the cluster. A contribution of |RM| of 10 rad m⁻² corresponding to limiting values expected from the |RM| contributions of the Milky Way (Pakmor et al. 2018) is not included (see also Anderson et al. 2022 who find a galactic foreground contribution of 13 rad m⁻²).

In Fig. 15, we present the results of our RM mock observations for the SKA at three different redshifts z = 6.2, 4.6, 3.5. The RM projections have sidelengths and projection depths of 2R₂₀₀. We indicate the virial radius (solid line) and the center as an inner circle with radius 20 kpc (dashed line). A resolution of 0.04″ enables the observation of the center of the simulated cluster, where the proto-BCG is active, at z = 6.2. This results in strong signals of ≈10⁴ rad m⁻² in the very center and in lower signals ≈10² rad m⁻² in a filamentary structure towards the top of the image. As redshift decreases, the signal strength grows while remaining nearly constant at the center. The increasing cluster size at lower redshifts allows for a more detailed study of the magnetic field’s morphology and the identification of merging galaxies, such as those near the virial radius in the top right of the image. Faraday RM observations conducted by LOFAR of the BCG of the spiderweb galaxy cluster at z = 2.2 demonstrate how the central AGN magnetizes the BCG’s circumgalactic medium (Anderson et al. 2022). Our mock observations agree with their actual observational RMs of ∼1500 rad m⁻² in the source frame.

Fig. 15. Square Kilometre Array mock Faraday RM observations of Halo3 for three different redshifts. The panels have sidelengths 40 kpc (left) and 2R200(z) (middle and right), and a projection depth of 2R200. The virial radius R200(z) is indicated with the solid circle. The central part of radius 20 kpc, indicative of the BCG, is shown with the dashed circle. Observations of the magnetized (proto-)BCG are possible as early as z = 6.2, where the cluster center shows high RM values, reflecting the interaction between AGN feedback and the magnetized ICM. From z = 4.6 onwards, the entire cluster center and the magnetic field morphology can be studied with a dense enough grid of background and embedded polarized sources. At z ≤ 3.5, the magnetic field can be studied out to larger radii using Faraday rotation.

However, the effectiveness of such experiments is constrained by the density of background sources available for RM measurements. For SKA, we expect to reach 60 sources/deg² (Jung et al. 2023), which translates to 0.003, 0.011, 0.027 sources per image for our z = 6.4, 4.6, 3.5 mock observations respectively. Nevertheless, radio emission can also arise from the cluster-intrinsic sources in form of radio lobes from AGNs (van Weeren et al. 2019). Furthermore, it has been suggested that radio halos exhibit a polarized fraction of 15–30% in the central regions of clusters, where the magnetic field is strong, particularly when observed at higher frequencies (Heald et al. 2020). Hence, instead of using a sparsely populated set of point sources to deduce the Faraday RM properties, we can use these diffuse and extended radio sources that continuously illuminate the cluster with polarized radio emission.

6. Discussion and conclusions

6.1. Comparison to observations and previous work

Verifying these findings through observations is challenging. Magnetic field strength estimates from radio brightness rely on the uncertain equipartition assumption (Ruszkowski & Pfrommer 2023), while estimates from Faraday RM data require a precise electron density along the line of sight. Additionally, fluctuating magnetic field directions can cause cancellation effects. However, this issue can be addressed with high-quality RM data, which enables the derivation of the magnetic power spectrum using second-order Faraday RM statistics (Vogt & Enßlin 2005; Kuchar & Enßlin 2011). Despite the large uncertainties, radio observations have already provided some support for our findings. Clusters at approximately z ≈ 1 show magnetic field strengths comparable to those in nearby clusters (Di Gennaro et al. 2021, 2023). Additionally, observations of high-redshift protoclusters with redshifts in the range 4.25 > z > 2.2 have revealed a magnetized ICM surrounding the BCG (Anderson et al. 2022; Cordun et al. 2023; Chapman et al. 2024). Notably, Faraday RM measurements of a protocluster at z = 2.2 yield RM values similar to those obtained from our mock observations of the radiative cluster at z = 3.5 (cf. Fig. 15 and Anderson et al. 2022). Moreover, we anticipate that the future capability of SKA will enable RMs for massive proto-BCGs as early as z = 6.4, thus providing an opportunity to test our predictions of early magnetic field growth.

However, the discrepancy between the magnetic Prandtl number in our simulations (P_m ≈ 1) and the real Prandtl numbers in the ICM (P_m ≈ 10²⁹; see, e.g., Tribble 1993; Kunz et al. 2022, Chapter 2.6.3) makes our simulation setup very different from reality. This discrepancy mainly influences the scale of magnetic field folding. Nevertheless, it has been shown that a P_m ≈ 1 dynamo behaves comparably to a large Prandtl number dynamo in growing the magnetic field, as the growth rate depends on the eddy turnover rates at the viscous scales (which is shown in, e.g., Schekochihin et al. 2004).

Our simulation methodology and findings are consistent with those of Marinacci et al. (2015), who compared magnetic field evolution of a radiative simulation to a non-radiative simulation in a cosmological box. We expanded on their work with our analysis by employing a higher resolution, an updated galaxy formation model, and an extended analysis. Consistent with our findings, they reported faster growth and saturation in the radiative simulation, although they did not analyze the physical reason for it. Importantly, their amplification and saturation occurred at lower redshifts, indicating lower growth rates, which can be attributed to their lower resolution. Both our findings and those of Marinacci et al. (2015) underscore the significance of considering galaxy formation physics when investigating magnetic fields in galaxy clusters.

Steinwandel et al. (2022) conducted a cosmological, non-radiative simulation of a cluster similar to our Halo3 but with half its mass. Despite using a resolution ten times better in mass than ours, they reported lower growth rates, with their magnetic field saturating at around z ≈ 2. In contrast, Steinwandel et al. (2024) simulated the same cluster at a resolution 20 times higher in mass than ours, finding that the magnetic field saturated at around z = 4, which is slightly earlier than in our radiative simulation results. This suggests that including radiative physics can accelerate magnetic field growth, compensating for lower resolution to a certain extent, but not indefinitely.

In an ideal scenario of a small-scale dynamo, the saturation is dependent on the magnetic Prandtl number, with Pr_m ≫ 1 yielding higher (43.8 − 1.3%) saturation values of the magnetic-to-turbulent energy density and Pr_m ≪ 1 yielding lower (2.43 − 0.135%) saturation values (Schober et al. 2015). Our simulations, however, exhibit saturation levels of ∼10% of the total kinetic energy in the center of mass frame in the central region. These values align with observations (Botteon et al. 2022), and they are also consistent with simulations of isotropic homogeneous turbulence conducted by Beresnyak (2012), where saturation values of approximately 4–5% were found. Similar agreement is observed in the simulations by Vazza et al. (2014), who explored the small-scale dynamo in galaxy clusters and reported an efficiency of kinetic turbulent energy conversion to magnetic energy of the order of 0.3, which is still within the same order of magnitude as our results, albeit slightly higher. Also Steinwandel et al. (2024) find saturation values of 0.1 in their halo, which is similar in mass and size to our halo. In contrast, simulations of disk galaxies demonstrate much higher saturation values, approaching nearly 10–100% (Pakmor et al. 2020; Whittingham et al. 2021; Pfrommer et al. 2022). As Pakmor et al. (2024) showed, these high saturation values are linked to disk galaxies, whereas dwarf galaxies with more spherical shapes exhibit lower saturation levels. They attribute this difference to the presence of a mean-field rotational dynamo in disks, which is not present in objects without ordered rotation.

Our results indicate that the BCG in the cluster center and merging galaxies play a pivotal role in enriching the ICM with magnetically pre-processed gas. However, the specific causes of gas loss from galaxies remain unclear. Our findings reveal that most galaxies entering the cluster arrive without self-bound gas. Moreover, we propose that AGN feedback may not significantly expel gas before z = 4. Stellar feedback appears to assume a more prominent role at higher redshifts in the initial galaxies, but whether it generates sufficiently robust stellar feedback driven winds extending beyond the galaxies’ virial radii remains uncertain. Ram pressure stripping is likely to be a prevalent phenomenon in both the host cluster and pre-merging substructures. Our findings are consistent with those of Kurinchi-Vendhan et al. (2024), who explored the causes of galaxy quenching in the TNG simulation. Their results suggest that AGN feedback and ram pressure stripping contribute to galaxy quenching, with most galaxies in cluster vicinities being quenched at around z ≈ 3. Due to various uncertainties in feedback mechanisms and high-redshift galaxy gas dynamics, which cannot be fully addressed within current galaxy formation models, we conclude that understanding the precise mechanisms responsible for galaxy quenching may not be necessary for concluding how the ICM is magnetized.

6.2. Conclusions

In this study, we have analyzed two types of cosmological galaxy cluster simulations. The first, the ‘radiative’ simulation, includes the galaxy formation model used in the IllustrisTNG simulations, while the ‘non-radiative’ simulation focuses only on gravity and MHD. We examined the evolution of the magnetic field strength in four clusters to account for sample variance and ensure robustness. In the remainder of our paper, we focused on the most massive cluster, Halo3, for detailed analysis, though the same qualitative results were found in the other clusters. Our analysis shows that the radiative simulation exhibits faster magnetic field amplification and saturation, as seen in Fig. 2. We also identified a connection between magnetic field evolution in the cluster and in galaxies. Our key findings are as follows:

• Cluster magnetic fields are intrinsically linked to galaxy clusters and their formation process. Initially, the magnetic field experiences amplification in the first galaxies that are formed in the cluster, which eventually merge to form the BCG. Subsequently, both the cluster mass and the magnetic field grow with the accretion of galaxies (Fig. 5).

• Magnetic fields grow fast within the forming galaxies at early times (z ≈ 8) due to a fluctuating small-scale dynamo that is excited by supersonic compressible (Burgers) turbulence induced by gravitational infall of cold gas (Fig. 7). We have the following evidence supporting the existence of such a dynamo: 1) The magnetic field of the ISM grows exponentially in time (Fig. 2); (2) a power spectrum analysis revealed a Kazantsev scaling in the kinematic phase of the dynamo on large scales, and after saturation of the small-scale field, we observed the growth of the magnetic coherence length (Fig. 6); and 3) the fit exponential growth rates agree with the theoretically expected rates of Burgers turbulence when we adopt our numerical Reynolds number (Fig. 10). This growth is already sufficient to match observational inferences and could be even faster for realistic Reynolds numbers in the ISM. Eventually, this leads to a saturated magnetic field in the ISM at around z ≈ 4.

• The magnetized gas in the BCG is transported to the surrounding ICM via turbulence that is injected by stellar feedback, the interaction with accreted gas, and (to a lesser extent) AGN feedback (Figs. 8 and 7).

• Later on, merging galaxies are responsible for enriching the ICM with the galactic magnetic field by means of galactic outflows driven by SN and AGN feedback as well as ram pressure stripping. Pure adiabatic expansion of the magnetic field in star formation regions to ICM densities could explain a large fraction of the ICM magnetic field strength (see Fig. 9). However, this overestimates the measured field strengths, as this simplified procedure neglects mixing losses and magnetic fields trapped in stars or retained in the ISM. This demonstrates that “pre-processing” (Fujita 2004), which is commonly associated with metal enrichment, likely also plays a significant role in the magnetic enrichment of the ICM in merging structures, which also manifests as a metallicity-magnetic field strength correlation in the ICM (Fig. 12).

• Cluster mergers inject incompressible (Kolmogorov) turbulence (Fig. 6) that triggers a fluctuating small-scale dynamo in the entire ICM at lower redshifts (z ≲ 5.5). This dynamo utilizes the pre-magnetized plasma expelled from galaxies as seed fields for further amplification, as we demonstrate by the following points: 1) A power spectrum analysis of the kinetic and magnetic energy density in the ICM revealed a Kazantsev scaling on large scales and showed an increasing magnetic coherence length after saturation on small scales (Fig. 11); (2) this dynamo is able to maintain an approximately constant volume-averaged field strength until z = 0 (Fig. 2) and counteracts the expected decrease of magnetic field strength over time as a result of mixing of accreted gas with smaller magnetization, giving rise to a constant magnetic-to-kinetic energy ratio after z ≈ 3.5 (Fig. 13); and 3) the correlation of B and gas density follows the relation of self-gravitating clouds and cluster observations (Fig. 12).

• Our results suggest that magnetic field strengths in clusters at redshift z = 3.5 should be comparable to those in nearby clusters (Fig. 2). Moreover, in the central galaxy, magnetic field strengths should be larger than 10 μG as early as z = 5.5 (Figures 2 and 5).

• We demonstrated that the magnetic field observed today in the weakly collisional ICM is always amplified on collisional scales: first, in the collisional ISM during the epoch of BCG assembly at z ≈ 8 and later on in the ICM on the magnetic coherence scale, which is always larger than the particle mean free path (Fig. 14). Hence, this justifies the use of MHD for studying the cluster dynamo, and it facilitates and by and large bypasses the extraordinary challenge to explain magnetic field growth in a weakly collisional plasma from the tiny seed values (St-Onge et al. 2020; Zhou et al. 2024).

A major outcome of our study is the realization that the magnetic field within the cluster experiences its initial amplification at high redshifts within the protocluster BCG. This amplification is attributed to compressive turbulence generated by gravitational infall of cold gas, which excites a small-scale dynamo. While previous research has theoretically explored how compressive turbulence can enhance magnetic fields in galaxies, this phenomenon has not yet been seen in cosmological simulations of galaxy clusters. Therefore, it is essential to rigorously test this outcome and determine its robustness for different galaxy formation models in the future.

Theoretical studies of the magnetic field evolution in galaxy clusters align with the forthcoming or ongoing telescope surveys. The James Webb Space Telescope can illuminate the formation and morphology of high-redshift protoclusters (Forrest et al. 2023; Muñoz et al. 2023). The recently launched X-Ray Imaging and Spectroscopy Mission is primed to study ICM turbulence in nearby clusters with exquisite detail (de Vries et al. 2023; Miller et al. 2023). Radio surveys such as the LOFAR Two-metre Sky Survey (Shimwell et al. 2017; Di Gennaro et al. 2021, 2023) and the upcoming SKA (Gitti et al. 2018) promise to explore the interaction between cosmic rays and magnetic fields within the ICM across a large range of redshifts. Hence, the motivation for theoretical studies employing simulations remains compelling.

Acknowledgments

We thank the referee for a careful reading of our manuscript and the valuable suggestions to improve it. We acknowledge support by the European Research Council under ERC-AdG grant PICOGAL-101019746 and by the German Science Foundation (DFG) through the Research Unit FOR-5195. TB gratefully acknowledges funding from the European Union’s Horizon Europe research and innovation programme under the Marie Skłodowska-Curie grant agreement No 101106080. JW acknowledges support by the German Science Foundation (DFG) under grant 444932369. RW acknowledges funding of a Leibniz Junior Research Group (project number J131/2022). The authors gratefully acknowledge the Gauss Centre for Supercomputing (GCS) for providing computing time on the GCS Supercomputer SuperMUC-NG at the Leibniz Supercomputing Centre (LRZ) in Garching, Germany, under project pn68cu. Software: We use Paicos (Berlok et al. 2024) to load and process the data. Paicos is a parallelized object-oriented Python package for analysis of simulations performed with AREPO https://github.com/tberlok/paicos. We also used astropy (Astropy Collaboration 2013, 2018, 2022), matplotlib (Hunter 2007), h5py (Collette 2013), cython (Behnel et al. 2011), scipy (Virtanen et al. 2020) and numpy (van der Walt et al. 2011).

References

Appendix A: Spherical profiles and Fourier transforms

A.1. Spherical profiles

The TNG model only produces cool core clusters at high redshift z ≳ 3 and does not allow for massive cool core clusters at low redshifts (Pakmor et al. 2023). This is because the kinetic energy injection in the AGN radio mode is quickly dissipated at the cluster centers, which leads to an increase of the central cluster entropy well above 30 keV cm², which is the cool core entropy threshold (Cavagnolo et al. 2009). Potentially, a more realistic AGN jet modeling (Ehlert et al. 2023; Weinberger et al. 2023) might leave the cool cores intact while still self-regulating star formation. By contrast, non-radiative simulations consistently produce non-cool cores due to the lack of cooling. The underlying gas density profiles also shape the resulting magnetic field strength profiles. To account for this behavior, we examine the co-evolution of magnetic field, density and metallicity. Figure A.1 shows spherically averaged radial profiles of Halo3 for density, magnetic field strength, entropy, and metallicity in radiative (upper panel) and non-radiative (lower panel) simulations for 6 different redshifts. In the second panel column of Fig. A.1, we present spherically averaged, volume-weighted root mean square magnetic field strength profiles. We fit single- or double-exponential profiles to our simulated magnetic profiles, depending on the requirement at a specific redshift:

$$\begin{array}{cc}\hfill B(r)& =B_0{\mathrm{e}}^{-k_{0}r},\hfill \end{array}$$(A.1)

$$\begin{array}{cc}\hfill B(r)& =B_1{\mathrm{e}}^{-k_{1}r}+B_2{\mathrm{e}}^{-k_{2}r},\hfill \end{array}$$(A.2)

where B_i (i ∈ {0, 1, 2}) are the magnetic field strengths at r = 0 and k_i⁻¹ are the characteristic scales of the corresponding exponential profiles and have units of kpc. We fit both of our simulations at z = 0 and z = 4.5 in the range 0.003R₂₀₀ < r < R₂₀₀ and report the fit values in Table A.1. Note that at z = 0, in the saturated dynamo regime, the magnetic field profile scales with the square root of the density profile (see Sect. 4.4.2), which could explain the slight deviations between the fit and measured curves. In the non-radiative run, the growth is smooth, while in the radiative simulation, magnetic fields increase exponentially, with exceptionally high strengths at the center at z = 3.5 (these were not seen previously when averaging over a larger volume).

The right panel of Fig. A.1 displays metallicity profiles exhibiting a similar trend: elevated values in the center at high redshift and a smoother distribution at lower redshifts. The transformation from cuspy to smooth occurs around z = 3.5. At z = 0, the distribution becomes even flatter, in comparison to the magnetic field profile.

Gas transport mechanisms are still a subject of investigation. In our model, stellar and AGN feedback contribute (especially in destroying the cool core), but mergers appear to be the primary drivers, injecting substantial energy that stirs the gas. Mergers introduce pre-enriched gas to the outskirts, flattening the magnetic and metallicity profiles (one merger can be seen as an extended bump in the magnetic field and metallicity profiles at z = 3.5 after which the profiles are flattened). The magnetic profiles experience a steeper increase in the outskirts in comparison to the metallicity profiles. This shows how the dynamo acts on top of the pre-enriched gas. The similarity between the density and magnetic profiles, i.e., the enhanced cores at high redshifts and the subsequent flattening, shows the contributions of adiabatic compression to the magnetic evolution. Also here, the magnetic field experiences a stronger increase in the outskirts, which is because of the small-scale dynamo that additionally amplifies the field.

Our study of pre-enrichment aligns with previous research, emphasizing the significant role of mergers in shaping these profiles. Notably, Vogelsberger et al. (2018), employing the same TNG galaxy formation model, report analogous findings. They find that by z = 0, 80% of all metal in the ICM are accreted.

Fig. A.1. We present spherically averaged radial profiles of density, magnetic field, entropy and metallicity for the radiative (upper panel) and the non-radiative simulation (lower panel) of Halo3. We show 6 different redshifts spanning a range z = 0 − 6.5 and a radius 10−3R200 < r < R200. We fit double or single exponential profiles to all magnetic field profiles shown (see Table A.1) and depict two of them for both simulations, respectively (for z = 0 and 4.5, using dashed lines). The radiative simulation shows a similar evolution in all profiles: at high redshifts, core values are very high. With decreasing redshift, mergers and feedback processes distribute the enhanced core values over larger radii. Both, metallicity and magnetic fields are results of pre-processing from galaxies. The profiles of the magnetic field strength experience a stronger increase at the outer radii in comparison to the metallicity profiles, which shows that a dynamo is operating on top of the pre-processed gas.

A.2. Fourier transform and power spectra

In order to compute the power spectrum of a (double) exponential profile, we need to first take the Fourier transform of a single exponential profile:

$$\begin{array}{cc}\hfill \mathcal{F}\left({\mathrm{e}}^{-k_{0}r}\right)(k)& ={\displaystyle {\int }_{-\infty }^{\infty }}{\mathrm{d}}^{3}x{\mathrm{e}}^{-k_0\sqrt{x^2+y^2+z^2}}{\mathrm{e}}^{-i\mathit{k}\mathit{x}}\hfill \\ \hfill & =4\pi {\displaystyle {\int }_0^{\infty }}\mathrm{d}r{r}^2{\mathrm{e}}^{-k_{0}r}\frac{\sin \left(kr\right)}{\mathit{kr}}=\frac{8\pi k_0}{{[k^2+k_0^2]}^2},\hfill \end{array}$$(A.3)

where $r=\sqrt{{\mathbf{x}}^2}$ and k = |k|. The 1D power spectrum P(k) is then obtained via (Pfrommer et al. 2022):

$$\begin{array}{cc}\hfill P(k)& =4\pi k^2{P}_{3\mathrm{D}}(k)\propto k^2{\left|\mathcal{F}\left(B_0{\mathrm{e}}^{-k_{0}r}\right)(k))\right|}^2=\frac{{\left(B_{0}8\pi k_{0}k\right)}^2}{{[k^2+k_0^2]}^4}\hfill \end{array}$$(A.4)

In the case of a double-exponential profile, the power spectrum is given by the absolute square of the Fourier transform of equation (A.2).

The numerically obtained power spectra (shown in Figs. 6 and 11) are computed by setting up a 3D regular grid onto which we interpolate the Voronoi data. We then apply a zero-padded fast Fourier transform of the three spatial components of $\sqrt{\rho }\mathbf{v}$ and $\mathbf{B}/\sqrt{4\pi }$, respectively. The power spectrum is obtained by taking the absolute square of the Fourier transformed components and then by averaging the power in N_bins = 1024 bins.

In the radiative simulations for 6.5 > z > 3.5, the radial gas density profile follows a double exponential form, which has an inner and outer characteristic scale. Figure A.2 shows that for z > 3.5, the inner characteristic scale (k₁⁻¹) of the double exponential profile is comparable to the magnetic coherence scale, indicating that the outer scale of magnetic turbulence could be impacted by that feature of the magnetic profile. A possible interpretation is that, in the inner regions of the cluster, the magnetic field grows efficiently through dynamo activity and saturates with a coherence length that is given by the average magnetic profile.

Fig. A.2. Comparison of the magnetic coherence scale (Equation 5) with the characteristic scales of a single and double exponential profile (see Table A.1) as a function of redshift. For z > 3.5, the inner characteristic scale (k1) of the double exponential profile in radiative simulations aligns with the magnetic coherence scale, suggesting a possible influence on the power spectra.

All Tables

All Figures

Fig. 1. Thin projections through the center of Halo3 displaying the root mean square magnetic field strength and vorticity, respectively, at z = 3 in the non-radiative (first row) and radiative (second row) simulations. The bottom two rows show the same but at z = 0. The virial radius R200 is indicated with a white circle. From left to right, the panels have sidelengths L = 4R200,  2R200,  0.5R200 and a thickness of L/10. At z = 3, the magnetic field in the radiative simulation is nearly saturated. The BCG appears as a bright spot with a strong magnetic field. An object infalling from the left has also developed a strong magnetic field that will contribute to the cluster’s final magnetic field. By contrast, at z = 3 the magnetic field in the non-radiative simulation remains weak. The infalling object is not visible, and there is no BCG. This underscores the magnetic field growth in galaxies and its transport to the ICM via merging halos.

In the text

Fig. 2. Evolution of magnetic field strength (derived from volume-weighted magnetic energy densities) in radiative (upper panels) and non-radiative (lower panels) simulations within R200. The left panels focus on Halo3; the right panels show the magnetic field strength across four halos. The magnetic field amplification is faster in radiative simulations. Its evolution in the entire ICM closely mirrors that of the ISM across all 4 halos To show this in more detail, we mark 4 phases of magnetic field evolution in the radiative simulation of Halo3: I) the collapse of the first galaxies, II) the amplification in the BCG, III) the amplification in the host cluster by mergers that contribute their magnetically pre-enriched gas to the host halo, IV) low redshift mergers that constantly inject turbulence, which drives a fluctuating small-scale dynamo.

In the text

Fig. 3. Resolution study of the evolution of the magnetic field strength averaged within R200 for the radiative (left) and non-radiative (right) simulations. We present the evolution for all four halos (shown in different colors) and compare two resolution levels: “zoom8” and “zoom12” (see text). The radiative simulation shows no significant differences between the two zoom levels for all four halos, reflecting the fact that the magnetic field is amplified predominantly in the ISM, which is well resolved in both resolutions. In the non-radiative simulation, significant differences emerge between the two zoom levels, with zoom8 simulations showing smaller growth rates and lower final field strengths for the two smaller clusters.

In the text

Fig. 4. Evolution of the gas-relaxation parameter 2Ekin, gas/|Epot|, which quantifies the dynamical state of the ICM for all four halos in the radiative (left) and non-radiative (right) simulations. In the radiative simulations, the high-redshift phase is dominated by frequent mergers and strong radiative gas cooling, which results in accretion and feedback events. In consequence, the ICM finds itself in an unrelaxed state in the forming protocluster. The non-radiative simulations display lower values of 2Ekin, gas/|Epot| at high redshift due to the lack of feedback and cooling. Mergers at late times can be seen as upwards fluctuations in this quantity in both simulation runs.

In the text

Fig. 5. Redshift evolution of the physical radius corresponding to 3R200 (≈9 Mpc) surrounding the cluster in the radiative simulation of Halo3. From left to right, we color code the spherically averaged fields of vorticity, metallicity, and magnetic field strength. The upper row shows the redshift range 20 > z > 0.1, highlighting the late time evolution dominated by mergers, while the bottom row shows the redshift range 11 > z > 4, focusing on the evolution in the central galaxy at higher redshifts. The dashed lines indicate the virial radius (R200) and the central cluster region (0.25R200). We keep the maximum radius in the lower panels at a fixed physical value of 150 kpc (see text). The lower panels show how vorticity and metallicity are enhanced in the center already at high redshifts, z > 9. The central magnetic field amplification sets in shortly afterwards at z ≈ 7. Turbulence from gravitational infall of cool gas amplifies the magnetic field and drives a small-scale dynamo. Stellar feedback additionally injects turbulence, distributes metals, and mixes magnetically enriched gas with less enriched gas. The upper panel highlights how merging substructures contribute their magnetically enriched gas to the host cluster. These merging structures also exhibit higher metallicity, indicating the connection between strong magnetic fields and the metal-rich ISM.

In the text

Fig. 6. Kinetic (dashed) and magnetic (solid) power spectra of Halo3. The left panels show the power spectra in the central regions within a fixed physical radius of R200(z = 9) = 20 kpc for four different redshifts (z = 4.5, 5.5, 7, 8.5) for both the radiative (upper) and the non-radiative (lower) simulations aimed at tracing the high-redshift dynamo in the BCG. The right panels show the ICM power spectra within 0.5R200(z), excluding the central 20 kpc, at six different redshifts (z = 0, 1.5, 2.5, 3.5, 4.5, 6.5) focusing on dynamo activity in the ICM within the radiative and non-radiative simulations. Dotted lines illustrate different scaling behaviors (see text). The left panels show that a small-scale dynamo (with the Kazantsev scaling ∝k3/2) driven by compressive Burgers turbulence (∝k−2) is only active in the central regions of the radiative simulation. By contrast, the right panels show the existence of a small-scale dynamo driven by incompressible turbulence (∝k−5/3) throughout the entire ICM within 0.5R200 in both simulations, albeit at a faster rate in the radiative simulation. The contribution of a Lorentzian profile, which we show for z = 4.5 in both simulations and additionally z = 0 in the radiative simulation (it scales as ∝k2 on scales larger than the scale height), at lower redshifts is minimal because the sharp cutoff at the characteristic exponential scale is located at rather large scales.

In the text

Fig. 7. Thin projections through the center of Halo3 in the radiative simulation with sidelength 2.2R200 at z = 6.2 and thickness 0.5R200. From top left to bottom right, the nine panels display the gas velocity in the cluster rest frame, sound speed, Mach number, star formation rate surface density, vorticity, the ratio of the free-fall time to the eddy-turnover time, particle number density, magnetic field strength, and metallicity. The velocity field is shown together with its absolute value in the first panel. Regions with cold, supersonic gaseous inflows create vorticity and display strong magnetic fields. Star formation in the high density regions creates stellar feedback driven winds that interact with the inflowing gas and produce vorticity. This vorticity efficiently mixes highly magnetized gas with weaker magnetized gas, distributing metals across the cluster.

In the text

Fig. 8. Energy injection rates into gas cells of Halo3 in the radiative simulation are shown. The panels differentiate between feedback launched from the BCG center (r < 20 kpc; left panel) and from the rest of the cluster, excluding the center (right panel). We present different energy injection rates that originate from stellar feedback (yellow) as well as thermal (purple) and kinetic (orange) AGN feedback. The injected thermal AGN feedback (dashed purple) is numerically reduced in star-forming cells where an effective equation of state is enforced, implying a rapid energy loss and a smaller effective thermal AGN feedback energy rate (solid purple). Radiative cooling losses (blue) are displayed for the ICM in the right panel. At high redshifts, both in the BCG and the entire cluster, the primary energy injection rate is due to stellar feedback. Subsequently, AGN feedback also becomes a significant contributor.

In the text

Fig. 9. Left panel: Magnetic energy density in the ISM (red). We measured the energy density in an expanding sphere with a radius equal to the virial radius at the time. When the magnetic energy density inside this sphere grew for the first time by a factor of 102 compared to the adiabatically compressed field right after the protocluster collapse, we fixed the radius at this redshift. Right panel: Magnetic energy density in the ICM (blue) compared to the magnetic energy density expected from expanding the ISM magnetic field purely adiabatically to the volume of the ICM (red dotted). This estimate even overshoots the total magnetic energy in the ICM. The lower measured ICM magnetic field strengths, compared to adiabatically expanded ISM fields, may result from depletion during star formation, retention in the ISM, or unaccounted mixing losses.

In the text

Fig. 10. Left: Time evolution of the magnetic energy density with fits of the exponential growth rates for the ISM in all four halos (dotted lines). The fits started when the magnetic energy density for the first time grew to a value at or above 4 × 10−7 erg cm−3 and ended when the magnetic energy density for the first time reached a value of at least 40 erg cm−3. We subtracted the contributions from adiabatic compression following the relation given in equation (7), where we fixed B0 and ρ0 to be the measured values after the cluster collapse, and εB(t), ρ(t) were measured in each snapshot. Right: Fit growth rates (circles) along with the expected growth rates (crosses) that include error bars as described in the text. Our analytic and measured growth rates agree within a factor of two, with the expected growth rates slightly larger than the measurements.

In the text

Fig. 11. Kinetic (dashed) and magnetic (solid) power spectra of Halo3 presented within a radius of 0.5R200 at six different cosmic times, t, that are equally spaced (3.6 Gyr < t < 13.7 Gyr with approximately 2 Gyr separation in between; this corresponds to 1.5 > z > 0) for both the radiative (left) and the non-radiative (right) simulations. These times correspond to our phase IV (defined in Fig. 2), where the small-scale dynamo operates in the ICM. This fluctuating dynamo amplifies the field ejected from galaxies and grows it in strength and correlation scale via turbulence injected by cluster merger events. This process also occurs in the non-radiative simulations, albeit at a much reduced level because in this model there is no galactic pre-processing of magnetic fields. Over time, the peak magnetic energy shifts to larger scales.

In the text

Fig. 12. Mass weighted histograms in the magnetic field versus gas density plane (upper two rows) in the ICM. Adiabatic compression of a flux-frozen magnetic field implies B ∝ ρ2/3, and a self-gravitating system with a saturated dynamo scales as B ∝ ρ1/2 (see black dashed lines). The radiative simulation starts to deviate early and first at high densities (corresponding to small scales) from the adiabatic relation, hinting toward dynamo activity at high redshifts. A rapid growth and early saturation is visible. By contrast, in the non-radiative simulation, magnetic amplification is significantly delayed and saturation is only reached at z = 0. The lower row shows the magnetic field vs. metallicity histogram in the ICM. The contours of constant density are indicated by brown-orange lines. These panels illustrate how the magnetic field is first amplified in high metallicity gas (in the center) and then, subsequently, in lower metallicity gas. The correlation between both quantities highlights the contributions of pre-enriched gas (in magnetic fields and metals) from merging galaxies in enriching the outskirts with both magnetic fields and metals.

In the text

Fig. 13. Ratios of the magnetic-to-kinetic energy density in the ICM calculated by averaging within spherical shells for the radiative simulation of Halo3. The progressively lighter colors indicate radial shells that move to larger radii. The magnetic energy saturates at around 10% in the innermost bin. The constant ratios at low redshifts indicate dynamo activity, while the lower values in the outer radial bins are likely due to a higher proportion of bulk kinetic in comparison to turbulent energy, as well as weaker magnetic fields in lower-density gas (see text).

In the text

Fig. 14. Comparison of the particle mean free path in the fully ionized ICM (equation (11)) and the magnetic coherence length (equation (5), which is close to the peak of the magnetic power spectrum from Fig. 6). This shows that by the time exponential amplification in the ICM begins (at z = 5.5), the magnetic coherence length has already significantly increased as a result of the dynamo activity in the ISM and the successive transport of the magnetic field into the ICM by means of galactic winds or ram pressure stripping, which further expanded the coherence length. This demonstrates that the dynamo in the ICM operates in the fully collisional regime throughout this process, rendering our MHD approximation valid.

In the text

Fig. 15. Square Kilometre Array mock Faraday RM observations of Halo3 for three different redshifts. The panels have sidelengths 40 kpc (left) and 2R200(z) (middle and right), and a projection depth of 2R200. The virial radius R200(z) is indicated with the solid circle. The central part of radius 20 kpc, indicative of the BCG, is shown with the dashed circle. Observations of the magnetized (proto-)BCG are possible as early as z = 6.2, where the cluster center shows high RM values, reflecting the interaction between AGN feedback and the magnetized ICM. From z = 4.6 onwards, the entire cluster center and the magnetic field morphology can be studied with a dense enough grid of background and embedded polarized sources. At z ≤ 3.5, the magnetic field can be studied out to larger radii using Faraday rotation.

In the text

Fig. A.1. We present spherically averaged radial profiles of density, magnetic field, entropy and metallicity for the radiative (upper panel) and the non-radiative simulation (lower panel) of Halo3. We show 6 different redshifts spanning a range z = 0 − 6.5 and a radius 10−3R200 < r < R200. We fit double or single exponential profiles to all magnetic field profiles shown (see Table A.1) and depict two of them for both simulations, respectively (for z = 0 and 4.5, using dashed lines). The radiative simulation shows a similar evolution in all profiles: at high redshifts, core values are very high. With decreasing redshift, mergers and feedback processes distribute the enhanced core values over larger radii. Both, metallicity and magnetic fields are results of pre-processing from galaxies. The profiles of the magnetic field strength experience a stronger increase at the outer radii in comparison to the metallicity profiles, which shows that a dynamo is operating on top of the pre-processed gas.

In the text

	Fig. A.2. Comparison of the magnetic coherence scale (Equation 5) with the characteristic scales of a single and double exponential profile (see Table A.1) as a function of redshift. For z > 3.5, the inner characteristic scale (k1) of the double exponential profile in radiative simulations aligns with the magnetic coherence scale, suggesting a possible influence on the power spectra.
In the text