The AGORA High-resolution Galaxy Simulations Comparison Project. II. Isolated Disk Test

28 Pages·2016·13.32 MB·English

DRAFTVERSIONNOVEMBER1,2016
PreprinttypesetusingLATEXstyleemulateapjv.12/16/11
THEAGORAHIGH-RESOLUTIONGALAXYSIMULATIONSCOMPARISONPROJECT.II:ISOLATEDDISKTEST
JI-HOONKIM1,2,3,4,OSCARAGERTZ5,6,ROMAINTEYSSIER7,MICHAELJ.BUTLER8,†,‡,DANIELCEVERINO9,†,‡,JUN-HWAN
CHOI10,†,‡,ROBERTFELDMANN7,11,†,BENW.KELLER12,†,‡,ALESSANDROLUPI13,†,‡,THOMASQUINN14,†,‡,YVESREVAZ15,†,‡,
SPENCERWALLACE14,†,NICKOLAYY.GNEDIN16,17,18,‡,SAMUELN.LEITNER19,‡,SIJINGSHEN20,‡,BRITTOND.SMITH21,‡,ROBERT
THOMPSON22,‡,MATTHEWJ.TURK23,‡,TOMABEL1,2,KENZAS.ARRAKI24,SAMANTHAM.BENINCASA12,SUKANYA
CHAKRABARTI25,COLINDEGRAF20,AVISHAIDEKEL26,NATHANJ.GOLDBAUM22,PHILIPF.HOPKINS3,CAMERONB.HUMMELS3,
6 ANATOLYKLYPIN24,HUILI27,PIEROMADAU28,13,NIRMANDELKER29,26,LUCIOMAYER7,KENTARONAGAMINE30,31,SARAH
1 NICKERSON7,BRIANW.O’SHEA32,JOELR.PRIMACK33,SANTIROCA-FA`BREGA26,VADIMSEMENOV17,IKKOHSHIMIZU30,
0 CHRISTINEM.SIMPSON34,KEITATODOROKI35,JAMESW.WADSLEY12,ANDJOHNH.WISE36FORTHEAGORACOLLABORATION37
2 1KavliInstituteforParticleAstrophysicsandCosmology,SLACNationalAcceleratorLaboratory,MenloPark,CA94025,USA
2DepartmentofPhysics,StanfordUniversity,Stanford,CA94305,USA
t 3DepartmentofAstronomy,CaliforniaInstituteofTechnology,Pasadena,CA91125,USA
c
O 4EinsteinFellow,[email protected]
5DepartmentofPhysics,UniversityofSurrey,Guildford,Surrey,GU27XH,UnitedKingdom
8 6LundObservatory,DepartmentofAstronomyandTheoreticalPhysics,LundUniversity,SE-22100Lund,Sweden
7CentreforTheoreticalAstrophysicsandCosmology,InstituteforComputationalScience,UniversityofZurich,Zurich,8057,Switzerland
2
8Max-Planck-Institutfu¨rAstronomie,D-69117Heidelberg,Germany

9Zentrumfu¨rAstronomiederUniversita¨tHeidelberg,Institutfu¨rTheoretischeAstrophysik,69120Heidelberg,Germany
]
A 10DepartmentofAstronomy,UniversityofTexas,Austin,TX78712,USA
11DepartmentofAstronomy,UniversityofCaliforniaatBerkeley,Berkeley,CA94720,USA
G 12DepartmentofPhysicsandAstronomy,McMasterUniversity,Hamilton,ONL8S4M1,Canada
. 13Institutd’AstrophysiquedeParis,SorbonneUniversites,UPMCUnivParis6etCNRS,75014Paris,France
h 14DepartmentofAstronomy,UniversityofWashington,Seattle,WA98195,USA
p 15InstituteofPhysics,Laboratoired’Astrophysique,E´colePolytechniqueFe´de´raledeLausanne,CH-1015Lausanne,Switzerland
- 16ParticleAstrophysicsCenter,FermiNationalAcceleratorLaboratory,Batavia,IL60510,USA
o
17DepartmentofAstronomyandAstrophysics,UniversityofChicago,Chicago,IL60637,USA
r
t 18KavliInstituteforCosmologicalPhysics,UniversityofChicago,Chicago,IL60637,USA
s 19DepartmentofAstronomy,UniversityofMaryland,CollegePark,MD20742,USA
a
20KavliInstituteforCosmology,UniversityofCambridge,Cambridge,CB30HA,UnitedKingdom
[
21InstituteforAstronomy,UniversityofEdinburgh,RoyalObservatory,Edinburgh,EH93HJ,UnitedKingdom
2 22NationalCenterforSupercomputingApplications,UniversityofIllinois,Urbana,IL61801,USA
v 23SchoolofInformationSciences,DepartmentofAstronomy,UniversityofIllinois,Urbana,IL61801,USA
6 24DepartmentofAstronomy,NewMexicoStateUniversity,LasCruces,NM88001,USA
6 25SchoolofPhysicsandAstronomy,RochesterInstituteofTechnology,Rochester,NY14623,USA
0 26CenterforAstrophysicsandPlanetaryScience,RacahInstituteofPhysics,TheHebrewUniversity,Jerusalem,91904,Israel
27DepartmentofAstronomy,UniversityofMichigan,AnnArbor,MI48109,USA
3
28DepartmentofAstronomyandAstrophysics,UniversityofCaliforniaatSantaCruz,SantaCruz,CA95064,USA
0
29DepartmentofAstronomy,YaleUniversity,NewHaven,CT06520,USA
.
0 30DepartmentofEarthandSpaceScience,GraduateSchoolofScience,OsakaUniversity,Toyonaka,Osaka,560-0043,Japan
1 31DepartmentofPhysicsandAstronomy,UniversityofNevada,LasVegas,NV89154,USA
6 32DepartmentofComputationalMathematics,ScienceandEngineering,DepartmentofPhysicsandAstronomy,NationalSuperconductingCyclotron
1 Laboratory,MichiganStateUniversity,Lansing,MI48824,USA
: 33DepartmentofPhysics,UniversityofCaliforniaatSantaCruz,SantaCruz,CA95064,USA
v 34HeidelbergerInstitutfu¨rTheoretischeStudien,69118Heidelberg,Germany
i 35DepartmentofPhysicsandAstronomy,UniversityofKansas,Lawrence,KS66045,USA
X
36CenterforRelativisticAstrophysics,SchoolofPhysics,GeorgiaInstituteofTechnology,Atlanta,GA30332,USA
r 37http://www.AGORAsimulations.org/
a †Theseauthors,inalphabeticalorder,contributedtothearticlebyleadingtheeffortwithineachcodegroupstoperformandanalyzesimulations.and
‡Theseauthors,inalphabeticalorder,contributedtothearticlebydevelopingGRACKLEandimplementingitsinterfacewithparticipatingcodes.
DraftversionNovember1,2016
ABSTRACT
Using an isolated Milky Way-mass galaxy simulation, we compare results from 9 state-of-the-art gravito-
hydrodynamicscodes widely used in the numerical community. We utilize the infrastructure we have built
for the AGORA High-resolution Galaxy Simulations Comparison Project. This includes the common disk
initial conditions, common physics models (e.g., radiative cooling and UV background by the standardized
package GRACKLE) and common analysis toolkit yt, all of which are publicly available. Subgrid physics
models such as Jeans pressure floor, star formation, supernova feedback energy, and metal production are
carefullyconstrainedacrosscodeplatforms. With numericalaccuracythatresolvesthedisk scale height,we
findthatthecodesoverallagreewellwithoneanotherinmanydimensionsincluding: gasandstellarsurface
densities, rotation curves, velocity dispersions, density and temperature distribution functions, disk vertical
heights, stellar clumps, star formation rates, and Kennicutt-Schmidt relations. Quantities such as velocity
dispersions are very robust (agreementwithin a few tens of percent at all radii) while measures like newly-
2 J. KIMETAL. FORTHEAGORACOLLABORATION
formed stellar clump mass functions show more significant variation (difference by up to a factor of ∼3).
Systematic differences exist, for example, between mesh-based and particle-based codes in the low density
region,andbetweenmorediffusiveandlessdiffusiveschemesinthehighdensitytailofthedensitydistribution.
Yetintrinsiccodedifferencesaregenerallysmallcomparedtothevariationsinnumericalimplementationsof
thecommonsubgridphysicssuchassupernovafeedback.Ourexperimentreassuresthat,ifadequatelydesigned
in accordance with our proposedcommon parameters, results of a modernhigh-resolutiongalaxyformation
simulationaremoresensitivetoinputphysicsthantointrinsicdifferencesinnumericalschemes.
Keywords:cosmology:theory–galaxies:formation–galaxies:evolution–galaxies:kinematicsanddynamics
–ISM:structure–methods:numerical–hydrodynamics
1. INTRODUCTION Through a multi-platform approach from the beginning, we
strive to improve all our codes by “increasing the level of
Decadesof strenuouseffortsbycomputationalastrophysi-
realism and predictive power of galaxy simulations and the
cistshavepropellednumericalexperimentstobecomeoneof
understandingofthefeedbackprocessesthatregulategalaxy
the most widely used tools in theorizing how galaxies form
metabolism,”(Kimetal.2014)andbydoingsotofindsolu-
in the Universe. Numerical experiments are often the only
tionstolong-standingproblemsingalaxyformation.Because
meanstoputourtheoryto atest, theresultofwhichwe can
theinterplaybetweennumericalresolutionandsubgridmod-
comparewithobservationaldatato validatethemodel’sfea-
elingsofstellarphysicsiscrucialingalaxy-scalesimulations,
sibility. Sincethesuccessofgalaxyformationtheoryispred-
we require thatsimulationsbe designedwith state-of-the-art
icatedonrobustnumericalexperiments,itisonlyreasonable
resolution,.100pc,whichiscurrentlyallowedwithinreal-
that we apply the same scientific standard of reproducibility
isticcomputationalcostbounds.
togalaxyformationsimulations. Inotherwords,itshouldbe
In the Project’s flagship paper, Kim et al. (2014), we ex-
considered as a fundamentalprinciple that researchers must
plained the philosophy behind the Project and detailed the
notestablish findingsfroma singlenumericalexperimentas
publiclyavailableProjectinfrastructurewehaveputtogether.
scientific knowledge. Only after the result is reproducedin-
We also described the proof-of-concept test, in which we
dependentlybyotherresearchersandproventobenotaniso-
field-testedourinfrastructurewithadarkmatter-onlycosmo-
latedincidencecanwebuildanyconclusivetheoryabouthow
logicalzoom-insimulation finding a robust convergencebe-
galaxiesactuallyformintheUniverse.
tween participating codes. More than 140 researchers from
However, the task of replicating galaxy simulations, or
over60academicinstitutionsworldwidehavesinceagreedto
equivalently, comparing simulations between codes, has not
take part in the Collaboration, many of whom having been
receivedhigh priority.1 Instead, the task is consideredcom-
actively engaged in working groups and sub-projects.2 The
plex and time-consuming because one needs to ensure that
cohort of numerical codes participating in the Project cur-
identicalphysicsisusedinanidenticalinitialconditionwith
rently include, but are not limited to in future studies: the
identical runtime settings. It is sometimes perceived as te-
Lagrangian smoothed particle hydrodynamics codes (SPH;
dious and unrewarding for early-career researchers. In fact,
Gingold & Monaghan 1977; Lucy 1977; Monaghan 1992)
thelackofreproducibilitychecksisnotuniquetothefieldof
numerical galaxy formation (e.g., Open Science Collabora- CHANGA, GADGET, GASOLINE, and GEAR, and the Eule-
rianadaptivemeshrefinementcodes(AMR;Berger&Oliger
tion2015;NatureSurvey2016). Anditscauseisnotsimply
anunwillingnessofonlycomputationalastrophysicists,either 1984; Berger & Colella 1989) ART-I, ART-II, ENZO, and
(Everett&Earp2015). Rather,addressingthesystem(orthe RAMSES, and the mesh-free finite-volume Godunov code
lackthereof)whichchecksthereproducibilityofsimulations GIZMO(seeSection5forinformationoneachcode).
In this second report of our continuing endeavor, we use
wouldrequireacollectiveactionbytheentirecommunity. It
an isolated Milky Way-mass galaxy simulation to compare
cannot be simply about asking individual researchers to re-
9 widely used state-of-the-artgravito-hydrodynamicscodes.
leasetheirdatadumps,butitshouldbeaboutbuildingasys-
As in all comparison studies in AGORA, the participating
temthatincentivizessimulationspublishedinareproducible
codes share the common initial condition (i.e., generated by
manner. Itshouldalsobeaboutassemblinganinfrastructure
thatreducesthecostofreproducibilitychecks,onwhichsim- MAKEDISK;seeSection2),commonphysicsmodels(e.g.,ra-
diativecoolingandUVbackgroundprovidedbythestandard-
ulations are verified routinely and effortlessly (Nosek et al.
2015;Begley&Ioannidis2015). izedpackageGRACKLE;seeSection3.1;Smith&etal.2016;
TheAGORAHigh-resolutionGalaxySimulationsCompar- Bryan et al. 2014; Kim et al. 2014),3 and common analysis
ison Project (Assembling Galaxies Of Resolved Anatomy) is platform(i.e,yttoolkit;Turketal.2011).4 Weadoptspatial
the collective response by the numerical galaxy formation resolutionof80pc thatresolvesthescale heightofthe disk.
communitytosuchachallenge.Sinceitsfirstmeetingin2012 This helps the codes to be less dependent on phenomeno-
at the University of California at Santa Cruz, AGORA Col- logical prescriptions of sub-resolution processes which are
laborationhas aimed to comparegalaxy-scalenumericalex- inevitablyintroducedin low-resolution(> kpc) simulations.
perimentsonavarietyofcodeplatformswithstate-of-the-art Asmoderngalaxyformationsimulationswithstate-of-the-art
resolution.Oursharedgoalistoensurethatphysicalassump- resolutionand physics prescriptionsbecome more and more
tionsareresponsibleforanysuccessingalaxyformationsim- computationallyexpensive,itistimelythatwecomparehigh-
ulations, rather than artifacts of particular implementations.
2SeetheProjectwebsiteathttp://www.AGORAsimulations.org/for
1Codecomparisonsintheastrophysicalcommunityhavepreviouslybeen more information on the Project including its membership, and its task-
undertaken, albeit with simplified physics in a different scale (e.g., Frenk orientedandscience-orientedworkinggroups.
etal.1999;O’Sheaetal.2005),orfocusingonlyonhydrodynamicssolvers 3Thewebsiteishttp://grackle.readthedocs.org/.
(e.g.,Agertzetal.2007;Taskeretal.2008). 4Thewebsiteishttp://yt-project.org/.
AGORACOMPARISON. II:ISOLATEDDISK 3
Table1
InitialConditionCharacteristics
Darkmatterhalo Stellardisk Gasdisk Stellarbulge
Densityprofile Navarroetal.(1997) Exponential Exponential Hernquist(1990)
Structuralproperties M200=1.074×1012M⊙,vc,200=150kms−1, Md,⋆=3.438×1010M⊙, Md,gas=8.593×109M⊙, Mb,⋆=4.297×109M⊙,
R200=205.5kpc,c=10,l =0.04 rd=3.432kpc,zd=0.1rd fgas=0.2 Mb,⋆/Md=0.1
Numberofparticles 105 105 105 1.25×104
Particlemass mDM=1.254×107M⊙ m⋆,IC=3.437×105M⊙ mgas,IC=8.593×104M⊙ m⋆,IC=3.437×105M⊙
resolutionisolateddisksimulationstocheckhowsuccessfully profile(M /M =0.1;Hernquist1990). The disk or bulge
b,⋆ d
thesegalaxiesarereproducedbytheirpeers.5 Readersshould starsintheICdonotcontributetothefeedbackbudget.
notethatourintentionisnottoidentifya“correct”or“incor- AmongthethreeresolutionchoicesprovidedinKimetal.
rect” code, butto focusinstead onjuxtaposingthe codesfor (2014)hereweemploya“low-resolution”ICwhichhas105
physicalinsightsandlearn howmuchscatter one shouldex- particleseach forthe halo, the stellar disk, andthe gasdisk,
pectamongmodernnumericaltoolsin the field (see Section and1.25×104particlesforthebulge.Theinitialgastemper-
7formorediscussion). ature in the disk is set to 104 K, not to the specific internal
The remainder of this article is organized as follows. In energy computed by MAKEDISK. The initial metal fraction
Section2weexplaintheisolateddiskinitialconditionusedin in the gas disk is 0.02041.9 When the gasdisk is initialized
thestudy.Thecommoninputphysicsandruntimeparameters
onmesh-basedcodes(ART-I,ART-II,ENZO,andRAMSES),
required in the participating codes are discussed in Section
insteadofusingtheparticlesprovidedbyMAKEDISK,were-
3 and 4, respectively. Then Section 5 describes 9 hydrody-
quirethattheparticipantsuseananalyticdensityprofileof
namicscodesthatparticipatedin thiscomparison. Section6
presentstheresultsofourcomparisonfocusingonsimilarities r d,gas(r,z)=r 0e−r/rd·e−|z|/zd (1)
and discrepancies discovered in various multi-dimensional
analyses.FinallyinSection7wesummarizeourfindingsand with r =M /(4p r2z ), where r is the cylindricalradius
0 d,gas d d
conclude the paper with remarks on future work. We will andzistheverticalheightfromthediskplane.Tosetupadisk
also stress the importance of collaborative and reproducible in a centrifugalequilibriumwe also ask thatthe participants
research in the numerical galaxy formation community the utilizetherotationalvelocityprofilebinnedfromanactualgas
AGORAProjectstrivestopromote. particle distribution within |z|<z .10 In mesh-based codes,
d
weadditionallyincludeauniformlylowdensitygashalowith
2. INITIALCONDITION n =10−6cm−3 and zero initial velocity, since they cannot
H
In this section we describe the Milky Way-mass isolated have cells with zero density. The halo is initially set to 106
initial condition (IC) we adopt in this study. While this IC K and zero metallicity. Note thatthis gaseoushalo doesnot
ispartofasetofdiskICsgeneratedforAGORAsimulations existinparticle-basedcodes(SPHcodesorGIZMO).
that were first introduced in Section 2.2 of the Project flag- It is worth to note one point about the common disk IC
shippaper(Kimetal.2014),webrieflyexplainitsimportant adoptedhere. Asreadersmayfindin Section6.6,in ourex-
structuralpropertiesforcompleteness.6 periment that includes radiative cooling, star formation and
ThediskgalaxyICwithpropertiescharacteristicofMilky feedback (“Sim-SFF”; see Section 3), only ∼109M addi-
⊙
Way-mass galaxiesat redshift z∼1 is generated with a pri- tionalstars form in 500 Myr in all codeson average. When
vatelysharedversionofMAKEDISK(Springeletal.2005).7,8 comparedwiththestellar andgascomponentspresentinthe
The IC has the following components(see also Table 1 and IC, this means only a ∼3% increase in stellar mass, and a
Figure 1): (1) a dark matter halo with M200 = 1.074× ∼12%decreaseingasmass. AsthediscussioninSection6
1012M , R =205.5 kpc and circular velocity of v = shouldmakeclear,giventhisrelativelysmallchangeinstellar
⊙ 200 c,200
150 kms−1 that follows the Navarro-Frenk-White (NFW; mass,itisexpectedthatthestellarfeedbackinourexperiment
Navarro et al. 1997) profile with concentration parameter willbeveryinefficient.
c=10andspinparameterl =0.04,(2)anexponentialdisk
withM =4.297×1010M ,scalelengthr =3.432kpcand
d ⊙ d
s2c0a%leghaesiginhtmzdas=s (0i..e1.,rdftha=t iMs com/pMosed=o0f.28)0,%(3)staarssteallnadr but9toT1h.i5s7f6ra1ctZio⊙nfaolrvGalRuAeC0K.0L2E0v421.1coorrreasbpoovned.sIttiosb1ecZa⊙usfeorthGeRsoAlaCrKmLeEtavl2li.c0-,
gas d,gas d ityunit Z⊙wasupdatedfrom0.02041to0.01295inGRACKLEv2.1. Since
bulgewithMb,⋆=4.297×109M⊙thatfollowstheHernquist cooling rates pre-tabulated by CLOUDYare at1 Z⊙, not at specific metal
fraction value, the cooling rates in GRACKLE’sequilibrium cooling mode
willdifferdependingonwhichGRACKLEversionisadopted(seeSection3.1
5Comparisonsofcosmologicalzoom-insimulationsarealsointhemaking formoreon GRACKLE). Forexample, inthecurrent study, thecodes us-
totesttherobustnessofthecodesuiteover13.8Gyrofevolution. Seethe ingGRACKLEv2.1(CHANGA,GASOLINE,GADGET-3,andGIZMO)show
Project’sflagshippaper,(Kimetal.2014),formoreinformation. slightly enhanced cooling rates thantheones using GRACKLEv2.0orbe-
6ThepublicDropboxlinkishttp://goo.gl/8JzbIJ. low(ART-I,ART-II,ENZO,RAMSES,andGEAR). Generallyspeaking,ini-
7MAKEDISKisanearlierrealizationofacodesimilartoGALIC(Yurin& tialgasmetallicity shouldbesetupsothatitisconsistentwiththechosen
Springel2014).GALICispubliclyavailable,anditswebsiteishttp://www.h- GRACKLEversioninterfacingwiththecode. Wereferinterestedreadersto
its.org/tap-software-en/galic-code/. theGRACKLEv2.1releasenoteathttps://goo.gl/BNRfwJ.
8WhileMilkyWay’s fgasis∼10%,typicalgalaxieswiththeMilkyWay 10 This actual initial velocity profile, provided in vcirc SPH.dat in
stellarmassatz∼0have fgas∼20%(Catinellaetal.2010). Inthisregard, ourpublicDropboxlink,isdifferentfromthefilevcirc.datproducedby
onecansaythatwemodelamoretypicalgalaxyatz∼0thanMilkyWay. MAKEDISKitself.Thedifferenceis∼5%inthecentralfewkpc.
4 J. KIMETAL. FORTHEAGORACOLLABORATION
Figure1. The0MyrsnapshotsoftheisolatedMilkyWay-massgalaxysimulationsby9participatingcodes.Diskgassurfacedensitiesina30kpcbox,edge-on
(top)andface-on(bottom),producedwiththecommonanalysistoolkityt. Forvisualizationsoftheparticle-basedcodeshereafter(Figures1-3,14-15,32,34,
35)–butnotinanyotheranalysesexceptthesefigures–ytusesanin-memoryoctreeonwhichgasparticlesaredepositedusingsmoothingkernels.Comparing
0Myrsnapshots–dumpedimmediatelyaftereachcodereadsintheIC–istochecktheexactidentityofICsinterpreted byeachcode. SeeSection6.1for
moreinformationonthisfigure,andSection5fordescriptionsofparticipatingcodesinthiscomparison. Thefullcolorversionofthisfigureisavailableinthe
electronicedition.Thehigh-resolutionversionsofthisfigureandarticleareavailableattheProjectwebsite,http://www.AGORAsimulations.org/.
3. COMMONPHYSICS:SIM-NOSFANDSIM-SFF codes)orthefinestcell(mesh-basedcodes)againstunphysi-
calcollapseandtoavoidartificialfragmentationduetounre-
Wenowdescribethecommonphysicsemployedinourex-
solved pressure gradient(Trueloveet al. 1997;Robertson &
perimentincludingequilibriumgascooling,metagalacticUV
Kravtsov2008). In practice, it is achievedbyenforcingthat
background, star formation, and energy and metal yields by
thelocalJeanslengthl besufficientlyresolvedwiththe
supernovae. Notethatthecommonphysicsadoptedhereisa Jeans
finestresolutionelementsatalltimes. Thatis,
variationof the commonphysicsmodelrecommendedin all
AGORA simulations by default; see Section 3 of the Project l =N D x (2)
flagshippaper(Kimetal.2014).Forthepresentstudyallpar- Jeans Jeans
ticipatingcodegroupsareaskedtoruntwosimulationsstart- whereD x=80pcistheadoptedspatialresolution(finestcell
ingfromtheidenticalIC:(1)“Sim-noSF”withradiativegas sizeorsofteninglength;seeSection4.1)andN =4isthe
Jeans
coolingbutwithoutstarformationorfeedback,and(2)“Sim- JeansnumberadoptedfromTrueloveetal.(1997).Thisgives
SFF”withradiativecooling,starformationandfeedback. In therequiredpressurefloorvalueas
Section 3.1, we first list the gas physicsthat are commonin
bothSim-noSFandSim-SFF.TheninSection3.2,thesubgrid 1
P = N2 Gr 2 D x2 (3)
prescriptionsofstellarphysicsforSim-SFFareexplained. Jeans gp Jeans gas
3.1. GasPhysics: RadiativeCooling,UVBackground,and where G is the gravitational constant, g =5/3 is the adia-
PressureFloor batic index, and r gas is the gas density. Note that NJeans is
not necessarily equal to the parameter controlling the pres-
Therateatwhichthegasinourgalaxyradiativelycoolsis
sure support in each code. For actual parameter choices for
determined by AGORA’s standard chemistry and cooling li-
selected codes, see Appendix A. For implementations using
braryGRACKLE(Smith&etal.2016;Bryanetal.2014;Kim
polytropesinART-II andRAMSES,seeSections5.2and5.4,
etal. 2014).3 Forthisstudy,theequilibriumcoolingversion
respectively.
ofGRACKLEisinterfacedwitheachparticipatingcode,either
viaGRACKLE’soriginalinterfaceorviaN.Gnedin’sauxiliary
API.11 In the chosen equilibrium cooling mode, GRACKLE 3.2. StellarPhysics: StarFormation,andEnergy,Mass,
MetalYieldsFromCore-CollapseSupernovae
follows tabulated cooling rates pre-computed by the pho-
toionizationcode CLOUDY (Ferlandetal.2013).12 Thepre- Inadditiontothegasphysicsdescribedintheprevioussec-
computedlook-uptablealsoincludesmetalcoolingratesfor tion, Sim-SFF incorporates subgrid models for star forma-
solarabundances,1 Z , asafunctionofgasnumberdensity tionandsupernovafeedback. First,aparcelofgasabovethe
⊙
and temperature. These metal cooling rates are then scaled thresholdnH,thres=10cm−3=r gas,thres/mHproducesstarsat
linearlywithmetallicitywhichisfollowedinoursimulations aratethatfollowsthelocalSchmidtlawas
as a separate passive scalar.13 We also adopt metagalactic dr e r
⋆ ⋆ gas
UVbackgroundradiationatz=0byHaardt&Madau(2012) = (4)
dt t
providedby GRACKLE. Forthe differencebetweenthe cho- ff
senUVbackgroundmodelandpreviouscalculationssuchas wherer isthestellardensity,t =(3p /(32Gr ))1/2 isthe
Haardt&Madau(1996)orFaucher-Gigue`reetal.(2009),we local fre⋆e-fall time, and e =1f%f is the star fograms ation effi-
referthereaderstoSection3.3ofKimetal.(2014). ciencyperfree-falltime. W⋆ecautionthattheparametere is
Lastly, a non-thermal Jeans pressure floor is applied to ⋆
notnecessarily equalto the star formationefficiencyparam-
stabilize the scales of the smoothing length (particle-based
eter foundin each code (e.g., for CHANGA and GASOLINE,
see AppendixB). Fora new star particleto spawn, itshould
11Thewebsiteishttps://bitbucket.org/gnedin/agora api/. have at least the mass of a gas particle in the IC, m =
12Thewebsiteishttp://www.nublado.org/. gas,IC
13See,however,footnote9onhowadifferentversionofGRACKLEmay 8.593×104M⊙. NotethatnH,thres adoptedinthisexperiment
affectthecoolingratesforthegaswiththesamemetalfraction(butnotthe isforthisparticularrunonly,andrepresentswheretheJeans
samemetallicityinterpretedbyCLOUDY). polytropeintersectswitha typicalT−r equationofstate in
A
G
Figure2. The500MyrcompositeofgassurfacedensitiesfromSim-noSFwithradiativegascoolingbutwithoutstarformationorsupernovafeedback. Eachframeiscenteredonthegalacticcenter–location O
R
ofmaximumgasdensitywithin1kpcfromthecenterofgasmass. Forvisualizationsoftheparticle-basedcodeshereafter(Figures1-3,14-15,32,34,35)–butnotinanyotheranalysesexceptthesefigures– A
ytusesanin-memoryoctreeonwhichgasparticlesaredepositedusingsmoothingkernels. SeeSection5fordescriptionsofparticipatingcodesinthiscomparison,andSection6.1foradetailedexplanationof C
thisfigure.ComparewithFigure14.Simulationsperformedby:DanielCeverino(ART-I),RobertFeldmann(ART-II),MikeButler(ENZO),RomainTeyssier(RAMSES),SpencerWallace(CHANGA),BenKeller O
(GASOLINE),Jun-HwanChoi(GADGET-3),YvesRevaz(GEAR),andAlessandroLupi(GIZMO). Thefullcolorversionofthisfigureisavailableintheelectronicedition. Thehigh-resolutionversionsofthis M
figureandarticleareavailableattheProjectwebsite,http://www.AGORAsimulations.org/. PA
R
I
S
O
N
.
I
I
:
I
S
O
L
A
T
E
D
D
I
S
K
Figure3. SameasFigure2butforSim-SFFwithstarformationandfeedback. SeeSection6.1foradetailedexplanationofthisfigure. SeealsoSection3.2forthecommonstarformationprescriptionandthe
guidelineforsupernovafeedback,andSection5fortheexactdepositschemeofthermalfeedbackenergyimplementedineachcode.ComparewithFigures15,21,29,32,34,and35.
5
6 J. KIMETAL. FORTHEAGORACOLLABORATION
ourdisks(seeSections5.2and5.4).14 4.2. MinimumHydrodynamicalSmoothingLength
Newstarparticlesinjectenergy,mass,andmetalsbackinto Forparticle-basedcodes(includingGIZMO;seeSection5.9
theinterstellarmedium(ISM)throughcore-collapse(TypeII) andfootnote28),werequirethatthehydrodynamicalsmooth-
supernovae. AssumingtheAGORAstandardChabrier(2003) inglengthsforcollisionalparticlesdonotdropbelow20%of
initial mass function (IMF) and that stars with masses be- the gravitational softening lengths. Unlike the gravitational
tween8and40M⊙ explodeasTypeIIsupernovae,oneType softeningkernel,exactsmoothingkernelchoicesdifferfrom
IIsupernovaoccursperevery91 M⊙stellarmassformed(see code to code, and are detailed for each of the particle-based
Section 3.5 of Kim et al. 2014). With the AGORA recom- codesinSection5. Wealsoreferthereadersinterestedinthe
mendedfittingformulaeEqs. (4)-(6)ofKimetal.(2014)and actualparameterchoicestoAppendixCagain.
the assumed IMF, this single burst is found to release 2.63
M⊙ of metals15 and14.8 M⊙ ofgas(includingmetals). Per 4.3. RefinementStrategy
every91 M stellarmass,thesemetalandmassareinstanta-
⊙ We recommend to mesh-based code groups that a cell is
neouslydepositedintoitssurroundingafteradelaytimeof5
splitinto8childcellsoncethecellcontainsmoremassthan
Myr,alongwithanetthermalenergyof1051ergs. m =8.593×104M (1 gas particle mass in the IC of
We note that exact deposit schemes for energy, mass, and gas,IC ⊙
particle-basedcodes), or 8 collisionlessparticles (disk/bulge
metalsareleftateachparticipant’sdiscretion. We donotin-
starparticlesintheICwithm =3.437×105M ,ordark
tendtooverlyspecifyasinglecommondepositschemewhich ⋆,IC ⊙
willneedtobeinevitablydifferentfromonecodetoanother matterparticleswithmDM=1.254×107M⊙).Thismakesthe
gridsberefinedinafashionsimilartotheLagrangianbehav-
(e.g.,betweenmesh-basedcodesandparticle-basedcodes),as
iorofparticle-basedcodes,andkeepstheratioofcollisionless
we arguedin Section 3.8 of the Project flagship paper (Kim
particlenumberstogascellsapproximatelyunityonaverage.
etal. 2014). Nevertheless,forall mesh-basedcodes(ART-I,
However,exactrefinementstrategiesdifferslightlyfromcode
ART-II,ENZO,andRAMSES),thesamestrategywaschosen:
to code, and are detailed for each of the mesh-based codes
thermalenergy,mass,andmetalsareaddedtothecellwherea
in Section 5.16 We continue to refine the grids down to the
5Myroldstarparticlesitsatthetimeofexplosion,andtothis
cell only. For particle-based codes (CHANGA, GASOLINE, resolutionlimitD x=80pc(seeSection4.1)wherethenon-
thermalpressurefloorkicksin(seeSection3.1).
GADGET-3,GEAR,andGIZMO),eachcode’sdepositscheme
isdiscussedindetailinSection5. InfutureAGORAprojects,
5. PARTICIPATINGCODES
we plan to calibrate different feedback schemes against ob-
servations and against one another. We refer the readers to In this section we introduce the 9 gravito-hydrodynamics
Section7formorediscussiononthisfuturework. codestakingpartinthistest,focusinginparticularonhydro-
dynamicssolvers,refinementschemesformesh-basedcodes
4. COMMONRUNTIMEPARAMETERS (ART-I, ART-II, ENZO, and RAMSES), and supernovafeed-
Here we review the runtime parameters each group is re- back implementations for particle-based codes (CHANGA,
quiredtoadopt,suchasgravitationalsofteningandhydrody- GASOLINE, GADGET-3, GEAR, and GIZMO). We leave out
namicsmoothinglengthsforparticle-basedcodesandrefine- details that are commonly adopted across platforms such as
mentthresholdsformesh-basedcodes. gas cooling (Section 3.1) or star formation (Section 3.2), or
thatwereincludedintheAGORAflagshippapersuchasgrav-
4.1. GravitationalSofteningLengthandFinestMeshSize itational dynamics(Section 5 of Kim et al. 2014). We also
For all codes the gas mass resolution in hydrodynamics pointoutthatthecodesinvolvedinfutureAGORAstudiesare
needs to be set as close as possible to m = 8.593× notnecessarilylimitedtotheonesdescribedherein.
gas,IC
104M . Assumingthatwewishtoresolveaself-gravitating
⊙ 5.1. ART-I
clumpwith64oftheseresolutionelements,thecorresponding
Jeanslengthscalebecomes In ART-I, differential equations of fluid dynamics are in-
tegrated using a shock-capturing Eulerian method described
64m 1/3 in Khokhlov (1998). It uses a 2nd-order accurate Godunov
l =2 gas,IC =348.7pc, (5)
Jeans (cid:18)(4p /3)n (cid:19) solver(Godunov1959)thatevaluatesEulerianfluxesbysolv-
H,thres
ing the Riemann problem at every cell interface (Colella &
andtherefore,fromEq. (2)wechooseaspatialresolutionof Glaz1985). LeftandrightstatesoftheRiemannproblemare
80pc. Thisvalueisusedasthefinestcellsize D x formesh- obtained by piecewise linear interpolation (van Leer 1979).
basedcodes,andasthegravitationalsofteninglengthegravfor IncontrasttootherversionsofART(seeSection5.2.1ofKim
particle-basedcodes. Forallparticle-basedcodestakingpart et al. 2014), ART-I with distinctive star formation and feed-
in the presentstudy, gravity is softenedaccordingto the cu- back recipes (e.g., Ceverino & Klypin 2009; Ceverino et al.
bicsplinekernel(e.g.,Eq. (A1)ofHernquist&Katz 1989). 2014)havebeendevelopedbyA.Klypinandcollaborators.
Forreadersinterestedintheactualparameterchoices,inAp- Theoctree-based,multi-leveladaptivemeshallowsusersto
pendixCweexaminethemeaningsofrelevantparametersin controlthegridstructureattheindividualcelllevel. Forthis
differentparticle-basedcodes. comparison,the ART-I groupusesa 1283 rootgridcovering
a (1.304 Mpc)3 box, then achieves an ∼ 80 pc cell size at
14 AsnotedinKimetal.(2014),starformationprescription parameters
suchasnH,thresore⋆,theinitialmassofstarparticles,andthestochasticity maximum7levelsofrefinement. Themassthresholdsabove
ofstarformation,areallhighlydependentonnumericalresolution.Anideal-
izedtestlikethedisksimulationpresentedhereareessentialtotuneupsuch 16Givendifferencesinrefinementmachineriesamongmesh-basedcodes
parametersforcomputationallyexpensivecosmologicalsimulations. itisimpractical,ifnotimpossible,toimposeanexactlyidenticalrefinement
15 Perunitstellarmassformed,thetotalfractionalejectedmetalmasses criterionacrossallcodes.Weinsteadadoptatrial-and-errorapproachwithin
(oxygenandironcombined)isMZ=2.09MO+1.06MFe=2.09×0.0133+ theguidelinepresentedinSection4.3,whichresultedinallmesh-codeseven-
1.06×0.0011=2.9%. tuallyconvergingtoasimilaroverallgridstructure.
AGORACOMPARISON. II:ISOLATEDDISK 7
Figure4. Cylindrically-binnedgassurfacedensityprofilesat500Myrfor Figure5. SameasFigure4butforSim-SFFwithstarformationandfeed-
Sim-noSFwithoutstarformationorfeedback.Thecylindricalradiusisfrom back.ComparewithFigures22and27.
thegalactic center–location ofmaximumgasdensitywithin1kpcfrom
thecenterofgasmass. Inallanalysesforparticle-basedcodeshereafter–
exceptthegraphicalvisualizationssuchasFigures1-3–rawparticlefields
areused,nottheinterpolatedorsmoothedfieldsconstructedinyt.Shownin
thebottompanelisthefractionaldeviationfromthemeanoftheseprofiles.
SeeSection6.1formoreinformationonthisfigure.They-axisrangeofthe
toppaneliskeptidenticalamongFigures4-7and22foreasiercomparison.
whichacellisadaptivelyrefinedintoanoctof8childcellsare forde-refinementprovidedtheyalso satisfy a numberof ad-
m =8.593×104M andm =3.437×105M forgas ditionalconstraints. Finally, cellsare refined(de-refined)by
gas,IC ⊙ ⋆,IC ⊙
and collisionlessparticles, respectively.17 For the supernova splittingtheminto8(bymerging8childcells).
feedbackschemetodepositthethermalenergyadoptedbyall The pressure floor implemented in ART-II affects cells at
thehighestlevelofrefinementbymodifyingthegaspressure
mesh-basedcodes,wereferthereaderstoSection3.2.
valuesthatentertheRiemannsolver(i.e.,nottheactualpres-
sureortemperaturefields)with
5.2. ART-II
ART-II solves the gravito-hydrodynamicsequations using Pcell=max(PJeans, Pgas) (6)
aparticle-mesh+EulerianAMRapproach. ART-II features =max(nHkBTJeans, Pgas) (7)
MPI parallelization for distributed memory machines, flexi-
where P is the value entering the Riemann solver, P is
ble time-stepping hierarchy, and a variety of unique physics cell gas
the gaspressurefield in the simulation, k is the Boltzmann
modules(e.g.,Gnedin& Kravtsov2011;Agertzet al.2013) B
constant, n =r /m , and T =T(n /n ) with T =
developedbyN.Gnedin,A.Kravtsovandcollaborators. H gas H Jeans J H H,J J
Forthepresentstudy,startingfromauniform1283rootgrid 1800KandnH,J=8cm−3. Thispolytropechoiceisdesigned
covering (1.311 Mpc)3, cells are refined up to 7 additional tomatchthecommonprescriptionEq.(3)withNJeans≃4.For
thesupernovafeedbackschemetodepositthethermalenergy
levels to reach the finest size of 80 pc. Spherical regionsof
adoptedbyallmesh-basedcodes,seeSection3.2.
4 (6, 10) rootgrid cells radius aroundthe box center are al-
ways refined to at least 3 (2, 1) additional levels relative to 5.3. ENZO
therootgrid.The(de-)refinementprocedureconsistsofthree
ENZO is a block-structured adaptive mesh code, devel-
steps. First, cells are marked for refinementif the gas mass
opedbyanopen-source,community-drivenapproach(Bryan
in the cell exceeds0.6m =0.6×8.593×104M , or if
gas,IC ⊙ & Norman 1997; O’Shea et al. 2004; Bryan et al. 2014).18
the cell contains 2 or more dark matter and/or star particles
Among a variety of solver choices, for this comparison the
that were present in the IC. We then use a diffusion step to
3rd-orderaccurate piecewise parabolic method (PPM) is se-
also mark neighboringcells for refinement and thus smooth
lectedto reconstructthe leftandrightstates ofthe Godunov
theshapeoftheregionstoberefined. Bycontrast,cellswith
problem (Colella & Woodward 1984; Bryan et al. 1995),
gasmassesbelow0.2m orwithoutparticlesaremarked
gas,IC alongwithaHarten-Lax-vanLeerwithContact(HLLC)Rie-
mann solver (Toro et al. 1994).A maximum 30% of the re-
17WenotethatART-IandENZOcannotrefinecellsbyparticlenumbers, quiredCourant-Friedrichs-Lewy(CFL)timestepisusedtoad-
butonlybyparticlemasses. Bycontrast, inthereportedruns, ART-IIand
RAMSESrefinecellsbyparticlenumbers.Therefinementcriteriaarechosen vancefluidelements;i.e.,CFLsafetyfactor=0.3.Inaddition
toensureanagreementamongmesh-basedcodesinoverallgridstructures.
Seealsofootnote16. 18Thewebsiteishttp://enzo-project.org/.
8 J. KIMETAL. FORTHEAGORACOLLABORATION
Figure6. Vertically-binned gas surface density profiles at 500 Myr for Figure7. SameasFigure6butforSim-SFFwithstarformationandfeed-
Sim-noSFwithout starformationorfeedback. Theheightistheabsolute back.
verticaldistancefromthex−ydiskplanecenteredonthegalacticcenter,
|zi−zcenter|.Showninthebottompanelisthefractionaldeviationfromthe
meanoftheseprofiles. They-axisrangeofthetoppaneliskeptidentical
amongFigures4-7and22foreasiercomparison.
to solving the conservation equations for mass, momentum m = 8.593×104M , or if the number of dark matter
gas,IC ⊙
andenergy,theequationforinternalenergyisalso solvedin and/orstarparticlesthatarepresentintheICexceeds8.20 We
parallel,andtheconservativeornon-conservativeformulation then mark additionalcells by performinga mesh smoothing
isadaptivelyselectedbasedonalocalestimateoftheenergy operation,expandingtheinitialareabyonecellwidthinevery
truncation errors. This ensures that the gas temperature re- direction.Whennewcellsarecreatedoroldcellsdestroyed,
mainsphysical,eveninhighlysupersonicregions. density,momentumandinternalenergyareusedasaveraging
The ENZO group uses a 643 initial root grid covering a and interpolating variables, thereby preventing a grid point
(1.311 Mpc)3 simulationbox,thenachieves80pcresolution withspurioustemperature.
with maximum 8 levels of refinement. The mass thresholds InRAMSES,thegaspressurefieldincludesthenon-thermal
abovewhichacellisrefinedbyfactorsoftwoineachaxisare pressure support term given by a temperature polytrope
mgas,IC=8.593×104M⊙ and8m⋆,IC =8×3.437×105M⊙ TJeans =m TJ(n′H/nH,J) with meanmolecularweight m , n′H=
for gas and collisionless particles, respectively.17 The non- r gasXH/mH, TJ =1800K, nH,J =8cm−3 and XH =0.76.21
thermalpressurefloorEq. (3)isusedtomodifythegaspres- AsinART-II,thispolytropeapproximatelymatchesthecom-
sureinsidetheRiemannsolver,butnottoaltertheactualgas monpressuresupportprescriptionEq.(3).Newlycreatedstar
energyfield. Forthesupernovafeedbackschemeadoptedby particlesinRAMSEShaveafixedmassofmgas,IC,buttheyare
allmesh-basedcodes,wereferthereaderstoSection3.2. spawnedwithaPoissonprobabilitydistributionwhoseparam-
eters are designed to mimic the local Schmidt law, Eq. (4).
5.4. RAMSES Lastly,forthecommonsupernovafeedbackschemeadopted
RAMSES is an octree-based adaptive mesh code featuring byallmesh-basedcodes,wereferthereaderstoSection3.2.
an unsplit 2nd-order accurate Monotone Upstream-centered
Scheme for ConservationLaws (MUSCL) Godunovscheme 5.5. CHANGA
forthe gaseouscomponent(Teyssier 2002).19 For thiscom- CHANGAisareimplimentationofGASOLINE(seeSection
parison,RAMSESgroupusesaidealgasequationofstatewith 5.6) in the CHARM++ runtime system.22 CHARM++ (Kale
g =5/3, along with the HLLC Riemann solver (Toro et al. & Krishnan1996)23 enablesthe overlapof computationand
1994) and the MinMod slope limiter (Roe 1986). The CFL
communicationand providesadaptiveload balancinginfras-
safety factor for controlling the time step is set to 0.5. The
tructure,allowingCHANGAtoscaletohundredsofthousands
dualenergyformalismadoptedinENZOsimulations(Section
5.3)isalsousedinRAMSESruns. 20 Readers shouldnoticesubtledifferences hereinrefinementstrategies
Forthisstudy,startingfromauniform1283 rootgridcov- betweenRAMSESandothermesh-basedcodes.Newlyformedstarsarecon-
ering(320 kpc)3,cellsarerefinedupto5additionallevelsto sideredaspartofthebaryonicfluid,sotheydonotchangetheparticlerefine-
mentbasedsolelyoncollisionlessparticlesintheIC.
achievean ∼ 80 pc cell size. The refinementprocessworks
21Thismeansthatinordertoretrievegasinternalenergyortemperature
asfollows.First,newrefinementistriggeredonacell-by-cell (e.g.,Section6),thepressuresupporttermneedstobesubtractedoutfrom
basisifthebaryonicmass(gas+newlyformedstars)exceeds RAMSES’spressurefield,whichistheonlyfieldbeingtracked.
22Thewebsiteishttp://www-hpcc.astro.washington.edu/tools/changa.html.
19Thewebsiteishttp://www.itp.uzh.ch/∼teyssier/Site/RAMSES.html. 23Thewebsiteishttp://charm.cs.uiuc.edu/.
AGORACOMPARISON. II:ISOLATEDDISK 9
Figure8. Cylindrically-binned, mass-weighted averages of gas vertical Figure9. SameasFigure8butforSim-SFFwithstarformationandfeed-
heightsforSim-noSFwithoutstarformationorfeedback. Theheightisthe back.
absoluteverticaldistancefromthex−ydiskplanecenteredonthegalactic
center, |zi−zcenter|. Showninthebottompanelisthefractionaldeviation
fromthemeanoftheseprofiles. SeeSection6.1foranexplanationonhow
thisfigureismade.
ofprocessorcores(Menonetal. 2015). The hydrodynamics 5.6. GASOLINE
in CHANGA closely followsthat of GASOLINE. SPH forces GASOLINEisamassivelyparallelSPHcode,firstdescribed
arecalculatedusingthemethodofRitchie&Thomas(2001),
inWadsleyetal.(2004),thathassubsequentlybeenupdated
andenergyisdiffusedusingtheschemeofShenetal.(2010),
with modern SPH features. It contains a subgrid model for
bothof whichprovidinga moreaccuratetreatmentofmulti-
turbulentmixingofmetalsandenergy(e.g.,Shenetal.2010),
phaseISM.Timestepsaredeterminedbytheminimumofan
a timestep limiter by Saitoh & Makino (2009) (see Section
accelerationandaCFLcriterion. Furthermore,thetimesteps
5.5),andageometricdensityestimatorforSPHforceexpres-
of neighbors are kept within a factor of 2 of each other as
sions(seeSection2.4ofKelleretal.2014,foralatestdetailed
in Saitoh & Makino (2009) in order to accurately integrate
descriptionofthecodeanditsperformance).
highlysupersonicflows.
For the current work, the GASOLINE group uses a Wend-
Forthiswork,ak-thnearestneighboralgorithmisusedto
land C4 smoothing kernel (Dehnen & Aly 2012) with
find the N = 64 nearestneighborswhich are smoothed
smooth N = 200 neighbors.25 The same feedback scheme as
withtheWendlandC4kernel(Dehnen&Aly2012)todeter- smooth
CHANGA’s(Section5.5)isimplemented,smoothedwiththe
minehydrodynamicproperties. Unlikeconventionalversions
WendlandC4kernelover64neighbors(notN =200)to
of CHANGA or GASOLINE, the supernovae thermal energy, smooth
better match the amount of mass heated by feedback events
mass, and metalsare directly distributedto the 64 neighbor-
with other particle-based codes. Gas particles that receive
inggasparticles.24 Gasparticlesthatareneighborsofparti-
feedbackcomputetheirrequiredCFLtimestepatthetimestep
clesthatwillexplodeasasupernovaintheirnexttimestepare
prior to receiving feedback, which helps to prevent numeri-
putontimestepssuitablefortheirpost-supernovathermalen-
cal instability and overcooling. GRACKLE coolingis imple-
ergy,preventingthemfrombeingonamuchsmallertimestep
mented by applying a half timestep of cooling, then a full
required in the CFL condition. GRACKLE cooling is imple- timestep of external PdV heating, followed by a final half
mentedbutit doesnotself-consistentlyaccountfor thePdV
timestepofcooling,asinCHANGA(Section5.5).
work or other external sources of energy, a requirement for
CHANGAandGASOLINE’senergyintegration.Therefore,we
5.7. GADGET-3
splittheenergyintegrationintoahalftimestepofGRACKLE
cooling,thenafulltimestepofPdV heating,andfinallyasec- GADGET-3isanupdatedversionof GADGET-2,acosmo-
ondhalftimestepofcooling. logical tree-particle-mesh (TPM) SPH code that was origi-
nallydevelopedbyV.Springel(Springeletal.2001;Springel
2005).26 GADGET-3 hasimportantupdatesfrom GADGET-
24 Thefeedback prescription usedinthis experiment needed areimple-
mentationofthefeedbackroutinenormallyusedinCHANGAandGASOLINE
(e.g.,Stinsonetal.2006). Inparticular, inprevious work, thesupernovae 25 NotethedifferenceinNsmooth fromCHANGAinSection5.5. Dehnen
ratedeterminedbythestellarageandIMFisconvertedtoanenergyinjec- &Aly(2012)showedthatthiskernelcanuselargerneighbornumberswith-
tion“rate”whichisthenincorporatedintothethermalenergyintegrationof outthepairinginstabilitywhichmayeffectivelyremoveresolution,andthat
neighboringgasparticles. Bycontrast,forthisstudy,asupernovaeventoc- doingthisimprovesperformanceonanumberofbasichydrodynamicstests.
cursinstantaneously,makingtherateanill-definedquantityintheirexisting 26 The website is http://www.h-its.org/tap-software-en/gadget-code/ or
energyintegrationmachinery. http://www.mpa-garching.mpg.de/gadget/.
10 J. KIMETAL. FOR THEAGORACOLLABORATION
Figure10. Gasrotationvelocitycurvesat500MyrforSim-noSFwithout Figure11. Same as Figure 10 but for Sim-SFF with star formation and
starformationorfeedback. Thecylindrical radiusandrotational velocity feedback.ComparewithFigure24.
arewithrespecttothegalacticcenter–locationofmaximumgasdensity
within1kpcfromthecenterofgasmass.Showninthebottompanelisthe
fractionaldeviationfromthemeanoftheseprofiles. SeeSection6.2fora
detailedexplanationonhowthisfigureismade.
2, such as domain decomposition and dynamic tree recon- (Durier&DallaVecchia2012),artificialviscosity(Monaghan
structionwhichmayslightlyaltertheN-bodydynamics. The & Gingold1983)supplementedwith the Balsara switch (f
ij
GADGET-3 code used in this comparison is a modified ver- from Balsara 1995), and particle-based time-dependent vis-
sionoftheoriginalGADGET-3 byK.Nagamineandhiscol- cositycoefficient(Rosswogetal.2000).
laborators, which includes pressure-entropy formulation by For this study, the standard cubic spline smoothingkernel
Hopkins (2013), time-dependentartificial viscosity, variable (Monaghan & Lattanzio 1985) in GADGET-2 is used with
smoothing lengths, among others (e.g., Choi & Nagamine N = 50. The feedback energy, mass, and metal injec-
ngb
2012;Thompsonetal.2014;Aoyamaetal.2016). tionintotheISMisimplementedfollowingthestandardSPH
Forthepresentstudy,the GADGET-3groupadoptsaquin- scheme. The implementationcomprisesthe followingsteps.
tic spline smoothing kernel (Morris 1996) with N = 64. Every time a star particle explodes, we first find the nearest
ngb
Theimplementationofsupernovafeedbackisbasedonanup- gasparticles,accordingtotheweightednumberofneighbors
datedversionofTodoroki(2014)thatlargelyfollowsaSedov- asdefinedinSpringel&Hernquist(2002).Adesirednumber
Taylor blast wave method outlined in Stinson et al. (2006, ofneighborsN =50isused.Thenweinjectthermalenergy
ngb
2013,butnottheircoolingshutoffmodel). Theexactmodel andyieldsintotheneighboringgasparticles,weightedbythe
used in the currentwork is fully described in Aoyama et al. SPHsplinekernel. GRACKLE coolingisperformedafterthe
(2016), but, in brief, the implementation comprises the fol- kick step, once gasparticleshave eventuallyreceivedsuper-
lowing steps. Every time a star particle explodes, we com- novafeedbackenergyandoncethe size of thenexttimestep
putethe“shockradius”basedonChevalier(1974)andMcKee is known. Theadiabatic cooling/heatingis first applied, and
& Ostriker (1977), andthen find the gasparticles within the thentheradiativeoneprovidedbyGRACKLE.
radius. We then inject thermal energy and metal yields into
theidentifiedgasparticleswithintheshockradius,weighted 5.9. GIZMO
by the SPH spline kernel. Finally, we note that the results
GIZMO (Hopkins2015)isanewmesh-freeGodunovcode
of this version of GADGET-3 is not representative of all the based on discrete tracers, aimed at capturingthe advantages
GADGET-3codesinthecommunity,becausesomeofthere- ofbothLagrangianandEuleriantechniques.27 Thenumerical
sultsarestronglydependentonthedetailedimplementations
scheme implemented in GIZMO, initially proposed by Lan-
ofbaryonicphysics,suchasstarformationandfeedback.
son&Vila(2008),followstheimplementationofGaburov&
Nitadori (2011) and relies on the discretization of the Euler
5.8. GEAR equationsof hydrodynamicsamong a set of discrete tracers.
GEAR is a self-consistent, fully parallelized, chemo- Unlike in the moving mesh technique, where the volume is
dynamical tree SPH code (Revaz & Jablonka 2012) which partitioned by a Voronoi tessellation, GIZMO distributes the
is built on the publicly available GADGET-2 code (see Sec- volumefractionassignedtothetracersthroughakernelfunc-
tion 5.7; Springel 2005). The simulations reported here are tion. Forthe currentwork, GADGET’sstandardcubic spline
run with the improvements dicussed in Revaz et al. (2016), smoothingkernelis used with Nngb =32. Note, unlike SPH
includingthepressure-entropyformulationproposedbyHop-
kins (2013), individual and adaptive time-stepping schemes 27Thewebsiteishttp://www.tapir.caltech.edu/∼phopkins/Site/GIZMO.html.