© The Authors 2024

1. Introduction

Magnetic fields are observed in a wide range of systems throughout the Universe. They are present in stars (e.g., Hubrig et al. 2011), including the Sun (e.g., Scherrer et al. 1977; Schmassmann et al. 2018) and the interstellar medium of various types of galaxies in the local Universe (e.g., Beck & Wielebinski 2013) and at intermediate redshifts (e.g., Widrow 2002; Bernet et al. 2008). Radio observations even reveal the existence of magnetic fields of a few μG in the intracluster medium (ICM; e.g., Vogt & Enßlin 2003; Bonafede et al. 2010). This corresponds to a magnetic energy density that is close to the energy density of turbulent motions. Incidentally, magnetic fields on the order of a few μG have been observed in diffuse radio emission in galaxy clusters up to redshift z ≃ 0.7 (Di Gennaro et al. 2021). The origin of such strong magnetic fields, in the largest gravitationally bound objects of the Universe, galaxy clusters, remains unclear.

The seeds for the ICM magnetic field could have already been generated in the very early Universe, as predicted from certain models of inflation (e.g., Fujita et al. 2015; Adshead et al. 2016; Talebian et al. 2020) and cosmological phase transitions (Quashnock et al. 1989; Törnkvist 1998; Ellis et al. 2019). Cosmological seed magnetic fields could have been further amplified by a macroscopic quantum effect called the chiral dynamo (e.g., Joyce & Shaposhnikov 1997; Brandenburg et al. 2017; Schober et al. 2022). Due to baryon overdensity, the Biermann battery term (Biermann 1950) provides another mechanism for seed field generation before recombination (e.g., Naoz & Narayan 2013). However, observations of the cosmic microwave background provide strict upper limits for the comoving field strength of cosmological seed fields between 10⁻⁹ (Planck Collaboration XIX 2016) and 10⁻¹¹ G (Jedamzik & Saveliev 2019) on length scales above approximately 1 Mpc. Alternatively, seed fields could be generated through astrophysical mechanisms such as the Biermann battery during structure formation (e.g., Kulsrud et al. 1997a; Ryu et al. 1998) or reionization (e.g., Subramanian et al. 1994; Gnedin et al. 2000). Yet, astrophysical magnetogenesis produces magnetic fields that are too weak to explain the observed values in the ICM. Hence, an efficient amplification process is needed.

One commonly discussed candidate for magnetic field amplification is the small-scale turbulent dynamo (Schekochihin & Cowley 2006; Miniati & Beresnyak 2015; Domínguez-Fernàndez et al. 2019). This is a magnetohydrodynamical (MHD) phenomenon, in which turbulence stretches, twists, and folds the magnetic field on small spatial scales. In the kinematic stage, the magnetic energy spectrum exhibits a Kazantsev-like power law (Kazantsev 1968) with a peak at the resistive scale k_η (Schekochihin et al. 2004; Brandenburg et al. 2023). For typical plasma parameters of the ICM, this corresponds to k_η ≈ (10⁴ km)⁻¹ (Schekochihin & Cowley 2006). However, observations of the Faraday rotation measure in the ICM reveal a typical reversal scale of the magnetic field of a few kiloparsecs (e.g., Bonafede et al. 2010 for the Coma cluster). This could indicate that the small-scale dynamo is in the nonlinear evolutionary stage, where magnetic energy shifts to larger spatial scales (Schekochihin et al. 2002; Beresnyak 2012; Schleicher et al. 2013). Whether or not this scenario can explain the existence of μG magnetic fields in galaxy clusters depends on the properties of ICM turbulence across cosmic times.

In hierarchical structure formation, galaxy clusters have gradually formed through a series of dark matter halo mergers and gas accretion (e.g., Mo et al. 2010). Such mergers are some of the most energetic events in the Universe and can drive turbulence. The hierarchical assembly of dark matter halos leading to the formation of galaxy clusters has been extensively studied, both analytically through the extended Press-Schechter theory (Press & Schechter 1974; Cole et al. 2000) and via cosmological simulations (for example the Millennium Simulation; Springel et al. 2005) and IllustrisTNG (e.g., Marinacci et al. 2018). Such tools can be employed to study the evolution of turbulence and therefore magnetic field amplification in the ICM. In fact, the growth rate of the classical small-scale dynamo in its kinematic stage increases with increasing Reynolds number Re (Kazantsev 1968), which measures the ratio between inertial and viscous forces. While in galaxies Re reaches values of approximately 10¹⁰ (e.g., Schober et al. 2013) and therefore allows us to explain the generation of the strong random component of galactic magnetic fields by a small-scale dynamo, it is estimated to be only on the order of a few tens in the ICM of local clusters (assuming the standard Spitzer viscosity; e.g. Schekochihin & Cowley 2006 and references therein). With such a low value of Re, amplifying the seed magnetic fields to equipartition with turbulent gas motions through a classical small-scale dynamo might not be possible.

Another complication of the classic small-scale dynamo scenario arises from the fact that MHD theory might not be suitable for describing the ICM. The applicability of MHD depends on the collisionality of the plasma. The ion-ion collision frequency ν_ii is inversely proportional to the plasma temperature T as ν_ii ∝ T^−3/2 (Fitzpatrick 2014). With typical virial temperatures of 10⁷ to 10⁸ K generated by mergers of dark matter halos, ${\nu }_{\mathrm{ii}}^{-1}/t_{\mathrm{clust}}\approx {10}^{-2}$ (e.g., Kravtsov & Borgani 2012), where t_clust is the typical dynamical time of a cluster. This suggests that the ICM plasma is weakly collisional. In weakly collisional plasmas (WCPs), the first magnetic moment of a particle, $\mu \equiv m{v}_{\perp }^2/(2B)$, is an adiabatic invariant, meaning that whenever a change in the magnetic field B is produced, the perpendicular velocity v_⊥ of a particle with mass m adjusts in order to conserve μ. This produces a bias in the thermal motions of the particles, leading to pressure anisotropies along and perpendicular to the local direction of the magnetic field. The main issue is that, in a μ-conserving environment, dynamo action is impossible. Indeed, when the magnetic field grows, a corresponding change in the anisotropic temperature should occur, but there is not enough thermal energy in the plasma to be redistributed that way (Helander et al. 2016). Besides, pressure anisotropies trigger microscale, fast-growing kinetic instabilities that break the μ-invariance, therefore allowing the dynamo process to kick in. The most commonly encountered instabilities in the literature are the mirror and the firehose ones¹ (e.g., Chandrasekhar et al. 1958; Parker 1958; Rudakov & Sagdeev 1961; Gary 1992; Southwood & Kivelson 1993; Hellinger 2007; Kunz et al. 2014; Melville et al. 2016). Therefore, in WCPs, discarding kinetic effects from macroscale dynamics is no longer possible, which constitutes an enormous theoretical and numerical challenge.

In the past years, an extensive amount of simulations and theoretical modeling of WCPs have been conducted for different magnetization states. For instance, Santos-Lima et al. (2014) performed simulations of the turbulent dynamo based on the double-adiabatic Chew-Goldberger-Low equations (Chew et al. 1956), in which they implemented an anomalous collisionality aiming to mimic the particle-scattering produced by firehose and mirror instabilities. Furthermore, Rincon et al. (2016) used low-resolution hybrid-kinetic simulations of a turbulent dynamo in a collisionless, unmagnetized plasma to demonstrate that a dynamo effect is effectively possible in such a system. Hybrid particle-in-cell simulations of the turbulent dynamo in the “kinematic” and nonlinear phases, as well as in the saturated state of a collisionless magnetized plasma, have been performed by St-Onge & Kunz (2018), where they calculated the effective collisionality produced by the microscale kinetic instabilities. Finally, St-Onge et al. (2020) performed extensive simulations of a WCP based on MHD equations supplemented by a field-parallel Braginskii viscous stress. A significant result from these two studies is that the parallel viscosity “knows” about the direction and strength of the magnetic field, which basically leads to a dependency of the Reynolds number on the magnetic field, meaning that the Reynolds number Re is replaced by an effective version Re_eff(B). In the magnetized “fluid” regime (corresponding to magnetic field strengths above a few nG), St-Onge et al. (2020) argued that the effective Reynolds number scales as Re_eff ≈ β²ℳ⁴, where β (thermal pressure to magnetic pressure ratio) and ℳ are the plasma parameter and the Mach number, respectively. On the other hand, Rincon et al. (2016) found that in the unmagnetized (collisionless) regime, the magnetic field evolves into a folded structure with a magnetic power spectrum close to the high magnetic Prandtl number MHD case. However, the form of this relation in the magnetized “kinetic” regime is still an open question. Despite advancements in simulations and modeling, investigating the dynamics of magnetic fields over cosmological timescales on which galaxy clusters form remains a challenging task. This is due to the vast range of values of ρ_i/L₀ (with ρ_i being the ion Larmor radius and L₀ the fluid-scale stretching) that must be covered, going from seed-field values (which can be as low as 10⁻²⁰ G) to values close to equipartition in the ICM (approximately 10⁻⁶ G). This range covers both the kinetic and fluid scales, making it computationally infeasible at present.

In this work, we circumvented the difficulties related to the resolution of cosmological simulations by using a semi-analytical approach based on merger trees. These algorithms function by tracing the merger history of dark matter halos and are used either as a compact representation of the merging and accretion rates of dark matter halos in simulations (e.g., Tweed et al. 2009) or as statistical tools to study the formation and evolution of cosmic large-scale structures (e.g., Dvorkin & Rephaeli 2011; Gómez et al. 2022). In particular, Cole et al. (2000) developed the GALFORM model, which is based on the extended Press-Schechter theory (Snyder et al. 1997) to study the evolution of galaxy formation. This model was later extended by Parkinson et al. (2008), who modified the conditional mass function from the latter in order to match the statistical predictions from the Millennium simulations (Springel et al. 2005). Our objective is to generate merger trees based on the modified GALFORM algorithm and to study the evolution of magnetic fields in galaxy clusters by implementing a semi-analytical model of the dynamo based on the results from state-of-the-art simulations of WCPs.

The paper is organized as follows. In Sect. 2, we present the global theoretical framework of the different magnetization regimes in WCPs, along with the physics of the pressure-anisotropy-driven microscale instabilities. In Sect. 3, we present the modified GALFORM algorithm and the construction of our semi-analytical model of the turbulent dynamo. Our results are discussed in Sect. 4, and our conclusions are presented in Sect. 6.

2. Theory of weakly collisional plasmas

To understand how initial seed fields are amplified by dynamo processes in WCPs, we must divide the dynamo problem into different magnetization regimes, and we deal with them separately.

2.1. The unmagnetized regime (UMR)

The UMR corresponds to the situation where the ion cyclotron frequency Ω_i is much smaller than any other dynamical frequency ω, and the Larmor radius ρ_i is much larger than any other spatial scale l_clust of the system. In other words, we have Ω_i/ω ≪ 1 and ρ_i/l_clust ≫ 1. This corresponds to magnetic field amplitudes of

$$\begin{array}{c}\hfill B\lesssim {10}^{-18}\left(\frac{n_{\mathrm{e}}}{{10}^{-3}{\mathrm{cm}}^{-3}}\right){\left(\frac{T}{5\mathrm{keV}}\right)}^{-3/2}\mathrm{G},\end{array}$$(1)

where n_e is the electron density. Equation (1) indicates that magnetic fields on the order of 10⁻¹⁸ G are strong enough to drive the ICM into a magnetized regime (St-Onge et al. 2020).

Fully kinetic numerical simulations of the Vlasov equation were conducted by Rincon et al. (2016) and show that in an initially unmagnetized plasma, magnetic field amplification does occur in stochastically driven, nonrelativistic flows. Moreover, they have shown that the magnetic field evolves with a power spectrum that resembles its high-Prandtl-number MHD counterpart. This is a solid justification that turbulent folding of the magnetic field is a main process that amplifies magnetic energy during merger events during the formation of a galaxy cluster.

2.2. The magnetized regime (MR)

As soon as ρ_i ≈ l_clust and Ω_i ≈ ω, the system enters the magnetized regime, and the new scaling relations will drastically change the magnetic field amplification problem. Magnetization along with the condition ρ_i, α/λ_mfp, α ≪ 1, where λ_mfp, α stands for the mean free path of particle species α, leads to a pressure tensor $\overline{\mathit{p}}$ with anisotropic components with respect to the local direction of the magnetic field $\widehat{\mathit{b}}\equiv \mathit{B}/B$, which can be expressed as

$$\begin{array}{c}\hfill \overline{\mathit{p}}=p_{\perp }\overline{\mathit{I}}+(p_{\Vert }-p_{\perp })\widehat{\mathit{b}}\widehat{\mathit{b}},\end{array}$$(2)

where $\overline{\mathit{I}}$ is the unity tensor, and p_∥ and p_⊥ are, respectively, the pressure along and perpendicular to the local direction of the magnetic field. The pressure anisotropy term p_∥ − p_⊥ plays a major role in the stability of a WCP, as explained below. By taking the first velocity moment of the Vlasov equation in the guiding-center limit (Kulsrud 1983), we obtain the following modified Euler momentum equation:

$$\begin{array}{c}\hfill mn\frac{\mathrm{d}\mathit{v}}{\mathrm{d}t}=-\mathrm{\nabla }(p_{\perp }+\frac{B^2}{8\pi })+\mathrm{\nabla }·[\widehat{\mathit{b}}\widehat{\mathit{b}}(p_{\perp }-p_{\Vert })+\frac{B^2}{4\pi }],\end{array}$$(3)

where m, n, and v, respectively, are the proton mass, density, and velocity, and where d/dt ≡ ∂_t + (v ⋅ ∇) is the Lagrangian derivative along the trajectory of a fluid element. Equation (3) requires an additional closure relation constraining the evolution of the anisotropic pressure. Considering the guiding-center limit of the Vlasov equation (Kulsrud 1983), we obtain the double-adiabatic equation of state, also known as the Chew-Goldberger-Low (CGL) equations (Chew et al. 1956). This is written as

$$\begin{array}{c}\hfill \frac{d}{\mathrm{d}t}\left(\frac{p_{\perp }}{\mathit{nB}}\right)=\frac{d}{\mathrm{d}t}\left(\frac{p_{\Vert }{B}^2}{n^3}\right)=0.\end{array}$$(4)

In the limit of the CGL equations, any local change of the magnetic field creates a corresponding local pressure anisotropy (Rincon 2019):

$$\begin{array}{c}\hfill \mathrm{\Delta }\equiv \frac{p_{\perp }-p_{\Vert }}{p},\end{array}$$(5)

where p ≡ (p_∥ + 2p_⊥)/3. However, μ-conservation constitutes in itself a fundamental problem regarding dynamo effects. Indeed, it has been pointed out by Kulsrud et al. (1997b) that any increase in the magnetic field amplitude has to cause the same order-of-magnitude increase in the perpendicular thermal energy of the plasma, which is energetically impossible in a high-β plasma stirred by subsonic turbulent motions. In other words, magnetic-field amplification is impossible in a μ-conserved, magnetized plasma (Helander et al. 2016; Santos-Lima et al. 2014). For a dynamo action to kick in, μ-conservation has to be broken. In the Braginkskii limit (Braginskii 1965), where the collision frequency is large compared to the inverse of the timescale on which the magnetic field strength changes, we obtain the following approximation:

$$\begin{array}{c}\hfill \mathrm{\Delta }\simeq \frac{1}{{\nu }_{\mathrm{ii}}}\frac{1}{B}\frac{\mathrm{d}B}{\mathrm{d}t}=\frac{1}{{\nu }_{\mathrm{ii}}}\widehat{\mathit{b}}\widehat{\mathit{b}}:\mathrm{\nabla }\mathit{u}.\end{array}$$(6)

We note that the second equality in the equation above comes from the classical ideal induction equation and the assumption of incompressibility of the flow. Equation (6) implies that any change in the magnetic field will automatically create corresponding changes in the pressure anisotropies. Whenever Δ grows, the system is driven away from stability, and plasma instabilities are triggered, the well-known of which are the mirror and firehose instabilities (Chandrasekhar et al. 1958; Parker 1958; Rudakov & Sagdeev 1961; Gary 1992; Southwood & Kivelson 1993; Hellinger 2007; Melville et al. 2016). The interested reader is invited to consult the extensive aforementioned literature for more details about those instabilities, as we only list basic concepts here. The mirror instability is triggered when Δ > 1/β, and the firehose one is triggered when Δ ≲ −2/β. All those instabilities develop on scales close to the ion Larmor radius at a frequency close to the gyrofrequency (e.g., Melville et al. 2016). In other words, they are much faster compared to the fluid-scale dynamics of a WCP. The effect of those instabilities on the global dynamics of the magnetic field is still unclear, and, hence, the mechanism that allows those kinetic-scale instabilities to pin the system at the stability threshold during the nonlinear stage remains elusive. Two scenarios are found in the literature. There is a possibility that mirror and firehose fluctuations screen particles from variations of the magnetic-field amplitude. For instance, in the firehose case, the stretching rate of the magnetic field by turbulent eddies is reduced, which could also be seen as an enhancement of the plasma viscosity (Mogavero & Schekochihin 2014). Another possible effect of the kinetic-scale instabilities is to scatter particles (Kunz et al. 2014). This particle-scattering process acts as an effective reduction of the plasma viscosity. Overall, both processes are observed in simulations (Kunz et al. 2014; Riquelme et al. 2015). The transition between those two processes is still not clear, although the screening mechanism seems to be predominant over the particle-scattering one for small magnetic field variations δB/B ≪ 1 (Rincon 2019).

Given that effects from pressure-anisotropy-induced kinetic instabilities are much faster compared to the dynamical timescales of the system considered, a common numerical method to prevent the system from crossing the stability thresholds employs the so-called hard-wall limiters (Sharma et al. 2006), which are

$$\begin{array}{c}\hfill p_{\perp }-p_{\Vert }=\{\begin{array}{cc}max(-\frac{B^2}{4\pi },3{\mu }_{\mathrm{B}}\widehat{\mathit{b}}\widehat{\mathit{b}}:\mathrm{\nabla }\mathit{u})\hfill & \mathrm{firehose}\mathrm{limit}\hfill \\ min(\frac{B^2}{8\pi },3{\mu }_{\mathrm{B}}\widehat{\mathit{b}}\widehat{\mathit{b}}:\mathrm{\nabla }\mathit{u})\hfill & \mathrm{mirror}\mathrm{limit}\hfill \end{array},\end{array}$$(7)

with μ_B being the Braginskii viscosity (Braginskii 1965). Pinpointing the system at the instability threshold using Eq. (7) implies a modified collisionality of the form ${\nu }_{\mathrm{ii},\mathrm{eff}}\propto \beta \widehat{\mathit{b}}\widehat{\mathit{b}}:\mathrm{\nabla }\mathit{u}$, which leads to an effective Reynolds number that depends on the magnetic field through the β-parameter (St-Onge et al. 2020) as

$$\begin{array}{c}\hfill {\mathrm{Re}}_{\mathrm{eff}}\approx {\beta }^2{\left(\frac{v_{\mathrm{turb}}}{v_{\mathrm{th},i}}\right)}^4.\end{array}$$(8)

However, such limiters can only be applied when ν_ii, eff ≲ Ω_i, which corresponds approximately to magnetic field amplitudes above

$$\begin{array}{c}\hfill B\gtrsim 6{\left(\frac{n_{\mathrm{e}}}{{10}^{-3}{\mathrm{cm}}^{-3}}\right)}^{2/5}{\left(\frac{T}{5\mathrm{keV}}\right)}^{1/2}{\left(\frac{\mathcal{M}}{0.2}\right)}^{3/5}{\left(\frac{L_0}{100\mathrm{kpc}}\right)}^{-1/5}\mathrm{nG}.\end{array}$$(9)

To sum up, in order to model the magnetic field evolution from seed-field values to equipartition-level values observed in low-redshift galaxy clusters, we have to consider three different magnetization regimes (St-Onge et al. 2020): the unmagnetized regime, where condition (1) holds, the magnetized “kinetic” regime, in which hard-wall limiters (7) cannot be applied, and the magnetized “fluid” regime, where condition (9) holds and pressure anisotropies are regulated by kinetic instabilities. Therefore, modeling the evolution of the magnetic field during the entire formation of a galaxy cluster requires covering a wide range of values for the ratio of the Larmor radius over the typical length scale of the system, which is far beyond current numerical capabilities.

3. Model for magnetic field evolution in clusters

3.1. Merger tree code and general strategy

In the standard model of cosmology, all structures form from initial primordial density fluctuations. Merger trees constitute a valuable tool for modeling the statistical evolution of dark matter halos. We used the output of such tools, that is the history of the assembly of various dark matter halos, to model the cosmological evolution of plasma parameters (temperature, gas density, turbulence nature, etc.). This allowed us to calculate the growth rate of the small-scale dynamo and therefore enabled us to study the evolution of the magnetic field over cosmological timescales. This is currently unfeasible through numerical simulations, given that the computational power required to model the dynamics of kinetic instabilities surpasses our current capabilities.

To generate dark matter merger trees, we used the modified GALFORM code described in Parkinson et al. (2008). This code is a modification of the GALFORM semi-analytical model developed by Cole et al. (2000). The latter describes cosmic structure formation by employing the mass function from the extended Press-Schechter theory, also implementing star formation feedback, chemical evolution, dust, starburst phases, and so on. For a more detailed description of the GALFORM algorithm, the interested reader is invited to consult the original papers or see Appendix A. In this study, we used the cosmological parameters obtained from Planck (Planck Collaboration VI 2020), that is Ω_m = 0.315, Ω_Λ = 0.685, and H₀ = 67.4 km s⁻¹ Mpc⁻¹. We generated merger trees with a final mass of M_z = 0 = 10¹⁵ M_⊙ and considered values of the maximum redshift of z_max = 1, 1.5, 2, …, 5. We employed a constant sampling in the redshift space of Δz = z_max/300, which easily satisfies the condition P ≪ 1, where P is the average number of progenitors in a given mass range defined by Eq. (A.4). The mass resolution is set to M_res = 10¹² M_⊙. Finally, we generated a total of N_tree = 10³ merger trees for each set of parameters. In Appendix B, we discuss the effect of the number of merger trees on the evolution of various physical quantities. It is shown in Fig. B.1 that increasing the number of merger trees above N_tree ≈ 500 has very little impact on the final curves.

3.2. Radial profiles

For each ith subhalo $H_{z_{\alpha }}^{(i)}$ at a given redshift z_α, we assume that the dark matter distribution, ρ_DM, is determined by the Navarro-Frenk-White profile (Navarro et al. 1996), which is given by

$$\begin{array}{c}\hfill {\rho }_{\mathrm{DM}}(r)=\frac{{\rho }_{\mathrm{s}}}{x{(1+x)}^2},\end{array}$$(10)

where ρ_s is an integration constant and x ≡ r/r_s, where r is the radial coordinate. The scale radius r_s is specified by the concentration parameter c as r_vir ≡ cr_s. Therefore, determining c and ρ_s for each subhalo gives us the dark matter profile directly. We determine c following a similar approach to that used in Dvorkin & Rephaeli (2011), where energy conservation was imposed after multiple mergers or simple matter accretion. A detailed description of this process is presented in Appendix E. We further used the result reported by Ostriker et al. (2005), which states that the gas density and temperature distributions, assuming a hydrostatic equilibrium of a polytropic gas in a dark matter halo, are given by

$$\begin{array}{c}\hfill {\rho }_{\mathrm{g}}(r)={\rho }_0{[1-\frac{\mathrm{\Lambda }}{1+\xi }(1-\frac{ln(1+r/r_{\mathrm{s}})}{r/r_{\mathrm{s}}})]}^{\xi },\end{array}$$(11)

and

$$\begin{array}{c}\hfill T(r)=T_0[1-\frac{\mathrm{\Lambda }}{1+\xi }(1-\frac{ln(1+r/r_{\mathrm{s}})}{r/r_{\mathrm{s}}})],\end{array}$$(12)

where

$$\begin{array}{c}\hfill \mathrm{\Lambda }\equiv \frac{4\pi G{\rho }_{\mathrm{s}}{r}_{\mathrm{s}}^2\mu m_{\mathrm{p}}}{k_{\mathrm{B}}{T}_0},\end{array}$$(13)

with μm_p being the mean molecular weight, ξ ≡ (Γ − 1)⁻¹ the polytropic index, m_p the proton mass, and k_B the Boltzmann constant. We adopted μ = 1/2 and Γ = 1.2. We first determined the constant T₀ by imposing that the temperature at the virial radius should equal the virial temperature of a truncated, single isothermal sphere, which yields

$$\begin{array}{c}\hfill T(r_{\mathrm{vir}})\equiv T_{\mathrm{vir}}=\frac{\mu m_{\mathrm{p}}}{k_{\mathrm{B}}}{v}_{\mathrm{vir}}^2\simeq 3.6\times {10}^5\mathrm{K}{\left(\frac{v_{\mathrm{vir}}}{100\mathrm{km}{\mathrm{s}}^{-1}}\right)}^2,\end{array}$$(14)

where the virial velocity is approximately given by

$$\begin{array}{c}\hfill v_{\mathrm{vir}}^2\approx \frac{\mathit{GM}}{r_{\mathrm{vir}}}·\end{array}$$(15)

Then, the constant ρ₀ is calculated by integration:

$$\begin{array}{c}\hfill {\displaystyle {\int }_0^{r_{\mathrm{vir}}}}4\pi r^2{\rho }_{\mathrm{g}}(r)\mathrm{d}r=M_{\mathrm{b}},\end{array}$$(16)

where M_b is the total mass of baryons in each dark matter subhalo. We assume that M_b/M_DM = 0.1 (e.g., McCarthy et al. 2007).

Our model of turbulent velocity is based on the work of Shi et al. (2018). They investigated the evolution of turbulence in the intracluster medium after major mergers using a set of non-radiative hydrodynamical cosmological simulations and highlighted that the turbulent velocity increases with radius². They also analyzed the radial dependence of the turbulent power spectrum and highlighted that the turbulent injection scale is independent of the radius. Therefore, we assumed the turbulent velocity of each subhalo to be

$$\begin{array}{c}\hfill v_{\mathrm{turb}}(r)\equiv v_0{(1+\frac{r}{r_{\mathrm{vir}}})}^{1/2},\end{array}$$(17)

where v₀ is an integration constant. We calculated the latter making use of the work by Vazza et al. (2012). The turbulent injection length is assumed to be of the form

$$\begin{array}{c}\hfill L_0\equiv \frac{1}{{\alpha }_0}{r}_{\mathrm{vir}},\end{array}$$(18)

where α₀ is a free parameter, and it allows us us to study the effect of various injection scales on our model. The different values tested are, respectively, α₀ = 5, 10, 20, and 50.³

3.3. Plasma parameters deduced from averaged profiles

Now we assign a unique scalar value to each subhalo $H_{z_{\alpha }}^{(i)}$ by averaging over the one-dimensional, spherically symmetric radial profiles described in the previous section. For a given physical quantity Φ ∈ (T, c, ρ_g, v_turb), we calculate the following average ⟨Φ⟩_ρg, weighted by the gas density distribution:

$$\begin{array}{c}\hfill {⟨\mathrm{\Phi }⟩}_{{\rho }_{\mathrm{g}}}\equiv \frac{{\int }_0^{r_{\mathrm{vir}}}\mathrm{\Phi }(r){\rho }_{\mathrm{g}}(r)\mathrm{d}r}{{\int }_0^{r_{\mathrm{vir}}}{\rho }_{\mathrm{g}}(r)\mathrm{d}r}·\end{array}$$(19)

Here, r_vir is the virial radius calculated for each subhalo by the relation

$$\begin{array}{c}\hfill M_{\mathrm{vir}}=\frac{4}{3}\pi r_{\mathrm{vir}}^3{\mathrm{\Delta }}_{\mathrm{c}}{\rho }_{\mathrm{c}},\end{array}$$(20)

where M_vir is the virial mall of the subhalo, ρ_c is the critical density of the Universe at a given redshift, and Δ_c is the overdensity constant. We adopted Δ_c = 200. Then, at redshift z_α, all values resulting from Eq. (19) are again averaged with the mass of each subhalo. All in all, each merger tree is attributed a one-dimensional, time-evolving profile given by

$$\begin{array}{c}\hfill \overline{{⟨\mathrm{\Phi }⟩}_{{\rho }_{\mathrm{g}}}}(z_{\alpha })\equiv \frac{{\sum }_{i\in z_{\alpha }}{M}_i{⟨\mathrm{\Phi }⟩}_{{\rho }_{\mathrm{g}}}}{{\sum }_{i\in z_{\alpha }}{M}_i}·\end{array}$$(21)

The index i ∈ z_α indicates that the sum is performed over all subhalos present at redshift z_α.

This approach results in a more compact representation of the evolution of the different plasma parameters. Indeed, we want to establish a global trend that does not take into account the numerous dynamic phenomena operating within the ICM, such as active galactic nuclei jet cooling flows or galactic wakes, and so on (e.g., Subramanian et al. 2006, and references therein). It also allows us to bypass the difficult task of implementing a radial profile for the magnetic field in each subhalo. For instance, we would need to assume how the magnetic energy of progenitors is distributed after a merging event in the final halo or assume the resulting distribution. Finally, we average the N_trees = 10³ distributions of all physical parameters by fitting all values at a given redshift with a skew-normal distribution. A detailed explanation of this procedure is given in Appendix C.

From the averaged values, $\overline{{⟨\mathrm{\Phi }⟩}_{{\rho }_{\mathrm{g}}}}$, we deduced more plasma parameters while considering a completely ionized plasma made of protons and electrons. The thermal velocity of particle α, where α refers either electrons or ions, is given by

$$\begin{array}{c}\hfill v_{\mathrm{th},\alpha }\equiv \sqrt{\frac{2k_{\mathrm{B}}T}{m_{\alpha }}}·\end{array}$$(22)

The ion-ion collision frequency is calculated assuming momentum exchange only between protons, because of the ratio of one thousand between proton and electron masses. The derivation can be found, for instance, in Fitzpatrick (2014), and leads to

$$\begin{array}{c}\hfill {\nu }_{\mathrm{ii}}=\frac{8}{3}{\left(\frac{m_{\mathrm{e}}}{m_i}\right)}^{1/2}\frac{{\pi }^{3/2}{e}^4{n}_{\mathrm{e}}log\mathrm{\Lambda }}{\sqrt{m_{\mathrm{e}}}}\frac{1}{{\left(k_{\mathrm{B}}T\right)}^{3/2}}·\end{array}$$(23)

The parallel viscosity, which only damps motions parallel to the local direction of the magnetic field, derived in Braginskii (1965), is given by

$$\begin{array}{c}\hfill {\mu }_{\Vert }\equiv v_{\mathrm{th},\alpha }{\lambda }_{\mathrm{mfp},\alpha }.\end{array}$$(24)

The Reynolds number is defined as

$$\begin{array}{c}\hfill \mathrm{Re}\equiv \frac{v_{\mathrm{turb}}{L}_0}{{\mu }_{\Vert }},\end{array}$$(25)

where L₀ is the injection scale of turbulence. Finally, the equipartition field strength is given by

$$\begin{array}{c}\hfill B_{\mathrm{equ}}=\sqrt{4\pi {\rho }_{\mathrm{g}}{v}_{\mathrm{turb}}^2}.\end{array}$$(26)

3.4. Magnetic fields

It is mentioned in Sect. 2 that in the magnetized “fluid” regime, pressure anisotropies could be well regulated by kinetic-scale instabilities, such as the firehose and mirror instabilities. This leads to a relation Re_eff = Re_eff(B), namely that given by Eq. (8). On the other hand, the latter does not hold in the magnetized “kinetic” regime. The instabilities do not scatter particles fast enough to keep the pressure anisotropies near marginal stability. Therefore, the exact dependency of the effective Reynolds number on the magnetic field in this regime is still an open question.

To circumvent this issue, we adopted the following strategy. We constructed an effective Reynolds number for the three different magnetization states. In the unmagnetized regime, the Reynolds number is given by its classical definition, that is Eq. (25). In the magnetized “fluid” regime, we adopted expression (8). However, due to the uncertain nature of the Reynolds number in the magnetized “kinetic” regime, we constructed three different models (denoted L, M, and U, respectively), as shown in Fig. 1. This approach was inspired by the qualitative illustration presented in Fig. 6.1 of St-Onge (2019). Our Model L considers a classical Reynolds number until it reaches the value of the magnetized “fluid” regime. In Model M, a linear interpolation between the unmagnetized and magnetized “fluid” regimes are assumed. Finally, Model U considers a constant value throughout the magnetized “kinetic” regime, equal to the maximum Reynolds number in the fluid regime. By considering these different models, we aim to study certain characteristics of dynamos operating in weakly collisional plasmas, such as the time required to reach equipartition with the turbulent kinetic energy. It should be noted that Model L only constitutes a lower limit on the effective Reynolds number in the situation where the effect of kinetic instabilities is to reduce the effective parallel viscosity by modifying the scattering of particles. In the case where the system is brought back to an equilibrium state by reducing the stretching rate of the magnetic field (and hence increasing the effective viscosity), the value of the effective Reynolds number would decrease. Model U does not represent an upper limit on the Reynolds number either. We established these models arbitrarily to study a wide range of values of the effective Reynolds number.

Fig. 1. Representation of different models for the effective Reynolds number for the following values: ne = 10−3 cm−3, Re = 100, T = 108 K, vturb = 300 km s−1, vth, i = 1000 km s−1, L0 = 100 kpc, and Bequ = 10−6 G. In the unmagnetized regime, the solid black line is the Reynolds number that is not modified by the regulation of pressure anisotropies and is given in Eq. (25). Between the limiting magnetic field strengths given in Eqs. (1) and (9), we construct the three different models, Models U, M, and L. The black dotted line corresponds to results obtained from simulations of weakly collisional plasmas presented in St-Onge et al. (2020).

Next we had to establish the equation that governs the evolution and dynamics of the magnetic field. In this work, we assumed that the magnetic field is amplified by the turbulent dynamo, which consists of the stretching, twisting, and folding of the magnetic field lines by turbulent motions (Kazantsev 1968). If we consider the framework of Kolmogorov turbulence, the rate-of-strain tensor is dominated by motions at the viscous scale (Schekochihin & Cowley 2006), which means |∇u|≈(v_turb/L₀)Re^1/2. On the other hand, in a plasma where a Braginskii-type viscosity is assumed, parallel motions to the magnetic field would be damped and only the parallel component of the rate-of-strain tensor $\widehat{\mathit{b}}\widehat{\mathit{b}}:\mathrm{\nabla }\mathit{u}$ is relevant. Therefore, in using Eq. (6) but replacing Re by Re_eff, we obtain

$$\begin{array}{c}\hfill \frac{1}{B}\frac{\mathrm{d}B}{\mathrm{d}t}\approx \frac{v_{\mathrm{turb}}}{L_0}{\mathrm{Re}}_{\mathrm{eff}}^{1/2},\end{array}$$(27)

where Re_eff depends on the magnetization state and is given in Fig. 1. We solved Eq. (27) for all magnetization states. Once, the magnetic field approaches equipartition with the turbulent velocity field, nonlinear effects should set in and lead to saturation. However, modeling these nonlinear effects is beyond the scope of this work, and we instead simply assumed that the exponential growth ends when the magnetic field reaches the equipartition field strength B_equ.

All in all, we constructed N_tree = 10³ merger trees with N_z = 300 steps between z = 0 and z_max. Once plasma parameters are calculated and averaged according to Eqs. (19) and (21), Eq. (27) is solved between two redshift steps. Our choice of time resolution between two given redshift steps is explained in detail in Appendix D.

4. Results

4.1. Basic plasma parameters from 1D profiles

Figure 2 showcases the temporal evolution of various plasma parameters (Sect. 3.3) for different z_max values associated with merger trees generated using the modified GALFORM algorithm. Initially, we observe a similarity in the trends of the temperature and turbulent velocity curves. The algorithm demonstrates a decline in these curves, by approximately 10 to 20% of the total runtime, followed by a local minimum and a subsequent steady increase until it reaches z = 0. Additionally, we notice that higher values of z_max lead to steeper local minima. While the underlying cause of these curve shapes remains unclear, we can provide some explanations. Firstly, it is important to note that the merger trees have been constructed with a uniform mass resolution of M_res = 10¹² M_⊙, representing a thousandth of the final cluster’s mass. However, the GALFORM algorithm incorporates two stop conditions: reaching a predetermined value of z_max, or the subhaloes reaching the resolution mass during the mass-splitting procedure, regardless of the redshift. Consequently, a higher z_max corresponds to fewer subhaloes in the tree near this redshift. Nonetheless, deviations from the mean values of temperature and turbulent velocity at z = 0 remain minimal. Furthermore, the concentration parameter c is computed at each node of the tree by considering energy conservation, accounting for accreted matter. When the merger tree extends to higher redshift values, the accretion of additional matter can substantially alter the final value of the concentration parameter, which could explain why c values of the curves with 1 ≲ z_max ≲ 2 at z = 0 are higher. However, it is observed from Fig. 2 that variations in c values at z = 0 do not seem to impact the electron density curves. Remarkably, the c value converges to approximately c ≃ 1.2 for z ≳ 2.5, in relatively good agreement with the concentration–mass relation derived in Biviano et al. (2017).

Fig. 2. Evolution of various plasma quantities as a function of redshift and for different values of zmax. Each curve represents the average value of the quantity over Ntree = 103 merger trees (see Appendix C) for which a redshift-dependent, one-dimensional profile has been calculated according to Eqs. (19) and (21). The age of the Universe tU at each redshift is plotted on the second x-axis. In order to assess the validity of our model, we added observed values for the temperature and the thermal electron density (see the detailed list in the main text).

In the preceding paragraph, several concerns regarding the validity of our magnetic field model were raised, as numerical artifacts arise for redshifts close to z_max. As previously stated, both the temperature and turbulent velocity curves exhibit a decreasing trend that precedes a global minimum and then displays a constant increase until z = 0. However, it has been observed that this trend strongly relies on the initial conditions set by z_max, resulting in a numerical effect in our model. In order to mitigate the potential impact of these artifacts on our conclusions, we decided to impose constraints on the value of z_max. Figure 2 clearly indicates that all curves converge toward a common value for redshift values z ≳ 2. Therefore, we opted to exclude merger trees with z_max ≤ 2 for further analysis. Although selecting the merger tree with the highest z_max seems to be a plausible option, it is worth noting that the number of dark matter subhaloes decreases as z_max increases due to the way the GALFORM algorithm operates. Therefore, from now on, we only use the merger tree with z_max = 4 to implement our magnetic field amplification model. Additionally, we initiate the dynamo process at redshift values equal to or below z_start ≤ 2, ensuring that the influence of the mass resolution remains marginal.

Finally, in order to assess the robustness of our model, we plot observational values of various physical quantities for nine different clusters in Fig. 2. It should be noted that this selection of clusters does not include only clusters in the same dynamical state, some are undergoing merging processes, while others are in a relaxed state. For instance, Miranda et al. (2008) combined X-ray observations and strong lensing analysis of RX J1347.5−1145, and revealed a complex structure. They suggested a merger scenario between dark matter subclumps that would account for discrepancies with mass estimates from the virial theorem. Moreover, Barrena et al. (2002) provided evidence of a major collision between 1E 0657−56 and a subcluster, and Maurogordato et al. (2008) provided optical observations suggesting that Abell 2163 had undergone a recent (≈0.5 Gyr) collision. For the temperature, we plotted results from X-ray observations for the Coma cluster, Abell 1795, PKS 0745−19, SW J1557+3530, and SW J0847+1331 from Moretti et al. (2011) and references therein. We also plotted the temperature from X-ray observations from Wallbank et al. (2022) for Abell 2163, Abell 754, RX J1347.5−1145, and 1E 0657−56. Our model results in an average temperature at z = 0, which is in good agreement with most clusters of the observational sample. We also present the thermal electron density values for the Coma cluster, both in its central region and its periphery, at r = r₂₀₀. These values were obtained from the X-ray morphology Churazov et al. (2021) using the data from SRG/eROSITA. This enables us to assert that the typical plasma values expected at low redshift are overall in accordance with our model.

4.2. Magnetic field amplification

In this section, we present the magnetic field amplification obtained from our model that is based on an effective Reynolds number. It is worth recalling that we integrated our dynamo model into the merger tree with z_max = 4 and initiated the amplification process at various redshifts between z_start = 2 and z_start = 0.1.

4.2.1. Example of magnetic field evolution

Figure 3 shows an example of the magnetic field evolution for the three effective Reynolds number models, detailed in Sect. 3.4, for α₀ = 20 (see Eq. (18) for the model of the turbulent forcing scale), z_start = 1.5, and B₀ = 10⁻²⁰ G. The gray curves represent the magnetic field evolution of all individual N_tree = 10³ merger trees. In red, we show the skew-normal average median with a 90% confidence interval (see Appendix C). We note that the estimated error during the exponential phase of Model L is larger than that of the other models. This can be partially explained by the fact that, in Model L, kinetic instabilities only modify the Reynolds number once the magnetic field amplitude reaches the order of nG. Compared to Models M and U, in Model L, the Reynolds number remains at a value several orders of magnitude lower for a longer time, allowing the various curves to spread more until the Reynolds number increases drastically, when B exceeds the threshold given by Eq. (9). In all three cases, Fig. 3 demonstrates that once an amplitude of a few nG is achieved, the amplification of the magnetic field toward equipartition is nearly instantaneous.

Fig. 3. Evolution of magnetic field for α0 = 20 and zstart = 1.5 and for the three models for the effective Reynolds number. The gray curves are the magnetic field evolution of all Ntree = 103 merger trees. The red curve indicates the average obtained by fitting all curves at each redshift value with a skew-normal distribution (see Appendix C). The colored area shows the 90% confidence interval for such a distribution. The age of the Universe tU at each redshift is plotted on the second x-axis.

Figure 4 illustrates the time evolution of the magnetic field for the three distinct effective Reynolds number models, detailed in Sect. 3.4, for a specific value of α₀ = 20, starting at redshifts of z_start = 0.5, 1, and 1.5.

Fig. 4. Evolution of magnetic field amplitude over redshift for different models of effective Reynolds number. The curves and shaded areas correspond, respectively, to the average and the 90% confidence interval given by the averaging process described in Appendix C. The age of the Universe tU at each redshift is plotted on the second x-axis. Each row corresponds to a different model, and each column corresponds to a different starting redshift of the dynamo. Model L (top row) has the lowest growth rate, while Models M and U (middle and bottom rows, respectively) show faster growth rates. The figure also shows that the dynamo was much faster in the past (higher redshift values).

Figure 4 illustrates the same evolution of the magnetic field as Fig. 3, but this time for z_start = 0.5, 1, 1.5 and for initial seed field values ranging from log₁₀(B [G]) = −20 to −10. Firstly, it is evident that for Models M and U, the curves coincide for B₀ ≳ 10⁻¹⁶ G and, in particular, the time required to reach equipartition with the turbulent velocity field is independent of B₀. Therefore, we chose to omit the error bars in our subsequent analyses described below. It also appears that the growth rate of the magnetic field becomes faster as the redshift increases. This effect can only be justified by the evolution of various plasma parameters incorporated in Eq. (27) for the growth rate.

4.2.2. Time to equipartition

Figure 5 shows the values of redshift z_equ at which the magnetic energy density reaches at least 90% of the kinetic energy density of turbulence, for the three models of Re_eff, and for different values of α₀, B₀, and z_start. Figure 6 shows the minimum redshift z_min at which the dynamo has to start such that equipartition with the turbulent velocity field can be reached, for the three effective Reynolds number models and for various values of α₀. Firstly, it can be observed that equipartition can be reached in all the different configurations that we implemented. Moreover, it is clear that Models M and U are indistinguishable for log₁₀(B₀ [G]) ≳ −17. Furthermore, below these values, the differences between their respective curves are not significant.

Fig. 5. Redshift zequ at which magnetic field reaches at least 90% of equipartition with turbulent velocity field; for different values of initial magnetic seed field. Each row corresponds to the different models of the effective Reynolds number in the magnetized “kinetic” regime (Sect. 3.4). Each column corresponds to a different turbulent injection scale.

The important conclusions that can be drawn from these results are as follows. Firstly, assuming that the modification of effective ion-ion collisionality is the sole physical quantity modified by kinetic instabilities, it can reasonably be accepted that Model L represents a lower limit of z_min. Additionally, even though Model U does not represent an upper limit to the effective Reynolds number, it can be assumed that any faster model would not result in any changes to the values presented in Fig. 6. In summary, any effective model for the Reynolds number lying between the curves associated with Models M and U would lead to an explosive dynamo so fast that their evolution could not be distinguished with our resolution in redshift. On the other hand, considering a model between the curves L and M would indeed lead to changes in the values of z_min compared to those presented in Fig. 6. Finally, we are limited in the conclusions we can draw due to our simple limitation in modeling the dynamics of the turbulent field, particularly its injection scale. If α₀ were to evolve, especially through certain mechanisms of turbulent decay, our results would obviously be modified. We leave a more comprehensive modeling of the turbulent velocity field for future work.

Fig. 6. Minimum redshift zmin from which a dynamo has to start to take the magnetic field B0 up to equipartition with the turbulent velocity field at z = 0; for different values of the initial magnetic seed field and for each model for the effective Reynolds number. Each panel corresponds to a specific value of the injection scale defined by Eq. (18).

5. Discussion

Using merger trees has enabled us to study the evolution of certain plasma quantities, including the magnetic field, without relying on computationally expensive simulations. However, such an approach presents a number of uncertainties. Firstly, we only used the modified GALFORM algorithm, as it was calibrated to the mass function of the Millennium simulations. We do not know (yet) if our conclusions would change when considering another merger tree algorithm.

When applying the merger tree algorithm, we have to pay special attention to numerical artifacts. In particular, we constrained the value of z_max (redshift at which the merging process starts) to z_max = 4, so that our results are not influenced by numerical effects related to mass resolution (see Fig. 2). Furthermore, we decided to start the dynamo process from z_start ≃ 2 or below. Indeed, below z_start = 2, all curves with z_max ≥ 4 are approximately identical. Therefore, our ability to study dynamo processes at z ≳ 2 depends on the model for the concentration parameter considered.

Also, the effective Reynolds number has to be chosen carefully, mainly because Models L and U do not constitute a lower and upper bound per se. This may be susceptible to introducing biases in the minimum redshift, z_min, at which the magnetic field amplification needs to start in order to reach equipartition at z = 0. However, our results have shown that Models M and U produce nearly identical results. This is not surprising, because even if Models M and U are not the same, they both differ from the classical Reynolds number by (many) orders of magnitude.

Furthermore, confining the evolution of the turbulent injection scale appears to be one of the few necessary conditions to constrain the evolution of the magnetic field. It should be noted that other sources of turbulence are likely to exist during certain phases of the formation of a cluster, potentially with different forcing length scales. Other sources include, for example, subcluster and galactic wakes (e.g., Subramanian et al. 2006) and jets from active galactic nuclei (e.g., Fujita et al. 2020).

Finally, for the scenario presented in Fig. 3, for instance, the magnetic field strength increases by 14 orders of magnitude, while the temperature increases by one order of magnitude. In this specific example, the β parameter changes from β ≈ 10²⁸ at z = 1.5 to β ≈ 10² at z = 0. Therefore, we do not expect direct feedback on cluster evolution from the magnetic field through magnetic pressure. However, even at times where β ≫ 1, the strength and direction of the viscosity are affected by the instabilities. This implies that indirect feedback from the magnetic field on the flow can occur, leading to effects that are not captured by classical hydrodynamic simulations (see also the discussion in St-Onge & Kunz 2018).

6. Conclusions

In this work, we used the modified GALFORM code (Parkinson et al. 2008) to investigate the amplification of magnetic fields in the intracluster medium during the formation of galaxy clusters by successive mergers of dark matter halos. We focused on clusters with a typical mass of M = 10¹⁵ M_⊙.

To this end, we implemented one-dimensional profiles for the distribution of dark matter, baryonic gas, temperature, and turbulent velocity, assuming spherical symmetry and following the methodology proposed by Dvorkin & Rephaeli (2011) and Johnson et al. (2021). For a given merger tree realization, the profiles of all subhalos at a given redshift were averaged in order to calculate relevant plasma parameters as a function of redshift. This process was performed for N_tree = 10³ merger trees to establish statistical properties. We compared the temperature at z = 0 obtained from our model with observed values from a sample of selected clusters presented in Moretti et al. (2011) and Wallbank et al. (2022). We also compared the electron density values of our model with SRG/eROSITA observations of the Coma cluster of Churazov et al. (2021). Overall, we found our model to be in good agreement with the observed values from a sample of low-redshift clusters.

We then constructed a semi-analytical model for the amplification of the magnetic field that incorporates the effects of pressure-anisotropy-driven kinetic instabilities (firehose and mirror) for different magnetization regimes (Melville et al. 2016; Rincon et al. 2016; St-Onge & Kunz 2018; St-Onge et al. 2020). Next, we estimated the growth rate of the magnetic field for different values of the strength of the magnetic seed field and forcing length scales of turbulence. This allowed us to determine the redshift at which equipartition is reached in different scenarios. In practice, we define equipartition as the time when the magnetic energy density reaches at least 90% of the turbulent kinetic energy density. Specifically, for an injection scale of L₀ = r_vir/20, where r_vir is the virial radius of the dark matter halo, the magnetic field evolves in the following way: for our slowest dynamo (effect Reynolds number Model L), the dynamo amplification for seed field values of log₁₀(B₀ [G]) = −20 and −10 has to start at least from, respectively, z_start ≃ 1 and 0.2 to reach equipartition until the present day. For the other faster scenarios (Models M and U), those values are reduced to z_start ≃ 0.3 and ≃0.1, respectively. Although our Model U cannot be considered an upper limit for the effective Reynolds number, our Model L can constitute a lower limit, provided that the increase of ion-ion collisionality is the sole effect created by the kinetic instabilities to bring the system back to marginal stability.

Overall, our study demonstrates that merger trees can be a valuable tool for constraining the evolution of magnetic fields over cosmological times, especially in galaxy clusters. In particular, they allow us to probe kinetic plasma effects in the history of the ICM that are inaccessible with state-of-the-art cosmological simulations. Although our approach is not as accurate as fully kinetic plasma simulations, it could potentially constitute a computationally efficient alternative for constraining the dependence of the effective Reynolds number on the magnetic field in the magnetic “kinetic” regime when combined with future and more resolved radio observations.

Acknowledgments

We are very grateful to Matthew Kunz and Jonathan Squire for providing us with all the necessary ingredients for modeling magnetic fields in the intracluster medium, during the Nordita program “Magnetic Field Evolution in Low Density or Strongly Stratified Plasmas” in 2022. Undoubtedly, our work would not have been possible without their contribution to the physics of weakly collisional plasmas. We further thank Abhijit B. Bendre for his numerous comments on this study and the referee for a very constructive report. Finally, we acknowledge the support from the Swiss National Science Foundation under Grant No. 185863.

References

Appendix A: GALFORM & modified GALFORM models

The GALFORM model (Cole et al. 2000) is based on the extended Press-Schechter theory (Press & Schechter 1974), which considers the following conditional mass function:

$$\begin{array}{cc}\hfill f(M_1|M_2)dln{M}_1=& \sqrt{\frac{2}{\pi }}\frac{{\sigma }_1^2({\delta }_1-{\delta }_2)}{{({\sigma }^2-{\sigma }_2^2)}^{3/2}}\hfill \\ \hfill & exp(-\frac{1}{2}\frac{{({\delta }_1-{\delta }_2)}^2}{{\sigma }_1^2-{\sigma }_2^2})\left|\frac{dln\sigma }{dln{M}_1}\right|dln{M}_1,\hfill \end{array}$$(A.1)

where f(M₁|M₂) represents the fraction of mass from halos of mass M₂ at redshift z₂ that is contained in progenitor halos of mass M₁ at earlier redshift z₁. The parameter δ_i denotes the linear overdensity threshold for spherical collapse at redshift z_i. The quantity σ(M) denotes the rms linear density fluctuation extrapolated at z = 0 in spheres containing mass M. For simplicity, we adopted σ_i(M)≡σ_i. Taking the limit z₁ → z₂, we obtain

$$\begin{array}{cc}& {\frac{\mathit{df}}{d{z}_1}|}_{z_1=z_2}dln{M}_{1}d{z}_1\hfill \\ \hfill & =\sqrt{\frac{2}{\pi }}\frac{{\sigma }_1^2}{{({\sigma }_1^2-{\sigma }_2^2)}^{3/2}}\frac{d{\delta }_1}{d{z}_1}\left|\frac{dln{\sigma }_1}{dln{M}_1}\right|dln{M}_{1}d{z}_1,\hfill \end{array}$$(A.2)

which allows us to find the mean number of halos of mass M₁ into which a halo of mass M₂ splits when we increase redshift upwards with a redshift step dz₁:

$$\begin{array}{c}\hfill \frac{\mathit{dN}}{d{M}_1}=\frac{1}{M_1}\frac{\mathit{df}}{d{z}_1}\frac{M_2}{M_1}d{z}_1.\end{array}$$(A.3)

Assuming a resolution mass M_res, we can also determine the average number of progenitors with masses M₁ in the interval M_res < M₁ < M₂/2, which is given by

$$\begin{array}{c}\hfill P={\displaystyle {\int }_{M_{\mathrm{res}}}^{M_2/2}\frac{\mathit{dN}}{d{M}_1}d{M}_1,}\end{array}$$(A.4)

along with the fraction of mass of the final object in progenitors below M_res, given by

$$\begin{array}{c}\hfill F={\displaystyle {\int }_0^{M_{\mathrm{res}}}\frac{\mathit{dN}}{d{M}_1}\frac{M_1}{M_2}d{M}_1.}\end{array}$$(A.5)

The GALFORM algorithm operates as follows. We start with the redshift and the mass of the final halo. Then we consider a redshift step dz₁ that satisfies P ≪ 1 (meaning that a halo is unlikely to have more than two progenitors at redshift z + dz). A random number R ∈ [0, 1] is generated. If R > P the main halo is not split, and its mass is reduced to M₂(1 − F) to account for the accreted mass. If R ≤ P, then a random mass M_res > M₁ > M₂/2 is generated to produce two new halos with masses of M₁ and M₂(1 − F)−M₁. The same process is repeated until the maximum redshift or the resolution mass is reached. However, there are some issues with such an approach. One of them is that the conditional mass function (A.1) does not match what is found in N-body simulations. Yet, statistical properties produced by the GALFORM algorithm have similar trends with those of merger trees constructed from high-resolution N-body simulations, with an error on mass progenitors that increases with redshift. In that spirit, Parkinson et al. (2008) replaced the initial statistics Eq. (A.3) with

$$\begin{array}{c}\hfill \frac{\mathit{dN}}{d{M}_1}\to \frac{\mathit{dN}}{d{M}_1}G({\sigma }_1/{\sigma }_2,{\sigma }_2/{\sigma }_2),\end{array}$$(A.6)

where G is a “perturbative” function that matches the generated merger tree with simulations. They also made the following assumption:

$$\begin{array}{c}\hfill G({\sigma }_1/{\sigma }_2,{\sigma }_2/{\sigma }_2)=G_0{\left(\frac{{\sigma }_1}{{\sigma }_2}\right)}^{{\gamma }_1}{\left(\frac{{\delta }_2}{{\sigma }_2}\right)}^{{\gamma }_2},\end{array}$$(A.7)

for which they obtained G₀ = 0.57, γ₁ = 0.38 and γ₂ = −0.01. Finally, they obtained halo abundances consistent with the Sheth-Tormen mass function (see Parkinson et al. 2008, and references therein). In our work, we used the same values for G₀, γ₁, and γ₂.

Appendix B: Effect number of generated trees

Figure B.1 shows the effect of the number of generated merger trees, on the skew-normal average (Sect. C) calculated for the temperature, the baryonic gas density, and the turbulent velocity at three different redshifts. The solid horizontal lines correspond to the value calculated with N_tree = 10³ merger trees. It is evident that convergence is reached, even for N_tree ≳ 10^2.5 merger trees. This justifies our choice of setting N_tree = 10³.

Fig. B.1. Skew-normal average of temperature, gas density, and turbulent velocity at three different values of redshift as a function of the number of generated merger trees. The solid lines represent the values of the physical quantities for Ntree = 103.

Appendix C: Skew-normal averaging process

Here, we present our method for averaging the time evolution curves of all N_tree = 10³ merger trees for all physical quantities. Initially, one might consider using the conventional arithmetic average to estimate the overall behavior of these curves at each specific redshift. However, this approach presents a significant flaw when analyzing the evolution of the magnetic field. By examining the left panel of Fig. 3, it becomes evident that at approximately z ≃ 0.75, while a few magnetic field curves have reached equipartition, the majority are still in the phase of exponential growth, with values between B ≈ 10⁻¹⁴ G and B ≈ 10⁻¹² G. Consequently, applying the classical arithmetic average to this dataset would yield a value close to equipartition, which does not accurately represent the actual trend exhibited by all the curves.

We decided to adopt the following strategy. At each given redshift, we fit the distribution of the values of all curves with a skew-normal distribution. However, special attention has to be paid to the magnetic field, as shown in Fig. C.1. Here, we show the histograms of different quantities, where the y-axis is normalized so that the total area of the histogram equals one. The top panel shows the temperature distribution at redshift z = 0.29, and the middle one shows the distribution of the magnetic field at the same redshift. If we were to fit the whole magnetic field sample, a skew-normal distribution would not necessarily be suitable as two distinct maxima appear. Therefore, we adopted the following convention. Whenever the number of trees that have reached equipartition is dominant, we fit the region delimited by dashed lines on the middle panel, which is shown on the bottom panel. However, if this number is smaller than the trees outside this region, we fit the trees in the latter. Finally, we note that the situation where the two histograms, distinguishing the equipartition magnetic field from other values, have the same number of occurrences is extremely rare.

Fig. C.1. Illustration of our method using a skew-normal distribution to average different physical quantities at a given redshift, for many merger trees. (top) Distribution of the temperature values from all merger trees at redshift z = 0.75. The blue line represents the skew-normal distribution used to fit the data. The dashed line and the blue-shaded line, respectively, are the distribution’s median (that is used as an estimator for the average of all curves) and its corresponding 90% confidence interval. (middle) Distribution of magnetic field values at this same redshift. The two dashed black lines indicate the region that is fit with the skew-normal distribution. (bottom) Fitting within the region in the middle panel. Those three panels were computed with Reeff Model L, for log10(B0 [G]) = −20, α0 = 20, and zstart = 1.5.

Appendix D: Time resolution of the dynamo model

As described in Sect. 3.4, we solved Eq. (27) between every redshift step of our merger tree (we recall that N_z = 300). Redshift is converted to time using the cosmology calculator of Wright (2006) with the cosmological parameters given in Sect. 3.1. However, considering only the two values of time between two redshift steps is not efficient enough to highlight the complexity of the different stages of the dynamo, such as nonlinear effects when the magnetic field is approaching equipartition. Therefore, we created a discretized time grid, between two given redshift steps, whose time step is kept constant (meaning that the number of time grid cells between two different consecutive reshifts will vary). The constant time step is given as follows. We can see that Eq. (27) is of the form B(t)∝e^Γ_efft, where ${\mathrm{\Gamma }}_{\mathrm{eff}}\equiv (v_{\mathrm{turb}}/L_0){\mathrm{Re}}_{\mathrm{eff}}^{1/2}$ is the growth rate of the magnetic field. However, given our construction for the effective Reynolds number in Sect. 3.4, it is obvious that Re_eff ≥ Re, which would yield a smaller timescale of the growth rate than with the classical Reynolds number. Therefore, in order to avoid over-discretizing our grid, we considered the timescale coming from the classical Reynolds number, Γ⁻¹, to be a minimum timescale that should not be exceeded. Therefore, the constant value Δt of all grid cells N_Δt contained in our discretized time between two redshift steps is given by

$$\begin{array}{c}\hfill \mathrm{\Delta }t\equiv \frac{{\mathrm{\Gamma }}^{-1}}{{\alpha }_B}=\frac{1}{{\alpha }_B}{\left(\frac{v_{\mathrm{turb}}}{L_0}{\mathrm{Re}}^{1/2}\right)}^{-1},\end{array}$$(D.1)

where α_B is a constant value.

We determine the latter by plotting the evolution of the magnetic field, for a fixed set of merger-tree parameters, and for different values of α_B. Figure D.1 shows such results for B₀ = 10⁻¹⁵ G and α₀ = 20. It appears that when α_B ≈ 10^3.5, all curves converge. Therefore, in our work, we adopted a constant value of α_B = 10⁴. We note that we checked that this value of α_B does not produce other results for all other parameters of B₀ and α₀.

Fig. D.1. Magnetic field evolution for different values of αB. Each column corresponds to a different value of redshift at which the dynamo starts. Each row corresponds to a different model for the effective Reynolds number described in Sect. 3.4.

Appendix E: Estimation of the c parameter

We describe our method to determine the concentration parameter in every subhalo of the generated merger trees here. We started by initiating a value of c for subhalos at z_max, with

$$\begin{array}{c}\hfill c(M)=4.67{\left(\frac{M}{{10}^{14}{h}^{-1}{M}_{\odot }}\right)}^{-0.11},\end{array}$$(E.1)

which was proposed by Neto et al. (2007) and where h is the Hubble constant in units of 100 km s⁻¹ Mpc⁻¹ (see Springel et al. 2005). This choice of the initial condition given in Eq. E.1 is somehow arbitrary, but after trying different values of c for the earliest halos, we found that the final result does not depend on this particular choice. We calculated the total energy of every subhalo using the same expressions adopted in Johnson et al. (2021), which is the sum of its kinetic and potential energy, given, respectively, by

$$\begin{array}{c}\hfill T=2\pi [r_{\mathrm{vir}}^3\rho (r_{\mathrm{vir}}){\sigma }^2(r_{\mathrm{vir}})+{\displaystyle {\int }_0^{r_{\mathrm{vir}}}}GM(r)\rho (r)r\mathrm{d}r]\end{array}$$(E.2)

and

$$\begin{array}{c}\hfill W=-\frac{G}{2}{\displaystyle [{\int }_0^{r_{\mathrm{vir}}}\frac{M^2(r)}{r^2}\mathrm{d}r+\frac{M^2(r_{\mathrm{vir}})}{r_{\mathrm{vir}}}].}\end{array}$$(E.3)

Here, σ(r) is the velocity dispersion that is determined by the Jeans equation:

$$\begin{array}{c}\hfill \frac{d(\rho {\sigma }^2)}{\mathit{dr}}=-\rho (r)\frac{GM(r)}{r^2}·\end{array}$$(E.4)

Therefore, each halo was attributed a total energy of E_h = T + W, which depends on parameters ρ_s and r_s. When a subhalo was the result of N merging subhalos, we computed the energy E_h, f = T_f + W_f of the resulting halo, which we equalized with the energy of the system of all progenitors. We expressed the total energy E_N of a system of N subhalos using the same strategy presented in Dvorkin & Rephaeli (2011), which is

$$\begin{array}{c}\hfill E_N={\displaystyle \sum _{i=1}^N{E}_{\mathrm{h},i}+\sum _{i\ne j}{U}_{\mathit{ij}}+\sum _{i\ne j}{U}_{\mathrm{acc},ij},}\end{array}$$(E.5)

where

$$\begin{array}{c}\hfill U_{\mathit{ij}}\equiv -\frac{G{M}_i{M}_j}{d_{0,ij}}\end{array}$$(E.6)

is the gravitational energy of two halos at the largest point of their separation, with d_0, ij being the distance at which subhalos i and j are bound, and where

$$\begin{array}{c}\hfill U_{\mathrm{acc},ij}\equiv -\frac{G(M_i+M_j)\mathrm{\Delta }M}{r_{\mathrm{vir},\mathrm{f}}}\end{array}$$(E.7)

is the energy component of the accreted matter; r_vir, f stands for the virial radius of the resulting halo, and ΔM ≡ M_f − ∑_iM_i is the mass that was not resolved during the merging event, below the value of M_res. We used the same value of d₀ as the one adopted in Dvorkin & Rephaeli (2011), which is d_0, ij = 5(r_vir, i + r_vir, j). Our strategy was then to compute the energy for each subhalo (whether it has one or multiple progenitors) and find the value of c that matches the energy conservation with its progenitors. Once c was calculated, the parameter ρ_s was obtained by integrating the dark matter density to match the mass of each subhalo.

All Figures

Fig. 1. Representation of different models for the effective Reynolds number for the following values: ne = 10−3 cm−3, Re = 100, T = 108 K, vturb = 300 km s−1, vth, i = 1000 km s−1, L0 = 100 kpc, and Bequ = 10−6 G. In the unmagnetized regime, the solid black line is the Reynolds number that is not modified by the regulation of pressure anisotropies and is given in Eq. (25). Between the limiting magnetic field strengths given in Eqs. (1) and (9), we construct the three different models, Models U, M, and L. The black dotted line corresponds to results obtained from simulations of weakly collisional plasmas presented in St-Onge et al. (2020).

In the text

Fig. 2. Evolution of various plasma quantities as a function of redshift and for different values of zmax. Each curve represents the average value of the quantity over Ntree = 103 merger trees (see Appendix C) for which a redshift-dependent, one-dimensional profile has been calculated according to Eqs. (19) and (21). The age of the Universe tU at each redshift is plotted on the second x-axis. In order to assess the validity of our model, we added observed values for the temperature and the thermal electron density (see the detailed list in the main text).

In the text

Fig. 3. Evolution of magnetic field for α0 = 20 and zstart = 1.5 and for the three models for the effective Reynolds number. The gray curves are the magnetic field evolution of all Ntree = 103 merger trees. The red curve indicates the average obtained by fitting all curves at each redshift value with a skew-normal distribution (see Appendix C). The colored area shows the 90% confidence interval for such a distribution. The age of the Universe tU at each redshift is plotted on the second x-axis.

In the text

Fig. 4. Evolution of magnetic field amplitude over redshift for different models of effective Reynolds number. The curves and shaded areas correspond, respectively, to the average and the 90% confidence interval given by the averaging process described in Appendix C. The age of the Universe tU at each redshift is plotted on the second x-axis. Each row corresponds to a different model, and each column corresponds to a different starting redshift of the dynamo. Model L (top row) has the lowest growth rate, while Models M and U (middle and bottom rows, respectively) show faster growth rates. The figure also shows that the dynamo was much faster in the past (higher redshift values).

In the text

	Fig. 5. Redshift zequ at which magnetic field reaches at least 90% of equipartition with turbulent velocity field; for different values of initial magnetic seed field. Each row corresponds to the different models of the effective Reynolds number in the magnetized “kinetic” regime (Sect. 3.4). Each column corresponds to a different turbulent injection scale.
In the text

	Fig. 6. Minimum redshift zmin from which a dynamo has to start to take the magnetic field B0 up to equipartition with the turbulent velocity field at z = 0; for different values of the initial magnetic seed field and for each model for the effective Reynolds number. Each panel corresponds to a specific value of the injection scale defined by Eq. (18).
In the text

	Fig. B.1. Skew-normal average of temperature, gas density, and turbulent velocity at three different values of redshift as a function of the number of generated merger trees. The solid lines represent the values of the physical quantities for Ntree = 103.
In the text

Fig. C.1. Illustration of our method using a skew-normal distribution to average different physical quantities at a given redshift, for many merger trees. (top) Distribution of the temperature values from all merger trees at redshift z = 0.75. The blue line represents the skew-normal distribution used to fit the data. The dashed line and the blue-shaded line, respectively, are the distribution’s median (that is used as an estimator for the average of all curves) and its corresponding 90% confidence interval. (middle) Distribution of magnetic field values at this same redshift. The two dashed black lines indicate the region that is fit with the skew-normal distribution. (bottom) Fitting within the region in the middle panel. Those three panels were computed with Reeff Model L, for log10(B0 [G]) = −20, α0 = 20, and zstart = 1.5.

In the text

	Fig. D.1. Magnetic field evolution for different values of αB. Each column corresponds to a different value of redshift at which the dynamo starts. Each row corresponds to a different model for the effective Reynolds number described in Sect. 3.4.
In the text