Moving mesh cosmology: the hydrodynamics of galaxy formation
Debora Sijacki, Mark Vogelsberger, Dusan Keres, Volker Springel, Lars Hernquist
MMon. Not. R. Astron. Soc. , 000–000 (0000) Printed 25 September 2018 (MN L A TEX style file v2.2)
Moving mesh cosmology: the hydrodynamics of galaxyformation
Debora Sijacki (cid:63) , Mark Vogelsberger , Duˇsan Kereˇs , † , Volker Springel , , and LarsHernquist Harvard-Smithsonian Center for Astrophysics, 60 Garden Street, Cambridge, MA, 02138, USA Department of Astronomy and Theoretical Astrophysics Center, University of California, Berkeley, CA 94720-3411, USA Department of Physics, Center for Astrophysics and Space Sciences, University of California at San Diego,9500 Gilman Drive, La Jolla, CA 92093, USA Heidelberg Institute for Theoretical Studies, Schloss-Wolfsbrunnenweg 35, 69118 Heidelberg, Germany Zentrum f¨ur Astronomie der Universit¨at Heidelberg, ARI, M¨onchhofstr. 12-14, 69120 Heidelberg, Germany
25 September 2018
ABSTRACT
We present a detailed comparison between the well-known smoothed particle hydro-dynamics (SPH) code
GADGET and the new moving-mesh code
AREPO on a numberof hydrodynamical test problems. Through a variety of numerical experiments withincreasing complexity we establish a clear link between simple test problems withknown analytic solutions and systematic numerical effects seen in cosmological sim-ulations of galaxy formation. Our tests demonstrate deficiencies of the SPH methodin several sectors. These accuracy problems not only manifest themselves in idealizedhydrodynamical tests, but also propagate to more realistic simulation setups of galaxyformation, ultimately affecting local and global gas properties in the full cosmologi-cal framework, as highlighted in companion papers by Vogelsberger et al. (2011) andKeres et al. (2011). We find that an inadequate treatment of fluid instabilities in
GAD-GET suppresses entropy generation by mixing, underestimates vorticity generation incurved shocks and prevents efficient gas stripping from infalling substructures. More-over, in idealized tests of inside-out disk formation, the convergence rate of gas disksizes is much slower in
GADGET due to spurious angular momentum transport. Insimulations where we follow the interaction between a forming central disk and orbit-ing substructures in a massive halo, the final disk morphology is strikingly differentin the two codes. In
AREPO , gas from infalling substructures is readily depleted andincorporated into the host halo atmosphere, facilitating the formation of an extendedcentral disk. Conversely, gaseous sub-clumps are more coherent in
GADGET simula-tions, morphologically transforming the central disk as they impact it. The numericalartifacts of the SPH solver are particularly severe for poorly resolved flows, and thusinevitably affect cosmological simulations due to their inherently hierarchical nature.Taken together, our numerical experiments clearly demonstrate that
AREPO deliversa physically more reliable solution.
Key words: methods: numerical – cosmology: theory – cosmology: galaxy formation
Numerical simulations have become an indispensable tool forstudying astrophysical phenomena. Over the last decade, thenumerical accuracy and fidelity of simulation methods have (cid:63)
Hubble Fellow. E-mail: [email protected] † Hubble Fellow. undergone drastic improvements, both in terms of resolvabledynamic range and of the complexity of the physical pro-cesses that are now routinely incorporated into the codes.This remarkable progress in numerical techniques coupledwith the ever increasing power of high performance com-puting platforms has led to major advances in a number oftopics, such as star formation (Klessen et al. 2009), accre-tion disk dynamics and associated jet phenomena (Gammieet al. 2003; De Villiers & Hawley 2003), supernova explo- c (cid:13) a r X i v : . [ a s t r o - ph . C O ] J un Sijacki et al. sions (Janka et al. 2007), black hole coalescence (Pretorius2007), and cosmic structure formation (Springel et al. 2005).Clearly, an in-depth understanding of many astrophys-ical processes, particularly those that are inherently non-linear, relies crucially on the accuracy and realism of thenumerical approach. While the latter is largely determinedby the degree of adequateness, comprehensiveness and con-sistency of the physical model assumed, the former dependssensitively on the discretization scheme adopted for theequations, which with increasing resolution should convergeto the correct solution. Nonetheless, even for some elemen-tary physical problems, with known analytic solutions, sim-ulation methods can sometimes produce poorly convergedresults, or even converge to a solution which is, however,different from the expected result (Springel 2010b).In the context of hydrodynamical cosmological simula-tions, an important example is given by the seminal work ofFrenk et al. (1999) (the Santa Barbara Comparison Project).This study performed a detailed comparison between twelvedifferent simulation codes that tracked the non-radiativeevolution of a forming galaxy cluster in a cold dark mat-ter (CDM) cosmology. The initial conditions were gener-ated independently by each group from a provided lineartheory density or displacement field. Very different numeri-cal methods, which can be broadly classified into smoothedparticle hydrodynamics (SPH) and mesh-based techniques,employing also different gravity solvers and different effec-tive resolutions, were then used by the groups to follow theformation and evolution of the target object. Reassuringly,the Santa Barbara Comparison Project has shown that theglobal properties of the simulated object, both in terms ofdark matter and gas distribution, were in reasonable agree-ment among the codes. However, detailed gas properties ofthe Santa Barbara cluster, notably in the central region,exhibited a much poorer level of consistency. In particular,Frenk et al. (1999) pointed out that there is a systematic dis-crepancy in the central entropy profiles predicted by SPHand mesh-based codes, with the former producing power-law entropy profiles all the way to the centre and the latteryielding some form of entropy core.Even though progress in numerical techniques has led tosubstantial improvements in current simulation codes com-pared to those considered in Frenk et al. (1999), the sys-tematic difference in the predicted central cluster entropystill persists today, and its origin has not been fully under-stood thus far. Agertz et al. (2007) have performed a setof numerical experiments with SPH and mesh codes, wherethey have simulated the evolution of a cold, dense blob in ahot windtunnel. By comparing the outcomes from differentmethods (and at various resolutions) against the character-istic blob disruption timescale expected analytically, theyconclude that in SPH codes the development of fluid in-stabilities can be numerically hampered. The suppressionof dynamical fluid instabilities, such as Kelvin-Helmholtz,Rayleigh-Taylor, and Richtmyer-Meshkov instabilities, leadsto less efficient mixing of fluid elements with different specificentropy, and hence also inhibits entropy generation throughmixing in the simulated system. In simulations of collidingisolated galaxy clusters, Mitchell et al. (2009) have shownthat the central entropy profiles obtained with the SPH code
GADGET and the mesh code
FLASH show a similar level of discrepancy as found by Frenk et al. (1999), which they at-tribute to the different levels of mixing.There have been a number of attempts to improve thedescription of mixing in SPH by modifying the standardset of discretized equations (see e.g. Price 2008; Wadsleyet al. 2008; Heß & Springel 2010) or by incorporating a Rie-mann solver in place of the artificial viscosity (e.g. Inutsuka2002; Cha & Whitworth 2003; Murante et al. 2011). Forexample, Price (2008) has shown that the introduction ofan artificial thermal conductivity term can improve the be-havior of fluid elements at contact discontinuities, which inturn leads to better developed Kelvin-Helmholtz instabili-ties. More generally, Wadsley et al. (2008) suggested that aphysical modeling of heat diffusion is necessary when simu-lating high Reynolds number flows (both for SPH and mesh-based methods), and that by applying such an approach aflat entropy core is likely a more correct solution for non-radiative galaxy cluster simulations.Differences in the hydrodynamical solver between SPHand mesh-based codes are possibly not the only reason forthe systematically different central entropy profiles in theSanta Barbara Comparison Project. As mentioned above,the groups involved in this study did not use identical initialconditions, and did not perform their simulations at equalnumerical resolutions and with identical gravity solvers,which opens up the possibility of additional sources for dis-crepancies. A more recent code comparison study by Heit-mann et al. (2008) evolved uniform, dark matter only cos-mological boxes in a ΛCDM universe with 10 different codes,starting from the same initial conditions. They showed thatfor large systems there is a reassuring agreement in the halomass functions between the codes, as well as in the internalstructures of haloes in the outer regions. However, the studyby Heitmann et al. (2008) also revealed significant discrepan-cies for small halos, demonstrating that the typical root gridresolution commonly adopted in adaptive mesh refinement(AMR) codes for simulations of cosmic structure formationis overly coarse and leads to a suppression of low mass haloformation (see also O’Shea et al. 2005). This emphasizes theneed to simultaneously strive for high accuracy both in thegravity solver and in the hydrodynamics. Additionally to thehigh accuracy of the gravity and hydro solver, it is also ofprime importance to adopt sufficiently high mass and spatialresolution for the studied problem at hand. For example, ifbesides pure hydrodynamics other physical processes, suchas star formation and associated feedback, are modelled incosmological simulations, it is necessary to resolve all star-forming haloes with a sufficiently large number of resolutionelements to reach convergent results (Springel & Hernquist2003b).In the present work we perform a detailed comparisonstudy between the widely used SPH code
GADGET and thenovel moving-mesh code
AREPO on a number of hydrody-namical test problems with increasing levels of complexity.We have adopted a combination of existing test problemsand newly devised numerical experiments (cid:63) in order to pro-vide a clear link between the results of our test problems and (cid:63) (cid:13)000
AREPO on a number of hydrody-namical test problems with increasing levels of complexity.We have adopted a combination of existing test problemsand newly devised numerical experiments (cid:63) in order to pro-vide a clear link between the results of our test problems and (cid:63) (cid:13)000 , 000–000 he hydrodynamics of galaxy formation those of full cosmological simulations, which are discussedin detail in our companion papers (Vogelsberger et al. 2011,hereafter Paper I) and (Keres et al. 2011, hereafter Paper II).Our work thus contributes to the understanding of the signif-icant differences in baryon properties found in cosmologicalsimulations between mesh-based and SPH techniques, bothat a global level (see Paper I) and at the scale of individualgalaxies (see Paper II).A unique advantage of our comparison study lies in thefact that simulations with GADGET and
AREPO can bestarted from identical initial conditions and that both codesemploy the same gravity solver. Note that in case of
AREPO gravitational softenings for the gas can be kept fixed as is thecase of
GADGET or can be computed adaptively determinedby the cell size. In the numerical experiments with gas self-gravity we explore both fixed and adaptive gas softening,where in the case of adaptive softenings a floor equal to the
GADGET gas softenings is set. We find that these differentchoices of gravitational softenings in
AREPO do not affectour findings. This allow us to isolate cleanly how the dif-ferent hydro solvers used by
GADGET and
AREPO impactgas properties. We first consider elementary hydrodynami-cal numerical tests such as a strong 1D Sod shock tube test,a 2D implosion test (which is widely used for benchmark-ing mesh codes, but has rarely been considered in SPH),and the “blob” test (Agertz et al. 2007). We then perform anumber of isolated or merging halo simulations in the non-radiative regime, aimed at understanding how shocks andfluid instabilities affect their gaseous atmospheres. Finally,we study the differences between
GADGET and
AREPO inradiative simulations, where we follow inside-out disk forma-tion and interactions between the central disks and orbitingsubstructures.This paper is organized as follows. In Section 2, we pro-vide a brief overview of the numerical codes used and thetypes of simulations performed. Section 3 represents the coreof the paper where all of our numerical experiments are dis-cussed. Finally, we summarize our findings in Section 4.
GADGET
GADGET is a massively parallel TreePM-SPH code widelyused in numerical astrophysics. In this study we adopt thelatest
GADGET-3 version. A detailed description of an ear-lier version can be found in Springel (2005).
GADGET is fullyadaptive in time and space, and in its entropy formulationfor SPH (Springel & Hernquist 2002) manifestly conservesboth energy and entropy in the absence of artificial viscosity.Gravitational forces are computed with an octtree method(Barnes & Hut 1986; Hernquist 1987). To speed up the com-putation, long-range forces can be optionally evaluated witha PM method, with the tree being restricted to short-rangeforces only.In all our tests, we adopt as standard value for the artifi-cial viscosity strength α = 1 .
0, and we use 64 neighbours forkernel interpolation in three dimensional simulation, unlesswe specifically vary these parameters to assess their effect.
As mentioned in the introduction, there have been a num-ber of recent proposals to improve the standard SPH imple-mentation in various ways, for example by invoking a time-dependent artificial viscosity (Morris & Monaghan 1997;Dolag et al. 2005), a modified density estimate (e.g. Ritchie& Thomas 2001; Heß & Springel 2010; Saitoh & Makino2012), a decoupling of the hot and cold neighbours in mul-tiphase flows (Marri & White 2003; Okamoto et al. 2003),an artificial thermal conductivity (Price 2008), a modifiedequation of motion (Heß & Springel 2010; Abel 2011), an ex-plicit modelling of mixing (Wadsley et al. 2008), an enlargedneighbour number combined with a different kernel shape(Read et al. 2010), or a replacement of the artificial viscosityby a Riemann solver (e.g. Inutsuka 2002; Cha & Whitworth2003; Murante et al. 2011). While some of these modifica-tions of the standard SPH formalism deliver more accurateresults in targeted numerical experiments, no consensus hasyet emerged whether any of these approaches (or a combina-tion thereof) is sufficiently robust for cosmological applica-tions and leads to universally more accurate results in galaxyformation simulations. Therefore, in this work we focus onthe traditional SPH implementation rather than on the var-ious possible modifications proposed recently. We note thatthis standard formulation of SPH has also been typicallyemployed in state-of-the art cosmological SPH calculations(e.g. Crain et al. 2009; Di Matteo et al. 2012). Furthermore,as discussed in Paper I, there are other, generic issues withSPH that have not been resolved by any of the above mod-ifications. These issues ultimately mean that SPH does notcurrently have a formal convergence condition, which alsocomplicates rigorous evaluations of variants of the standardSPH algorithm.
AREPO
AREPO (Springel 2010a) is a new massively parallel sim-ulation code, which uses the same gravity solver as
GAD-GET (augmented with the possibility of adaptive gravita-tional softenings for the gas), but employs a completelydifferent method for the evolution of the fluid. It adoptsa second-order accurate finite volume technique, where thesolution of the Euler equations is computed by an unsplitGodunov method equipped with an exact Riemann solver.Throughout we use the default choice of the slope limiterin
AREPO which prevents the linear reconstruction to over-or undershoot the maximum/minimum values of neighbour-ing cells, as described in detail in Springel (2010a). Unlikestandard finite-volume codes used in numerical astrophysics,
AREPO solves the equations on an unstructured Voronoimesh, which is allowed to freely move with the fluid. Theresulting quasi-Lagrangian nature of
AREPO automaticallyguarantees spatial adaptivity and greatly reduces numericaldiffusivity even in the presence of large bulk flows. Com-pared to standard Eulerian mesh codes,
AREPO has the ad-vantage of being fully Galilean invariant (as is
GADGET aswell), it is less prone to advection errors and over-mixing,and preserves contact discontinuities better. c (cid:13) , 000–000 Sijacki et al.
In this study, we employ a mesh regularization method † in AREPO by default, based on a Lloyd algorithm (for de-tails see Springel 2010a), which ensures that the geometriccentre of each Voronoi cell is sufficiently close to the cell’smesh-generating point to ensure good accuracy of the spa-tial reconstruction. In cosmological simulations (see PaperI), an alternative regularization criterion has proven to beadvantageous, based on the maximum opening angle underwhich a cell face is seen from the mesh-generating point.We have checked for a number of test runs presented in thisstudy that this alternative regularization method does notaffect any of the results described here.For most of the simulations presented here we do notuse the possibility of mesh refinement and de-refinement op-erations (see Springel 2010a), except in our numerical ex-periments with star formation presented in Section 3.4.3,where we employ it for verification of our findings. Meshde-/refinements are used to constrain the mass of cells to asmall range around a target value (equal to the gas parti-cle mass in the matching
GADGET run). This restricts thespectrum of star particle masses which are generated fromgaseous cells, and thus assures that N-body heating effectsare minimized. Also, as our default choice we use the en-ergy formulation of
AREPO . We verified for each numericalexperiment that there is no significant spurious transfer ofkinetic into thermal energy.
We perform both radiative and non-radiative hydrodynam-ical simulations, where in the former case gas is representedby a primordial mixture of hydrogen and helium in an op-tically thin limit (Katz et al. 1996). In simulations withradiative cooling, we employ a subresolution multi-phasemodel for star formation and associated supernova feedback(Springel & Hernquist 2003a). We slightly modify the be-haviour of this model for gas elements which are hot butalready above the density threshold for star formation, byallowing them to settle quickly onto the effective equationof state: if their newly estimated temperature from radiativecooling would fall below the temperature of the multiphasemedium, we set it equal to the multiphase temperature. Thischange has been introduced to make the sub-grid star for-mation module consistent with its current implementation in
AREPO . We also perform some simulations with the subreso-lution star formation model where spawning of star particlesis intentionally prevented, but the cold, dense gas above thedensity threshold for star formation is still governed by theeffective equation of state. These simulations are particu-larly useful for understanding the development of dynami-cal instabilities between cold and hot media, and have theadvantage of not being prone to fragmentation which mightaffect pure cooling runs.Note that we deliberately consider only very simplifiedbaryonic physics implementations in all presented numeri-cal experiments, as detailed above. This is for two reasons. † Note that in the presence of mesh regularization the Galileaninvariant nature of the AREPO code is not violated.
First, we need to make sure that gas cooling and star forma-tion processes are treated on an as equal footing as possiblein order to allow a meaningful comparison of the two codes.Thus, we adopt a simple multi-phase model for star forma-tion which is largely insensitive to the detailed structureof the gas on very small scales where gas joins the inter-stellar medium and is described by a relatively stiff effec-tive equation of state. Second, our aim is to isolate in anas clean manner as possible the differences between simu-lated systems with
GADGET and
AREPO stemming fromthe discretization of the fluid equations alone. This goalis given precedence over trying to reproduce observationalfindings through a more sophisticated modelling of addi-tional physics. While it is likely that such additional physi-cal processes will change the properties of some of our simu-lated systems by possibly different degrees in
GADGET and
AREPO , it is of significant interest in its own right to dis-entangle numerical uncertainties from uncertainties in thephysical modelling of star formation and associated feed-back processes.
For simplicity, many of the numerical experiments presentedin this study have been evolved without gas self-gravity. Forsimulations where self-gravity of the gas is nevertheless in-cluded, we note that
GADGET employs constant gravita-tional softening for gas particles in the manner of Hernquist& Katz (1989), while
AREPO uses either the same constantgravitational softening or an adaptive softening determinedby the cell size with a minimum softening value set to thefixed softening of the matching
GADGET run. We have veri-fied that this does not lead to substantial differences for anyof the tests presented in this study.
Even though
GADGET and
AREPO calculations use thesame gravity solver and can be initiated from identical ini-tial conditions, due to the very different nature of the hydrosolver it is not straightforward to define unambiguously acomparison strategy at the same or equivalent hydro reso-lution. In this study we choose to keep the number of res-olution elements the same in the both codes, i.e. to havethe same number of SPH particles as Voronoi cells, corre-sponding to a comparable mass resolution in the gas. Thisallows us to adopt indeed the same initial conditions and toalso have similar mass resolution in the stellar componentin those simulations where star formation is included. Fur-thermore, the CPU costs of the two codes are then roughlycomparable, as discussed in Paper I.We note, however, that this choice implies that the “ef-fective” spatial resolution of
GADGET is lower than thatof the matching moving-mesh calculation, given that fluidproperties in SPH are evaluated by kernel averaging overa somewhat larger number of neighbours than needed in
AREPO for its stencil of neighbouring cells, for example forgradient estimates. For this reason, we perform all of ournumerical experiments at a number of different resolutions,which also help us to gain some insight into the convergenceproperties of the two codes. Nonetheless, it is important to c (cid:13)000
AREPO for its stencil of neighbouring cells, for example forgradient estimates. For this reason, we perform all of ournumerical experiments at a number of different resolutions,which also help us to gain some insight into the convergenceproperties of the two codes. Nonetheless, it is important to c (cid:13)000 , 000–000 he hydrodynamics of galaxy formation stress that the convergence properties of the SPH methodare still not well understood. For example, one would ulti-mately require that the number of neighbours be increasedwith increasing total particle number (see e.g. Rasio 2000;Read et al. 2010; Robinson & Monaghan 2011, and discus-sion in Paper I), but the appropriate scaling of the neighbournumber with increasing resolution is currently unknown. Itis common practice, which we adopt as well, to simply al-ways keep the number of neighbours fixed when the totalparticle number is varied, even though the discretized rep-resentation of the density field is not guaranteed to convergeto its underlying smooth distribution with this choice. For all the tests presented in this study (except for the“blob” experiment of Agertz et al. 2007, for which we takepublicly available initial conditions) we generate the initialconditions in terms of simple particle/point distributions.Depending on the simulated problem at hand, we either pop-ulate the simulated domain with particles uniformly spacedon a regular lattice or we adopt Poisson sampling of a den-sity field. For simulations with
AREPO , the positions of par-ticles in the initial conditions define mesh generating points.All purely hydrodynamical numerical experiments have beenperformed with periodic initial conditions. In the case of sim-ulated problems where we follow the evolution of an isolatedgaseous halo or merging haloes, we select vacuum boundaryconditions for
GADGET . In simulations with
AREPO , it ishowever necessary to enclose the simulation in a finite vol-ume in order to make the tessellation with a mesh well posed.If needed, we thus add a low resolution background grid toour initial conditions, which is chosen sufficiently large asto encompass the spatial extent of our simulated objects atall times. The procedure adopted for this background gridgeneration, which also reduces the Poisson noise in the ini-tial mesh geometry, is described in detail in Section 9.4 ofSpringel (2010a).
As an introductory problem, we consider a strong shock witha Mach number M = 6 . L x = 10 . N = 200 particles (cells) havebeen placed on a regular grid such that for x < . P = 30 . ρ = 1 .
0, while for x ≥ . P = 0 .
14 and ρ = 0 . γ = 1 .
4, and assume that the fluid isinitially at rest. In the test run with
GADGET , we adopt astandard value of the artificial viscosity equal to α = 1 . N ngb by setting it to 5, 7,11, or 15, appropriate for the 1D nature of the test.In Figure 1, we show the gas density, velocity, en-tropy (i.e. P/ρ γ ) and pressure at time 0 .
25 for
GADGET with N ngb = 5 (left-hand panels), and for AREPO (middle panels). Blue symbols give the values of individual parti-cles/cells, dashed red lines represent the initial conditions,while dotted red lines are the analytic solution to this Rie-mann problem. It can be seen that in both
GADGET and
AREPO the post-shock properties of the fluid are capturedwell, but the shock and the contact discontinuity are signifi-cantly broader in
GADGET which also shows a characteristic’pressure-blip’ at the contact discontinuity (see also Springel2010b). Contrary to what one may perhaps suppose, thepost-shock oscillations present in
GADGET are not causedby inadequate artificial viscosity, but are instead induced byan inaccurate treatment of the sharp contact discontinuityof the initial conditions. To demonstrate this point, we haveperformed exactly the same shock tube test but this timesmoothing the initial contact discontinuity over 5 particlesin density and internal energy with a Hann window func-tion, so as to reduce its sharpness. Green cross symbols inthe third row of Figure 1 (see also the right-hand panelswhere we zoom into the region around the shock) illustratehow the gas entropy is affected by this choice of smoothedinitial conditions. Post-shock oscillations and the “spike” inentropy at the contact discontinuity are greatly reduced.It is interesting to note that also in the case of
AREPO the “spike” in the entropy at the contact discontinuity iscured by our smoothed initial conditions. This feature isabsent in the results of standard grid codes, as we have veri-fied by running this test with a static mesh option in
AREPO (and by running it with
ENZO ), which tends to broaden con-tact discontinuities (depending on the advection speed) andthus largely wash-out this feature. This can be seen in thelower right-hand panel of Figure 1 where we show the resultsfrom the
AREPO run with a static mesh with magenta trian-gle symbols. The much lower numerical diffusivity of
AREPO for contact discontinuities preserves the initial start-up fea-ture to much higher accuracy, but at the same time this canlead to a larger ‘wall heating’ effect (Rider 2000) than instatic grid codes, which tend to smooth out at some levelthe initial start-up errors at the contact discontinuity (seealso the description of the Noh problem in Springel 2010a).When we increase the number of neighbors from N ngb =5 (as shown in Figure 1) to 15 in the shock tube tests per-formed with GADGET (but keeping the total number of par-ticles constant), the spatial region over which the contactdiscontinuity and the shock are broadened increases pro-gressively, but the pressure jump between the post-shockand pre-shock gas outside of the broadened region remainthe same, yielding consistent Mach numbers. Based on thesetests we later discuss in more detail the consequences ofshock broadening in
GADGET when simulating the radialinfall of gas into static dark matter haloes in Section 3.3.2.
While one-dimensional shock tube tests are an essential ba-sic benchmark for hydrodynamical code performance, theyare far less demanding than multi-dimensional flow wherecomplex interactions of non-planar shocks may occur, as isthe case in realistic structure formation simulations. Suchtests have however not been examined widely in the SPH lit-erature thus far. We hence perform the so-called “implosion c (cid:13) , 000–000 Sijacki et al. x x v x x P / (cid:0) x P GADGET x x v x x P / (cid:0) x P AREPO x P / (cid:0) x P / (cid:0) GADGETAREPO
Figure 1.
One dimensional shock tube problem with Mach number M = 6 . t = 0 .
25. Left-hand panels: GADGET with standardartificial viscosity α = 1 . N ngb = 5. Middle panels: AREPO with moving mesh. Dashed red lines denote initial conditions, whiledotted red lines represent the analytic solution. In the third row, green cross symbols illustrate gas entropy in a shock tube of the samestrength but with initially smoothed contact discontinuity adopting a Hann window. Right-hand panels: zoom into the entropy profilearound the shock - here magenta triangle symbols in the lower panel are for the AREPO run with a static mesh. test” in two dimensions (Hui et al. 1999) ‡ . To set up this testproblem, we select a computational domain L x = L y = 0 . N x = N y = 200 particles or cells, respectively (we havealso run higher resolution versions with N x = N y = 400and 800). The initial pressure and density are P = 1 . ρ = 1 . x + y > .
15, while P = 0 .
14 and ρ = 0 .
125 other-wise. The adiabatic index is γ = 1 . GADGET adopting anartificial viscosity of α = 1 . N ngb = 22 smoothingneighbours (suitable for this 2D test). While previously thistest has been considered using reflective boundary condi-tions, we here adopt periodic boundary conditions due totheir more straightforward implementation in SPH codes.In Figure 2, we show density maps for simulations per-formed with GADGET (left-hand panels) and with
AREPO (right-hand panels), at times t = 0 .
1, 0 .
3, 0 . . N x = N y = 200. The complex, evolving gas density ‡ See also “Comparison of Several difference schemeson 1D and 2D Test problems for the Euler equa-tions” by Liska, R. and Wendroff, B., on ∼ liska/CompareEuler/compare8/ and ∼ jstone/Athena/tests/ . structure is caused by the continuous interaction of shocksthroughout the computational domain. Initially, due to thediscontinuity in density and pressure along x + y = 0 .
15, aplanar mild shock front develops perpendicular to the x = y diagonal traveling towards the origin. Given that we haveadopted periodic boundary conditions, the gas will interacton all four sides of the simulated box, and in particular, in-teracting shock waves in the lower left-hand corner resultin the formation of a narrow jet along the x = y diago-nal (which is clearly visible only in the AREPO run). As thetraveling shocks accelerate the fluid in the regions of densitydiscontinuity a Richtmyer-Meshkov instability develops, asmanifested by “mushroom-like” features in the gas distribu-tion.While the global gas density distribution in Figure 2agrees reasonably well in
GADGET and
AREPO runs (i.e.with respect to the shape of the regions with different den-sities and the magnitude of the density differences), there areseveral significant differences: i) shocks and contact discon-tinuities are much sharper in AREPO , in agreement with ourfindings in Section 3.1.1; ii) the density distribution in
GAD-GET appears not only smoothed out, but it is also noisier,as evidenced by the graininess of the density maps, which iscaused by intrinsic noise in multidimensional flows in SPH c (cid:13)000
GAD-GET appears not only smoothed out, but it is also noisier,as evidenced by the graininess of the density maps, which iscaused by intrinsic noise in multidimensional flows in SPH c (cid:13)000 , 000–000 he hydrodynamics of galaxy formation GADGET y ρ g a s AREPO y ρ g a s GADGET y ρ g a s AREPO y ρ g a s GADGET y ρ g a s AREPO y ρ g a s GADGET x y ρ g a s AREPO x y ρ g a s Figure 2.
Implosion test in 2D at times t = 0 .
1, 0 .
3, 0 . .
7. Left-hand panels: GADGET ( α = 1 . N ngb = 22). Right-hand panels:AREPO moving-mesh. For each row, the density scale is the same in the left-hand and right-hand panels, covering the following densityrange: ρ gas = 1 . − . (cid:13) , 000–000 Sijacki et al.
GADGET: 200 x y V o r ti c it y GADGET: 800 x y V o r ti c it y AREPO: 200 (FD) x y V o r ti c it y AREPO: 200 (Internal) x y V o r ti c it y Figure 3.
Vorticity maps of the implosion test in 2D at time t = 0 .
3. Upper panels: GADGET runs at two different resolutions with N x = N y = 200 (left-hand panel) and N x = N y = 800 (right-hand panel). For both resolutions, vorticity maps are computed by finitedifferencing the velocity field. Lower panels: vorticity maps of the AREPO run with N x = N y = 200 computed by finite differencingthe velocity field (left-hand panel), and by calculating vorticity in the code based on a discretized curl-operator directly applied to theVoronoi cells (right-hand panel). (Springel 2010b); iii) most strikingly, perhaps, is the quitedifferent appearance of the low density gas in the bottomleft corner – a narrow, extended, dense jet along the di-agonal is largely absent in GADGET (due to the broaden-ing of the contact discontinuity) and Richtmyer-Meshkovinstabilities along discontinuities are suppressed. We haveverified that higher resolution simulations with
GADGET ( N x = N y = 400 and 800) result in sharper shocks andcontact discontinuities (with a somewhat feeble jet form-ing), but they cannot satisfactorily cure the absence ofwell-developed Richtmyer-Meshkov instabilities. This is inagreement with previous studies (e.g. Agertz et al. 2007;Springel 2010b) indicating fundamental limitations of theSPH method in its widely used form. We also note that themoving mesh in AREPO does an excellent job of maintain-ing symmetry across the x = y diagonal. Furthermore, thelevel of numerical diffusion of discontinuities is very low, asindicated by the narrowness and length of the jet.In addition to the density fields it is also instructiveto analyze the vorticity distribution in the implosion test, as shown in Figure 3. To construct vorticity maps we ei-ther i) first compute spatially adaptive velocity field mapsfor each component on a uniform grid (in our case only x and y components) and then finite difference them to obtainthe curl or ii) we compute spatially adaptive vorticity maps,where the curl has been evaluated directly in the code. For GADGET runs, vorticity maps computed with method i) areshown in the upper panels of Figure 3 for N x = N y = 200and N x = N y = 800. Using method ii) vorticity maps lookessentially the same. In the case of AREPO (shown in thelower panels for the same run with N x = N y = 200) thevorticity maps are also very similar when computed withmethod i) or ii) . However, regardless of the exact details ofthe vorticity map generation there are marked differencesin the vorticity field between GADGET and
AREPO . Forthe N x = N y = 200 run, the vorticity map in GADGET islargely featureless and noisy, with artificial suppression ofvorticity generation at locations where surfaces of constantdensity and of constant pressure are not aligned (i.e. baro-clinic source term ∝ ∇ ρ × ∇ p ). Strikingly, even in the case c (cid:13)000
AREPO . Forthe N x = N y = 200 run, the vorticity map in GADGET islargely featureless and noisy, with artificial suppression ofvorticity generation at locations where surfaces of constantdensity and of constant pressure are not aligned (i.e. baro-clinic source term ∝ ∇ ρ × ∇ p ). Strikingly, even in the case c (cid:13)000 , 000–000 he hydrodynamics of galaxy formation of the high resolution GADGET run with N x = N y = 800(upper right-hand panel) vorticity generation in the regionswith large baroclinic term is largely suppressed compared tothe AREPO simulation, even though the latter has a muchsmaller number of resolution elements. It is also worth not-ing that the overall magnitude of the vorticity is relativelyhigh in
GADGET with respect to
AREPO in the regionswhich should be characterized by having low vorticity. Thisis due to the noise alluded to above, visible both in the den-sity and vorticity maps. This noise is inherently present inmultidimensional SPH simulations. It is caused by inaccu-rate pressure gradient estimates and is particularly evidentin the subsonic regime (for further details see e.g. Springel2010b; Bauer & Springel 2012). The jitter in gas velocitiescaused by the noisy gradient estimates introduces an artifi-cial vorticity ’floor’ throughout the simulated volume of thebox. The poor treatment of vorticity in
GADGET has alsodirect consequences for the effectiveness of the widely usedBalsara switch (Balsara 1995) for artificial viscosity whichrelies on evaluation of velocity divergence and curl.Hence, from the results of the “implosion test” we con-clude that while
GADGET can accurately capture the mainfluid properties in the case of dynamically interacting shockswith complicated geometries, there are systematic biases indetailed aspects, related to the development of fluid insta-bilities and to the level of fluid mixing at different entropies.These unwanted features lead to the suppression of angu-lar momentum transport by vortices and to the damping ofturbulence in the wake of curved shocks (see also Bauer &Springel 2012).
Our findings from the implosion test (Section 3.1.2) indicatethat there are fundamental differences in the treatment offluid instabilities between
GADGET and
AREPO . To investi-gate this issue in detail, we first consider the so-called “blob”test as proposed by Agertz et al. (2007), which has been an-alyzed with many different codes by now (Agertz et al. 2007;Heß & Springel 2010; Murante et al. 2011). The idea behindthis test is to simulate the evolution of a dense cold blobimmersed in a hot windtunnel, mimicking in a simplifiedway the motion of a satellite galaxy through the intraclus-ter medium (see also Heß & Springel 2011). If the relativevelocity between the blob and the surrounding medium is su-personic, a bow shock will develop in front of the blob. Addi-tionally, dynamical instabilities (mostly of Kelvin-Helmholtzand Rayleigh-Taylor type) will grow in the subsonic part ofthe flow between the bow shock and the surface of the blob(see Agertz et al. 2007, for a detailed description). Theseinstabilities will greatly influence the evolution of the blob,leading to its eventual disintegration on the characteristicKelvin-Helmholtz timescale set by the initial conditions.To simulate this problem, we adopt the initial con-ditions used in the original Agertz et al. (2007) paper,which are publicly available § . The simulated domain con-sists of a periodic box with extensions L x = 2000 kpc, § L y = 2000 kpc, L z = 8000 kpc, and the blob is initiallyplaced at x c = y c = z c = 1000 kpc. The blob has aradius R blob = 197 kpc, and it is ten times colder anddenser than the surrounding medium, such that pressureequilibrium is ensured. The density and temperature of theexternal medium is ρ medium = 3 . × − M (cid:12) kpc − and T medium = 10 K, respectively, and it is moving with aconstant velocity v medium = 1000 km s − along the z -axis.The adiabatic index is set to γ = 5 /
3. Following Agertzet al. (2007), we define the characteristic Kelvin-Helmholtztimescale as t KH = 1 . t cr , where t cr is the blob crushingtime defined as t cr = 2 R blob ( ρ blob /ρ medium ) / /v medium . Forthese initial conditions, the characteristic Kelvin-Helmholtztimescale is t KH ∼ .
98 Gyr.We have performed the blob test at three different res-olutions with
GADGET and
AREPO . The resolutions usedare: 32 × × × × × × y = 100 kpc) centred around the blob position.Left-hand panels illustrate GADGET runs at times t = t KH , t = 2 t KH and t = 3 t KH , while the right-hand panels arefor the simulations with AREPO in the moving mesh mode.The position and the shape of the bow shock are rather sim-ilar between
GADGET and
AREPO runs, especially at earlytimes, but the shock is broader and less crisp in
GADGET .In the run with
GADGET , the blob acquires a cap-like ap-pearance caused by the internal shock which compresses it,and it undergoes continuous ablation due to the low pres-sure region which forms in the wake of the blob (Agertzet al. 2007). While initially the blob evolution is similar in
AREPO , after t = t KH well developed Kelvin-Helmholtz andRayleigh-Taylor instabilities lead to an efficient shredding ofthe blob, mixing it with the external medium. In Figure 5we also show projected surface density maps centred aroundthe blob position for the 32 × ×
128 and 64 × × t = 3 t KH .We quantify the time evolution of the blob in Fig-ure 6. Here we show the remaining blob mass fraction asa function of time, where the material associated with theblob is selected such as to satisfy: T < . T medium and ρ > . ρ blob , as has been done in the previous studies.We compare GADGET and moving-mesh
AREPO simula-tions, performed at three different resolutions (as indicatedon the legend). While for our highest resolution runs thereis a broad agreement between
GADGET and
AREPO for t < t KH , blob mass fractions are systematically differentafterwards. This is exactly when the transition in the massloss rate from the ablation-dominated to the fluid instabil-ity dominated regime occurs. In the latter regime, AREPO clearly delivers a physically more trustworthy solution. Themoving-mesh
AREPO simulation agrees qualitatively verywell with the outcome of Eulerian grid codes studied inAgertz et al. (2007). However, as noted by Springel (2011),there is a small but systematic difference with the moving-mesh code delivering slightly higher remaining blob massfor t > . × t KH . As expected, the GADGET simulationresults agree well with the findings of Agertz et al. (2007)for the other SPH codes. This indicates that the inaccura- c (cid:13) , 000–000 Sijacki et al. t = t KH GADGET: 128x128x512 z [ kpc ] x [ kp c ] -6 -5 -4 Σ g a s -6 -5 -4 t = t KH AREPO: 128x128x512 z [ kpc ] x [ kp c ] -6 -5 -4 Σ g a s -6 -5 -4 t = 2 t KH GADGET: 128x128x512 z [ kpc ] x [ kp c ] -6 -5 -4 Σ g a s -6 -5 -4 t = 2 t KH AREPO: 128x128x512 z [ kpc ] x [ kp c ] -6 -5 -4 Σ g a s -6 -5 -4 t = 3 t KH GADGET: 128x128x512 z [ kpc ] x [ kp c ] -6 -5 -4 Σ g a s -6 -5 -4 t = 3 t KH AREPO: 128x128x512 z [ kpc ] x [ kp c ] -6 -5 -4 Σ g a s -6 -5 -4 Figure 4.
Projected surface density maps in units of [M (cid:12) kpc − ] at times t = t KH (top row), t = 2 t KH (middle row) and t = 3 t KH (bottom row) for GADGET (left-hand panels; α = 1 . N ngb = 64) and AREPO (right-hand panels; moving mesh). The thickness ofthe slices is ∆ y = 100 kpc, and they are centred on y c = 1000 kpc. Dynamical fluid instabilities are responsible for blob shredding onthe characteristic Kelvin-Helmholtz timescale in the case of AREPO, but they are artificially suppressed in the run with GADGET,prolonging the survival time of the blob. cies we find for GADGET are inherent to the standard SPHmethod, and have prompted a number of works suggestingpossible improvements to the SPH method (e.g. Price 2008;Wadsley et al. 2008; Heß & Springel 2010; Saitoh & Makino2012, see also Section 2.1.2).
To explore more directly the impact of shocks and fluid in-stabilities in the context of structure formation, we have de-vised a number of idealized test problems involving isolatedhalo models. Here we briefly describe how we set up and ver-ify hydrostatic equilibrium configuration of the initial con-ditions, which form the backbone for a series of numericalexperiments discussed in Sections 3.3 and 3.4.The initial conditions are constructed by populatingstatic or live dark matter potentials with gas particles, whose positions are drawn randomly. For the dark matter distribu-tion we assume a Hernquist profile (Hernquist 1990) so thatwe can easily generate self-consistent models when we uselive haloes. The gas density profile traces the dark matterat large radii but is slightly softened in the centre, i.e. ρ gas ( r ) = M vir (2 πa )( x + x )( x + 1) , (1)where M vir is the virial mass of the system, a is the Hern-quist scale length parameter, x = r/a , and x = 0 .
01 isthe softening scale length parameter for the gaseous halo.For the simulations without any net rotation, the initial gasvelocities are set to zero, while for the simulation with non-vanishing angular momentum we assign gas velocities suchthat the halo is characterized by the dimensionless spin pa-rameter λ = J | E | / GM / , (2) c (cid:13)000
01 isthe softening scale length parameter for the gaseous halo.For the simulations without any net rotation, the initial gasvelocities are set to zero, while for the simulation with non-vanishing angular momentum we assign gas velocities suchthat the halo is characterized by the dimensionless spin pa-rameter λ = J | E | / GM / , (2) c (cid:13)000 , 000–000 he hydrodynamics of galaxy formation t = 3 t KH GADGET: 32x32x128 z [ kpc ] x [ kp c ] -6 -5 -4 Σ g a s -6 -5 -4 t = 3 t KH AREPO: 32x32x128 z [ kpc ] x [ kp c ] -6 -5 -4 Σ g a s -6 -5 -4 t = 3 t KH GADGET: 64x64x256 z [ kpc ] x [ kp c ] -6 -5 -4 Σ g a s -6 -5 -4 t = 3 t KH AREPO: 64x64x256 z [ kpc ] x [ kp c ] -6 -5 -4 Σ g a s -6 -5 -4 Figure 5.
Projected surface density maps in units of [M (cid:12) kpc − ] at t = 3 t KH for the low (top row) and intermediate resolution simulation(bottom row) with GADGET and AREPO. The thickness of the slices is ∆ y = 100 kpc and they are centred on y c = 1000 kpc. Even atour lowest resolution, AREPO captures the dynamical evolution of the cold blob much more accurately. t / t KH M c l ( t ) / M c l ( ) GADGET: 32x32x128GADGET: 64x64x256GADGET: 128x128x512AREPO: 32x32x128AREPO: 64x64x256AREPO: 128x128x512
Figure 6.
The remaining blob mass fraction as a function of timein units of t KH . GADGET and AREPO results at three differentresolutions are shown, as indicated in the legend. where J represents the angular momentum, E is the totalenergy of the halo, and we assume solid body rotation.For validation purposes we evolve isolated haloes with GADGET and
AREPO for 2 .
45 Gyr with gas self-gravity andno radiative losses. These test runs confirm that the gas isin very good hydrostatic equilibrium within the dark mat-ter potential. The differences in gas density, temperature and entropy between
GADGET and
AREPO are within a few per-cent throughout the whole halo at the final time in the caseof static dark matter haloes (see Figure A1 in Appendix A;for live haloes see Section 3.5).We also observe that even though the initial gas ve-locities are zero (in the case with λ = 0) some small gasvelocities develop over time ( σ gas , ∼
30 km s − ). This isprimarily caused by the initial Poisson sampling of gas po-sitions, which implies an initial state that is not perfectlyrelaxed. While this numerical artifact could be avoided byexplicitly relaxing the initial conditions, we note that theseresidual gas velocities do not have any significant bearingon our results: the total gas kinetic energy E kin is less than0 .
005 of either the total potential or the total internal gas en-ergy at the final time, as shown in Figure A2 in Appendix A.It is, however, interesting to note that while the total σ gas isvery similar between GADGET and
AREPO , there are somesystematic differences in the radial profiles of σ gas , whichare caused by dissipation of gas motions on different spatialscales (see bottom panel of Figure A1). We now analyze how differences in shock capturing between
GADGET and
AREPO affect the radial infall of gas in darkmatter haloes, a problem of direct cosmological interest.For this purpose we intentionally adopt a set-up as sim-ple as possible in order to isolate differences between thecodes driven by the shock treatment only. We initially setup gas in hydrostatic equilibrium and at rest within a staticHernquist potential with mass M vir = 10 M (cid:12) , scale length a = 176 kpc, and a gas fraction of f gas = 0 .
17. We consider c (cid:13) , 000–000 Sijacki et al. ρ g a s ( r ) [ M O • k p c - ] Softened Hernquist Halo Softened Hernquist Halo Softened Hernquist Halo Softened Hernquist Halo Softened Hernquist Halo Softened Hernquist Halo T [ K ] r [ kpc ] S ~ T / ρ / Figure 7.
Radial profiles of gas density, temperature and entropyat time t = 0 . N gas = 10 and r soft = 14 . N gas = 10 and r soft = 6 . N gas = 10 and r soft = 3 . r [ kpc ] S ~ T / ρ / Analytic Hernquist Halo Analytic Hernquist Halo Analytic Hernquist Halo Analytic Hernquist Halo Analytic Hernquist Halo Analytic Hernquist Halo
Figure 8.
Radial profiles of gas entropy at time t = 0 . N gas = 10 (dotted lines), N gas = 10 (dashed lines) and N gas = 10 (continuous lines),which are largely overlapping. These runs have been computedby adopting an analytic Hernquist dark matter halo without anysoftening. The black vertical dotted line denotes the virial radiusof the system. both an analytic gravitational potential and a potential witha centrally softened core ¶ . We introduce a small modifica-tion in the codes such that gas self-gravity is switched off –the gas only feels the static dark matter potential and hydro-dynamical forces in a purely non-radiative regime. We thenartificially reduce the internal energies of gas particles/cellsso as to displace the gas from the equilibrium solution. Thenewly assigned gas temperature is 4 . × K, equal for allresolution elements. This test is hence analogous to the well-studied Evrard collapse (Evrard 1988), but it is even simplerin nature because we intentionally neglect gas self-gravity.As a consequence of the dramatic reduction in its tem-perature, the gas will suddenly lose pressure support andradially free-fall towards the centre. As the gas collapses, aradial shock develops in the centre, steepening while it prop-agates outwards and ploughing trough the remainder of theouter material which is still falling in. The Mach number ofthe shock varies over the range ∼ −
8, well matched tothe 1D shock tube test problem described in Section 3.1.1.Finally, as the shock propagates beyond the virial radiusof the halo, the gas will reach a new hydrostatic equilib-rium solution within the static dark matter potential. InFigure 7, we show radial profiles of gas density, temperatureand entropy computed at time t = 0 . GADGET , N gb = 64, α = 1 .
0: bluelines;
AREPO : red lines) we perform three runs at differentresolutions: N gas = 10 and r soft = 14 . N gas = 10 and r soft = 6 . N gas = 10 and r soft = 3 . r soft values indicate the spatial scale over which we smooth theanalytic dark matter potential. They do not necessarily rep- ¶ We modify the analytic Hernquist potential by convolving itwith the spline softened potential in the centre, as we describe indetail in Section 3.5. c (cid:13)000
AREPO : red lines) we perform three runs at differentresolutions: N gas = 10 and r soft = 14 . N gas = 10 and r soft = 6 . N gas = 10 and r soft = 3 . r soft values indicate the spatial scale over which we smooth theanalytic dark matter potential. They do not necessarily rep- ¶ We modify the analytic Hernquist potential by convolving itwith the spline softened potential in the centre, as we describe indetail in Section 3.5. c (cid:13)000 , 000–000 he hydrodynamics of galaxy formation resent the minimum spatial resolution of these simulations,given that the gas is not self-gravitating and that the darkmatter halo is rigid. The smoothing lengths of gas particlesin the central region are of the order of r hsml ∼ GAD-GET for the lowest resolution run with N gas = 10 , while thetypical central cell sizes in AREPO are a few hundred pc.From Figure 7 it can be seen that there are systematicdifferences in the gas properties predicted by the simulationswith
GADGET and
AREPO . In particular, the gas entropydistribution is broader in
GADGET both in pre- and post-shock gas in all three simulations, while for the two lowerresolution runs there is a slight mismatch in the exact po-sition of the shock front and in its strength (similar to thefindings for the Evrard collapse in Springel 2010a), which isminimized in the case of N gas = 10 particles. Note that thedifferences in shock front position between different resolu-tion simulations for a given code are not driven by resolu-tion effects but by different spatial softening of the centralpotential, which affects the gas collapse in the innermost re-gions. We have checked this explicitly by performing runswith N gas = 10 , and 10 for both codes, but this timesimulating cold gas collapse within an analytic Hernquistpotential, as shown in Figure 8. In this case the shock prop-erties for different resolution runs are almost identical for agiven code, indicating that the simulation with N gas = 10 resolution elements is in principle sufficient for capturingthe shock position accurately. Nonetheless, regardless of theresolution, differences between GADGET and
AREPO in thepre- and post-shock gas persist also in the case of analyticHernquist potentials.The differences in entropy content of pre- and post-shock gas in
GADGET and
AREPO are in part due to thelarge shock broadening in SPH, as discussed in Sections 3.1.1and 3.1.2. In fact, if we adopt N ngb = 32 (which is consideredthe minimum number still permissable in three-dimensionalsimulations) instead of N ngb = 64, the gas entropy profile in GADGET becomes less broad, but is still not as sharp as in
AREPO , indicating that simulations with a larger number ofparticles are needed in
GADGET than in
AREPO to recovershock features with the same accuracy.Moreover, there is another numerical effect leading tospatially different entropy generation in
GADGET : in theconverging subsonic part of the flow, artificial viscosity(as implemented in
GADGET ) leads to artificial dissipationwhich increases the entropy in the pre-shock gas. While thisfeature is clearly visible in Figure 40 of Springel (2010a) forthe case of the Evrard collapse, here we see that it also en-larges the central entropy in
GADGET . The reason for thisis the following: as soon as the gas is brought out of equi-librium and starts free-falling towards the centre, the gasentropy will be boosted in the central region due to an ac-tive artificial viscosity in the converging flow, creating anentropy bump that extends up to several kpc (or even sev-eral tens of kpc in the case of potentials with large cores)away from the centre, even though the shock has not fullyformed yet at this point. This artificial entropy generationis much smaller in
AREPO . As the shock forms and propa-gates outwards, it will lead to additional physical dissipationof much higher magnitude, bringing the entropy profiles of
GADGET and
AREPO into better agreement. Interestingly,in the case of the Evrard collapse, the initial difference incentral entropy profiles is minimized with time due to the gas self-gravity (as we checked explicitly by running a sim-ulation with exactly the same gas configuration but withgas self-gravity and without static dark matter potential),while it persists in our test runs even when the new equilib-rium solution of the system is reached. The central entropyis higher in
GADGET by a factor ∼ . ∼ . GADGET and
AREPO canbe aggravated in the case of non self-gravitating gas.Note, however, that the difference in the central entropyprofiles between the two codes goes in the opposite directionto what is found in non-radiative simulations of hierarchi-cally forming galaxy clusters, where the central entropy ishigher in mesh-based calculations. This indicates that accre-tion shocks during cosmological structure formation do notseem to be the likely cause of this central entropy discrep-ancy.
We now further increase the complexity of the problem byconsidering two gaseous spheres instead of one, collapsinginto one common static dark matter halo placed in-betweenthe two spheres. Each sphere has an initial spatial displace-ment from the centre of the halo. This test is similar inspirit to a number of previous works which analyzed colli-sions of two galaxy clusters in isolation (see e.g. Ricker &Sarazin 2001; Ritchie & Thomas 2002; McCarthy et al. 2007;Springel & Farrar 2007; Mitchell et al. 2009; ZuHone 2011),but here we devise a cleaner set-up in order to minimizeadditional possible numerical effects (e.g. gravitational N-body heating, differences due to gas self-gravity, etc.). Asin Section 3.3.2, we simulate a static analytic Hernquistdark matter halo (with exactly the same parameters) butinstead of one cold gaseous sphere we generate two identi-cal cold spheres separated by 1 . x -axis, andwe again neglect any radiative losses and gas self-gravity.The gaseous spheres are constructed in the same way as inSection 3.3.2 i.e. gas is first set up to be in hydrostatic equi-librium within a static Hernquist dark matter potential, andthen its internal energy is reduced to 4 . × K. For eachcode we perform runs with different gas particle numbers,i.e. N gas = 10 , N gas = 10 and N gas = 10 per sphere.Under the gravitational pull from the central dark mat-ter potential the two cold spheres collapse towards its centreand violently collide. The interesting aspect of this problemis that radial symmetry is broken and the gas interactionresults in much more complicated shock geometries. This isillustrated in Figure 9, where we show a time sequence ofprojected gas density maps (first two columns) and mass-weighted entropy maps (last two columns) for runs with N gas = 2 × performed with the two codes. While the GADGET simulation shows qualitatively similar gas struc-tures, the detailed properties substantially differ.Initially, as the spheres start to fall in, the gas is com-pressed in the centre of the halo, generating a spherical over-density, which is somewhat broader and less peaked in
GAD-GET , largely due to a poorer effective spatial resolution andnon-negligible artificial viscosity (see Section 3.3.2). Also,during this initial stage more entropy is produced in the c (cid:13) , 000–000 Sijacki et al. -300 -200 -100 0 100 200 300-300-200-1000100200300-300 -200 -100 0 100 200 300-300-200-1000100200300
GADGET -300 -200 -100 0 100 200 300 x [ kpc ] -300-200-1000100200300 y [ kp c ] -300 -200 -100 0 100 200 300-300-200-1000100200300-300 -200 -100 0 100 200 300-300-200-1000100200300 Σ g a s -300 -200 -100 0 100 200 300-300-200-1000100200300-300 -200 -100 0 100 200 300-300-200-1000100200300 AREPO -300 -200 -100 0 100 200 300 x [ kpc ] -300-200-1000100200300 -300 -200 -100 0 100 200 300-300-200-1000100200300-300 -200 -100 0 100 200 300-300-200-1000100200300 Σ g a s -300 -200 -100 0 100 200 300-300-200-1000100200300-300 -200 -100 0 100 200 300-300-200-1000100200300 GADGET -300 -200 -100 0 100 200 300 x [ kpc ] -300-200-1000100200300 y [ kp c ] -300 -200 -100 0 100 200 300-300-200-1000100200300-300 -200 -100 0 100 200 300-300-200-1000100200300 S g a s -300 -200 -100 0 100 200 300-300-200-1000100200300-300 -200 -100 0 100 200 300-300-200-1000100200300 AREPO -300 -200 -100 0 100 200 300 x [ kpc ] -300-200-1000100200300 -300 -200 -100 0 100 200 300-300-200-1000100200300-300 -200 -100 0 100 200 300-300-200-1000100200300 S g a s -300 -200 -100 0 100 200 300-300-200-1000100200300-300 -200 -100 0 100 200 300-300-200-1000100200300 GADGET -300 -200 -100 0 100 200 300 x [ kpc ] -300-200-1000100200300 y [ kp c ] -300 -200 -100 0 100 200 300-300-200-1000100200300-300 -200 -100 0 100 200 300-300-200-1000100200300 Σ g a s -300 -200 -100 0 100 200 300-300-200-1000100200300-300 -200 -100 0 100 200 300-300-200-1000100200300 AREPO -300 -200 -100 0 100 200 300 x [ kpc ] -300-200-1000100200300 -300 -200 -100 0 100 200 300-300-200-1000100200300-300 -200 -100 0 100 200 300-300-200-1000100200300 Σ g a s -300 -200 -100 0 100 200 300-300-200-1000100200300-300 -200 -100 0 100 200 300-300-200-1000100200300 GADGET -300 -200 -100 0 100 200 300 x [ kpc ] -300-200-1000100200300 y [ kp c ] -300 -200 -100 0 100 200 300-300-200-1000100200300-300 -200 -100 0 100 200 300-300-200-1000100200300 S g a s -300 -200 -100 0 100 200 300-300-200-1000100200300-300 -200 -100 0 100 200 300-300-200-1000100200300 AREPO -300 -200 -100 0 100 200 300 x [ kpc ] -300-200-1000100200300 -300 -200 -100 0 100 200 300-300-200-1000100200300-300 -200 -100 0 100 200 300-300-200-1000100200300 S g a s -300 -200 -100 0 100 200 300-300-200-1000100200300-300 -200 -100 0 100 200 300-300-200-1000100200300 GADGET -300 -200 -100 0 100 200 300 x [ kpc ] -300-200-1000100200300 y [ kp c ] -300 -200 -100 0 100 200 300-300-200-1000100200300-300 -200 -100 0 100 200 300-300-200-1000100200300 Σ g a s -300 -200 -100 0 100 200 300-300-200-1000100200300-300 -200 -100 0 100 200 300-300-200-1000100200300 AREPO -300 -200 -100 0 100 200 300 x [ kpc ] -300-200-1000100200300 -300 -200 -100 0 100 200 300-300-200-1000100200300-300 -200 -100 0 100 200 300-300-200-1000100200300 Σ g a s -300 -200 -100 0 100 200 300-300-200-1000100200300-300 -200 -100 0 100 200 300-300-200-1000100200300 GADGET -300 -200 -100 0 100 200 300 x [ kpc ] -300-200-1000100200300 y [ kp c ] -300 -200 -100 0 100 200 300-300-200-1000100200300-300 -200 -100 0 100 200 300-300-200-1000100200300 S g a s -300 -200 -100 0 100 200 300-300-200-1000100200300-300 -200 -100 0 100 200 300-300-200-1000100200300 AREPO -300 -200 -100 0 100 200 300 x [ kpc ] -300-200-1000100200300 -300 -200 -100 0 100 200 300-300-200-1000100200300-300 -200 -100 0 100 200 300-300-200-1000100200300 S g a s -300 -200 -100 0 100 200 300-300-200-1000100200300-300 -200 -100 0 100 200 300-300-200-1000100200300 GADGET -300 -200 -100 0 100 200 300 x [ kpc ] -300-200-1000100200300 y [ kp c ] -300 -200 -100 0 100 200 300-300-200-1000100200300-300 -200 -100 0 100 200 300-300-200-1000100200300 Σ g a s -300 -200 -100 0 100 200 300-300-200-1000100200300-300 -200 -100 0 100 200 300-300-200-1000100200300 AREPO -300 -200 -100 0 100 200 300 x [ kpc ] -300-200-1000100200300 -300 -200 -100 0 100 200 300-300-200-1000100200300-300 -200 -100 0 100 200 300-300-200-1000100200300 Σ g a s -300 -200 -100 0 100 200 300-300-200-1000100200300-300 -200 -100 0 100 200 300-300-200-1000100200300 GADGET -300 -200 -100 0 100 200 300 x [ kpc ] -300-200-1000100200300 y [ kp c ] -300 -200 -100 0 100 200 300-300-200-1000100200300-300 -200 -100 0 100 200 300-300-200-1000100200300 S g a s -300 -200 -100 0 100 200 300-300-200-1000100200300-300 -200 -100 0 100 200 300-300-200-1000100200300 AREPO -300 -200 -100 0 100 200 300 x [ kpc ] -300-200-1000100200300 -300 -200 -100 0 100 200 300-300-200-1000100200300-300 -200 -100 0 100 200 300-300-200-1000100200300 S g a s Figure 9.
Projected surface density maps (first two columns; in units of [M (cid:12) kpc − ]) and mass-weighted entropy maps (last two columns;in internal units) for GADGET and AREPO simulations, showing the collision of two gaseous spheres with N gas = 2 × at times t = 0 .
54 Gyr (first row), t = 1 . t = 2 . t = 15 Gyr (fourth row). The plotted spatial domain is0 . × . × central region in our SPH calculation. As more gas falls in,a shock develops which rapidly assumes a cocoon-like ge-ometry elongated perpendicular to the direction of collapse(see top panels of Figure 9). From the gas density maps itcan be seen that the shock front is narrower and sharperin AREPO , whereas in
GADGET it has a more splotchy-like appearance, caused by kernel averaging. High entropyplumes propagating outwards along the y -axis are clearlyvisible in the right-hand panels, where the central entropyin GADGET within ∼
250 kpc is still higher than in themoving mesh simulation. As the cocoon propagates againstthe infalling material, gas in the very centre is pushed per-pendicular to the x -axis, generating a dense sheet-like region (see second row of Figure 9). Dynamical fluid instabilities atthe boundary of this dense region induce typical mushroom-like morphologies (cap-like in projection). However, even forour highest resolution simulations with N gas = 2 × par-ticles, there are some marked differences between the twocodes. Mushroom-like features (corresponding to the red-orange colours in the density maps) originating at the veryboundary of the dense central region are more coherent in GADGET , while in
AREPO they break up and mix more ef-ficiently with the surrounding medium. This also leads tothe more efficient mixing of different entropy gas in the verycore in the moving mesh simulation, as can be seen from theentropy maps. c (cid:13)000
AREPO they break up and mix more ef-ficiently with the surrounding medium. This also leads tothe more efficient mixing of different entropy gas in the verycore in the moving mesh simulation, as can be seen from theentropy maps. c (cid:13)000 , 000–000 he hydrodynamics of galaxy formation
10 100 1000 r [ kpc ] S ~ T / ρ / Figure 10.
Radial entropy profiles at t ∼
15 Gyr from the start ofthe simulation when the system has reached equilibrium. For eachcode (GADGET: blue lines; AREPO: red lines) three simulationswith increasing gas particle number are shown: N gas = 2 × (dotted lines), N gas = 2 × (dashed lines), N gas = 2 × (con-tinuous lines). The vertical dotted line indicates the virial radiusof the underlying dark matter halo. While for N gas ≥ × theentropy profiles seem converged for each code, they converge to avery different result. The differences in the fluid properties at the early stagesof the simulated system as described above are, however,not the only reason why the thermodynamic properties ofthe gas are systematically discrepant between the two codeswhen a new equilibrium state is reached. With time, denseshells of gas completely disperse by mixing with the sur-rounding material in
AREPO , while in the case of
GAD-GET filaments and blobs of dense gas survive and gradu-ally sink back to the centre (see third row of Figure 9).This buoyantly-driven deposition of low entropy materialin
GADGET causes even larger differences between the fi-nal entropy distributions, which are illustrated in Figure 10(see also bottom row of Figure 9). The radial entropy pro-files are computed once the system has reached hydrostaticequilibrium at time t ∼
15 Gyr from the start of the simula-tion. The blue lines denote
GADGET results at three differ-ent resolutions (dotted lines: N gas = 2 × , dashed lines: N gas = 2 × , solid lines: N gas = 2 × ), while the redlines are for the runs with AREPO . For N gas ≥ × ,both codes seem to produce converged entropy profiles, butthey converge to very different results. While in GADGET the entropy profiles steadily decrease towards the centre,in the moving mesh code a large entropy core is produced.This systematic difference in the central entropy profiles isin good agreement with the previous study by Mitchell et al.(2009) of idealized major merger simulations of haloes in thenon-radiative regime.
We now devise a numerical experiment to capture the evo-lution of cold, dense blobs in a more realistic setting, ratherthan a uniform density windtunnel (see Section 3.2.1). Forthis purpose we consider our default isolated halo with astatic, analytic Hernquist profile for the dark matter com-ponent, gas in hydrostatic equilibrium, and we include gasself-gravity, but neglect any radiative losses. We additionallypopulate this halo with 10 blobs, with the following prop-erties: the gas pressure within blobs is set to be 0 .
01 of themaximum intracluster pressure within 800 kpc radius fromthe centre; the radius of each blob is R blob = 20 kpc; thetotal mass of all blobs is 20% of the total intracluster gasmass, i.e. M blobs = 0 . f gas M halo ; the adiabatic index is thesame for the blobs and for the surrounding gas, γ = 5 / v r (only inward radial ve-locity), v θ , and v φ , assuming a random distribution for eachvelocity component starting from a characteristic velocityvalue of 200 km s − . Blob velocity values range then from ∼
230 km s − to ∼
510 km s − , with the average velocity ofall blobs being v mean ∼
400 km s − . In this way, individualblobs will not reach the cluster centre all at the same time,giving them more realistic orbits than purely radial ones. Ini-tially, the blobs are roughly in pressure equilibrium with thesurrounding gaseous halo. We perform this numerical exper-iment at two different resolutions: N gas = 10 , N blob = 10 (low resolution simulation), and N gas = 10 , N blob = 10 (higher resolution simulation).In Figure 11, we show the time evolution of the denseblobs moving through the isolated halo, for the higher reso-lution run. In the top panels, the projected surface densitymap is plotted at the initial time, where both the intraclus-ter gas structure and the blob properties are exactly thesame in the two codes. Initially, for t ≤ GADGET and
AREPO ,and they also have very similar morphologies. As the blobsstart approaching the inner cluster region, well-defined bowshocks develop ahead of each blob and ram-pressure strip-ping ablates the blobs. Additionally, the Kelvin-Helmholtzand Rayleigh-Taylor instabilities arise which tend to dis-rupt the blobs on a characteristic timescale of several Gyrs,as discussed in Section 3.2.1 (but note that here gas is self-gravitating). The lower rate of ram pressure stripping andsuppression of fluid instabilities in
GADGET has a signifi-cant effect not only on blob morphologies, but also on theirorbits. At t = 1 .
37 Gyr, as illustrated in the middle panels,it is still possible to cross-identify each blob in
GADGET with the respective blob in the
AREPO run, but the blobsin
GADGET have lost less material and thus appear denserand some of them are closer to the centre.This different loss of blob material leads to systemati-cally diverging orbits due to higher buoyancy and dynam-ical friction forces acting on each blob in
GADGET . Thiscan be clearly seen in the bottom panels where the blobsin
GADGET are markedly more concentrated and have es- c (cid:13) , 000–000 Sijacki et al. -1000 -500 0 500 1000-1000-50005001000-1000 -500 0 500 1000-1000-50005001000
GADGET -1000 -500 0 500 1000 x [ kpc ] -1000-50005001000 y [ kp c ] -1000 -500 0 500 1000-1000-50005001000 Σ g a s
50 100 150 20050 100 150 200 -1000 -500 0 500 1000-1000-50005001000-1000 -500 0 500 1000-1000-50005001000
AREPO -1000 -500 0 500 1000 x [ kpc ] -1000-50005001000 y [ kp c ] -1000 -500 0 500 1000-1000-50005001000 Σ g a s
50 100 150 20050 100 150 200 -1000 -500 0 500 1000-1000-50005001000-1000 -500 0 500 1000-1000-50005001000
GADGET -1000 -500 0 500 1000 x [ kpc ] -1000-50005001000 y [ kp c ] -1000 -500 0 500 1000-1000-50005001000 Σ g a s
50 100 150 20050 100 150 200 -1000 -500 0 500 1000-1000-50005001000-1000 -500 0 500 1000-1000-50005001000
AREPO -1000 -500 0 500 1000 x [ kpc ] -1000-50005001000 y [ kp c ] -1000 -500 0 500 1000-1000-50005001000 Σ g a s
50 100 150 20050 100 150 200 -1000 -500 0 500 1000-1000-50005001000-1000 -500 0 500 1000-1000-50005001000
GADGET -1000 -500 0 500 1000 x [ kpc ] -1000-50005001000 y [ kp c ] -1000 -500 0 500 1000-1000-50005001000 Σ g a s
50 100 150 20050 100 150 200 -1000 -500 0 500 1000-1000-50005001000-1000 -500 0 500 1000-1000-50005001000
AREPO -1000 -500 0 500 1000 x [ kpc ] -1000-50005001000 y [ kp c ] -1000 -500 0 500 1000-1000-50005001000 Σ g a s
50 100 150 20050 100 150 200
Figure 11.
Projected surface density maps in units of [M (cid:12) kpc − ] at the initial time (top panels), at t = 1 .
37 Gyr (middle panels), andat t = 2 .
65 Gyr (bottom panels) for GADGET ( α = 1 . N ngb = 64) and AREPO. The plotted spatial domain is 2 × × R = 755 kpc. While in the beginning the orbits and the morphologies of theblobs are very similar in the two codes, for t ≥ (cid:13) , 000–000 he hydrodynamics of galaxy formation -1000 -500 0 500 1000-1000-50005001000-1000 -500 0 500 1000-1000-50005001000 GADGET -1000 -500 0 500 1000 x [ kpc ] -1000-50005001000 y [ kp c ] -1000 -500 0 500 1000-1000-50005001000 -1000 -500 0 500 1000-1000-50005001000-1000 -500 0 500 1000-1000-50005001000 AREPO -1000 -500 0 500 1000 x [ kpc ] -1000-50005001000 y [ kp c ] -1000 -500 0 500 1000-1000-50005001000 Figure 12.
Spatial distribution of gas at t = 1 .
37 Gyr that was contained in the blobs at time t = 0. The plotted spatial domain is2 × × -400 -200 0 200 400-400-2000200400 -400 -200 0 200 400-400-2000200400 GADGET -400 -200 0 200 400 x [ kpc ] -400-2000200400 y [ kp c ] -400 -200 0 200 400-400-2000200400 V o r ti c it y -400 -200 0 200 400-400-2000200400 -400 -200 0 200 400-400-2000200400 AREPO -400 -200 0 200 400 x [ kpc ] -400-2000200400 y [ kp c ] -400 -200 0 200 400-400-2000200400 V o r ti c it y Figure 13.
Projected mass-weighted vorticity maps in units of [km s − kpc − ] (absolute value of the z -component only) for a run withGADGET (left-hand panel) and with AREPO (right-hand panel). The maps are computed at time t = 1 .
37 Gyr from the start of thesimulation. The thickness of the projected region is ∆ z = 1Mpc. In the wake of the bow shocks produced by the moving blobs turbulenceis generated. The spatial extent of the turbulent wakes is significantly larger in AREPO. sentially all fallen to the innermost cluster region, while in AREPO they are at much larger cluster-centric distances be-ing gradually eroded. Moreover, given that the blobs in themoving mesh calculation are ablated more efficiently andthat therefore the dynamical friction exerted on the blobsis lower, they have higher velocities when passing at the pericentre and thus can reach larger distances after the firstpassage, as visible in the bottom panels of Figure 11. Asthe gas is self-gravitating, the blobs are also subject to tidalstripping when passing close to the innermost regions whichcontributes to the ablation of the blob material.To illustrate more clearly the differences in mass loss of c (cid:13) , 000–000 Sijacki et al. the blobs in
GADGET and
AREPO , in Figure 12 we showprojected maps of gas material at t = 1 .
37 Gyr (correspond-ing to the middle panels of Figure 11) which initially be-longed to the blobs. To compute these maps in the case of
GADGET , we show integrated mass along the line-of-sight(∆ z = 2 Mpc) for all particles initially contained in theblobs. In the case of AREPO , we use a tracer field to fol-low spreading of fluid elements initially within the blobs:each cell is characterized by an additional scalar field
Tracer which at t = 0 is equal to 1 for all blob cells and 0 else-where. The tracer field essentially evolves as a dye cast onthe moving fluid, and at some time t > Tracer value in-dicates the mass fraction of the material which initially wasin the blobs for each cell. In the right-hand panel of Fig-ure 12, we plot the projected density-weighted tracer field.The dynamical range is the same in both panels, rangingfrom the maximum of the projected quantity to 10 − of thismaximum value, using a logarithmic colour mapping. As an-ticipated, cold dense blobs in GADGET lose much less ma-terial while in the run with
AREPO the lost blob material issignificantly more diffuse. In the wake of the infalling blobsprominent tails develop, extending up to several 100 kpc.The mass deposition of blob material is clearly much morespatially extended in the moving-mesh code and occurs overlarger cluster-centric distances, allowing fluids with differententropies to intermingle and to affect host halo propertieson larger scales.The formation of bow shocks in front of the movingblobs also implies that vorticity will be generated due tothe baroclinic term. To explore this issue in more detail, inFigure 13 we show vorticity maps for the run with
GAD-GET (left-hand panel) and for the simulation with
AREPO (right-hand panel) at t = 1 .
37 Gyr (corresponding to themiddle panels of Figure 11). The vorticity maps are con-structed by first evaluating projected mass-weighted veloc-ity maps for the x and y -components, and then by takingthe absolute value of the z -component of the curl. Note thatwhile the positions and the structure of the bow shocks arerather similar in the two codes, as expected, there is a strik-ing difference between the spatial extent of the high vor-ticity regions produced in the wake of the bow shocks (de-noted with green colours), where turbulent motions shouldbe generated. For example, focusing on the right-most blobcentred on x ∼
250 kpc and y ∼ −
100 kpc in projection, themean projected vorticity value in its wake is a factor of ∼ GADGET willhave an impact on the level of turbulence injected by curvedshocks that are associated with galaxy or subhalo motionsthrough the intracluster medium (ICM), producing a biasin the amount of non-thermal pressure support in galaxyclusters (see also e.g. Dolag et al. 2005; Vazza et al. 2011;Iapichino et al. 2011).As the blobs orbit at larger cluster-centric distances in
AREPO for a longer time and mix with the surroundingmedium more efficiently than is the case for
GADGET , theylead to higher entropy generation over a wider range of radii.This is shown in Figure 14, where we compute radial entropyprofiles for all intracluster gas including the blobs, at the fi-nal time t ∼
10 Gyr when the system has reached hydrostaticequilibrium. Both in the low and high resolution runs with
AREPO the gas entropy profile is significantly higher in the r [ kpc ] S ~ T / ρ / Figure 14.
Radial entropy profiles at time t ∼
10 Gyr after thestart of the simulation, at which point all the blobs have been dis-rupted and the system has reached a new equilibrium. For eachcode (GADGET: blue lines; AREPO: red lines) two resolutionsimulations are shown: N gas = 10 and N blob = 10 (dashedlines), N gas = 10 and N blob = 10 (continuous lines). Greycurves with the same line-styles represent the initial entropy pro-files of these simulations. Additional green line shows AREPOsimulation with N gas = 10 and N blob = 10 , but with the fixedgravitational softening as in GADGET. AREPO runs with adap-tive and fixed softenings for the gas produce very similar entropyprofiles. Vertical lines with the same line-styles denote the gravi-tational softenings and vertical dotted line the virial radius of thesystem. The interaction of the moving dense blobs with the intra-cluster medium leaves systematically different imprints on the gasentropy profiles when simulated with AREPO or with GADGET. inner regions up to r ∼
100 kpc. Instead, in
GADGET , theentropy content of the dense blobs is increased less as theysink towards the centre, such that they settle on a loweradiabat, corresponding to smaller radii. We also note thatthese systematic differences in entropy profiles are entirelydue to the different hydro solvers employed by
GADGET and
AREPO , and not due to the different choice of gravita-tional softenings for gas (fixed in
GADGET and adaptive in
AREPO , but with a floor set equal to the
GADGET value).In fact, in Figure 14, the green line shows the entropy profileobtained from an identical low resolution
AREPO simulationwhere instead of an adaptive a constant gravitational soft-enings is used for the gas, in the same way as in
GADGET .It is also instructive to analyze how the entropy profilesof the intracluster medium evolve with time. In Figure 14 weshow radial entropy profiles at the initial time (grey lines)for both resolutions. The initial entropy profile of the runwith a higher particle number extends further inwards dueto the better spatial resolution, while some differences inthe outer regions are due to different positions of the blobswhich are drawn randomly. The kink in the initial entropyprofile for r ≥
500 kpc is caused by the blobs which popu-late this region and have low entropy content with respectto the surrounding ICM. During the first 0 . c (cid:13)000
500 kpc is caused by the blobs which popu-late this region and have low entropy content with respectto the surrounding ICM. During the first 0 . c (cid:13)000 , 000–000 he hydrodynamics of galaxy formation tropy profiles in GADGET and
AREPO do not change no-ticeably, remaining nearly identical to each other and to thevalues prescribed by the initial conditions. At 1 Gyr, thegas entropy starts rising in the very centre in the run with
AREPO (it is higher by a factor of ≤ GADGET before the central entropy rises by thesame amount. At that time, however, the central entropyprofile in
AREPO is already about a factor of 6 higher thanit was initially, and this systematic difference in the centralentropy values persists until the final simulated time.At time t ∼ GADGET and
ENZO codes show that infalling satellites in
GADGET sink to smaller radii and are characterized by lowerentropy content than is the case for
ENZO , confirming theimportance of the physical processes discussed here in thefull cosmological setting. This indicates that at least partof the systematic difference in central entropy values foundbetween SPH and grid codes in the Santa Barbara clustercomparison project (Frenk et al. 1999) is due to the differenttreatment of stripping and mixing of the cold material fromthe infalling satellites.
We now investigate the time evolution of an isolated 10 M (cid:12) mass halo subject to radiative cooling and star formation.We first simulate a halo without any net rotation and com-pare gas and stellar properties of this system between thetwo codes. The results illustrated here are for haloes withstatic dark matter potentials, but we find very similar resultsif we consider live haloes instead.In Figure 15 we show star formation rates as a functionof time for our simulated object consisting of N gas = 10 or N gas = 10 particles/cells, calculated either with GADGET or AREPO . Over the whole simulated time-span of 10 Gyrthe star formation rates are very similar in the two codes.This indicates that not only do the implemented gas coolingrates match very well, but also that the sub-grid model forstar formation and supernova feedback is consistent betweenthe codes, even at a relatively low numerical resolution (asimilar conclusion is reached for mergers of isolated galaxiesby Hayward et al. 2012).To demonstrate more clearly that gas radiative cool-ing and star formation proceed in a very similar mannerin our isolated haloes in the absence of any net rotation,in Figure 16 we plot the total mass of cold baryons (starsand gas with entropy < in internal units) as a functionof time for GADGET (blue lines) and
AREPO (red lines),at two different resolutions, i.e. N part = 10 (dotted lines)and N part = 10 . For the same number of SPH particles ascells used in AREPO the amount of cold baryons matches to t [Gyr] S FR [ M O • / y r ] GADGET LR GADGET HR AREPO LR AREPO HR
Figure 15.
Time evolution of the star formation rate of a M vir = 10 M (cid:12) isolated halo which radiatively cools and hasnegligible spin. Illustrated are low ( N gas = 10 ; dotted lines) andhigh ( N gas = 10 ; continuous lines) resolution simulations withGADGET (blue) and AREPO (red). The agreement in star for-mation rates is very good over the whole simulated time-span. t [Gyr] M c o l d [ M O • ] ∆ M c o l d Figure 16.
Time evolution of the cold baryonic mass (stars andgas with entropy < in internal units) for a M vir = 10 M (cid:12) isolated halo which radiatively cools and has negligible spin. Illus-trated are low ( N gas = 10 ; dotted lines) and high ( N gas = 10 ;continuous lines) resolution simulations with GADGET (blue)and AREPO (red). In the inset, the ratio of M cold values foundin AREPO and GADGET is shown for both numerical resolu-tions. The agreement of the amount of cold baryons formed isextremely good.c (cid:13) , 000–000 Sijacki et al. log [ r gas / kpc] l og [ S = u / ρ / / i n t e r n a l un it s ] t = 5Gyr GADGET AREPO Figure 17.
2D histogram of the gas entropy as a function of radialdistance from the centre for GADGET (blue) and AREPO (red)at time t = 5 Gyr for N part = 10 (dotted lines) and N part = 10 (continuous lines). The gaseous halo is allowed to radiatively coolbut there is no net rotation. The vertical dotted lines indicate thegravitational softening used in the GADGET run, which corre-spond to the floor values of the adaptive softenings in the movingmesh calculations. within a few percent between the two codes, as is evidentfrom the inset in the plot where we show the ratio of M cold values found with our moving-mesh code and with GAD-GET . In the case of live haloes we find that the total massof cold baryons exhibits the same level of agreement.Furthermore, in Figure 17 we show a 2D histogram ofintracluster gas entropies as a function of cluster-centric dis-tance, after 5 Gyr from the start of the simulation for thelow resolution and high resolution runs. The distribution ofgas entropies is almost identical in the two codes for r > r soft (note that below r soft the simulation results are not trust-worthy due to the limited gravitational resolution on thesescales). It can be seen that there is some difference in thegas entropy close to the virial radius of the system, which iscaused by different boundary conditions (vacuum for GAD-GET and a uniform low resolution grid for
AREPO ).These results are in line with the expectation that gascooling and condensation in this simulated system is de-termined entirely by the gas properties at the cooling ra-dius (Bertschinger 1989; White & Frenk 1991; Hernquist &Springel 2003) which corresponds very closely between thesimulation codes.
Even though gas cooling and star formation proceed in aremarkably similar way in
GADGET and the moving-meshcode for haloes with vanishing spins, this is not guaranteedto remain the case once some degree of rotation is included.We therefore simulate exactly the same isolated haloes asin the previous section, but imposing a certain level of gas t [Gyr] M c o l d [ M O • ] ∆ M c o l d Figure 18.
Time evolution of the cold baryonic mass (starsand gas with entropy < in internal units) for a M vir =10 M (cid:12) isolated halo which radiatively cools and has a spin of λ = 0 .
4. Illustrated are low ( N gas = 10 ; dotted lines) and high( N gas = 10 ; continuous lines) resolution simulations with GAD-GET (blue) and AREPO (red). In the inset, the ratio of M cold values found in AREPO and GADGET is shown for both nu-merical resolutions. The agreement between the amount of coldbaryons formed is poorer for N gas = 10 , but improves for thehigh resolution run, with M cold in AREPO being ∼ − rotation within the static dark matter potential. To highlightthe effect, we use a large spin parameter equal to λ = 0 . < in inter-nal units) for the rotating haloes simulated with GADGET (blue curves) and
AREPO (red curves) at different resolu-tions. Comparing M cold with our findings from Figure 16 fornon-rotating haloes indicates that overall less gas cools fromthe hot phase if the gas spins. This effect is not surprisinggiven our initial conditions. Even though the gas densityand temperature distribution are initially identical, in thesimulations where there is considerable spin the gas will besubject to centrifugal accelerations, preventing it from col-lapsing radially. The (partial) centrifugal support will tendto reduce the gas densities and hence the cooling rates.More importantly, from Figure 18 it can be seen thatthere is poorer agreement in the amount of cold baryonsbetween GADGET and
AREPO for a rotating gaseous halo.The discrepancy is larger for the low resolution run with N gas = 10 , where the final M cold value after 10 Gyr is about30% higher in the moving-mesh code. The reason for thisdiscrepancy is twofold: at low resolution, GADGET under-estimates the cooling rate somewhat, while
AREPO overes-timates it. The SPH result turns out to be more stable atpoor resolution than the mesh-based calculation. Here it isadvantageous for SPH that even for few particles a clearlydefined phase boundary between cold and hot gas is main-tained (in fact, a ‘pressure blip’ in SPH leads to a sampling c (cid:13)000
AREPO overes-timates it. The SPH result turns out to be more stable atpoor resolution than the mesh-based calculation. Here it isadvantageous for SPH that even for few particles a clearlydefined phase boundary between cold and hot gas is main-tained (in fact, a ‘pressure blip’ in SPH leads to a sampling c (cid:13)000 , 000–000 he hydrodynamics of galaxy formation log [ r gas / kpc] l og [ S = u / ρ / / i n t e r n a l un it s ] LOW RES t = 5Gyr GADGET AREPO log [ r gas / kpc] l og [ S = u / ρ / / i n t e r n a l un it s ] HIGH RES t = 5Gyr GADGET AREPO
Figure 19.
2D histograms of gas entropy as a function of radial distance from the centre for GADGET and AREPO at t = 5 Gyr forthe 10 M (cid:12) mass halo which radiatively cools and spins at low ( N gas = 10 ; left-hand panel) and high ( N gas = 10 ; right-hand panel)resolution. The vertical dotted line is the gravitational softening used in the GADGET run, which correspond to the floor value of theadaptive softenings in the moving mesh calculations. gap at this boundary) whereas in AREPO this boundaryis blurred at low resolution, causing slightly elevated cool-ing. This trend is also confirmed by the gas radial veloci-ties which are least negative in the low resolution
GADGET run and most negative in the low resolution
AREPO simula-tion. The radial velocities systematically differ within a few100 kpc and throughout most of the simulated time-span.For N gas = 10 , the radial velocities obtained with the twocodes are much closer, which also translates into a bettermatch of the amount of cold baryons, with M cold higher by ∼ −
15% in
AREPO .At least part of the remaining difference in M cold be-tween GADGET and
AREPO is driven by the interactionbetween the gas which is cooling towards the centre andthe cold gas that is already in the disk. This is shown inFigure 19, where we plot 2D histograms of the gas entropyas a function of cluster-centric distance at t = 5 Gyr forthe low (left-hand panel) and high resolution simulations(right-hand panel). In the moving-mesh code the gas sur-rounding the central disk has a range of entropy values, withsome fluid elements exhibiting very large entropies due tothe shock in the immediate vicinity of the disk. Conversely,in the runs with GADGET and especially at low resolution,there is a clear gap between cold material in the disk andhotter gas in the halo (see e.g. Agertz et al. 2007; Springel2010b). Similarly, as described in Section 3.3.2, SPH par-ticles are affected by artificial viscosity in the convergingpart of the flow which can slightly offset cooling losses. Notethat for N gas = 10 the difference in gas entropy structurearound the disk is lower between the codes, as expected,given that unwanted artificial viscosity effects are reducedand that the gap between the cold and hot gas phase due torepulsive pressure forces is smaller.From Figure 19 is it also evident that the extent of the cold disk is different between the codes, especially at lowresolutions. To quantify this important effect, in Table 1we summarize the main properties of the forming disk. For N gas = 10 the half-mass radius of the gaseous disk in GAD-GET is almost a factor of two smaller than in the higherresolution run, while in the moving mesh code R gas , HM is70% of the value we obtain with N gas = 10 . This indi-cates that the convergence rate of the gas disk size is slowerin the case of GADGET , due to spurious transfer of angu-lar momentum from the cold to the hot phase (Okamotoet al. 2003) and due to the artificial viscosity in the case ofpoorly resolved disks (note, however, that the total angu-lar momentum is manifestly conserved in
GADGET ). In thecase of the stellar disks they are essentially not resolved inour low resolution runs ( R stars , HM ∼ r soft = 14 kpc), whilethe half-mass radius is ∼
40% higher in
AREPO at higherresolution. Contrary to
GADGET , total angular momentumconservation is not automatically guaranteed in the movingmesh code, particularly for disks resolved with only a smallnumber of cells. Nonetheless, it is reassuring that R gas , HM increases with higher resolution in AREPO (attesting thatspurious angular momentum transport inwards is probablysmall) and that the value of R gas , HM obtained with bothcodes for N gas = 10 is relatively close. This indicates thatthe differences in the disk sizes in this numerical experi-ment are at least partly driven by resolution effects, and inparticular by the slow convergence rate of GADGET sim-ulations. Comparing the resolution requirements to obtaingood agreement between the codes here with Section 3.4we note that in the case of cooling gaseous haloes with asignificant spin more resolution elements are needed to fol-low accurately the thermo-dynamical interaction betweenthe central cold disk and the inflowing hotter gas. c (cid:13) , 000–000 Sijacki et al. code N gas R gas , max R gas , HM M gas , HM R stars , HM M stars , HM M cold [kpc] [kpc] [10 M (cid:12) ] [kpc] [10 M (cid:12) ] [10 M (cid:12) ]GADGET 10 . . . . . . . . . . . . . . . . . . . . . . . . Table 1.
Gaseous and stellar properties of the disk at t = 5 Gyr which forms in an isolated 10 M (cid:12) halo which radiatively cools androtates. For simulations with two different resolutions with GADGET and AREPO, we list the maximum radial extent of the gaseousdisk, its half-mass radius and the total mass within the half-mass radius, in the third, fourth and fifth columns, respectively. Stellarhalf-mass radius and the total stellar mass enclosed within it are shown in the sixth and seventh columns. In the eighth column the totalmass of cold baryons is given, being equal to 2 × ( M gas , HM + M stars , HM ). All gas particles/cells with entropy less than 10 (in internalunits) are assumed to be part of the disk, while for the stellar disk we consider all star particles that form in the simulated volume.code N gas N blob R gas , max R gas , HM M gas , HM R stars , HM M stars , HM M cold [kpc] [kpc] [10 M (cid:12) ] [kpc] [10 M (cid:12) ] [10 M (cid:12) ]GADGET 10 . . . . . . . . . . . . . . . / / . . . . . . . . . . . . . . . . . . . . . . / / . . . . . . . Table 2.
Gaseous and stellar properties of the disk at t = 5 Gyr that formed in an isolated 10 M (cid:12) halo which radiatively cools, rotatesand contains 10 orbiting substructures. For simulations with three different resolutions with GADGET and AREPO we list the maximumradial extent of the gaseous disk, its half-mass radius and the total mass within the half-mass radius, in the fourth, fifth and sixth columns,respectively. Stellar half-mass radius and the total stellar mass enclosed within it are shown in the seventh and eighth columns. In theninth column the total mass of cold baryons is given, being equal to 2 × ( M gas , HM + M stars , HM ). All gas particles/cells with entropyless than 10 (in internal units) and within R gas , max are assumed to be part of the disk, while for the stellar disk we consider all starparticles that form in the simulated volume (including the blobs). We also list the gaseous disk properties in two additional simulations,denoted by EOS, performed at intermediate resolution. In these simulations we adopt the standard sub-grid model for star formation,with the dense, cold gas lying on the effective equation-of-state, but we prevent any spawning of star particles. Here we consider the evolution of cold, dense blobs embed-ded in a galaxy cluster as in Section 3.3.4, but now thewhole system is allowed to radiatively cool. To add an addi-tional layer of realism the blobs are constructed to be similarto cosmological substructures: they are equipped with theirown dark matter halo, and stars may form out of their gasduring the simulated time-span. More specifically, each blobis represented by a live Hernquist dark matter halo and gasin hydrostatic equilibrium. The virial mass of each blob is2 × M (cid:12) , the scale length parameter is a = 41 . f gas = 0 .
17. As before, we popu-late our default 10 M (cid:12) halo (that has a static dark matterpotential) with 10 identical substructures. The procedure forassigning blob positions and velocities is the same as in Sec-tion 3.3.4, but here we use a different random number seedwhich leads to the different initial positions and velocitieswith respect to Section 3.3.4. The positions of blob centresare in the range of ∼
650 to ∼
750 kpc, while the character-istic blob velocities range from 200 to 500 km s − . Also, inSection 3.3.4 the gas in the halo is initially at rest whilehere the gas has considerable angular momentum, whichcontributes to the relative velocity between the blobs and the gas. We simulate the evolution of this system at threedifferent resolutions: N gas = 10 , N blob = 10 (low resolu-tion), N gas = 10 , N blob = 10 (intermediate resolution),and N gas = 10 , N blob = 10 (high resolution), with gasself-gravity, and with cooling and the sub-grid model forstar formation.Additionally, we perform two simulations at intermedi-ate resolution where the standard sub-grid model for starformation is included as well, but where we simply pre-vent any star particles from being spawned out of the star-forming phase, denoted by EOS. This simulation set-up isparticularly useful for following the thermodynamical evo-lution of cold and hot gas for many Gyrs without the densecold gas being subject to fragmentation that would likelyoccur in pure cooling runs. For numerical experiments with AREPO , as a default choice, we have not used mesh refine-ment and de-refinement. However, we have performed ex-tra runs at intermediate resolution where we adopt a de-/refinement strategy to limit the mass range of gas cells(within a factor of 2 of the gas particle mass in the matching
GADGET run) which automatically imposes a narrow rangeof stellar masses as well. For these runs we have increasedthe number of gas cells/particles in each blob to 2 . × ,such that the gas particle/cell mass in blobs is exactly iden- c (cid:13)000
GADGET run) which automatically imposes a narrow rangeof stellar masses as well. For these runs we have increasedthe number of gas cells/particles in each blob to 2 . × ,such that the gas particle/cell mass in blobs is exactly iden- c (cid:13)000 , 000–000 he hydrodynamics of galaxy formation tical to the one in the parent halo (which is optimal forour de-/refinement method). These test runs recover veryclosely all of our results with the default set-up, confirmingthat N-body heating effects (e.g. due to the spectrum of starparticle masses) are not very important here.In Figure 20, we show projected surface density maps at t = 0 . t = 2 . t = 6 Gyr for simulations wherewe use the sub-grid model for star formation, but prevent thespawning of star particles. Initially the evolution of the blobsproceeds in a very similar fashion in simulations with GAD-GET and
AREPO . However, already after less than a Gyr(see top panels of Figure 20) blobs in the moving mesh codeare much more affected by ram pressure and dynamical fluidinstabilities, which cause efficient gas stripping. As the blobsreach the very inner regions, they interact with the formingdisk. In
GADGET simulations, the blobs have a significantlymore damaging effect on the disk, simply because they areless gas depleted. This is clearly visible in the central pan-els of Figure 20, where in
GADGET feeble spiral featuresare present, while in
AREPO the disk is more extended andspiral arms are well developed. Additionally, in the moving-mesh simulation a large fraction of the stripped material isdeposited in the main halo, and is gradually accreted ontothe central disk, promoting its growth. Blobs in
GADGET are more coherent and eventually lose their angular momen-tum due to hydro-/dynamical friction. As they merge withthe central disk, an ellipsoidal structure is formed (see bot-tom panels of Figure 20). Note that due to the static darkmatter halo the dynamical friction force arises only due tothe gas in the main halo. Thus, we expect the blobs willlose their angular momentum on an even faster time-scalein simulations with live dark matter haloes, as is the casefor cosmological runs. After 6 Gyrs, the difference between
GADGET and
AREPO simulations is significant: while in thelatter case an extended disk forms, with gaseous spiral armsreaching up to 60 kpc away from the centre, in
GADGET weare left with a flattened, amorphous blob.In numerical experiments where we do not suppress starparticles from forming, we find very similar results. In Fig-ure 21, we show projected stellar density maps at t = 6 Gyrusing our intermediate resolution simulations (which bothin resolution and in time match the bottom panels of Fig-ure 20). Clearly, the stellar distribution in GADGET is morecentrally concentrated, forming a small disk. On the con-trary, a well-defined extended stellar disk assembles in themoving-mesh calculation, endowed with a central bar. Wequantify the properties of gaseous and stellar disk in Table 2.The convergence rate of the gas disk size in
GADGET is evenslower than in the simulations without blobs, because addi-tionally to the spurious transfer of the angular momentumthe damaging effect of the blobs is the highest in the low res-olution runs, where gas stripping is largely suppressed. Thisleads to larger systematic difference in disk sizes betweenthe two codes with resolution, but for the highest resolutionruns
AREPO disks are only about ∼
30% larger.The interaction of cold blobs with the surroundingmedium influences the global star formation rates as well. InFigure 22, we compare the star formation histories in simula-tions without blobs (left-hand panel; see Section 3.4.2) withthe numerical tests where haloes are populated with 10 sub-structures (right-hand panel). In the case without blobs, thestar formation rates are somewhat higher in
AREPO espe- cially for the low resolution run due to the more efficientgas cooling and thus larger amounts of cold gas availablefor star formation. For the high resolution run the differ-ences in star formation rates between the codes are smaller,given that difference in the amount of cold gas between thecodes amount to ∼ t > AREPO blobs dramatically drops, and at ∼ R blob (which we varied fromthe scale length parameter to the virial radius) and for thesurrounding medium density we take the typical halo den-sity at the position of the blobs. With these assumptionswe obtain typical values for t KH in the range of 1 − R blob . These t KH values are comparable to thetimescale on which star formation within the blobs is sup-pressed in AREPO . In
GADGET simulations, star formationproceeds in the blobs even until 9 Gyrs, albeit at a progres-sively reduced rate. For t >
AREPO simulations have higher star for-mation rates over many Gyrs, which are typically larger bya factor of 2 than the
GADGET results. Even in our highestresolution moving mesh simulation the star forming disk hasroughly twice as large amount of cold gas, a difference whichoriginates from the interaction of hotter infalling gas withthe gas in the disk, as described in Section 3.4.2. However,also the star-forming gas is distributed over a larger area. Infact, while within the half-mass radius of the
GADGET diskthe star formation rate is ∼
40% higher in the moving meshrun, outside of it is a factor ∼ . AREPO star formation rate. The reason whycold, star forming gas in
GADGET is filling a smaller area inthe disk and is more confined to the dense arm segments andblobs (the so called “string-of-pears”, which can be also seenin the middle panel of Figure 20), is due to the SPH surfacetension originating at the interface between cold and hot me-dia in relative motion. In the right-hand panel of Figure 22we also show an identical intermediate resolution
AREPO run, but with fixed gas gravitational softening, the same asin the matching
GADGET run. It can be seen that the choiceof gas gravitational softening in
AREPO does not affect ourresults in any significant way: both the evolution of the SFRin the central disk and in the orbiting blobs remains verysimilar, indicating that the differences that we see betweenthe two codes are indeed entirely driven by hydrodynamicaleffects.These findings have immediate consequences for more c (cid:13) , 000–000 Sijacki et al. -200 -100 0 100 200-200-1000100200 -200 -100 0 100 200-200-1000100200
GADGET -200 -100 0 100 200 x [ kpc ] -200-1000100200 y [ kp c ] -200 -100 0 100 200 -200-1000100200-200 -100 0 100 200-200-1000100200 Σ g a s
50 100 150 20050 100 150 200 -200 -100 0 100 200-200-1000100200 -200 -100 0 100 200-200-1000100200
AREPO -200 -100 0 100 200 x [ kpc ] -200-1000100200 y [ kp c ] -200 -100 0 100 200 -200-1000100200-200 -100 0 100 200-200-1000100200 Σ g a s
50 100 150 20050 100 150 200 -200 -100 0 100 200-200-1000100200 -200 -100 0 100 200-200-1000100200
GADGET -200 -100 0 100 200 x [ kpc ] -200-1000100200 y [ kp c ] -200 -100 0 100 200 -200-1000100200-200 -100 0 100 200-200-1000100200 Σ g a s
50 100 150 20050 100 150 200 -200 -100 0 100 200-200-1000100200 -200 -100 0 100 200-200-1000100200
AREPO -200 -100 0 100 200 x [ kpc ] -200-1000100200 y [ kp c ] -200 -100 0 100 200 -200-1000100200-200 -100 0 100 200-200-1000100200 Σ g a s
50 100 150 20050 100 150 200 -200 -100 0 100 200-200-1000100200 -200 -100 0 100 200-200-1000100200
GADGET -200 -100 0 100 200 x [ kpc ] -200-1000100200 y [ kp c ] -200 -100 0 100 200 -200-1000100200-200 -100 0 100 200-200-1000100200 Σ g a s
50 100 150 20050 100 150 200 -200 -100 0 100 200-200-1000100200 -200 -100 0 100 200-200-1000100200
AREPO -200 -100 0 100 200 x [ kpc ] -200-1000100200 y [ kp c ] -200 -100 0 100 200 -200-1000100200-200 -100 0 100 200-200-1000100200 Σ g a s
50 100 150 20050 100 150 200
Figure 20.
Projected surface density maps in units of [M (cid:12) kpc − ] at times t = 0 . t = 2 . t = 6 Gyr (bottom panels) for a M vir = 10 M (cid:12) isolated halo which rotates, radiatively cools, and has 10 orbiting substructures. Thethickness of the projection is ∆ z = 2 Mpc. Although we have used our sub-grid model for star formation in these simulations, spawningof star particles has been intentionally prevented here. Gas stripping from the orbiting blobs is found to be very different in the twonumerical techniques. In GADGET, blobs largely survive, and when they interact with the central disk they tend to disrupt it, whereasin AREPO, the blobs are more easily shredded and have a less damaging effect on the forming disk. In fact, some of the stripped blobmaterial ends up contributing to the extended gaseous disk. c (cid:13)000
Projected surface density maps in units of [M (cid:12) kpc − ] at times t = 0 . t = 2 . t = 6 Gyr (bottom panels) for a M vir = 10 M (cid:12) isolated halo which rotates, radiatively cools, and has 10 orbiting substructures. Thethickness of the projection is ∆ z = 2 Mpc. Although we have used our sub-grid model for star formation in these simulations, spawningof star particles has been intentionally prevented here. Gas stripping from the orbiting blobs is found to be very different in the twonumerical techniques. In GADGET, blobs largely survive, and when they interact with the central disk they tend to disrupt it, whereasin AREPO, the blobs are more easily shredded and have a less damaging effect on the forming disk. In fact, some of the stripped blobmaterial ends up contributing to the extended gaseous disk. c (cid:13)000 , 000–000 he hydrodynamics of galaxy formation -200 -100 0 100 200-200-1000100200 -200 -100 0 100 200-200-1000100200 GADGET -200 -100 0 100 200 x [ kpc ] -200-1000100200 y [ kp c ] -200 -100 0 100 200 -200-1000100200-200 -100 0 100 200-200-1000100200 Σ s t a r
50 100 150 20050 100 150 200 -200 -100 0 100 200-200-1000100200 -200 -100 0 100 200-200-1000100200
AREPO -200 -100 0 100 200 x [ kpc ] -200-1000100200 y [ kp c ] -200 -100 0 100 200 -200-1000100200-200 -100 0 100 200-200-1000100200 Σ s t a r
50 100 150 20050 100 150 200
Figure 21.
Projected stellar density maps in units of [M (cid:12) kpc − ] at t = 6 Gyr for a M vir = 10 M (cid:12) isolated halo which rotates,radiatively cools, forms stars, and has 10 orbiting substructures. The maps have been constructed from our intermediate resolution runsand the thickness of the projection is ∆ z = 2 Mpc. While in GADGET a centrally concentrated disk forms that is surrounded by a feeblethick structure, in the simulation with AREPO, a well-defined extended stellar disk assembles with a central bar. t [Gyr] S FR [ M O • / y r ] GADGET LR GADGET HR AREPO LR AREPO HR NO BLOBS 0 2 4 6 8 10 t [Gyr] S FR [ M O • / y r ] WITH BLOBS GADGET MR AREPO MR
SFR DISK SFR BLOBS
Figure 22.
Time evolution of the star formation rate of a M vir = 10 M (cid:12) isolated halo which radiatively cools and rotates. In theleft-hand panel we show simulations without blobs, while in the right-hand panel the haloes contain 10 orbiting substructures. In theright-hand panel, the curves give results from simulations with three different resolutions, with GADGET and AREPO, using N gas = 10 (LR), N gas = 10 (MR), and N gas = 10 (HR) particles/cells. In the left-hand panel, only low and high resolution simulations are shown.In the simulations without blobs, the star formation rates are somewhat higher in AREPO, but the relative difference increases when weinclude the blobs. Higher star formation rates in AREPO emanate from the central disk which contains a greater amount of cold, star-forming gas which fills a larger area. In the inset plot, we compute the total star formation rate occurring within the blobs (dot-dashedand continuous lines) and the one coming from the disk (triple dot-dashed and continuous lines) for our intermediate and high resolutionsimulations with GADGET (blue curves) and AREPO (red curves). In the right-hand panel we also show intermediate resolution AREPOrun with the fixed gravitational softening for the gas (green lines) which produces very similar SFRs (both in the disk and in the blobs)as the simulation with the adaptive softening.c (cid:13) , 000–000 Sijacki et al. realistic astrophysical situations. For example, Puchweinet al. (2010) have simulated a high-resolution sample ofgalaxy clusters with
GADGET finding that up to 30% of intr-acluster stars form in dense cold blobs – remnants of infallingsatellites. Suppressed dynamical instabilities in
GADGET enhance the probability of survival for these dense blobs,which can then serve as sites of star formation, thus possi-bly over-predicting the number of intracluster stars. In fact,if we compute the total mass of cold baryons M cold (starsplus gas above the density threshold for star formation) intwo cosmological simulations (for further details see Paper I)where in one we prevent the spawning of stars and in theother we allow it (the simulations are otherwise identical), M cold is found to be very similar in GADGET regardlessof whether stars are formed or not. Instead, in the simu-lations with
AREPO , we find that M cold is systematicallyreduced at lower redshifts if the spawning of stars is pre-vented. This demonstrates explicitly that very low entropymaterial formed in GADGET cannot easily be shredded andmixed with higher entropy gas (at least in the absence of ad-ditional feedback processes such as galactic winds or blackhole heating), so that M cold is preserved (a similar conclu-sion has been reached independently by Heß & Springel2011). In contrast, in the simulations with AREPO if starformation is switched-off some of the cold material can bestripped out of galaxies due to dynamical instabilities, re-turning it to diffuse form and lowering M cold . Importantly,this also implies that the number of stars formed in AREPO will be much more sensitive to the characteristic timescalefor star formation and to the numerical resolution (neededto resolve fluid instabilities) than is the case for
GADGET .In a recent paper by Agertz et al. (2011) it has been shownthat lower star formation efficiency leads to the productionof larger disks in cosmological simulations. While this resultis in line with our findings, it is important to recognize thatthe actual cause is quite different: in the study by Agertzet al. (2011), differences in the physical modelling of starformation in cosmological simulations affect the disk sizes,while here the cause lies in different accuracies of the hydrosolvers involved.
We now investigate potential artificial gas heating due tothe Poisson noise induced by the finite number of particlesin dark matter haloes. Steinmetz & White (1997) outlinedthe analytic theory of this effect, and confirmed it with non-radiative and radiative numerical experiments that quan-tified gravitational N-body noise present in structure for-mation simulations. For an equilibrium system they defineda characteristic N-body heating timescale which is propor-tional to the cube power of the dark matter velocity disper-sion and inversely proportional to the dark matter particlemass and dark matter density. Due to this inverse propor-tionality to the dark matter density, N-body heating is ex-pected to be strongest in the innermost regions of haloes.From their analysis it follows that the dark matter parti-cle mass adopted in numerical simulations should be lowerthan a critical mass which is of order of few times 10 M (cid:12) for galaxy clusters (see their equation [10]), otherwise theradiative cooling losses may be overwhelmed by spuriousheating. ρ g a s ( r ) [ M O • k p c - ] T [ K ] r [ kpc ] S ~ T / ρ / Figure 23.
Radial profiles of gas density, temperature and en-tropy at time t ∼
10 Gyr for GADGET (blue lines) and forAREPO (red lines). For each code we show two different reso-lution runs: N gas = 10 and r soft = 14 . N gas = 10 and r soft = 3 . (cid:13)000
10 Gyr for GADGET (blue lines) and forAREPO (red lines). For each code we show two different reso-lution runs: N gas = 10 and r soft = 14 . N gas = 10 and r soft = 3 . (cid:13)000 , 000–000 he hydrodynamics of galaxy formation R a d i a l A cce l e r a ti on s Hernquist S p li n e - s o f t e n e d
1 2 3 r [ kpc ] -0.10-0.050.000.050.10 S i m u l a ti on / F it -
1 2 3
Figure 24.
Comparison of fitting functions for spline-softenedgravitational accelerations (equation 4; blue dashed lines) withactual gravitational accelerations obtained from very high reso-lution runs with matching r soft values (green continuous lines).The limiting cases of spline-softened accelerations for small andlarge radii are illustrated with dotted red lines. The bottom panelshows the residuals of our fitting functions over the relevant rangeof radii. Since nowadays the typical mass resolution of hydro-dynamical cosmological simulations has improved dramati-cally, with simulated galaxy clusters containing several times10 up to 10 particles in zoom-in runs (Dolag et al. 2009;Sijacki et al. 2009; Vazza et al. 2010), possible numericalartifacts due to a grainy dark matter distribution are rarelyaddressed in the literature (but see e.g. Kay et al. 2000;Borgani et al. 2006; Vazza 2011). Recently, Springel (2010a)has pointed out that mesh based codes are potentially moreseverely affected by gravitational N-body heating than SPHcodes, due to their better ability to detect weak shocks,which in this case is an unwanted feature.Here we perform a number of idealized test problemsaimed at explicitly addressing the N-body heating problem,and in particular to understand whether there are system-atic differences between GADGET and
AREPO in this re-spect. Note that intentionally all previous tests (except forSection 3.4.3) have been performed either without any darkmatter component or with a static dark matter potential toavoid such possible Poisson-noise imprints.We first consider our standard isolated halo, where wereplace the static dark matter potential with a live dark mat-ter halo in which dark matter particles have a characteris-tic velocity dispersion for the given mass of the system (butthere is no net rotation). As in the case of the rigid potential,we populate the halo with gas particles in hydrostatic equi-librium which are initially at rest. We include self-gravityof the gas component and evolve the system non-radiativelyfor 10 Gyr.In Figure 23, we show radial profiles of gas density, tem-perature and entropy at the final time for both codes attwo resolutions: N gas = N DM = 10 particles (dotted lines) and N gas = N DM = 10 particles (continuous lines). Thecentre of the system is defined by the position of the mostbound particle. For cluster-centric distances greater than thegravitational softening value, r soft , the entropy profiles of GADGET and
AREPO agree to within ∼
20% for the lowresolution runs, and to better than 10% for the high res-olution simulation. For r < r soft , the differences betweenthe codes in the entropy profile amount up to a factor of 2,with the gas entropy in
AREPO being systematically higher.An extra test simulation performed with
AREPO in its dualentropy mode (where entropy instead of energy is explic-itly conserved for cells whose Riemann problems with ad-jacent cells all have Mach number less than 1 .
1) resultedin identical radial profiles as for the
AREPO simulation inthe standard energy mode. This clearly indicates that thesmall entropy core which develops in
AREPO is not causedby weak shocks generated by the grainy dark matter poten-tial. Rather, central gas particles/cells acquire small scalevelocities when simulated with live haloes. In the case of themoving-mesh code, these velocities lead to fluid mixing andtemperature equilibration in the innermost regions, while in
GADGET such mixing does not occur. Note, however, that inthe gravitating systems, r soft represents the minimum spa-tial scale below which the simulation results are not trust-worthy due to discreteness effects. Therefore, provided thatgravitational softening lengths are chosen conservatively, thegas properties of galaxy clusters in an equilibrium configura-tion simulated with GADGET and
AREPO match very well.We now consider a more challenging problem, where wesimulate the infall of cold gas into dark matter haloes. Thisis the same problem discussed in Section 3.3.2, but herewe compare outcomes from live versus static dark matterhaloes. To match numerical experiments with rigid and livehaloes as closely as possible we adopt the following proce-dure.For live haloes, the motion of dark matter particles isprevented in the code to avoid local deformations of thegravitational potential which will not be present in the staticcase. Effectively, in this way the distribution of dark matteris “frozen”, and live haloes are coarse grained representa-tions of static haloes.For static haloes, we convolve the analytic Hernquistpotential with a spline softened potential in the centre, suchthat gas particles feel exactly the same gravitational accel-eration as in the case of live haloes, where the central poten-tial needs to be softened to minimize effects from two-bodyencounters. For the Hernquist potential, gravitational accel-erations are given by g ( r ) = − GM vir ( r + a ) . (3)We modify the gravitational acceleration felt by gas particlesin the code, viz. g ( r ) = − GM vir (cid:0) πr h M vir (cid:1) − + ( r + a ) exp( − h /r ) , (4)where h and h are free coefficients (to first approximation h ∼ h ∼ r soft ). We determine the value of h and h forthree different resolution runs by fitting equation (4) to thegravitational acceleration of a live “frozen” halo with 10 dark matter particles simulated three times assuming r soft c (cid:13) , 000–000 Sijacki et al. ρ g a s ( r ) [ M O • k p c - ] GADGET Static Halo Live Halo GADGET Static Halo Live Halo GADGET Static Halo Live Halo GADGET Static Halo Live Halo GADGET Static Halo Live Halo GADGET Static Halo Live Halo T [ K ] r [ kpc ] S ~ T / ρ / ρ g a s ( r ) [ M O • k p c - ] AREPO Static Halo Live Halo AREPO Static Halo Live Halo AREPO Static Halo Live Halo AREPO Static Halo Live Halo AREPO Static Halo Live Halo AREPO Static Halo Live Halo AREPO Static Halo Live Halo T [ K ] r [ kpc ] S ~ T / ρ / Figure 25.
Radial profiles of gas density, temperature and entropy at time t = 0 . N DM = 10 and r soft = 14 . N DM = 10 and r soft = 6 . N DM = 10 and r soft = 3 . N gas = 10 is kept fixed. Inthe right-hand panel we also show the entropy profile obtained with the moving-mesh code in the entropy mode for our highest resolutionsimulation (continuous orange line). The vertical black lines with the same style indicate the softening scales of the dark matter potential.The black vertical dot-dashed lines denote the virial radius of the system. c (cid:13)000
Radial profiles of gas density, temperature and entropy at time t = 0 . N DM = 10 and r soft = 14 . N DM = 10 and r soft = 6 . N DM = 10 and r soft = 3 . N gas = 10 is kept fixed. Inthe right-hand panel we also show the entropy profile obtained with the moving-mesh code in the entropy mode for our highest resolutionsimulation (continuous orange line). The vertical black lines with the same style indicate the softening scales of the dark matter potential.The black vertical dot-dashed lines denote the virial radius of the system. c (cid:13)000 , 000–000 he hydrodynamics of galaxy formation values appropriate for the numerical resolutions we want toinvestigate.In Figure 24 we illustrate the outcome of this method.In the upper panel, the green lines denote the radial gravi-tational acceleration of the live “frozen” halo with 10 darkmatter particles, adopting r soft = 3, 6 . N gas = 10 (corresponding to r soft , gas = 3 kpcfor GADGET , while gravitational softenings are adaptivein
AREPO but with a floor of 3 kpc), while we increasethe number of dark matter particles from N DM = 10 ( r soft , DM = 14 kpc), to 10 ( r soft , DM = 6 . ( r soft , DM = 3 kpc). For each live halo we run a match-ing static halo simulation where the number of gas parti-cles/cells is kept the same, i.e. N gas = 10 . In all runs, gasself-gravity is neglected and there are no radiative losses.In close encounters of gas and dark matter, the effectivegravitational softening is the maximum between r soft , gas and r soft , DM , such that r soft , DM is the relevant length scale of theproblem.In Figure 25 we show radial profiles of gas density, tem-perature and entropy at time t = 0 . GADGET , while the right-hand pan-els show the results with
AREPO . Due to gravitational N-body heating there are systematic differences in the centralgas properties. The central gas density is lower, while thecentral gas temperature and entropy are higher for the livehaloes, due to spurious transfer of energy from dark matterto gas, which heats the gas and makes it expand. In general,N-body heating effects are largest for N DM = 10 and small-est for N DM = 10 , and they are confined to spatial regionswithin r soft , DM , which is reassuring. We can see by compar-ing the left-hand to the right-hand panels that AREPO isindeed much more sensitive to spurious gravitational heat-ing, as discussed in Springel (2010a), with the central en-tropy boosted by two orders of magnitude. In the right-handpanel of Figure 25, we also show radial entropy profile ob-tained with
AREPO in the dual entropy mode (where en-tropy instead of energy is explicitly conserved for cells withall Mach numbers less than 1 .
1) with N DM = 10 particles(orange continuous line). This indicates that at least partof the central entropy core is generated by the weak shockswhich gas cells experience when moving through the grainydark matter potential.At the final time t = 2 .
45 Gyr, when the system hasreached an equilibrium state, the difference between thematching live and static halo runs for r > r soft , DM is ∼ ∼ ∼ GADGET , while it is somewhat higher in
AREPO , i.e. ∼ − r < r soft the ratios of central entropy values are ≤ ≤ ≤ GADGET and ≤ ≤
150 and ≤ AREPO , for low, intermediate and high resolution runs.Even though the amount of artificial heating is significantlylarger in
AREPO than it is for
GADGET , the departuresbetween static and live halo radial profiles always occurwithin r soft , DM for both codes. Provided that gravitationalsoftening values for the dark matter component are chosencautiously, our numerical experiments hence indicate thatthe systematic discrepancy in central entropy values foundbetween SPH and mesh-based codes for the Santa Bar-bara cluster comparison project (Frenk et al. 1999; Springel2010a) is unlikely to be due to effects from a grainy darkmatter potential during cosmological structure formation. In this study we have carried out a detailed comparison be-tween the SPH code
GADGET and the new moving-meshcode
AREPO on a number of hydrodynamical test problems,which are crucial for understanding cosmological simulationsof galaxy formation. In a purely hydrodynamical regimewithout gas self-gravity or an external gravitational poten-tial we have first carried out a set of numerical experimentspreviously considered in the literature, some of which haverarely been shown for SPH codes, such as the 2D implosiontest. We have then focused on idealized non-radiative galaxycluster simulations, specifically aimed towards benchmark-ing differences in hydro solvers for problems with shocks andfluid instabilities. In simulations where radiative losses wereincluded we have analyzed the amount of baryons which coolfrom the hot halo atmospheres in
GADGET and
AREPO ,both for rotating and non-rotating haloes. In the formercase, we have also studied how the central baryonic disksform, and how orbiting substructures affect the disk mor-phology and the star formation rate. Finally, we have con-structed special test problems designed to gauge the effectof gravitational N-body heating on the gaseous properties ofthe haloes. Our main conclusions are as follows: • While post-shock fluid properties are captured wellboth in
GADGET and
AREPO in the case of 1D shock tubetests with high Mach numbers, the shocks are significantlybroader in
GADGET , and substantial post-shock oscillationsdevelop, largely because of an inadequate treatment of theinitial contact discontinuity.
AREPO in the moving-meshmode preserves the contact discontinuity much more accu-rately, but it is broadened if we employ a static mesh, dueto larger numerical diffusivity in this case. Note that themore accurate treatment of the contact discontinuity in themoving mesh code can lead to a larger ‘wall heating’ effect(Rider 2000) than in static grid codes, which tend to wash-out at some level the initial start-up errors at the contactdiscontinuity (see also description of the Noh problem inSpringel 2010a). • For shocks with complicated geometries in multi-dimensions, differences between
GADGET and
AREPO aremore striking. Even though the global fluid properties re-main similar, the sampling of the fluid properties is muchnoisier in
GADGET , and the development of dynamical fluidinstabilities is inhibited. These problems are not specific to
GADGET , but are inherent to the standard SPH method as c (cid:13) , 000–000 Sijacki et al. a whole (see also Agertz et al. 2007; Springel 2010b, and dis-cussion in Section 2.1.2). The suppression of fluid instabili-ties reduces the amount of entropy generation by mixing andartificially prolongs the lifetime of gaseous structures whichare moving through a medium with a different density. Ad-ditionally, vorticity generation in the wake of curved shocksdue to the baroclinic source term is largely suppressed in
GADGET , which directly impacts angular momentum trans-fer by vortices and the level of generated turbulence. • These fundamental differences between
GADGET and
AREPO identified in simple hydrodynamical test problemsaffect the properties of gas in more realistic, cosmologicallymotivated simulations as well. Specifically, in non-radiativeidealized simulations of merging galaxy clusters and in simu-lations of isolated haloes with orbiting gaseous substructureswe find that in
AREPO : i) an entropy core is produced in thecentre due to more efficient fluid mixing, ii) the gas strippingrate from the orbiting substructures is larger, and iii) morevorticity is produced in the wake of curved shocks. Thesefindings are in line with results of previous works which sim-ulated similar problems with AMR codes (e.g. Mitchell et al.2009; Vazza et al. 2011; Vazza 2011). Moreover, unphysicaldissipation of shocks and subsonic turbulence in GADGET ,as shown by Bauer & Springel (2012) (see also Section 4.2of Paper I), leads to the heating of the halo outskirts ratherthan of the central regions. • In radiative simulations of isolated haloes without anynet rotation, gas cooling and condensation into stars pro-ceeds in a very similar fashion in
GADGET and
AREPO ,given that the gas properties at the cooling radius matchclosely. However, for spinning haloes there is a net differ-ence in the total amount of cold baryons, which is higherin
AREPO , and this difference persists at a level of about10 −
15% in our highest resolution runs. Baryonic diskswhich form due to dissipative collapse of rotating gas aresystematically larger in
AREPO at low resolution, and theconvergence rate of the gas disk sizes is higher. • In numerical experiments where we follow the interac-tion between a forming central disk and 10 orbiting substruc-tures in an isolated halo which radiatively cools, the finaldisk morphology is significantly different. While in
AREPO an extended disk is produced with well developed spiral armsand a central bar, in
GADGET the disk is more centrally con-centrated and amorphous. Orbiting substructures are muchmore efficiently stripped of their gas content in the moving-mesh calculation and the material is incorporated into thehost halo atmosphere. Instead, in
GADGET gaseous sub-structures are more coherent, thus they lose more angu-lar momentum from hydro-/dynamical friction, and whenpassing through the central disk they induce morphologicaltransformations. • While star formation is more readily truncated in in-falling substructures due to gas stripping, extended gaseousdisks in
AREPO have significantly larger star formation ratesfor many Gyrs than is the case for
GADGET . • Due to its better ability to detect weak shocks,
AREPO is more sensitive to gravitational N-body heating (Springel2010a). This is confirmed by our specifically designed nu-merical experiments, which allow us to quantify the magni-tude of spurious gas heating for simulations of haloes witha finite number of dark matter particles with respect to theanalytic dark matter potentials. Nonetheless, the spatial ex- tent of this artificial heating is reassuringly constrained tolie within the gravitational softening length for typical set-ups, which is a scale below which simulation results are nottrustworthy due to resolution effects anyway.The numerical experiments presented in this studyclearly demonstrate that several important shortcomings ofthe SPH solver not only affect idealized test problems butare equally detrimental in more realistic setups relevant forstructure formation, and ultimately in full cosmological sim-ulations. This is especially the case because of: i) the compli-cated flows involving multiphase media which are the normin cosmological simulations, and ii) the hierarchical natureof structure formation where low mass systems are alwayspoorly resolved. In both of these regimes the hydrodynami-cal solver of the standard SPH method exhibits the largestinaccuracies. On the other hand, our numerical tests confirmand significantly extend the findings of Springel (2010a) thatthe moving-mesh code AREPO delivers a physically more ac-curate representation of the evolution of inviscid gases. Notethat in the current work we deliberately kept the modellingof the baryon physics at a very simple level, so as to isolatethe differences between the hydro solvers in as clean a man-ner as possible. While inclusion of more realistic feedbackmechanisms, such as supernovae winds and AGN heating,will likely modify the properties of the simulated galaxiessignificantly, it is of prime importance to disentangle nu-merical inaccuracies of the hydro solver from the uncertain-ties of the feedback physics modelling. It is hence clear thatcosmological simulations with
AREPO have the potential toprovide a much more realistic description of structure for-mation in the Universe (see also Papers I and II), somethingwe will explore in more depth in forthcoming studies.
ACKNOWLEDGEMENTS
We would like to thank Andrew MacFadyen and DanielEisenstein for very useful discussions and suggestions onthe topic, and Manfred Kitzbichler and Diego Mu˜noz forcarefully reading the manuscript. We would like to thankthe anonymous referee for many constructive suggestions,which helped to improve the presentation of the results.DS acknowledges NASA Hubble Fellowship through grantHST-HF-51282.01-A. DK acknowledges NASA Hubble Fel-lowship through grant HST-HF-51276.01-A. The computa-tions in this paper were performed on the Odyssey clustersupported by the FAS Science Division Research ComputingGroup at Harvard University.
REFERENCES
Abel T., 2011, MNRAS, 413, 271Agertz O., Moore B., Stadel J., Potter D., Miniati F., ReadJ., Mayer L., Gawryszczak A., Kravtsov A., Nordlund˚A., Pearce F., Quilis V., Rudd D., Springel V., Stone J.,Tasker E., Teyssier R., Wadsley J., Walder R., 2007, MN-RAS, 380, 963Agertz O., Teyssier R., Moore B., 2011, MNRAS, 410, 1391Balsara D. S., 1995, Journal of Computational Physics, 121,357Barnes J., Hut P., 1986, Nature, 324, 446 c (cid:13)000
Abel T., 2011, MNRAS, 413, 271Agertz O., Moore B., Stadel J., Potter D., Miniati F., ReadJ., Mayer L., Gawryszczak A., Kravtsov A., Nordlund˚A., Pearce F., Quilis V., Rudd D., Springel V., Stone J.,Tasker E., Teyssier R., Wadsley J., Walder R., 2007, MN-RAS, 380, 963Agertz O., Teyssier R., Moore B., 2011, MNRAS, 410, 1391Balsara D. S., 1995, Journal of Computational Physics, 121,357Barnes J., Hut P., 1986, Nature, 324, 446 c (cid:13)000 , 000–000 he hydrodynamics of galaxy formation Bauer A., Springel V., 2012, MNRAS, p. 3102Bertschinger E., 1989, ApJ, 340, 666Borgani S., Dolag K., Murante G., Cheng L.-M., SpringelV., Diaferio A., Moscardini L., Tormen G., Tornatore L.,Tozzi P., 2006, MNRAS, 367, 1641Cha S.-H., Whitworth A. P., 2003, MNRAS, 340, 73Crain R. A., Theuns T., Dalla Vecchia C., Eke V. R., FrenkC. S., Jenkins A., Kay S. T., Peacock J. A., Pearce F. R.,Schaye J., Springel V., Thomas P. A., White S. D. M.,Wiersma R. P. C., 2009, MNRAS, 399, 1773De Villiers J.-P., Hawley J. F., 2003, ApJ, 589, 458Di Matteo T., Khandai N., DeGraf C., Feng Y., CroftR. A. C., Lopez J., Springel V., 2012, ApJ, 745, L29Dolag K., Borgani S., Murante G., Springel V., 2009, MN-RAS, 399, 497Dolag K., Vazza F., Brunetti G., Tormen G., 2005, MN-RAS, 364, 753Evrard A. E., 1988, MNRAS, 235, 911Frenk C. S., White S. D. M., Bode P., Bond J. R., BryanG. L., Cen R., Couchman H. M. P., Evrard A. E. e. a.,1999, ApJ, 525, 554Gammie C. F., McKinney J. C., T´oth G., 2003, ApJ, 589,444Hayward C., et al., 2012, in preparationHeitmann K., Luki´c Z., Fasel P., Habib S., Warren M. S.,White M., Ahrens J., Ankeny L., Armstrong R., O’SheaB., Ricker P. M., Springel V., Stadel J., Trac H., 2008,Computational Science and Discovery, 1, 015003Hernquist L., 1987, ApJS, 64, 715Hernquist L., 1990, ApJ, 356, 359Hernquist L., Katz N., 1989, ApJS, 70, 419Hernquist L., Springel V., 2003, MNRAS, 341, 1253Heß S., Springel 2011, MNRAS, submittedHeß S., Springel V., 2010, MNRAS, 406, 2289Hui W. H., Li P. Y., Li Z. W., 1999, Journal of Computa-tional Physics, 153, 596Iapichino L., Schmidt W., Niemeyer J. C., Merklein J.,2011, MNRAS, 414, 2297Inutsuka S.-I., 2002, Journal of Computational Physics,179, 238Janka H.-T., Langanke K., Marek A., Mart´ınez-Pinedo G.,M¨uller B., 2007, Phys. Rep., 442, 38Katz N., Weinberg D. H., Hernquist L., 1996, ApJS, 105,19Kay S. T., Pearce F. R., Jenkins A., Frenk C. S., WhiteS. D. M., Thomas P. A., Couchman H. M. P., 2000, MN-RAS, 316, 374Keres D., et al., 2011, MNRAS, submitted, astro-ph:1109.4638Klessen R. S., Krumholz M. R., Heitsch F., 2009, SpecialIssue on Computational Astrophysics, Eds. Lucio Mayer,Advanced Science Letters, 4, 258-285 (2011)Marri S., White S. D. M., 2003, MNRAS, 345, 561McCarthy I. G., Bower R. G., Balogh M. L., Voit G. M.,Pearce F. R., Theuns T., Babul A., Lacey C. G., FrenkC. S., 2007, MNRAS, 376, 497Mitchell N. L., McCarthy I. G., Bower R. G., Theuns T.,Crain R. A., 2009, MNRAS, 395, 180Morris J. P., Monaghan J. J., 1997, Journal of Computa-tional Physics, 136, 41Murante G., Borgani S., Brunino R., Cha S.-H., 2011, MN-RAS, accepted, astro-ph: 1105.1344 Okamoto T., Jenkins A., Eke V. R., Quilis V., Frenk C. S.,2003, MNRAS, 345, 429O’Shea B. W., Nagamine K., Springel V., Hernquist L.,Norman M. L., 2005, ApJS, 160, 1Pretorius F., 2007, ”Relativistic Objects in Compact Bi-naries: From Birth to Coalescence”, Editor: Colpi et al.,Pulisher: Springer Verlag, Canopus Publishing Limited,eprint arXiv:0710.1338Price D. J., 2008, Journal of Computational Physics, 2271,10040Puchwein E., Springel V., Sijacki D., Dolag K., 2010, MN-RAS, 406, 936Rasio F. A., 2000, Progress of Theoretical Physics Supple-ment, 138, 609Read J. I., Hayfield T., Agertz O., 2010, MNRAS, 405, 1513Ricker P. M., Sarazin C. L., 2001, ApJ, 561, 621Rider W. J., 2000, Journal of Computational Physics, 162,395Ritchie B. W., Thomas P. A., 2001, MNRAS, 323, 743Ritchie B. W., Thomas P. A., 2002, MNRAS, 329, 675Robinson M., Monaghan J. J., 2011, arXiv:1101.2240Saitoh T. R., Makino J., 2012, ApJ submitted, arXiv:1202.4277Sijacki D., Springel V., Haehnelt M. G., 2009, MNRAS,400, 100Springel V., 2005, MNRAS, 364, 1105Springel V., 2010a, MNRAS, 401, 791Springel V., 2010b, ARA&A, 48, 391Springel V., 2011, Invited review for the volume ”Tessel-lations in the Sciences: Virtues, Techniques and Applica-tions of Geometric Tilings”, eds. R. van de Weijgaert, G.Vegter, J. Ritzerveld and V. Icke, astro-ph: 1109.2218Springel V., Farrar G. R., 2007, MNRAS, 380, 911Springel V., Hernquist L., 2002, MNRAS, 333, 649Springel V., Hernquist L., 2003a, MNRAS, 339, 289Springel V., Hernquist L., 2003b, MNRAS, 339, 312Springel V., White S. D. M., Jenkins A., Frenk C. S.,Yoshida N., Gao L., Navarro J., Thacker R., Croton D.,Helly J., Peacock J. A., Cole S., Thomas P., CouchmanH., Evrard A., Colberg J., Pearce F., 2005, Nature, 435,629Steinmetz M., White S. D. M., 1997, MNRAS, 288, 545Vazza F., 2011, MNRAS, 410, 461Vazza F., Brunetti G., Gheller C., Brunino R., 2010, NewAstronomy, 15, 695Vazza F., Brunetti G., Gheller C., Brunino R., Br¨uggenM., 2011, A&A, 529, A17Vazza F., Dolag K., Ryu D., Brunetti G., Gheller C., KangH., Pfrommer C., 2011, MNRAS, submitted, astro-ph:1106.2159Vogelsberger M., Sijacki D., Keres D., Springel V., Hern-quist L., 2011, MNRAS, submitted, astro-ph: 1109.1281Wadsley J. W., Veeravalli G., Couchman H. M. P., 2008,MNRAS, 387, 427White S. D. M., Frenk C. S., 1991, ApJ, 379, 52ZuHone J. A., 2011, ApJ, 728, 54 c (cid:13) , 000–000 Sijacki et al. ρ g a s ( r ) [ M O • k p c - ] T [ K ] r [ kpc ] σ g a s , D [ k m / s ] Figure A1.
Radial profiles of gas density, temperature and1D velocity dispersion at time t = 0 .
05 Gyr (dotted lines) and t = 2 .
45 Gyr (continuous lines) for GADGET (blue lines) andAREPO (red lines) with N gas = 10 in a static dark matter halowith a Hernquist profile. Vertical black lines indicate the soften-ing scales of the gas equal to r soft = 3 . t [ Gyr ] E k i n / E t o t , | E po t | / E t o t , E i n t / E t o t Figure A2.
Time evolution of the ratio of gas kinetic (dottedlines), potential (dashed lines), and internal (continuous lines) tothe total gas energy in GADGET (blue) and AREPO (red) sim-ulations of isolated halos in hydrostatic equilibrium. Gas is self-gravitating and represented by N gas = 10 resolution elementsand is evolved within a static Hernquist dark matter halo. APPENDIX A: ISOLATED HALOES INHYDROSTATIC EQUILIBRIUM
Here we show how gaseous spheres placed in hydrostaticequilibrium within a static Hernquist dark matter haloevolve with time to determine the level of accuracy of ourinitial conditions. The general set-up is described in detailin Section 3.3.1. Specifically for the figures presented herewe assume that the gas is self-gravitating and that it hasinitially zero velocity.In Figure A1, we show gas density, temperature and 1Dvelocity dispersion radial profiles at t = 0 .
05 Gyr (dottedlines) and t = 2 .
45 Gyr (continuous lines) for
GADGET (bluelines) and
AREPO (red lines) with N gas = 10 resolution el-ements. As discussed in Section 3.3.1, due to the Poissonsampling of the gas positions in the initial conditions theyare not perfectly relaxed, which leads to the developmentof small-scale random motions, as evidenced by a non-zerogas velocity dispersion. Over time these residual gas veloci-ties are dissipated, as can be seen from the bottom panel ofFigure A1. This leads to a slight readjustment of the tem-perature distribution, which is very similar in both codes.Note, however, that the kinetic energy of the random gasmotions is at all times a very small fraction of both the po-tential and the internal energy, as illustrated in Figure A2,amounting to only about 0 .
2% of the total gas energy after2 .
45 Gyr. c (cid:13)000