Cosmic Reionization On Computers I. Design and Calibration of Simulations
DD raft version O ctober
1, 2018
Preprint typeset using L A TEX style emulateapj v. 5 / / COSMIC REIONIZATION ON COMPUTERS I. DESIGN AND CALIBRATION OF SIMULATIONS N ickolay Y. G nedin
Draft version October 1, 2018
ABSTRACTCosmic Reionization On Computers (CROC) is a long-term program of numerical simulations of cosmicreionization. Its goal is to model fully self-consistently (albeit not necessarily from the first principles) allrelevant physics, from radiative transfer to gas dynamics and star formation, in simulation volumes of up to100 comoving Mpc, and with spatial resolution approaching 100 pc in physical units. In this method paper wedescribe our numerical method, the design of simulations, and the calibration of numerical parameters. Usingseveral sets (ensembles) of simulations in 20 h − Mpc and 40 h − Mpc boxes with spatial resolution reaching125 pc at z =
6, we are able to match the observed galaxy UV luminosity functions at all redshifts between6 and 10, as well as obtain reasonable agreement with the observational measurements of the Gunn-Petersonoptical depth at z < Subject headings: cosmology: theory – cosmology: large-scale structure of universe – galaxies: formation –galaxies: intergalactic medium – methods: numerical INTRODUCTION
Study of cosmic reionization has been highlighted by thelast decadal survey as one of the most promising areas of as-trophysical research in the current decade. Progress in thisarea directly influences many other fields of astrophysics,from thermal evolution of the Lyman- α forest to propertiesof early galaxies.Because the observational constraints on reionization arelimited, theoretical modeling, including numerical simula-tions, play a relatively larger part in reionization studies thanin many other fields of modern astrophysics. Historically,simulations of reionization were mostly confined to two op-posite limits: simulations of small spatial volumes with de-tailed treatment of relevant physics, or large volume simula-tions with simplified physical modeling (Iliev et al. 2009a;Aubert & Teyssier 2010; Friedrich et al. 2011; Ahn et al.2012; Shapiro et al. 2012, and Trac & Gnedin (2009) for areview of the earlier work). Both approaches su ff er from seri-ous limitations. Small box simulations can model individualionizing sources with su ffi cient physical detail, but fail to ac-count for the large-scale correlations between them. Largebox simulations include these correlations, but, by ignoringgas dynamics, are not able to model ionizing sources self-consistently. The inability of the simulations to include allrelevant scales resulted in a recent surge in semi-numericaland purely analytical approximate methods (Furlanetto et al.2004; Furlanetto & Oh 2005; Choudhury & Ferrara 2005;Furlanetto et al. 2006; Choudhury & Ferrara 2006; Zahn et al.2007a; Mesinger & Furlanetto 2007; Alvarez & Abel 2007;Shull & Venkatesan 2008; Zahn et al. 2011; Mitra et al. 2011;Venkatesan & Benson 2011; Mesinger et al. 2011; Kuhlen& Faucher-Gigu`ere 2012; Alvarez & Abel 2012; Mitra et al.2012; Zhou et al. 2013; Battaglia et al. 2013; Robertson et al.2013; Kaurov & Gnedin 2013; Sobacchi & Mesinger 2014).That’s where Moore’s Law comes to the rescue. The un- Particle Astrophysics Center, Fermi National Accelerator Laboratory,Batavia, IL 60510, USA; [email protected] Kavli Institute for Cosmological Physics, The University of Chicago,Chicago, IL 60637 USA; Department of Astronomy & Astrophysics, The University ofChicago, Chicago, IL 60637 USA relenting exponential increase in the supercomputing powermeans that sooner or later the gap between small- and large-box simulations is going to be bridged. In fact, this time isnow - the new generation of supercomputing platforms thathave recently been and are planned to be deployed in the US ,the so-called “peta-scale” platforms (since they get close to orexceed 10 floating-point-operations per second), are partic-ularly suitable for large-scale simulations of reionization thattreat fully self-consistently the radiative transfer of ionizingradiation and gas dynamics.Taking advantage of this technological progress, we havestarted a Cosmic Reionization On Computers (CROC) projectthat aims, over the course of several years, to produce nu-merical simulations of reionization that model fully self-consistently (albeit not necessarily from the first principles)all relevant physics, from radiative transfer to gas dynamicsand star formation, in simulation volumes of up to 100 co-moving Mpc and with spatial resolution approaching 100 pcin physical units.In this first paper in a series, we focus primarily on the tech-nical aspects of our simulations, such as the description of thenumerical method, simulation design, and the calibration ofsimulation parameters. We present the original scientific re-sults from our simulations in the subsequent publications. NUMERICAL TOOLS
Our main simulation tool is the Adaptive Refinement Tree(ART) code (Kravtsov 1999; Kravtsov et al. 2002; Rudd et al.2008). The ART code is an implementation of the AdaptiveMesh Refinement (AMR) technique with the Fully ThreadedTree data structure (Khokhlov 1998). It includes a wide rangeof physical processes that make it particularly suitable formodeling cosmic reionization. Specifically, the current ver-sion of the code includes the following physical ingredients(in addition to standard ingredients of gravity, dark matter,and gas dynamics).
Cooling and Heating of hydrogen and helium is computed For example, “Stampede” at Texas Advanced Computing Center,“Kraken” at Oak Ridge National Lab, “Hopper” and “Edison” at Livermore-Berkeley Lab, “Mira” at Argonne National Lab, “Blue Waters” at NCSA,etc. a r X i v : . [ a s t r o - ph . C O ] J u l “on the fly”, taking into account all relevant processes ina time-dependent manner, without any assumptions of pho-toionization or collisional equilibrium. Abundance of heavyelements is tracked self-consistently in ART, but, in a mostgeneral case, computing the full dependence of the coolingand heating functions on the incident radiation field is a com-plex task in itself, and cannot be currently implemented ex-actly in cosmological simulations, unless the radiation field isconstant in space (Kravtsov 2003; Wiersma et al. 2009).Since the latter is not a reasonable assumption during reion-ization or in the ISM of galaxies, ART uses an approximatemethod of Gnedin & Hollon (2012) that allows to compute themetallicity-dependent part of the cooling and heating func-tions for an arbitrary time-dependent and spatially-variableradiation field. ART, thus, is able to account for several physi-cal e ff ects that are missed in most other cosmological simula-tions codes, such as suppression of cooling in strong radiationfields, dependence of the LTE temperature on the radiationspectrum, etc (see Gnedin & Hollon (2012) for some repre-sentative examples). Radiative Transfer of ionizing and ultraviolet radiation iscurrently implemented in ART using the Optically Thin Vari-able Eddington Tensor (OTVET) approximation of Gnedin& Abel (2001). While OTVET is an approximation, it hasbeen extensively tested against exact schemes (Iliev et al.2006a, 2009b). The Iliev et al. tests underscored one un-desirable feature of the original ART implementation of theOTVET method - excessive numerical di ff usion around ion-ization fronts. The implementation of the OTVET scheme inART was substantially revised after those tests, and the cur-rently used approach eliminates numerical di ff usion almostcompletely; a full description of our current implementationof OTVET is presented in Appendix C. Thus, OTVET re-mains a highly suitable method for modeling cosmic reion-ization (see Gnedin & Abel 2001, for detailed discussion ofthe limitations and inaccuracies of OTVET).In our simulations, we include ionizing radiation from starsfully self-consistently (in a time-dependent and spatially-inhomogeneous manner), because it is the primary driver ofthe reionization process. Other sources of ionizing radiation(quasars, recombination radiation from helium that can ion-ize hydrogen, bremsstrahlung, etc) we only include in thecosmic background, because these sources are either weaklyclustered (helium recombination radiation) or too rare to sig-nificantly a ff ect the radiation field in a typical region of theuniverse (bright quasars). Both components - the radiationfrom local sources and the radiation from distant sources (i.e.cosmic background) - are treated separately in ART, and thencombined together to derive a single solution of the radiativetransfer equation (see Appendix C). The advantage of this ap-proach is that it allows to account for the contribution of raresources (like quasars) to the cosmic background without ac-tually requiring an impractically large simulation volume.Contributions from helium recombination andbremsstrahlung can be easily computed exactly. Ourmodel for the quasar contribution is presented in AppendixA.Of course, the cosmic background is only important if themean free path of ionizing photons is su ffi ciently large, sothat the radiation field from distant sources is comparable toor above the radiation field from local sources at a typical lo-cation in the universe.Since we are running several independent realizations foreach set of numerical parameters, the post-reionization evolu- tion of the IGM would not be captured correctly if we com-puted the mean free path for the cosmic background fromwithin one simulation box - periodic boundary conditions willextend that box over the whole universe, whereas it is sup-posed to represent just one sub-volume of the universe, andonly the full set of independent realizations should be treatedas a numerical model for the whole universe. Hence, we usethe fit from Songaila & Cowie (2010) to account for LymanLimit absorptions in the cosmic background; the backgroundis then still subject to local absorptions inside shielded re-gions, as captured by the radiative transfer solver in Equation(C7). Radiation from local ionizing sources is absorbed fullyself-consistently with the actually simulated gas distributionin the box (Eq. C4). Molecular Hydrogen chemistry (both gas-phase and dust-based) can be followed in complete detail in ART (Gnedin &Kravtsov 2011). However, since spatial resolution of our sim-ulations ( (cid:38)
100 pc) is too coarse to resolve the scale heights ofgalactic disks, it would make little sense to use the full molec-ular chemistry module in this work. Instead, we use the fittingformulae of Gnedin & Draine (2014, in preparation), derivedfrom a large set of small-volume, high resolution simulations,to reliably account for the environmental dependence of themolecular gas on such ISM properties as dust-to-gas ratio orlocal interstellar radiation field. These fitting formulae aresimilar to the ones presented in Gnedin & Kravtsov (2011),but they also account for the overlap of damping wings ofseparate absorption lines in the Lyman-Werner band at highmolecular column densities.
Star Formation cannot yet be modeled from the first prin-ciples in cosmological simulations, and needs to be imple-mented with a phenomenological “sub-grid” model. In thelast several years an important observational advance hasbeen made in understanding star formation on galactic scales.Both, local (Leroy et al. 2008; Bigiel et al. 2008; Bolatto et al.2011; Bigiel et al. 2011; Leroy et al. 2012, 2013) and interme-diate redshift (Genzel et al. 2010; Daddi et al. 2010; Tacconiet al. 2013) observational studies find that the star formationrate surface density on several-hundred-pc scales correlateswell, and approximately linear, with the surface density ofmolecular gas. We use this observed correlation to define ourstar formation recipe in an entirely empirical manner, Σ SFR = Σ H τ SF , (1)where Σ SFR is the star formation rate surface density, Σ H isthe surface density of the molecular gas (including the contri-bution of helium), and τ SF is the molecular gas depletion timescale. We ignore the slightly sub- or super-linear slopes some-times found in observations, since with our resolution we areonly able to resolve a modest range of surface densities wherethe di ff erence between an exactly linear and a slightly non-linear slopes is negligible.The currently most widely accepted viewpoint is that thedepletion time scale τ SF ≈ − τ SF = . ff ect of varying this parameter on our results be-low. Stellar Feedback is implemented in our simulations with thecurrent industry standard “blastwave” or “delayed cooling”model (Stinson et al. 2009; Governato et al. 2010; Agertz et al.2011; Brook et al. 2012; Agertz et al. 2013a; Stinson et al.2013). While this model is purely phenomenological, it isknown to reproduce many of the observed properties of realgalaxies well, and, hence, is an appropriate numerical tool atpresent. The delay time-scale τ BW is a parameter that we varyas a part of the simulation calibration procedure. Our fiducialvalue of τ BW =
10 Myr is consistent with the usage of thisfeedback model in the field.
Ionizing Radiation from Stars is the dominant contributor tothe global reionization process. The exact amount and spec-trum of that radiation depend on stellar IMF and on local ab-sorption inside the galaxy (usually quantified by “escape frac-tions”). For our modeled stars we use a fixed Kroupa IMF; theshape of the ionizing spectrum is adopted from Starburts99modeling (Leitherer et al. 1999) and is plotted in Fig. 4 of Ri-cotti et al. (2002a). The total UV and ionizing luminosities ofa single-age stellar population with mass m ∗ and metallicity Z ∗ can be computed with Starburst99; we fit numerical resultswith the following formula: L ion = (cid:15) UV . × − m ∗ c Z . ∗ (1 + . Z ∗ ) f ( t ) , where f ( t ) is such that (cid:90) t f ( t ) dt = x (cid:16) . + x (cid:17) + x (cid:0) . + x (cid:1) , and x = t / (3 Myr). At late times ( t (cid:29)
10 Myr) the ionizingemissivity from a single-age stellar population falls o ff withtime more rapidly than UV light. Our fit behaves in betweenthe ionizing and UV emissivities, since our OTVET imple-mentation requires the same time-dependence of the sourcefunction for all radiation bands; that ansatz causes at most afew percent error.The parameter (cid:15) UV is unity for the unattenuated stellar out-put. However, in a numerical simulation with finite spatialresolution some of the absorptions are not accounted for. Forexample, some of ionizing photons are absorbed in the parentmolecular cloud from which stars form, further absorptionsoccur in the atomic ISM on scales below the e ff ective resolu-tion of the radiative transfer solver. To account for all of theseunresolved photon losses, we include the (cid:15) UV factor and treatit as a free parameters of our model.Both, our star formation model and the model for ionizingemissivity ignore the contribution of Pop III stars. This is,necessarily, a simplification. Recent studies of the transitionfrom Pop III to Pop II star formation modes (Ricotti et al.2002a,b; Wise & Abel 2008; Wise & Cen 2009; Wise et al.2012; Muratov et al. 2013a,b) demonstrated, that the transi-tion is rapid and occurs early on, hence the contribution ofPop III stars to the global reionization budget is small. How-ever, if in the future these studies are shown to be incorrect,our simulations should then be considered as a strict lowerlimit that only accounts for ionizing radiation from the cur-rently observed stars and quasars.Another source of ionizing photons that we do not includeis annihilation radiation from dark matter. We leave exploringthis reionization source to future work. DESIGN OF THE SIMULATIONS
Resolution Requirements
The spatial resolution is constrained by the range of scaleson which our star formation model is valid, hence we do notvary it in this paper. The mass resolution is a simulation pa-rameter, and, hence, needs to be set carefully. There are sev- eral physical arguments that can be used to pick up importantmass scales, but the eventual choice of the mass resolutionmust be done in a convergence study.The most strict requirement is that all halos that can coolvia atomic hydrogen line cooling ( M tot (cid:38) (1 − × M (cid:12) )need to be resolved (Gnedin & Ostriker 1997; Gnedin 2000;Iliev et al. 2006b; Trac & Cen 2007; Zahn et al. 2007b; Mc-Quinn et al. 2007; Trac et al. 2008; Pawlik et al. 2009). Thatrequirement may be an overkill, though, since dwarf galax-ies in the Local Group contain only a few stars (e.g. Kravtsov2009, and references therein) - if these galaxies are typical ofdwarf galaxies at high redshift, the mass resolution require-ments may be relaxed.We present the detailed analysis of the mass convergenceproperties of our simulations in the Appendix D. The op-timal resolution for our simulation, measured as the total“equivalent particle mass” M (the sum of masses of a sin-gle dark matter particle and an average mass of a baryoniccell in the absence of refinement - or the mass of a singledark matter particle compensated by the ratio of Ω M / Ω DM ),would be about 10 M (cid:12) , which is equivalent of resolving a,say, 20 h − Mpc box (in comoving units) with 1024 dark mat-ter particles. As Appendix D shows, such a simulation wouldaccount for 80% of all ionizing photons.Unfortunately, such a resolution is still too computation-ally expensive at present, so as our fiducial mass resolutionwe adopt eight times coarser resolution of 7 × M (cid:12) (resolv-ing 20 h − Mpc box with 512 dark matter particles). Withthis coarser resolution we account for about 55% of all ioniz-ing photons and for 70% of ionizing emissivity at z =
6. Atpresent such precision is su ffi cient to match the existing ob-servational data; by the time the JWST significantly expandsobservational samples, faster supercomputing platforms willallow us to increase our mass resoluton to the equivalent of1024 particles in a 20 h − Mpc box.Because in ART the number of dark matter particlesthroughout the simulation remains fixed, while cells of theadaptive mesh are created and destroyed dynamically, it ismore convenient to quantify the resolution in terms of thenumber of dark matter particles. Each simulation starts withthe same number of adaptive mesh cells as dark mater parti-cles, to ensure the consistency of the mass resolution in twomain gravitating components. As the simulation proceeds, thenumber of cells usually grows with time, so that by the end ofthe simulation the number of cells is a factor of several higherthan the number of dark matter particles.Our fiducial simulation series is presented in Table 1. Eachof the simulations in the series is run with additional 6 levelsof refinement (except B20HR, which is run with 5 levels ofrefinement, to maintain the same spatial resolution as otherruns), achieving the same cell size of 125 pc in proper unitsat z = z = ∼ h − Mpc and40 h − Mpc boxes. A simulation with the 80 h − Mpc is cur-rently feasible to complete on the largest available machines,but it is su ffi ciently computationally expensive (requiring 20-30 million CPU hours depending on the platform); with thecomputational resources available to us, we will only be ableto a ff ord one per year beginning with 2014. Running Sets of Simulations
Even our largest “B80” simulations, with the 80 h − Mpcbox size, will be only marginally large enough for obtainingconvergent results on the distribution of sizes of ionized re-gions or on observable properties of high redshift galaxies.Thus, in order to extend the reach of our simulations, we runsets of independent realizations of initial conditions for eachparticular choice of the box size and simulation parameters,accounting for the cosmic variance on the scale of the boxsize.This variance is commonly referred to as the “DC mode”.The ART code supports the DC mode in arbitrary cosmologywithout any approximations, following the method describedin Gnedin et al. (2011). In fact, as has been shown in thatpaper, if a set of independent realizations with a given box sizeproperly accounts for the DC mode, it becomes statisticallyequivalent to a simulation with a several times larger box.Table 2 lists all sets of simulations that we present in this pa-per. We have performed several sets of 20 h − Mpc that serveas our primary parameter exploration data. Most of the sci-ence results we extract from two sets of 3 independent real-izations each of the 40 h − Mpc box. We highlight in boldfacethe simulation sets that we consider “fiducial” - in the nextsection we justify that particular choice.In order to distinguish individual realizations in each sim-ulation set, we label them with capital letters. For exam-ple, the fiducial 20 h − Mpc set (B20.sf1.uv2.bw10) includessix independent realizations A-F. A set B20.sf1.uv2.bw40only includes 3 simulations D, E, and F, which have iden-tical initial conditions to the simulations D, E, and F fromthe fiducial set B20.sf1.uv2.bw10. Hence, we have the abilityto both compare simulations with identical initial conditionsbut varied physical parameters (like B20.sf1.uv2.bw10.D andB20.sf1.uv2.bw40.D) and simulations with identical physicalparameters but di ff erent realizations of initial conditions (likeB20.sf1.uv2.bw10.A and B20.sf1.uv2.bw10.B).We also use comparison between 20 h − Mpc and40 h − Mpc boxes as a rudimentary convergence test.Since every computed physical quantity may have di ff erentconvergence requirements, we do not discuss numericalconvergence in a separate sub-section, but rather include suchdiscussion together with the calibration for each simulationparameter. TABLE 1S imulation P arameters Run Box size Resolution Number of(comoving) (proper at z =
6) DM particlesFiducial Series“B20” 20 h − Mpc 125 pc 512 “B40” 40 h − Mpc 125 pc 1024 “B80” 80 h − Mpc 125 pc 2048 Convergence Study“B20LR” 20 h − Mpc 125 pc 256 “B20HR” 20 h − Mpc 125 pc 1024 TABLE 2C ompleted S imulation S ets Set Id τ SF (cid:15) UV τ BW Stopping Number of(Gyr) (Myr) redshift realizations20 h − Mpc boxesB20.sf1.uv1.bw10 1.5 0.1 10 5 6 [A-F]
B20.sf1.uv2.bw10 h − Mpc boxes
B40.sf1.uv1.bw10 CALIBRATION OF THE SIMULATIONS
Star Formation Model
One of the largest existing observational data sets on thesources of reionization are the UV luminosity functions ofhigh redshift galaxies. Matching them would validate our starformation and feedback models.As our final simulation data, we combine the fidu-cial 20 h − Mpc and 40 h − Mpc sets (B20.sf1.uv2.bw10 andB40.sf1.uv1.bw10). Six independent realizations of the20 h − Mpc box are equivalent, by volume, to 0.75 of a sin-gle 40 h − Mpc box, increasing our fiducial set from 3 to 3.7540 h − Mpc boxes.In order to predict UV luminosities of our model galax-ies, we use the Flexible Spectral Population Synthesis (FSPS)code (Conroy et al. 2010; Conroy & Gunn 2010). One com-plication in computing stellar luminosities in far UV is aproper account of cosmic dust. A fully self-consistent dustmodel would require a complex and computationally expen-sive ray-tracing through the simulated galaxies, and is beyondthe scope of this first paper. Instead, we adopt a simple butreasonable dust obscuration model, which we delegate to theappendix, as it has only modest e ff ect on our results (and noe ff ect at all at z > z = z = . A simple star forma-tion model of Equation (1) is able to reproduce the observa-tions for all z >
5. The agreement becomes worse at lowluminosities at z =
5, and that is not particularly surprising- since our simulations maintain the fixed resolution in co-moving units, the spatial resolution degrades as simulationsevolve, and the feedback model becomes progressively lessaccurate, especially in the low mass galaxies.The sensitivity of our star formation model to the numericalparameters τ SF and τ BW is explored in Figure 2. To compare The exact values of redshifts are matched to Hubble filters used in obser-vations: z = . , . , . , . , . , .
23 22 21 20 19 18 17 16 15 M -10 -9 -8 -7 -6 -5 -4 -3 -2 -1 − z φ ( M ) [ M p c − m ag − ] z ≈ z ≈ z ≈ z ≈ z ≈ z ≈ F ig . 1.— Ultraviolet galaxy luminosity functions from a combination ofsimulation sets B20.sf1.uv2.bw10 and B40.sf1.uv1.bw10 at 6 di ff erent red-shifts. Lines with bands show the average luminosity functions with the rmsvariation over our e ff ective 3.75 40 h − Mpc boxes (the error of the mean is,respectively, a factor of √ . = .
23 22 21 20 19 18 17 16 15 M -10 -9 -8 -7 -6 -5 -4 -3 -2 -1 − z φ ( M ) [ M p c − m ag − ] z ≈ z ≈ z ≈ z ≈ z ≈ τ SF =1 . ,τ BW =10 Myr τ SF =1 . ,τ BW =40 Myr τ SF =0 .
75 Gyr ,τ BW =10 Myr τ SF =0 .
75 Gyr ,τ BW =40 Myr F ig . 2.— Ultraviolet galaxy luminosity functions for simulation sets withvaried parameters of the star formation model. The time-scale for delayedcooling τ BW and the depletion time τ SF (i.e. star formation e ffi ciency) doa ff ect the computed luminosity functions, but not a directly proportional waydue to well-known self-regulation of star formation. apples and apples, we use only 3 independent realizations D-F from our fiducial set B20.sf1.uv2.bw10, and compare themwith simulation sets B20.sf1.uv2.bw40, B20.sf2.uv2.bw10,
22 21 20 19 18 17 M β UV z =6 z =7 z =8 z =9 F ig . 3.— Ultraviolet continuum slopes of galactic SEDs for our simu-lated galaxies as a function of UV luminosity (average - lines, rms - shaded).Squares with error-bars are observational measurements from Bouwens et al.(2013); z = z >
5, and the reason forthe discrepancy is currently unclear. and B20.sf2.uv2.bw40 (which only included 3 simulations D-F each, and were only continued to z = .
7, to save com-putational resources). A longer time-scale τ BW for delayedcooling does have a significant e ff ect on the simulated galax-ies, but not a directly proportional one - luminosity functionsfor the set B20.sf1.uv2.bw40 match the fiducial set very wellif shifted horizontally by about 0.5 magnitude, which corre-sponds to only a factor of 1.6. The sensitivity of our star for-mation model to the star formation time-scale τ SF (or, equiv-alently, to star formation e ffi ciency) is stronger - a change in τ SF by a factor of 2 makes a similar 0.5 magnitude change ofthe simulated luminosity functions, but the e ff ect is still sub-stantially sub-linear. Finally, since the two parameters haveopposite e ff ects, a higher star formation e ffi ciency (lower τ SF )can always be compensated by a stronger feedback (longer τ BW ), as is illustrated by the set B20.sf2.uv2.bw40. That re-sult should not be surprising at all - it is well established thatstellar feedback “self-regulates” star formation on kiloparsecscales (Schaye et al. 2009; Agertz et al. 2011; Hopkins et al.2011; Agertz et al. 2013b; Hopkins et al. 2013). The sensi-tivity of our star formation model to numerical parameters issomewhat higher than is usually found at lower redshifts, re-flecting the fact that complete self-regulation takes some timeto get established; that is also consistent with prior work (e.g.Schaye et al. 2009; Agertz et al. 2011).Another observational constraint on our star formationmodel is o ff ered by the observed UV continuum slopes ofhigh redshift galaxies. To measure the UV continuum slopesof model galaxies, we compute the monochromatic fluxes at1300Å, 1400Å, 1500Å, 1600Å, and 1700Å for our fiducialB20.sf1.uv2.bw10 set , and fit a power-law relation F λ ∝ λ β using a simple least-squares fit in a log-log plane. In Figure3 we show the average UV continuum slopes measured fromour simulations and the rms scatter around them, as well asobservational measurements at z = −
7. Our simulationsmatch the observed slopes at z = z = z = Such calculations would be too computationally expensive for all ofB40.sf1.uv1.bw10 simulations, and are not really needed.
The latter discrepancy is somewhat surprising, since bymatching the whole evolution of luminosity functions, oursimulations reproduce not only star formation rates at a givenredshift, but the whole prior star formation histories of galax-ies. Never-the-less, it is possible to come up with severalpotential reasons for the discrepancy: our dust obscurationmodel (Appendix B) may be over-simplistic (since it is in-dependent of the galaxy luminosity), the FSPS code may notbe su ffi ciently accurate, as may be observational, broad-bandbased diagnostics for the true spectral slope; finally, our starformation and stellar feedback models may be too crude. Forall these reasons, we leave resolving this apparent discrepancyto future work. Ionizing Emissivity
Our star formation model, specified by parameters τ SF and τ BW , appears to work reasonably well. The last remaining pa-rameter in the simulations that needs to be calibrated is theescape fraction up to the simulation resolution limit (cid:15) UV . Thisparameter can be calibrated with the observed spectra of high-redshift quasars - ionizing emissivity of our sources is propor-tional to (cid:15) UV , and, therefore, the whole process of reionizationand its aftermath - the Ly α forest at z = − ff ected by (cid:15) UV .In order to model absorption spectra of high-redshiftquasars, we compute synthetic Ly α spectra along 1000 totallyrandom lines of sight at several simulation snapshots. Sincesimulations use periodic boundary conditions, we extend eachline of sight to twice the simulation box size - random linesof sight in a periodic universe are only weakly a ff ected by theartificial periodicity if they do not extend beyond twice thebox length. For the 20 h − Mpc box that corresponds to about ∆ z = .
15 at z =
6, which is a typical redshift resolution foraveraging properties of Ly α spectra in observed spectra (Fanet al. 2006b).In Figure 4 we show the evolution of the mean Gunn-Peterson optical depth and the volume weighted mean HI frac-tion in four of our simulation sets: three 20 h − Mpc boxeswith (cid:15) UV = .
1, 0.2, and 0.4, and our fiducial 40 h − Mpc box.There are several conclusions that can be drawn from that fig-ure.Firstly, changing ionizing emissivity by a factor of 2changes the moment of overlap of ionizing bubbles, indi-cated by the sharp drop in volume weighted mean HI fraction(Gnedin 2000, 2004; Gnedin & Fan 2006), by about ∆ z ≈ . α forestare best matched by (cid:15) UV = . h − Mpc - the fact that outfiducial sets B20.sf1.uv2.bw10 and B40.sf1.uv1.bw10 havedi ff erent values of the (cid:15) UV parameter, but similar reioniza-tion histories, indicates that even 3 independent realizationsof a 40 h − Mpc are not enough to obtain an accurate predic-tion for the mean Gunn-Peterson optical depth or the volumeweighted mean neutral fraction.Finally, our simulations do not match the observationalpoints at z (cid:38) α forest at eachredshift interval way better than the observational measure-ments, and yet they are still far from convergence. Hence, theobservational values are not the converged results either, andthe discrepancy between the observational constraints and our τ G P z -5 -4 -3 -2 -1 › X H I fi B ,† UV =0 . B ,† UV =0 . B ,† UV =0 . B ,† UV =0 . B ,† UV =0 . F ig . 4.— Gunn-Peterson optical depth (top) and volume-weighted HI frac-tion (bottom) vs redshift for four simulation sets: B20.sf1.uv1.bw10 (blue),B20.sf1.uv2.bw10 (green), B20.sf1.uv4.bw10 (red), and B40.sf1.uv1.bw10(black). Solid lines show the average quantities, while the shaded regionsfor fiducial sets B20.sf1.uv2.bw10 and B40.sf1.uv1.bw10 span the limits ofvariation among 6 / τ G P z -5 -4 -3 -2 -1 › X H I fi F ig . 5.— Average Gunn-Peterson optical depth (top) and volume-weightedHI fraction (bottom) vs redshift for 6 independent realizations of the fiducialset B20.sf1.uv2.bw10. Data points are from Fan et al. (2006b). simulations does not necessarily imply a major failure of ourphysical model, but may also be due to just cosmic variance.Large cosmic variance is further illustrated in Figure 5,where we show the average Gunn-Peterson optical depth andthe volume weighted mean HI fraction for all 6 independentrealization of our fiducial B20.sf1.uv2.bw10 set. As one cansee, the variations from realization to realization are very large(recall, that the 20 h − Mpc box is well-matched to the redshiftbin of observational measurement of ∆ z ≈ . τ GP and the gas density, whichmakes τ GP extremely sensitive to outliers. For example, con-sider N lines of sight at some redshift during reionization, ofwhich only one has non-zero transmitted flux F . The aver-age flux over N lines is F / N , and the corresponding Gunn-Peterson optical depth is τ GP = − ln ( F / N ) = τ + ln N , where τ is the Gunn-Peterson optical depth along that oneline of sight with partial absorption. In observations, typicallyonly a few lines of sight contribute to a given redshift bin,i.e. N = −
5. Hence, a single line of sight with τ = z = .
2) results in the measured“average” of τ GP ≈ . − .
6, only slightly above τ .Based on this reasoning, we use only post-reionizationLy α data ( z <
6) for calibrating our simulations, and relyon the derived volume-weighted mean HI fraction from Fanet al. (2006b) as preferred calibration data, even if the Gunn-Peterson optical depth is a directly observed quantity - (cid:104) X HI (cid:105) V is, e ff ectively, a convolution over the whole distribution func-tion of observed Ly α fluxes in individual spectral pixels, justlike τ GP , but it is more sensitive to typical values of the trans-mitted flux, while τ GP is heavily weighted towards the tail ofthe distribution. In that sense it is a less biased, even if lessdirect, observational probe.Hence, our preferred value for the parameter (cid:15) UV is between0.1 and 0.2 (with 0.2 matching observations better); the com-putational expense of simulations (and lack of convergencebetween 20 h − Mpc and 40 h − Mpc boxes) prevents us fromactually fitting for the value of (cid:15) UV to higher precision. CONCLUSIONS
Cosmic Reionization On Computers (CROC) project is along-term simulation campaign for modeling the process ofcosmic reionization in su ffi ciently large simulations volumes (above 100 Mpc) with detailed physical modeling and spatialresolution better than 0.5 kpc (simulation cell size of less than150 pc).A simple model of star formation and feedback, based onthe linear Kennicutt-Schmidt relation in the molecular gas anda widely used “delayed cooling” or “blastwave” feedback, isable to reproduce the observed galaxy UV luminosity func-tions in the whole redshift range z = −
10 with a value forthe molecular gas depletion time of τ SF = . z ∼ z = − (cid:15) UV parameter, a quantity thatdescribes photon losses on scales unresolved in our simula-tions, of (cid:15) UV = . − . α forest in the spectra of SDSS quasars at z <
6. An even betterconsistency can be achieved by fitting the simulations to thedata, albeit at (presently unrealistically) large computationalexpense.The observed increase in the Gunn-Peterson optical depthat z > ∆ z ≈
1. Since theobservational constraints have even less statistical power thanour simulations, they have not yet converged on the true evo-lution of the average Gunn-Peterson optical depth, and may,therefore, be significantly biased.We are grateful to Andrea Ferrara and Matthew McQuinnfor valuable comments and suggestions that significantly im-proved the original manuscript.Simulations used in this work have been performed onthe Joint Fermilab - KICP cluster “Fulla” at Fermilab, onthe University of Chicago Research Computing Center clus-ter “Midway”, and on National Energy Research Supercom-puting Center (NERSC) supercomputers “Hopper” and “Edi-son”.
APPENDIXMODEL FOR QUASAR SOURCES
For the intrinsic quasar SED we use our own fit to the Richards et al. (2006) model, ν L ν ∝ . × − (1 + ( h ν/
300 eV)) . e − h ν/
500 keV + . / h ν ) . (cid:16) + (9 eV / h ν ) (cid:17) . + . e − h ν/ . . The evolution of the bolometric quasar luminosity function has been determined by Hopkins et al. (2007). We fit the quasarbolometric luminosity density as a function of redshift as L QSO ( z ) = . L (cid:12) Mpc (cid:18)(cid:104) e . z − (cid:105) / + (cid:104) e − . z (cid:105) / (cid:19) − , and we use a k -correction of k ion = . DUST OBSCURATION MODEL
In this work we adopt a simple, but reasonable dust obscuration model for our simulated galaxies, in which the dust opticaldepth at wavelength λ is estimated as τ λ = σ D ( λ ) Z ∗ N e ff , (B1)where N e ff is a parameter to be calibrated from the observational data, Z ∗ is the stellar metallicity in solar units, and σ D ( λ ) is thedust cross-section at solar metallicity. For the specific extinction law, we take the SMC dust (Weingartner & Draine 2001), since
22 21 20 19 18 17 M τ z =5 z =6 z =7 F ig . 6.— Dust optical depth at λ = ff erent colors correspond toredshifts in the rainbow order. most of our simulated galaxies at z = − σ D ( λ ) = . × − cm (cid:32) λ (cid:33) . . Figure 6 shows the observational constraints on the dust abundance in z = z = N e ff ≈ × cm − . At z > N e ff : N e ff = × cm − min (1 , max (0 , . − z )))( N e ff = × cm − for z ≤ N e ff = z ≥
8, and is linearly interpolated in between to avoid discontinuities). With thatansatz the dust obscuration at z = OPTICALLY THIN VARIABLE EDDINGTON TENSOR APPROXIMATION IN THE ART CODE
Two-field Ansatz for the Radiation Field
Consider a radiative transfer equation in the expanding universe, in comoving reference frame, ac ∂ J ν ∂ t + (cid:126) n ∂ J ν ∂(cid:126) x − aHc (cid:32) ν ∂ J ν ∂ν − J ν (cid:33) = − k ν J ν + S ν , (C1)where J ν ( t , (cid:126) x , (cid:126) n ) is the radiation specific intensity, k ν is the absorption coe ffi cient (per unit comoving distance), and S ν is thesource function (in appropriate units).In the following, it is assumed that A1: all sources have the same spectral shape, i.e. the frequency dependence of S ν can be factored out, and that A2: sources are isotropic ( S ν does not depend on (cid:126) n ).Both these assumptions can be relaxed, if necessary.With these assumptions, we can write S ν ( t , (cid:126) x ) = L ν ρ S ( t , (cid:126) x ) . For example, for stellar sources, ρ S ( t , (cid:126) x ) can be the mass density of massive stars.In cosmological simulations, it is convenient to replace specific intensity J ν with two separate functions, f ν and g ν , as follows: J ν = ¯ J ν f ν + L ν ( g ν − ¯ g ν f ν ) , (C2)where ¯ J ν is the spatially and angle averaged specific intensity (i.e. cosmic background), which satisfies the following equation, ac ∂ ¯ J ν ∂ t − aHc (cid:32) ν ∂ ¯ J ν ∂ν − J ν (cid:33) = − ¯ k ν ¯ J ν + ¯ S ν , where we defined the mean absorption coe ffi cient as radiation-field-weighted,¯ k ν ≡ (cid:104) k ν J ν (cid:105) ¯ J ν . Analogously, ¯ g ν is the spatially and angle averaged g ν .The reason to impose such an ansatz is to simplify the frequency dependence of the radiation field - in principle, one wouldneed to follow several hundred radiation fields, one for each frequency bins, in order to compute accurately most of ionizationand other chemical rates. This is not practical, obviously. Hence, the goal of ansatz (C2) is to concentrate most of the frequencydependence in the pre-factors, and hope that both f ν and g ν depend on the frequency only moderately.For example, one can imagine that the radiation field shining on a particular place in space is a combination of a contributionof the cosmic background (perhaps, attenuated by additional local absorption) and the radiation from nearby sources (perhaps,also attenuated by additional local absorption). In that case f ν and g ν would only have to account for the local absorption, and inplaces where the local absorption is negligible, would become completely frequency independent.The complete details of our approach to modeling the frequency dependence are given in Gnedin & Abel (2001). Here we justmention that we represent the full frequency dependence of f ν and g ν as f ν = f exp − (cid:88) j q j ( ν ) w j , where q j ( ν ) are a set of pre-determined basis functions (not necessarily mutually orthogonal), and w j are weights, which need tobe computed by the RT solver (the expression for g ν is completely analogous). In the present implementation we choose q j ( ν )to be photo-ionization cross-sections for HI, HeI, and HeII (i.e. the sum over j includes just 3 terms), and weights w j and thenormalization f (or g ) are derived by sampling f ν (or g ν ) at HI, HeI, and HeII ionization thresholds as well in the opticallythin limit. The advantage of this approach is that in the most likely physical situation - when the opacity at a given location isdominated by a single system along a line of sight to a radiation source (or to the radiation background in case of g ν ) - the fullfrequency dependence of the incident radiation field is recovered exactly.The specific form of ansatz (C2) is dictated by the need to satisfy the consistency condition, since J ν now depends on its ownaverage; averaging both sides of Equation (C2) over space and angles gives¯ J ν = ¯ J ν ¯ f ν + L ν (¯ g ν − ¯ g ν ¯ f ν ) , which can be satisfied if ¯ f ν = f ν and g ν remain always non-negative, which is a desirable property for a numerical implementation (as it avoids numerical lossof precision).Equation (C2) is nothing more than an ansatz, we replaced one unknown function with two, hence we can impose a conditionon these two functions. The condition we impose is that terms with f ν and g ν cancel out separately. In the following we dropthe frequency subscripts for brevity, all quantities except ρ S remain functions of frequency. In addition, we adopt a Newtonianapproximation for f and g , omitting terms with the Hubble parameter for them, but retaining cosmological terms for ¯ J ν (seeGnedin & Abel 2001, for more detailed description of this approximation). With these simplifications, one obtains f (cid:104) − ¯ k ¯ J + ¯ S (cid:105) − L f ac ∂ ¯ g ∂ t + ( ¯ J − L ¯ g ) D fdl + L Dgdl == − ( ¯ J − L ¯ g ) k f − Lkg + L ρ S , (C3)where we use a shorthand Ddl ≡ ac ∂∂ t + (cid:126) n ∂∂(cid:126) x for the derivative along the light cone.The condition we impose is then Dgdl = − kg + ρ S , (C4)which, in Newtonian limit (ignoring terms with 1 / c ) has a simple solution, g ( t , (cid:126) x , (cid:126) n ) = (cid:90) ∞ dl ρ S ( t , (cid:126) x + (cid:126) nl ) e − τ ( (cid:126) x ,(cid:126) x + (cid:126) nl ) , (C5)where τ ( (cid:126) x , (cid:126) x ) is the optical depth between points (cid:126) x and (cid:126) x , τ ( (cid:126) x , (cid:126) x ) = | (cid:126) x − (cid:126) x | (cid:90) ds k (cid:0) (cid:126) x + s ( (cid:126) x − (cid:126) x ) (cid:1) . An even more familiar form is the angle average of g , (cid:104) g (cid:105) ( t , (cid:126) x ) = π (cid:90) d x (cid:48) ρ S ( (cid:126) x (cid:48) )( (cid:126) x − (cid:126) x (cid:48) ) e − τ ( (cid:126) x ,(cid:126) x (cid:48) ) , ρ S / (4 π r ) over all sources, diminished by the opacity between the source and the current location.In particular, g is manifestly positive everywhere in the computational domain.Using Equation (C4) in (C3), we find ( ¯ J − L ¯ g ) (cid:34) D fdl + k f (cid:35) = f (cid:104) ¯ k ¯ J − L ¯ ρ S (cid:105) + L f ac ∂ ¯ g ∂ t . Averaging Equation (C4) over space and angle results in ac ∂ ¯ g ∂ t = −(cid:104) kg (cid:105) + ¯ ρ S . Combining the last two equations together, we find( ¯ J − L ¯ g ) (cid:34) D fdl + k f (cid:35) = f (cid:104) ¯ k ¯ J − L (cid:104) kg (cid:105) (cid:105) . (C6)Finally, we can use Equation (C2) to compute the average absorption,¯ k ¯ J ≡ (cid:104) kJ (cid:105) = ( ¯ J − L ¯ g ) (cid:104) k f (cid:105) + L (cid:104) kg (cid:105) . Substituting (cid:104) kg (cid:105) from the above equation into Equation (C6), we obtain the final equation for the function f , D fdl = − k f + f (cid:104) k f (cid:105) . (C7)Equations (C4) and (C7) are our master equations for the cosmological radiative transfer. At this point no approximations havebeen made except for the assumptions A1 and A2 above.One undesirable property of Equation (C7) is that it numerically unstable. To see that, we can average it over space and angle, ac ∂ ¯ f ∂ t = (cid:104) k f (cid:105) ( ¯ f − . (C8)Value ¯ f = ∂ ¯ f ∂ t > f > ∂ ¯ f ∂ t < f <
1. To circumvent thisproblem, we multiply the last term in equation (C7) by a function q ( ¯ f ), D fdl = − k f + q ( ¯ f ) f (cid:104) k f (cid:105) , (C9)where q (1) =
1. It is easy to show that, if q (cid:48) (1) < −
1, then Equation (C9) is numerically stable.A fiducial choice for q is q ( x ) = x − , but, in the future, other forms for that function need to be explored. OTVET Approximation in the Two-field Ansatz
In Gnedin & Abel (2001) it is shown how to derive a single di ff usion-like equation for the angle average of fields f and g .Namely, if F ν ( t , (cid:126) x ) ≡ π (cid:90) d Ω f ν ( t , (cid:126) x , (cid:126) n ) , G ν ( t , (cid:126) x ) ≡ π (cid:90) d Ω g ν ( t , (cid:126) x , (cid:126) n ) , then (again omitting the frequency dependence for brevity) ∂ G ∂ξ = ∂∂ x j k ∂ Gh i jG ∂ x i − kG + ρ S , (C10) ∂ F ∂ξ = ∂∂ x j k ∂ Fh i jF ∂ x i − kF + q ( ¯ F ) F (cid:104) kF (cid:105) , (C11)where d ξ = ˆ c dt / (2 a ), ˆ c ≤ c is the “reduced speed of light” (Gnedin & Abel 2001), and averaging in Equation (C11) is done overthe space (obviously, ¯ f = ¯ F ).We choose Eddington tensors di ff erently in Equations (C10) and (C11): for a cosmological simulation in a periodic box, h i jG is chosen as the optically thin Eddington tensor from all sources inside a periodic box, while the Eddington tensor for thebackground radiation is taken to be isotropic, h i jF = δ i j / Elliptic Solver for OTVET Di ff usion-like Equation Consider OTVET-type equation for some function E ( ξ, (cid:126) x ) and some tensor h i j ( ξ, (cid:126) x ), ∂ E ∂ξ = ∂∂ x j (cid:32) k ∂ Eh i j ∂ x i (cid:33) − kE + s , (C12)where k ( ξ, (cid:126) x ) is the absorption coe ffi cient and s ( ξ, (cid:126) x ) is the source term. This equation is discretized in some fashion in spaceon a set of indicies { I } , where I may may an SPH particle number, a set of indicies ( i , j , k ) on a regular mesh, or any otherdiscretization. We assume that the discretization is such that all quantities E I , h i jI , k I , and s I are co-located in space on the sameset of resolution elements { I } .The discretization scheme is associated with a spatial scale ∆ x (cell size on a regular mesh, SPH kernel size, etc). Each cycleof the OTVET solver (for example, a time-step of a hydro scheme) consists of the set of consecutive iterations, which we labelas E ( n ) I , where n = , , , ... .The second-order term in Equation (C12) is discretized as ∂∂ x j (cid:32) k ∂ Eh i j ∂ x i (cid:33)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) I ≈ ∆ x ˆ D I [ E ; a ] , where a I ≡ k I ∆ x is the dimensionless absorption coe ffi cient and ˆ D I [ E ; a ] is a linear operator on the set of all values E I ,ˆ D I [ E ; a ] = (cid:88) J w I , J E J , with the sum being over all indicies J . Dimensionless weights w I , J depend on various a I as appropriate, and most of w I , J are zerobecause the second order term only has a finite support.For example, for a regular mesh with I = ( i , j , k ),ˆ D I [ E ; a ] = U xi + , j , k − U xi , j , k a i + / , j , k + U xi − , j , k − U xi , j , k a i − / , j , k + U yi , j + , k − U yi , j , k a i , j + / , k + U yi , j − , k − U yi , j , k a i , j − / , k + U zi , j , k + − U zi , j , k a i , j , k + / + U zi , j , k − − U zi , j , k a i , j , k − / , where a i + / , j , k = ( a i , j , k + a i + , j , k ) / + (cid:15) , etc, and a small o ff set (cid:15) = − is added to avoid division by zero (we have verified itmakes no e ff ect on the actual solution). Fluxes U j are defined as U xi , j , k = E i , j , k h xxi , j , k + (cid:16) E i , j + , k h xyi , j + , k + E i , j − , k h xyi , j − , k + E i , j , k + h xzi , j , k + + E i , j , k − h xzi , j , k − (cid:17) , U yi , j , k = E i , j , k h yyi , j , k + (cid:16) E i + , j , k h yxi + , j , k + E i − , j , k h yxi − , j , k + E i , j , k + h yzi , j , k + + E i , j , k − h yzi , j , k − (cid:17) , U zi , j , k = E i , j , k h zzi , j , k + (cid:16) E i + , j , k h zxi + , j , k + E i − , j , k h zxi − , j , k + E i , j + , k h zyi , j + , k + E i , j − , k h zyi , j − , k (cid:17) , etc.Let us define two new, iteration-independent, discretized quantities, A I ≡ γ + γ (cid:0) β a I − w I , I (cid:1) ( w I , I < D I [ E ; a ]) and B I ≡ s I ∆ x − (1 − β ) a I E (0) I , where α , β , and γ are constants.Then one iteration of the OTVET elliptic solver consists of computing d ( n ) I = ˆ D I [ E n ; a ] − β a I E ni + B I , for all I , followed by updating E I as E ( n + I = E ( n ) I + α A I d ( n ) I . Numerical stability requires α <
1. A set of values that works particularly well is α = . ,β = . ,γ = . The left panel of Figure 7 shows the error in the propagation of the ionization front as a function of time in the Test ffi cient to achieve 2% precision in the ionization front evolution, and mere 10 iterations2 F ig . 7.— Left: Error in the position of the ionization front r I as a function of time in the Iliev et al. (2006a) Test ffi cient to achieve the convergence. Right: Comparison of the old and new OTVET implementations in Test1 of Iliev et al. (2006a). In the new implementation the I-front is sharper and it location is more accurate. This images are to be compared with Figure 6 of Ilievet al. (2006a).
22 20 18 16 14 12 10 M -10 -9 -8 -7 -6 -5 -4 -3 -2 -1 − z φ ( M ) [ M p c − m ag − ] z ≈ z ≈ z ≈ z ≈ z ≈ M =9 × M fl M =7 × M fl M h [M fl ] -9 -8 -7 -6 -5 -4 -3 -2 -1 − z ˙ ρ ∗ ( > M h ) [ M fl / y r / M p c ] M → M =9 × M fl M =7 × M fl M =5 × M fl z F ig . 8.— Left: Galaxy luminosity functions for runs with variable mass resolution: B20.sf1.uv2.bw10 (dash-dotted lines, total equivalent particle mass M = × M (cid:12) ) and B20HR.sf1.uv2.bw10 (dashed lines, total equivalent particle mass M = × M (cid:12) ). Right: cumulative star formation rate density in halos vsthe lowest halo mass for the same two runs, plus the lowest resolution run B20LR.sf1.uv2.bw10 (dotted lines, total equivalent particle mass M = × M (cid:12) ).Solid lines show the extrapolation of these 3 simulations to the limit M →
0. The insert plots the degree of convergence in the total star formation rate forB20.sf1.uv2.bw10 and B20HR.sf1.uv2.bw10 runs. already give 5% precision. The right panel of the same figure shows the improvement in tracking the ionization front with thenew scheme. It should be compared with Figure 6 of Iliev et al. (2006a) - now the quality of OTVET solution approaches that ofthe best ray-tracing schemes.
MASS CONVERGENCE TEST
In order to test mass convergence in our simulations, we compare 3 20 h − Mpc runs with di ff erent particle counts (and, hence,mass resolutions): the lowest resolution B20LR.sf1.uv2.bw10 (256 particles, total equivalent particle mass M = × M (cid:12) ),the fiducial B20.sf1.uv2.bw10 (512 particles, total equivalent particle mass M = × M (cid:12) ), and the highest resolutionB20HR.sf1.uv2.bw10 (1024 particles, total equivalent particle mass M = × M (cid:12) ). Figure 8 shows the luminosity functionsat z = −
10 (except for B20LR.sf1.uv2.bw10) on the left panel and the cumulative star formation rate density on the right panel.In numerical analysis it is customary to extrapolate the results at finite resolution in order to estimate the convergence. Weextrapolate the cumulative star formation rate density ˙ ρ ∗ ( < M h | M ) as a function of the simulation resolution M to the limit3 M → ρ ∗ ( < M h | M ) = ˙ ρ ∗ ( < M h |
0) exp( − α ( M h ) M β ( M h )) (D1)and fitting three parameters ˙ ρ ∗ ( < M h | α ( M h ), and β ( M h ) from the three simulation values for M = × M (cid:12) , 7 × M (cid:12) ,and 5 × M (cid:12) in each bin of halo mass M h . Admittedly, the adopted functional form is somewhat arbitrary; we checked severaldi ff erent functional forms and found this one to provide the best fits to the simulation results. It can also be justified by noticingthat extrapolation is often more precise in logarithmic space, and ansatz (D1) is a simple power-law extrapolation in the log space.With such extrapolation, we can estimate the completeness of our simulations for the total star formation rate as a function ofredshift. The highest resolution B20HR.sf1.uv2.bw10 accounts for 90% of all star formation at z = z =
10. Ourfiducial runs include 70% of all star formation at z = z =
10. For the total, integrated overtime production of ionizing radiation, the highest resolution run B20HR.sf1.uv2.bw10 misses 20% of all ionizing photons, whilethe fiducial runs miss 45%.
REFERENCESAgertz, O., Kravtsov, A. V., Leitner, S. N., & Gnedin, N. Y. 2013a, ApJ,770, 25—. 2013b, ApJ, 770, 25Agertz, O., Teyssier, R., & Moore, B. 2011, MNRAS, 410, 1391Ahn, K., Iliev, I. T., Shapiro, P. R., Mellema, G., Koda, J., & Mao, Y. 2012,ApJ, 756, L16Alvarez, M. A. & Abel, T. 2007, MNRAS, 380, L30—. 2012, ApJ, 747, 126Aubert, D. & Teyssier, R. 2010, ApJ, 724, 244Battaglia, N., Trac, H., Cen, R., & Loeb, A. 2013, ApJ, 776, 81Becker, R. H., Fan, X., White, R. L., Strauss, M. A., Narayanan, V. K.,Lupton, R. H., Gunn, J. E., Annis, J., Bahcall, N. A., Brinkmann, J.,Connolly, A. J., Csabai, I., Czarapata, P. C., Doi, M., Heckman, T. M.,Hennessy, G. S., Ivezi´c, ˇZ., Knapp, G. R., Lamb, D. Q., McKay, T. A.,Munn, J. A., Nash, T., Nichol, R., Pier, J. R., Richards, G. T., Schneider,D. P., Stoughton, C., Szalay, A. S., Thakar, A. R., & York, D. G. 2001,AJ, 122, 2850Bigiel, F., Leroy, A., Walter, F., Brinks, E., de Blok, W. J. G., Madore, B., &Thornley, M. D. 2008, AJ, 136, 2846Bigiel, F., Leroy, A. K., Walter, F., Brinks, E., de Blok, W. J. G., Kramer, C.,Rix, H. W., Schruba, A., Schuster, K., Usero, A., & Wiesemeyer, H. W.2011, ApJ, 730, L13 + Bolatto, A. D., Leroy, A. K., Jameson, K., Ostriker, E., Gordon, K., Lawton,B., Stanimirovi´c, S., Israel, F. P., Madden, S. C., Hony, S., Sandstrom,K. M., Bot, C., Rubio, M., Winkler, P. F., Roman-Duval, J., van Loon,J. T., Oliveira, J. M., & Indebetouw, R. 2011, ApJ, 741, 12Bouwens, R. J., Illingworth, G. D., Franx, M., Chary, R.-R., Meurer, G. R.,Conselice, C. J., Ford, H., Giavalisco, M., & van Dokkum, P. 2009, ApJ,705, 936Bouwens, R. J., Illingworth, G. D., Franx, M., & Ford, H. 2007, ApJ, 670,928Bouwens, R. J., Illingworth, G. D., Oesch, P. A., Labb´e, I., Trenti, M., vanDokkum, P., Franx, M., Stiavelli, M., Carollo, C. M., Magee, D., &Gonzalez, V. 2011, ApJ, 737, 90Bouwens, R. J., Illingworth, G. D., Oesch, P. A., Labbe, I., van Dokkum,P. G., Trenti, M., Franx, M., Smit, R., Gonzalez, V., & Magee, D. 2013,ArXiv e-printsBowler, R. A. A., Dunlop, J. S., McLure, R. J., Rogers, A. B., McCracken,H. J., Milvang-Jensen, B., Furusawa, H., Fynbo, J. P. U., Taniguchi, Y.,Afonso, J., Bremer, M. N., & Le Fevre, O. 2013, ArXiv e-printsBradley, L. D., Trenti, M., Oesch, P. A., Stiavelli, M., Treu, T., Bouwens,R. J., Shull, J. M., Holwerda, B. W., & Pirzkal, N. 2012, ApJ, 760, 108Brook, C. B., Stinson, G. S., Gibson, B. K., Kawata, D., House, E. L.,Miranda, M. S., Macci`o, A. V., Pilkington, K., Roˇskar, R., Wadsley, J., &Quinn, T. R. 2012, MNRAS, 426, 690Choudhury, T. R. & Ferrara, A. 2005, MNRAS, 361, 577—. 2006, MNRAS, 371, L55Conroy, C. & Gunn, J. E. 2010, ApJ, 712, 833Conroy, C., White, M., & Gunn, J. E. 2010, ApJ, 708, 58Daddi, E., Bournaud, F., Walter, F., Dannerbauer, H., Carilli, C. L.,Dickinson, M., Elbaz, D., Morrison, G. E., Riechers, D., Onodera, M.,Salmi, F., Krips, M., & Stern, D. 2010, ApJ, 713, 686Fan, X., Carilli, C. L., & Keating, B. 2006a, ARA&A, 44, 415Fan, X., Strauss, M. A., Becker, R. H., White, R. L., Gunn, J. E., Knapp,G. R., Richards, G. T., Schneider, D. P., Brinkmann, J., & Fukugita, M.2006b, AJ, 132, 117Friedrich, M. M., Mellema, G., Alvarez, M. A., Shapiro, P. R., & Iliev, I. T.2011, MNRAS, 413, 1353 Furlanetto, S. R., Hernquist, L., & Zaldarriaga, M. 2004, MNRAS, 354, 695Furlanetto, S. R., McQuinn, M., & Hernquist, L. 2006, MNRAS, 365, 115Furlanetto, S. R. & Oh, S. P. 2005, MNRAS, 363, 1031Genzel, R., Tacconi, L. J., Gracia-Carpio, J., Sternberg, A., Cooper, M. C.,Shapiro, K., Bolatto, A., Bouch´e, N., Bournaud, F., Burkert, A., Combes,F., Comerford, J., Cox, P., Davis, M., Schreiber, N. M. F., Garcia-Burillo,S., Lutz, D., Naab, T., Neri, R., Omont, A., Shapley, A., & Weiner, B.2010, MNRAS, 407, 2091Gnedin, N. Y. 2000, ApJ, 535, 530—. 2004, ApJ, 610, 9Gnedin, N. Y. & Abel, T. 2001, New Astronomy, 6, 437Gnedin, N. Y. & Fan, X. 2006, ApJ, 648, 1Gnedin, N. Y. & Hollon, N. 2012, ApJS, 202, 13Gnedin, N. Y. & Kravtsov, A. V. 2011, ApJ, 728, 88Gnedin, N. Y., Kravtsov, A. V., & Chen, H.-W. 2008, ApJ, 672, 765Gnedin, N. Y., Kravtsov, A. V., & Rudd, D. H. 2011, ApJS, 194, 46Gnedin, N. Y. & Ostriker, J. P. 1997, ApJ, 486, 581Governato, F., Brook, C., Mayer, L., Brooks, A., Rhee, G., Wadsley, J.,Jonsson, P., Willman, B., Stinson, G., Quinn, T., & Madau, P. 2010,Nature, 463, 203Hopkins, P. F., Keres, D., Onorbe, J., Faucher-Giguere, C.-A., Quataert, E.,Murray, N., & Bullock, J. S. 2013, ArXiv e-printsHopkins, P. F., Quataert, E., & Murray, N. 2011, MNRAS, 417, 950Hopkins, P. F., Richards, G. T., & Hernquist, L. 2007, ApJ, 654, 731Iliev, I. T., Ciardi, B., Alvarez, M. A., Maselli, A., Ferrara, A., Gnedin,N. Y., Mellema, G., Nakamoto, T., Norman, M. L., Razoumov, A. O.,Rijkhorst, E.-J., Ritzerveld, J., Shapiro, P. R., Susa, H., Umemura, M., &Whalen, D. J. 2006a, MNRAS, 371, 1057Iliev, I. T., Mellema, G., Pen, U.-L., Merz, H., Shapiro, P. R., & Alvarez,M. A. 2006b, MNRAS, 369, 1625Iliev, I. T., Pen, U.-L., McDonald, P., Shapiro, P. R., Mellema, G., &Alvarez, M. A. 2009a, Ap&SS, 320, 39Iliev, I. T., Whalen, D., Mellema, G., Ahn, K., Baek, S., Gnedin, N. Y.,Kravtsov, A. V., Norman, M., Raicevic, M., Reynolds, D. R., Sato, D.,Shapiro, P. R., Semelin, B., Smidt, J., Susa, H., Theuns, T., & Umemura,M. 2009b, MNRAS, 400, 1283Kaurov, A. A. & Gnedin, N. Y. 2013, ApJ, 771, 35Khokhlov, A. M. 1998, JCPh, 143, 519Kravtsov, A. V. 1999, PhD thesis, New Mexico State University—. 2003, ApJ, 590, L1—. 2009, ArXiv:0906.3295Kravtsov, A. V., Klypin, A., & Ho ff man, Y. 2002, ApJ, 571, 563Kuhlen, M. & Faucher-Gigu`ere, C.-A. 2012, MNRAS, 423, 862Leitherer, C., Schaerer, D., Goldader, J. D., Delgado, R. M. G., Robert, C.,Kune, D. F., de Mello, D. F., Devost, D., & Heckman, T. M. 1999, ApJS,123, 3Leroy, A. K., Bigiel, F., de Blok, W. J. G., Boissier, S., Bolatto, A., Brinks,E., Madore, B., Munoz-Mateos, J.-C., Murphy, E., Sandstrom, K.,Schruba, A., & Walter, F. 2012, AJ, 144, 3Leroy, A. K., Walter, F., Brinks, E., Bigiel, F., de Blok, W. J. G., Madore, B.,& Thornley, M. D. 2008, AJ, 136, 2782Leroy, A. K., Walter, F., Sandstrom, K., Schruba, A., Munoz-Mateos, J.-C.,Bigiel, F., Bolatto, A., Brinks, E., de Blok, W. J. G., Meidt, S., Rix, H.-W.,Rosolowsky, E., Schinnerer, E., Schuster, K.-F., & Usero, A. 2013, AJ,146, 19McQuinn, M., Lidz, A., Zahn, O., Dutta, S., Hernquist, L., & Zaldarriaga,M. 2007, MNRAS, 377, 1043Mesinger, A. & Furlanetto, S. 2007, ApJ, 669, 663Mesinger, A., Furlanetto, S., & Cen, R. 2011, MNRAS, 411, 955 Mitra, S., Choudhury, T. R., & Ferrara, A. 2011, MNRAS, 413, 1569—. 2012, MNRAS, 419, 1480Muratov, A. L., Gnedin, O. Y., Gnedin, N. Y., & Zemp, M. 2013a, ApJ, 772,106—. 2013b, ApJ, 773, 19Oesch, P. A., Bouwens, R. J., Illingworth, G. D., Gonzalez, V., Trenti, M.,van Dokkum, P. G., Franx, M., Labb´e, I., Carollo, C. M., & Magee, D.2012, ApJ, 759, 135Oesch, P. A., Bouwens, R. J., Illingworth, G. D., Labb´e, I., Franx, M., vanDokkum, P. G., Trenti, M., Stiavelli, M., Gonzalez, V., & Magee, D.2013a, ApJ, 773, 75Oesch, P. A., Bouwens, R. J., Illingworth, G. D., Labbe, I., Smit, R., vanDokkum, P. G., Momcheva, I., Ashby, M. L. N., Fazio, G. G., Huang, J.,Willner, S. P., Gonzalez, V., Magee, D., Trenti, M., Brammer, G. B.,Skelton, R. E., & Spitler, L. R. 2013b, ArXiv e-printsPawlik, A. H., Schaye, J., & van Scherpenzeel, E. 2009, MNRAS, 394, 1812Richards, G. T., Lacy, M., Storrie-Lombardi, L. J., Hall, P. B., Gallagher,S. C., Hines, D. C., Fan, X., Papovich, C., Vanden Berk, D. E., Trammell,G. B., Schneider, D. P., Vestergaard, M., York, D. G., Jester, S., Anderson,S. F., Budav´ari, T., & Szalay, A. S. 2006, ApJS, 166, 470Ricotti, M., Gnedin, N. Y., & Shull, J. M. 2002a, ApJ, 575, 33—. 2002b, ApJ, 575, 49Robertson, B. E., Furlanetto, S. R., Schneider, E., Charlot, S., Ellis, R. S.,Stark, D. P., McLure, R. J., Dunlop, J. S., Koekemoer, A., Schenker,M. A., Ouchi, M., Ono, Y., Curtis-Lake, E., Rogers, A. B., Bowler,R. A. A., & Cirasuolo, M. 2013, ApJ, 768, 71Rudd, D. H., Zentner, A. R., & Kravtsov, A. V. 2008, ApJ, 672, 19Schaye, J., Vecchia, C. D., Booth, C. M., Wiersma, R. P. C., Theuns, T.,Haas, M. R., Bertone, S., Du ff y, A. R., McCarthy, I. G., & van de Voort,F. 2009, MNRAS, 1888Schenker, M. A., Robertson, B. E., Ellis, R. S., Ono, Y., McLure, R. J.,Dunlop, J. S., Koekemoer, A., Bowler, R. A. A., Ouchi, M., Curtis-Lake,E., Rogers, A. B., Schneider, E., Charlot, S., Stark, D. P., Furlanetto,S. R., & Cirasuolo, M. 2013, ApJ, 768, 196 Shapiro, P. R., Iliev, I. T., Mellema, G., Ahn, K., Mao, Y., Friedrich, M.,Datta, K., Park, H., Komatsu, E., Fernandez, E., Koda, J., Bovill, M., &Pen, U.-L. 2012, in American Institute of Physics Conference Series, Vol.1480, American Institute of Physics Conference Series, ed. M. Umemura& K. Omukai, 248–260Shull, J. M. & Venkatesan, A. 2008, ApJ, 685, 1Sobacchi, E. & Mesinger, A. 2014, ArXiv e-printsSongaila, A. & Cowie, L. L. 2010, ApJ, 721, 1448Stinson, G. S., Brook, C., Macci`o, A. V., Wadsley, J., Quinn, T. R., &Couchman, H. M. P. 2013, MNRAS, 428, 129Stinson, G. S., Dalcanton, J. J., Quinn, T., Gogarten, S. M., Kaufmann, T., &Wadsley, J. 2009, MNRAS, 395, 1455Tacconi, L. J., Neri, R., Genzel, R., Combes, F., Bolatto, A., Cooper, M. C.,Wuyts, S., Bournaud, F., Burkert, A., Comerford, J., Cox, P., Davis, M.,F¨orster Schreiber, N. M., Garc´ıa-Burillo, S., Gracia-Carpio, J., Lutz, D.,Naab, T., Newman, S., Omont, A., Saintonge, A., Shapiro Gri ffiffi