Comparing Simulations of AGN Feedback
Mark L. A. Richardson, Evan Scannapieco, Julien Devriendt, Adrianne Slyz, Robert J. Thacker, Yohan Dubois, James Wurster, Joseph Silk
AAccepted to ApJ
Preprint typeset using L A TEX style emulateapj v. 5/2/11
COMPARING SIMULATIONS OF AGN FEEDBACK
Mark L. A. Richardson , Evan Scannapieco , Julien Devriendt , Adrianne Slyz , Robert J. Thacker , YohanDubois , James Wurster , Joseph Silk Accepted to ApJ
ABSTRACTWe perform adaptive mesh refinement (AMR) and smoothed particle hydrodynamics (SPH) cos-mological zoom simulations of a region around a forming galaxy cluster, comparing the ability of themethods to handle successively more complex baryonic physics. In the simplest, non-radiative case,the two methods are in good agreement with each other, but the SPH simulations generate centralcores with slightly lower entropies and virial shocks at slightly larger radii, consistent with what hasbeen seen in previous studies. The inclusion of radiative cooling, star formation, and stellar feedbackleads to much larger differences between the two methods. Most dramatically, at z = 5 , rapid coolingin the AMR case moves the accretion shock well within the virial radius, while this shock remainsnear the virial radius in the SPH case, due to excess heating, coupled with poorer capturing of theshock width. On the other hand, the addition of feedback from active galactic nuclei (AGN) to thesimulations results in much better agreement between the methods. In this case both simulationsdisplay halo gas entropies of 100 keV cm , similar decrements in the star-formation rate, and a dropin the halo baryon content of roughly 30%. This is consistent with AGN growth being self-regulated,regardless of the numerical method. However, the simulations with AGN feedback continue to differ inaspects that are not self-regulated, such that in SPH a larger volume of gas is impacted by feedback,and the cluster still has a lower entropy central core. INTRODUCTION
In the hierarchical cold dark matter and dark energy(ΛCDM) model for galaxy formation, matter condensesinto small clumps that then merge to create increasinglymassive objects over time. This model has provided sev-eral predictions that are in excellent agreement with ob-servations (e.g., Spergel et al. 2007; Larson et al. 2011).For star formation to begin in ΛCDM, the temperatureof gas contained within condensed dark matter ‘halos’must cool sufficiently to allow the formation of galaxies.Because larger galaxies have more gravitational compres-sion, and hence a higher temperature, they might be ex-pected to take longer to cool and form stars, with thelargest galaxies only now reaching significant star forma-tion rates (SFRs).Observational studies of star formation rates, on theother hand, have found several surprising trends: thecosmic star formation rate density reaches a peak at z (cid:39)
2, 5-3 billion years after the Big Bang and thendecreases until today (Hopkins & Beacom, 2006; Karimet al. 2011), SFRs peaked earlier in more massive galax-ies and more recently in smaller galaxies (e.g. Guzman Sub-department of Astrophysics, University of Oxford, Ke-ble Road, Oxford, OX1 3RH School of Earth and Space Exploration, Arizona State Uni-versity, Tempe, AZ 85287 Observatoire de Lyon, UMR 5574, 9 avenue Charles Andr´e,F-69561 Saint Genis Laval, France Department of Astronomy & Physics, Saint Mary’s Univer-sity, Halifax, NS, B3L 3C3, Canada Sorbonne Universit´es, UPMC Univ Paris 06, UMR 7095, In-stitut d’Astrophysique de Paris, F-75014, Paris, France CNRS, UMR 7095, Institut d’Astrophysique de Paris, F-75014, Paris, France Monash Centre for Astrophysics and School of Physics andAstronomy, Monash University, Victoria, 3800, Australia Department of Physics and Astronomy, The Johns HopkinsUniversity Homewood Campus, Baltimore, MD 21218, USA et al. 1997; Brinchmann & Ellis 2000), and today SFRsare lower in more massive galaxies (Heavens et al. 2004;Panter et al. 2007). These trends, known collectively as‘downsizing’ (e.g. Cowie et al. 1996; Bauer et al. 2005;Panter et al. 2007; Karim 2011), are clearly at odds withthe naive predictions of the ΛCDM model.A similar such discrepancy involves the properties ofgalaxy clusters, as constrained by high-resolution X-rayand radio observations such as those from the ChandraObservatory and the Very Large Array. While many clus-ters appear to be quiescent, about a third show strongpeaks in their central X-ray surface brightness distribu-tions, indicating that their gas is cooling rapidly (e.g.Fabian & Nulsen 1977; Nulsen et al. 1982; Stewart etal. 1984; Fabian 1994; Tamura et al. 2001; Cavagnolo etal. 2009). However, this cooling is neither accompaniedby strong star formation nor a significant fraction of gascolder than 1 keV (e.g. Peterson et al. 2001; Raffertyet al. 2006; McNamara & Nulsen 2007). Instead galaxyformation is halted by an unknown energy source (e.g.Croton et al. 2006).A prominent theory is that energetic feedback fromactive galactic nuclei (AGN) is required to explain thesetwo discrepancies (e.g., Scannapieco & Oh 2004; Springelet al. 2005; Thacker et al. 2006; Dunn & Fabian 2006;Sijacki et al. 2007; Booth & Schaye 2009; Dubois etal. 2013; Martizzi et al. 2013). AGN are among themost energetic objects in the Universe, characterizedby their extremely luminous cores powered by the in-fall of gas from a relativistic accretion disk (e.g., Rees1984) onto a supermassive black hole (SMBH) with mass M BH > M (cid:12) . AGN are associated with two modes offeedback into their environments. The kinematic or radiomode is associated with collimated relativistic jets, andlow and radiatively inefficient accretion rates (e.g., Fal-cke & Biermann, 1999; Sambruna et al. 2000; Merloni & a r X i v : . [ a s t r o - ph . GA ] M a y Heinz 2007), while the quasar or wind mode is associatedwith isotropic energy deposit, and high and radiativelyefficient accretion (e.g., Silk & Rees 1998). Dependingon the feedback mode, it has been shown that AGN canprovide the energy needed to maintain the hot ICM (e.g.,Dunn & Fabian 2006), with kinematic feedback creatingthe large buoyant bubbles (e.g., Dunn et al. 2006). Itis also clear that this feedback can hamper cooling ofgalactic halo gas, preferentially reducing the SFR firstin large halos at early times, and then smaller halos atlate times (e.g., Scannapieco & Oh 2004; Scannapieco etal. 2005). AGN also act to remove, via accretion andheating, the low specific angular momentum (sAM) gasfrom the central region of its halo, gas that would oth-erwise form stars. This leads to a net increase in thegalactic sAM. However, AGN can also increase galacticsAM where feedback has heated gas and reduced the ofpossibly high angular momentum gas into the galacticdisk (e.g., Dubois et al. 2013; Genel et al. 2015; Nelsonet al. 2015), .Unfortunately, AGN feedback is extremely difficult tosimulate as its effects span several orders of magnitude,originating on sub-parsec scales, and impacting kilopar-sec and even megaparsec scales. Thus numerical methodsmust implement a subgrid prescriptions for injecting thefeedback model if they are resolving cluster or cosmo-logical scales. The nature of this feedback is also highlydebated, with different studies focusing on different in-put mechanisms for the feedback energy, and differentenvironments in which to study its effects.The first models of AGN feedback were limited tosmoothed particle hydrodynamics (SPH) simulations.Springel et al. (2005) used the SPH code
GADGET-2 (Springel 2005) to study individual and merging galax-ies with black holes (BHs) and feedback, injected intothe simulation thermally and isotropically. They foundthat the feedback energy was sufficient to regulate theBH growth. Thacker et al. (2006b; 2009) carried outSPH simulations using the
HYDRA (Couchman et al. 1995;Thacker & Couchman 2006) code with an isotropic kine-matic outflow model for AGN feedback. Their resultsreproduced the antihierarchical turnoff in the quasar lu-minosity function, as well as the spatial distribution ofquasars on both small and large scales. However, the im-pact of feedback was significantly less than predicted byanalogous semianalytic models, a difference that couldbe traced to in-shock cooling as occurred in they SPHsimulation. Sijacki et al. (2007) introduced a dual modeprescription for their AGN simulations in
GADGET-2 , re-sulting in better agreement between the simulated galaxystellar mass density and observations than previous mod-els. Booth & Schaye (2009) performed SPH simulationsusing the
GADGET III code with a modified version of theAGN prescription of Springel et al. (2005), and foundthat the BHs greatly suppressed the star formation inhigh-mass galaxies, self-regulating their feedback suchthat they agreed with the M BH - σ v relation, a tight cor-relation between BH mass, M BH and the velocity disper-sion of the host galaxy’s bulge, σ v (e.g., Ferrarese & Mer-ritt 2000; Gebhardt et al. 2000; Tremaine et al. 2002).Recently, Planelles et al. (2014), Le Brun et al. (2014),and Pike et al. (2014) used variants of GADGET III and
GADGET-2 to study the impact of AGN on a range ofgalaxy masses and clusters, and how AGN models can be tuned to produce better agreement between simulatedand observed galaxy groups and clusters. Finally, Baraiet al. (2013) and Wurster & Thacker (2013) each usedSPH simulations (
GADGET III and
HYDRA , respectively)to consider the impact of various AGN feedback modelson isolated galaxies and mergers, showing some successin replicating the M BH - σ v relation.All such SPH simulations, which use Lagrangianschemes to solve the equations of fluid dynamics, are veryefficient at resolving dense structures. However, SPH isnot without its shortcomings. As particles moves alongwith the mass, SPH is ill-equipped to resolve the low-density environment surrounding galaxies and clusters,although it is through this medium that the feedbackinteracts with the surrounding structure. Furthermore,traditional, also called standard, SPH has difficulties ac-curately modeling shocks and mixing (e.g., Morris 1996;Marri & White 2003; Agertz et al. 2007; Hopkins 2013),which are essential when studying the impact of feedbackon surrounding structure. Fortunately, there is on-goingeffort and success in reformulating the SPH method toovercome the mixing and shock issues discussed above(e.g., Price 2012; Hopkins 2013), although they have notyet become the standard practice.Thus, recently adaptive mesh refinement (AMR) meth-ods are gaining interest, as AMR is a shock-capturingmethod, it can set high spatial resolution in any regionof the simulation volume and resolves shocks with onlya few cells. Dubois et al. (2010) introduced kinematicradio-mode feedback in RAMSES simulations, and then inDubois et al. (2013) used both radio and quasar mod-els of AGN feedback in the AMR code
RAMSES (Teyssieret al. 2002) to study the growth of a galaxy cluster ina full cosmological simulation. Their study focused onthe accretion history of the cluster SMBH and the effectof feedback on the gas content and temperature. Theyfound that only with AGN feedback were they able togreatly heat the gas and affect its ability to accrete ontothe central galaxy, thus limiting the overall SFR andaccretion on to the SMBH, while being in good agree-ment with the empirical M BH - σ v relation. Martizzi et al.(2013) also used RAMSES to simulate an isolated clusterhalo, studying how quickly cluster gas heated by AGNfeedback cools back into the central region.Aside from using different hydrodynamical methods,these collection of studies are quite diverse in how theymodel the deposited feedback energy and its impact onthe environment. Thus the role of the simulation methodin determining the conclusions of these studies is diffi-cult to disentangle in the absence of a comparison thatemploys the same physical model across different simu-lation codes. The importance of comparisons betweensimulation techniques has also been highlighted by theusefulness of such studies in models of structure forma-tion without AGN feedback. For example, the SantaBarbara Cluster Comparison Project (Frenk et al. 1999)compared the results of twelve numerical codes using thesame initial conditions to study the virialization of a mas-sive galaxy cluster, without including feedback or radia-tive cooling, finding that agreement was best for prop-erties of the dark matter and worst for the total X-rayluminosity. Similar comparisons with higher resolutionwere done in Voit et al. (2005) and then Mitchell et al.(2009), focusing on the central entropy profiles of clus-ters, and they showed that the core entropy was lowerin SPH than in AMR due to under-mixing in SPH, andto a lesser degree, over-mixing in AMR. A detailed com-putational comparison between AMR and SPH methodsfor standard hydrodynamic turbulence was performed byAgertz et al. (2007), which demonstrated the inherentdifficulties of modeling shear layers with standard SPH.This work also highlighted the effect of steep density gra-dients in standard SPH simulations, where an effectivesurface tension between the two phases would lead to lim-ited mixing, overcooling, and angular momentum trans-port (Kaufmann et al. 2007).Code comparisons continue to be important. For theAquila project, Scannapieco et al. (2012) studied theformation of a Milky Way-size galaxy, wherein 13 differ-ent numerical methods were used, starting from the sameinitial conditions, with no attempt to use identical sub-grid models (gas cooling and the formation and feedbackof stars and AGN). They found that while the varietyof different gas physics led to a large span in the phys-ical characteristics of the final galaxy, it was inconsis-tent with observations of real galaxies. They also foundthat gas cooled more efficiently in grid codes, leading tohigher star formation rates. In contrast to the Aquilaproject, the now underway AGORA project is combin-ing the works of 95 scientists to use a variety of numeri-cal codes with as identical as possible implementations ofvarious baryonic physics to produce more observationallyconsistent galaxies (Kim et al. 2014). There has been asuite of work comparing galactic and cosmological sim-ulations from the recently introduced moving-mesh code
AREPO , with the SPH code
GADGET III (e.g., Springel2010; Bauer & Springel 2012; Kereˇs et al. 2012; Sijackiet al. 2012; Torrey et al. 2012a; Vogelsberger et al. 2012;Nelson et al. 2013), which only recently included AGNfeedback (Hayward et al. 2014). These have furtherdemonstrated where shortcomings exist in the ability ofSPH codes to mix merging material. In particular, Kereˇset al. (2012) found that with radiative gas cooling
AREPO had either the same or lower central entropy profiles de-pending on the halo mass, which is opposite the resultsof grid codes in non-radiative simulations (e.g., Voit etal. 2005; Mitchell et al. 2009). Finally, the nIFTy Cos-mology workshop has also lead to simulation comparisonpapers that include AMR and SPH codes, and the mov-ing mesh code
AREPO (Sembolini et al. 2015b, 2015a).For non-radiative simulations they reproduced the radialprofile results of Mitchell et al. (2009). By introduc-ing cooling and AGN feedback, at z = 0 they see betterqualitative agreement in the halo profiles, but with largerscatter.In this work, we wish to continue the effort to com-pare the results from different codes as they simulate thecosmologically consistent formation of a cluster environ-ment, including AGN feedback. We perform two simula-tion suites, one with AMR and one with standard SPH,from the same initial conditions, studying the impact ofdifferent subgrid physics, including cooling, star forma-tion, stellar and AGN feedback models. Note that weattempt to implement nearly identical sub-grid baryonicphysics models, using the same parameter values in thesemodels to emphasize the role of the numerical method.We compare the ability of these two numerical methodsto model the evolution of the cluster environment and its response to these subgrid models, including the gas tem-perature, SFRs, and gas content. We stress that theseresults are only applicable to standard SPH implemen-tations, and future work comparing with non-standardSPH implementations is required. In a companion paper(Richardson et al. 2016 in prep) we compare the AMRand SPH impact of AGN feedback on the characteristicsof halos ranging from 10 to 10 . M (cid:12) .The structure of this paper is as follows. In § § § NUMERICAL METHODS
Simulations were conducted with either the AMR code
RAMSES (Teyssier 2002) or the SPH code
HYDRA (Couch-man et al. 1995). The initial conditions were generatedusing the mpgrafic (Prunet et al. 2008) package whichcreates a realization of the density fluctuations on a gridaccording to the desired power spectrum, and uses theZel’dovich approximation to calculate the correspondingparticle velocities. Our initial conditions were generatedat a redshift of z = 43 .
2, centered on a region in which acluster halo with virial mass M vir = 2 × M (cid:12) formsby z = 0. We assumed a ΛCDM cosmology with cosmo-logical parameters (Ω ∆ , Ω M , Ω b , σ , h ) = (0.73, 0.27,0.044, 0.8, 0.7) from the 7-year WMAP (Komatsu etal. 2011). Our simulations were carried out in a 100 h − Mpc comoving box with periodic boundaries andwere run to z = 3. In both particle and grid simulationswe assumed a zoomed-in realization of this box, with aspherical high-resolution region 25 h − Mpc in diameter,an effective dark matter particle number of 1024 , anda mass resolution of 8 . × M (cid:12) . This high resolutionregion was selected to contain all particles found in thehalo by z = 1, and thus constitutes a conservative esti-mate of the necessary high-resolution region for z = 3.Outside of the high resolution region, we uniformly de-creased the dark matter particle resolution until reachingan effective particle number of 64 , for AMR, and 128 ,for SPH, in the outer regions of the box. The dark mat-ter particle initial conditions were identical between thegrid and particle simulations.We carried out the grid simulations with the RAMSES code (Teyssier 2002), which uses an unsplit second-orderGodunov scheme for evolving the Euler equations forthe gas.
RAMSES variables are cell-centered and inter-polated to the cell faces for flux calculations, which arethen used with a Harten-Lax-van Leer-Contact Riemannsolver (van Leer 1979; Einfeldt 1988). For calculatingthe gravitational potential and accelerations, collisionlessstar, black hole, and dark matter particles had their massmapped to the grid with a Cloud-in-Cell (CIC) scheme(Birdsall & Fuss, 1997). The CIC method was also usedfor comparing gas and particle densities for star forma-tion and sink particle generation, discussed further be-low. Gas within the high-resolution region was refinedusing a semi-Lagrangian technique. When more than 8dark matter particles were in a cell, or when the baryondensity in a cell was 8 times more than the cosmic aver-age, the cell was split into 8, doubling the spatial reso-lution. We aimed for a fixed maximum physical resolu-tion of ∆ x min = 545 pc, where the maximum refinementlevel was increased with increasing cosmic scale factor,with increments occurring at z (cid:39)
39, 19, 9, and 4. Thusthe spatial resolution at any one time varied from 435physical pc to 870 physical pc. To avoid over-resolvingthe dark matter in dense gas regions, we set the maxi-mum level to map the dark matter particles into cells at l max , DM = 15. For comparison, at z = 4, the gas wasrefined up to l max , g = 16 levels of resolution. Thus from z = 4 to 3 the dark matter had an effective softeninglength of 2 physical kpc, given by twice the width of gridcell at l = 15.We carried out the particle-based simulations with the HYDRA code (Couchman et al. 1995; Thacker & Couch-man 2006), which uses an adaptive particle-particle,particle-mesh method (Couchman 1991) to calculatethe gravitational forces and an SPH method (Gingold& Marigold 1977; Lucy 1977) to calculate the hydrody-namic forces. This is a standard implementation of SPHwithout a forced conservation of entropy (e.g. Springel2005), and does not include methods for better captur-ing contact discontinuities (e.g. Hopkins et al. 2013).We used the same initial conditions for dark matter asfor the
RAMSES run, and overlaid the gas particles ontothe dark matter particle positions, which then trace thegas density field. Note,
HYDRA employs a modification ofthe kernel gradient to prevent the formation of the pairinstability (Thomas & Couchman 1992), and where pres-sure is important we confirm that gas and DM particlesseparate within a few time steps.
HYDRA uses the S2 grav-itational softening length (Hockney & Eastwood 1981), (cid:15) , for mediating close encounter scattering events, withthe minimum smoothing length given by h min = (cid:15)/ h min = 2∆ x min = 1090 physical pc, giving thesame softening length for dark matter particles in thetwo codes at late redshifts. Although the spatial resolu-tion of a grid or particle method is not exactly equal tothese two quantities, we found that the star formationhistories were sufficiently similar when relating these pa-rameters in this way. We discuss this further in § K, and metal fine-structure cooling following Rosen & Bregman (1995) atcooler temperatures. The gas temperature was not al-lowed to drop below T = 500 K via radiative and metalline cooling, although adiabatic cooling below this limitwas permitted. After z = 8 . , heating from an ultravio-let background was modeled following Haardt & Madau(1996). For simplicity and consistency between codes,the metallicity was set at a constant value of a third solarfor the entirety of the simulation, as typically observed inthe intracluster medium (Loewenstein 2004). We notethat a recently identified error in the cooling prescrip-tion in RAMSES shows that its metal cooling tables also include the contribution from hydrogen and helium, thusin regions where H and He are the main source of cooling(e.g., high temperature synchrotron tail) we were addingan extra third of cooling (since this extra term scales withthe metallicity). However, for consistency, the same cool-ing tables were used in the
HYDRA simulation, thus thiserror does not account for any differences between thetwo methods.
Star Formation & Stellar Feedback
In both types of simulations, star formation followeda Schmidt-Kennicut law (Schmidt 1955; Kennicut 1998),with the star formation rate given as dM ∗ dt = c sf M g t ff , (1)where c sf is the star forming efficiency, M g is the gasmass, and t ff ≡ (32 Gρ/ π ) / is the gravitational free-fall time in the vicinity of the star-forming region, with G the gravitational constant and ρ the gas mass density ina given resolution element. We set c sf = 0 . , consistentwith observations of giant molecular clouds (Krumholz& Tan 2007). Star formation was implemented wherethe gas density was above a hydrogen number densitythreshold, sufficient to overcome the local hydrodynamicpressure. For these resolutions, we set this threshold tobe n ∗ = 0 .
05 hydrogen atoms cm − . We also incorpo-rate a baryon overdensity threshold, comparing the localdensity of gas and stars to dark matter. This is neces-sary to prevent star formation from proceeding in cosmicfilaments at high redshift. For SPH, if the total baryonmass within two smoothing lengths of a gas particle wasmore than 0.25 that of the enclosed dark matter, then theparticle was permitted to proceed with the star forma-tion prescription discussed below. For AMR, we checkedon a cell-by-cell case, and required a stricter thresholdthan SPH to ensure the gas in the smaller volume wasindeed dominating the local mass. We allowed this tovary with time since the dark matter softening transi-tioned to 1 kpc, and ensured agreement with the globalstar formation rate in the SPH simulation in the zoom-region of the box. This resulted in the gas overdensitythreshold ranging from 0.25 to 20 times the dark mattermass in a given cell. We assumed that cold gas above thedensity threshold belonged to the multiphase interstellarmedium (ISM), which we could not resolve. As such weemployed a polytropic equation of state for such gas with T = T ( n/n ∗ ) κ − , with polytropic index κ = 4 /
3, andISM temperature T = 500 K. The density is also consis-tent with particle simulations of similar resolution (e.g.,McCarthy et al. 2010; Scannapieco et al. 2012 and refer-ences therein; Sijacki et al. 2012; Hayward et al. 2014).The implementation of star formation in RAMSES is de-scribed in Rasera & Teyssier (2006). We defined a unitstellar mass, m ∗ ,R = ∆ x n ∗ m p /X H , where m p is theproton mass, and X H = 0 .
76 is the hydrogen mass frac-tion. Thus the unit stellar mass is the mass in a cell at thethreshold density, and is therefore the minimum stellarmass. If a cell’s gas density was above the star-formationthreshold density and overdensity then we used Equa-tion (1) with M g = ∆ x ρ to determine the amount ofstellar mass expected to be created in the next time step.Comparing this expected mass with the unit stellar mass Fig. 1.—
Star formation density history up to z = 3 in thezoom region, comparing HYDRA and
RAMSES . Black squares showthe global star formation rate density in the high-resolution zoomregion for
RAMSES while blue circles show the results for
HYDRA .Note the delayed star formation for
RAMSES runs near z = 4. Themaximum refinement level increases at this redshift, followed by astrong increase in star formation as we better resolve the centraldensity peak of star-forming clumps. The best agreement is above z ≈ z ≈
3, when the resolutions were most comparable. yielded an expected number of stars to be formed, whichwas used as the expectation value for a random integerdrawn from a Poisson distribution. A star particle wasthen generated with the same velocity as the cell and amass equal to the random number times the unit stel-lar mass, and the cell’s mass was reduced by the starparticle’s mass.The implementation of star formation in
HYDRA is de-scribed in detail in Thacker & Couchman (2000). Eachgas particle accumulated a stellar component followingEquation (1), where M was that gas particle’s mass, and ρ was the local gas density at the particle. Once theaccumulated mass was equal to the star particle mass,which we set to be half the high-resolution gas particlemass, a star particle was formed, and the gas particle’smass was equivalently reduced. Until the time at which astar was formed, the gas dynamics used the total gas par-ticle mass, not just the non-stellar component. The onlyexception to this was that the mass in Equation (1) onlyincluded the non-stellar component of the gas particle.Once a gas particle that already made one star particlehad converted 80% of its remaining mass into stars, thenthat entire gas particle was turned into a second starparticle.In Figure 1 we compare the history of star formationfor the two codes in a 25 h − comoving Mpc sphere cen-tered on the region of interest. The varying resolution ofthe RAMSES runs is manifested by the delay in star forma-tion shortly before an increase in maximum refinementlevel, occurring at roughly z = 9 and z = 4. While for h min = 2∆ x there was fairly good agreement before usingthe overdensity threshold discussed above, the inclusionof this threshold results in stars forming only in galaxycores.We assumed 10% of the stellar mass formed was contained in high-mass stars that contributed feedbackthrough type II supernovae (SNe), with 10 erg of en-ergy released for 10 M (cid:12) of high-mass stars (i.e. per 100M (cid:12) of total stars). We thus used a total efficiency of (cid:15) fb = 5 × erg g − (Sommer-Larsen et al. 1999).This energy was deposited immediately into the vicinityof the formed star particle. Although the typical lifetimeof a high-mass star was on order tens of time-steps, oursingle particle is representative of several giant molecu-lar clouds, and we expect little issue to arise from inject-ing this energy without a delay at those kpc resolutions(see Wurster & Thacker 2013). The implementation ofthe stellar feedback in RAMSES is described in Dubois &Teyssier (2008). This energy was deposited kinemati-cally with a radius of a single cell around the star par-ticle. The injected mass, momentum, and energy areconsistent with a Sedov blast wave solution. The imple-mentation of the stellar feedback in
HYDRA is described inThacker & Couchman (2000). Since particle positionsare not isotropic, we injected this energy thermally, usinga kernel weighting for gas particles within the star par-ticle’s former smoothing length. Post-feedback coolinguses a multiphase description of the gas density, whichis much lower than the cold gas. In this model, pres-sure equilibrium is assumed, and the resulting decreasein density is determined by the net energy increase. Inthis way, the method is similar to other delayed coolingschemes (see for example Gerritsen & Icke 1997). Thismultiphase density then gradually returns to the particledensity, set by a half-life time of 1 Myr, consistent withthe time for the blast wave to reach the cooling radius fordensities near our star-formation threshold and feedbackenergy (Blondin et al. 1998).
Black Holes & AGN Feedback
Black holes (BHs) were modeled as sink particles inboth types of simulations, and they were formed wherethe local gas and stellar density were both above the starformation criteria, n ∗ = 0 .
05 cm − . All BHs had a seedmass of 8 × M (cid:12) , and to ensure only one BH was madeper galaxy, they were only allowed to form in locationsat least 30 comoving kpc from all other BHs. In RAMSES ,this particle was given the same momentum as its cell oforigin, and the cell’s mass was reduced by the seed mass.In
HYDRA , the BHs were spawned at the same locationand with the same momentum as the source gas particle,and the gas mass was reduced by the seed mass. Each BHparticle in
HYDRA had a smoothing length that overlappedwith roughly 60 neighboring gas particles. Black holescould merge when they were within 4 resolution units inAMR, or two smoothing lengths in SPH.The black holes accreted gas following the Bondi-Hoyle-Littleton rate (Bondi 1952), as described inDubois et al. (2013) and Wurster & Thacker (2013).The accretion rate is given by dM BH dt = 4 πα G M ¯ ρ (¯ c + ¯ u ) / , (2)where M BH is the black hole mass, ¯ ρ is the local averagegas density, ¯ c s is the local average sound speed, ¯ u is thelocal average gas speed, and α is a dimensionless boostfactor with α = max[1,( n/n ∗ ) ] (Booth & Schaye 2009),which accounts for our inability to resolve the cold high TABLE 1Simulation Summary
Runs Code Boxsize ∆x h min2 M DM3 M g3 M ∗ AMR-
RAMSES
100 550 - 80 - 0.2
SPH-
HYDRA
100 - 1100 80 16 8
Runs Cool/Reion. n ∗ (cm − ) c sf (cid:15) fb (erg g − ) (cid:15) r (cid:15) f M BH , s ( M (cid:12) ) -NC N - - - - - - -FID
Y 0.05 0.1 5 × - - - -QSO Y 0.05 0.1 5 × × comoving Mpc h − , physical pc, M (cid:12) density ISM gas around the BH. We set the maximumaccretion rate to be the Eddington accretion rate, dM Edd dt = 4 π GM BH m p (cid:15) r σ t c , (3)where σ t is the Thompson cross section, c is the speedof light, and (cid:15) r is the radiative efficiency, set to 0.1 forthe Shakura & Sunyaev (1973) model of accretion ontoa Schwarzschild BH. In RAMSES , each BH particle has acloud of sensors that sample the surrounding gas to de-termine the average gas quantities. In
HYDRA we useda kernel-weighted average of the gas particles that over-lap with the BHs smoothing length. Since mass is dis-cretized in
HYDRA , we used an internal and a dynamicBH mass. The internal mass was incremented by the ac-creted amount, and once this mass was larger than thedynamical mass by half a gas particle mass, the closestgas particle was accreted and the dynamical mass was in-creased by this amount. This dynamical mass was usedto calculate the gravitational effect of the BH, while theinternal mass was used to determine accretion and feed-back properties.The BH particles were advected similar to the darkmatter particles. To model the effect of gas on theBHs, we included a drag force, thus avoiding spuri-ous oscillations of the BHs about their local potentialminimum. This dynamical friction is set to F DF = f gas παρ ( GM BH / ¯ c s ) , where f gas is a factor whose valueis depends on the local Mach number, with 0 < f gas < HYDRA , the dragforce is calculated using the smoothed gas values of theBH’s neighboring particles.In this work, we have employed the quasar-mode AGNfeedback following Dubois et al. (2013). At every stepthe thermal quasar feedback mode injected energy into asphere of radius ∆ x centered on the BH, at an injectionrate of ˙ E AGN = (cid:15) f (cid:15) r ˙ M BH c , where (cid:15) f is a free parameterset to 0.15 to reproduce the M BH − M b , M BH − σ b , andBH density in the local universe (see Dubois et al. 2012).The implementation of this feedback mode in HYDRA wasbuilt upon the existing work by Wurster & Thacker(2013), which heats gas particles within two smoothinglengths. For both methods, the feedback was depositedwith a uniform specific energy inside the bubble radius.For each simulation code, we ran three simulations (seeTable 1). One simulation was carried out without cool-ing, reionization, star formation, or AGN, so as to bestunderstand, from the ground up, how the two meth-ods compare. This simulation is labelled
AMR-NC and
SPH-NC (for no cooling) in
RAMSES and
HYDRA , respec-tively. A second simulation includes cooling, reioniza-tion, star formation and SNe feedback. This simulation is labelled
AMR-FID and
SPH-FID (for fiducial) in
RAMSES and
HYDRA , respectively. A final simulation includes theseprocesses, but additionally tracks the evolution of AGNand their associated feedback. This is labelled
AMR-QSO and
SPH-QSO in RAMSES and
HYDRA , respectively.For the
QSO run, to make the two implementations assimilar as possible, we used a temperature ceiling of 10 K, and in the SPH run, energy was deposited kernel-weighted to the gas particles within the BH smoothinglength. Regardless of the simulation method, if the AGNfeedback energy at a given step would heat its environ-ment above the temperature ceiling, then the excess en-ergy was saved for the following step, and accretion wasstalled until this excess energy was administered. Withthis work we wish to highlight the dependence of clusterenvironment and galaxy gas evolution on feedback imple-mentations and how this evolution is dependent on thenumerical method. RESULTS
We break up our results into four main sections, as wefocus on the cluster environment and its gas, and how thecluster responds to the addition of increasingly complexgas physics in AMR and SPH simulations. We beginby looking at the characteristics of the gas in the clusterhalo and sequentially increase the physics included in thesimulations. A summary of the different simulations isgiven in Table 1. We first discuss the Non-Cooling runs,and then compare these simulations with our Fiducialruns, which include radiative cooling, star formation, andstellar feedback. We then introduce our AGN simulationsand compare these with our fiducial runs. For these firstthree sections we focus in great detail at snapshots at z =5 and z = 3 , and we end our results section by lookingat the continuous evolution of the cluster gas, stellar andblack hole components between these two snapshots. Non-cooling Runs
We first ran the two NC simulations, which did not in-clude gas cooling, star formation, reionization, or blackholes. These two simulations are crucial for ensuringthe gravity solvers are in good agreement and that theartificial viscosity implemented in HYDRA produces con-sisten shock heating. They also allow us to compare oursimulations directly with other non-radiative work (e.g.,Voit et al. 2005; Mitchell et al. 2009). Unlike the morecomplex cases, we only ran these simulations to z = 5,which given their simplicity was sufficiently advanced tobe able to compare with the more complex runs. Virialquantities, taken at the radius within which the averagedensity is 200 times denser than the critical density, arepresented for all simulations in Table 2.The left column of Figure 2 compares the density of TABLE 2Cluster Characteristics z Run M tot ( M (cid:12) ) r vir (kpc) T vir (K) S vira (keV cm ) (cid:16) M bar M tot (cid:17) / (cid:16) Ω b Ω m (cid:17) AMR-NC . × . ×
17 0.92
AMR-FID . × . ×
17 1.1
AMR-QSO . × . ×
15 0.68
SPH-NC . × . ×
18 0.92
SPH-FID . × . ×
18 0.92
SPH-QSO . × . ×
17 0.684
AMR-FID . × . ×
44 1.0
AMR-QSO . × . ×
40 0.61
SPH-FID . × . ×
44 0.92
SPH-QSO . × . ×
41 0.683
AMR-FID . × . ×
128 1.0
AMR-QSO . × . ×
123 0.74
SPH-FID . × . ×
126 0.92
SPH-QSO . × . ×
123 0.68 a Following Voit et al. (2005), but per H atom instead of electrons. -29 -27 -25 -23 logDensity (g cm − ) logTemperature (K) logTemperature (K) (cid:18)(cid:17)(cid:19)(cid:16)(cid:18)(cid:17)(cid:19)(cid:16) (cid:18)(cid:17)(cid:19)(cid:16)(cid:18)(cid:17)(cid:19)(cid:16) (cid:18)(cid:17)(cid:19)(cid:16)(cid:18)(cid:17)(cid:19)(cid:16) z-Proj z-Proj z-Slice z =5 A M R - N C S P H - N C Fig. 2.—
Comparisons of density and temperature for
AMR-NC (top) and
SPH-NC (bottom) out to 8 virial radii ( r vir = 70 physical kpc,indicated by a white circle) at z = 5. Each image is thus 1.1 physical Mpc across centered on the halo. Left (Middle): Projections along the z -axis showing the gas density (density-weighted temperature). Right: Slices of the same region for temperature highlighting the positionof the virial shock. the gas out to 8 virial radii surrounding the cluster at z = 5 from the two simulations (top: AMR, bottom:SPH). Unsurprisingly, there is significant correspondencebetween the RAMSES and the
HYDRA plots. However we note significantly more dense gas present in the SPHrun. We attribute this to a combination of the differ-ent gravitational calculations on small scales betweenthe two methods, and numerical dissipation of entropyin dense regions in SPH. On small scales
HYDRA usesa particle-particle gravity solver, while
RAMSES uses thestandard adaptive particle mesh. Particle-particle meth-ods have better short-range force resolution, which givesmore clumps (e.g., O’Shea et al. 2005). Additionally,
HYDRA uses an energy-conserving, rather than entropy-conserving, implementation of SPH. Therefore, in denseregions some entropy is dissipated numerically, furtherdifferentiating these clumps in SPH. The larger numberof dense knots in classic SPH is well established, and hasbeen discussed in detail (e.g., Frenk et al. 1999; Kauf-mann et al. 2006; Power et al. 2014).The middle and right columns of Figure 2 show thetemperature in projection and slice, respectively, of thegas for the same region. The pixelation seen in the SPHslice is a result of the mapping used to produce the im-age with the same visualization software as for the AMRimages (Turk et al. 2011, http://yt-project.org/). Inprojection we see better agreement between the two sim-ulations, while the slices highlight the different extentof the virial shocks. The SPH run has a larger shockradius that is spread over a wider spatial extent, con-sistent with the artificial viscosity injecting entropy at asomewhat larger radii than in AMR. This earlier shockheating causes more high-entropy gas in SPH since thedensities are very similar at the virial radius. This in-creased entropy may also be related to the poor abilityof SPH to model subsonic turbulence, as described inBauer & Springel (2012). Gas accretion into clustersis dependent on subsonic gas, which in AMR correctlycascades to smaller scales but in SPH thermalizes nearthe driving scale. It is this, combined with heating bythe artificial viscosity, that leads to an expanded virialshock.Figure 3 presents radial profiles of the gas density andentropy, here expressed as S ≡ k B T /n / , where k B isthe Boltzmann constant, and n H is the hydrogen num-ber density, for all gas, out to eight times the virial ra-dius ( r vir = 70 physical kpc). Virial quantities are anno-tated with dotted grey lines, to aid in comparisons withworks that look at multiple halos, switching to scale-freeunits normalized by these virial quantities (e.g, Voit etal. 2005; Mitchell et al. 2009; Kereˇs et al. 2012). Also in-cluded in the left plots of this figure are volume-averageddensity and entropy (taken as the volume-averaged tem-perature divided by the volume-averaged density to the2/3 power), which for direct comparison between AMRand SPH, we normalize by the virial quantities.For 0.2 r vir < r < r vir , the average entropy and den-sity profiles agree very well, and they are consistent withprevious comparisons looking at a range of mass scales(Frenk et al. 1999; Voit et al. 2005; Mitchell et al. 2009;Richardson et al. 2013; Sembolini et al. 2015). In partic-ular, the entropy profile matches the observations of Voitet al. (2015), where S ∝ r . . The average AMR gas isslightly more dense through the virial shock, which maybe due to SPH smoothing out the shock width. Out-side of the virial radius, the average density is again inagreement between the two methods, although the en-tropy is slightly increased in SPH near the virial radius.Sembolini et al. (2015) saw a similar small increase inthe entropy at the virial radius in their HYDRA simula-tion compared to their AMR comparison run. As we discussed above, this is due to higher temperature at thevirial radius.Within 0.2 r vir , on the other hand, the SPH gas den-sity is higher than in AMR, and the entropy reaches alower value than AMR before plateauing. The higher,more centrally peaked density in SPH is seen in Figure 2in the center of the cluster. This is also consistent withprevious comparisons (Frenk et al. 1999; Voit et al. 2005;Mitchell et al. 2009; Richardson et al. 2013; Semboliniet al. 2015), and Mitchell et al. (2009) demonstratedthat the lower entropy in the central region of SPH sim-ulations is due to a reduced amount of mixing, with thelow-entropy particles instead sinking to the center. Itis unclear, however, just how much over-mixing occursin AMR simulations, thus the true central profile is ex-pected to be slightly lower than in AMR. However, thedensity discrepancy is not sufficient to explain the en-tropy difference, thus the SPH gas is slightly cooler inthe center.Looking at the full distribution of density beyond thevirial radius, we see that the density in HYDRA extendsto higher values than in
RAMSES , consistent with satelliteobjects also collapsing to slightly higher densities in SPH,and visible in Figure 2. The entropy distribution has awider distribution of entropy in the center of the SPHsimulations, consistent with less mixing of high and lowentropy particles.Figure 4 shows two-dimensional entropy-density andtemperature-density distributions functions for gaswithin the virial radius, along with one-dimensional dis-tribution functions for each of these quantity. In theSPH simulation, high-density gas is found at lower en-tropy and even higher densities than in the AMR case(compare with Figure 3). This gas is at the center of thecluster, where the entropy is numerically dissipated, andnot accurately reinjected through mixing (Mitchell et al.2009) or from large scales (Bauer & Springel 2011). Thislarger fraction of high-density, low-entropy gas in SPHis also seen for halo substructure in the nIFTy compari-son (Sembolini et al. 2015A). A larger fraction of gas isfound at very high entropy in SPH, corresponding to lowdensity, hotter gas at the virial radius, as we discussedabove. Besides the very high density gas in SPH, the dis-tribution of density and temperature in SPH and AMRare very consistent.In general, the subtle differences that appear betweenthe simulations are not surprising due to the difficultyin implementing an ad hoc artificial viscosity in SPH,the energy-conserving implementation of the fluid equa-tions in SPH, and the tendency for SPH to undermix andAMR to overmix. We proceed, aware that these smalldifferences may compound when including cooling, starformation, and feedback.
Fiducial Runs
Next we look at the
FID simulations, which compareAMR and SPH with the inclusion of cooling, reioniza-tion, star formation, and stellar feedback. We first beginby comparing these results with the z = 5 NC results, andthen compare the FID results in more detail at z = 3.The halo virial quantities are listed in Table 2, whichshows that the inclusion of cooling has led to a highergas fraction than in the AMR-NC simulations, with the
AMR-FID values now slightly exceeding the cosmic aver-
Fig. 3.—
Radial profile plots at z = 5 showing the gas entropy (top) and gas density (bottom) vs radius for AMR-NC (middle) and
SPH-NC (right) out to 8 virial radii from the cluster. 0.2, 0.7, and 1.0 virial radii, and the virial density and entropy are indicated by thegrey dotted lines (see Table 2). The virial entropy is calculated following Voit et al. 2005, except we use the hydrogen density instead ofthe electron density. Color corresponds to log gas mass probability distribution fraction, where the integral of this quantity over a givenplot is unity. The solid black (dashed blue) lines in the middle (right) plots demark the volume-weighted average density and entropy(taken as the average temperature over average density to the 2/3 power) with radius for all gas in the AMR (SPH) simulation. The leftplots compare the average values of AMR vs SPH, using the same color and line scheme as the middle and right plots, scaled by the virialquantities so that they are directly comparable.
Fig. 4.—
Phase plots at z = 5 showing the gas entropy (top)and temperature (bottom) vs density for AMR-NC (left) and
SPH-NC (right) within the virial radius. Color corresponds to log gasmass fraction. The horizontal dotted line demarks the virial en-tropy and temperature of the cluster. On the left edge in linearunits are one-dimensional probability distributions of the log en-tropy (top) and log temperature (bottom) comparing the relativedistribution of gas mass for both AMR (black solid line) and SPH(blue dashed line), while the bottom edge shows in linear units theone-dimensional probability distribution of log gas density, showntwice to facilitate comparison with the above phase plots. age.In Figure 5 we show projections of the gas density for the
AMR-FID and
SPH-FID simulations at z = 5 . In the
FID simulations the gas can cool, resulting in more con-densed structures. Thus the filaments are thinner, andthe galaxies are collapsed to thin disks. On these scales indensity we see little impact of star formation or feedback.In comparison, we see an amplification of the differencesseen in Figure 2, with more clumps in SPH than in AMR.In AMR, the gas filaments are smooth ribbons of nearuniform density gas with large galaxies residing withintheir nodes. In SPH, these filaments are much more in-homogeneous, housing many more small clumps that areless dense than AMR galaxies, but more dense than thesurrounding filament gas. The clumpiness in
SPH-NC isnow compounded in
SPH-FID by the fact that this denser,lower-entropy gas has shorter cooling times, leading toquicker fragmentation times for the filament as a whole.Inside the virial radius the dense gas is completely frag-mented into individual parcels. The larger clumps agreebetween AMR and SPH, and are cospatial with clumpsin dark matter. However, the additional clumps found inSPH, which are seen in other studies of classic SPH (e.g.,Frenk et al. 1999; Kaufmann et al. 2006; Power et al.2014), are not cospatial with dark matter clumps and aredue to artificial dissipation. We stress however that asthe densest clumps fragment along the polytrope, theyare forced to increase in temperature. Now that the fill-ing factor of dense gas has decreased, we expect that theaddition of AGN feedback will be more efficient in SPHand better able to blow away the more tenuous ambientgas, as we discuss in detail in § -29 -27 -25 -23 logDensity (g cm − ) logTemperature (K) logTemperature (K) (cid:18)(cid:17)(cid:19)(cid:16)(cid:18)(cid:17)(cid:19)(cid:16) (cid:18)(cid:17)(cid:19)(cid:16)(cid:18)(cid:17)(cid:19)(cid:16) (cid:18)(cid:17)(cid:19)(cid:16)(cid:18)(cid:17)(cid:19)(cid:16) z-Proj z-Proj z-Slice z =5
10 5 0 5 10 x (kpc) y ( k p c ) T e m p e r a t u r e ( K )
10 5 0 5 10 x (kpc) y ( k p c ) -26 -25 -24 -23 -22 -21 -20 D e n s i t y ‡ g c m ·
10 5 0 5 10 x (kpc) y ( k p c ) T e m p e r a t u r e ( K )
10 5 0 5 10 x (kpc) y ( k p c ) -26 -25 -24 -23 -22 -21 -20 D e n s i t y ‡ g c m · A M R - F I D S P H - F I D Fig. 5.—
Same as Figure 2 but for the
FID simulations, measuring 1.1 physical Mpc across. The inserts in the projection plots are slicesof density (left) and temperature (middle) focused on the galaxy scale, measuring 28 physical kpc across, and reach out to 0.2 r vir . These differences are visibile in the temperature projec-tions and slices shown in the middle and right panels ofFigure 5, and are much more apparent than in the NCcomparison runs shown in Figure 2. Thus the increasedpost-shock temperature and slightly wider shock in SPH(e.g., Hubber et al. 2013) leads to a lower cooling rate.We discuss the expected cooling times below.Figure 6 presents the profile diagram of the gas at z = 5out to 8 virial radii, to compare with Figure 3. Here wehave split the average quantities into a hot componentabove 50,000 K, and a total component. The hot gasis more directly relatable to other work (e.g., Kereˇs etal. 2012) and with our NC runs, and less susceptibleto the different clumping behavior. Cooling in the fidu-cial case also leads to a large population of gas at verylow entropies and high densities, which is seen in boththe AMR and SPH simulations. However in the AMRrun there is a two-phase medium, with cold and densematerial found at the same radial distances as the hot,tenuous gas. In the SPH run, on the other hand, thehot and cold phases are segregated, with cold materialfound almost exclusively in the center, surrounded by ahot diffuse region. Thus, the addition of cooling has am-plified the ability of low entropy gas in SPH to sink to the center of the halo, and the two methods yield verydifferent results within 0.2 r vir .Beyond this radius, out to 0.7 r vir , there is better agree-ment of the hot gas between the two methods, whileAMR is better able to model cold streams with entropyof roughly 10 − keV cm . Beyond 0.7 r vir and out toroughly the virial radius, the entropy differences are evenstronger than in the NC runs. Here the post-shock gascan efficiently cool in the AMR run, leading to the shockradius lying within the virial radius. In the SPH run, onthe other hand, the shock radius is at the virial radius.This is not what is expected from theory. At z = 5,this system has a dynamical time of roughly 0.11 H − while the cooling time is approximately 0.03 H − . Thusthe post-shock gas should cool quicker than the typicalgrowth time of the halo, and we expect the virial shock tolie within the virial radius (e.g. White & Rees 1978; Birn-boim & Dekel 2003). In our HYDRA run, no star particleshave been made within the virial radius of this clusterby z = 5, thus this larger virial shock and smaller fillingfactor of cold gas is not due to a difference in the SPHSNe feedback. Instead, this appears to be due to theexcess heating near the virial radius seen in the SPH-NC run, coupled to a reduced cooling rate as the gas shock is1
Fig. 6.—
Same as Figure 3, but for the
FID simulations. The red (magenta) lines in the middle (right) plots demark the volume weightedaverage density with radius for gas hotter than 50,000 K in the AMR (SPH) simulation, respectively. The left plots compare the averagevalues, scaled by the virial values, of AMR vs SPH, using the same color and line scheme as the middle and right plots. broadened in SPH. This lower cooling at the virial radiusin SPH then leads to an inflated post-shock region. Thisis similar to the comparison of the projected tempera-ture of a galaxy halo taken from the SPH code
GADGET with the moving mesh code
AREPO reported in Nelson etal. (2013), in which the SPH virial shock was located atlarger radii than in the moving-mesh case.
Fig. 7.—
Same as in Figure 4, but for the
FID simulations. Thevertical dotted line marks the star-formation density threshold atroughly 10 − g cm − . Phase diagrams of the gas at z = 5 out to the virialradii are given in Figure 7. These plots show many newfeatures not seen in the the NC results shown in Fig-ure 4. The medium is now found mostly at 10,000 - 20,000 K, where the cooling function has a local mini-mum. As gas slowly cools through this regime, it be-comes denser and therefore has lower entropy. At densi-ties above 3 × − g cm − the gas cools more efficiently,until cooling to the enforced polytope, where the tem-perature is forced to scale with ρ / . Finally, post-shockgas is heated to just above the virial temperature, andcools inefficiently, except at high densities. The SPH-NC gas at high temperatures that extends to a densityof ρ (cid:39) − g cm − is able to cool in the SPH-FID simulation, forming lines extending down to the cooler, T (cid:39) K regime. Finally, given that gas was foundat slightly higher densities in
SPH-NC than in
AMR-NC ,we naively expected the cooling rate for the
SPH-FID gasto be faster. However, since the polytrope gas extendsto higher densities and temperatures in AMR, SPH ap-pears to prevent gas in the polytrope from moving tohigher densities. This may be due to undermixing, as dis-cussed before. Additionally, AMR creates star particlesof smaller mass than in SPH, thus AMR is more quicklyremoving pressure support from the densest regions, pos-sibly leading to the buildup of more high-density poly-trope gas. A future study of star particle mass and starformation rate is needed to better understand this effect.Next we carry out the same analysis at z = 3 , whichwe will also use to compare with the QSO runs, in theredshift regime where AGN feedback becomes more im-portant. Density and temperature projections of the gasat this redshift are given in Figure 8. In density, thevarious accretion filaments have coalesced into a morecondensed cosmic web that is less volume filling, withdenser gas in the cosmic nodes. The temperature pro-jections also show an increased virial temperature, thatextends to a larger radius, and the filaments are alsoencased within post-shock heated gas. The virial shockradius in SPH and AMR are both now at or beyondthevirial radius at this redshift, since both the cooling time2 -29 -27 -25 -23 logDensity (g cm − ) logTemperature (K) logTemperature (K) (cid:18)(cid:17)(cid:19)(cid:16)(cid:18)(cid:17)(cid:19)(cid:16) (cid:18)(cid:17)(cid:19)(cid:16)(cid:18)(cid:17)(cid:19)(cid:16) (cid:18)(cid:17)(cid:19)(cid:16)(cid:18)(cid:17)(cid:19)(cid:16) z-Proj z-Proj z-Slice z =3
40 20 0 20 40 x (kpc) y ( k p c ) T e m p e r a t u r e ( K )
40 20 0 20 40 x (kpc) y ( k p c ) -26 -25 -24 -23 -22 -21 -20 D e n s i t y ‡ g c m ·
40 20 0 20 40 x (kpc) y ( k p c ) T e m p e r a t u r e ( K )
40 20 0 20 40 x (kpc) y ( k p c ) -26 -25 -24 -23 -22 -21 -20 D e n s i t y ‡ g c m · A M R - F I D S P H - F I D Fig. 8.—
Same as Figure 5 but for z = 3, where now r vir = 220 physical kpc. The images thus measure about 3.5 Mpc across centeredon the halo. and the dynamical time are roughly 0.11 H − . In AMRthe filaments clearly remain relatively cold, while this ismuch more difficult to see in SPH. Looking at the tem-perature slices (right), indeed the filaments inside theshock-heated gas are cool in SPH, but they are com-prised of smaller, more fragmented, gas clumps than inAMR, while the post shock gas is much more extended.The discrepant behavior, that is the overheating and un-dercooling at the shock, that explained the hotter gas inSPH at z = 5 at the virial radius has had a runawayeffect, and by z = 3 its impact is even more extreme. Inthe galaxy-scale inset, it is clear that there is cold gas inboth AMR and SPH, within which stars form, but thisgas is more ordered in AMR, consistent with a higherspecific angular momentum.In the profile plots, presented in Figure 9, we see thesame qualitative behavior in entropy and density forAMR and SPH as at z = 5. A similar trend in theentropy of the hot gas is visible in the volume-averageprofiles, in particular at small radii ( r < . r vir ) theSPH entropy is lower than in AMR, at intermediate radii(0 . r vir < r < r vir ) the SPH entropy is larger than inAMR, at r (cid:39) r vir the entropies are in good agreement,and out to 2 r vir the SPH entropy is again higher, with temperatures near the virial temperature. However, theratio of AMR and SPH entropies are larger than theirvalues at z = 5 , thus the behavior of the two codes iseven less consistent. The densities, on the other hand,show the same level of consistency between the two codesat z = 3 as at z = 5.Thus at both z = 5 and z = 3 with radiative cooling,star formation, and reionization, the overall trend is forstandard SPH to have lower entropy cluster cores andlarger entropy at larger radii compared with AMR. Thisis consistent with other comparisons between standardSPH and AMR for non-radiative simulations (e.g., Voitet al. 2005; Mitchell et al. 2009), but oddly, this is notwhat is seen in comparisons between standard SPH andthe moving mesh code AREPO . In fact, Kereˇs et al. (2012)saw that for intermediate mass halos similar to our halothat the moving mesh simulations have even lower en-tropy values in the core and higher entropy values nearthe virial radius than the SPH runs. The authors arguedthat the moving mesh was capturing a cooling flow, andthus SPH was not capturing sufficient cooling in the gas.While this may be the case, the authors also discuss themixing of low entropy gas at the center of the clusterin
AREPO which is not captured in
GADGET . Yet Mitchell3
Fig. 9.—
Same as Figure 6, but for z = 3. Fig. 10.—
Same as in Figure 7, but for z = 3. et al. (2009) compared SPH to Eulerian simulations toshow that such mixing injects heat into the central re-gion, leading to a higher, not lower entropy value. WhileMitchell et al. (2009) did not include cooling, we still findthe same behavior in our AMR simulations with cooling,suggesting that dissipative heating is in fact sufficient tooffset the cooling in the center. This is even with ourassumed constant metallicity of one third solar, and anoverestimate of the cooling from H and He (see § S is above 100 keV cm ; Oh & Benson2003). In AMR, this gas instead only reaches S = 10 keV cm . The temperature profiles again are in good agree-ment, although the post-virial shock gas does still ex-tended to slightly hotter temperatures. Thus SPH locksup more gas in the diffuse, high-entropy phase, whichwill have an impact on the amount of star-forming gas inthe cluster. Finally, in SPH there is a small feature ex-tending from the star-forming ISM to hotter, denser gas.This is a post-SNe feedback region, where gas is instantlymoved to higher temperatures in the phase diagram, andthen has its cooling artificially delayed. This results inthe heated gas first expanding adiabatically, dropping indensity and temperature, until as its cooling ramps up itcools more quickly and its entropy drops. AGN Feedback Runs
We now look at the
QSO simulations to see how theinclusion of AGN feedback impacts the halo and galaxygas, and how AMR and SPH compare in this context.We first begin by comparing these results with the z = 5 NC and FID results, followed by the
FID results at z = 3 . The halo virial quantities for these simulations arelisted in Table 2. The inclusion of AGN feedback hasled to a much lower gas fraction than in the
FID simu-lations, which is now consistent between the AMR andSPH simulations. As we show below, with AGN feedbackthe virial shock is now beyond the virial radius in bothAMR and SPH, and it is for this reason that the bary-onic fraction has dropped, as gas spends a longer time inthe outer halo before cooling and falling within the virialradius.In the left column of Figure 11, we show large-scale(out to 8 r vir ) projections of the gas density for the AMR-QSO (top) and
SPH-QSO (bottom) simulations at z = 5 , to compare with Figures 2 and 5. The introduction ofAGN feedback results in very little apparent impact onthe density distribution at this redshift. This is mostlydue to only a 30% drop in overall gas fraction, which isdifficult to see with the scale spanning 6 orders of mag-4 -29 -27 -25 -23 logDensity (g cm − ) logTemperature (K) logTemperature (K) (cid:18)(cid:17)(cid:19)(cid:16)(cid:18)(cid:17)(cid:19)(cid:16) (cid:18)(cid:17)(cid:19)(cid:16)(cid:18)(cid:17)(cid:19)(cid:16) (cid:18)(cid:17)(cid:19)(cid:16)(cid:18)(cid:17)(cid:19)(cid:16) z-Proj z-Proj z-Slice z =5
10 5 0 5 10 x (kpc) y ( k p c ) T e m p e r a t u r e ( K )
10 5 0 5 10 x (kpc) y ( k p c ) -26 -25 -24 -23 -22 -21 -20 D e n s i t y ‡ g c m ·
10 5 0 5 10 x (kpc) y ( k p c ) T e m p e r a t u r e ( K )
10 5 0 5 10 x (kpc) y ( k p c ) -26 -25 -24 -23 -22 -21 -20 D e n s i t y ‡ g c m · A M R - Q S O S P H - Q S O Fig. 11.—
Same as Figure 5 but for the
QSO simulations. nitude. However, under close comparison with Figure 5,we see that the densest regions in the
QSO runs havethe largest decrement in density compared with the
FID runs, while the filaments inside the cluster are somewhatdenser. In AMR we see that the filaments are slightlypushed off-center, and are thinner. In SPH, there is alower number of collapsed clumps, suggesting that thefeedback is offsetting the numerical dissipation. Finally,the diffuse halo medium inside the virial radius is slightlydenser with AGN. Thus, while by z = 5 AGN feedbackis not effective in moving large amounts of matter out ofthe cluster environment, it does reduce the central dens-est peaks and delay gas accretion into the halo.The middle and right columns of Figure 11 present pro-jections and slices of the gas temperature, respectively.While AGN feedback does not affect the gas density, itclearly affects the gas temperature. In both simulationtypes, the AGN feedback results in a larger volume of hotgas than in the FID runs, with an increase in the hot gasvolume-filling factor of roughly 10, and a virial shock ra-dius well beyond the virial radius. However, note that inboth
FID and
QSO simulations the amount of heated gasis consistently greater for the SPH simulations. Giventhat more gas is heated in SPH, it is thus surprising thatthe temperature slices reveal that the two methods give a consistent picture of the halo gas (within the white cir-cle). The halo virial radii are 66 physical kpc and 70physical kpc for
AMR-QSO and
SPH-QSO , respectively,and within the halo the gas is nearly uniform at the virialtemperature of 8 × K. Although the AMR filamentsare more continuous than in SPH, similar to what is seenin the
FID simulations, the cold filaments in
AMR-QSO truncate at a larger radius than in
AMR-FID and pointslightly away from the center of the cluster. The fila-ments’ truncation radius is signified by a second shocklocated at the virial radius. Thus, although it is notvisible in the large-scale projections, the AGN is impact-ing filament gas on small scales in the AMR simulation,which is also seen in higher resolution runs in Dubois etal. (2013). We do not see a similar impact on the fila-ments in SPH however, and this is consistent with whatis seen in Di Matteo et al. (2012).We interpret the similarity of the AMR and SPH halogas as follows: although the injection scale of the AGNfeedback for the two methods is not identical, the AGNeventually heats the halo gas sufficiently to self-regulateits accretion, and this self-regulating gas configurationoccurs at a “critical entropy” (Oh & Benson 2003; Scan-napieco & Oh 2004). We find that this critical entropy islargely code independent. However, it is clear that SPH5
Fig. 12.—
Same as Figure 6 but for the
QSO simulations. causes more collateral gas heating, such that the gas atvery large radii is heated to 10 K in the process of in-creasing the halo gas to the same critical entropy. Thisbecomes much clearer at z = 3 , as discussed below.Beyond the halo, we see that the intergalactic mediumin the SPH simulation is significantly cooler than in theAMR simulation. This difference is unfortunately causedby an erroneous switch that turned off reionization heat-ing for gas below 500 K in our SPH-QSO simulation.However, the shock increases the temperature in AMRby three orders of magnitude, and therefore the fact thatthe IGM gas is artificially cooled results in little impacton the halo itself.The impact of AGN feedback on the gas profiles ofthe cluster environment at z = 5 are presented in Fig-ure 12. Even though the temperature projections showa strong signature of AGN feedback in both simulations,the average profiles in this figure are largely similar to thefiducial cases shown in Figure 6. However, in the QSO simulations, the distribution of gas at very high entropieshas increased, extending roughly an order of magnitudehigher than in Figure 6, with the largest increases occur-ring at very small radii ( r <
20 physical kpc) and verylarge radii ( r >
100 physical kpc). Additionally, the in-clusion of AGN has resulted in less low entropy gas at10-20 kpc in the AMR simulation, while in the SPH sim-ulation there is less low-entropy gas (
S < )at all radii, but more mid-entropy gas ( S (cid:39) )at 5-20 physical kpc. Finally, by comparing the QSO re-sults with the non-radiative profile given in Voit et al.(2005), we see that the inclusion of AGN feedback hasresulted in a higher average entropy value at all radii forAMR, and a higher average entropy for the central gasin the SPH case.AGN feedback acts to diminish the SPH and AMRcentral density values and average values, bringing thetwo methods into better agreement. However, the hotSPH gas extends to higher densities, indicative of a very
Fig. 13.—
Same as in Figure 7, but for the
QSO simulations. recent feedback event in which the gas has not yet had anopportunity to expand. This central, hot, feature is visi-ble in the insets shown in Figure 11, which also show thepresence of more cool gas in SPH, explaining the lowerentropy profile. Overall there is much better agreementbetween the two methods than in the
FID runs, wherethe AGN heating results in similarly increased entropyprofiles necessary for self-regulating their accretion.In Figure 13 we show the phase diagrams of the z = 5gas within the virial radius. The one-dimensional PDFplots of the entropy make it much clearer that theamount of high-entropy gas has increased by includ-ing AGN feedback. In AMR-QSO , much more gas has S ≥
100 keV cm than in AMR-FID , and this highentropy gas is found over a range of densities, from ρ = 3 × − g cm − to the star-formation threshold6 -29 -27 -25 -23 logDensity (g cm − ) logTemperature (K) logTemperature (K) (cid:18)(cid:17)(cid:19)(cid:16)(cid:18)(cid:17)(cid:19)(cid:16) (cid:18)(cid:17)(cid:19)(cid:16)(cid:18)(cid:17)(cid:19)(cid:16) (cid:18)(cid:17)(cid:19)(cid:16)(cid:18)(cid:17)(cid:19)(cid:16) z-Proj z-Proj z-Slice z =3
40 20 0 20 40 x (kpc) y ( k p c ) T e m p e r a t u r e ( K )
40 20 0 20 40 x (kpc) y ( k p c ) -26 -25 -24 -23 -22 -21 -20 D e n s i t y ‡ g c m ·
40 20 0 20 40 x (kpc) y ( k p c ) T e m p e r a t u r e ( K )
40 20 0 20 40 x (kpc) y ( k p c ) -26 -25 -24 -23 -22 -21 -20 D e n s i t y ‡ g c m · A M R - Q S O S P H - Q S O Fig. 14.—
Same as Figure 11 but for z = 3. density of 10 − g cm − . In SPH, there is a smaller in-crease in gas with S ≥
100 keV cm , and most of thisis at low densities ≈ × − g cm − . Thus, withinthe halo AGN are heating gas at a range of densities inAMR, preventing it from cooling to very high densities,while in SPH, AGN are more efficient at heating the sur-rounding diffuse gas. This increase in high entropy, hotgas is also shown in the one-dimensional PDF plots oftemperature, which illustrate that in the AMR case theAGN heat less gas but this gas is heated to higher en-tropies and temperatures than in the SPH case. This ispartially due to the way the feedback energy is injectedinto the simulation, with AMR distributing the energyto the neighboring cells of the BH particle, and SPH dis-tributing the energy over the BH particle’s smoothinglength.Finally, as shown in the profile plots, in SPH the AGNis causing significant heating of dense gas in the poly-trope. This gas occupies the high-density region of theSPH plots, where the very dense but also hot gas is un-physical. Yet it appears to be a long-lived feature. Thisoccurs when a sink particle resides in a cold, compact,dense clump (such as seen in the inset of Figure 5) butit is not sufficiently massive for Eddington-limited feed-back to overcome the clump’s binding energy. Instead, the AGN must grow while keeping the clump at a hottemperature, until finally dumping sufficient energy intothe clump to overcome gravity.Moving to z = 3, projections of gas density and pro-jections and slices of gas temperature out to 8 r vir areshown in Figure 14. We see the same small effects asat z = 5. Thus, even at later times when AGN have alarger impact, thermal feedback does not lead to large-scale redistribution of gas. The temperature projectionsnow present significantly more extended hot gas, which,if we consider all gas heated above the typical IGM tem-perature of 10 K, is more volume filling in SPH. How-ever, if we consider the region with projected tempera-ture above 10 . K, then we again see that there is fairlygood agreement between the two methods. This regionin SPH is only slightly larger and reaches temperaturesonly slightly hotter than in AMR. Thus we find thatthe methods produce a halo gas temperature distribu-tion that is much more consistent than in the z = 3 FID runs and the
QSO runs at z = 5. By z = 3 the AGNhas heated the gas to roughly the same entropy in bothmethods, the entropy required for self-regulation.This agreement on intermediate and large-scales aside,the inserts, showing the gas density and temperature inthe inner 0.2 r vir , do show a few differences. Thus while7 Fig. 15.—
Same as Figure 12, but at z = 3. the AGN feedback acts to make the bulk of the halomaterial consistent thermodynamically, on smaller scalesthe gas is impacted differently. In SPH there is a morecentrally collapsed, hot, gas clump, which is broader indensity than in temperature, and whose outskirts aredenser than its interior. This region is in the processof a feedback-driven expansion, leading to a shock frontwith higher density and therefore increased cooling. Thisexpansion phase is made clearer in the gas profiles shownin Figure 15. In SPH the gas density peaks outside of thecenter at roughly 4 physical kpc, where entropy drops.We have confirmed that this is not due to an improperchoice of halo center.On intermediate and large scales in Figure 15 there isagain better agreement between the two methods. Out-side of the central clump, the AGN has led to a muchlarger amount of gas at high entropy at large radii, com-pared with the FID simulations in Figure 9. Within thevirial radius, the average entropy and density profilesare even more consistent than in the
FID and z = 5 QSO runs. Given that in SPH there is a central clump withhigher density, we would expect to see lower entropy inthe core of the SPH halo. However, the gas has finallybeen heated sufficiently to drive an outflow, giving a cen-tral peak in the SPH entropy. We have looked at z = 3 . z = 3 is a recent event, and in general lower entropycores in SPH are more common.Finally, in Figure 16 we compare the phase plots at z = 3 for the cluster. We can see the expanding centralclump in the SPH diagram, where a large plume of isopy-cnic gas has continuously been heated via AGN feedbackuntil reaching a sufficiently high temperature that thepressure can unbind the clump. However, note that al-though this unphysical region appears quite large, it con-tains only a small fraction of the gas in the halo. Instead,a greater amount of gas is found near 10 − g cm − in Fig. 16.—
Same as Figure 13, but at z = 3. both methods. The AGN feedback has led to similaramounts of gas in the hot diffuse halo, although in AMRthere is more gas held at 10 K just below the star for-mation density threshold, while this gas is found abovethis threshold in SPH. This is due to the clumping na-ture of the cold gas in SPH, where hot gas is less ableto disrupt these clumps through instabilities. In AMR,on the other hand, the dense gas forms streams insteadof clumps, with larger surface area that is more suscep-tible to turbulence. Thus in AMR the AGN feedback isbetter able to dissolve and/or physically move the veryinner filamentary gas into the warmer, diffuse medium,consistent with what has been seen in previous work (e.g.,Dubois et al. 2013; Nelson et al. 2015).In summary, the inclusion of AGN has lead to a moreconsistent halo environment for the two methods. This8
Fig. 17.—
Top Left: Evolution of the total halo mass of the cluster. The
AMR-FID ( SPH-FID ) simulation is presented with a solid black(blue) line where we mark the measured values in an output by squares (circles). The
AMR-QSO ( SPH-QSO ) simulation is presented witha dashed green (red) line. The magenta dotted line gives an exponential mass law, M = M f e − α ( z − z f ) = 2 . × e − . z − , followingWechsler et al. (2002). Bottom Left: Evolution of the halo mass accretion histories, given by d ln M/dz . Line and point styles match theabove plot, with the magenta line giving a constant value of 1.15, consistent with the exponential parameterization for M ( z ). Middle:Evolution of the dark matter and baryon components. Right: Evolution of the gas fraction normalized by the cosmic mean baryon fraction. is suggestive of a self-regulation scenario, which becomesstronger at later times. The result is that the AGN heatsthe halo environment to sufficiently high entropies so thatit turns off further gas accretion. It is very interestingthat, even though the two methods cannot probe pre-cisely the same spatial scales, the resulting entropy isvery similar. On large scales in SPH, however, the re-sult of AGN feedback heating the halo gas to this self-regulatory temperature is that a much larger volume ofgas is also impacted. This would have interesting im-plications for surveys of the thermal Sunyaev-Zel’dovicheffect (e.g., Spacek et al 2016). Evolution of Gas, Stars and Black Holes
Halo Gas
Finally, we turn to a more detailed study of the halo’sevolution. In Figure 17 we summarize the mass historyof the halo, taken from a range of outputs from roughly z (cid:39) −
3. The total halo mass evolution is consistentwith pure exponential growth, similar to that seen fora range of halos in Wechsler et al. (2002). Using anoffset at z = 3, so that the leading constant is roughlythe mass of the halo at the end of our simulations, weget an equation of the form M tot ( z ) = M ( z = 3) e − α ( z − , (4)where we find M ( z = 3) = 2 . × M (cid:12) , and α = 1 . α value is consistent with the most massive objectsconsidered in Wechsler et al. (2002). By breaking upthe mass into components, we see the dark matter massis more consistent between the AMR and SPH runs thanthe baryons. The AMR-FID run has a higher baryon frac-tion than the cosmic mean at all times, but this fractiondrops with time. This is due to the virial temperatureincreasing as the cluster grows, leading to longer cool- ing times at further radii, slowly depleting the amountof baryons that fall inside the halo. In the
SPH-FID runwe see similar behavior, with the gas fraction droppingin time. However, the gas fraction starts out at roughlythe cosmic mean, and then drops more quickly than inAMR. This difference is consistent with a higher post-shock temperature in SPH that is at a larger radius thanin AMR. The inclusion of AGN feedback has little im-pact on the parameterization of the halo’s growth. Bylooking at the effective slope, we can see there is a bitmore spread in the
QSO runs, with a slightly shallower α at large redshifts that becomes steeper at lower redshifts.In the evolution of the halo mass, only the baryonsappear affected by AGN feedback, resulting in a dropof roughly 30% in the amount of baryons in the clus-ter. With AGN feedback, there is more scatter in theevolution of the baryon fraction, but the AMR and SPHresults are more consistent than in the FID runs. As wehave seen, the AGN leads to roughly the same halo gasentropy, thus better agreement in the baryon gas fractionis expected.The specific angular momentum (sAM), j , and spin pa-rameter, λ = j/ ( √ r vir v vir ) (e.g., Bullock et al. 2001),are presented in Figure 18, where v vir is the virial ve-locity. The sAM of dark matter in the outer halo( r > . r vir ) grows at a smooth rate over the course ofthe simulation. Except at very early times, the growthof sAM is unaffected by AGN feedback. By normalizingby the virial radius and velocity it is clear that λ for theouter halo is roughly constant in time with a value near0.034, consistent with tidal torque theory (e.g., White1984; Bett et al. 2007). Accreting dark matter slowlyexchanges its sAM with the halo, leading to a slowerrate of growth of sAM in the central 0 . r vir , correspond-ing to the galactic scale. Even on small scales the growthof sAM is not impacted by AGN feedback. As seen in9 Fig. 18.—
Top two rows: Evolution of the specific angular mo-mentum (top row) and spin parameter (second row) of the clus-ter dark matter (left) and gas (right) outside 0.1 r vir . Bottomtwo rows: Evolution of the specific angular momentum (third row)and spin parameter (bottom row) of the cluster dark matter (left)and gas (right) within 0.1 r vir . An approximate fit to the AMR-FID ( SPH-FID ) simulation is presented with solid black (blue) lineswhere we mark the measured values in an output by squares (cir-cles). An approximate fit to the
AMR-QSO ( SPH-QSO ) simulationis presented with dashed green (red) lines.
Danovich et al. (2015), the dark matter sAM on smallscales is roughly an order of magnitude below the accret-ing dark matter value.The gas sAM behaves quite differently. Before z = 4the sAM is roughly constant in the outer halo, until at z = 4 it has roughly the same sAM as dark matter, atwhich point it grows faster than the dark matter sAM.By z = 3 the gas sAM is a factor of 3 above that ofdark matter. At very early times, the lower gas fractionsseen previously has led to lower sAM in the QSO sim-ulations. This suggests that at early time much of theangular momentum generation is via accretion of highsAM gas. However, by z = 4 the outer gas sAM hasthe same value in all simulations, which suggests thata significant amount of the torque generating the angu-lar momentum is gravitational in nature, consistent withtidal torque theory. The spin parameter is also roughlyconstant for gas in the outer halo, staying a factor ofroughly two above the dark matter, except at z = 4.Considering we are presenting the angular momentum ofall gas, and not just cold gas, it is encouraging to seethe agreement with Danovich et al. (2015), who see that cold gas in the outer halo has spin parameters of roughlythree times that of dark matter. A more detailed studyof the angular momentum of the different gas phases isbeyond the scope of this work.At small radii, however, there are significant differ-ences between SPH-FID and
AMR-FID , where the SPHsAM is consistently lower than in AMR. While SPH ex-plicitly conserves angular momentum in the absence ofartificial viscosity, Okamoto et al. (2003) demonstratedthat due to its difficulty in modeling the layer betweendifferent fluids, standard SPH can transfer angular mo-mentum from the dense disk to the hot diffuse halo. Ad-ditionally, Kaufmann et al. (2007) showed that disksin low-resolution SPH simulation, of comparable resolu-tion to this and other cosmological simulations, unphysi-cally transport angular momentum to the outer medium.However, AMR does not explicitly conserve angular mo-mentum, its grid has a preferred direction, and it suffersfrom advection errors. Thus AMR also suffers from spu-rious angular momentum dissipation. Yet, we see herethat the central 0.1 r vir has an order of magnitude moreangular momentum in AMR than in SPH. Thus, the SPHresolution issues may be playing the dominant role. Ad-ditionally, it is possibly that the more collimated filamentstreams in AMR incur more sAM from tidal torques sinceit has a high quadrupole moment (e.g., Danovich 2015),which is also weakly supported at high z and large radii.The inclusion of AGN feedback also leads to markedlydifferent behavior in the inner region for the two meth-ods. In SPH, the inclusion of feedback has led to a com-bination of accretion of low angular momentum gas, andremoval of this low AM gas from the central region. Thisleads to a net increase in the central SPH sAM. However,in AMR the sAM drops when we include AGN feedback.The two mechanisms leading to an increased sAM in SPHare also operating in AMR. Thus the difference is thatgas heated via AGN feedback moves the cold stream fila-ments outwards, inhibiting gas accretion from these outerregions and transfering the filamentary gas’ angular mo-mentum through shocks into the outer halo. While thesAM is increasing for all simulations, the spin parame-ters are dropping. Thus the angular momentum buildupin the central region is not keeping pace with the growthof the cluster, regardless of feedback processes. Star formation
In Figure 19 we present the halo stellar mass and starformation rates for the four simulations containing stars.We fit the halo stellar mass functions with exponentialpower laws, similar to the total halo mass fits in § z = 3, itis always the case that the AMR simulations have morestars than in the SPH simulations. This is mostly duethe gas reaching higher densities in AMR, as shown inFigure 7. Since the star formation rate scales with ρ . ,the higher density gas in AMR leads to a much largerstellar masses. The inclusion of AGN feedback leads to asignificantly reduced amount of stars, due to removal oflow angular momentum gas, and heating of surroundinghalo gas. Unsurprisingly, AGN have the biggest impactin SPH, where there is an average increase in the centralangular momentum and more large-scale heating, whichreduces the star formation rate in surrounding structuresoutside the virial radius, and thus a lower amount of0 Fig. 19.—
Top: Evolution of the halo stellar mass of the clus-ter. The
AMR-FID ( SPH-FID ) simulation is presented with a solidblack (blue) line where we mark the measured values in a outputby squares (circles). The
AMR-QSO ( SPH-QSO ) simulation is pre-sented with a dashed green (red) line. Lines are fitted exponentiallaws of the same form as Equation (4). Bottom: Evolution of thehalo star formation rate, taken as the change in halo stellar massbetween adjacent outputs. Lines, given by the derivative of theexponential fits in the top plot, follow the same color scheme. stellar mass is accreted.
TABLE 3
Halo stellar mass fitting parameters
Run M ∗ α ∗ AMR-FID 9 . × M (cid:12) . × M (cid:12) . × M (cid:12) . × M (cid:12) The average SFR is provided in Figure 19 as well, takenas the change in halo stellar mass divided by the changein time between snapshots. We include the parametricfits, equivalent to the derivative of the exponential fits,given by dM ∗ dt = α ∗ M ∗ H (1 + z ) (cid:112) Ω M (1 + z ) + Ω Λ , = 14 . (cid:12) yr − (cid:16) α ∗ (cid:17) (cid:18) M ∗ M (cid:12) (cid:19) (cid:18) h . (cid:19) × (1 + z ) (cid:112) Ω M (1 + z ) + Ω Λ , (5)where α ∗ is the exponential slope parameter, and M ∗ isthe stellar mass at a particle redshift, which follows anexponential growth fit equivalent to Equation (4), whoseparameters are given in Table 3. The star formation rateis consistently lower in SPH, but grows faster than inAMR. Since the rate is lower, there is more gas presentat later times from which to form stars in SPH. For the FID runs, this is consistent with a longer cooling time inSPH. In the
QSO runs, the SPH feedback has a largerimpact on large scales than in AMR, thus the SFR islower at high redshifts. However, at late times the halo
Fig. 20.—
Top: Evolution of the sum of the mass of the halo’sthree most massive black holes. The
AMR-QSO ( SPH-QSO ) simu-lation is presented with a dashed green (red) line, where we markthe measured values in an output by squares (circles). Middle:Evolution of the total growth rate of the three most massive haloblack holes, taken as the change in mass between adjacent outputs,and thus includes mass growth due to merging. Lines follow thesame color scheme as above. The dark green long-dashed (darkred dotted) line indicates the sum of the Eddington accretion ratesfor these 3
AMR-QSO ( SPH-QSO ) black holes. Bottom: Same asmiddle, but normalized by the Eddington accretion rate sum. gas becomes consistent between the two methods, whichresults in the SFRs approaching better agreement.
Cluster black hole growth
In Figure 20, we show the evolution of the sum of themass of the halo’s three most massive black holes, andtheir net growth rate, both in absolute terms and normal-ized by the sum of their Eddington accretion rate, eachgiven in Equation (3). Since the star formation rate islower in the SPH runs, it takes longer before the stellardensity reaches the threshold to trigger the formation ofa BH. Thus, the SPH BHs form later and are less massiveat very high redshifts. Since BHs are generated at thelimit of density resolution, this could be a large sourceof discrepancy between the two codes. It is thus very in-teresting that the halo environment is so similar at latetimes. Due to the delayed formation, the SPH BHs to-tal growth rate is near Eddington, until the BH massesapproaches that of their AMR counterparts, at whichpoint the growth rates are in good agreement, consistentwithin a factor of two in absolute terms, or roughly thesame value relative to their Eddington rate. This is con-sistent with the AGNs entering the self regulation phase.At z < . z = 3 in Figure 15. During the period ofbest agreement between the two methods, 3 . < z < CONCLUSIONS
To better understand the impact of AGN feedback onthe formation and evolution of a large cluster and howthis evolution can be biased by numerical effects, we havesimulated the growth of a cluster from identical initialconditions in two different numerical methods. Using theAMR code
RAMSES and the SPH code
HYDRA we have at-tempted to match their radiative cooling, star formation,stellar feedback, and their black hole formation, growth,and energetic feedback processes, with the caveat thatthe fundamental differences between the methods willmake modeling even simple structure formation not nec-essarily identical. By comparing these simulations withsuccessively more complex baryonic physics, we have ob-served the following key points: • Regardless of the treatment of baryonic processes,SPH consistently has a lower central entropy pro-file than AMR, with the sole exception being di-rectly after an energetic feedback event. Whilethis has been seen before for non-radiative simu-lations, comparisons with moving mesh codes havesuggested that SPH codes may have higher centralentropy profiles. • Gas that is heated by the virial shock can efficientlycool at high redshift in simulations, although thiscooling may not be captured numerically if an ar-tificial viscosity is employed or if the shock widthis not well resolved. Future work studying the res-olution dependence of this feature will shed furtherlight on this issue. • AGN feedback reduces the baryonic fraction of ha-los by roughly 30%, regardless of the numericalmethod. The baryon fraction is highly dependent on the location of the virial shock, as it sets thebottleneck for subsequent baryonic accretion. • AGN feedback leads to better agreement betweenthe two methods on the thermal state of the halogas, consistent with a picture of self-regulation. Re-gardless of the numerical method, the AGN willaccrete matter until moving most of the gas to en-tropies of roughly 100 keV cm , at which point sub-sequent accretion diminishes. • We see hints of AGN feedback impacting the ter-mination point and orientation of filamentary coldflows in AMR, acting to push these streams beyondthe impact radius of the AGN. In SPH, the fila-ments are more discrete, allowing the AGN heatedgas to escape around the flows, with little impacton their location or orientation. This may also ex-plain the decrease in angular momentum of the cen-tral halo in AMR. • AGN feedback leads to a reduction in the total starformation rate of the cluster halo, by up to an orderof magnitude by z = 3. Future work comparing thetwo methods at lower redshift is needed.In summary, AGN clearly play an important role inthe evolution and regulation of cluster growth. Theirpossible observational impact is becoming clearer as sur-veys become larger, and hydrodynamic simulations be-come more complex. Further work exploring the detailedphysical implications of AGN feedback and its interactionwith the environments is essential for understanding thecosmic history of the Universe.This research has been supported by the Balzan Foun-dation via the University of Oxford. M.L.A.R. was sup-ported by NSF grant AST11-03608 and the NationalScience and Engineering Research Council of Canada.E.S. was also supported by the National Science Founda-tion under grant AST11-03608 and NASA theory grantsNNX09AD106 and NNX15AK826. R.J.T. is supportedby a Discovery Grant from NSERC, the Canada Foun-dation for Innovation, the Nova Scotia Research and In-novation Trust, and the Canada Research Chairs Pro-gram. Simulations were conducted with the CFI-NSRITfunded St Mary’s Computational Astrophysics Labora-tory, and the High Performance Computing environmentat the University of Leicester. M.L.A.R would also liketo thank Debora Sijacki for very insightful discussion,and the International Balzan Prize Foundation for con-tributing to the breadth of this work. REFERENCESAgertz, O., Moore, B., Stadel, J., et al. 2007, MNRAS, 380, 963Aubert, D., Pichon, C., & Colombi, S. 2004, MNRAS, 352, 376Barai, P., Viel, M., Borgani, S., et al. 2013, MNRAS, 430, 3213Bauer, F. E., Fabian, A. C., Sanders, J. S., Allen, S. W., &Johnstone, R. M. 2005, MNRAS, 359, 1481Bauer, A., & Springel, V. 2012, MNRAS, 423, 2558Benson, A. J., Bower, R. G., Frenk, C. S., et al. 2003, ApJ, 599,38Bett, P., Eke, V., Frenk, C. S., et al. 2007, MNRAS, 376, 215Bauer A., & Springel V. 2012, MNRAS, 423, 2558Birdsall, C. K., & Fuss, D. 1997, JCoPh, 135, 141 Birnboim, Y., & Dekel, A. 2003, MNRAS, 345, 349Blondin, J. M., Wright, E. B., Borkowski, K. J., & Reynolds, S.P. 1998, ApJ, 500, 342Bondi, H. 1952, MNRAS, 112, 195Booth, C. M., & Schaye, J. 2009, MNRAS, 398, 53Brinchmann, J. & Ellis, R. 2000, ApJ, 536, L77Bullock, J. S., Dekel, A., Kolatt, T. S., et al. 2001, ApJ, 555, 240Chapon, D., Mayer, L., & Teyssier, R. 2013, MNRAS, 429, 3114Cole, S., Norberg, P., Baugh, C., et al. 2001, MNRAS, 326, 255Couchman H. M. P. 1991, ApJ, 368, L232