nIFTy galaxy cluster simulations II: radiative models
We have simulated the formation of a massive galaxy cluster (M = 1.110) in a CDM universe using 10 different codes (RAMSES, 2 incarnations of AREPO and 7 of Gadget), modeling hydrodynamics with full radiative subgrid physics. These codes include Smoothed-Particle Hydrodynamics (SPH), spanning traditional and advanced SPH schemes, adaptive mesh and moving mesh codes. Our goal is to study the consistency between simulated clusters modeled with different radiative physical implementations - such as cooling, star formation and AGN feedback. We compare images of the cluster at , global properties such as mass, and radial profiles of various dynamical and thermodynamical quantities. We find that, with respect to non-radiative simulations, dark matter is more centrally concentrated, the extent not simply depending on the presence/absence of AGN feedback. The scatter in global quantities is substantially higher than for non-radiative runs. Intriguingly, adding radiative physics seems to have washed away the marked code-based differences present in the entropy profile seen for non-radiative simulations in Sembolini et al. (2015): radiative physics + classic SPH can produce entropy cores. Furthermore, the inclusion/absence of AGN feedback is not the dividing line -as in the case of describing the stellar content- for whether a code produces an unrealistic temperature inversion and a falling central entropy profile. However, AGN feedback does strongly affect the overall stellar distribution, limiting the effect of overcooling and reducing sensibly the stellar fraction.
keywords:methods: numerical – galaxies: haloes – cosmology: theory – dark matter
This paper is a continuation of the nIFTy cluster comparison project Sembolini et al. (2015): a study of the latest state-of-the-art hydrodynamical codes using simulated galaxy clusters as a testbed for theories of galaxy formation. Simulations are indispensable tools in the interpretation of astronomical observations of these objects (see for instance Borgani & Kravtsov 2011). Although early -body simulations only modeled the gravitational evolution of collisionless effects of dark matter (White 1976; Fall 1978; Aarseth, Turner & Gott 1979), these were vital for interpreting galaxy surveys and unveiling the cosmic web of the large scale structure of the Universe. The focus of modern simulations (see e.g. Vogelsberger et al. 2014; Schaye et al. 2015) has shifted to modeling galaxy formation in a cosmological context , incorporating the key physical processes that govern galaxy formation and the intra-cluster medium (ICM).
The details of the physical processes that are part and parcel of building a galaxy remain uncertain. Naturally, these processes include the conversion of gas to stars and the feedback of energy and metals from supernovae into the surrounding medium (see e.g. Voit 2005 for a review of the radiative processes which govern the evolution of the baryonic component). Galaxy clusters offer an ideal testbed for the study of these processes and their complex interplay, precisely because their enormous size encompasses a wide range of relevant scales. As mentioned before, the goal of modern simulations is now focused on modeling galaxy formation, incorporating the key physical processes that drive galaxy formation - such as the cooling of a collisional gaseous component (e.g. Pearce et al. 2000; Muanwong et al. 2001; Davé, Katz & Weinberg 2002; Kay et al. 2004; Nagai, Kravtsov & Vikhlinin 2007; Wiersma, Schaye & Smith 2009), the birth of stars from cool overdense gas (e.g. Springel & Hernquist 2003; Schaye & Dalla Vecchia 2008), the growth of black holes (Di Matteo, Springel & Hernquist, 2005), and the injection of energy into the inter-stellar medium by supernovae (e.g. Metzler & Evrard 1994; Borgani et al. 2004; Davé, Oppenheimer & Sivanandam 2008 ;Dalla Vecchia & Schaye 2012) and powerful AGN outflows (e.g. Thacker, Scannapieco & Couchman, 2006; Sijacki et al., 2007; Puchwein, Sijacki & Springel, 2008; Sijacki et al., 2008; Booth & Schaye, 2009; Steinborn et al., 2015). These processes span an enormous dynamic range, both spatial and temporal, from the sub-pc scales of black hole growth to the accretion of gas on Mpc scales from the cosmic web.
One of the main issues with radiative simulations of galaxy clusters is that they tend to convert a large fraction of gas into stars. Observationally, only 10-15 per cent of the baryon component of clusters is expected to be in the stellar phase (Gonzalez, Zaritsky & Zabludoff, 2007), but radiative runs which only include stellar feedback are affected by overcooling and usually convert too large a fraction of the gas (above 30 per cent) inside the cluster virial radius into stars (Borgani & Kravtsov, 2011).
Recent work on hydrodyamic simulations has identified AGN feedback as a suitable candidate for overcoming this problem (e.g. Puchwein, Sijacki & Springel, 2008; Puchwein et al., 2010; Fabjan et al., 2010; McCarthy et al., 2010; Battaglia et al., 2012; Martizzi et al., 2012; Planelles et al., 2013; Pike et al., 2014; Le Brun et al., 2014). Heating from AGN occurs via the release of energy during accretion of the ICM gas onto a supermassive black hole hosted by the central cluster galaxy: this energy is sufficiently high to remove gas from the inner regions of clusters. At the same time, AGN heating may also be able to explain the lack of gas in the central region of dynamically relaxed clusters (the “cool core” clusters). Pre-ejection of gas by AGN in the high-redshift progenitors of present-day clusters may also be crucial (McCarthy et al., 2011).
This is where the nIFTy cluster comparison project comes in, building on a long history of important comparison studies of simulated clusters (e.g. the Santa Barbara project, Frenk et al. 1999, hereafter SB99) as well as galaxies (e.g. the Aquila project - Scannapieco et al. 2012 - and the AGORA project - Kim et al. 2014). All codes and subgrid modules attempt to model the key processes of galaxy formation. In our first paper, Sembolini et al. (2015) (hereafter S15), we addressed a well known issue, first highlighted in SB99: mesh-based and traditional SPH codes produced galaxy cluster entropy profiles that were not in agreement. Grid based codes displayed a constant entropy core whereas traditional SPH codes produces profiles that continued to fall all the way towards the centre. The latter behavior was due to the artificial surface tension and the associated lack of multiphase fluid mixing in classic SPH (e.g. Wadsley, Veeravalli & Couchman, 2008; Mitchell et al., 2009). Modern SPH codes attempted to address the lack of mixing through a variety of means: artificial conduction (Price, 2008; Valdarnini, 2012) and pressure-entropy formulations (Ritchie & Thomas, 2001; Saitoh & Makino, 2013; Hopkins, 2013). In S15, we clearly showed that modern SPH is able to create clusters with flat entropy cores that are indistinguishable from those generated by mesh-based codes.
Here we tackle the subgrid physics implemented in a variety of state-of-the-art codes. We extend the analysis presented in S15 by performing simulations of the same cluster with full physics runs where codes have radiative mechanisms describing gas cooling, star formation, supernova feedback, black hole accretion and AGN feedback. We used 10 different codes (RAMSES, 2 incarnations of AREPO, 7 of Gadget), allowing each method to choose their favorite radiative processes modeled by subgrid physics. This allows us to study how the different mechanisms, especially star formation and AGN feedback, influence the properties of simulated clusters. We examine the overall cluster environment and we focus our analysis on revisiting the gas entropy profiles.
The rest of this paper is organized as follows: in Section 2 we briefly describe the codes used and the subgrid physics adopted by each code along with a brief description of the data set. We then discuss our results in Sections 3-5: starting with an overview of the bulk properties of the cluster and the effect of radiative physics (Section 3); followed by the dark matter distribution (Section 4); we continue our analysis by studying the baryon distribution (Section 5): in Section 5.1 we describe key properties of the gas component such as the temperature, entropy and gas fraction, concluding our analysis by presenting the code-to-code differences in the distribution of stars (Section 5.2). We report our conclusions in Section 6.
2 The Codes
The initial nIFTy cluster comparison project, as presented in Sembolini et al. (2015), included 13 codes – RAMSES, ART, AREPO, Hydra and 9 variants of the Gadget code. In this study, we consider the subset of these codes in which full radiative subgrid physics has been included. A comprehensive summary of the approach taken to solve the hydrodynamic equations in each of these codes can be found in S15; here we provide a brief recap of this summary, with a focus on a description of the sub-grid physics implemented in each code. Table 1 lists the codes included in this work and their basic characteristics (the definition of modern and classic SPH codes, as well as that of grid-based and moving-mesh codes, is provided in Section 2 of S15).
|Grid-based||RAMSES||Y||Y||Ramses-AGN||Teyssier et al. (2011)|
|Moving mesh||AREPO||Y||Y||Arepo-IL||Vogelsberger et al. (2013, 2014)|
|G3-PESPH||Y||N||Huang et al. (in prep.)|
|G3-Magneticum||Y||Y||Hirschmann et al. (2014)|
|Classic SPH||G3-Music||Y||N||G3-Music||Sembolini et al. (2013)|
|G2-MusicPI||Piontek & Steinmetz (2011)|
|G3-OWLS||Y||Y||Schaye et al. (2010)|
|G2-X||Y||Y||Pike et al. (2014)|
2.1 Mesh-based Codes
Ramses (Teyssier, Perret)
Ramses is an adaptive mesh refinement code. For fluid dynamics a directionally unsplit, second order Godunov scheme with the HLLC Riemann solver is used. The N-body solver is an adaptive particle mesh code, for which the Poisson equation is solved using the multi-grid technique. The grid is adaptively refined on a cell-by-cell basis, following a quasi-Lagrangian refinement strategy whereby a cell is refined into 8 smaller new cells if its dark matter or baryonic mass grows by more than a factor of eight. Time integration is performed using an adaptive, level-by-level, time stepping strategy. Parallel computing is based on the MPI library, with a domain decomposition set by the Peano-Hilbert space filling curve.
Cooling & Heating: Gas cooling and heating is performed assuming coronal equilibrium with a modification of the Haardt & Madau (1996) UV background and a self-shielding recipe based on Aubert & Teyssier (2010), with an exponential cut-off of the radiation flux with critical density H/cm. All Hydrogen and Helium cooling and heating processes are included following Katz, Weinberg & Hernquist (1996). Metal cooling is added using the Sutherland & Dopita (1993) metal-only cooling function at solar metallicity, multiplied by the local metallicity of the gas in solar units. In this particular project, we use also a temperature floor K to prevent spurious fragmentation of our relatively poorly resolved galactic discs.
Star Formation: Star formation is implemented as a stochastic process using a local Schmidt law, as in Rasera & Teyssier (2006). The density threshold for star formation was set to H/cc, and the local star formation efficiency per gas free fall time was set to 5 per cent.
Stellar Population Properties & Chemistry: Each star particle is treated as a single stellar population (SSP) with a Chabrier (2003) IMF. Mass and metal return to the gas phase by core collapse supernovae only. A single average metal specie is followed during this process and advected in the gas as a passive scalar, to be used as an indicator of the gas metallicity in the cooling function.
Stellar Feedback: In this project, no feedback processes related to the stellar population are used.
SMBH Growth & AGN Feedback: The formation of SMBH particle is allowed using the sink particle technique described in Teyssier et al. (2011). When the gas density is larger than the star formation density threshold, a boost in the Bondi accretion rate is allowed, using the boost function proposed by Booth & Schaye (2009). The SMBH accretion rate is never allowed to exceed the instantaneous Eddington limit. SMBH particles are evolved using a direct gravity solver, to obtain a more accurate treatment of their orbital evolution. SMBH particles more massive than are allowed to merge if their relative velocity is smaller than their pairwise scale velocity. Less massive SMBH particles, on the other hand, are merged as soon as they fall within 4 cells from another SMBH particle. The AGN feedback used is a simple thermal energy dump with of specific energy, multiplied by the instantaneous SMBH accretion rate. .
Here we use two different versions of AREPO: one including AGN feedback (Arepo-IL) and one not including it (Arepo-SH).
AREPO uses a Godunov scheme on an unstructured moving Voronoi mesh; mesh cells move (roughly) with the fluid. The main difference between AREPO and traditional Eulerian AMR codes is that AREPO is almost Lagrangian and Galilean invariant by construction. The main difference between AREPO and SPH codes (see next subsection) is that the hydrodynamic equations are solved with a finite-volume Godunov scheme. The version of AREPO used in this study conserves total energy in the Godunov scheme, rather than the entropy-energy formalism described in Springel (2010). Detailed descriptions of the galaxy formation models implemented in AREPO can be found in Vogelsberger et al. (2013) and Vogelsberger et al. (2014), but the key features can be summarized as follows (hereafter we describe the features of Arepo-IL, the radiative models used for Arepo-SH are the same than G3-Music, and are listed later in this section).
Cooling & Heating: Gas cooling takes the metal abundance into account. The metal cooling rate is computed for solar composition gas and scaled to the total metallicity of the cell. Photoionization and photoheating are followed based on the homogeneous UV background model of Faucher-Giguère et al. (2009) and the self-shielding prescription of Rahmati et al. (2013). In addition to the homogeneous UV background, the ionizing UV emission of nearby AGN is taken into account.
Star Formation: The formation of stars is followed with a multi- phase model of the interstellar medium which is based on Springel & Hernquist (2003) (hereafter SH03) but includes a modified effective equation of state above the star formation threshold, i.e. above a hydrogen number density of 0.13 cm
Stellar Population Properties & Chemistry: Each star particle is treated as a single stellar population (SSP) with a Chabrier (2003) IMF. Mass and metal return to the gas phase by AGB stars, core collapse supernovae and Type Ia supernovae is taken into account. Nine elements are followed during this process (H, He, C, N, O, Ne, Mg, Si, Fe).
Stellar Feedback: Feedback by core collapse supernovae is implicitly invoked by the multiphase star formation model. In addition, we include a kinetic wind model in which the wind velocity scales with the local dark matter velocity dispersion ( 3.7) The mass-loading is determined by the available energy which is assumed to be 1.0910erg per core collapse supernova. Wind particles are decoupled from the hydrodynamics until they fall below a specific density threshold or exceed a maximum travel time. This ensures that they can escape form the dense interstellar medium.
SMBH Growth & AGN Feedback: SMBHs are treated as collisionless sink particles. Particles with a mass of 10M are seeded into halos once they exceed a mass of 510M. The BHs subsequently grow by Bondi-Hoyle accretion with a boost factor of = 100. The Eddington limit on the accretion rate is enforced in addition. AGN are assumed to be in the quasar mode for accretion rates larger than 5 per cent of the Eddington rate. In this case 1 per cent of the accreted rest mass energy is thermally injected into nearby gas. For accretion rates smaller than 5 per cent of the Eddington rate, AGN are in the radio mode in which 7 per cent of the accreted rest mass energy is thermally injected into spherical bubbles (similar to Sijacki et al. 2007). Full details of the black hole model are given in Sijacki et al. (2014).
2.2 SPH Codes
Gadget3-X (Murante, Borgani, Beck)
G3-X code is a development of the non-public Gadget3 code. It includes an improved SPH scheme, described in Beck et al. (2015). Main changes with respect to the standard Gadget3 hydro are: (i) an artificial conduction term that largely improves the SPH capability of following gas-dynamical instabilities and mixing processes; (ii) a higher-order kernel (Wendland C4) to better describe discontinuities and reduce clumpiness instability; (iii) a time-dependent artificial viscosity term to minimize viscosity away from shock regions. Both pure hydrodynamical and hydro/gravitational tests on the performance of our improved SPH are presented in Beck et al. (2015).
Cooling & Heating: Gas cooling is computed for an optically thin gas and takes into account the contribution of metals, using the procedure of Wiersma, Schaye & Smith (2009), while a uniform UV background is included following the procedure of Haardt & Madau (2001).
Star Formation: Star formation is implemented as in Tornatore et al. (2007) , and follows the star formation algorithm of SH03 – gas particles above a given density threshold are treated as multi-phase. The effective model of SH03 describes a self-regulated, equilibrium inter-stellar medium and provides a star formation rate that depends upon the gas density only, given the model parameters.
Stellar Population Properties & Chemistry: Each star particle is considered to be a single stellar population (SSP). We follow the evolution of each SSP, assuming the Chabrier (2003) IMF. We account for metals produced in the SNeIa, SNeII and AGB phases, and follow 15 chemical species. Star particles are stochastically spawned from parent gas particles as in SH03, and get their chemical composition of their parent gas. Stellar lifetimes are from Padovani & Matteucci (1993); metal yields from Woosley & Weaver (1995) for SNeII, Thielemann et al. (2003) for SNeIa, and van den Hoek & Groenewegen (1997) for AGB stars.
Stellar Feedback: SNeII release energy into their surroundings, but this only sets the hot gas phase temperature and, as a consequence, the average SPH temperature of gas particles. Supernova feedback is therefore modeled as kinetic and the prescription of SH03 is followed (i.e. energy-driven scheme with a fixed wind velocity of 350 km/s, wind particles decoupled from surrounding gas for a period of 30 Myr or until ambient gas density drops below 0.5 times the multiphase density threshold).
SMBH Growth & AGN Feedback: AGN feedback, follows the model described in Steinborn et al. (2015). Nevertheless, while this model includes a Bondi-Hoyle like gas accretion (Eddington limited) onto SMBH, distinguishing the cold and the hot component (their Eq. 19), here we only consider the cold accretion, using a fudge-factor in the Bondi-Hoyle formula. In other words, . The radiative efficiency is variable, and it is evaluated using the model of (Churazov et al., 2005). Such a model outputs separately the AGN mechanical and radiative power as a function of the SMBH mass and the accretion rate; however, here we sum up these powers and give the resulting energy to the surrounding gas, in form of purely thermal energy. We set the efficiency of AGN feedback/gas coupling to .
We tuned the parameters of our new hydro scheme using the tests presented in Beck et al. (2015), and those of the AGN model for reproducing observational scaling relations between SMBH mass and stellar mass of the host galaxies. We note that we did not make any attempt to tune parameters to reproduce any of the observational properties of the ICM. First results on the application of this code to simulations of galaxy clusters, including the reproduction of the Cool Core/ Non-Cool Core dichotomy, can be found in Rasia et al. (2015).
A black hole (BH) of mass 510 is seeded at the centre of each friends-of-friends (FoF) group whose mass exceeds 2.510 and which does not already contain a BH.
Gadget3-PESPH (February, Davé, Katz, Huang)
This version of Gadget uses the pressure-entropy SPH formulation of Hopkins (2013) with a 128 neighbour HOCTS(n=5) kernel and the time-dependent artificial viscosity scheme of Morris & Monaghan (1997).
Cooling & Heating: Radiative cooling using primordial abundances is modeled as described in Katz, Weinberg & Hernquist (1996), with additional cooling from metal lines assuming photo-ionization equilibrium following Wiersma, Schaye & Smith (2009). A Haardt & Madau (2001) uniform ionizing UV background is assumed.
Star Formation: Star formation follows the approach set out in SH03, where a gas particle above a density threshold of cm is modeled as a fraction of cold clouds embedded in a warm ionized medium, following McKee & Ostriker (1977). The star formation rate obeys the Schmidt (1959) law and is proportional to , with the star formation timescale scaled to match the =0 Kennicutt (1998) relation. In addition, the heuristic model of Rafieferantsoa et al. (2014), tuned to reproduce the exponential truncation of the stellar mass function, is used to quench star formation in massive galaxies. A quenching probability , which depends on the velocity dispersion of the galaxy, determines whether or not star formation is stopped in a given galaxy; if it is stopped, each gas particle eligible for star formation first has its quenching probability assessed, and if it is selected for quenching then it is heated to 50 times the galaxy virial temperature, which unbinds it from the galaxy.
Stellar Population Properties & Chemistry: Each star particle is treated as a single stellar population with a Chabrier (2003) IMF throughout. Metal enrichment from SNeIa, SNeII and AGB stars are tracked, while 4 elements – C, O, Si and Fe – are also tracked individually, as described by Oppenheimer & Davé (2008).
Stellar Feedback: Supernova feedback is assumed to drive galactic outflows, which are implemented using a Monte Carlo approach analogous to that used in the star formation prescription. Outflows are directly tied to the star formation rate, using the relation SFR, where is the outflow mass loading factor. The probability for a gas particle to spawn a star particle is calculated from the subgrid model described above, and the probability to be launched in a wind is times the star formation probability. If the particle is selected to be launched, it is given a velocity boost of in the direction of , where and are the particle instantaneous velocity and acceleration, respectively. This is a highly constrained heuristic model for galactic outflows, described in detail in Davé et al. (2013), which utilizes outflows scalings expected for momentum-driven winds in sizable galaxies (km s), and energy-driven scalings in dwarf galaxies. In particular, the mass loading factor (i.e. the mass outflow rate in units of the star formation rate) is for galaxies with velocity dispersion km s , and for km s.
SMBH Growth & AGN Feedback: These processes are not included.
G3-Magneticum is an advanced version of Gadget3. In this version, a higher order kernel based on the bias-corrected, sixth-order Wendland kernel (Dehnen & Aly, 2012) with 295 neighbors is included. The code also incorporates a low viscosity scheme to track turbulence as original described in Dolag et al. (2005) with improvements following Beck et al. (2015). Gradients are computed with high-order scheme (Price, 2012) and thermal conduction is modeled isotropically at 1/20th of the Spitzer rate (Dolag et al., 2004). The simulation is run with a time-step limiting particle wake-up algorithm (Pakmor et al., 2012). The models adopted for cooling, star formation and stellar feedback are the same that in G3-X, but with different parameters.
Cooling & Heating: The simulation allows for radiative cooling according to (Wiersma, Schaye & Smith, 2009) and heating from a uniform time-dependent ultraviolet background (Haardt & Madau, 2001). The contributions to cooling from each one of 11 elements (H, He, C, N, O, Ne, Mg, Si, S, Ca, Fe) have been pre-computed using the publicly available CLOUDY photoionization code (Ferland et al., 1998) for an optically thin gas in (photo-)ionization equilibrium.
Star Formation: We model the interstellar medium (ISM) by using a subresolution model for the multiphase ISM of Springel & Hernquist (2003). In this model, the ISM is treated as a two-phase medium, in which clouds of cold gas form by cooling of hot gas, and are embedded in the hot gas phase assuming pressure equilibrium whenever gas particles are above a given threshold density.
Stellar Population Properties & Chemistry: We include a detailed model of chemical evolution according to Tornatore et al. (2007). Metals are produced by SNII, by supernovae type Ia (SNIa) and by intermediate and low-mass stars in the asymptotic giant branch (AGB). Metals and energy are released by stars of different masses, by properly accounting for mass-dependent life-times (with a lifetime function according to Padovani & Matteucci 1993), the metallicity-dependent stellar yields by Woosley & Weaver (1995) for SNII, the yields by van den Hoek & Groenewegen (1997) for AGB stars, and the yields by Thielemann et al. (2003) for SNIa. Stars of different masses are initially distributed according to a Chabrier initial mass function (IMF; Chabrier 2003).
Stellar Feedback: The hot gas within the multiphase model describing the ISM is heated by supernovae and can evaporate the cold clouds. A certain fraction of massive stars (10 per cent) is assumed to explode as supernovae type II (SNII). The released energy by SNII (10 erg) triggers galactic winds with a mass loading rate proportional to the star formation rate (SFR) with a resulting wind velocity of = 350 km/s.
SMBH Growth & AGN Feedback: Our simulations include prescriptions for the growth of black holes and the feedback from active galactic nuclei (AGN) based on the model of Springel et al. (2005a) and Di Matteo, Springel & Hernquist (2005) with the same modifications as in Fabjan et al. (2010) and some new, minor changes as described below. The accretion onto black holes and the associated feedback adopts a sub-resolution model. Black holes can grow in mass by either accreting gas from their environments, or merging with other black holes. The gas accretion rate is estimated by the Bondi-Hoyle Lyttleton approximation, (Hoyle & Lyttleton 1939; Bondi & Hoyle 1944; Bondi 1952). The black hole accretion is always limited to the Eddington rate and a characteristic boost factor of 100 is applied as only the accretion to large scale is captured. Unlike in Springel et al. (2005a), in which a selected gas particle contributes to accretion with all its mass, we include the possibility for a gas particle to accrete only with a fraction (1/4) of its original mass. A fraction r = 0.1 of the accreted mass is converted into energy, and a fraction f = 0.1 of this energy is then thermally coupled with gas within the smoothing length of the BH, weighted using the same SPH kernel used for the hydrodynamics. Following Sijacki et al. (2007), when the accretion rate drops below a given threshold, it is assumed that there is a transition from a ÒquasarÓ mode to a ÒradioÓ mode of AGN feedback, and the feedback efficiency is enhanced by a factor of 4. In contrast to Springel et al. (2005a), we modify the mass growth of the BH by taking into account the feedback, e.g., . Other more technical modifications on the BH dynamics with respect to the original implementation have been included. We refer the reader to Dolag et al. (2015) & Hirschmann et al. (2014) for more details, where we also demonstrate that the bulk properties of the AGN population within the simulation are quite similar to the observed AGN properties.
Gadget3-OWLS (McCarthy, Schaye)
This is a heavily modified version of Gadget3 using a classic entropy-conserving SPH formulation with a 40 neighbor M3 kernel.
Cooling & Heating: Radiative cooling rates for the gas are computed on an element-by-element basis by interpolating within pre-computed tables (generated with the CLOUDY code; cf. Ferland et al. 2013) that contain cooling rates as a function of density, temperature and redshift calculated in the presence of the cosmic microwave background and photoionization from a Haardt & Madau (2001) ionizing UV/X-ray background (further details in Wiersma, Schaye & Smith, 2009).
Star Formation: Star formation follows the prescription of Schaye & Dalla Vecchia (2008) – gas with densities exceeding the critical density for the onset of the thermogravitational instability is expected to be multiphase and to form stars (Schaye, 2004). Because the simulations lack both the physics and numerical resolution to model the cold interstellar gas phase, an effective equation of state (EOS) is imposed with pressure for densities where . As described in Schaye & Dalla Vecchia (2008), gas on the effective EOS is allowed to form stars at a pressure-dependent rate that reproduces the observed Kennicutt-Schmidt law (Kennicutt, 1998) by construction.
Stellar Population Properties & Chemistry: The ejection of metals by massive- (SNeII and stellar winds) and intermediate-mass stars (SNeIa, AGB stars) is included following the prescription of Wiersma et al. (2009). A set of 11 individual elements are followed (H, He, C, Ca, N, O, Ne, Mg, S, Si and Fe), which represent all the important species for computing radiative cooling rates.
Stellar Feedback: Feedback is modeled as a kinetic wind (Dalla Vecchia & Schaye, 2008) with a wind velocity and a mass loading , which corresponds to using approximately 40 per cent of the total energy available from SNe for the adopted Chabrier (2003) IMF. This choice of parameters results in a good match to the peak of the SFR history of the universe (Schaye et al., 2010).
SMBH Growth & AGN Feedback: Each black hole can grow either via mergers with other black holes within the softening length or via Eddington-limited gas accretion, the rate of which is calculated using the Bondi-Hoyle formula with a modified efficiency, setting as in Booth & Schaye (2009). The black hole is forced to sit on the local potential minimum, to suppress spurious gravitational scattering (Springel et al., 2005b). Feedback is done by storing up the accretion energy (assuming , ) until at least one particle can be heated to a fixed temperature of (Booth & Schaye, 2009). A friends-of-friends algorithm is run on the fly and FOF haloes with at least 100 dark matter particles (and that do not yet have a black hole particle) are seeded with a black hole particle. The initial mass of this particle is set to 10 times the (initial) gas mass.
Gadget2-X (Kay, Newton)
This is a modified version of the original Gadget2 Tree-PM code that uses the classic entropy-conserving SPH formulation with a 40 neighbor M3 kernel. A detailed description of the code can be found in Pike et al. (2014), but its key features can be summarized as follows.
Cooling & Heating: Cooling follows the prescription of Thomas & Couchman (1992) – a gas particle is assumed to radiate isochorically over the duration of its timestep. Collisional ionization equilibrium is assumed and the cooling functions of Sutherland & Dopita (1993) are used, with the metallicity =0 to ignore the increase in cooling rate due to heavy elements. Photo-heating rates are not included but the gas is heated to a minimum at and .
Star Formation: Star formation follows the method of Schaye & Dalla Vecchia (2008); it assumes an equation of state for the gas with , with an effective adiabatic index of for constant Jeans mass. Gas is converted to stars at a rate given by the Kennicutt-Schmidt relation (Schmidt, 1959; Kennicutt, 1998), assuming a disc mass fraction . The conversion is done stochastically on a particle-by-particle basis so the gas and star particles have the same mass.
Stellar Population Properties & Chemistry: Each star particle is assumed to be a single stellar population with a Salpeter (1955) IMF.
Stellar Feedback: A prompt thermal Type II SNe feedback model is used. This assumes that a fixed number, , of gas particles are heated to a fixed temperature, , with values of and =10 chosen to match observed hot gas and star fractions (cf. Pike et al., 2014). Heated gas is allowed to interact hydrodynamically with its surroundings and radiate.
SMBH Growth & AGN Feedback: A variation on the Booth & Schaye (2009) model is used. Black holes are seeded in friends-of-friends (FOF) haloes with more than 50 particles at , at the position of the most bound star or gas particle, which is replaced with a black hole particle. The gravitational mass of the replaced particle is unchanged but an internal mass of is adopted for the calculation of feedback. Each black hole can grow either via mergers with other black holes within the softening length or via Eddington-limited gas accretion, the rate of which is calculated using the Bondi-Hoyle formula with a modified efficiency, setting as in Booth & Schaye (2009). The black hole is forced to sit on the local potential minimum, to suppress spurious gravitational scattering. Feedback is done by storing up the accretion energy (assuming , ) until at least one particle can be heated to a fixed temperature of . This high temperature was chosen for high-mass clusters to match their observed pressure profiles – a lower temperature causes too much gas to accumulate in cluster cores because there is insufficient entropy to escape to larger radius).
Gadget3-MUSIC (Yepes, Sembolini)
This is the original code adopted for MUSIC-2 dataset (Sembolini et al., 2013), simulated using a modified version of the Gadget3 Tree-PM code that uses classic entropy-conserving SPH formulation with a 40 neighbor M3 kernel. The basic SH03 model was used, the key features of which can be summarized as follows. In this work we also present Gadget2-MUSIC, an alternative version of MUSIC performed using the radiative feedbacks described in Piontek & Steinmetz (2011) (G2-MusicPI since now on).
Cooling & Heating: Radiative cooling is assumed for a gas of primordial composition, with no metallicity dependence, and the effects of a background homogeneous UV ionizing field is assumed, following Haardt & Madau (2001).
Star Formation: The SH03 model is implemented.
Stellar Population Properties & Chemistry: A Salpeter (1955) IMF is assumed, with a slope of -1.35 and upper and lower mass limits of and respectively.
Stellar Feedback: This has both a thermal and a kinetic mode; thermal feedback evaporates the cold phase within SPH particles and increases the temperature of the hot phase, while kinetic feedback is modeled as a stochastic wind (as in SH03) – gas mass is lost due to galactic winds at a rate , which is proportional to the star formation rate , such that , with . SPH particles near the star forming region will be subjected to enter in the wind in an stochastic way. Those particles impacted upon by the wind will be given an isotropic velocity kick of = 400 km/s and will freely travel without feeling pressure forces up to 20 kpc distance from their original positions
SMBH Growth & AGN Feedback: These processes are not included.
Colour & line style scheme
In all the radial plots below we distinguish codes including AGN feedback from codes which only include stellar feedback. The first group is identified by dashed lines and the second one by solid lines. Each code is identified by a different color. In all the plots, the codes are ordered by decreasing gas fraction at from left to right (or top to bottom).
2.3 The Data
We use zoom simulations of clusters produced with a variety of codes running full physics (FP) models, building upon the dark matter only and non-radiative simulations of S15. The initial conditions for our zoom simulations were drawn from the MUSIC-2 cluster catalog (Sembolini et al. 2013; Sembolini et al. 2014; Biffi et al. 2014)
The mass of a gas element naturally varies in our mesh codes. Star particle masses varies from code to code depending on how many generations of stars a gas element produces and the mass of the gas element being converted into a star particle.
3 Bulk Properties
Before we focus on the various components of our simulated clusters, we analyze the impact that the different subgrid models adopted in full physics simulations (FP) have on the bulk properties of the cluster.
As already mentioned in Section 1, one of the main goals of modern simulations is to give a description of the baryonic (galaxies and ICM) component of clusters which succeeds in reproducing observational results. We therefore start our analysis by testing how the different codes used in this work compare with measurements of the gas and stellar components as provided by observations. We show in Figure 1 the values of as calculated at , the radius enclosing = 500 times the critical density (the gas fraction with respect to the total mass of the cluster) against those of (the star fraction) evaluated at the same overdensity. The green area indicates the range of values allowed by observations; as observational results still do not agree (e.g. Gonzalez et al. 2013 invokes higher gas fractions for massive clusters with respect to previous results, see Section 5 for a more detailed discussion), we set very non-restrictive limits to the extreme permitted values: 0.11 0.174 (the value of the cosmic ratio according to WMAP7) and 0.005 0.03. We see that most of the codes not including AGN feedback show values of the stellar fraction which have been ruled out by observations, although they are able to reproduce the gas content. In this work we do not use an observational approach to estimate baryonic masses (e.g. measuring the gas fractions from synthetic X-rays observations), but we estimate the masses by simply counting the number of particles inside a fixed radius.
Figure 2 shows a selection of global properties calculated within , the radius enclosing 200 times the critical density: radius, mass, mass-weighted gas temperature, gas and stellar fractions, shape parameters (here we report the values of the minor semi-axes, and , normalized to that of the major semi-axis, ) and the one dimensional velocity dispersion, . The first feature is that the scatter in FP simulations is higher than in the non-radiative (NR) case (see S15). The mean values for the total mass, radius, shape (with the exception in this case of Ramses-AGN) and DM velocity dispersion are extremely close to those in the non-radiative runs and still have very low scatter (less than 2 per cent).
More importantly, pronounced differences lie in the baryonic sector. The temperature (4.3 keV, corresponding approximately to 510K) is 20 per cent higher in FP simulations than in NR models (3.7 keV) and has a scatter around 5 per cent compared to that of 2 per cent registered in the NR comparison. The gas fraction is lower than what was found in the non-radiative case (as some of the gas has been converted to stars), especially for the codes which do not include AGN feedback. The overall fractions show significant scatter: and a code-to-code scatter of 30-40 per cent; the discrepancies are more dramatic for the stellar component, where varies between . The total baryon fraction () shows a more moderate scatter (around 10 per cent) and most of the codes show values around 0.16, very close to the cosmic ratio (here we adopt the value -used for our simulations- of 0.174 reported using WMAP7+BAO+SNI data by Komatsu et al. 2011). Ramses-AGN is the outlier, showing a baryon fraction that is slightly larger than the cosmic ratio ( 0.18). Interestingly, we observe a trend in the AGN codes, from Ramses-AGN to G3-OWLS the temperature tends to increase and at the same time, the gas fraction tends to decrease. This may suggest a variation in feedback strength from left-to-right (as more and more gas is expelled, the remaining gas is hotter).
Figure 3 shows how the main global cluster properties reported in Figure 2 changed in full physics simulations with respect to the NR runs reported in S15. The quantities that exhibit less scatter (e.g. mass and radius) are, as expected, also the ones whose values were basically unchanged with respect to the NR models, with differences lower than 1 per cent (only for Ramses-AGN some of these values are 5 per cent higher than its NR version) and scatters between 1 and 3 per cent. The temperature and gas fraction, which depend only on the baryon component and are therefore more affected by radiative processes, exhibit higher differences: as the gas is heated by the different energy injection mechanisms included in the FP simulations, temperatures are on average 10 per cent higher (with the only exception of Ramses-AGN, which registers a temperature a few per cent lower than its NR model) with a scatter of 7 per cent. Furthermore, as part of the baryon component is now converted into stars, the gas fraction is now substantially lower: we find a median value of 15 per cent and a scatter of 13 per cent. On the other hand, the methods with the lowest portion of baryons converted into stars (see Section 5.2), such as Ramses-AGN and G3-X, show a gas fraction very close to the value registered for the corresponding NR version. The total baryon fraction is either almost unaltered or 5-10 per cent lower than in the NR case for almost all the codes.
4 Dark matter
A visual comparison of the density field centered on the cluster at is presented in Figure 4 and density profiles are shown in Figure 5. Although all the codes successfully recover the same object and its main features (e.g. the position of the main subhalo, which in the maps is located at 7 o’clock close to , except for Ramses-AGN, which seems to have a slightly different merger phase), the dark matter distribution differs significantly more than what was found in S15 for the dark matter-only and non-radiative models.
These differences in the dark matter distribution arise in response to the baryons. As baryons cool they can pull in dark matter with an effect similar to adiabatic contraction (Eggen, Lynden-Bell & Sandage 1962; Zel’dovich et al. 1980). This contraction may look surprising at first sight as dark matter dominates the mass budget of the cluster, exceeding baryonic matter by a factor of 6. However, the gravitational field in the central regions of a halo is dominated by stars, which formed from the condensations of cooling baryons. The amount of the contraction was studied for the first time in cosmological simulations by Gnedin et al. (2004) (and recently revisited by Capela, Pshirkov & Tinyakov 2014). These studies indicated that cooling and star formation can produce clusters and galaxies with central dark matter densities that are an order of magnitude higher than analogues in non-radiative runs. Duffy et al. (2010) studied the effects of feedback from star formation and AGN, finding large variations and much less contraction when AGN feedback is included.
Of greater significance is the variety in the dark matter distributions, most easily seen in the radial profiles of Figure 5. The first notable systematic is that codes which exhibit a stronger contraction are those which do not include AGN feedback (G3-Music, G2-MusicPI, Arepo-SH), with the exception of G3-PESPH. These codes have inner regions (kpc) with densities a factor of 2 higher than the other codes. Many studies show that simulations of clusters that lack a physical mechanism to stop the central cooling of the gas are affected by the problem of overcooling (e.g. Suginohara & Ostriker 1998; Lewis et al. 2000; Tornatore et al. 2003; Nagai & Kravtsov 2004). These codes have a notably higher fraction of the baryons in the form of cold gas and stars within the virial radius than inferred from observations, 30-50 per cent vs 10-20 per cent, and are expected to produce more stars (see Section 5.2 for a more detailed discussion). This picture fits with their higher dark matter concentrations.
Codes that include AGN feedback do not have such a pronounced contraction, with dark matter profiles similar to that reported for NR runs (see Figure 2 in S15). The interesting exception noted before is G3-PESPH, which has a profile similar to G2-X and G3-X. Among the AGN codes, Arepo-IL experiences the smallest contraction, a factor of 2 less than the other codes. As the contraction is related to the star formation efficiency, it is no surprise to find that Arepo-IL is one of the codes with the fewest stars (see Figure 13 in Section 5.2).
The profiles not only show systematic differences, the code-to-code scatter in full physics simulations is considerably higher (up to a factor of 5 between the two different versions of AREPO at the center of the halo) than that observed in the DM-only and non-radiative runs (see Figures 1, 2 and A1 of S15), where differences never exceeded 20 per cent. This scatter occurs primarily in the central regions. The cluster outskirts show a scatter of 10 per cent. The large difference between the two different versions of AREPO confirms how the dark matter distribution depends on the subgrid physics adopted, and in particular by how energy is injected into the gas reservoir.
A visual comparison of the gas density field centered on the cluster at is presented in Figure 6. There is a substantial amount of variation in the central gas density, with some methods (Arepo-IL, G3-X, G2-MusicPI, G3-OWLS) having significantly larger extended nuclear regions. Some codes appear to show numerous small dense gas clumps in the cluster outskirts, especially those including AGN feedback: in this case AGN prevents gas from cooling and forming stars, and therefore more gas is left in these substructures. More significantly, we observe that different subgrid physics applied to the same code (AREPO) produces very different gas environments.
Figure 7 allows a visual comparison of stellar density distributions. The projected stellar densities appear to show even more variation. Ramses-AGN, Arepo-IL and G3-X have dense stellar objects whereas Arepo-SH, G3-Music and G2-MusicPI have significantly more extended stellar distribution. G3-OWLS also has an extended intra-cluster stellar halo but also has numerous stellar concentrations. Moreover, features of the gas distribution do not map to features in the stellar distribution, i.e., an extended gas distribution does not necessarily produce an extended stellar distribution. For instance, both G2-X and Ramses-AGN show a very high gas concentration in the core, but the latter produces a much more limited star distribution.
The gas differences seen in Figure 6 are also evident in the radial gas density profiles presented in Figure 8. The code-to-code scatter in the central regions is 40 per cent and decreases in the outskirts of the cluster. The outliers are G3-PESPH, which produces the lowest central density in the core (a factor of times smaller), and Ramses-AGN, which has the highest. In the outskirts the differences among codes are much more contained at overdensities lower than 2500 (although Ramses-AGN, Arepo-IL and G3-X show slightly higher gas densities). Interestingly, we also notice that G3-Music and Arepo-SH, which adopt the same star formation model (SH03), show very similar gas fraction profiles in the outskirts. The scatter is generally higher than in the non-radiative case (see Figure 6 of S15). As anticipated visually by Figure 6, the same hydrodynamics code with different subgrid physics produces different gas distributions (e.g. AREPO). Furthermore, in the behavior of the gas density there is not a clear distinction between grid-based and modern SPH codes on the one hand and classic SPH on another hand as highlighted in the NR case (Figure 6 of S15).
We next show in Figure 9 the radial mass-weighted temperature profiles, defined as:
where and are the mass and temperature of the gas particles/cells. The code-to-code scatter is large, especially at the center of the cluster. G3-OWLS, G3-Magneticum, Ramses-AGN and G2-X show a central temperature inversion, similar to that observed in non-radiative, classic SPH simulations: the inner temperature is times smaller than the peak value, which here is keV. In particular, G3-Magneticum shows a very sharp temperature inversion at R 0.1Mpc: this effect is probably due to overcooling, as a large portion of the gas in the core is converted into a massive gaseous BCG. In contrast, all the other codes display rising profiles going towards the the core, a behavior that is observed in modern SPH and mesh-based non-radiative simulations (see Figure 7 of S15). The typical peak temperature for this cluster in these codes is keV. The outlier amongst the codes with no temperature inversion is Arepo-SH, which has an inner temperature exceeding 20 keV. Intriguingly, pronounced differences between codes including and not including AGN feedback are not visible. It is also interesting that some classic SPH codes (such as G3-Music), which in the non-radiative simulations produce a central temperature inversion, now produce monotonically rising temperature profiles (in agreement with Rasia et al. 2014, which pointed out that radiative processes decrease the tension in temperature profiles between classic SPH and adaptive-mesh codes).
We combine the gas density and temperature to produce the radial gas entropy profiles shown in Figure 10, where we adopt the definition of entropy commonly used in the observational X-ray literature:
where is the number density of free electrons of the gas. We observe that the differences between modern and classic SPH methods that had been displayed for the non-radiative case (see Figure 8 of S15) have been washed away to a certain extent with the inclusion of radiative subgrid physics. Radiative processes dominate the effect that different treatments of artificial viscosity and entropy dissipation have on the entropy profile. That is not to say that codes produce the same profile. Codes with temperature inversions (G3-Magneticum, G2-X) still stand out. However, the key result is that classic SPH codes such as G3-Music and G3-OWLS no longer produce declining entropy profiles with decreasing radius: they now exhibit an almost-flat entropy core. The other classic SPH code, G2-X, still displays a falling entropy inner profile. Subgrid physics can wash away the differences between classic SPH and mesh codes. Interestingly, the modern SPH code G3-PESPH, which produced a falling inner entropy profile more similar to classic SPH in NR simulations than to other modern SPH methods, is now indistinguishable from the Arepo-SH entropy profile. We also note that the introduction of radiative physics in the mesh code AREPO has pushed the entropy profile in the opposite direction. In non-radiative runs, AREPO produces flat entropy cores but it now has a shallow slope in both its subgrid versions. The grid-based code RAMSES shows an almost flat entropy core, although significantly lower than some classic SPH codes such as G3-Music. Another key result is that AGN feedback does not seem to play a dominant role in governing the entropy profile (e.g. the G3-Music and G3-X entropy profiles are similar). In general, codes produce an “almost-flat” central entropy profile, matching the observed overall X-ray profiles (see e.g. Walker et al. 2012). X-ray observations show flat entropy cores for NCC clusters and declining entropy profiles for CC clusters (see also Rasia et al. 2015 for a discussion on the effect of AGN artificial diffusion on the entropy profiles and their relative importance in establishing the cool-coreness of clusters).
A natural follow-up question to ask is whether similar (dis)agreement between codes is seen for the gas fraction (see Figure 11):
Another key result of S15 was that classic SPH codes typically have very baryon rich cores, (R0.1h) 0.4, whereas newer SPH schemes, mesh codes and AREPO produce cores with . Full physics simulations contain little gas in the central regions as a result of star formation, regardless of the code used. Ramses-AGN, Arepo-IL, G2-X and G3-X show a gas fraction that is significantly higher than for the other codes at and, in the case of Ramses-AGN, it exceeds the cosmic ratio outside : as shown in Figure 3, its value at is even higher than in the NR case.
The key systematic difference between codes arises from AGN feedback, which produces in Ramses-AGN, Arepo-IL and G3-X the most evident effect (with the last two showing very similar results). AGN feedback increases gas fractions throughout the cluster with respect to radiative runs with no AGN, especially outside . G3-OWLS is the only code including AGN feedback which has baryon fractions similar to codes with SN feedback only. This difference is in stark contrast to the non-radiative simulations, where with a moderate scatter.
Given the systematic differences presented here, a natural question to ask is which code+subgrid physics is in reasonable agreement with observations of the cluster environment, especially with the aim of using the gas fractions of simulated clusters for cosmological purposes. As pointed out by various studies on the gas fraction of galaxy clusters based on X-ray observations (e.g. LaRoque et al. 2006; Ettori et al. 2010; Zhang et al. 2010; Maughan 2014), gas is expected to account for around 11-12 per cent of the total mass at , which corresponds to approximately 65-70 per cent of the cosmic ratio.
Some AGN codes are in tension with these observations and have 90 per cent (which corresponds to more than 15 per cent of gas with respect to the total mass). All the other codes are largely compatible with these results, and these include include methods both with and without AGN feedback. Moving inward to smaller radii, Zhang et al. (2010) and Vikhlinin et al. (2009) find lower values (around 9-10 per cent) at , in keeping with the general trend of falling gas fractions seen in all simulations (whether NR or not). These values are achived in our comparison by the same set of codes that were found to be in agreement with observational results at .
Nevertheless, in a recent work Gonzalez et al. (2013) suggested that massive clusters may have a higher gas content than what was reported by most of observational studies, estimating a gas fraction around 14 per cent for 210M: these results would support the high gas fraction obtained by codes including AGN, such as Ramses-AGN, Arepo-IL and G3-X. This is supported also by Pratt et al. (2009), which suggests at the same overdensity 14 per cent for massive clusters, measuring values up to 16 per cent for individual clusters.
Here we do not examine the stellar component in detail, i.e., the properties of the galaxies, but defer such an analysis to a companion paper, (Elahi et al., in prep.). Instead we only focus on the overall stellar profiles presented in Figures 12-13. These figures show that the stellar distribution does not extend as far as the gas or dark matter distributions and that galaxies dominate the baryonic content of the central regions. As before, however, the profiles show major code-to-code scatter and systematic differences, and generally a clear separation between codes including and including AGN feedback, with one notable exception, G3-PESPH. The profiles of the star density are shown in Figure 12. All the codes which only include stellar feedback and not AGN show very concentrated stellar densities, around a factor of 5 larger than those of the codes which do include AGN. G2-MusicPI is the code with the highest stellar density within .
Unlike the gas densities, the disagreement does not vanish at the cluster outskirts: gas density profiles are mainly determined by gravity in the outskirts, while star formation is determined by local cooling/feedback. The residuals are flat and non-zero out to well past , as shown in the upper panel of Figure 12; at there is still an order of magnitude difference between the code with the highest stellar density (Arepo-SH) and that with the lowest (Ramses-AGN).
Similarly to the case of the gas component, we define the star fraction as:
and we show the profiles in Figure 13.
Most codes (actually all but Arepo-IL, and to a lesser extent Ramses-AGN and G3-X) have stellar dominated central regions ( 1). The importance of AGN feedback in preventing overcooling is indicated by the fact that only codes without AGN feedback typically have a factor of 2-3 larger than the rest, not only in the cluster core but also in the outskirts. In fact, at the codes including AGN feedback already have 20 per cent, while for the others the star component still accounts for around 40 per cent of the cosmic ratio. At all the codes with AGN feedback show values of below 10 per cent, while codes with only stellar feedback have a mean value of 30 per cent. Ramses-AGN, Arepo-IL and G3-X are the codes which most efficiently reduce star formation, showing 10 per cent already at . Interestingly, G3-PESPH is again an outlier in the codes that do not include AGN feedback, having similar profiles to G3-X and G2-X. Amongst all the codes, Arepo-IL is the only one which does not show a monotonically-falling star fraction, exhibiting a small inversion in the cluster core (this may be due to an offset between the BCG and the cluster center). In general, AGN feedback decreases the stellar fraction by a factor of 80-100 per cent.
Measurements of the stellar mass of galaxy clusters from observations still do not agree: for massive clusters, Giodini et al. (2009) and more recently Gonzalez et al. (2013) reported a star fraction between 1 and 2 per cent (corresponding to about 5-10 per cent of the cosmic ratio for WMAP5-WMAP7 cosmologies), while Sanderson et al. 2013 give a value closer to 3 per cent (15-20 per cent of the cosmic baryon fraction). In spite of these discrepancies, observations seem to agree that in massive clusters of galaxies the star mass does not exceed 3 per cent of the total cluster mass.
Previous works on the baryon contents of hydrodynamical simulations of galaxy clusters have pointed out that methods only including stellar feedback produce an excess of stars (see for instance Sembolini et al. 2013) and star fractions not lower than 5 per cent, while codes which take advantage of AGN feedback are able to reproduce stellar masses compatible with observations (e.g. Planelles et al. 2013). Other detailed comparisons between hydrodynamic simulations and observations can be found in McCarthy et al. (2010) and Le Brun et al. (2014). Our results confirm this trend, as all the codes with AGN feedback succeed in recovering stellar fractions below 3 per cent at . Interestingly, G3-PESPH is the only code with only stellar feedback which reports a stellar fraction more similar to that of the methods including AGN feedback. This can be explained considering the wind model adopted by G3-PESPH, which strongly suppresses low-mass galaxies using high mass loading. This results in a slower buildup of massive galaxy progenitors at early epochs and less dry merger growth within cluster environments at later epochs (see e.g Oppenheimer et al. 2010; Davé, Oppenheimer & Finlator 2011; Davé, Finlator & Oppenheimer 2011).
Nevertheless, as discussed in the previous section, some codes with AGN feedback show an excess of gas mass in their outskirts which makes their values for the gas fraction incompatible with observational results. Surprisingly, codes that include cooling, star formation, and SN feedback are in better agreement with these observations than some of those that also include AGNs, which can give baryon fractions that are too high.
It is also interesting that the codes that recover the most realistic results of the gas and stellar fractions are in general those which have been previously calibrated with observations: for instance, in G3-OWLS the AGN heating temperature has been tuned in order to synthetically reproduce X-ray, SZ and optical cluster properties matching with observations (see Le Brun et al. 2014); G2-X calibrated its AGN model to match the pressure profiles measured by Planck (Planck Collaboration et al., 2013), and tuned the SN feedback parameters to get reasonable agreement with the gas and star fractions (see Pike et al. 2014).
6 Summary & conclusions
This work is the second paper of the nIFTy cluster comparison project series. In the first nIFTy paper, 13 different codes have been used to simulate the same massive cluster, describing the baryon component only by means of non-radiative hydrodynamics: we showed that modern SPH codes are able to reproduce the same results as grid-based codes - same gas density and temperature profiles and a large constant entropy core (S15).
Here, we have studied how cluster properties and code-to-code discrepancies change when the the realism of the description of the baryon component is improved by adding radiative mechanisms - such as cooling, star formation, supernovae feedback, black hole accretion and AGN feedback. We have investigated the performance of 10 modern astrophysical simulation codes - RAMSES, 2 versions of AREPO and 7 versions of Gadget with different SPH implementations. All the simulations have been run using a common set of parameters (e.g. time step accuracy, gravitational softening, dimension of the particle mesh) adopted for S15, but allowing each method to choose the radiative processes modeled by subgrid prescriptions.
We find that - in contrast with what we reported for non-radiative comparison- the differences between classic SPH, modern SPH and grid-based codes are now washed away by the differences in the subgrid physics. The main discrepancies are between codes which include AGN feedback and those which only consider stellar feedback. For instance, the two versions of AREPO show significantly different results in the gas and star fraction. Nevertheless, AGN feedback does not always play a dominant role: in particular, entropy profiles do not seem to be sensitive to the inclusion of AGN feedback. Nevertheless, the addition of radiative models seems to drastically change the entropy cores produced by different codes.
Our main results can be summarized as follows:
Global properties of the cluster -such as the total mass or shape- as calculated at are recovered by all the codes with little scatter (less than 2 per cent) and values extremely close to the non-radiative case. The discrepancies are more evident when we consider baryon properties: temperatures are on average 20 per cent higher than the NR runs with a scatter of 5 per cent. The star fraction - which is the global property most strongly dependent on the chosen subgrid physics- has a scatter larger than 60 per cent.
Although all the codes (except for Ramses-AGN in some cases) agree well on main the features of the cluster, the dark matter distribution appears to have larger scatter amongst the codes than in the of case of the dark matter-only and non- radiative models. This happens in the central regions, where the gravitational field of the cluster is dominated by stars. The dark matter is pulled in by cooling baryons and stars. The codes which do not include AGN feedback, which are those with the highest cooling rates and star fraction in the center, are therefore the ones with the highest dark matter concentrations in the core.
The gas density profiles show a larger scatter than in the non-radiative case. In this case, we do not observe any difference between grid-based, modern and classic SPH codes; similarly, AGN feedback does not seem to play a dominant role.
Temperatures are higher than in the non-radiative case and have a large scatter. Some of the codes show a central temperature inversion, similar to that observed in non-radiative, classic SPH simulations; other codes have a temperature which behaves monotonically, as for modern SPH codes in the non-radiative case. Interestingly, in the full physics case some codes which in their non-radiative version were exhibiting a central inversion of the temperature now show a monotonically decreasing profile (e.g. G3-Music).
Entropy profiles are strongly affected by radiative processes and they present a completely different scenario than in the non-radiative case. The differences between classic SPH codes , which showed an entropy profile falling towards the cluster center, and grid-based and modern SPH codes, which showed a flat entropy core, have now disappeared. Most of codes produce an Õalmost-flatÕ central entropy profile, matching the overall X-ray profiles observed. AGN feedback is not necessary to flatten entropy profiles.
As expected, codes including AGN are able to limit the problem of overcooling and produce a star content compatible with observations. Codes with only stellar feedback show extremely star-dense cores and an excess of stars also in the outskirts.
AGN feedback seems also to increase the fraction of gas at all the cluster radii, especially in the outskirts. Comparison with observational studies of the gas fraction at and show that some of the codes including AGN produce an excess of gas of around 15-20 per cent. Codes which have been previously calibrated with observations generally get more realistic results for the gas and star fractions.
The next papers of the nIFTy cluster comparison series will investigate in more detail how the different codes and physical mechanisms adopted describe a wide range of cluster properties: an upcoming work (Elahi et al., in preparation) will study more deeply the properties of the star component and of substructures; another work (Cui et al., in preparation) will focus on the differences between dark-matter only, non radiative and full physics runs. Subsequent papers will look at the recovery of cluster properties such as X-ray temperature and Sunyaev-ZelÕdovich profiles, gravitational lensing, the cluster outskirts and hydrostatic-mass bias, all of which will add to our understanding of how consistently the results of different codes can inform our understanding of galaxy cluster properties.
The authors would like to express special thanks to the Instituto de Fisica Teorica (IFT-UAM/CSIC in Madrid) for its hospitality and support, via the Centro de Excelencia Severo Ochoa Program under Grant No. SEV-2012-0249, during the three week workshop “nIFTy Cosmology” where this work developed. We further acknowledge the financial support of the University of Western 2014 Australia Research Collaboration Award for “Fast Approximate Synthetic Universes for the SKA”, the ARC Centre of Excellence for All Sky Astrophysics (CAASTRO) grant number CE110001020, and the two ARC Discovery Projects DP130100117 and DP140100198. We also recognize support from the Universidad Autonoma de Madrid (UAM) for the workshop infrastructure.
GY and FS acknowledge support from MINECO (Spain ) through the grant AYA 2012-31101. GY also thanks the Red Espaola de Supercomputacin for granting the computing time in the Marenostrum Supercomputer at BSC. AK is supported by the Ministerio de Economía y Competitividad (MINECO) in Spain through grant AYA2012-31101 as well as the Consolider-Ingenio 2010 Programme of the Spanish Ministerio de Ciencia e Innovación (MICINN) under grant MultiDark CSD2009-00064. He also acknowledges support from the Australian Research Council (ARC) grants DP130100117 and DP140100198. He further thanks Radio Dept for lesser matters.. CP acknowledges support of the Australian Research Council (ARC) through Future Fellowship FT130100041 and Discovery Project DP140100198. WC and CP acknowledge support of ARC DP130100117. EP acknowledges support by the Kavli Foundation and the ERC grant ÒThe Emergence of Structure during the epoch of ReionizationÓ. STK acknowledges support from STFC through grant ST/L000768/1. RJT acknowledges support via a Discovery Grant from NSERC and the Canada Research Chairs program. Simulations were run on the CFI-NSRIT funded Saint Mary’s Computational Astrophysics Laboratory. SB & GM acknowledge support from the PRIN-MIUR 2012 Grant “The Evolution of Cosmic Baryons” funded by the Italian Minister of University and Research, by the PRIN-INAF 2012 Grant “Multi-scale Simulations of Cosmic Structures”, by the INFN INDARK Grant and by the “Consorzio per la Fisica di Trieste”. IGM acknowledges support from a STFC Advanced Fellowship. DN, KN, and EL are supported in part by NSF AST-1009811, NASA ATP NNX11AE07G, NASA Chandra grants GO213004B and TM4-15007X, the Research Corporation, and by the facilities and staff of the Yale University Faculty of Arts and Sciences High Performance Computing Center. PJE is supported by the SSimPL program and the Sydney Institute for Astronomy (SIfA), DP130100117. JIR acknowledges support from SNF grant PP00P2 128540/1. CDV acknowledges financial support from the Spanish Ministry of Economy and Competitiveness (MINECO) through the 2011 Severo Ochoa Program MINECO SEV-2011-0187 and the AYA2013-46886-P grant. AMB is supported by the DFG Research Unit 1254 ’Magnetisation of Interstellar and Intergalactic Media’ and by the DFG Cluster of Excellence ’Origin and Structure of the Universe’. RDAN acknowledges the support received from the Jim Buckee Fellowship. The AREPO simulations were performed with resources awarded through STFCs DiRAC initiative. The authors thank Volger Springel for helpful discussions and for making AREPO and the original GADGET version available for this project. JS acknowledges support from the European Research Council under the European UnionÕs Seventh Framework Programme (FP7/2007-2013)/ERC grant agreement 278594-GasAroundGalaxies. G3-PESPH Simulations were carried out using resources at the Center for High Performance Computing in Cape Town, South Africa
The authors contributed to this paper in the following ways: FS, GY, FRP, AK, CP, STK & WC formed the core team that organized and analyzed the data, made the plots and wrote the paper. AK, GY & FRP organized the nIFTy workshop at which this program was initiated. GY supplied the initial conditions. PJE assisted with the analysis. All the other authors, as listed in Section 2 performed the simulations using their codes. All authors had the opportunity to proof read and comment on the paper.
The simulation used for this paper has been run on Marenostrum supercomputer and is publicly available at the MUSIC website
This research has made use of NASA’s Astrophysics Data System (ADS) and the arXiv preprint server.
- pagerange: nIFTy galaxy cluster simulations II: radiative models–References
- pubyear: 2015
- Specifically, it is cluster 19 of MUSIC-2 dataset; all the initial conditions of MUSIC clusters are available at http://music.ft.uam.es
- A dark-matter only simulation containing 2048 particles in a (1Gpc) cube performed using ART (Kravtsov, Klypin & Khokhlov, 1997) at the NASA Ames Research centre Prada et al. (2012)
- Aarseth S. J., Turner E. L., Gott, III J. R., 1979, ApJ, 228, 664
- Aubert D., Teyssier R., 2010, ApJ, 724, 244
- Battaglia N., Bond J. R., Pfrommer C., Sievers J. L., 2012, ApJ, 758, 74
- Beck A. M. et al., 2015, ArXiv e-prints
- Biffi V., Sembolini F., De Petris M., Valdarnini R., Yepes G., Gottlöber S., 2014, MNRAS, 439, 588
- Bondi H., 1952, MNRAS, 112, 195
- Bondi H., Hoyle F., 1944, MNRAS, 104, 273
- Booth C. M., Schaye J., 2009, MNRAS, 398, 53
- Borgani S., Kravtsov A., 2011, Advanced Science Letters, 4, 204
- Borgani S. et al., 2004, MNRAS, 348, 1078
- Capela F., Pshirkov M., Tinyakov P., 2014, Phys. Rev. D, 90, 083507
- Chabrier G., 2003, PASP, 115, 763
- Churazov E., Sazonov S., Sunyaev R., Forman W., Jones C., Böhringer H., 2005, MNRAS, 363, L91
- Dalla Vecchia C., Schaye J., 2008, MNRAS, 387, 1431
- Dalla Vecchia C., Schaye J., 2012, MNRAS, 426, 140
- Davé R., Finlator K., Oppenheimer B. D., 2011, MNRAS, 416, 1354
- Davé R., Katz N., Oppenheimer B. D., Kollmeier J. A., Weinberg D. H., 2013, MNRAS, 434, 2645
- Davé R., Katz N., Weinberg D. H., 2002, ApJ, 579, 23
- Davé R., Oppenheimer B. D., Finlator K., 2011, MNRAS, 415, 11
- Davé R., Oppenheimer B. D., Sivanandam S., 2008, MNRAS, 391, 110
- Dehnen W., Aly H., 2012, MNRAS, 425, 1068
- Di Matteo T., Springel V., Hernquist L., 2005, Nature, 433, 604
- Dolag K., Gaensler B. M., Beck A. M., Beck M. C., 2015, MNRAS, 451, 4277
- Dolag K., Jubelgas M., Springel V., Borgani S., Rasia E., 2004, ApJ, 606, L97
- Dolag K., Vazza F., Brunetti G., Tormen G., 2005, MNRAS, 364, 753
- Duffy A. R., Schaye J., Kay S. T., Dalla Vecchia C., Battye R. A., Booth C. M., 2010, MNRAS, 405, 2161
- Eggen O. J., Lynden-Bell D., Sandage A. R., 1962, ApJ, 136, 748
- Ettori S., Gastaldello F., Leccardi A., Molendi S., Rossetti M., Buote D., Meneghetti M., 2010, A&A, 524, A68
- Fabjan D., Borgani S., Tornatore L., Saro A., Murante G., Dolag K., 2010, MNRAS, 401, 1670
- Fall S. M., 1978, MNRAS, 185, 165
- Faucher-Giguère C.-A., Lidz A., Zaldarriaga M., Hernquist L., 2009, ApJ, 703, 1416
- Ferland G. J., Korista K. T., Verner D. A., Ferguson J. W., Kingdon J. B., Verner E. M., 1998, PASP, 110, 761
- Ferland G. J. et al., 2013, Rev. Mexicana Astron. Astrofis., 49, 137
- Frenk et al. C. S., 1999, ApJ, 525, 554
- Gill S. P. D., Knebe A., Gibson B. K., 2004, MNRAS, 351, 399
- Giodini S. et al., 2009, ApJ, 703, 982
- Gnedin O. Y., Kravtsov A. V., Klypin A. A., Nagai D., 2004, ApJ, 616, 16
- Gonzalez A. H., Sivanandam S., Zabludoff A. I., Zaritsky D., 2013, ApJ, 778, 14
- Gonzalez A. H., Zaritsky D., Zabludoff A. I., 2007, ApJ, 666, 147
- Haardt F., Madau P., 1996, ApJ, 461, 20
- Haardt F., Madau P., 2001, in Clusters of Galaxies and the High Redshift Universe Observed in X-rays, Neumann D. M., Tran J. T. V., eds., p. 64
- Hirschmann M., Dolag K., Saro A., Bachmann L., Borgani S., Burkert A., 2014, MNRAS, 442, 2304
- Hopkins P. F., 2013, MNRAS, 428, 2840
- Hoyle F., Lyttleton R. A., 1939, Proceedings of the Cambridge Philosophical Society, 35, 592
- Katz N., Weinberg D. H., Hernquist L., 1996, ApJS, 105, 19
- Kay S. T., Thomas P. A., Jenkins A., Pearce F. R., 2004, MNRAS, 355, 1091
- Kennicutt, Jr. R. C., 1998, ApJ, 498, 541
- Kim J.-h. et al., 2014, ApJS, 210, 14
- Knollmann S. R., Knebe A., 2009, ApJS, 182, 608
- Komatsu E. et al., 2011, ApJS, 192, 18
- Kravtsov A. V., Klypin A. A., Khokhlov A. M., 1997, ApJS, 111, 73
- LaRoque S. J., Bonamente M., Carlstrom J. E., Joy M. K., Nagai D., Reese E. D., Dawson K. S., 2006, ApJ, 652, 917
- Le Brun A. M. C., McCarthy I. G., Schaye J., Ponman T. J., 2014, MNRAS, 441, 1270
- Lewis G. F., Babul A., Katz N., Quinn T., Hernquist L., Weinberg D. H., 2000, ApJ, 536, 623
- Martizzi D., Teyssier R., Moore B., Wentz T., 2012, MNRAS, 422, 3081
- Maughan B. J., 2014, MNRAS, 437, 1171
- McCarthy I. G., Schaye J., Bower R. G., Ponman T. J., Booth C. M., Dalla Vecchia C., Springel V., 2011, MNRAS, 412, 1965
- McCarthy I. G. et al., 2010, MNRAS, 406, 822
- McKee C. F., Ostriker J. P., 1977, ApJ, 218, 148
- Metzler C. A., Evrard A. E., 1994, ApJ, 437, 564
- Mitchell N. L., McCarthy I. G., Bower R. G., Theuns T., Crain R. A., 2009, MNRAS, 395, 180
- Morris J. P., Monaghan J. J., 1997, Journal of Computational Physics, 136, 41
- Muanwong O., Thomas P. A., Kay S. T., Pearce F. R., Couchman H. M. P., 2001, ApJ, 552, L27
- Nagai D., Kravtsov A. V., 2004, in IAU Colloq. 195: Outskirts of Galaxy Clusters: Intense Life in the Suburbs, Diaferio A., ed., pp. 296–298
- Nagai D., Kravtsov A. V., Vikhlinin A., 2007, ApJ, 668, 1
- Oppenheimer B. D., Davé R., 2008, MNRAS, 387, 577
- Oppenheimer B. D., Davé R., Kereš D., Fardal M., Katz N., Kollmeier J. A., Weinberg D. H., 2010, MNRAS, 406, 2325
- Padovani P., Matteucci F., 1993, ApJ, 416, 26
- Pakmor R., Edelmann P., Röpke F. K., Hillebrandt W., 2012, MNRAS, 424, 2222
- Pearce F. R., Thomas P. A., Couchman H. M. P., Edge A. C., 2000, MNRAS, 317, 1029
- Pike S. R., Kay S. T., Newton R. D. A., Thomas P. A., Jenkins A., 2014, MNRAS, 445, 1774
- Piontek F., Steinmetz M., 2011, MNRAS, 410, 2625
- Planck Collaboration et al., 2013, A&A, 550, A131
- Planelles S., Borgani S., Dolag K., Ettori S., Fabjan D., Murante G., Tornatore L., 2013, MNRAS, 431, 1487
- Prada F., Klypin A. A., Cuesta A. J., Betancort-Rijo J. E., Primack J., 2012, MNRAS, 423, 3018
- Pratt G. W., Croston J. H., Arnaud M., Böhringer H., 2009, A&A, 498, 361
- Price D. J., 2008, J. Comp. Phys., 227, 10040
- Price D. J., 2012, Journal of Computational Physics, 231, 759
- Puchwein E., Sijacki D., Springel V., 2008, ApJ, 687, L53
- Puchwein E., Springel V., Sijacki D., Dolag K., 2010, MNRAS, 406, 936
- Rafieferantsoa M., Davé R., Anglés-Alcazar D., Katz N., Kollmeier J. A., Oppenheimer B. D., 2014, ArXiv e-prints
- Rahmati A., Pawlik A. H., Raievic M., Schaye J., 2013, MNRAS, 430, 2427
- Rasera Y., Teyssier R., 2006, A&A, 445, 1
- Rasia E. et al., 2015, ArXiv e-prints
- Rasia E. et al., 2014, ApJ, 791, 96
- Ritchie B. W., Thomas P. A., 2001, MNRAS, 323, 743
- Saitoh T. R., Makino J., 2013, ApJ, 768, 44
- Salpeter E. E., 1955, ApJ, 121, 161
- Sanderson A. J. R., O’Sullivan E., Ponman T. J., Gonzalez A. H., Sivanandam S., Zabludoff A. I., Zaritsky D., 2013, MNRAS, 429, 3288
- Scannapieco C., Wadepuhl M., Parry O. H., Navarro J. F., Jenkins A., Springel V., Teyssier R., Carlson E. e. a., 2012, MNRAS, 423, 1726
- Schaye J., 2004, ApJ, 609, 667
- Schaye J. et al., 2015, MNRAS, 446, 521
- Schaye J., Dalla Vecchia C., 2008, MNRAS, 383, 1210
- Schaye J. et al., 2010, MNRAS, 402, 1536
- Schmidt M., 1959, ApJ, 129, 243
- Sembolini F., De Petris M., Yepes G., Foschi E., Lamagna L., Gottlöber S., 2014, MNRAS, 440, 3520
- Sembolini F., Yepes G., De Petris M., Gottlöber S., Lamagna L., Comis B., 2013, MNRAS, 429, 323
- Sembolini F. et al., 2015, ArXiv e-prints
- Sijacki D., Pfrommer C., Springel V., Enßlin T. A., 2008, MNRAS, 387, 1403
- Sijacki D., Springel V., Di Matteo T., Hernquist L., 2007, MNRAS, 380, 877
- Sijacki D., Vogelsberger M., Genel S., Springel V., Torrey P., Snyder G., Nelson D., Hernquist L., 2014, ArXiv e-prints
- Springel V., 2010, MNRAS, 401, 791
- Springel V., Hernquist L., 2003, MNRAS, 339, 289
- Springel V. et al., 2005a, Nature, 435, 629
- Springel V. et al., 2005b, Nature, 435, 629
- Steinborn L. K., Dolag K., Hirschmann M., Prieto M. A., Remus R.-S., 2015, MNRAS, 448, 1504
- Suginohara T., Ostriker J. P., 1998, ApJ, 507, 16
- Sutherland R. S., Dopita M. A., 1993, ApJS, 88, 253
- Teyssier R., Moore B., Martizzi D., Dubois Y., Mayer L., 2011, MNRAS, 414, 195
- Thacker R. J., Scannapieco E., Couchman H. M. P., 2006, ApJ, 653, 86
- Thielemann F.-K. et al., 2003, Nuclear Physics A, 718, 139
- Thomas P. A., Couchman H. M. P., 1992, MNRAS, 257, 11
- Tornatore L., Borgani S., Dolag K., Matteucci F., 2007, MNRAS, 382, 1050
- Tornatore L., Borgani S., Springel V., Matteucci F., Menci N., Murante G., 2003, MNRAS, 342, 1025
- Valdarnini R., 2012, A&A, 546, A45
- van den Hoek L. B., Groenewegen M. A. T., 1997, A&AS, 123, 305
- Vikhlinin A. et al., 2009, ApJ, 692, 1033
- Vogelsberger M., Genel S., Sijacki D., Torrey P., Springel V., Hernquist L., 2013, MNRAS, 436, 3031
- Vogelsberger M. et al., 2014, Nature, 509, 177
- Voit G. M., 2005, Reviews of Modern Physics, 77, 207
- Wadsley J. W., Veeravalli G., Couchman H. M. P., 2008, MNRAS, 387, 427
- Walker S. A., Fabian A. C., Sanders J. S., George M. R., 2012, MNRAS, 427, L45
- White S. D. M., 1976, MNRAS, 177, 717
- Wiersma R. P. C., Schaye J., Smith B. D., 2009, MNRAS, 393, 99
- Wiersma R. P. C., Schaye J., Theuns T., Dalla Vecchia C., Tornatore L., 2009, MNRAS, 399, 574
- Woosley S. E., Weaver T. A., 1995, ApJS, 101, 181
- Zel’dovich Y. B., Klypin A. A., Khlopov M. Y., Checketkin, V M., 1980, Sov. J. nucl. Phys., 31, 664
- Zhang Y.-Y. et al., 2010, ApJ, 711, 1033