OPEN ACCESS

On the Small-scale Turbulent Dynamo in the Intracluster Medium: A Comparison to Dynamo Theory*

Ulrich P. Steinwandel¹ orcid, Ludwig M. Böss2,3, Klaus Dolag2,4 and Harald Lesch²

Published 2022 July 08 • © 2022. The Author(s). Published by the American Astronomical Society.

The Astrophysical Journal, Volume 933, Number 2

Author e-mails: usteinwandel@flatironinstitute.org

Author affiliations
¹ Center for Computational Astrophysics, Flatiron Institute, 162 5th Avenue, New York, NY 10010, USA
² University Observatory Munich, Scheinerstraße 1, D-81679 Munich, Germany
³ Excellence Cluster ORIGINS, Boltzmannstraße 2, D-85748 Garching, Germany
⁴ Max Planck Institute for Astrophysics, Karl-Schwarzschildstraße 1, D-85748 Garching, Germany

Dates
Received 2021 August 18; Revised 2022 April 19; Accepted 2022 May 16

Citation
Steinwandel et al. 2022 APJ 933 131

Unified Astronomy Thesaurus concepts
Galaxy clusters (584); Cosmic magnetic fields theory (321); Extragalactic magnetic fields (507); Magnetohydrodynamical simulations (1966); Intracluster medium (858); Magnetic fields (994)

DOI: https://doi.org/10.3847/1538-4357/ac715c

Abstract

We present non-radiative, cosmological zoom-in simulations of galaxy-cluster formation with magnetic fields and (anisotropic) thermal conduction of one massive galaxy cluster with Mvir ∼ 2 × 10¹⁵M⊙ at z ∼ 0. We run the cluster on three resolution levels (1×, 10×, 25×), starting with an effective mass resolution of 2 × 10⁸M⊙, subsequently increasing the particle number to reach 4 × 10⁶M⊙. The maximum spatial resolution obtained in the simulations is limited by the gravitational softening reaching ϵ = 1.0 kpc at the highest resolution level, allowing one to resolve the hierarchical assembly of the structures in fine detail. All simulations presented are carried out with the SPMHD code gadget3 with an updated SPMHD prescription. The primary focus of this paper is to investigate magnetic field amplification in the intracluster medium. We show that the main amplification mechanism is the small-scale turbulent dynamo in the limit of reconnection diffusion. In our two highest resolution models we start to resolve the magnetic field amplification driven by the dynamo and we explicitly quantify this with the magnetic power spectra and the curvature of the magnetic field lines, consistent with dynamo theory. Furthermore, we investigate the ∇ · B = 0 constraint within our simulations and show that we achieve comparable results to state-of-the-art AMR or moving-mesh techniques, used in codes such as enzo and arepo. Our results show for the first time in a cosmological simulation of a galaxy cluster that dynamo action can be resolved with modern numerical Lagrangian magnetohydrodynamic methods, a study that is currently missing in the literature.

Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

1. Introduction

Magnetic fields are observed across all scales within the universe, from the interstellar medium (ISM) on the scales of molecular clouds (e.g., Heiles & Crutcher 2005; Crutcher 2012; Clark et al. 2014, 2019; Hu et al. 2019; Sullivan et al. 2021) over the large-scale field structure in galaxies (e.g., Basu & Roy 2013; Beck & Krause 2005; Greaves et al. 2000; Jones et al. 2020; Lacki & Beck 2013; Robishaw et al. 2008; Tabatabaei et al. 2008; Watson & Wyld 2001) to the intracluster medium (ICM; e.g., Bonafede et al. 2009, 2010, 2013; Böhringer et al. 2016; Clarke et al. 2001; Hu et al. 2020).

While there are significant differences in the structure of the magnetic field within the ISM, the intergalactic medium, and the ICM there seems to be some observational evidence that the magnetic field strength within these very different components is in the order of a few microgauss (see, e.g., Beck 2015; Crutcher 2012; van Weeren et al. 2019; for reviews on magnetic fields on galaxy, ISM, and ICM scales) and seems to be correlated on galaxy and even galaxy-cluster scales.

Specifically, for galaxy clusters one can determine the magnetic field strength from the Faraday rotation measurement (RM) of radio galaxies that are located in the foreground and background of the cluster of interest (e.g., Brentjens & de Bruyn 2005; Burn 1966; Clarke 2004; Murgia et al. 2004; van Weeren et al. 2019). The RM measurements of the Coma-galaxy cluster are of the highest quality as they can be obtained from measurements of the magnetic field from seven individual radio galaxies located in the central part of the cluster (see Bonafede et al. 2010, for the details). Additionally, the magnetic field of the Coma cluster could be further constrained over seven more radio galaxies detected around the location of the infalling galaxy group NGC 4839 (Bonafede et al. 2013). From these studies, one can derive a value on the order of a few microgauss for the magnetic field in the Coma cluster. However, the picture may change for cool-core clusters where higher magnetic fields of the order of a few tens of microgauss are observed (see van Weeren et al. 2019, for a more detailed review of the observed magnetic field strengths in different galaxy clusters). Specifically, Vogt & Enßlin (2003) point to a cluster magnetic field of around 12 μG in the cluster Hydra as a conservative estimate.

One can place a lower limit on the magnetic field strength in galaxy clusters, which is given by the absence of inverse Compton emission that should originate from photons that are scattered on cosmic-ray electrons in the cluster environment. This should infer a spectrum of hard X-ray emission, which could clearly be distinguished from the background thermal Bremsstrahlung spectra (e.g., Rephaeli 1979; Rephaeli et al. 1994; Sarazin & Kempner 2000). However, modeling the inverse Compton emission is not straightforward (see van Weeren et al. 2019, and references therein). Despite, the difficulty in the modeling the hard X-ray spectrum of inverse Compton radiation, one can derive an upper limit from which one can infer a lower limit on the magnetic field strength in galaxy clusters ranging from 0.1–0.5 μG in the Coma cluster (e.g., Rossetti & Molendi 2004), the Bullet-Cluster (e.g., Wik et al. 2014) and Abell 2163 (e.g., Sugawara et al. 2009; Ota et al. 2014).

While there seems to be some observational consensus on the magnetic field strengths in galaxy clusters, the situation for the origin of these magnetic fields is less clear. The underlying problem is that the magnetic field amplification in the ICM is supposedly driven by the turbulence that is injected by shocks during the structure formation process (e.g., Miniati et al. 2001; Iapichino & Brüggen 2012; Iapichino et al. 2013, 2017). The involvement of ICM MHD turbulence makes this problem particularly difficult to control, as turbulence is not well understood in numerical simulations to begin with and is thus even less as the driver of magnetic field amplification (e.g., Donnert et al. 2018, for a very detailed discussion of this problem).

Despite the somewhat tedious process of understanding turbulence in numerical simulations one can draw a two-sided picture of the magnetization of the ICM. As already pointed out in the classical picture of turbulent amplification of magnetic fields in the ICM, the magnetic field could be amplified during the structure formation process when the largest structures (galaxy clusters) assemble in the universe. In this framework, the idea is that the magnetic field is amplified due to turbulence in the ICM driven by cosmic accretion that drives strong shocks in the ICM. While previous research indicates the presence of very strong shocks in the ICM, they are very rare with a volume filling factor of less than 0.1% for shocks with a Mach number larger than 10 (e.g., Miniati et al. 2001; Vazza et al. 2011). This renders the ICM (on average) as a subsonic medium, where the turbulent kinetic energy is orders of magnitude smaller than the thermal energy. However, we note that the total kinetic energy (bulk motion plus turbulence) is typically in equipartition with the thermal component and the galaxies can easily move at supersonic speed with respect to the ICM. However, the argument is about the volume filling factor of sub versus supersonic turbulence and since high Mach number shocks are extremely rare (even in the outskirts of the ICM) the ICM can be considered as formally subsonic. Hence, magnetic field amplification in the ICM can be described beautifully by the dynamo theory first developed by Kazantsev and Kraichnan (Kazantsev 1968; Kraichnan & Nagarajan 1967), which has been further advanced by several authors since then (e.g., Boldyrev & Cattaneo 2004; Kazantsev et al. 1985; Kulsrud & Anderson 1992; Kulsrud et al. 1997; Ruzmaikin et al. 1988; Subramanian & Barrow 2002; Xu & Lazarian 2020; Zel’dovich 1965, 1970, 1983).

This theory self-consistently describes the amplification of magnetic fields via subsonic turbulence due to stretching, twisting, and subsequent folding of magnetic field lines and can be tested through means of the magnetic power spectrum and the distribution of the curvature of magnetic field lines derived from high-resolution numerical simulations (e.g., Schekochihin et al. 2004; Porter et al. 2015). This theory has been widely applied on scales of galaxies (e.g., Wang & Abel 2009; Kotarba et al. 2009; Beck et al. 2012; Pakmor & Springel 2013; Martin-Alvarez et al. 2020; Marinacci et al. 2015; Marinacci & Vogelsberger 2016; Pakmor et al. 2017, 2020; Rieder & Teyssier 2016, 2017a, 2017b; Steinwandel et al. 2019, 2020a, 2020b; Tricco et al. 2016). While the concept of turbulent magnetic field amplification works well on galaxy scales, the ideal application for this theoretical framework is magnetic field amplification in the ICM due to the subsonic nature of ICM turbulence.

Already in early numerical simulations of magnetic fields in galaxy clusters it has been pointed out that the amplification via the small-scale turbulent dynamo is quite likely without showing direct evidence of the process (Dolag et al. 1999, 2001, 2002, 2005; Brüggen et al. 2005; Dubois & Teyssier 2008; Ryu et al. 2008; Vazza et al. 2014) and recently the first efforts have been made to show direct evidence of an acting small-scale-turbulent dynamo on the scales of galaxy clusters in Eulerian codes (e.g., Vazza et al. 2018; Roh et al. 2019) by directly comparing to the statistics, which is enforced by dynamo theory and high-resolution numerical simulations of the small-scale turbulent dynamo in the ICM regime (e.g., Schekochihin et al. 2004; Porter et al. 2015).

Generally, the idea of the small-scale turbulent dynamo is that the magnetic field is generated on the scales of the turbulent eddies and is thus of scale-free nature. That means the stretching, twisting, and folding of field lines can occur on parsec scales in the ISM or on megaparsec scales in the ICM with a growth rate that is proportional to the eddy turnover time with respect to the ambient medium. The field is then propagated to larger scales via an inverse turbulence cascade when the magnetic energy density reaches equipartition with the turbulent energy density stored in the smallest eddies. This process leads to an increase of the power stored in the magnetic field on larger scales with a subsequent decrease on the largest scales that is predicted by the evolution of the energy spectra in Kazantsev–Kraichnan theory.

In this paper, we study the buildup of the magnetic field in numerical simulations based on the theory of the small-scale turbulent dynamo in the framework of a Lagrangian simulation framework that is currently missing in the literature on magnetic field amplification in galaxy clusters.

This paper is structured as follows. In Section 2, we discuss the theoretical background that is needed to understand the basics of the theory of small-scale turbulent dynamo. In Section 3, we discuss the basics of the numerical algorithms we use to carry out our simulations, including the handling of nonideal MHD. In Section 4, we discuss the details of the simulation suite along with the initial conditions and the adopted naming conventions. In Section 5, we present the results of the simulations and carry out the analysis that is needed to study the turbulent dynamo in the ICM. In Section 6, we summarize our findings, conclude our results, and comment on model limitations and future work. Furthermore, we discuss model variations in the Appendix.

2. Short Overview of Small-scale Turbulent Dynamo Theory

Magnetic field amplification in the ICM is most likely driven by turbulence that is injected during structure formation. The idea is that tiny magnetic fields of which the origin is still under debate (e.g., Biermann 1950; Demozzi et al. 2009; Matarrese et al. 2005; Gnedin et al. 2000; Rees 1987, 1994, 2005, 2006) are amplified to large-scale coherent fields via stretching, twisting, and folding of magnetic field lines. This process is limited by the magnetic tension force that makes every stretch-twist-fold process inherently more difficult and finally saturates the dynamo when the turbulent kinetic energy in the smallest eddies is in equipartition with the magnetic field energy. Over the past 60 yr several authors have continuously refined the theoretical understanding of turbulent magnetic field amplification (e.g., Boldyrev & Cattaneo 2004; Kraichnan & Nagarajan 1967; Kazantsev 1968; Kazantsev et al. 1985; Kulsrud & Anderson 1992; Kulsrud et al. 1997; Ruzmaikin et al. 1988; Subramanian & Barrow 2002; Zel’dovich 1965, 1970, 1983; Xu & Lazarian 2020, to just name a few). The picture of turbulent dynamo amplification fits perfectly to the use case of magnetic field amplification in galaxy clusters as the ICM is very hot and therefore is characterized by the subsonic turbulent energy cascade that is in very good agreement with the Kolmogorov theory of turbulence (Kolmogorov 1941).

Every theory of magnetic field amplification starts with the induction equation in the continuum limit of MHD:

∂B∂t=∇×(∇×B)+ηΔB,(1)

where η is the magnetic diffusivity. A quite intuitive approach toward an understanding of the magnetic field structure developed by the small-scale turbulent dynamo can be derived by statistically studying fluctuations in the velocity and magnetic field. A general vector field can always be Fourier decomposed. Thus, this means the velocity field can be written as v = ∫d³kkvk × eik·r and the magnetic field can be written as B = ∫d³kkBk × eik·r. The ultimate goal of the Fourier analysis of these fluctuations is to derive the distribution of the power in the magnetic field that is given via

Emag〈B2〉8π=∫PM(k)dk.(2)

One can find the time derivative of PM(k) (see Kulsrud & Zweibel 2008, for the details of the derivation):

∂PM(k)∂t=∫K(k,k0)M(k0)dk0−2βk2PM(k),(3)

which describes the time evolution of the magnetic power spectrum as a function of structure function K and turbulent resistivity β. The combination of Equation (2) and Equation (3) yields

dEmagdt=2γEmag,(4)

where γ denotes the growth rate of the dynamo. From this one can straightforwardly see that the magnetic field strength is increased by a factor of 2 as a function of the eddy turnover time. This is consistent with the stretching, twisting, and folding as it is assumed to occur on the small-scale turbulent dynamo. Furthermore, it is worth noting that in this prescription the growth rate is then directly correlated with the eddy turnover rate of the smallest eddies and energy is carried to larger scales by an inverse turbulent cascade. In the kinematic regime, one can find PM(k) by solving

∂PM(k)∂k=γ5k2∂2PM(k)∂k2−2k∂PM(k)∂k+6PM(k)−2k2λresPM(k),(5)

where λres is the resistivity. This differential equation can be solved with standard methods and one obtains

PM(k,t)∝e3/4γtk3/2,(6)

which directly indicates the exponential growth of modes with k3/2. Therefore, the small-scale turbulent dynamo can be clearly identified over the shape of the magnetic energy spectra. However, in practice power spectra can be tricky to measure in Lagrangian codes and we suggest also considering the imprint of the dynamo on the magnetic tension force as a quantitative prior to measuring dynamo action in numerical simulations as previously suggested by various groups on galaxy and ICM scales (e.g., Schekochihin et al. 2004; Vazza et al. 2018; Steinwandel et al. 2019).

3. Numerical Method

We carry out the simulations presented in this paper with the Tree-SPMHD-Code p-gadget3, which is the developers’ version of the Tree-SPH-Code p-gadget2 (Springel 2005). We use a modern implementation of SPH, as presented in Beck et al. (2016) that includes time-dependent artificial viscosity and conduction and employs higher-order kernel functions, given as the well-studied Wendland functions (Wendland 1995, 2004; Dehnen & Aly 2012). However, as thermal conduction is very important in the ICM we do not use the time-dependent artificial conduction implementation within the simulations but rather use the physical conduction implementation first presented in Jubelgas et al. (2004) and later updated by Arth et al. (2014) to a conjugate gradient solver for improved convergence and stability of the scheme. The conjugate gradient solver employed is similar to the one developed in Petkova & Springel (2009) for use in galaxy formation simulations that include direct radiative transfer. We run all the simulations with physical conduction and 1/20 of the canonical Spitzer value. This is motivated by a number of research works that have adopted this value based on high-resolution simulations of plasma instabilities, which in turn can suppress the heat flux for thermal conduction (e.g., McCourt et al. 2012; Berlok et al. 2021). However, we note that the exact suppression factor of the Spitzer value is still heavily debated in the literature. Nevertheless, we note that our value of 1/20 is well within the limits of previous studies such as in Kannan et al. (2017). Our version of p-gadget3 further includes magnetic fields and magnetic dissipation as presented in Dolag & Stasyszyn (2009). While we will use isotropic thermal conduction in our default runs in the main paper we carry out some additional runs utilizing anisotropic thermal conduction following the prescription in Arth et al. (2014). We note that the usage of anisotropic thermal conduction increases the computational cost of a simulation by roughly 20%. In the following, we briefly discuss the specifics of the underlying SPH equations.

3.1. Kernel Function and Density Estimate

The current SPH scheme implemented in our code is based on the density-entropy formulation of SPH, which means that we smooth the density field in the following fashion:

ρi=∑jmjWij(xij,hi),(7)

where hi is the smoothing length. The summation carried out in Equation (7) is computed over the neighboring particles within the kernel Wij(xij, hi):

Wij(xij,hi)=1hi3w(q).(8)

In our simulations, the kernel function is used with 295 neighboring particles. The function w(q) is given by

w(q)=136564π1−q61+8q+25q2+32q3,(9)

for q < 1. For q > 1, we set w(q) to zero and it is known as the Wendland C6 kernel function (Wendland 1995, 2004; Dehnen & Aly 2012).

3.2. Equation of Motion in SPH and SPMHD

The equation of motion (EOM) for SPH can conveniently be derived from a discrete Lagrangian as presented in Price (2012) by using the physical principle of least action. This has the distinct advantage that the derived formulation is conserving energy, momentum, and angular momentum by construction. This leads to the SPH formulation of the EOM in the pure hydrodynamic case:

dvidt=−∑jmjficoPjρj2∂Wij(hi)∂ri+fjcoPjρj2∂Wij(hj)∂ri.(10)

with fjco given by

fjco=1+hj3ρj∂ρj∂hj−1.(11)

A similar argument can be made for the SPMHD case leading to the SPH formulation of the MHD EOM given as

dvidt=−∑jmjficoPi+12μ0Bi2ρi2∇iWij(hi)+fjcoPj+12μ0Bj2ρj2∇iWij(hj)+1μ0∑jmjficoBi[Bi·∇iWij(hi)]ρi2+fjcoBj[Bj·∇iWij(hj)]ρj2.(12)

The presence of the magnetic field alters the EOM in several ways. First, the presence of the magnetic field leads to an additional pressure component, apart from the thermal pressure within the fluid. This additional pressure component scales as B². The fact that the magnetic pressure component scales as B² is crucial to establish pressure equipartition with the thermal pressure relatively quickly. Second, the term on the left-hand side of Equation (12) arises due to the divergence cleaning constraint ∇ · B = 0. This term is problematic because it breaks the symmetry of the underlying Lagrangian in a way that the system is no longer invariant under the rotation of the system. Thus, from Noether’s theorem one can easily see that the SPMHD equations do not strictly conserve angular momentum.

3.3. Formulation of the Induction Equation in Nonideal MHD with Effective ηm

Furthermore, how the magnetic field is influenced by the EOM it is not only interesting, but how the magnetic field itself evolves with time is also interesting. The evolution of the magnetic field is generally given by the induction equation that takes the form:

∂B∂t=∇×v×B+∇×ηm(∇×B),(13)

which can be reformulated as

dBdt=(B·∇)v−B(∇·v)+ηmΔB.(14)

To this day, most MHD simulations of galaxies or galaxy clusters drop the last term of Equation (14) and very few simulations include these terms (e.g., Kotarba et al. 2011; Bonafede et al. 2011; Steinwandel et al. 2019, 2020a). However, this term is crucial for modeling the plasma in an accurate fashion⁵ and thus we include it in our cosmological galaxy-cluster simulations. The parameter ηm is hereby an effective diffusion parameter that is comprised of the contribution due to thermal conduction by ηCoulomb that is related to the thermal conductivity σ and the turbulent diffusion coefficient of the plasma ηturb. While the exact value of the diffusion coefficient in the ICM is under debate and several processes yield different limits (e.g., Strong et al. 2007; Lesch & Hanasz 2003; Schlickeiser et al. 1987; Schuecker et al. 2004; Maier et al. 2009; Rebusco et al. 2006) we use a moderate value of ∼ ηturb ∼ 2 × 10²⁷ cm² s⁻¹. This is in good agreement with the classical Spitzer model for the ICM (Spitzer 1956). For a more detailed discussion on the choice of ηm we refer the reader to Section 4.2 of Bonafede et al. (2011).

Following Dolag & Stasyszyn (2009) the diffusion term in the induction equation takes the form

∂B∂t=ηmρiHa2∑j=mjρj(Bi−Bj)(ri−rj)∣(ri−rj)∣·∇iWi,(15)

where W is the kernel and Bi and Bj are the magnetic field vectors at the positions ri and rj. Ha² is the internal scaling applied in gadget needed for correct unit conversion from the comoving field B to the physical field B/a². However, as the magnetic field is dissipated this introduces an additional entropy term for the thermal plasma, manifesting in the rate of change of entropy:

dAdt=−γ−12μ0ρiγ−1∑jmjρjBi−Bj2(ri−rj)∣(ri−rj)∣·∇iWi.(16)

Further implementation details can be found in Dolag & Stasyszyn (2009).

3.4. Thermal Conduction

Thermal conduction is believed to be a process of importance in the ICM. We will briefly describe the physical process as well as the numerical implementation of the process into our simulation code for the isotropic and the anisotropic cases. We note that while we run most of the simulations with isotropic conduction in the presence of magnetic fields we run one simulation with anisotropic conduction and show the results in Appendix C. We note that at our adopted 5% of the Spitzer value the difference between the isotropic case is only marginal and we adopted isotropic thermal conduction to save computational cost. We adopt thermal conduction in our simulations, since we attempt a bottom-up approach for modeling the ICM where we first want to study the plasma astrophysical aspects of the ICM, including magnetic fields and thermal anisotropic conduction (in this paper), as well as anisotropic viscosity and cosmic rays (in follow-up work), before we switch to more complex physical models for cooling, star formation, and active galactic nuclei (AGN) feedback. This approach is driven by first understanding first the plasma astrophysics of the ICM before complicating the picture with ad hoc subgrid models to study galaxy formation physics.

3.4.1. Isotropic Case

Thermal conduction is the physical process that describes heat transfer via scattering of free electrons. Thus, in order to properly work one needs a high ionization fraction of the underlying plasma. In the ICM this is straightforwardly achieved as the ICM has virial temperatures that can easily reach 10⁸ K. One can follow Spitzer (1956) to get the heat flux as a function of the gradient of the temperature distribution via

Q=−κ∇T,(17)

where κ is the conduction coefficient. In the classical Spitzer case, one makes the assumption of an idealized Lorentzian gas for which one can find the canonical Spizter value given by

κSp=202π3/2kBTe5/2me1/2e4ZlnΛ,(18)

where kB is the Boltzmann constant, me is the electron mass, Te is the electron temperature, e is the elementary charge, Z is the average number of protons in the plasma, and ln Λ is the Coulomb logarithm. As the temperature of clusters is very high the conduction coefficient shows a strong dependence on the temperature one can infer a strong contribution of heat conduction to the dynamical processes driven in the ICM. In a realistic plasma one typically does not reach the Spitzer regime and conduction is suppressed. This is strongly dependent on the average number of protons present in the plasma (see Spitzer & Härm 1953), which results in a suppression factor of 0.3 (Spitzer 1956) under the assumption of a primordial distribution of the gas, which is a good first-order assumption for cosmological simulations. However, one often applies the following parameterization of κ in cosmological simulations:

κ=1.31·neλekBkBTeme1/2,(19)

which can be refactored to obtain

κ=4.6·1013Te108K5/240lnΛergscmK,(20)

where ne is the electron number density and λe is their mean free path. From this one can directly see why the heat transport is dominated by electrons or why the heat transport by protons in the plasma is subdominant, as the conductivity scales with the inverse mass. Thus, heat conduction is dominated by the electron population of the plasma.⁶ We therefore neglect any contribution of protons and furthermore make the assumption that the Coulomb logarithm is constant with lnΛ=37.8. This picture is incomplete as for now we have assumed that the length scale that is given via

lT=T∣∇T∣(21)

is much larger than the mean free path of the electrons. Strictly speaking one can only make this assumption for a higher density plasma. However, in a lower density plasma such as the ICM where the temperature gradient and the mean free path of the electrons are roughly of the same order, transporting energy by conduction is limited by the low number of interaction rates in the plasma. Thus, we are in the conduction-saturation limit of Cowie & McKee (1977), who computed the limited heat flux for a low-density plasma

Qsat=0.4nekBT2kBTπme1/2.(22)

Now, one can interpolate between Equation (21) and Equation (22) to obtain the total heat flux

Qtot=−κTlT+4.2λ∇T∣∇T∣.(23)

This is equivalent to a renormalized conduction coefficient

κ=κSp1+4.2λ/lT.(24)

Therefore, one can finally formulate the rate of change of the energy per unit mass:

dudt=−1ρ∇·Q=1ρ∇·(κ∇T).(25)

3.4.2. Anisotropic Case

In a magnetized plasma heat conduction is slightly more complicated because as the electrons are charged particles their scattering processes are influenced by the magnetic field structure. While electrons can move freely alongside magnetic field lines their motion perpendicular to the field lines is suppressed; hence, the term anisotropic conduction. The trajectory of electrons in the presence of magnetic fields is well studied and they gyrate around magnetic field lines with the Larmor frequency

ωg=eBmec,(26)

where c is the speed of light. This affects the movement of the electrons in the presence of magnetic fields as pointed out by Frank-Kamenetskii (1967). Following Braginskii (1965) one can subdivide the heat flux into three additive terms

Q=−κ∣∣∇∣∣T−κ⊥∇⊥T−κΛB×∇T,(27)

where the first two are referred to as the parallel and perpendicular components and the last term is called the “Hall term.” We drop this term as we will see that it will vanish once we start the discretization of our numerical scheme. However, the remaining contributions are extremely tedious to describe. To do so it is often useful to introduce the so-called diffusion coefficients D that are related to the conduction coefficient κ by κ∼DnekB. One can now distinguish between two cases for the perpendicular diffusion coefficient: (i) a diffusion coefficient that scales as B−2 and (ii) a diffusion coefficient that scales with B⁻².

In the former case, one can assume that the diffusion coefficient D can be described as D ≈ v²τ. Since electrons can move freely along the field lines we get D∣∣ = D. On the other hand, perpendicular to the field lines electrons can only move by interchanging cyclotron frequencies so that D⊥≈λv/(ωg2τ2) and one can straightforwardly determine that D⊥/D∣∣ ∝ B⁻² for ωgτ ≪ 1. If the gyro-radius is of the order of the mean free path, one finds the relation D⊥/D∣∣≈1/(1+ωg2τ). For typical values of the ICM, one obtains D⊥/D∣∣ ≈ 10⁻²⁸.

However, in practice, it is a bit more complicated because the transport process perpendicular to the field lines interacts with turbulent diffusion processes and possibly reconnection diffusion events, which make the interaction highly nonlinear, and experiments conducted in the laboratory indicate that the scaling is rather of order B⁻¹ than B⁻² given by Bohm diffusion (Guthrie et al. 1949).

Finally, we can write the heat flux in the case of anisotropic thermal conduction as

Q=−κ∇T·B∣∇T·B∣∇T=−κcos(θ)∇T.(28)

In practice, we split the conduction equation into two parts, which results in the rate of change of the energy per unit mass of the form

dudt=1ρ∇·κ∣∣Bˆ·∇TBˆ+κ⊥∇T−Bˆ∇TBˆ.(29)

3.4.3. Numerical Implementation

Now we need to discretize Equation (29) to obtain its SPH formulation. First, we rewrite Equation (29) by factoring (Bˆ·∇T)Bˆ, which yields

dudt=1ρ∇·(κ∣∣−κ⊥)Bˆ·∇TBˆ+κ⊥∇T.(30)

A straightforward way to solve this equation is to do an operator split and solve for the divergence and the temperature gradient in chained SPH loops. In practice, this has the disadvantage that it increases the computational cost and leads to increased noise in the solution as each loop will independently add partition noise to the final result. Thus, this is not the favored way of solving this problem. We will follow the methodology derived in Petkova & Springel (2009), who discretized a similar diffusion equation in the context of radiative transfer. In the following, we will derive the discrete form of the anisotropic conduction equation but only for the first term. The reason for this is that the second term in Equation (30) can be solved as described in Jubelgas et al. (2004) or Steinwandel et al. (2020b) and thus we refer the reader to those papers for a review of how to solve the second term in SPH. The first term takes the discrete form of

dudt1st=1ρ∑α,β∂∂xα(κ∣∣−κ⊥)BˆαBˆβ∂∂xβT.(31)

We note that α and β represent the components of a tensor of second order. We substitute the tensor components by Aαβ=(κ∣∣−κ⊥)BˆαBˆβ, which yields

dudt1st=12ρ∑α,β∂2AαβT∂xα∂xβ−T∂2Aαβ∂xα∂xβ+Aαβ∂2T∂xα∂xβ.(32)

Now we want to rewrite the second derivative present in Equation (32). This can be achieved by using the following identity for an arbitrary vector Y (e.g., Price 2012):

∑α∂2Yi∂xα2=2∫d3xj(Yj−Yi)xijT·∇iWij∣xij∣2.(33)

We can use Equation (33) to rewrite the right-hand side of Equation (32) to obtain

dudt1st=1ρi∫d3xixjT(Aj+Ai)(Tj−Ti)∣xij∣2∇iWij.(34)

Finally, we can write the integral on the right-hand side of Equation (34) by the SPH sum over the neighbors to get

dudt1st=μ(γ−1)kBρi∑jNngbmjρj·xijT(Aj+Ai)(uj−ui)∣xij∣2∇iWij,(35)

where μ is the molecular weight and γ is the adiabatic index, which we set to 5/3 in the whole simulation domain. We note that we essentially recover the same scaling of the anisotropic case with the temperature difference between SPH particles as we know it from Jubelgas et al. (2004) and Steinwandel et al. (2020b). The problem with Equation (35) is that there is a precondition for actually solving this for the tensor Ai + Aj, which needs to be positive definite to establish physical heat flux from hot to cold. This is related to the fact that A scales with the difference between κ∣∣ and κ⊥, which can be negative for strong anisotropies of the heat flux. Thus, technically the heat flux can be negative and flow from cold to hot, which is nonphysically. There are several methods to circumvent this. First, one could essentially take the approach of flux-limited diffusion. Second, only do the anisotropic part if the tensor is positive definite. Third, isotropize the tensor. We follow Petkova & Springel (2009) and do the latter by adding an artificial isotropic component to the heat flux tensor, which yields A→αA+1/3(1−α)tr(A)1. We follow Petkova & Springel (2009) and adopt α = 2/3. We note that we actually do not do a discretization of the Hall term since in our parameterization it is simply vanishing from the discretized equations.

Finally, we solve the actual differential equation. We do this by adopting the method from Petkova & Springel (2009) with the so-called conjugate gradient method, which requires an additional SPH loop but shows very accurate results in practice. The conjugate gradient method is in principle a numerical method used to solve a matrix inversion problem of the form

C·x=b.(36)

The method is implicit and iterates to a solution with a very high convergence order. However, it is numerically unstable, which means that it does not guarantee a solution for arbitrary particle distributions. However, the underlying idea of a conjugate gradient method is the following. The algorithm is supposed to monotonically approach the solution of each element of the inverted matrix by adopting a weight that is dependent on the residual of the previous iterations based on a good initial guess of the state of the physical system. In order for this to work, C needs to be real, positive definite. However, as we already use a correction to get an isotropic version of our heat flux tensor we can assume that both of these assumptions are valid. The final task that remains is to write Equation (35) in the form of a conjugate gradient,

uin+1=uin+∑jNngb,cij(uin+1−ujn+1),(37)

where we adopted dui/dt=Δu/Δt=(uin+1−uin)/Δt. The cij are computed as follows:

cij=−(γ−1)μkBmjΔtρiρjxijT∣xij∣2Ai+Aj∇iWij.(38)

This is the conjugate gradient version of our anisotropic thermal conduction equation (the first term). We note that the scheme can be numerically unstable. However, this can be resolved by switching to a biconjugate method (bicgstab), which applies convergence through a second direction (geometrically speaking).

We note that there are a lot of similar implementations that treat different parts of thermal conduction mostly in the realm of galaxy-cluster formation (Ruderman et al. 2000; Dolag et al. 2004; Schekochihin et al. 2008; Rasera & Chandran 2008; Parrish et al. 2009; Sharma et al. 2010; ZuHone et al. 2013; Suzuki et al. 2013; Komarov et al. 2014; ZuHone et al. 2015; Dubois & Commerćon 2016; Kannan et al. 2017; Yang & Reynolds 2016).

4. Initial Conditions and Simulations

We run a suite of six non-radiative (without cooling and star formation) cosmological zoom-in simulations of a galaxy cluster with a target mass of M₂₀₀ ∼ 2 × 10¹⁵ M⊙ to the target redshift z = 0. We run the simulations at three different mass resolutions, which we will refer to as 1×, 10×, and 25×. In this naming convention, the leading number indicates 1, 10, or 25 times the mass resolution that is achieved in the Magneticum high-resolution cosmological volume simulations (Hirschmann et al. 2014). We choose this setup for the following reasons. First, this allows us a detailed study of the plasma physics that is acting within the ICM without being polluted by the additional energy input of supernovae (SNe) or AGN. This enables us to investigate the magnetic field amplification in the ICM without the need to retune feedback parameters on the different resolution levels and directly sets the origin of the turbulence injected in the system to be driven by gravitational forces only. Thus, we can link any amplification of the magnetic field directly to the gravo-turbulence within the ICM without considering turbulence driven by stellar or AGN feedback. However, we will comment on the caveats of this in Section 6 in greater detail. Second, previous studies of the dynamo in the ICM that were obtained with the grid code enzo presented in Vazza et al. (2018) have been carried out with a similar non-radiative setup, which allows for a pristine comparison with their results.

4.1. Simulation Setup

The initial conditions for the cluster at hand are selected from a lower resolution volume with a box size of 1.0 Gpc and a base resolution of 1024³ dark matter particles leading to a particle resolution of around 10⁸ M⊙. We use a Wilkinson Microwave Anisotropy Probe 7 cosmology with Ω₀ = 0.24, ΩΛ = 0.76, Ωbaryon = 0.04, h = 0.72, and σ₈ = 0.8. We select the dark matter particles at z = 0 in cosmological volume at base resolution and trace them back to a resolution-dependent initial redshift using the code zic (Tormen et al. 1997). The code allows for arbitrary shapes of the high-resolution regions to avoid overhead by oversampling the high-resolution regions when it is simply tied to a sphere or an ellipsoid. The initial redshifts for the three resolution levels are as follows: zini = 70 (1×), 140 (10×), and 180 (25×).

The region that is selected for re-simulation is chosen to be large enough to avoid pollution of the target halo at z = 0 by lower resolution intruder particles that originate from the lower resolution large-scale structure in the gigaparsec volume. For each resolution level, one dark-matter-only test run has been carried out to ensure the quality of the initial conditions in which the gas particles have been coevolved as separate dark matter species. As the mass resolution and the applied force softening change on all of the three resolution levels, we sum the basic simulation parameters in Table 1. We note that the force softening in our lowest resolution run is 3.0 kpc which is still smaller than the spatial resolution of the highest resolution run of Vazza et al. (2018). In our highest resolution run, the force softening is pushed to 1 kpc.

Table 1
Overview of the Number of Particles, Mass Resolution, Gravitational Softening Lengths, and Physical Field Strength of the Three Targeted Resolution Levels

Particle Numbers
|  | 1× | 10× | 25×
Gas particles | Ng | ∼ 1 × 10⁶ | ∼ 1 × 10⁷ | ∼ 3 × 10⁷
Dark matter | Ndm | ∼ 1 × 10⁷ | ∼ 1 × 10⁷ | ∼ 3 × 10⁷
Mass resolution [M⊙]
Gas particles | mgas | ∼ 1.4 × 10⁸ | ∼ 1.4 × 10⁷ | ∼ 5.6 × 10⁶
Dark matter | mdm | ∼ 6.9 × 10⁸ | ∼ 6.9 × 10⁷ | ∼ 2.7 × 10⁷
Gravitational softening [kpc]
Gas particles | ϵgas | 3.0 | 1.4 | 1.0
Dark matter | ϵdm | 3.0 | 1.4 | 1.0
Seed field (physical)
|  | 5.02 × 10⁻¹¹ | 1.98 × 10⁻¹⁰ | 3.28 × 10⁻¹⁰

As we carry out all simulations in a non-radiative fashion we need an initial seed field as this is the only possible origin of the magnetic field in a scenario like this. For this, we choose a default value of Bseed = 10⁻¹⁴ G (comoving). As we start our simulations at different redshifts, this corresponds to a variation in the physical seed field of Binit,ph=Bseed·(1+zinit)2 for the different resolutions, which we sum in Table 1. We note that this is a rather conservative choice for the magnetic seed field and other simulation groups often take larger values (e.g., Vazza et al. 2018).

Furthermore, we will test different physical settings on our 10× model. For this purpose, we need to introduce naming conventions and model variation parameters.

4.2. Naming Conventions and Model Variations

The MHD simulations with the abovementioned default settings are referred to as 1×, 10×, and 25×. Throughout the paper, we will make quite limited use of our hydrodynamics-only models. Therefore, whenever we refer to one of the hydrodynamics-only simulations we will explicitly label it as the hydrodynamic realization of the model 1×, 10×, or 25×. The focus of this work is magnetic field amplification in the ICM and the hydrodynamics-only simulations are merely reference simulations used for the comparison with the MHD runs.⁷

Furthermore, we will test three important model variations on the 10× run and introduce new naming conventions for the specifics of the respective run. The three model variations are as follows:

First, we will investigate the variation of the magnetic seed field to a 10 times smaller and a 10 times higher value than the value we chose in our default setting (for 10× this is 1.98 × 10⁻¹⁰ G.), resulting in 1.98 × 10⁻¹¹ and 1.98 × 10⁻⁹ G. These two simulations are labeled as 10×-low-seed and 10×-high-seed and will be the subject of Appendix A.

Second, we will investigate how robust our results are on the choice of the numerical diffusion parameter. Our default setting for the diffusion parameter is ∼ ηturb =2 × 10²⁷ cm² s⁻¹ (Bonafede et al. 2011). In order to understand the dependence of the magnetic field growth on the diffusion parameter, we will also vary this to a 10 times lower value and a 10 times higher value of ηturb = 2 × 1026 and ≈2 × 10²⁸ cm² s⁻¹, respectively. For these simulations, we introduce the naming conventions 10 × low-eta and 10×-high-eta, which are the subject of Appendix B.

Third, we will carry out one additional run with anisotropic thermal conduction where we directly include the magnetic field structure in our thermal conduction solver via a biconjugate gradient solver (Arth et al. 2014). All other simulations are carried out with physical, but isotropic conduction, which is a potential caveat in the MHD case. However, isotropic conduction is already computationally expensive to solve and takes up roughly 15% of the computing time. Anisotropic conduction is even more demanding in terms of computational cost and the memory imprint of the code. Thus, we only carry out one simulation labeled as 10×-ani with the effect of an anisotropic physical conduction, which will be the subject of Appendix C.

We show an overview of all the simulations with the physics variations that we carried out in this work in Table 2.

Table 2
Overview of the Physics Variations Adopted throughout Our Different Simulations

Name | Nonideal MHD | η [cm² s⁻¹] | Thermal Conduction | Anisotropic Thermal Conduction | κ/κs
1× | ✓ | 2 × 10²⁷ | ✓ | × | 0.05
10× | ✓ | 2 × 10²⁷ | ✓ | × | 0.05
25× | ✓ | 2 × 10²⁷ | ✓ | × | 0.05
1×-NO | × | ⋯ | ✓ | × | 0.05
10×-NO | × | ⋯ | ✓ | × | 0.05
25×-NO | × | ⋯ | ✓ | × | 0.05
10×-low-seed | ✓ | 2 × 10²⁷ | ✓ | × | 0.05
10×-high-seed | ✓ | 2 × 10²⁷ | ✓ | × | 0.05
10×-low-eta | ✓ | 2 × 10²⁶ | ✓ | × | 0.05
10×-high-eta | ✓ | 2 × 10²⁸ | ✓ | × | 0.05
10×-ani | ✓ | 2 × 10²⁷ | ✓ | ✓ | 0.05

5. Results

In this section, we present the results of our simulations. First, we discuss the cosmological assembly of the structure in terms of halo mass and will investigate the general impact of the magnetic field on the structure formation process. This is followed by a detailed investigation of the buildup of the magnetic field from the initial redshift down to redshift z = 0. Finally, we will briefly discuss the impact of the divergence cleaning constraint on the structure of the ICM in our simulated galaxy cluster.

5.1. Cosmological Assembly and General Cluster Properties

First, we briefly discuss the cosmological assembly of our quintessential galaxy cluster to gauge that it reaches a halo mass of around M₂₀₀ ∼ 2 × 10¹⁵ M⊙ at redshift z = 0. We show this in Figure 1 for our 1× resolution simulations with and without magnetic field from redshift z = 4–0. The structure itself starts to form at a much higher redshift and has already assembled around ∼2 × 10¹³ M⊙ by redshift z = 4, which is around 1% of the mass that it will acquire by redshift z = 0. The cluster undergoes very rapid growth between redshift z = 4 and 1.5 from ∼2 × 10¹³ to ∼3 × 10¹⁴ M⊙, which corresponds to a growth rate of around 9.5 × 10¹³ M⊙ Gyr⁻¹. After that, the system transits into a phase of weaker growth until redshift z = 0.8 in which it doubles its mass. This is followed by a major merger at redshift z = 0.8 at which the system roughly acquires another 30% of its total mass up to that point pushing it just below the 10¹⁵ M⊙ mark and then finally transits to the regime of the continued growth of the system via smooth accretion. This is followed by another major merger at redshift z = 0.3. Past redshift z = 0.2, the system quietly assembles the rest of its mass until it reaches a final mass of ∼ 2 × 10¹⁵ M⊙ at redshift z = 0. We note that the reference run without the effects of the magnetic field (red line in Figure 1) follows the MHD run very closely with an error below the 1% margin for most of the evolution of the system. This gauges the expected very weak effect of the presence of the magnetic field on the large-scale assembly of the structure and shows that the magnetic field does not alter the behavior of the cosmological assembly of the structure. However, there is one exception to this, which is the slight delay in the first major merger of the system at around redshift z = 0.8, which slightly delays the merger and could potentially originate from the fact that the additional pressure component that is present as the magnetic field within the system slows down the collapse of the baryons into the dark matter halo, which indirectly slows down the assembly of the dark matter mass in the center of the halo.

Figure 1. We show the accretion history of the clusters for the MHD case (blue) and for the hydrodynamical (HD) case (red) for reference only for our lowest resolution run to gauge that magnetic fields have a weak effect on the structure formation process.

5.2. Morphology of the Cluster

We start the evaluation of our results by visualizing the key quantities of the cluster for our different runs in Figure 2. In the top row we show the model 1×, in the middle row we show the model 10×, and in the bottom row, we show the model 25×. The panels on the left show the gas surface density, the panels in the center show the temperature distribution of the cluster, and the panels on the right show the magnetic field strength in the three different models. The white dashed circle in the center of each panel indicates the virial radius of the cluster. In the 1× run we can see a clear lack of resolution, especially in the cluster outskirts beyond the virial radius of the system. This manifests as a vanishing substructure in the density distribution compared to the higher resolution models 10× and 25×. However, the largest difference between the 1× model and the 10× and 25× models can be seen in the temperature and magnetic field distributions. Visually it appears that there is more hot gas around 10⁸ K in the virial radius for the two higher resolution models 10× and 25×. Furthermore, we can identify a clear trend of an increase in the magnetic field strength within the virial radius by a factor of around 3 from the 1× model to the 10× and 25× models. In this context, we want to note that the particles with the maximum field strength within the simulation are located around the cluster center. For the 1× simulation, the particle with the maximum field strength has a value of ∼120 μG, for the 10× run we find ∼180 μG, and for the 25× run we find ∼240 μG. We note that there are very few particles on each resolution level that have similar magnetic field strength (it is around 10 particles for the 25× simulation, which have a field beyond 100 μG). We will discuss in Section 5.6 to which degree this behavior is driven by our nonzero divergence in all the simulations. Despite the slightly larger field strength in the runs 10× and 25× we want to point out the increase of magnetic field line structures that we can capture in the higher resolution runs 10× and 25× compared to the 1× run.

Figure 2. We show the projections of density (left), temperature (center), and magnetic field (right) for our runs of 1× (top), 10× (center), and 25× (bottom) at redshift z = 0. While the 10× and the 25× simulations show structural similarities, in the magnetic field and temperature structure, there are structural differences in the run 1×. Moreover, the central magnetic field increases by a factor of around 2.5 from 1–25×.

Moreover, we gauge the assembly of the cluster as a function of redshift for the run 25× for three different redshifts as shown in Figure 3. We show the same quantities as in Figure 2. The top row of Figure 3 shows the cluster at z = 2.331, the second row shows the cluster at a redshift of z = 1.180, and the third row is the same as the third row of Figure 2 showing the cluster at redshift z = 0. One can see that the central density of the cluster is continuously increasing as a function of redshift, which occurs due to subsequent merger and accretion events. We can see at redshift z = 2.331 that a massive structure is about to fall in from the top right. At z = 1.180, we can see smaller structures infalling from beyond the virial radius, which show extended tails of stripped gas. At redshift z = 0, the cluster evolves to a more relaxed state with a lower number of infalling objects. There is also a clear evolution in the temperature profiles, which we can see in the center panels from top to bottom where the cluster gas is strongly heated through its formation process down to redshift z = 0. The magnetic field structure is of particular interest as it is apparent from the evolution of the magnetic field on the right-hand side of Figure 3 from top to bottom that we can find a fully developed magnetic field with around a few microgauss already at redshift z = 2.331, which seems visually to decrease but occupies a larger volume as the system evolves toward redshift z = 1.180. At redshift z = 0, we find a fully developed field within the virial radius. Visually, the magnetic field amplification seems to be correlated with the turbulence that is injected via the structure formation process. We want to specifically point out that the field is stronger at around redshifts z = 2 and 1 compared to the field strength at redshift z = 0.

Figure 3. We show the projected density (left), temperature (center), and magnetic field (right) for three different points in time, at z ∼ 2 (top), z ∼ 1 (center), and z ∼ 0 (bottom) for our 25× run. We can clearly see that we have a fully developed magnetic field structure by redshift z ∼ 2.

5.3. Radial Evolution of the Cluster

Before we start the discussion on magnetic field amplification via the turbulent dynamo in the ICM we want to briefly report our results for radial profiles of central physical quantities at redshift z = 0. In Figure 4, we show the radial profiles out to a radius of 4 Mpc for the density (top left), the magnetic field (top right), the temperature (bottom left), and the pressure (bottom right) for the runs 1× (blue), 10× (red), and 25× (magenta). For all quantities, we find declining profiles as a function of the radius. As we are specifically interested in the magnetic field evolution of the cluster we note the most important findings regarding the radial trend of the magnetic field strength as a function of resolution. As the resolution increases from 1–10× and finally to 25× we find an increase in the central magnetic field from around 5 μG in the case of the 1× simulation over 9 μG in the 10× simulation to 14 μG in the 25× simulation. We compare our predicted magnetic field profiles from our simulations to the best fit to a β model from the observations of the magnetic field in the Coma-galaxy cluster (gray line in the top left panel of Figure 4) and the best fit obtained from 24 simulations of galaxy clusters at the same resolution, and then our 1× run from the same parent dark matter box from Bonafede et al. (2011) (black line in the top panel of Figure 4). While our results for the 1× run are in good agreement with respect to the central magnetic field value in the cluster compared to Coma observations and the simulations of Bonafede et al. (2011) our higher resolution models overpredict the central magnetic field value roughly by a factor of 2.5. We will investigate the origin of this behavior in greater detail in Appendices A and B by varying the magnetic diffusion constant and the initial seed-field strength. Even though our higher resolution simulations predict a central magnetic field strength that is higher than the observed values in the Coma cluster we note that state-of-the-art simulations with Eulerian gird codes typically predict values that just reach the microgauss regime and are around the same factor too low compared to observed values within the Coma cluster that report central field strengths of around 7 μG (Bonafede et al. 2010). Moreover, we note that the cluster that we simulated is not really a Coma-cluster analog as this is a system that is in equilibrium at redshift z = 0 and the Coma cluster is not (e.g., Lyskova et al. 2019). Furthermore, as noted above we do not vary the parameters of the seed field and diffusion constant, which might impact the radial magnetic field distribution at redshift z = 0. We chose to do this in our default simulation runs to obtain pristine conditions for our study of the galactic dynamo, which is the central subject of this paper. Last but not least, the too high central field could also be related to our nonvanishing divergence of the field. We will discuss the impact of the divergence cleaning constraint on the too high central magnetic field strengths in Section 5.6.

Figure 4. We show radial profiles of the density (top left), volume-weighted magnetic field strength (top right), temperature (bottom left), and thermal pressure (bottom right) for our three galaxy-cluster simulations 1× (blue), 10× (red), and 25× (magenta). For the magnetic field, we overplot the best fit from Bonafede et al. (2011) (black line) and observations of the magnetic field within the Coma-galaxy cluster from Bonafede et al. (2010) (gray line).

5.4. Amplification of the Magnetic Field

We start the discussion of magnetic field amplification in our galaxy-cluster zoom-in simulations by considering the time evolution of the magnetic field within 1 virial radius (Rvir) from redshift z = 04–0. We show this in Figure 5. The magnetic field increases exponentially from the initial seed-field value between redshift z = 4 and 2 to a sub-equipartition value of around 0.05 μG in the 1× simulation. In the higher resolution simulations 10× and 25× we find a very different shape of the growth of the field as a function of redshift. Here, the magnetic field in the cluster peaks at around redshift z = 2 at a value of a few microgauss. From that point in time, the field decreases toward redshift z = 0 and settles at around 1–2 μG within the virial radius. This behavior is consistent with other cosmological simulations of magnetic field amplification (see, e.g., Garaldi et al. 2021). While the exponential increase in the magnetic field strength could potentially be related to magnetic field amplification by a small-scale turbulent dynamo driven by subsonic turbulence in the ICM, it is impossible to determine this from the evolution of the magnetic field alone. However, we can still estimate the growth rate of the magnetic field for the different runs. Essentially, we find that all three models are initially consistent with the exponential growth of the form

B∝exp(γ·t).(39)

For the 1× run, we find that γ ≈ 0.7 Gyr⁻¹ while for the 10× and 25× runs we find γ ≈ 0.15 Gyr⁻¹. The former growth rate indicates unresolved dynamo action in the 1× simulation. Furthermore, we note that the increase in the field we observe toward higher redshift is roughly consistent with an increase in the magnetic field toward higher redshift following the relation

B=Bz=01+zm,(40)

with a power-law index of around m = 0.5. This is in relatively good agreement with the predictions made for ska by Krause et al. (2009). It is intrinsically complicated to identify dynamo action in numerical simulations of galaxy and galaxy-cluster formation. This is mainly due to the fact that the fundamentals of dynamo theory are built on top of the theory of turbulence, which is generally not very well understood in HD numerical simulations. Generally, dynamos work by converting (turbulent) kinetic energy into magnetic field energy on the scale of small turbulent eddies. This process is saturated once equipartition between turbulent kinetic energy and magnetic field energy is reached. The magnetic field energy can then be transported to the larger scales in the so-called inverse turbulent cascade. However, the amplification of tiny magnetic seed fields by turbulence competes with the dissipation of magnetic fields on the smallest scales. Only if the interplay between the dissipation of the magnetic field and transport of the magnetic field alongside its amplification is modeled correctly, the dynamo will transit from the linear growth regime, into the nonlinear regime and finally saturate. The crux in achieving this is to have enough resolution on small scales to capture magnetic field amplification by turbulence but also enough resolution on larger scales to model magnetic field transport toward larger structures. This has been the subject of MHD research in many Eulerian grid codes in recent years (e.g., Ryu et al. 2008; Beresnyak & Miniati 2016; Schekochihin et al. 2004; Cho et al. 2009; Porter et al. 2015; Vazza et al. 2018) but there is little to no work on magnetic field amplification in Lagrangian methods. This paper is explicitly targeted to close the gap between the state of research in studies of magnetic field amplification within the ICM that has been put forward in recent years with Eulerian codes. In the following, we will present evidence of an acting small-scale turbulent dynamo in the ICM of our simulated galaxy clusters and evaluate the resolution dependence of the process by directly comparing it to the dynamo theory that has been put forward by Kraichnan & Nagarajan (1967) and Kazantsev (1968) and has been refined by several authors since then (e.g., Zel’dovich 1983; Kazantsev et al. 1985; Kulsrud & Anderson 1992; Kulsrud et al. 1997; Subramanian & Barrow 2002; Xu & Lazarian 2020).

Figure 5. We show the time evolution of the total magnitude of the magnetic field as a function of redshift. This indicates exponential growth of the magnetic field early in the formation history of our galaxy-cluster simulations. We indicate the different resolution levels with red (1×), blue (10×), and magenta (25×). Bottom: we show the evolution of the magnetic field energy for our cluster for all resolution levels and only find a very weak dependence on the resolution.

First, one can study the structure of the cluster in the density-magnetic field phase space to gauge the dependence of the magnetic field on its environment. We show this in Figure 6 for our highest resolution simulation (25×) and two different redshifts, redshift z = 0.3 (left) and redshift z = 0 (right). The gas cells are selected within 1 Rvir around the center of the cluster, which has been identified with the subfind algorithm (Springel et al. 2001; Dolag et al. 2009). The color code indicates the cell mass and shows how much mass is contained in each state in 10¹⁰ M⊙. We deliberately choose these two points in time to distinguish between the linear and nonlinear dynamo regimes at redshifts z = 0.0 and 0.3, respectively. Within this time frame, the cluster transits from a turbulent merging epoch toward a dynamically relaxed system. At both redshifts, we can identify gas at low magnetic field strengths and lower gas densities that are in good agreement with the power-law scaling obtained from the flux-freezing regime of an ideal MHD of an adiabatically collapsing system (B∝ρ23). At larger densities and higher magnetic fields, we can identify a different scaling, specifically at redshift 0. We find excellent agreement with the saturated dynamo regime with the scaling B∝ρ12 in the framework of reconnection diffusion (see, e.g., Xu & Lazarian 2020). While the saturation regime is in very good agreement with our redshift z = 0 results, this is not the case at z = 0.3 when the system undergoes a merger event followed by rapid smooth accretion of gas mass toward the cluster center. While some gas at lower densities is still following the saturation regime, there is a clear deviation in the high-density tail, which is identifiable as flattening of the scaling followed by a kink within the distribution at a density around ne = 3 × 10⁻³ cm⁻³. This could be evidence of rapid diffusion of the magnetic field at the highest densities, which would be in agreement with the recent theory proposed by Xu & Lazarian (2020), who study the turbulent dynamo in the framework of reconnection diffusion under gravo-turbulence in cooling star-forming cores. While Xu & Lazarian (2020) point out that the theory they develop could be of paramount importance in the regime of the first star formation, it is apparent that the idea of gravo-turbulence is of importance on galaxy-cluster scales as well. Thus, first we want to point out that the physical systems of a gravitational collapsing star-forming core are quite different from the cosmological assembly of a galaxy cluster. In a star-forming core, the idea would be that cooling enhances the collapse as the heat generated by the collapse can efficiently be radiated away. This means, that the system is heavily driven out of equilibrium. However, galaxy clusters are (to first order) virialized and thus in equilibrium, especially if one considers non-radiative simulations of clusters. There is one exception to this, which is when the cluster is undergoing a merger process and the merger remnant continues to accrete material that indirectly mimics the situation for which Xu & Lazarian (2020) derive their dynamo model. Xu & Lazarian (2020) derive the following scaling for the nonlinear growth regime of the dynamo under gravitational collapse:

Bρ23∝ρ257−16.(41)

Xu & Lazarian (2020) compare their derived scaling to the results of simulations of magnetic field amplification of first star formation taken from Sur et al. (2012) and find good agreement of their scaling relation with collapsing star-forming structures. However, this is a scale-free problem and can easily be tested in the framework of our cosmological galaxy-cluster simulations.

Figure 6. We show the evolution of the density-magnetic field strength phase space at redshift 0.3 (left) at which the system undergoes a heavy major merger and at redshift z = 0 at which the system transits toward a relaxed state. We include all the gas within the virial radius Rvir. In both cases, we find excellent agreement with the adiabatic compression limit (B ∝ ρ2/3, blue line) at low magnetic field strengths and lower densities. This can be associated with the gas that is falling toward the center of the structure. However, in the regime of higher magnetic fields and higher densities, we find good agreement with the power-law scaling that is expected from a saturated turbulent dynamo (B ∝ ρ1/2) at redshift 0, while there appears to be some deviation from the saturated dynamo at redshift z = 0.3, which could hint at a nonlinear dynamo regime.

Sur et al. (2012) point out that dynamo amplification under gravitational collapse can better be quantified by evaluating B/ρ2/3 than just by evaluating the phase space of magnetic field strength and density in star formation simulations (they simulate the gravitational collapse of a Bonnor–Ebert sphere). However, the physics driving dynamo amplification in the regime of star formation is quite similar to the formation scenario of a galaxy cluster as it undergoes collapse in the dark matter potential and one can directly test the linear growth regime and the nonlinear growth regime discussed in Sur et al. (2012) and Xu & Lazarian (2020), respectively, in the fashion suggested by Sur et al. (2012) in the regime of the formation of a massive galaxy cluster.

This is evaluated in Figure 7 where we show the average of B/ρ2/3 as a function of the density in equal log bins for our 1× (blue), 10× (red), and 25× (magenta) simulations for redshift z = 0.3 (top) and z = 0 (bottom). In this context, we can identify the linear growth regime as the monotonically increasing part as a function of density. For the 1× (lowest resolution) simulation, we find a steeply decreasing part at high densities indicating an unresolved nonlinear growth regime at both redshifts of interest. While we find a similar situation for redshift z = 0 in our two high-resolution simulations 10× and 25×, this is different at a higher redshift of z = 0.3 where we can clearly identify the nonlinear growth regime of the dynamo following the scaling of Xu & Lazarian (2020). We note that the disagreement with Xu & Lazarian (2020) at redshift z = 0 can be explained by the fact that our cluster at hand is a dynamically relaxed system at that time. At z = 0.3, however, the system is strongly collapsing following a previous major merger, providing the ideal conditions for magnetic field amplification via the refined theory of Xu & Lazarian (2020).

Figure 7. We show the quantity B/ρ2/3 as a function of the density in equal log bins to show the agreement of the theory of the turbulent dynamo under the gravitational collapse derived by Xu & Lazarian (2020) with our cosmological galaxy-cluster simulations at redshift z = 0.3 at which the cluster undergoes gravitational collapse after a major merger(top). At redshift z = 0 when the cluster is dynamically relaxed we find deviations from the above scaling. The colors indicate our 1× (blue), 10× (red), and 25× (magenta) runs.

Moreover, we can directly quantify how well our results agree with the scalings derived by Xu & Lazarian (2020) by directly fitting the relevant part shown in Figure 7 between 0.002 and 0.01 cm⁻³ with a function of the form log(Bρ2/3)=αlogn. We report the best-fit parameters for the redshifts z = 0 and 0.3 in Table 3 for the runs 10× and 25×, since these are the only resolution levels that exhibit the reported behavior by Xu & Lazarian (2020). This confirms our findings, that we are not able to reproduce the reported scaling of Xu & Lazarian (2020) if the cluster is relaxed at redshift z = 0, but the cluster can be fairly well characterized by the theory of Xu & Lazarian (2020) when the cluster undergoes heavy merging at around redshift z = 0.3. However, we note that the 10× run is still overestimating the slope and only the 25× run is capable of reproducing a slope that is in good agreement with the theory of reconnection diffusion under gravitational collapse.

Table 3
Best-fit Values for α as Obtained by Our Simulations

Fitting Parameters
| 10× (z = 0) | 25× (z = 0) | 10× (z = 0.3) | 25× (z = 0.3)
α | −0.34 ± 0.02 | −0.29 ± 0.03 | −0.25 ± 0.03 | −0.16 ± 0.02

However, despite the agreement with Xu & Lazarian (2020) we already saw the indication for this behavior shown in the left panel of Figure 6 as the kink in the phase-space distribution at a density of roughly ne 3 × 10⁻³ cm⁻³, which one can interpret as dissipation of magnetic field in high-density regimes.

Apart from the density-magnetic field strength phase space, there is another way of identifying the small-scale-turbulent dynamo by evaluating the magnetic power spectra of the simulations. This has become a standard test for identifying an acting small-scale turbulent dynamo in numerical simulations, specifically in the ISM (e.g., Balsara et al. 2004; Schekochihin et al. 2004; Porter et al. 2015; Hennebelle & Iffrig 2014; Gent et al. 2021) but has recently also become quite popular on the scales of galaxies (e.g., Butsky et al. 2017; Martin-Alvarez et al. 2018, 2020; Pakmor et al. 2017; Rieder & Teyssier 2016, 2017a, 2017b; Steinwandel et al. 2019, 2020a) and galaxy clusters (e.g., Dubois & Teyssier 2008; Ryu et al. 2008; Vazza et al. 2018) for studying the turbulent dynamo. We show the magnetic power spectra for our three MHD simulations in Figure 8 for 1× (blue), 10× (red), and 25× (magenta). We note that these are simply the redshift z = 0 power spectra. For each simulation, we calculated the power spectra by binning the SPH data to a grid. The grid has a resolution of 128³ and is represented by a cube with a side length of 3 × Rvir. The power on each scale is then computed by evaluating the Fourier modes on the grid scale.

Figure 8. We show the magnetic power spectra at redshift z = 0 for all three resolution levels in blue (1×), red (10×), and magenta (25×) of our galaxy-cluster zoom-in simulations. We see very good agreement with the predicted slope from the dynamo theory by Kazantsev (1968) on large scales (golden line) in our 10× and 25× resolution runs. However, we note that the 1× model predicts a slightly steeper slope than the k3/2 slope from Kazantsev (1968), which is supposedly related to a lack of resolution. Furthermore, we overplot the k−5/3 and k⁻² slopes for reference on smaller scales. We find a steeper slope on smaller scales than predicted by these scaling, which is in accordance with small-scale simulations of the turbulent dynamo (see, e.g., Schekochihin et al. 2004; Porter et al. 2015).

One can clearly see that there is excellent agreement between the power spectra determined by this methodology and the power-law slopes predicted by Kazantsev (1968) for the simulations 10× and 25×. In these two simulations, we find the best fit for the power-law slope between 1.2 and 1.5. However, we note that there is some variation of this depending on the exact resolution of the grid we adopt for the underlying power spectrum. For the 1× run, we find a fit best-fit value for the slope of around 0.9, which is much lower than the slope of 1.5 predicted by the theory of Kazantsev (1968), indicating that the inverse cascade of the small-scale turbulent dynamo is only poorly resolved in the 1× run. This is probably introduced by the lack of resolution we have in the outer parts of the cluster in the 1× run compared to the 10× and 25× runs. Therefore, we suggest that studies regarding the dynamo on the scales of the ICM require a cell mass resolution of around 10⁷ M⊙, which corresponds to a spatial resolution of around 2 kpc. While we generally find little difference between our 10× and 25× models, we would advise future dynamo studies with Lagrangian methods to adopt our 25× resolution for converged results on the power spectra, which results in a mass resolution of around 4 × 10⁶ and a spatial resolution of roughly 3 kpc. This is roughly in line with the findings of Vazza et al. (2018) for the grid code enzo, who obtain self-consistent power spectra for their two highest resolution runs. However, we note that while the Kazantsev (1968) spectra are recovered they are not converged and with increased resolution we find a shift of the spectra to the right, indicating that the growth rate is still increasing when we move from 10–25×. Additionally, we would like to add a brief analysis based on the findings of Schober et al. (2015), who present in their Table 1 some scalings for the shape of the magnetic power spectra with under resolved dynamo action. First, we derive the critical magnetic Reynolds number, RMcrit, in our simulations for which we generally find a value of around 1000, depending on the exact numbers for gas density, magnetic field energy density, as well as an estimate of the turbulent length scale and the turbulent velocity. For our estimate, we assume a gas density of 0.01 cm⁻³, a magnetic field strength of 1 μG as well as a typical response time of a cluster of a few tens of megayears, based on the growth rate of the magnetic field. Furthermore, we assume that typical velocities of turbulence are around 100 km s⁻¹ and the turbulence length scale in the ICM might be a few hundred kiloparsecs. This yields (depending on the exact values plugged in) an effective magnetic Reynolds number between 1000 and 10,000 in our highest resolution run. Compared to Table 1 in Schober et al. (2015), this implies that we resolve the effective magnetic Reynolds number for compressive, intermittent, and compressible turbulence as outlined. Furthermore, we can compute that labeled by Schober et al. (2015) as k⋆ for which we generally find a value between 0.05 and 0.1, again dependent on the exact values of turbulence that one might assume for the ICM. Hereby, we note that our highest resolution runs roughly reproduce the peak value of the power spectrum of that obtained out of Kolmogorov turbulence in comparison to Table 1 in Schober et al. (2015) in the case of PM > 1, further ensuring that we can resolve dynamo action (not the growth rate) in simulations of massive galaxy clusters.

Furthermore, we investigate the power spectra at four earlier times than redshift z = 0 for our 25× run. We show the results in Figure 9 for z = 0 (black line), z = 0.3 (light black line), z = 0.6 (gray line), z = 1.1 (light gray line), and z = 2.3 (very light gray line). Essentially, we find similar results as in the redshift z = 0 case with the difference that there is less power stored in the magnetic field at higher redshift. Despite the lower power in the magnetic field, we are still able to recover a power spectrum that is in accordance with the theory of Kazantsev (1968) already at redshift z ∼ 1, after half of the cosmological evolution of the cluster itself. However, we note an exception to this at redshift z = 2.3 where we find more power in the magnetic field compared to the lower redshift spectra. This is roughly consistent with the peak in the time evolution of the magnetic field strength at around redshift z = 2 and generally also consistent with expectations for higher redshift clusters (e.g., Krause et al. 2009).

Figure 9. We show the magnetic power spectra for four different redshifts of our highest resolution run 25×. For all presented redshifts we find very good agreement with the predicted slope from Kazantsev (1968) on larger scales (gold line). We overplot the k−5/3 and k⁻² slopes for reference on smaller scales.

Most studies on the turbulent dynamo on galaxy and galaxy-cluster scales stop at the point where they achieve the scalings predicted from the phase-space structure (see Figure 6) and the power spectra (see Figures 8 and 9). Schekochihin et al. (2004) pointed out quite early that the power spectra alone might not suffice to clearly identify dynamo action. Thus, they suggest a different (stronger scaling) based on the curvature of the magnetic field lines given via

K=(B·∇)B∣B2∣,(42)

which can be rewritten as

K=1B212∇(B·B)−B×(∇×B).(43)

We calculate this quantity as an additional output field in the code on the fly. We show the relation of the magnetic field strength as a function of the curvature of the field lines shown in Figure 10 and note that while our calculations are slightly too steep in the higher resolution models 10× and 25× we recover the declining trend of the field strength with the curvature following roughly KB1/2=constant. Generally, this is a good sign as this indicates that the increasing field strength is counteracting the bending of field lines by magnetic tension. Thus, the bending of field lines is suppressed by magnetic tension and the dynamo saturates in the regime of a few microgauss, as expected from the theory of small-scale-turbulent dynamo. The fact that our results are too steep could be related to our slightly too high magnetic field strengths in the cluster center. Therefore, this could in our case be related to some limitations of the model, which we will discuss in detail in Section 6.2. Nevertheless, we raise two additional points about the curvature. First, this relation is mainly inferred from high-resolution plasma physics simulations without the presence of self-gravity. Thus, it is a priori not clear why a galaxy cluster would exactly follow this relation as the gravitational collapse of structure will add an additional imprint to the curvature relation. Moreover, we note that this is not so different from that of Vazza et al. (2018), who are the only other group who ever checked this relation, where they also find a slight deviation from the results of Schekochihin et al. (2004). We strongly suspect that the gravitational collapse of structure is responsible for the change in the magnetic curvature relation. However, we cannot prove this statement as this requires a detailed study of high-resolution plasma physics simulations like those of Schekochihin et al. (2004), which include the effect of self-gravity.

Figure 10. We show the relation between magnetic field strength and magnetic field line curvature for all of our simulation runs 1× (blue), 10× (red), and 25× (magenta). We compare our results to the scaling derived from the high-resolution idealized dynamo simulations of Schekochihin et al. (2004). We can recover the decreasing trend following KB1/2=constant (black line). However, we find a tilt and a slightly steeper slope than expected in our higher resolutions simulations for the models 10× and 25×.

5.5. Probability Distribution Function of the Magnetic Field

Another aspect that is interesting when it comes to magnetic fields in galaxy clusters is the probability distribution function (PDF) of the field strength within the cluster. In this context, there are two very important questions to answer. First, how does the magnetic field distribution change with the resolution of the simulation, and second, how does it change as a function of time within the cluster region?

In Figure 11, we show the magnetic field PDF at redshift z = 0 for our cluster at the three targeted resolution levels 1× (blue), 10× (red), and 25× (magenta). In the lowest resolution run at 1×, we can see that the magnetic field distribution peaks at a low value between 10⁻⁵ and 10⁻⁴ μG, even at redshift z = o and there is only a very small fraction of the volume that reaches a magnetic field up to a few microgauss. This picture changes in the higher resolution runs, where we can indicate a small peak between 1 and 10 μG, which shows that roughly 10% of the clusters’ total volume reaches significant magnetic field strengths that are in agreement with observed magnetic fields in galaxy clusters.

Figure 11. We show the volume-weighted PDF of the density for all three models, 1× (blue), 10× (red), and 25× (magenta), within the virial radius of our simulated galaxy clusters. For the lowest resolution model 1×, there is a small plateau around the microgauss regime. In the higher resolution runs at 10× and 25× this evolves toward an extended bump showing a strong magnetization of around 20% of the governed volume of the system. There is a sharp drop in the PDF for stronger magnetic field strengths above 10 μG. This indicates that only very few gas cells govern the regime of extreme magnetic field strengths. It is interesting to point out that at low resolution the PDF peaks at relatively lower field strengths of around 10⁻⁹, indicating very inefficient dynamo amplification at the lowest resolution level as the turbulence is not strong enough to significantly amplify the field into the microgauss regime.

Furthermore, we take the time evolution of the PDF of the magnetic field into account by displaying a time sequence of the PDF for our simulation 25× that is color coded against cosmic time starting from an evolutionary state that marks 1 Gyr. We show this in Figure 12. Initially, the magnetic field distribution peaks at a value below 10⁻⁴ μG, which is slightly higher than our initial seed-field value for this run, due to adiabatic compression of the gas during structure formation. At later times, one can identify a clear shift in the PDF from low magnetic field values to high magnetic field values where the peak is shifted the farthest to the right after around 4 Gyr of evolution, which is roughly consistent with the slight peak in the mean magnetic field strength that we could observe for our higher resolution runs shown in Figure 5 at around or slightly before a redshift of z = 2. The peak shifts further to the left again with decreasing redshift and peaks at a value of around 7 μG by redshift z = 0. Furthermore, we note that by combining Figures 11 and 12 the dynamo action can be identified by a transition of the volume-weighted PDF that peaks around the initial seed-field value; at high resolution, the field is redistributed to higher field strengths in the microgauss regime by the dynamo, while at low resolution the volume filling phase remains at the seed-field value, even at redshift z = 0 due to unresolved dynamo action.

Figure 12. We show a time evolution of the volume-weighted PDF of the magnetic field starting after 1 Gyr of the evolution of the cluster for our highest resolution run. Initially, we find a lot of the gas at very low field strengths around 10⁻¹⁰ as this corresponds to the initial value of our physical field. The dynamo action that is induced due to the turbulence introduced by structure formation then amplifies the magnetic field as time progresses. This redistributes the volume filling phase of the field to the higher magnetic field regime, which is generated by the dynamo with decreasing redshift. Furthermore, the PDF shows higher field values at redshift z ∼ 2 compared to z ∼ 0.

5.6. Divergence Cleaning Constraint

Finally, we discuss the divergence constraint that is of importance for MHD simulations with both particle and grid-based methods. While in grid-based codes it is possible to enforce the divergence constraint by using the constrained transport (CT) scheme in which one computes the cell-centered magnetic field from the electric field on the edges of each cell, it is not yet clear whether this can be done similarly in particle codes. Regarding a CT scheme, one should keep two constraints in mind. First, strictly speaking with a CT scheme one enforces that the numerical realization of a physical field is divergence-free, not the actual physical field (i.e., the field is divergence-free in the projected grid geometry). Second, and more importantly, the divergence error of the magnetic field in a simulation directly traces the accuracy of the integration of the induction equation within the simulation. Thus, if present, the nonzero divergence gives a direct estimate of the integration error of the underlying integration scheme used for the induction equation. This information is at least partially lost in a CT scheme. While the first point is rather particular, and strictly speaking, this is true for every numerical realization of any physical system, the second one is rather important as it effectively measures how well the numerical scheme is handling the complexity of the MHD equations. As pointed out above, it is not clear if such a scheme can be constructed fully in Lagrangian code and thus particle codes (such as ours) rely on divergence cleaning schemes to control the error introduced by a nonzero divergence. In our simulations, we use the eight-wave Powell et al. (1999) cleaning scheme.

In Figure 13, we show the column density (left), relative divergence (middle), and the magnetic field (right) for our simulation 25× at redshift z = 0. While we can see that the high divergence regions track the regions with high magnetic fields, the mean of the absolute value of the divergence stays reasonably low. Nevertheless, we want to point out that the nonzero divergence is the potential origin of the high central magnetic fields that we observe in our higher resolution runs at redshift z = 0, if we compare it to recent observations of galaxy-cluster magnetic fields.

Figure 13. The gas surface density can be visualized on the left, the divergence of the magnetic field (h denotes the smoothing length) in the center, and the magnetic field strength on the right for our highest resolution simulation 25×. We note that our highest magnetic field strengths in the center are associated with the highest divergence.

However, given the shape of our temperature and density profiles observed in all clusters we believe it is more likely that we simply operate in a regime of adiabatic compression rather than magnetic field amplification in a spurious dynamo. Furthermore, we investigate the divergence of the magnetic field in the most massive structure within our simulations within the virial radius Rvir as a function of the density. We show our results in Figure 14 for redshift z = 0. In both Figures discussed so far, we measure the divergence as the relative divergence:

∇·Brel=∇·BhB,(44)

where h is the smoothing length of the kernel in our SPH simulation. We solve for this quantity on the fly by computing:

∇·Brel,i=∇·Bihi+hjBi+Bj·Bihi.(45)

in our simulation code, to obtain the most accurate description of our divergence error within the simulations. We directly write this data into our simulation snapshots. As can be seen from Figure 14, our relative divergence error remains very low for our 25× resolution simulations, with the bulk of our particles showing a relative divergence error of around 10⁻⁹. This is a very good result and shows that our code is capable of handling the divergence constraints in astrophysical MHD simulations to a sufficient amount. However, we note that there is an extended tail of the relative divergence error reaching up toward 1 for single particles at very high densities. These also trace the particles with the highest magnetic field in our simulations showing a field strength of around 100 μG. This fact could be related to the fact that our higher resolution simulations predict slightly too high central magnetic fields within the cluster. Still, given the distribution of the relative divergence error in our simulations, we can rule out magnetic monopoles as the primary amplification mechanism of the magnetic field in our simulations.

Figure 14. We show the density-divergence phase space color coded against the mass in the bin each bin. Our simulations show excellent divergence cleaning properties. We show the innermost 2 Mpc of the most massive structure within our simulation. One can clearly see that the bulk of the material shows excellent divergence cleaning properties with a relative divergence between 10⁻¹¹ and 10⁻⁵. However, we note that single particles reach very high relative divergence around 10⁻¹, while there are outliers reaching up to 1. Despite this, the excellent relative divergence makes it very unlikely that the magnetic field is solely amplified by magnetic monopoles.

Figure 15. We show the time evolution of the mean value of the divergence for all three resolution levels of our galaxy-cluster simulation with 1× (blue), 10× (red), and 25× (magenta). We can clearly see that the divergence of the field remains very low until the cluster starts to grow a significant field around redshift z = 4. The divergence fluctuates around zero with peak values reaching 10⁻⁵, which is an excellent value compared to other galaxy-cluster zoom simulations of the same kind (e.g., Vazza et al. 2018).

Moreover, we want to gauge the evolution of the divergence constraint as a function of time by showing the density-divergence phase space for the three different redshifts shown in Figure 15, ranging from redshift z = 2 over redshift z = 1 to redshift z = 0 on the right. Finally, we note that the results for the divergence of the magnetic field that we obtain in our galaxy-cluster zoom-in simulations are excellent, even compared to the results obtained with grid codes like enzo, ramses, and arepo, which can, for example, be seen by comparing our divergence-density phase space shown in Figure 14 with the results obtained by Pakmor et al. (2020). The origin of this is of our improved divergence constraint compared to previous SPH simulations can be understood as follows. First, all of our simulations are carried out in a non-radiative fashion, which avoids the cooling-driven collapse of high-density regions within the ICM into galaxies, which typically would host the regions with the largest divergence in the simulation domain. However, even compared to the non-radiative simulations of Vazza et al. (2018) we obtain a result that shows around an order of magnitude lower divergence in the radial trend, which shows that our simple Powell-cleaning scheme (Powell et al. 1999) is sufficient to capture the emergence of the magnetic field in the ICM, at least in non-radiative simulations. Moreover, we note that the divergence is already heavily suppressed compared to older SPH simulations just by adopting a modern form of SPH that makes use of the higher-order kernel first discussed by Wendland (1995). Essentially, one can imagine that the use of a higher-order kernel in SPH is similar to increasing the convergence order of the code while decreasing the spatial resolution.

Second, we derive the relative divergence of the magnetic field in an SPH-like fashion in our simulation code following Equation (45), which is the correct way of obtaining the relative divergence in SPH simulations. If the divergence is derived from the particle data alone by weighting it with the smoothing length and normalizing it to the absolute value of the magnetic field it can be heavily over- or underestimated, depending on the exact averaging that is applied, which is usually very sensitive to the outliers in the distribution. Deriving the divergence following Equation (45) is not prone to outliers as it is a kernel-averaged quantity. Therefore, we propose that in all particle-based methods the divergence should always be calculated by Equation (45) to avoid confusion in future SPMHD simulations and provide a cleaner comparison to Eulerian codes that often calculate this quantity by a finite-difference technique (e.g., Vazza et al. 2018).

6. Conclusions

6.1. Summary

We present SPMHD simulations of a massive galaxy cluster with a total mass of M₂₀₀ ∼ 10¹⁵ M⊙ as a study of resolution on three different resolution levels of 1× (∼10⁸ M⊙ per cell), 10× (∼10⁷ M⊙ per cell), and 25× (∼4 × 10⁶ M⊙ per cell). We investigated the structure, morphology, and evolution of the cluster, focused on the amplification of the magnetic field via the small-scale turbulent dynamo and discussed the limitations. The main conclusions of this work are the following:

With increasing resolution the central magnetic field strength in the core of the cluster increases by a factor of ∼3 from the base resolution run at 1× toward the highest resolution run at 25×. We note that this is higher than results obtained with Eulerian methods (e.g., Vazza et al. 2018) and observations of the Coma-galaxy cluster (e.g., Bonafede et al. 2011), but is still not an unrealistic value for cool-core clusters.

We find a steep exponential increase of the magnetic field as a function of cosmic time that flattens at a sub-equipartition value at around redshift z = 2 for our lowest resolution simulation 1× from which it increases at a slower rate as the cluster reaches redshift z = 0 with a field that remains slightly below the microgauss regime. In the higher resolution runs, 10× and 25×, we find that the magnetic fields peak at around redshift z = 2 and the fields saturate at a value of around ∼2 μG in our 10× and 25× models while it stays slightly at sub-equipartition in our 1× model (when one compares the mean in the virial radius at redshift z = 0).

The field increase toward higher redshift is consistent with predictions for ska from Krause et al. (2009).

We find strong evidence that the magnetic field is amplified by the small-scale-turbulent dynamo in the ICM, driven by turbulence introduced by mergers, shocks, and cosmic accretion. For the first time, we were thus able to unravel the nonlinear regime of the dynamo driven by gravo-turbulence in agreement with the recent theoretical model, developed by Xu & Lazarian (2020). Furthermore, we show evidence of the dynamo by the magnetic power spectra that take the form predicted by Kazantsev (1968) and investigate the dependence between magnetic field strength and field line curvature and find good agreement with the results of Schekochihin et al. (2004) and Vazza et al. (2018).

Finally, we analyzed the behavior of the divergence constraint in our simulations and find that while the divergence of the field increases with increasing resolution (which is the potential origin of our slightly too large central magnetic field strengths in the cluster center) it is in good agreement with results presented with state-of-the-art moving-mesh codes like arepo (e.g., Pakmor et al. 2020) on galaxy scales and with state-of-the-art grid codes like enzo on cluster scales (e.g., Vazza et al. 2018).

6.2. Model Limitations

Finally, we want to briefly discuss the consequences of the limitations of our modeling. First, we carry out the simulations without the effects of cooling, star formation, and feedback. This might be an important restriction as one could imagine that the magnetic field is first amplified during the collapse of halos to proto-galaxies and later via galactic dynamos in the galaxies themselves, which could redistribute their intergalactic magnetic field to the CGM and subsequently the ICM via galactic outflows, driven by AGN and SNe. Specifically, the lack of AGN feedback could be related to our too high central magnetic field strengths in our higher resolution runs 10× and 25×. In that scenario, the high magnetic fields in the center could efficiently support accretion toward a central SMBH, centered in the brightest cluster galaxy that could have formed under cooling in the center of the cluster. The accretion of material and the subsequent outflow of material from the AGN could efficiently transport the magnetic field outward and contribute to the magnetization of the void regions surrounding the cluster. Throughout all the simulations we carried out, we build up an approximate microgauss magnetic field in the cluster center. The redistribution of the magnetic fields from the center to the outer parts of the cluster via AGN feedback could also potentially help in decreasing the fast drop in the radial profile of the magnetic field in the outer parts of the cluster.

However, the aim of our study is to investigate magnetic field amplification in the ICM by the turbulent dynamo, driven by turbulence injected via shocks during structure formation processes. In order to understand how magnetic fields are growing in a cosmological context on large scales, we deliberately ignored these effects and ran the simulations in a non-radiative fashion.

Moreover, we did not discuss the origin of the magnetic field in our simulations, for example, via the implementation of the Biermann battery or other battery processes driven by reionization (Garaldi et al. 2021). While the former could provide an interesting scenario for a self-consistent treatment of seed magnetic fields that are injected during the structure formation process by an offset between pressure and density gradient, the latter remains unimportant as long as we carry out the simulations in a non-radiative fashion. However, while we do not self-consistently implement a process like the Biermann battery we tested different magnetic field strengths as the initial seed for our 10× runs with an order of magnitude difference and find very similar results. We discuss this in more detail in Appendix A. Additionally, we note that while we run all the simulations effectively in a nonideal MHD limit we adopted a constant diffusion coefficient on all resolution scales. The change in the diffusion coefficient could have a potential effect on the magnetic field distribution, especially in the cluster center where we find that the magnetic field strength is a factor of 2.5 higher compared to observations of the Coma-galaxy cluster. We test the impact of a varying diffusion coefficient by increasing and decreasing by an order of magnitude in Appendix B. Furthermore, although we run MHD simulations we adopted isotropic conduction ignoring that magnetic fields can alter the thermal transport process in an anisotropic fashion. While we explicitly state this a caveat, we tested the impact of anisotropic thermal conduction on basic cluster properties as well as basic dynamo properties and found only little differences when we include the effect for our 10× and 25× simulations. We discuss some of our findings in Appendix C.

In addition, we also did not take into account other nonthermal effects like cosmic-ray protons and cosmic-ray electrons, which have a potentially important impact on the magnetic field in terms of amplification and structure (Buck et al. 2020; Butsky et al. 2020; Hopkins et al. 2020a, 2020b, 2020c, 2020d).

We would like to thank the anonymous referee for a great report that helped us understand more details on the subject of dynamo theory. U.P.S. is supported by the Simons Foundation through a Flatiron Research Fellowship (FRF) at the Center for Computational Astrophysics. The Flatiron Institute is supported by the Simons Foundation. U.P.S. acknowledges the computing time provided by the resources at the Flatiron Institute on the cluster rusty. U.P.S. acknowledges the computing time provided by the Leibniz Rechenzentrum (LRZ) of the Bayrische Akademie der Wissenschaften on the machine SuperMUC-NG (pn72bu). U.P.S. acknowledges the computing time provided by c2pap (pr27mi). K.D., L.M.B., and U.P.S. acknowledge the computing time provided by the Leibniz Rechenzentrum (LRZ) of the Bayrische Akademie der Wissenschaften on the machine SuperMUC-NG (pr86re). This research was supported by Excellence Cluster ORIGINS, which is funded by the Deutsche Forschungsgemeinschaft (DFG; German Research Foundation) under Germany’s Excellence Strategy—EXC-2094-390783311. K.D. acknowledges funding for the COMPLEX project from the European Research Council (ERC) under the European Unionʼs Horizon 2020 research and innovation program grant agreement ERC-2019-AdG 860744. We thank the LRZ in Garching for the supercomputing resources and for an energy mix that is 100% comprised of renewable energy resources (e.g., https://www.top500.org/news/germanys-most-powerful-supercomputer-comes-online/, https://www.lrz.de/wir/green-it_en/).

Software: We used the cosmological simulation code gadget3 (Springel 2005; Dolag & Stasyszyn 2009; Beck et al. 2016) to run the simulations and use the language julia (Bezanson et al. 2014) to perform the analysis⁸ based on the packages that can be found at: https://github.com/LudwigBoess.

Data Availability

The data will be made available based upon a reasonable request to the corresponding author.

Appendix A
Dependence on the Seed Field

The choice of the adopted seed field at high redshift is somewhat random, and in fact, it is not our preferred choice of initializing a magnetic field as we would much rather seed the field by SNe. Generally, we think that the SN-seeding scenario is physically better motivated than the choice of a pseudo-random seed field. However, as the runs are adiabatic, specifically to obtain a closer comparison to the work of Vazza et al. (2018), we do not have cooling or star formation and cannot use our SN-seeding approach that has been employed in earlier work (e.g., Beck et al. 2012; Steinwandel et al. 2019, 2020a, 2020b). Nevertheless, we can investigate and point out the key changes that we observe for a 10 times higher and a 10 times lower initial magnetic seed field. First, we note that the choice of the seed field only marginally affects the magnetic field growth in most quantities. Thus, we are mostly interested in the radial trend of the magnetic field, which we find to be too high compared to the observations of Coma (see Section 5.3). This is a known issue in Lagrangian methods (e.g., Donnert et al. 2018, for a review) and it has often been suspected that this is related to the nonzero divergence behavior in Lagrangian methods. Second, we are also interested in the magnetic field curvature relation as this represents the direct imprint of the dynamo on the structure of the magnetic field lines, while power spectra often take a very generic form and it remains unclear what drives their power-law behavior, as we discussed in Section 5.4. Therefore, these are the two quantities we want to focus on here. We show the radial profiles for density, magnetic field, temperature, and thermal pressure in Figure 16. There are two notable changes in the magnetic field structure that can be observed as the seed field is increased. The magnetic field grows roughly by a factor of 2 in the very center of the cluster with increasing seed field, while the variations in the outskirts of the cluster are marginal but we still find a slightly larger field. We note that while we find a relatively strong field for our 10×-high-seed run with the higher seed field but also note that we find the lowest temperature in the cluster center for this run. The drop in the temperature in the cluster center is from the adiabatic nature of this set of simulations, possibly due to the absence of an AGN that provides energy to heat the gas.

Figure 16. Same as Figure 4 but for the runs 10× (red), 10×-low-seed (orange), and 10×-high-seed (green) to demonstrate the effect of a variation in the initial seed field on the evolution of the cluster.

Therefore, we note that to increase the cluster temperature in the center alone would yield a drop in the central field by roughly half a dex.

In Figure 17, we show the curvature relation as a function of the different seed fields and note that we find only very minor changes compared to our default runs. Nevertheless, the run with the higher magnetic seed-field strength, 10×-high-seed, captures the slope of Schekochihin et al. (2004) better at larger magnetic field line curvature.

Figure 17. Same as Figure 10 but for the runs 10× (red), 10×-low-seed (orange), and 10×-high-seed (green) to demonstrate the effect of a variation in the initial seed field on the evolution of the cluster.

Appendix B
Dependence on the Magnetic Diffusivity

Furthermore, we tested the dependence of our simulation results by adopting different magnetic diffusivity constants ηm in our nonideal MHD prescription by increasing and decreasing the constant by a factor of 10. We show the radial profiles that we obtain for these runs in Figure 18. The differences are very minor, and in fact, they are contained within the model scatter at the percent regime. This means that if we run the exact same simulation on a different machine or a different number of nodes, we obtain a similar change in the radial magnetic field distribution as inferred from the change in the diffusion coefficient. This could be interpreted in the following way. The dynamo arises from the interplay between magnetic field amplification and dissipation (on small scales) or diffusion (to larger scales) of the magnetic field. The marginal change that we see by changing the diffusion coefficient tells us that the diffusion and dissipation of the field are of numerical nature rather than physical, which is consistent with other works that do not include the nonideal MHD term to begin with (almost every other galaxy-cluster simulation). However, there is an imprint of the change in the diffusion coefficient, which is most apparent by considering the curvature relation shown in Figure 19. Here, we can clearly see that the magnetic field is redistributed from the regime of smaller curvatures to larger ones, which essentially forces the curvature relation to become less steep at higher field line curvature.

Figure 18. Same as Figure 4 but for the runs 10× (red), 10×-low-eta (orange), and 10×-high-eta (green) to demonstrate the effect of a variation in the magnetic resistivity of the plasma.

Figure 19. Same as Figure 10 but for the runs 10× (red), 10×-low-eta (orange), and 10×-high-eta (green) to demonstrate the effect of a variation in the magnetic resistivity of the plasma.

Appendix C
Dependence on the Conduction Model

Finally, we carry out one additional test on our 10× model, which is the inclusion of anisotropic thermal conduction. Generally, we note that the run with anisotropic thermal conduction produces the largest field in the radial trend of all the test runs (see Figure 20). However, this can be somewhat understood by taking the temperature profile into account. We can clearly see that the run 10×-ani shows the lowest temperature in the cluster center and shows a consistent increase in the magnetic field compared to the decrease in temperature. Moreover, we find a very weak change in the curvature relation by the inclusion of anisotropic thermal conduction, which we show in Figure 21. As a side note, we want to point out that these simulations are computationally expensive because we carry them out with physical conduction. The conduction module takes roughly 20% of the computing time, which is a considerable computational effort.

Figure 20. Same as Figure 4 but for the runs 10× (red) and 10×-ani (orange) to demonstrate the effect of anisotropic thermal conduction in the plasma.

Figure 21. Same as Figure 10 but for the runs 10× (red) and 10×-ani (orange) to demonstrate the effect of anisotropic thermal conduction in the plasma.

In conclusion, we note that the origin of the large magnetic fields in the center that are too high by a factor of 2.5 remains an unresolved issue in particle codes. We carefully checked that the divergence can almost be ruled out at this point as we show almost two orders of magnitude better behavior compared to similar simulations that report magnetic fields that are a factor of 2 too low compared to observations in the Coma cluster. Therefore, we need a more dedicated numerical study of the origin of the larger fields in Lagrangian codes, which will include model variations that include more physics, such as Braginskii viscosity. Moreover, we need to test more in-depth the diffusion of the field and whether it is introduced by the cleaning scheme or the nonideal MHD prescription. This will be the subject of future work.

ORCID iDs

Ulrich P. Steinwandel orcid https://orcid.org/0000-0001-8867-5026

References

↑

Arth, A., Dolag, K., Beck, A. M., Petkova, M., & Lesch, H. 2014, arXiv:1412.6533
ADS

↑

Balsara, D. S., Kim, J., Mac Low, M.-M., & Mathews, G. J. 2004, ApJ, 617, 339
Crossref ADS

↑

Basu, A., & Roy, S. 2013, MNRAS, 433, 1675
Crossref ADS

↑

Beck, A. M., Lesch, H., Dolag, K., et al. 2012, MNRAS, 422, 2152
Crossref ADS

↑

Beck, A. M., Murante, G., Arth, A., et al. 2016, MNRAS, 455, 2110
Crossref ADS

↑

Beck, R. 2015, A&ARv, 24, 4
Crossref ADS

↑

Beck, R., & Krause, M. 2005, AN, 326, 414
Crossref ADS

↑

Beresnyak, A., & Miniati, F. 2016, ApJ, 817, 127
Crossref ADS

↑

Berlok, T., Quataert, E., Pessah, M. E., & Pfrommer, C. 2021, MNRAS, 504, 3435
Crossref ADS

↑

Bezanson, J., Edelman, A., Karpinski, S., & Shah, V. B. 2014, arXiv:1411.1607
ADS

↑

Biermann, L. 1950, ZNatA, 5, 65
ADS

↑

Böhringer, H., Chon, G., & Kronberg, P. P. 2016, A&A, 596, A22
Crossref ADS

↑

Boldyrev, S., & Cattaneo, F. 2004, PhRvL, 92, 144501
Crossref ADS

↑

Bonafede, A., Dolag, K., Stasyszyn, F., Murante, G., & Borgani, S. 2011, MNRAS, 418, 2234
Crossref ADS

↑

Bonafede, A., Feretti, L., Murgia, M., et al. 2010, A&A, 513, A30
Crossref ADS

↑

Bonafede, A., Vazza, F., Brüggen, M., et al. 2013, MNRAS, 433, 3208
Crossref ADS

↑

Bonafede, A., Feretti, L., Giovannini, G., et al. 2009, A&A, 503, 707
Crossref ADS

↑

Braginskii, S. I. 1965, RvPP, 1, 205
ADS

↑

Brentjens, M. A., & de Bruyn, A. G. 2005, A&A, 441, 1217
Crossref ADS

↑

Brüggen, M., Ruszkowski, M., Simionescu, A., Hoeft, M., & Dalla Vecchia, C. 2005, ApJL, 631, L21
Crossref ADS

↑

Buck, T., Pfrommer, C., Pakmor, R., Grand, R. J. J., & Springel, V. 2020, MNRAS, 497, 1712
Crossref ADS

↑

Burn, B. J. 1966, MNRAS, 133, 67
Crossref ADS

↑

Butsky, I., Zrake, J., Kim, J.-h., Yang, H.-I., & Abel, T. 2017, ApJ, 843, 113
Crossref ADS

↑

Butsky, I. S., Fielding, D. B., Hayward, C. C., et al. 2020, ApJ, 903, 77
Crossref ADS

↑

Cho, J., Vishniac, E. T., Beresnyak, A., Lazarian, A., & Ryu, D. 2009, ApJ, 693, 1449
Crossref ADS

↑

Clark, S. E., Heiles, C., & Robishaw, T. 2019, arXiv:1903.07671
ADS

↑

Clark, S. E., Peek, J. E. G., & Putman, M. E. 2014, ApJ, 789, 82
Crossref ADS

↑

Clarke, T. E. 2004, JKAS, 37, 337
Crossref ADS

↑

Clarke, T. E., Kronberg, P. P., & Böhringer, H. 2001, ApJL, 547, L111
Crossref ADS

↑

Cowie, L. L., & McKee, C. F. 1977, ApJ, 211, 135
Crossref ADS

↑

Crutcher, R. M. 2012, ARA&A, 50, 29
Crossref ADS

↑

Dehnen, W., & Aly, H. 2012, MNRAS, 425, 1068
Crossref ADS

↑

Demozzi, V., Mukhanov, V., & Rubinstein, H. 2009, JCAP, 2009, 025
Crossref ADS

↑

Dolag, K., Bartelmann, M., & Lesch, H. 1999, A&A, 348, 351
ADS

↑

Dolag, K., Bartelmann, M., & Lesch, H. 2002, A&A, 387, 383
Crossref ADS

↑

Dolag, K., Borgani, S., Murante, G., & Springel, V. 2009, MNRAS, 399, 497
Crossref ADS

↑

Dolag, K., Jubelgas, M., Springel, V., Borgani, S., & Rasia, E. 2004, ApJL, 606, L97
Crossref ADS

↑

Dolag, K., Schindler, S., Govoni, F., & Feretti, L. 2001, A&A, 378, 777
Crossref ADS

↑

Dolag, K., & Stasyszyn, F. 2009, MNRAS, 398, 1678
Crossref ADS

↑

Dolag, K., Vazza, F., Brunetti, G., & Tormen, G. 2005, MNRAS, 364, 753
Crossref ADS

↑

Donnert, J., Vazza, F., Brüggen, M., & ZuHone, J. 2018, SSRv, 214, 122
Crossref ADS

↑

Dubois, Y., & Commerçon, B. 2016, A&A, 585, A138
Crossref ADS

↑

Dubois, Y., & Teyssier, R. 2008, A&A, 482, L13
Crossref ADS

↑

Frank-Kamenetskii, D. A. 1967, Diffusion and Heat Transfer in Chemical Kinetics [in Russian], Nauka Moscow (Princeton, NJ: Princeton Univ. Press)

↑

Garaldi, E., Pakmor, R., & Springel, V. 2021, MNRAS, 502, 5726
Crossref ADS

↑

Gent, F. A., Mac Low, M.-M., Käpylä, M. J., & Singh, N. K. 2021, ApJL, 910, L15
Crossref ADS

↑

Gnedin, N. Y., Ferrara, A., & Zweibel, E. G. 2000, ApJ, 539, 505
Crossref ADS

↑

Greaves, J. S., Holland, W. S., Jenness, T., & Hawarden, T. G. 2000, Natur, 404, 732
Crossref ADS

↑

Guthrie, G. J., Norton, W. J., & Schomaker, V. 1949, JChPh, 17, 310
Crossref ADS

↑

Heiles, C., & Crutcher, R. 2005, in Cosmic Magnetic Fields, ed. R. Wielebinski, & R. Beck Vol. 664, (Berlin: Springer), 137
ADS

↑

Hennebelle, P., & Iffrig, O. 2014, A&A, 570, A81
Crossref ADS

↑

Hirschmann, M., Dolag, K., Saro, A., et al. 2014, MNRAS, 442, 2304
Crossref ADS

↑

Hopkins, P. F., Chan, T. K., Garrison-Kimmel, S., et al. 2020d, MNRAS, 492, 3465
Crossref ADS

↑

Hopkins, P. F., Chan, T. K., Ji, S., et al. 2020a, MNRAS, 501, 3640
Crossref ADS

↑

Hopkins, P. F., Chan, T. K., Squire, J., et al. 2020b, MNRAS, 501, 3663
Crossref ADS

↑

Hopkins, P. F., Squire, J., Chan, T. K., et al. 2020c, MNRAS, 501, 4184
Crossref ADS

↑

Hu, Y., Lazarian, A., Li, Y., Zhuravleva, I., & Gendron-Marsolais, M.-L. 2020, ApJ, 901, 162
Crossref ADS

↑

Hu, Y., Yuen, K. H., Lazarian, V., et al. 2019, NatAs, 3, 776
Crossref ADS

↑

Iapichino, L., & Brüggen, M. 2012, MNRAS, 423, 2781
Crossref ADS

↑

Iapichino, L., Federrath, C., & Klessen, R. S. 2017, MNRAS, 469, 3641
Crossref ADS

↑

Iapichino, L., Viel, M., & Borgani, S. 2013, MNRAS, 432, 2529
Crossref ADS

↑

Jones, T. J., Kim, J.-A., Dowell, C. D., et al. 2020, AJ, 160, 167
Crossref ADS

↑

Jubelgas, M., Springel, V., & Dolag, K. 2004, MNRAS, 351, 423
Crossref ADS

↑

Kannan, R., Vogelsberger, M., Pfrommer, C., et al. 2017, ApJL, 837, L18
Crossref ADS

↑

Kazantsev, A. P. 1968, JETP, 26, 1031
ADS

↑

Kazantsev, A. P., Ruzmaikin, A. A., & Sokolov, D. D. 1985, ZhETF, 88, 487
ADS

↑

Kolmogorov, A. 1941, DoSSR, 30, 301
ADS

↑

Komarov, S. V., Churazov, E. M., Schekochihin, A. A., & ZuHone, J. A. 2014, MNRAS, 440, 1153
Crossref ADS

↑

Kotarba, H., Lesch, H., Dolag, K., et al. 2009, MNRAS, 397, 733
Crossref ADS

↑

Kotarba, H., Lesch, H., Dolag, K., et al. 2011, MNRAS, 415, 3189
Crossref ADS

↑

Kraichnan, R. H., & Nagarajan, S. 1967, PhFl, 10, 859
Crossref ADS

↑

Krause, M., Alexander, P., Bolton, R., et al. 2009, MNRAS, 400, 646
Crossref ADS

↑

Kulsrud, R. M., & Anderson, S. W. 1992, ApJ, 396, 606
Crossref ADS

↑

Kulsrud, R. M., Cen, R., Ostriker, J. P., & Ryu, D. 1997, ApJ, 480, 481
Crossref ADS

↑

Kulsrud, R. M., & Zweibel, E. G. 2008, RPPh, 71, 046901
Crossref ADS

↑

Lacki, B. C., & Beck, R. 2013, MNRAS, 430, 3171
Crossref ADS

↑

Lesch, H., & Hanasz, M. 2003, A&A, 401, 809
Crossref ADS

↑

Lyskova, N., Churazov, E., Zhang, C., et al. 2019, MNRAS, 485, 2922
Crossref ADS

↑

Maier, A., Iapichino, L., Schmidt, W., & Niemeyer, J. C. 2009, ApJ, 707, 40
Crossref ADS

↑

Marinacci, F., & Vogelsberger, M. 2016, MNRAS, 456, L69
Crossref ADS

↑

Marinacci, F., Vogelsberger, M., Mocz, P., & Pakmor, R. 2015, MNRAS, 453, 3999
Crossref ADS

↑

Martin-Alvarez, S., Devriendt, J., Slyz, A., & Teyssier, R. 2018, MNRAS, 479, 3343
Crossref ADS

↑

Martin-Alvarez, S., Slyz, A., Devriendt, J., & Gómez-Guijarro, C. 2020, MNRAS, 495, 4475
Crossref ADS

↑

Matarrese, S., Mollerach, S., Notari, A., & Riotto, A. 2005, PhRvD, 71, 043502
Crossref ADS

↑

McCourt, M., Sharma, P., Quataert, E., & Parrish, I. J. 2012, MNRAS, 419, 3319
Crossref ADS

↑

Miniati, F., Ryu, D., Kang, H., & Jones, T. W. 2001, ApJ, 559, 59
Crossref ADS

↑

Murgia, M., Govoni, F., Feretti, L., et al. 2004, A&A, 424, 429
Crossref ADS

↑

Ota, N., Nagayoshi, K., Pratt, G. W., et al. 2014, A&A, 562, A60
Crossref ADS

↑

Pakmor, R., Gómez, F. A., Grand, R. J. J., et al. 2017, MNRAS, 469, 3185
Crossref ADS

↑

Pakmor, R., & Springel, V. 2013, MNRAS, 432, 176
Crossref ADS

↑

Pakmor, R., van de Voort, F., Bieri, R., et al. 2020, MNRAS, 498, 3125
Crossref ADS

↑

Parrish, I. J., Quataert, E., & Sharma, P. 2009, ApJ, 703, 96
Crossref ADS

↑

Petkova, M., & Springel, V. 2009, MNRAS, 396, 1383
Crossref ADS

↑

Porter, D. H., Jones, T. W., & Ryu, D. 2015, ApJ, 810, 93
Crossref ADS

↑

Powell, K. G., Roe, P. L., Linde, T. J., Gombosi, T. I., & De Zeeuw, D. L. 1999, JCoPh, 154, 284
Crossref ADS

↑

Price, D. J. 2012, JCoPh, 231, 759
Crossref ADS

↑

Rasera, Y., & Chandran, B. 2008, ApJ, 685, 105
Crossref ADS

↑

Rebusco, P., Churazov, E., Böhringer, H., & Forman, W. 2006, MNRAS, 372, 1840
Crossref ADS

↑

Rees, M. J. 1987, QJRAS, 28, 197
ADS

↑

Rees, M. J. 1994, in NATO ASI Ser., Ser. C, ed. D. Lynden-Bell Vol. 422, (Dordrecht: Kluwer), 155
ADS

↑

Rees, M. J. 2005, in Cosmic Magnetic Fields, ed. R. Wielebinski, & R. Beck Vol. 664, (Berlin: Springer), 1
ADS

↑

Rees, M. J. 2006, AN, 327, 395
Crossref ADS

↑

Rephaeli, Y. 1979, ApJ, 227, 364
Crossref ADS

↑

Rephaeli, Y., Ulmer, M., & Gruber, D. 1994, ApJ, 429, 554
Crossref ADS

↑

Rieder, M., & Teyssier, R. 2016, MNRAS, 457, 1722
Crossref ADS

↑

Rieder, M., & Teyssier, R. 2017a, MNRAS, 471, 2674
Crossref ADS

↑

Rieder, M., & Teyssier, R. 2017b, MNRAS, 472, 4368
Crossref ADS

↑

Robishaw, T., Quataert, E., & Heiles, C. 2008, ApJ, 680, 981
Crossref ADS

↑

Roh, S., Ryu, D., Kang, H., Ha, S., & Jang, H. 2019, ApJ, 883, 138
Crossref ADS

↑

Rossetti, M., & Molendi, S. 2004, A&A, 414, L41
Crossref ADS

↑

Ruderman, M. S., Oliver, R., Erdélyi, R., Ballester, J. L., & Goossens, M. 2000, A&A, 354, 261
ADS

↑

Ruzmaikin, A., Sokolov, D., & Shukurov, A. 1988, Natur, 336, 341
Crossref ADS

↑

Ryu, D., Kang, H., Cho, J., & Das, S. 2008, Sci, 320, 909
Crossref ADS

↑

Sarazin, C. L., & Kempner, J. C. 2000, ApJ, 533, 73
Crossref ADS

↑

Schekochihin, A. A., Cowley, S. C., Kulsrud, R. M., Rosin, M. S., & Heinemann, T. 2008, PhRvL, 100, 081301
Crossref ADS

↑

Schekochihin, A. A., Cowley, S. C., Taylor, S. F., Maron, J. L., & McWilliams, J. C. 2004, ApJ, 612, 276
Crossref ADS

↑

Schlickeiser, R., Sievers, A., & Thiemann, H. 1987, A&A, 182, 21
ADS

↑

Schober, J., Schleicher, D. R. G., Federrath, C., Bovino, S., & Klessen, R. S. 2015, PhRvE, 92, 023010
Crossref ADS

↑

Schuecker, P., Finoguenov, A., Miniati, F., Böhringer, H., & Briel, U. G. 2004, A&A, 426, 387
Crossref ADS

↑

Sharma, P., Parrish, I. J., & Quataert, E. 2010, ApJ, 720, 652
Crossref ADS

↑

Spitzer, L. J. 1956, ApJ, 124, 20
Crossref ADS

↑

Spitzer, L., & Härm, R. 1953, PhRv, 89, 977
Crossref ADS

↑

Springel, V. 2005, MNRAS, 364, 1105
Crossref ADS

↑

Springel, V., White, S. D. M., Tormen, G., & Kauffmann, G. 2001, MNRAS, 328, 726
Crossref ADS

↑

Steinwandel, U. P., Beck, M. C., Arth, A., et al. 2019, MNRAS, 483, 1008
Crossref ADS

↑

Steinwandel, U. P., Dolag, K., Lesch, H., et al. 2020a, MNRAS, 494, 4393
Crossref ADS

↑

Steinwandel, U. P., Moster, B. P., Naab, T., Hu, C.-Y., & Walch, S. 2020b, MNRAS, 495, 1035
Crossref ADS

↑

Strong, A. W., Moskalenko, I. V., & Ptuskin, V. S. 2007, ARNPS, 57, 285
Crossref ADS

↑

Subramanian, K., & Barrow, J. D. 2002, MNRAS, 335, L57
Crossref ADS

↑

Sugawara, C., Takizawa, M., & Nakazawa, K. 2009, PASJ, 61, 1293
Crossref ADS

↑

Sullivan, C. H., Fissel, L. M., King, P. K., et al. 2021, MNRAS, 503, 5006
Crossref ADS

↑

Sur, S., Federrath, C., Schleicher, D. R. G., Banerjee, R., & Klessen, R. S. 2012, MNRAS, 423, 3148
Crossref ADS

↑

Suzuki, K., Ogawa, T., Matsumoto, Y., & Matsumoto, R. 2013, ApJ, 768, 175
Crossref ADS

↑

Tabatabaei, F. S., Krause, M., Fletcher, A., & Beck, R. 2008, A&A, 490, 1005
Crossref ADS

↑

Tormen, G., Bouchet, F. R., & White, S. D. M. 1997, MNRAS, 286, 865
Crossref ADS

↑

Tricco, T. S., Price, D. J., & Federrath, C. 2016, MNRAS, 461, 1260
Crossref ADS

↑

van Weeren, R. J., de Gasperin, F., Akamatsu, H., et al. 2019, SSRv, 215, 16
Crossref ADS

↑

Vazza, F., Brüggen, M., Gheller, C., & Wang, P. 2014, MNRAS, 445, 3706
Crossref ADS

↑

Vazza, F., Brunetti, G., Brüggen, M., & Bonafede, A. 2018, MNRAS, 474, 1672
Crossref ADS

↑

Vazza, F., Dolag, K., Ryu, D., et al. 2011, MNRAS, 418, 960
Crossref ADS

↑

Vogt, C., & Enßlin, T. A. 2003, A&A, 412, 373
Crossref ADS

↑

Wang, P., & Abel, T. 2009, ApJ, 696, 96
Crossref ADS

↑

Watson, W. D., & Wyld, H. W. 2001, ApJL, 558, L55
Crossref ADS

↑

Wendland, H. 2004, Scattered Data Approximation, Vol. 17, (Cambridge: Cambridge Univ. Press)
Crossref

↑

Wendland, H. 1995, AdCoM, 4, 389
Crossref

↑

Wik, D. R., Hornstrup, A., Molendi, S., et al. 2014, ApJ, 792, 48
Crossref ADS

↑

Xu, S., & Lazarian, A. 2020, ApJ, 899, 115
Crossref ADS

↑

Yang, H. Y. K., & Reynolds, C. S. 2016, ApJ, 829, 90
Crossref ADS

↑

Zel’dovich, Y. B. 1965, JETP, 21, 656
ADS

↑

Zel’dovich, Y. B. 1970, SvA, 13, 608
ADS

↑

Zel’dovich, Y. B. 1983, Magnetic Fields in Astrophysics (New York, NY: Gordon and Breach)
ADS

↑

ZuHone, J. A., Kunz, M. W., Markevitch, M., Stone, J. M., & Biffi, V. 2015, ApJ, 798, 90
Crossref ADS

↑

ZuHone, J. A., Markevitch, M., Ruszkowski, M., & Lee, D. 2013, ApJ, 762, 69
Crossref ADS

* Released on August, 19th, 2021.

⁵ Essentially, there are no dynamos that properly work without some form of diffusion

⁶ The inverse mass is sometimes referred to as the mobility when rescaled with the mean free path of the particle.

⁷ The hydrodynamics models are used as a sanity check for the MHD models in terms of mass and the growth of the size of the objects (see Section 5.1).

⁸ https://docs.julialang.org