Cosmic Reionization On Computers: Numerical and Physical Convergence

0.83 MB·

by  Nickolay Y. Gnedin

DraftversionMarch2,2016
PreprinttypesetusingLATEXstyleemulateapjv.5/2/11
COSMICREIONIZATIONONCOMPUTERS:NUMERICALANDPHYSICALCONVERGENCE
NickolayY.Gnedin1,2,3
DraftversionMarch2,2016
ABSTRACT
InthispaperIshowthatsimulationsofreionizationperformedundertheCosmicReionizationOnComputers
(CROC)projectdoconvergeinspaceandmass, albeitratherslowly. Afullyconvergedsolution(foragiven
starformationandfeedbackmodel)canbedeterminedatalevelofprecisionofabout20%,butsuchasolution
6 is useless in practice, since achieving it in production-grade simulations would require a large set of runs at
1 variousmassandspatialresolutions,andcomputationalresourcesforsuchanundertakingarenotyetreadily
0 available.
2 Inordertomakeprogressintheinterim,Iintroduceaweakconvergencecorrectionfactorinthestarforma-
tion recipe, which allows one to approximate the fully converged solution with finite resolution simulations.
r
a The accuracy of weakly converged simulations approaches a comparable, ∼ 20% level of precision for star
M formation histories of individual galactic halos and other galactic properties that are directly related to star
formationrates,likestellarmassesandmetallicities. Yetotherpropertiesofmodelgalaxies,forexample,their

1 HImasses,arerecoveredintheweaklyconvergedrunsonlywithinafactoroftwo.

Subject headings: cosmology: theory – cosmology: large-scale structure of universe – galaxies: formation –
]
O galaxies: intergalacticmedium–methods: numerical
C
. 1. INTRODUCTION Apractical,althoughnotnecessarilymathematicallyrigor-
h oussolutiontothisdilemmaistoseek”weak”convergence4,
p Asthestudyofcosmicreionizationentersa(nothingshort
i.e. reduced dependence of the simulation results on spatial
- of a true) Renaissance, with numerous observational probes
o (from JWST and 30-meter class telescopes to 21 cm experi- and mass resolution, when the actual numerical model to
r be simulated is adjusted as the resolution changes. This is
ments)justabouttoincreasethevolumeandqualityofobser-
st vationaldatahundreds-fold,thedemandontheorytoremain precisely the goal I have in this paper: sufficiently accurate
a weakly converged simulations will serve as a robust theoret-
onparwiththeforthcomingdataforcesacarefulreexamina-
[ ical model for the interpretation of the future observational
tionoftheaccuracyofthecurrenttheoreticalmodeling.
2 At present, numerical simulations offer the most accurate data, and as the diversity, quality and volume of the data in-
crease, theweaklyconvergednumericalmodels arelikelyto
v andrealistictheoreticalmodelsofreionization. Hence,evalu-
improvefurtherwiththecontinuingeffort.
2 atingtheprecisionofmodernsimulationsisanimportantand
0 timely theoretical task, while we are all eagerly waiting for Inordertopresentthesubjectrigorously,Istartwiththedis-
8 thefloodofnewdata. cussionofstrongandweakconvergencein§2.AsIhavemen-
5 Results from computer simulations can be considered as tionedabove, theactualterminologyhasbeen introducedby
0 solutionstoactualphysicalequationsonlyiftheyarenumeri- Schayeetal.(2015),sothispresentationisnotoriginal,butit
. servesasausefulbasisforthefurtherexposition. In§3Ienu-
1 callyconverged.Inpractice,however,itisoftenexceptionally
0 difficultorevenplainlyimpossibletoreachwhattheEAGLE merateoriginalsimulationsusedinthispaper.Alloftheexist-
ingsimulationworkisunitedunderasingleproject,Cosmic
6 team(Schayeetal.2015)termed”strong”convergence-i.e.,
1 theeffectiveindependenceofsimulationresultsofspatialand ReionizationOnComputers(CROC),thathasbeendescribed
in Gnedin (2014). Section 4 presents both my methodology
: mass resolution (c.f. Springel & Hernquist 2003; Naab et al.
v anditsapplicationtoachievingnumericalconvergenceinthe
2007; Springel et al. 2008; Mayer et al. 2008; Schaye et al.
i globalstarformationhistoryoftheuniverseforCROCsimu-
X 2010a;Vogelsbergeretal.2014;Schayeetal.2015,foranin-
lations. Finally,§5describestheweaklyconvergednumerical
completelist). Withoutaproperconvergencestudyhowever,
r modelsthatwillserveasaphysicalbasisforthenextinstall-
a the results from a numerical simulation may remain suspi-
mentofCROCsimulations.
cioustobesimplynumericalartifacts.
In galaxy formation modeling the situation is exacerbated
2. NUMERICALANDPHYSICALCONVERGENCE
bythefactthatwecanneverformulatetheproblemtosolveas
asystemofbasicphysicalequations,andalwayshavetorely 2.1. StrongandWeakNumericalConvergence
onsubgridmodelsforstarformationandfeedback. Anysuch Generically,anyphysicalquantityextractedfromasimula-
modelisnecessarilyscale-dependent,soastrongconvergence tionneedstobetestedfornumericalconvergence-avaluefar
limitmaynotexistorbephysicallymeaningless(seeSchaye from convergence is meaningless, as it is dominated by nu-
etal.2015,foramoreelaborateddiscussion). mericaltruncationerrors. Inthespecificcaseofcosmological
simulations, the two most crucial numerical parameters are
1ParticleAstrophysicsCenter,FermiNationalAcceleratorLaboratory, spatial and mass resolution. The former is often fixed in ei-
Batavia,IL60510,USA;[email protected]
thercomovingorphysicalcoordinates. InthispaperIprimar-
2KavliInstituteforCosmologicalPhysics,TheUniversityofChicago, ily consider simulations whose spatial resolution ∆r is fixed
Chicago,IL60637USA;
3Department of Astronomy & Astrophysics, The University of
Chicago,Chicago,IL60637USA 4IamagainfollowingtheSchayeetal.(2015)terminologyhere.
2
inphysicalcoordinates,andthedetailsoftheimplementation would not be valid for very small ∆r. For example, an em-
ofthefixedphysicalresolutionintheARTcodearepresented piricallinearKennicutt-Schmidtrelationusedasastarforma-
intheAppendix(§A). tion recipe for CROC simulations is known to work well on
Themassresolutionisconvenienttoparametrizebyaquan- scales of several hundred parsecs (Leroy et al. 2008; Bigiel
tity M ,whichisequaltotheparticlemassinapureN-body et al. 2008; Genzel et al. 2010; Daddi et al. 2010; Bolatto
1
simulationwiththesamenumberofparticles,i.e etal.2011;Bigieletal.2011;Leroyetal.2012,2013;Tacconi
etal.2013), butfailsorchangesonscalesoftensofparsecs
Ω
M ≡ M M , (Onodera et al. 2010; Shetty et al. 2014; Evans et al. 2014;
1 Ω −Ω DM Casasolaetal.2015).
M B
An alternative to strong convergence is “weak conver-
whereM isthedarkmatterparticlemassinthesimulation.
DM gence”, when parameters of the physical model are adjusted
Forexample,asimulationwith20h−1Mpcboxsizeand10243
astheresolutionincreases,
particleswouldhaveM =9.2×105M forouradoptedvalues
1 (cid:12)
of cosmological parameters (ΩM = 0.3036, ΩB = 0.0479, Qˆ = lim Q(∆r,M1|P,pj = fJ(∆r,M1)). (3)
h=0.6814inaflatΛCDMcosmologicalmodel). ∆r,M1→0
In principle, the third important numerical parameter is
Weakconvergenceislessappealing,asthephysicalmodelis
the dynamic range of initial conditions. For example, one
“tweaked”eachtimeresolutionchanges,buttheultimatetest
couldimagineasimulationwith10243 particlesbutwithini-
iswhethertheconvergedvaluesofparameters,
tial conditions set up on a 5123 grid. In practice, however,
overwhelmingmajorityofallnumericalstudiesincosmology pˆ = lim f (∆r,M ),
j J 1
woulduseinitialconditionsona10243 gridfora10243 par- ∆r,M1→0
ticle simulation, and all CROC simulations are set the same
are physically meaningful. If they are, then there is noth-
way too. Hence, the parameter M quantifies not only the
1 ingwrongwithweakconvergence. Infact,theEAGLEteam
massresolutionofthesimulation,butalsothedynamicrange
(Schaye et al. 2015) argued that weak convergence is physi-
ofinitialconditions.
callymoremeaningfulthanformalstrongconvergence,since
Hence, any physical quantity Q extracted from the simu-
anysubgridmodelisphenomenological,and,hence,mustbe
lation(itcanbeliterallyanything,fromaverylocalproperty
scale-dependent.
likethedensityvalueinagivenlocationtoaglobalvalue,like
An alternative interpretation of weak convergence is that
theaveragestarformationratedensity)canbeconsideredto
converged values pˆ are the “correct” ones, but at finite res-
j
beafunctionofnumericalparameters, olution (∆r,M ) parameters may be tweaked to produce a
1
Q= Q(∆r,M ). better converged value, i.e. Qˆ may be close to the uncon-
1 verged value Q(∆r,M |P,p ) for some values of p even if
1 j j
Ifthecosmologicalsimulationswereproducedfromthe“first thefiniteresolutionvalueattheconvergedparametervalues,
principles”, that would be the complete dependence. How- Q(∆r,M |P,pˆ ), issufficientlyfarfromtheconvergedvalue
1 j
ever, while many physical processes are modeled from the Qˆ.
firstprinciplesinCROCsimulations,starformationandstel-
lar feedback are followed with phenomenological models;
2.2. PhysicalConvergence
different star formation models5 may result in different sim-
Just by itself, numerical convergence does not imply the
ulationpredictions, and, hence, simulatedquantitiesalsode-
pend on the physical model (I label it P) and its parameters truth. A numerically converged value may be wrong, if the
underlyingphysicalmodeliswrong. Hence,oneshouldalso
p ,
j
Q= Q(∆r,M |P,p ). (1) consider“physicalconvergence”ofasimulation,i.e. thesen-
1 j sitivityofQˆ totheadoptedphysicalmodel P. Ofcourse,itis
Ideally, thefullynumericallyconvergedvalue Qˆ wouldbe impossibletoconsideranysortofalimitfor P, butonecan
obtained as a limit of Equation (1) for infinitely fine spatial explorethedifferencebetweensimulationresultswithtwoor
andmassresolution,butwithphysicalmodelparameterfixed, moredifferentphysicalmodels,
Qˆ =∆r,lMim1→0Q(∆r,M1|P,pj), (2) δQˆ =∆r,lMim1→0[Q(∆r,M1|P1)−Q(∆r,M1|P2)].
where hereafter I use a hat symbolˆto label a numerically If several sufficiently different physical models give similar
convergedresult(notice,thatsofaritisnotspecifiedhowthe results,itseemsreasonabletoassignhigherconfidenceinthe
doublelimitisactuallytaken;Iwilldiscussthislater). Sucha simulatedvalues.
limit is commonly called “strong convergence” after Schaye ThephysicalmodelusedinCROCsimulationsisfullyde-
etal.(2015),inasensethataphysicalmodelisfullyindepen- scribed in Gnedin (2014). That paper also explored in suf-
dent of numerical resolution. Alas, as was mentioned in the ficient detail the dependence of the star formation model on
Introduction and discussed by Schaye et al. (2015), in prac- some of its parameters, most importantly the molecular hy-
tice this is rarely achievable or useful, as phenomenological drogen depletion time τSF and the feedback model delayed
models typically are only valid for a finite range of spatial cooling time τBW. There is a degeneracy between these
andmassscales. I.e.,evenifthelimitofEquation(2)exists, twoparameters(fasterstarformationcanbecompensatedby
it may make little physical sense, as the physical model P strongerfeedback),buttheremainingparametricfreedomcan
onlybeconstrainedbyobservationalmeasurements.
5 Hereafter, for the sake of brevity, I will use the term “star formation InthispaperIdonotexplorethevariationofthesetwopa-
model”tomeanboththestarformationrecipeandthemodelforthestellar rametersanyfurther,consideringthepreviousstudysufficient
feedback. for now. Hence, all of the results presented below are also
3
performedforthismassresolution. Thisisfurtherdiscussed
TABLE1
MassResolutionChart in§4.2.
Becauseenforcingconstantresolutioninphysicalunitsin-
volves substantial amount of spatial smoothing, and, hence,
Label M1 Nparticlesin Nparticlesin resolution loss, reference simulations using the “Full Chem-
10h−1Mpcbox 20h−1Mpcbox istry”modelarerunatconstantcomovingresolution.
Individual simulations are labeled in the following way: a
simulation label starts with B10 or B20 for the box size and
“LR” 5.9×107M(cid:12) 1283 then is followed by a two-letter identification of the simula-
“MR” 7.4×106M(cid:12) 2563 5123 tionmassresolutionasshowninTable1. Thethirdfieldindi-
“HR” 9.2×105M(cid:12) 5123 10243
“UR” 1.2×105M(cid:12) 10243 catesthephysicalmodelused(“FULL”for“FullChemistry”,
“GD14”forGD14fits,etc). Thefourthfieldgivesthespatial
resolution in parsecs (c.f. “R100” stands for ∆r = 100pc),
if the resolution is fixed in physical units, or the maximum
subjecttovariationofτSF,τBW,andanyothermodelparame- level of refinement (c.f. “L7” for LMAX = 7) for simulations
ternotexplicitlydiscussedhere. with constant comoving resolution. For reference, a simula-
An important component of the star formation model that tion with 10h−1Mpc, high mass resolution (5123 root grid)
has not yet been explored fully is the model for the forma- and L = 7 has spatial resolution in physical units of
MAX
tion of molecular hydrogen. In CROC simulations stars are ∆r=20(11/(1+z))pc.
only allowed to form in molecular gas, hence modeling the For example, B10.HR.FULL.L7 labels a simulation in a
molecular gas is the cornerstone of the adopted star forma- 10h−1Mpc box with 5123 particles, the “Full Chemistry”
tionmodel. ARTsupportstwomethodsformodelingmolec- model, and spatial resolution fixed in comoving units, with
ularhydrogen-(A)ahighlycomputationallyexpensive, full L = 7, while B10.MR.GD14+.R200 labels a simulation
MAX
non-equilibrium chemical network, including H2 formation ina10h−1Mpcboxwith2563 particles,the“GD14+”physi-
on dust, and detailed radiative transfer of ionizing and ultra- calmodel(§B),andfixedphysicalspatialresolutionof200pc.
violet radiation with accurate account for H2 self-shielding Simulations with progressively higher resolution but the
(Gnedin&Kravtsov2010;Gnedin&Draine2014),and(B)a sameboxsize(likeB10.LR,B10.MR,B10.HR,etc)startwith
simplifiedH2formationmodelbasedonfittingformulasfrom initialconditionsthatpreservethesamelarge-scalemodes,so
Gnedin & Draine (2014)6 and radiative transfer of ionizing that the same individual objects can be identified in simula-
radiationonly.Iwillcallthefirstphysicalmodel“FullChem- tions with different mass resolution. I.e., initial conditions
istry” and the second, simplified model “GD14 Fits”. It is for, say, B10.MRareobtainedfromtheinitialconditionsfor
the second, simplified model that has been so far used in all B10.LR by adding small scale power between k (LR) and
Ny
CROCsimulations,asitwouldbeimpracticaltousethe“Full k (MR) = 2 × k (LR), initial conditions for B10.HR are
Ny Ny
Chemistry”modelinlarge-scale,non-zoom-insimulations. furtherrefinedbyaddingsmallscalepowerbetweenk (MR)
Ny
There exists, however, one serious limitation of the GD14 andk (HR)=2×k (MR),etc,wherek (...)istheNyquist
Ny Ny Ny
fittingformulas:theywerecalibratedonlowerredshift(z∼3) frequencyofthecorrespondingsimulationbox.
simulations,anddonotaccountforthefactthatmolecularhy-
drogenformationtimemaybelong,evenlongerthantheage
4. CONVERGINGONTHEGLOBALSTARFORMATIONHISTORY
oftheuniverse. Atz ∼ 3thisisnotanissue,sincetheageof
theuniverseis2Gyr, comparabletothemolecularhydrogen Oneofthemostgenericquantitiesinacosmologicalsimu-
depletion time. At high redshifts, z (cid:38) 10, the effect of the lationisthegloballyaveragedstarformationratedensityasa
finiteageoftheuniverseis, however, substantial; inparticu- functionofcosmictime(Springel&Hernquist2003;Schaye
lar, at z > 10 original GD14 fits would produce unphysical etal.2010b), oftencalled“globalstarformationhistory”. If
globalstarformationrateinexcessof“FullChemistry”sim- anumericalsimulationisnotconvergedonthegloballyaver-
ulations. Hence,IintroduceacorrectiontotheGD14model, agedstarformationratedensity,littleelsewouldbeconverged
callingit“GD14+”. Thiscorrectionisdefinedandvalidated too.Ialsousethissimplequantitytodemonstrateindetailthe
inAppendixB. methodologyusedforexploringconvergence.
3. DESCRIPTIONOFTHESIMULATIONS
4.1. SpatialConvergenceatFixedMassResolution
A complete description of the physics followed in CROC
simulationsispresentedinGnedin(2014)andIdonotrepeat Figure 1 shows such a convergence study for CROC high
ithereforthesakeofbrevity. resolution(M1 = 9.2×105M(cid:12))simulations. Let’sfirstfocus
ForexploringtheconvergencepropertiesofCROCsimula- on colored lines. They show 10h−1Mpc box simulations at
tions,Iuse10h−1Mpcand20h−1Mpcboxes,aslargersimu- severalspatialresolutions(from∆r = 25pcto∆r = 400pc)
lations would be too expensive for a sufficiently dense sam- forthedefault“GD14+”model. Simulationresultsdifferfor
pling of numerical parameters. Simulations that use GD14 differentspatialresolutions,soevenattheresolutionof∆r =
fitting formulas (i.e., prototypes for production simulations) 25pcsimulationshavenotyetfullyconverged. Isitpossible
are run at fixed spatial resolutions, with values ∆r = 25, 50, toextrapolatesimulationresultstofindtheconvergedanswer?
100, 200, and 400pc, except for the “ultra-high” (UR) res- Asthedefaultextrapolationscheme,IadoptaTaylorseries
olution, which would be impractical to sample fully with 5 expansionofthelogofaquantity(inthiscase-theglobalstar
spatialresolutions;onlythelowestresolution,∆r=400pc,is formationratedensityatagivenredshift),
 
ma6tioTnhesifimtutilnagtiofonrsmwuiltahsthareefuolbltanionne-deqbuyilriubnrinuimngnaetwwiodrek.gridofgalaxyfor- Q(∆r)= Q0exp−(cid:88)n Cj∆rj, (4)
i=1
4
102 0.1
Full Chemistry
z=7
GD14+
]
3pc 10 3c] 10-2
M p
M
/ z=10
r /
y r
/ y
fl /
M 1 fl
[ M 10-3
3z) ∆r=400pc ∆r)[
+
1 ∆r=200pc (∗ Q=Q exp( A∆rB )
˙ρ(∗0.1 ∆∆rr==15000ppcc ˙ρ 10-4 Q=Q00exp(−−A∆r−B∆r2)
∆r=25pc Q=Q exp( A∆r B∆r2 C∆r3)
∆r 0 0 − − −
→
10-2 100 101 102 103
6 7 8 9 10 12 14 16
∆r[pc]
z
Fig.1.— Global star formation histories for several “HR” (M1 = 9.2× 0.1
105M(cid:12)) simulations. Various colored solid lines (from red to purple in
rainbow order) show actual simulation results with spatial resolutions of z=7
∆r =400,200,100,50pc,and25pcrespectivelyforourdefault“GD14+”
model,whileasolidblacklinegivesthenumericallyconvergedhistory.Asa
]
testofphysicalconvergence,thedashedblacklineshowsthefullyconverged 3c 10-2
reference“FullChemistry”HRrun. p
M
whereQ isthespatiallyconvergedvalue, / z=10
0 r
y
Q0 = lim Q(∆r). /fl
∆r→0 M 10-3
)[
With 5 numerically sampled values of Q, one can determine r
∆
Q0 anduptofirst4coefficientsCj (n = 4). Inpractice,how- (∗ Q=Q exp( A∆rB )
ever,sincesimulationresultsarealwayssomewhatnoisy,itis ˙ρ 0 −
preferablytouseasmallernumberofdegreesoffreedomthan 10-4 Q=Q0exp(−A∆r−B∆r2)
Q=Q exp( A∆r B∆r2 C∆r3)
thenumberofdatapoints,toavoidfittingnumericalnoise. 0 − − −
Since the extrapolating procedure does depend on the
0 100 200 300 400 500
adoptedfunctionformforthefittingfunction,Ialsoconsider ∆r[pc]
analternative,3-parameterpower-lawfunctionalform,
(cid:16) (cid:17) Fig.2.—Extrapolationsofthesimulatedresultsatz = 7andz = 10to
Q(∆r)= Q exp −A∆rB , (5) thenumericallyconvergedvaluesalongthespatialresolutionaxis.Blackcir-
0
clestracenumericalresultsattheparticularredshift,whilelinesshowvarious
where Q , A and B > 0 are fitting parameters. I also ex- extrapolations(seetextformoredetails). Colors(red,green,andblue)cor-
0 respondtocaseswhen5,4,and3highestresolutiondatapointsareusedfor
perimentedwithotherfittingfunctions, includingusingTay- extrapolation.Anumericallyconvergedresult(extrapolationto∆r→0)can
lor series expansion for the value itself, rather than its log bedeterminedtoabout10%precision. Thetwopanelsarethesamefigure
(Q(∆r) = Q −(cid:80)n C ∆rj), but these two functional forms withlogarithmicandlinearhorizontalaxes(thelog-scaledplotiseasierto
0 i=1 j see,buttheactualamountofinterpolationto∆r → 0requiredisclearerin
givenaboveseemtocapturemostofvariationsduetoadopted
thelinearspace).
parametrization.
As an example of the extrapolation technique, Figure 2 highredshifts,butdeterioratestoabout20%afterz≈10.One
shows cuts through Fig. 1 at two representative redshift val- cannot say a priori which of the two solutions is the “right”
ues, and results of extrapolating numerically sampled values one;whilethe“FullChemistry”modelincludesmorephysics,
with functions (4) and (5), using either all 5 data points or italsotakesthe“blastwave”feedbackmodelintotheresolu-
subsets of 4 or even 3 highest resolution points. Differences tion limit where it is not supposed to work well. Hence, the
betweenthesevariouschoicesserveasanestimateoftheun- differencebetweentheblacklinesinFig.1should,atpresent,
certainties in the extrapolation procedure, which in this case be treated as an estimate of the theoretical error due to the
doesnotexceed10%. assumed H formation model (and there must be other theo-
2
Using the functional form (4) with n = 3 (n = 2 for the reticalerrorsontopofthisduetootherphysicalassumptions).
“Full Chemistry model) as the fiducial one, I can now ex- The interpretation of this figure is as follows. It demon-
trapolate simulated results to the limit ∆r → 0 at each red- stratesthatsimulationswithafinitespatialresolution∆r > 0
shift, and these extrapolations are shown with black lines in dohaveaconvergencelimit∆r → 0;thelimitisslightlyde-
Fig.1. Physicalconvergencecannowbetestedbycomparing pendent on the adopted extrapolating function, but is robust
the fully converged “Full Chemistry” simulations (the black to within about 10%, which can be treated as a “systematic
dashedlineinFig.1)and“GD14+”runs(ablacksolidline). error”ofthenumericallyconvergedresult. Thisnumerically
The (level of) agreement between these three lines demon- converged limit can be tuned to be within the 20% variation
stratesthephysicalconvergenceoftheCROCsimulations;it fromthe”FullChemistry”model,thusallowingonetousea
isnot,ofcourse,achievedserendipitously,butsimplytheout- simplified, ”GD14+” model for molecular hydrogen forma-
come of choosing the right values for the model parameters, tion as a reasonable approximation to a full non-equilibrium
D = 0.05,inEquation(B3). Theagreementisverygoodat calculation.
c
5
102
M 0,∆r=100 pc
z=7 1→
10-2 M1→0,∆r=50 pc
M 0,∆r=25 pc
3Mpc] 3Mpc] 10 M11→→0,∆r→0
M/yr/fl10-3 z=10 M/yr/fl 1
)[ 3[
M1 z)
(∗10-4 Q=Q exp( A MB ) 1+ ∆r 0,M =6 107M
˙ρ QQ==QQ000eexxpp((−−−AA MM111)−, B2 Mpo12in)ts ˙ρ(∗0.1 ∆∆∆rrr→→→000,,,MMM111==790×××110065MMflflfl
10-5 → 1→
104 105 106 107 108 10-2
6 7 8 9 10 12 14 16
M [M ]
1 z
fl
Fig.4.— Fully (i.e. spatially and mass) converged global star formation
histories for CROC simulations. Red lines show the sequence of mass-
z=7 convergedvalues,non-solidlinestrackingthreefixedspatialresolution(out
10-2 of total 5), while the red solid line is the limit M1 → 0,∆r → 0. Blue
linesshowtheoppositeapproach,withthebluesolidlinetracingthelimit
3c] ∆r →0,M1 →0. Solidblueandredlinescoincide,demonstratingthesta-
p bilityofextrapolation.
M
r/ 10-3 z=10 onlygoodtoabout20%,butthatissufficientatpresent,since
y
/ thephysicalconvergenceisonlygoodtoatmostthatlevel.
fl
M
)[ 4.3. FullNumericalConvergence
1
M
(∗10-4 Q=Q exp( A MB ) Having determined the spatially converged results along
˙ρ 0 − 1 spatialormassdirection,onecannowundertakethefulllimit
Q=Q0exp(−A M1), 2 points ∆r → 0,M1 → 0 or M1 → 0,∆r → 0. The results of such
Q=Q0exp(−A M1−B M12) doubleextrapolationareplottedinFigure4,withthetwolim-
10-5 its coinciding almost perfectly. This is, of course, fully ex-
0 107 pected,astheorderoftakingthelimitsin∆r and M should
M [M ] 1
1 not matter. The agreement, however, is not completely triv-
fl
Fig.3.—Extrapolationsofthesimulatedresultsatz=7andz=10tothe ial, since extrapolations along both spatial and mass direc-
numericallyconvergedvaluesalongthemassresolutionaxis. Blackcircles tionsrequirefittingfunctionstonoisydata;ifthefitsbecome
trace numerical results at the particular redshift, while lines show various unstable, the two limits would not agree. The actual agree-
extrapolations(seetextformoredetails).Colors(redandblue)correspondto
mentbetweenthemindicatesthestabilityoftheextrapolated
caseswhen4and3highestresolutiondatapointsareusedforextrapolation.
Anumericallyconvergedresult(extrapolationtoM1→0)canbedetermined values.
toabout20%precision.Thetwopanelsarethesamefigurewithlogarithmic
andlinearhorizontalaxes(thelog-scaledplotiseasiertosee,buttheactual 5. WEAKNUMERICALCONVERGENCEONGALACTIC
amountofinterpolationtoM1→0requiredisclearerinthelinearspace). PROPERTIES
While finding the fully converged result is important, it is
4.2. MassConvergenceatFixedSpatialResolution
not practical to strive for strong convergence in production
A similar procedure can be used to explore mass conver- runs.ComputingthesolidlineinFig.4requiresproducing16
gence. Foratechnicalreason,theARTcoderequiresthesize separate simulations (of which highest mass resolution runs
of the root grid to be a power of two, so the mass resolu- dominatethetotalcomputationalexpense). Theultimategoal
tion can, at present, be only sampled in factors of 8. That ofaconvergencestudyis,then,toachieveweakconvergence-
limits the number of samples we can have in most cases to i.e.,tofindcombinationsofsimulationparametersthatallow
just 3 (LR, MR, and HR). It is hard to estimate the robust- to reproduce the fully converged result in simulations with
ness of the extrapolation procedure with just 3 samples, so finitespatialandmassresolutions.
one additional simulation with “ultra-high” (UR) resolution Inaddition,convergingjustontheglobalstarformationhis-
(M1 =1.2×105M(cid:12))hasbeencompleted. Atsuchanextreme toryisnotenough,sincereproducingtheglobalhistorydoes
mass resolution, it is only practical to run at the lowest spa- notguaranteethatthedetailedmorphologyofreionizationor
tialresolution,∆r =400pc. Hence,asmallcornerofthefull internal properties of galaxies are reproduced. Among these
resolutionsamplegrid(M = 1.2×105M ,∆r ≤ 200pc)re- propertiesthemostimportantformodelingreionizationisstar
1 (cid:12)
mains unexplored in this work. It can be filled in the future, formationhistoriesofindividualgalactichalos. Indeed,sim-
ascomputingpowerincreasesevenfurther. ulationswithvariousmassandspatialresolutionareexpected
Figure3servesasadirectanalogofFig.2,butnowshowing to converge well on the distribution of sufficiently massive
theextrapolationalongthemassdirection. Becausethenum- galactic halos, since modern N-body simulations are highly
ber of data samples is less than in the spatial extrapolation robust(Kimetal.2014). Hence, aslongasthemassresolu-
case, the variation between various extrapolating functional tion of a simulation is enough to resolve galaxies with halo
formsisgreater. Inthemassdirectiontheconvergedvalueis masses M (cid:38) 108M that produce almost all ionizing pho-
h (cid:12)
6
10 10
∆HRr,→ ∆0r,=M110→0p0c, SFE 2 1 11..12 Ionizing budget
1 MR, ∆r=100pc, SFE××3.5 0.1 z≈6 001...890
z 6 6 7 8 9 10
3c] 0.1 ≈ 3c] 10-2 z≈7 z
p p
M 10-2 z 7 M 10-3
yr/ ≈ yr/ z≈8
/ / 10-4
fl10-3 fl
M M z 9
)[ z≈8 )[ 10-5 ≈
Mh10-4 Mh
> > 10-6 z 10
( z 9 ( ≈
˙ρ∗10-5 ≈ ˙ρ∗10-7
z− z−
610 10-6 z 10 610 10-8 ∆HRr+→W0,M,1 ∆→r0=100pc, SFE 1.2
≈ 10-9 HR ×
HR+W , ∆r=200pc, SFE 1.5
10-7 10-10 MR+WHMRR, ∆r=100pc, SFE××4
MR+W , ∆r=200pc, SFE 4
MR ×
10-8 10-11
108 109 1010 1011 108 109 1010 1011
M [M ] M [M ]
h h
fl fl
Fig.5.—Cumulativestarformationratedensityasafunctionofhalomassat Fig.6.—Cumulativestarformationratedensityasafunctionofhalomass
severalredshiftsforthefullyconvergedresult(solidblacklines)andtypical atseveralredshiftsforthefullyconvergedresult(solidblacklines)andfour
finiteresolutionruns(redandbluedashedlines). Notethatwhiletheglobal weaklyconvergedfiniteresolutionsimulations. Withtheweakconvergence
SFRhistoryisapproximatelymatchedbythefiniteresolutionsimulations correctionthefiniteresolutionsimulationsrecoverthewholeshapeofρ˙∗(>
(valuesofρ˙∗atMh→0asafunctionofredshift),thedetailedstarformation Mh)atallhalomasseswithanacceptable(∼ 20%)precision. Theinsert
ratesinhalosofdifferentmassesaresignificantlyoff.Differentredshiftsare showstheratiosofρ˙∗(Mh → 0)fromfourweaklyconvergedsimulations
shiftedverticallyforclarity. andthefullyconvergedsolution,whichistheprecisionwithwhichweakly
convergedsimulationsrecoverthetotalstarformationrateand, hence, the
tons(Gnedin&Kaurov2014),onecanattempttoconstructa globalionizingbudget.
weaklyconvergedmodelthatwillrecoverstarformationhis-
ambiguity during merger events. In addition, it would be a
toriesofindividualgalactichalosand, hence, willreproduce
non-local modification of the star formation recipe. In order
thedetailedreionizationmorphologyofafullyconvergedso-
toavoidthesepurelytechnical,butstillformidablecomplica-
lution.
tions,IintroduceaweakconvergencecorrectionfactorW(n )
Figure5illustratesthispointwithanexampleoftwofinite H
thatisafunctionofthetotalgashydrogendensityn only. It
resolution simulations. It shows a fraction of the global star H
is,obviously,justachoice,andotherformsoftheconversion
formation rate density contributed to by galactic halos of a
factor can be also designed. Hence, the weak convergence
givenmass,inacumulativedistribution,becauseinthatcase
correctionfactorisnotunique;eveniftwodifferentformsof
thevalueofρ˙ at M →0issimplytheglobalstarformation
∗ h W provide similar level of convergence for some particular
history (i.e. asymptotic values of black lines at M → 0 are
h galactic property, like ρ˙ (> M ), they may result in different
simplyvaluesofsolidlinesinFig.4attheseredshifts). ∗ h
degreeofconvergenceforsomeothergalacticproperties. Itis
Black solid lines show the fully converged (M → 0,
1 unclear, however, howonewouldapproachfindingan”opti-
∆r → 0)solutions,obtainedwiththemethodologydescribed
mal”convergencefactor.
in the previous section, and color lines give two typical fi-
Hence, to be specific, I choose the following form for
niteresolutionsimulations,whosestarformationefficiencies
W(n ),
(SFE) are adjusted to match the global SFR history at z = 6 H q +q (n /n )
(isi.eb.ythceonasgtrreuecmtioenn)t.beWtwheileentahlel 3glloibnaelsSfoFrRMhhis→tori0esatfozr=th6e W(nH|q<,q>,nc)= <1+(>nH/Hnc)c , (6)
latterarewithin∼ 25%ofthefullyconvergedsolutionatall where q , q , and n are parameters. Extensive parameter
< > c
times,theoverallshapeofthedistributionisnotcapturedby search results in the following weak convergence correction
thefiniteresolutionsimulations. factorsfor“HighResolution”(M =9.2×105M )runs,
1 (cid:12)
In order to fix this deviation, I introduce in the star for-
mationrecipeamultiplicative”weakconvergencecorrection W (n )≡W(n |3,1,10cm−3), (7)
HR H H
factor”W,
ρ˙ =W×ρ˙(orig), and“MediumResolution”(M1 =7.4×106M(cid:12))runs,
∗ ∗
W (n )≡W(n |3,0.3,10cm−3). (8)
whereρ˙(orig) isthevalueofthestarformationrateatagiven MR H H
∗
locationproducesbytheoriginal,withoutthecorrectingfac- Both factors enhance star formation in low mass halos; the
tor,simulationrecipe. Ideally,onewouldintroduceacorrec- “MediumResolution”formalsoreducesstarformationratein
tionfactorthatdependsonthehalomass-inthatcaseitcan highestdensityregions,toforcealargervariationinρ˙ (> M )
∗ h
simple be computed as the ratio of black and colored lines with the halo mass. The actual values of the parameters q ,
<
from Fig. 5. However, in practice that would require track- q ,andn areratherrobust;n canbevariedbetween3cm−3
> c c
inghalosasthesimulationevolves,andwouldalsointroduce and 30cm−3 with little change, and q and q can also be
< >
7
obtainedbycomparingnumericconvergenceforothergalac-
10
ticproperties-sincethefactorsWarefullyfixedintheseruns,
Ionizing budget
1.2 allothergalacticpropertiesareuniquelypredictedinthesim-
1 1.1
1.0 ulations.
0.1 z≈6 00..896 7 8 9 10 In Figure 8 I show three other galactic properties for the
] z fully converged solution and four weakly converged simula-
3c 10-2 z 7 tionsdiscussedabove. Adecentlevelofconvergenceforstel-
p ≈
M lar masses and metallicities is not surprising, since converg-
r/ 10-3 ingonfullstarformationhistoriesofindividualgalactichalos
y z 8
M/fl10-4 ≈ ceoffnevcetirvgeeldyagsuawraenllt.eeHsothwaetvmero,sat ogfivsetnellsatarrpfroorpmerattiieosnarheiswtoerlyl
M)[h 10-5 z≈9 doesnotconstraintmanyothergalacticproperties,forexam-
ple, atomic hydrogen fraction. As the right panel of Fig. 8,
>
(∗10-6 z 10 thisquantityisconvergedmuchworse,atthelevelofafactor
˙ρ ≈ of2only.Hence,theweaklyconvergedsimulationsshouldbe
6z− 10-7 usedwithcarewhenmakingpredictionsfor,forexample,21
0
1 cm emission - they will predict correctly the emission from
10-8 ∆r→0,M1→0 thelarge-scaledistributionofionizedbubbles(sincestarfor-
B20.MR+W , ∆r=200pc
MR mation rates and, hence, ionizing luminosities of individual
10-9 B20.MR+WMR, ∆r=100pc galactichalos arerecovered), but willoverestimatethe addi-
B20.HR+W , ∆r=100pc
HR
10-10 tional contribution from galaxies themselves, which can be
108 109 1010 1011 importantat∼25−50%level(Kaurov&Gnedin2015).
Mh[M ] An alternative, and likely a better approach (suggested by
fl
thereferee)couldbetomodifytheactualgasdensitythaten-
Fig.7.—Cumulativestarformationratedensityasafunctionofhalomass
atseveralredshiftsforthefullyconvergedresult(solidblacklines)andthree tersstarformationmodelandradiativetransfersolver(butnot
weakly converged finite resolution simulations in 20h−1Mpc boxes (two thehydrosolver,whichneedstoremainstrictlyconservative).
“MediumResolution”setswith100pcand200pcspatialresolutionandasin- It may be possible then with a single factor to achieve weak
gle“HighResolution”runwith100pcspatialresolution).Asemi-transparent convergenceinseveralsimulatedquantities, suchasstarfor-
bandshowstheerrorinthemeancomputedfrom6randomrealizations.
mationrates,metallicities,HImasses,andeventheinterstel-
variedbyupto50%. larradiationfield.Findingsuchauniversalcorrectionwillbe,
Figure6nowshowsρ˙∗(> Mh)vs Mh forfourweaklycon- though,asignificanteffort.
verged simulations, “High Resolution” and “Medium Reso-
lution”runswith100pcand200pc(withtheappropriatead- 6. CONCLUSIONS
justmentintheSFE7).Theoverallshapesarenowreproduced Becausenumericalsimulationsofferthemostaccurateand
atanacceptablelevelofprecision,andtheglobalstarforma- realistic theoretical models of reionization at present, evalu-
tionhistory(akathetotalionizingbudget)arealsorecovered atingtheirprecisionisanimportanttheoreticaltask. Results
towithin10%for100pcrunsandwithin20%for200pcsim- fromcomputersimulationscanbeconsideredassolutionsto
ulations. actual physical equations only if they are numerically con-
Whilemyprimaryfocusinthispaperarespatialandmass verged. Inpractice,however,itisoftenexceptionallydifficult
resolution, the simulation box size is also an important nu- or even plainly impossible to reach ”strong” convergence -
merical parameter. In Figure 7 I show weak convergence i.e.,theeffectiveindependenceofsimulationresultsofspatial
for several simulations in 20h−1Mpc boxes (as compared to andmassresolution(seeeq.2)-andonehastoseekinstead
10h−1Mpc for all the runs presented above): two “Medium ”weak” convergence. With weak convergence the simulated
Resolution” (5123) simulation sets with 100pc and 200pc physicalmodelatfinitespatialandmassresolutionisadjusted
spatial resolution and a single “High Resolution” (10243) as the resolution changes (eq. 3). One can consider such an
run with 100pc spatial resolution. Each “Medium Resolu- adjustment as tuning a finite resolution simulation to repro-
tion” box was simulated 6 times in independent realizations ducetheresultsofthefullyconvergedsolution,acompletely
of initial conditions (similar to runs B20.sf1.uv2.bw10 from legitimateandsensibleapproachtotake.
Gnedin2014),andthefigureshowstheaverageoverall6re- In this paper I show that simulations of reionization
alizations and its error (for the 200pc case; the error for the performed under the Cosmic Reionization On Computers
100pcisessentiallyidentical). Simulationswithtwicelarger (CROC)projectdoconvergeinspaceandmass,albeitslower
box sizes achieve a similar level of weak convergence irre- than one hoped for. A fully converged solution can be ob-
spectively of their spatial or mass resolution, illustrating the tainedatalevelofabout20%precision,andwhichisalsoin-
approximate independence of the accuracy of weak conver- dependentoftheadoptedmodelformolecularhydrogenfor-
gencefromtheboxsize. mation. While this can be considered an important achieve-
Sinceweakconvergencecorrectionfactorsaretunedtore- ment, it is useless in practice, since populating the grid of
coverthefullyconvergedsolutionspecificallyforρ˙∗(> Mh), various resolution values for production grade simulation,
independenttestsofthequalityofthesemodifiedrunscanbe neededforaccurateextrapolationtoformallyinfiniteresolu-
tion,wouldrequireaninsaneamountofcomputationaltime.
7 Remember, thattheSFE(quantifiedbythegasdepletiontimeτSF)is Inordertomakeprogressintheinterim,beforesuchlarge
afreeparameterofthesesimulations,tobefixedbymatchingsomedata- computer allocations become possible, I introduce a weak
forexample,theobservedgalaxyUVluminosityfunctions. Hence,inthis
convergence correction factor in the star formation recipe,
paperIamprimaryconcernedwiththerelativestarformationratesinvarious
runs,withtheabsolutescaletobedeterminedlater,afterthefullyconverged which allows to approximate the fully converged solution
solutioniscomparedtotheobservations. with finite resolution simulations. The accuracy of weakly
8
1 1 1
0.1 z≈6 0.1 z≈6 0.1 z≈6
10-2 z≈7 10-2 z≈7 10-2 z≈7
ΩM/Ω)BhM111000---543 zz≈≈89 Z[Z]∗fl1100--43 zz≈≈89 ΩM/Ω)BhM1100--43 zz≈≈89
6z10M/(−∗1100--76 z≈10 6z10−1100--65 z≈10 6z10M/(−HI10-5 z≈10
1100--98 ∆HHMMRRrRR++→++WWWW0,HHMMMRRRR,,,,1 ∆∆∆∆→rrrr0====121200000000ppppcccc,,,, SSSSFFFFEEEE××××1144..25 1100--87 ∆HHMMRRrRR++→++WWWW0,HHMMMRRRR,,,,1 ∆∆∆∆→rrrr0====121200000000ppppcccc,,,, SSSSFFFFEEEE××××1144..25 1100--76 ∆HHMMRRrRR++→++WWWW0,HHMMMRRRR,,,,1 ∆∆∆∆→rrrr0====121200000000ppppcccc,,,, SSSSFFFFEEEE××××1144..25
10-11008 109 1010 1011 10-9108 109 1010 1011 10-8108 109 1010 1011
Mh[Mfl] Mh[Mfl] Mh[Mfl]
Fig.8.—TestsofnumericalconvergenceforM∗vsMH,Z∗vsMh,andMHIvsMhfortheweaklyconvergedrunsshowninFig.6.Agoodlevelofconvergence
forM∗andZ∗isnotthatsurprising,sinceconvergingonfullstarformationhistoriesofindividualgalactichaloseffectivelyguaranteesthat. Agreementwith
atomicgasfractions-atrulyindependenttest-comesoutsignificantlyworse.
convergedsimulationsapproachesacomparable,∼20%level vergedsimulationsbecomeavailable.
of precision for star formation histories of individual galac-
tichalosandothergalacticpropertiesthataredirectlyrelated
Iamgratefultotheanonymousrefereeforinsightfulcom-
to star formation rates, like stellar masses and metallicities.
ments and an important suggestion for the future improve-
Yetotherpropertiesofmodelgalaxies,forexample,theirHI
ment of the current work. Simulations used in this work
masses,arerecoveredintheweaklyconvergedrunstoonlya
have been performed on the Joint Fermilab - KICP cluster
factoroftwo.
“Fulla” at Fermilab, on the University of Chicago Research
Overall, the weakly converged solutions will serve as rea-
ComputingCentercluster“Midway”,andonNationalEnergy
sonable interim theoretical counterparts for the first obser-
Research Supercomputing Center (NERSC) supercomputers
vational data from JWST, until more expensive, better con-
“Hopper”and“Edison”.
APPENDIX
IMPLEMENTATIONOFFIXEDPROPERRESOLUTIONINTHEARTCODE
Gridcodes(includingAMR)commonlymaintainspatialresolutionfixedincomovingcoordinates,whileSPHcodesoftenhold
resolution fixed in proper coordinates. The latter is a desirable property, since most subgrid physical models work well over a
limitedrangeofspatialscales. Hence,inagridsimulationtheproperresolutionmay,withtime,redshiftawayfromthepreferred
range. Asimple,approximateimplementationoffixedproperresolutioninanAMRcodehasbeenrecentlypresentedbyRosˇkar
etal.(2014). Intheirschemeanewrefinementlevelisaddedeverytimetheproperresolutionredshiftstoomuchawayfroma
presetvalue.
A Rosˇkar et al. (2014) approach is very simple, but it has one undesirable feature: as a new AMR level is added, the actual
spatialresolutionofthesimulationjumpsbyafactorof2. Sincethenewlevelwasheldbackartificially,alargefractionofthe
parent level would have been refined were that level active; as the new level is released, a significant fraction of cells on the
parentlevelgetsrefinedtothenewlevelonaveryshorttimescale(free-falltime-scaleofthenewlevel). Because,aswehave
shown above, CROC simulations do not achieve strong convergence in spatial resolution, such a jump in resolution creates a
comparablejumpintheglobalstarformationrate. Sucha“universe-quake”behavioris,obviously,totallyunphysicalandeven
violatescausality.
Inordertoavoidunphysicalbehavior,adifferentschemeisimplementedintheARTcodethatincludesseveralcomponents.
Slowdown,butnotsuppressrefinementbelow∆r.: Norefinementlevelisheldbackartificially,buttherefinementcriterionis
modifiedsothatonanylevelLbelowthefixedproperresolution∆rthedarkmatterand/orgasmassinacellrequiredfor
refinementtoahigherlevelL+1isincreasedbyafactor(∆x /∆r)3,sothatonlythedensestcellswouldrefinebelowthe
L
fixedproperresolution.
SmooththesourcetermforthePoissonequation.: ThetotaldensitythatsourcesthePoissonequationissmoothedonascale
comparableto∆r. Namely,foragivenvalueof∆rtwoneighboringrefinementlevelsL andL = L +1arefoundsuch
d u d
that∆x ≤∆rand∆x >∆r. ThenthedensityonlevelL thatisusedinthePoissonequationiscomputedas
Lu Ld u
ρ˜ =(1−w)ρ +wρ
Lu Lu Ld
with w = log (∆r/∆x ) and ρ being the actual total density on the refinement level L. On all levels L ≤ L the actual
2 Lu L d
totaldensityρ isusedinthePoissionequation,andforalllevelsL> L thedensityoftheparentcellatlevelL isusedin
L u u
thePoissionequation. Thisprocedureisnotpreciselyaconvolutionwithagivenwindowfunction,butitis(a)linearand
(b)givesρ˜ =ρ for∆r=∆x ,soitservesasasmoothingprocedureforρ˜.
L L L
Maintainpressurefloorconsistentlywithgravity.: Thepressurefloorthatpreventsnumericalfragmentationinthegasismain-
tainedsothatthelocalJeanslengthisalwaysresolvedwithatleast4∆r.
9
107 107
z=5 z=7
M [1,3] 1011M M [1,3] 1011M
106 Mh∈[1,3]×1010Mfl 106 Mh∈[1,3]×1010Mfl
∆rh=∈100pc× fl ∆rh=∈100pc× fl
105 ∆r=200pc 105 ∆r=200pc
∆r=400pc ∆r=400pc
¯/ρ104 ¯/ρ104
ρ ρ
103 103
102 102
101 101
10-1 100 101 10-1 100 101
r(kpc) r(kpc)
Fig.9.—Averagedarkmatterdensityprofilesforhalos(inproperdistancecoordinates)intwomassbinsandforthreespatialresolutionsinB10.HRsimulations.
Verticaldottedlinesshowradiiof200pc,400pc,and800pc-arunwith∆r=200pchastherealspatialresolutionabouttwiceworse.Twopanelsshowz=5
andz=7todemonstratethatthespatialresolutionisindeedreasonablyredshiftindependent.
Averagedensityprofilesinhalomassbinsfortwotestsimulationswith∆r =100pc,∆r =200pc,and∆r =400pcareshown
inFigure9. Asonecansee,therealresolutionofthesimulationsisabout2×∆r,ascouldbeexpectedforagridcode,andthe
resolutionisessentiallythesameattwodifferentredshifts,i.e.itremainsconstantinpropercoordinates.
MOLECULARHYDROGENFORMATIONMODEL
The target spatial resolution for CROC simulations is 100−200pc, and a detailed chemical model for molecular hydrogen
formationwouldnotworkatsuchresolution(Gnedin&Kravtsov2011). Instead,itisappropriatetousethefittingformulasof
Gnedin&Draine(2014)(GD14)thatarespecificallydesignedtoworkonscalesofseveralhundredparsecs.
Thereexists,however,oneseriouslimitationoftheGD14fittingformulas: theywerecalibratedonlowerredshiftsimulations,
anddonotaccountforthefactthatmolecularhydrogenformationtimemaybelong, evenlongerthantheageoftheuniverse.
At z ∼ 3 this is not an issue, since the age of the universe is 2Gyr, comparable to the molecular hydrogen depletion time. At
highredshifts,z (cid:38) 10,theeffectofthefiniteageoftheuniverseis,however,substantial. ThisisillustratedinFigure10,which
comparesareference,“FullChemistry”simulation(runB10.HR.FULL.L7)withatestsimulationwiththesamemassresolution
(“HR”),butimplementingtheoriginalGD14fitsasamodelformolecularhydrogenabundance(runB10.HR.GD14.R100). At
z (cid:38) 11, GD14 fits predict star formation rates (and, hence, molecular hydrogen abundance) in excess of the “Full Chemistry”
model. Since the “Full Chemistry” solves the actual time-dependent chemical network of reactions, it properly accounts for
thetimeavailableformolecularhydrogenformationinlowdustabundanceorlowdensityenvironments. Hence, thehighstar
formationratereturnedbytheGD14testsimulationisunphysical.
Inordertocorrectthisunphysicalbehavior,ImodifytheoriginalGD14fitsasfollows. Themolecularhydrogenfractioninthe
102
]
3c 10
p
M
/
r
y
/
fl
M 1
[
3
)
z
+
1
( 0.1
∗
˙ρ Full Chemistry, L =7
MAX
GD14 Model, ∆r=100pc
GD14+ Model, ∆r=100pc
10-2
6 7 8 9 10 12 14 16
z
Fig.10.—Globalstarformationhistoriesofseveraltest“HR”(M1 =9.2×105M(cid:12))simulations. Ablacklineshowsthereference“FullChemistry”HRrun
with7levelsofrefinement(runB10.HR.FULL.L7).AredlinepresentsthesimulationthatimplementstheoriginalGD14fitsasamolecularhydrogenformation
modelwith100pcresolution(runB10.HR.GD14.R100).LargeexcessoftheGD14modeloverthereferencemodelatz(cid:38)11isamanifestationofthefailureof
GD14fitsatearlytimes.Abluelineshowsthecorrected“GD14+”modelthateliminatestheunphysicalbehavior(runB10.HR.GD14+.R100).
10
10 10
GD14 Model GD14+ Model
z=10 z=10
z=9 z=9
z=7 z=7
1 1
4 +
1 4
GD 0.1 D1 0.1
R G
R
10-2 10-2
L=2(∆r=640pc@z=10) L=2(∆r=640pc@z=10)
L=3(∆r=320pc@z=10) L=3(∆r=320pc@z=10)
L=4(∆r=160pc@z=10) L=4(∆r=160pc@z=10)
10-3 10-3
10-3 10-2 0.1 1 10 10-3 10-2 0.1 1 10
R R
EXACT EXACT
(rFeidg).a1n1d.—z =Co7m(pbaluries)onanodftthhreeaecdtuiffalelryensitmreusloatleudtiomnos.leLcuinlaers-tsoh-oatwommiecdgiaansrvaatiloueRs,≡whρiHle2/sρeHmIiv-terrasnussluthceenotnbeacnodmspshuotewd1u0si%ng-9G0%D1p4efirctesnattiltewroanregdes.hTifwtso,zpa=ne1l0s
showtheoriginalGD14fitsandthemodifiedGD14+modelfromEquation(B3).
gasiscomputedas
R˜
f = , (B1)
H2 1+R˜
whereR˜ istheratioofthemoleculartoatomicgas,andisparametrizedas
q
R˜ = R, (B2)
1+q
where R is given by unmodified GD14 fits, their Eqs. (8-10). In the limit q → ∞, this modified, “GD14+” fit reduces to the
originalGD14fittingformula.
Thereisnorigorouswaytocomputethecorrectionfactorq(shortofrunninga“FullChemistry”simulation),soithastobe
implementedwithamodel. Theaveragerateofformationofmolecularhydrogeninasimulationcellis
n˙
H2 =3.5×10−17cm3/s×D n C ,
n MW H d
H2
where D is the dust-to-gas ratio in the Milky Way units (i.e., in the solar neighborhood D = 1) andC is the clumping
MW MW d
factorthataccountsforthenumericallyunresolveddensitystrictureonsmallscales(sinceformationofH isatwo-bodyprocess
2
-fullerdetails, includingequations, aregiveninGnedin&Kravtsov2011;Gnedin&Draine2014). Aplausibleansatzforthe
suppressionfactorqis
n˙
q∼ H2τ ,
n c
H2
whereτ issometimescale.Fromgeneralconsideration,itishardtoguesswhatitshouldbe.Forexample,itcanbeproportional
c
totheageoftheuniverse,inwhichcaseqscalesasDnH;itcanalsobeproportionaltothefree-falltimetff,inwhichcaseqscales
asDn1/2;itcanalsoscalewithdensityasthecoolingtime,τ ∝1/n ,inwhichcaseqscalesasDwithnodensitydependence.
H c H
Inordertoexplorethischoices,Iusethe”FullChemistry”runB10.HR.FULL.L7andapplyGD14fittingformulatotheactual
simulateddata,averagedoverlessrefinedlevels(thus,mimickinglowerresolution).Suchacomparisonisshownintheleftpanel
ofFigure11(itcanbedirectlycomparedtoFig.5ofGD14). TheagreementbetweentheactualcalculationandtheGD14fitting
formula is not nearly as good as in the GD14 paper, and that is the manifestation of the finite value of τ . I explored all three
c
possibleansatzesfortheqfactordiscussedabove(q ∝ Dnα withα = 0, 1/2, and1), andthebestresultisobtainedforα = 0
H
case. Inthatcasethefactorqcanbeparametrizedas
(cid:32) D(cid:33)1/2
q= , (B3)
D
c