Formation of First Galaxies inside Density Peaks and Voids under the Influence of Dark Matter-Baryon Streaming Velocity, I: Initial Condition and Simulation Scheme
aa r X i v : . [ a s t r o - ph . C O ] N ov Draft version November 29, 2018
Typeset using L A TEX twocolumn style in AASTeX62
Formation of First Galaxies inside Density Peaks and Voids under the Influence of Dark Matter-Baryon StreamingVelocity, I: Initial Condition and Simulation Scheme
Kyungjin Ahn and Britton D. Smith Department of Earth Sciences, Chosun University, Gwangju 61452, Korea San Diego Supercomputer Center, University of California at San Diego, San Diego 92093, USA
Submitted to ApJABSTRACTWe present a systematic study of the cosmic variance that existed in the formation of first stars andgalaxies. We focus on the cosmic variance induced by the large-scale density and velocity environmentengraved at the epoch of recombination. The density environment is predominantly determined by thedark-matter overdensity, and the velocity environment by the dark matter-baryon streaming velocity.Toward this end, we introduce a new cosmological initial condition generator BCCOMICS, which solvesthe quasi-linear evolution of small-scale perturbations under the large-scale density and streaming-velocity environment and generates the initial condition for dark matter and baryons, as either particlesor grid data at a specific redshift. We also describe a scheme to simulate the formation of first galaxiesinside density peaks and voids, where a local environment is treated as a separate universe. Theresulting cosmic variance in the minihalo number density and the amount of cooling mass are presentedas an application. Density peaks become a site for enhanced formation of first galaxies, which competewith the negative effect from the dark matter-baryon streaming velocity on structure formation.
Keywords: cosmology: theory — dark ages, reionization, first stars — Galaxy: formation INTRODUCTIONFirst stars and first galaxies form out of the pri-mordial environment which is devoid of any metal,and they are categorized as Population III (Pop III)objects. For such objects to form, at least high den-sity, efficient gas-phase cooling and shielding from thecoolant-dissociating radiation field are necessary. Suchconditions can be first met inside minihalos whichare the first nonlinear objects in the universe. Pio-neering numerical simulations used to find that firststars are very massive ( M ∗ & − M ⊙ ) andform in isolation inside minihalos (Abel et al. 2002;Bromm et al. 2002; Yoshida et al. 2006), but later,higher-resolution simulations found that multiple for-mation of intermediate-mass ( M ∗ . M ⊙ ) starswere likely as well (Turk et al. 2009; Stacy et al. 2010;Greif et al. 2011a). Semi-analytical studies also findthat intermediate- or even small-mass Pop III stars may Corresponding author: Kyungjin [email protected] have formed out of the primordial environment (Omukai2001; Nagakura & Omukai 2005). Different physicalconditions inside minihalos likely lead to a wide spec-trum of the mass of Pop III stars (Hirano et al. 2014,2015). Direct observation of truly metal-free stars isyet to be made in future surveys. Current observa-tions of ultra-iron-poor stars in the Milky Way (e.g.Caffau et al. 2011; Howes et al. 2015; Jacobson et al.2015), due to the existence of other heavy elements,cannot be taken as a convincing evidence of small-massPop III stars.Because first stars and first galaxies form first in-side minihalos, varying physical conditions of minihalos(forming through the hydrogen-molecule cooling) wouldlead to a variation in the outcome of the first star for-mation. In addition to the minihalo-to-minihalo vari-ance (Hirano et al. 2014, 2015), a large-scale variance inthe streaming flow of baryons against dark matter par-ticles was found to impact the formation and evolutionof minihalos (Tseliakhovich & Hirata 2010, TH here-after). This effect can be parameterized by the large-scale streaming velocity V bc ≡ V b − V c , where V b and Ahn & Smith V c are the bulk velocities of baryons and cold dark mat-ter, respectively. In terms of the comoving wavenumber, k ≃ [10 − − is the range strongly affected, andresults in the overall suppression of the matter densityfluctuation in this scale (TH). The coherence length of V bc , below which V bc hardly varies, is set by the baryonacoustic oscillation (BAO) process during the epoch ofrecombination and is in the order of a few comovingMpc.Subsequent studies have focused on both cosmo-logical and astrophysical implications of the impactof V bc , and indeed have found quantitative changesand new physical phenomena. Both the amplitudeand the peak location of the BAO feature, when ob-served through galaxy surveys, will be affected mostlythrough the impact of V bc on the galaxy formation(Dalal et al. 2010; Yoo et al. 2011; Blazek et al. 2016;Lewandowski et al. 2015; Slepian & Eisenstein 2015;Schmidt 2016). A new type of modulation on the 21-cmintensity mapping due to V bc and the correspondingboost of the signal may allow easier high-redshift cos-mology (McQuinn & O’Leary 2012). The minimummass of halos that can host first stars will be increasedfrom that without V bc , or first star formation will bedelayed, due to the advection of gas across the darkmatter potential well, as found by numerical simulations(Greif et al. 2011b; Maio et al. 2011; Stacy et al. 2011;O’Leary & McQuinn 2012). The mismatch of CDM andbaryon overdensities may induce baryon-dominated ob-jects such as globular clusters (Naoz & Narayan 2014)or yield the cosmic variance in the gas content of halos(Fialkov et al. 2012). The increase of the effective Jeansmass may induce the formation of large-mass, direct-collapse black holes (Hirano et al. 2017; Schauer et al.2017; Regan & Downes 2018).Unfortunately, numerical simulations mentionedabove (Greif et al. 2011b; Maio et al. 2011; Stacy et al.2011; O’Leary & McQuinn 2012; Hirano et al. 2017;Schauer et al. 2017), except O’Leary & McQuinn (2012),are likely to underestimate the negative effect of V bc on star formation and other similar effects. This is be-cause these work generated the initial condition basedon a typical Boltzmann code (such as CAMB) andimposed a sudden V bc at some initial time t i . Thisprocedure would then completely underestimate thenegative effect that had existed before t i . Instead, oneshould calculate the cumulative effect of V bc from therecombination epoch to t i and then generate the initialcondition, which would then contain the negative effectthat had existed before t i . O’Leary & McQuinn (2012)followed this track by first calculating the evolutionof perturbations (given in terms of equations by TH) under the influence of V bc and then generating initialconditions, which were used for numerical simulationsof their nonlinear growth. As a result, they provideda new initial condition generator CICsASS (Cosmolog-ical Initial Conditions for AMR and SPH Simulations)that can be used to simulate the nonlinear evolution ofsmall-scale structures under a given V bc environment.In the meantime, improvements on the original for-mulation of TH was made by considering long-rangemodes that had been neglected in TH but found tobe of higher significance than V bc . It was found thatthe large-scale overdensity, or equivalently the velocitydivergence, would impact the small-scale density fluc-tuation more efficiently than V bc (Ahn 2016, Ahn16hereafter; Blazek et al. 2016; Schmidt 2016). Ahn16adopted the original peak-background split scheme ofTH but with additional large-scale mode contributions,approximating long-range (small wavenumber) modesas a uniform local patch with given overdensities (∆ c and ∆ b ), velocity divergences (Θ c and Θ b ), V bc , andthe temperature fluctuation ∆ T . Ahn16 found thatsmall-scale perturbations would evolve in a biased way:the higher ∆ c was, the faster δ would grow. Schmidt(2016) showed that the galaxy bias and subsequentlytheir clustering are more strongly affected by the large-scale differential overdensity (∆ bc ≡ ∆ c − ∆ b , the differ-ence between the CDM overdensity ∆ c and the baryonoverdensity ∆ b ) and the differential velocity divergence(Θ bc ≡ Θ c − Θ b , the difference between the CDM ve-locity divergence Θ c and the baryon velocity divergenceΘ b ) than V bc . Blazek et al. (2016) invoked a similarperturbative approach as in Schmidt (2016) includingall 1-loop contributions, and found that the BAO fea-tures (both the amplitude and the peak location) shouldbe affected.In order to study the growth of small-scale structuresunder the given large-scale environment, characterizedmainly of ∆, Θ and V bc , a most accurate way wouldbe to (1) generate an initial condition at redshift highenough to make perturbative calculations valid and (2)then numerically simulate their evolution to the non-linear regime. For (1), even CICsASS falls short of therequirement: initial conditions are adequately generatedonly for mean-density (∆ = 0) regions in CICsASS be-cause the original formulation of TH is used. Ahn16showed clearly that δ would grow faster under larger ∆environment, and thus the full set of equations treat-ing large-scale impact, given by Ahn16, should be inte-grated in order to generate an initial condition. Thisrequirement gets even stronger if one is interested in avery high-density or a low-density environment, becausethen the error from just using equations of TH instead irst Galaxies in Peaks and Voids (Baryon-CDM COs-Mological Initial Condition generator for Small scales)that fulfills this requirement. For (2), because we wantto cover any overdensity environment, the usual set-upwith zero overdensity and a periodic boundary condi-tion would not work. This instead is allowed by us-ing the well-known trick of treating the environment asa separate universe with local cosmological parameters,which are different from those of the global universe (e.g.Goldberg & Vogeley 2004).This paper is organized as follows. We describe theframework for BCCOMICS in Section 2, including thespatial variance of large-scale environment, how this en-ters the evolution equation of small-scale perturbations,and the schematics for generating initial conditions fromthe transfer functions that become anisotropic. In Sec-tion 3 we describe a strategy to identify an overdense(underdense) patch as a separate universe by defininglocal parameters. In Section 4, as an application of theschemes in Sections 2 and 3, we describe a suite of nu-merical simulations of high-redshift, small-scale struc-ture formation inside varying large-scale overdensity and V bc environments. In Section 5, we conclude the paperwith a summary. Unless noted otherwise, length scalesare expressed in comoving length units. INITIAL CONDITION: BCCOMICSIn order to generate initial conditions for this study,we need an initial condition generator that implementsenvironmental effects from both streaming-velocity anddensity. This requires first solving the perturbationequations for the evolution of dark-matter and baryoncomponents at the least. In the early universe, pertur-bative description of their evolution is justified when thescales of interest have not entered the nonlinear regimeyet.The perturbation theory for the evolution of small-scale structures under the influence of only the large-scale streaming velocity environment was formulated byTH, and there exists an initial condition generator thatimplements this theory (CICsASS, O’Leary & McQuinn2012). In terms of overdensity ( δ ≡ ( ρ − ¯ ρ ) / ¯ ρ , definedwith the local density ρ and the average density ¯ ρ ), pe-culiar velocity ( v , with its divergence θ ≡ (1 /a ) ∇ · v where a is the scale factor and ∇ is the gradient inthe comoving frame) and temperature fluctuation ( δ T ≡ ( T − ¯ T ) / ¯ T , with the local baryon temperature T and themean baryon temperature ¯ T ), they solve the following linearized perturbation equation in Fourier space: ∂δ c ∂t = − θ c ,∂θ c ∂t = − H Ω m ( f c δ c + f b δ b ) − Hθ c ,∂δ b ∂t = − ia − V bc · k δ b − θ b ,∂θ b ∂t = − ia − V bc · k θ b − H Ω m ( f c δ c + f b δ b ) − Hθ b + a − k B ¯ Tµm H k { δ T + δ b } ,∂δ T ∂t = 23 ∂δ b ∂t − x e ( t ) t γ a − ¯ T γ ¯ T δ T , (1)where the subscript c and b denote the cold dark matterand baryons, respectively, V bc ( ≡ V b − V c ) is the large-scale streaming velocity of baryons against dark matter, k is a given wavenumber, H is the Hubble parameter,¯ T γ is the average radiation temperature, k B is the Boltz-mann constant, µ is the mean molecular weight, m H isthe hydrogen mass, t γ = 1 . × years, and x e ( t ) isthe global ionized fraction . As long as one limits thedensity environment to that of the mean-density, CIC-sASS correctly solves equation (1) and provides initialconditions that are accurate to the linear order.However, in order to include the effect of the den-sity environment in addition to that of the streaming-velocity environment, equation (1) is no longer validand a more accurate description is required. In addi-tion to V bc , a given large-scale patch will have in gen-eral non-zero overdensity (∆), velocity divergence (Θ)and temperature fluctuation (∆ T ). This induces cou-pling of large-scale and small-scale fluctuations (denoted“mode-mode coupling” hereafter), and was formulatedby Ahn16 into the following perturbation equation: ∂δ c ∂t = − (1 + ∆ c ) θ c − Θ c δ c ,∂θ c ∂t = − H Ω m ( f c δ c + f b δ b ) − Hθ c ,∂δ b ∂t = − ia − V bc · k δ b − (1 + ∆ b ) θ b − Θ b δ b ,∂θ b ∂t = − ia − V bc · k θ b − H Ω m ( f c δ c + f b δ b ) − Hθ b + a − k B ¯ Tµm H k { (1 + ∆ b ) δ T + (1 + ∆ T ) δ b } ,∂δ T ∂t = 23 (cid:26) ∂δ b ∂t + ∂ ∆ b ∂t ( δ T − δ b ) + ∂δ b ∂t (∆ T − ∆ b ) (cid:27) The original work by TH, for simplicity, ignores δ T and as-sumes a spatially constant sound speed, and works in the baryon-rest frame. Equation (1) instead uses the CDM-rest frame. Ahn & Smith − x e ( t ) t γ a − ¯ T γ ¯ T δ T . (2)Equation (2) is the governing equation for the evolutionof small-scale perturbations ( δ c , δ b , θ c , θ b , δ T ) under theinfluence of local large-scale environment (∆ c , ∆ b , Θ c ,Θ b , ∆ T ), accurate to the leading-order. Equation (2)is developed under the CDM-rest frame, such that anobserver is moving at the bulk velocity of CDM parti-cles, V c , inside the patch. This equation is based on theassumption that the two scales are well separated, andthus this formalism is a peak-background split scheme.We consider the variation of the large-scale fluctuationsat a scale of 4 comoving Mpc, which is a natural choicebecause the streaming-velocity at recombination is co-herent at a few comoving Mpc (TH).We introduce a newly developed, cosmological ini-tial condition generator BCCOMICS (Baryon-Cold darkmatter COsMological Initial Condition generator forSmall-scales) that incorporates the full impact of large-scale environment on small-scale fluctuations. Obvi-ously, we developed BCCOMICS because there has notbeen any initial condition generator that fully incorpo-rates this effect. The large-scale environments are real-ized as a set of three-dimensional fields and their evo-lution is calculated properly (Section 2.1). BCCOMICSthen solves the perturbation equation (Equ. 2) and gen-erates three-dimensional fields of small-scale quantities(Section 2.2). Incorporating non-zero overdensity envi-ronment requires a careful reassignment of cosmologicalparameters and scaling laws, which will be described inSection 3 in detail.2.1. Evolution of large-scale variance
We need to understand how large-scale ( & a fewMpc) fluctuations vary in space, and how they evolvein time until a dedicated epoch for the initial conditionis reached. This is because the impact of the environ-ment is continuous in time as seen in equation (2). Forany given patch chosen out of many large-scale patcheswith varying physical properties, once we are able totrack its evolution, we obtain the full understanding ofa given environment and this can be fed into equation(2). This evolution, if in the linear regime, can be easilyobtained by Boltzmann solvers such as CAMB, whichsolves Eulerian linear evolution equations of all compo-nents in the standard ΛCDM cosmology. Nevertheless,it would be somewhat more efficient computationally(we will justify this below) and also intuitive to have anapproximation based on a simplified evolution equation.Indeed, the evolution of large-scale variables after re-combination can be well approximated by solving for thegrowth of 4 independent modes under a very simplified evolution equation (Ahn16) given by ∂ ∆ + ∂t = − Θ + ,∂ Θ + ∂t = − H Ω m ∆ + − H Θ + ,∂ ∆ − ∂t = − Θ − ,∂ Θ − ∂t = − H Θ − , (3)where ∆ + ≡ f c ∆ c + f b ∆ b , Θ + ≡ f c Θ c + f b Θ b , ∆ − ≡ ∆ c − ∆ b , and Θ − ≡ Θ c − Θ b with f c ≡ ¯ ρ c / (¯ ρ c + ¯ ρ b )and f b ≡ − f c . They are growing, decaying, com-pensated, and streaming modes, which are solutionsto these linear equations (equation 3). The growingand decaying modes compose the adiabatic perturba-tion, and the compensated and streaming modes com-pose the isocurvature perturbation. Note that theisocurvature perturbation can be sourced primordiallyby the cosmic inflation in some inflation models, andsecularly from BAO during the recombination epoch(Tseliakhovich & Hirata 2010; Barkana & Loeb 2011).We only consider the latter possibility, which is self-consistently calculated by CAMB.Equation (3) is in this simple form because (1) fluctu-ation in the radiation component quickly decays in timeafter recombination, (2) the large-scale ( & a few Mpc)variance is likely to be free from the mode-mode cou-pling and (3) the pressure term a − k B ¯ Tµm H k { ∆ T + ∆ b } is negligible due to smallness of k . Then the temporalevolution of these variables for a given patch is given bya linear combination of these modes. Finally, we evolve∆ T passively, after solving for the evolution of { ∆ c , ∆ b ,Θ c , Θ b } beforehand, by integrating the following equa-tion (identical to equation 8 of Naoz & Barkana 2005): ∂ ∆ T ∂t = 23 ∂ ∆ b ∂t + x e ( t ) t γ a − (cid:26) ∆ γ (cid:18) T γ T − (cid:19) − ∆ T ¯ T γ ¯ T (cid:27) , (4)where we use the output transfer function from CAMBfor the radiation fluctuation ∆ γ ( z ), and the initial tem-perature fluctuation ∆ T ( z i ) is fixed by the scheme byNaoz & Barkana (2005):∆ T ( z i ) = ∆ T γ (cid:18) − T ¯ T γ (cid:19) + t γ x e a (cid:18) ∂ ∆ b ∂t − ∂ ∆ T γ ∂t (cid:19) , (5)where all the terms of the right hand side are evaluatedat z i . Note that equations (3-5) can be all solved in k -space, and an r -space map can be obtained from k -space quantity by randomization (Equ. 12 of Ahn16)and Fourier transformation.We note that it is important to consider all the 4modes in evolving equation (3), because the motion of irst Galaxies in Peaks and Voids r ≡ ρ r /ρ crit ) at least in the global evolu-tion of Ω m , which makes Ω m = 1 for quite long afterrecombination. Cosmological N-body+hydrodynamicscodes such as Enzo and Gadget used not to include theradiation component. For Enzo at least, therefore, wehave now implemented non-zero Ω r and Enzo correctlycalculates the corresponding Ω m ( z ) and the Hubble pa-rameter H ( z ) . Otherwise, a spurious effect in structureformation simulation will occur if one applies the trans-fer function from CAMB, which of course considers theradiation component, to N-body+hydro codes which donot have the radiation component.In practice, the following steps are performed for cal-culating the evolution of a patch. We refer readers tosection 2.2 and Appendix of Ahn16 for details. First,at z re , we use the transfer function from a widely-usedBoltzmann solver CAMB. By convolving the transferfunction with a Gaussian random seed, 3D maps of { ∆ c ,∆ b , Θ c , Θ b , V c , V b } , or { ∆ + , ∆ − , Θ + , Θ − , V c , V b } at z re can be realized. Second, we numerically solve for thetemporal evolutions of { ∆ + , ∆ − , Θ + , Θ − , V bc } . Evo-lution of each of these mode variables is fully describedby a single corresponding numerical solution (let’s say F ( z )), because equation (3) is linear. Third, at any z , one simply multiplies the fluctuation value at z re to F ( z ) and obtain the evolved value at z . For example, if F ∆ + ( z ) is the solution of ∆ + with normalization conven-tion F ∆ + ( z re ) = 1, ∆ + ( z ) = ∆ + ( z re ) F ∆ + ( z ). Fourth,as for ∆ T , one can use either the actual evolution (ob-tained by integrating equation 4) or a fitting formulagiven by equation (30) of Ahn16: thanks to the quickcoupling of ∆ T to ∆ b after z ≃
500 the fitting formulais given in terms of ∆ b , and smallness of ∆ T ( . − )before z ≃
500 allows us to simply set ∆ T = 0 when z & The up-to-date development version of Enzo, which is down-loadable from http://enzo-project.org, reflects this implementa-tion.
An important physical intuition can be seen by com-paring spatial maps and correlations of these variablesat z ≃ z re and a target redshift, e.g., z i = 200. Figure 1and Figure 2 show maps and histograms at z = 1000and z = 200, respectively. As seen in Figure 1(a), V cb ≡ V c − V b (arrows) is dominated by V c andthus show convergence into peaks and divergence fromvoids. Comparing Figure 1(a) and 1(b), we find that at z = 1000 ∆ c dominates over ∆ b in amplitude, and thesetwo are very weakly correlated, as also seen in Fig. 1(e).This is due to the tight coupling of baryons to photonsduring recombination, which can also be seen in the mapof ∆ T (Fig. 1c) showing the characteristic feature of thebaryon acoustic oscillation (BAO). CDM moves mainlythrough its self-gravity, producing a very tight correla-tion between ∆ c and Θ c (Fig. 1d). V bc and ∆ c are notcorrelated at all (Fig. 1f).The main difference between the two epochs lies in thebehavior of baryon fluctuations. At z = 200, aside fromthe amplitude, maps of ∆ c , ∆ b and ∆ T are almost indis-tinguishable (Figs. 2a-2c). This is because the baryonfluctuation is now governed predominantly by gravityrather than by photon fluctuation. The correlation be-tween ∆ c and Θ c is even tighter (Fig. 2d) than it is at z = 1000, and now ∆ b is strongly correlated with ∆ c (Fig. 2e). The latter fact is in stark contrast with theloose correlation seen at z = 1000 (Fig. 1e). Becauseof this, one may assume that baryons are “locked” intoCDM such that ρ b /ρ c = Ω b / Ω c in all patches, and con-sider ∆ c as the only important density environment after z = 200. This is indeed an approximation appropriatefor the study of first-galaxy formation, and allows oneto use a periodic boundary condition without worryingthe net inflow of baryons to the simulation box. Wetherefore take this approximation and isolate the simu-lation box with baryon and CDM contents fixed in thispaper. However, we note that the full degree of sucha “locking” has not happened yet, because ∆ b / ∆ c = 1but ∆ b / ∆ c ≃ .
5. This makes ρ b /ρ c slightly off fromthe cosmic abundance Ω b / Ω c , even though the small-ness of ∆ b and ∆ c at z = 200 make this approximationacceptable. Nevertheless, in time, the value of ∆ b / ∆ c gradually approaches 1, but still with some variance overthe large scale ( k . . − ) (Fig. 3). This is indeeda result of BAO. For a precision cosmology with galax-ies or the intergalactic medium (IGM), for example, oneshould therefore consider other large-scale variables aswell as ∆ c . V bc and ∆ c still remain uncorrelated (Fig.2f).2.2. Anisotropic transfer function and real-spacefluctuations
Ahn & Smith (a) (b) (c)(d) (e) (f)
Figure 1. (a) Map of CDM overdensity ∆ c (colored cells) and the streaming velocity V cb ≡ V c − V b (arrows) on a slice of200 Mpc are containing 50 cells. This is an arbitrarily chosen part of the periodic box with an actual volume of 604 Mpc .(b) Map of baryon overdensity ∆ b on the same slice. (c) Map of temperature overdensity ∆ T on the same slice. (d) Distributionof ∆ c and the CDM velocity divergence Θ c , with colors representing the number of cells in sampling bins. (e) Distribution of∆ c and ∆ b . (f) Distribution of ∆ c and V bc ≡ | V cb | . All panels are based on quantities at z = 1000. All the transfer functions (let us denote it by T s ( k , z ; V bc , { ∆ L } ) for small-scale component s un-der the influence of large-scale component L ; { ∆ L } is the abbreviation for all large-scale fluctuations) ofsmall-scale perturbations ( { δ s } = { δ c , δ b , θ c , θ b , δ T } )are anisotropic due to the impact of V bc , as expectedfrom equation (2). This makes it necessary to calculate T s ( k , z ; V bc , { ∆ L } ) at every wavenumber k . Or, wecan use an implied dependence T s ( k , z ; V bc , { ∆ L } ) = T s ( k, µ, z ; V bc , { ∆ L } ), where µ ≡ k · V bc / ( kV bc ). In T l depends on several environmental variables, in principle.Nevertheless, strong correlation at z = 1000 between ∆ c and Θ c allows us to remove the dependence on Θ c . In addition, we be-lieve that the dominant effect on small-scale perturbations comesmainly from V bc and ∆ c after the motion of baryons are approxi-mately synchronized with that of CDM. Therefore, in this paperwe do not investigate the dependence on other variables, ∆ b , Θ b and ∆ T , in our simple case studies (Section 4). practice, we first integrate equation (2) and generate T s ( k, µ, z ; V bc , { ∆ L } ) with dozens of logarithmic sam-ples k and a uniformly gridded samples of µ , for agiven patch with a specific set of large-scale fluctua-tions { ∆ L } and V bc . (1) This forms a 2D table of T s in terms of { k, µ m } for a given patch, where sub-script m is the integer index for the discretized valuesof µ , or µ m . (2) Then, for a given k and V bc , weperform a 2D interpolation at the actual set of points { k, µ } and obtain T s ( k , z ; V bc , { ∆ L } ). Once we ob-tain T s ( k , z ; V bc , { ∆ L } ), with proper normalization(to fluctuations at z = 1000), the power spectrum ofquantity l becomes P s ( k , z ; V bc , { ∆ L } ) = P s ( k, z =1000) T s ( k , z ; V bc , { ∆ L } ). (3) We apply two Gaussianrandom seeds G and G to obtain real and imaginaryparts of the k -space field, where both G and G have irst Galaxies in Peaks and Voids (a) (b) (c)(d) (e) (f) Figure 2.
Same as Figure 1, except that all panels are now based on quantities at z = 200. − .
01 1 100 k (Mpc − )0 . . . ∆ b / ∆ c z=10z=50z=200 Figure 3.
Overdensity ratio ∆ b / ∆ c calculated by the linearBoltzmann solver at different wavenumbers k and redshifts z . mean 0 and variance 1:Re( δ s ( k )) = G N (cid:18) P s ( k )2 V box (cid:19) / sign [ T s ( k )] , Im( δ s ( k )) = G N (cid:18) P s ( k )2 V box (cid:19) / sign [ T s ( k )] . (6) (4) Of course, real-space variables are all in real num-bers, thus we use the condition δ ∗ s ( k ) = δ s ( − k ) afterfilling only 1/2 of the allowed k -space ( ∗ denotes thecomplex conjugate). In addition, all monopole terms( k = 0) are assumed zero. Removing monopole termsallows one to use the usual periodic boundary condition,and we will explain in Section 3 how this becomes pos-sible even in the presence of non-zero ∆ c and Θ c . (5)Finally, we take the Fourier transform of δ s ( k ) and ob-tain δ s ( x ). Vectorization of the baryon velocity field isperformed through the relation v ( k ) = − ( ia k /k ) θ ( k )assuming a curl-free velocity field.The whole process (1) - (5) is straightforwardfor generating uniform-grid quantities. For particlequantities (of CDM particles; of baryon particles ifsmoothed particle hydrodynamics is used), we followthe usual Lagrangian perturbation theory (LPT), re-lating the displacement vector to the Eulerian overden-sity (Bouchet et al. 1992). Because we calculate theEulerian overdensity based on the linear-order Boltz-mann solver CAMB, we restrict our calculation to the1st-order Lagrangian perturbation theory (1LPT). A Ahn & Smith particle is displaced from its Lagrangian point q toits Eulerian point x by the displacement vector Ψ , as x = q + Ψ . To linear order, δ ( q ) = −∇ · Ψ ( q ) , (7)which can be again cast into the k -space quantities Ψ = ( i k /k ) δ (8)and allows using T δ ( k ) for δ ( k ). It is easy to gener-ate the real-space displacement field through the Fouriertransformation of Ψ ( k ) if displacing particles from cubi-cally and uniformly spaced positions; in case of a glass,an interpolation of this uniform-space Ψ( q ) field ontoglassy positions is further required. For the particlevelocity field, consequently, its generation first takes thesteps for the uniform-grid data, and interpolation of thisuniform-grid velocity field onto the displaced positionsare performed. In practice, however, such interpolationbecomes unnecessary if Ψ is very small compared to thelength scale ∼ | Ψ / ∇ · Ψ | . HOW TO SIMULATE DENSITY PEAKS ANDVOIDSBecause BCCOMICS solves equation (2) which in-cludes possible couplings between large-scale (in termsof { ∆ l } ) and small-scale modes (in terms of { δ s } ), BC-COMICS provides so far the most accurate (and only)initial condition regarding the large-scale environmen-tal effect on small-scale fluctuations. To take advantageof this fact and simulate structure formation in over-dense and underdense regions, however, we need to takefurther steps than those required for a mean-density en-vironment. In this section, we lay out a strategy forachieving this goal.3.1. Local Hubble parameter and the Friedmannequation
We will take an overdense (underdense) patch as a sep-arate universe. There can be other ways of simulatingan overdense patch, such as taking a much larger volumeas a simulation box and zooming into a patch of interestwith e.g. nested grids (e.g. O’Shea et al. 2015). Never-theless, this scheme has a merit of allowing the periodicboundary condition, because an overdense (underdense)patch in our universe is now treated as a mean-densitypatch in a universe with different cosmological parame-ters.Because the patch will detach from the global Hub-ble flow, almost all relevant cosmological parametersshould be redefined. Anything redefined in such a sep-arate universe will be called “local” and the relevant symbol will be capped by ˜. The mass conservation ofsuch a patch first defines the local Hubble parameter(Goldberg & Vogeley 2004):˜ H = H − ˙∆3(1 + ∆) , (9)where ∆ is the matter overdensity of a given patch, ˙ = ddt with cosmic time t , and we assume that the net influxesof CDM and baryons are both zero on this patch in thelocal viewpoint, or ˜Θ c = ˜Θ b = 0, in order to incor-porate the usual periodic box condition for simulation.Nevertheless, because we take a patch in a CDM-restframe, there will be a bulk flow of baryons with velocity˜ V b = V bc . Note also that ˜∆ c = ˜∆ b = 0 by definition,because we take a patch as a separate mean-density uni-verse.The local Friedmann equation can be written in vari-ous forms: 1 = ˜Ω m + ˜Ω r + ˜Ω Λ + ˜Ω K , (10)˜ H = ˜ H h ˜Ω m, ˜ a − + ˜Ω r, ˜ a − + ˜Ω Λ , + ˜Ω K, ˜ a − i / , (11)˜ H = ˜ H i " ˜Ω m,i (cid:18) ˜ a ˜ a i (cid:19) − + ˜Ω r,i (cid:18) ˜ a ˜ a i (cid:19) − + ˜Ω Λ ,i + ˜Ω K,i (cid:18) ˜ a ˜ a i (cid:19) − / , (12)where ˜ a = 1 convention is used in equation (11), andthe subscript i refers to an initial time with the initialscale factor ˜ a i . Equation (9) gives˜ H i H i = 1 − ( ˙∆) i H i (1 + ∆ i ) . (13)The initial local cosmological parameters become˜Ω m,i = ˜ ρ m,i ˜ ρ crit ,i = (1 + ∆ i ) ρ m,i ˜ ρ crit ,i = (1 + ∆ i ) ρ m,i ρ crit ,i ρ crit ,i ˜ ρ crit ,i = (1 + ∆ i )Ω m,i ˜ H i H i ! − , (14)˜Ω r,i = ˜ ρ r,i ˜ ρ crit ,i = ρ r,i ˜ ρ crit ,i = Ω r,i ˜ H i H i ! − , (15)˜Ω Λ ,i = ˜ ρ Λ ,i ˜ ρ crit ,i = ρ Λ ,i ˜ ρ crit ,i = Ω Λ ,i ˜ H i H i ! − , (16)˜Ω K,i = 1 − (cid:16) ˜Ω m,i + ˜Ω r,i + ˜Ω Λ ,i (cid:17) , (17)where the global (flat ΛCDM universe) Ω i values foreach component is of course given by, e.g.Ω m,i = ρ m,i ρ crit ,i = Ω m, a − i Ω m, a − i + Ω r, a − i + Ω Λ , . (18) irst Galaxies in Peaks and Voids Local cosmological parameters and redshiftmapping
We take Enzo as our model simulation code, and showhow cosmological parameters are assigned when an over-dense (underdense) patch is simulated. Application toother codes will be similar. Because Enzo parameterfile requires those at the “present”, such as
Cosmolo-gyOmegaMatterNow , we need to get the local values ofthese. However, we can arbitrarily define the “present”,only if the local scale factor ˜ a at that time, ˜ a , is nor-malized to 1. This convention is used in Enzo. We notethat the following assignment procedure is provided asa separate function code in BCCOMICS.Let us take some global redshift z ′ and scale factor a ′ = 1 / (1 + z ′ ) to be those of the local present. We thenneed to connect this information to the actual evolutionof ˜ a ( t ). We need to integrate d ˜ a/dt (Equ. 12) d ˜ ad ( tH i ) = ˜ a i ˜ H i H i " ˜Ω m,i (cid:18) ˜ a ˜ a i (cid:19) − + ˜Ω r,i (cid:18) ˜ a ˜ a i (cid:19) − + ˜Ω Λ ,i (cid:18) ˜ a ˜ a i (cid:19) + ˜Ω K,i / , (19)and then form a { tH i , ˜ a ( t ) } table. When performingnumerical integration in practice, we start from sometime ( t s ) deep inside the radiation-dominated epoch,with t s H i ≪
1, and then use the analytical expressionduring this epoch˜ a s = ˜ a i s t s H i ˜ H i H i q ˜Ω r,i , (20)where we set ˜ a i = a i “temporarily”. We find that t s H i = 10 − provides an excellent accuracy for post-recombination epoch, and we use ODE45 (Runge-Kutta4th order equivalent) of Matlab and gnu octave, or sim-ply the Simpson’s rule with appropriate time-binning.In either way we can easily obtain an accuracy of lessthan 10 − at all t .We then sample specific global time tH i ’s from thestarting (global) scale factor a i to the final (global)scale factor a ′ , in terms of the global time variable tH i = { t H i , t H i , ..., t N H i } . This becomes the set of N epochs for data output. The corresponding local scalefactor of the patch will be ˜ a = { ˜ a , ˜ a , ..., ˜ a N } , whichcan be obtained from the { tH i , ˜ a ( t ) } table. ˜ a N = a ′ ingeneral. And for convenience, we match ˜ a N to the localpresent. We thus rescale ˜ a to ˜˜ a = n ˜ a ˜ a N , ˜ a ˜ a N , ..., o , andthe corresponding redshifts become ˜ z = { ˜ z , ˜ z , ..., } = (cid:26)(cid:16) ˜ a ˜ a N (cid:17) − − , (cid:16) ˜ a ˜ a N (cid:17) − − , ..., − (cid:27) . This list is as-signed to CosmologyOutputRedshift in Enzo. When necessary, one can always map ˜ z to z through the { tH i , ˜ a ( t ) , a ( t ) } table.We also need to change several other cosmological pa-rameters. To obtain CosmologyOmegaMatterNow , e.g.,we calculate it by˜Ω m, = ˜Ω m,i (cid:18) ˜ a N ˜ a i (cid:19) − , " ˜Ω m,i (cid:18) ˜ a N ˜ a i (cid:19) − + ˜Ω Λ ,i + ˜Ω r,i (cid:18) ˜ a N ˜ a i (cid:19) − + ˜Ω K,i (cid:18) ˜ a N ˜ a i (cid:19) − , (21)and follow similar steps with correct power-law of(˜ a N / ˜ a i ) in the numerator for all other Cosmolo-gyOmega*Now ’s, where *= { CDM , Baryon , Matter , Lambda , Radiation } . Non-zero curvature term for non-zero ∆ m should be taken in carefully, and in case of(recent-version) Enzo, this term is calculated internallyinstead of being accepted as an input parameter. Cos-mologyHubbleConstantNow , which is the present Hubbleconstant in units of 100 km s − Mpc − , becomes˜ h = ˜ H i
100 km / s / Mpc " ˜Ω m,i (cid:18) ˜ a N ˜ a i (cid:19) − + ˜Ω Λ ,i + ˜Ω r,i (cid:18) ˜ a N ˜ a i (cid:19) − + ˜Ω K,i (cid:18) ˜ a N ˜ a i (cid:19) − / . (22)In addition, because the comoving box size is the properlength at “present”, this needs to be reassigned too.First, let us denote the comoving length of the mean-density box by L . Then the proper length of thebox at ˜ a i is L ˜ a i . The proper length of the box at˜ a N , which is the comoving length of the box, is then L ˜ a i (˜ a N / ˜ a i ) = L ˜ a N . This is then assigned to Cosmolo-gyComovingBoxSize .3.3.
Scaling laws and halo identification
Even when peaks and voids are treated as a separateuniverse, an observer there can measure quantities basedon the global properties of the Universe. For example,an observer that estimated the local matter content ˜Ω m through an observation inside a small volume (e.g. onlyinside a 4 Mpc patch) will realize that the global valueΩ m has been only underestimated after enlarging thesurvey volume. Therefore, scaling laws for time andlength are required.The scaling laws are given trivially. If the local lengthand time scales are ˜ L and ˜ t , the global length scales L and t are given by L = ˜ L a ˜ a (23)and T = ˜ T t ˜ t , (24)0 Ahn & Smith respectively, where t and ˜ t are ages of the actual universeand the local universe (patch), respectively.A few obvious but important applications are immi-nent. If one were to restrict the volume of a galaxysurvey to e.g. 200 Mpc centered at the observer, therewould occur a danger of wrongfully measuring the BAOscale (true value of ∼
150 Mpc) as 150 (˜ a/a ) Mpc where˜ a is the local scale factor of the 200 Mpc volume.Similarly, the cosmological wavenumber of fluctuationsshould be scale properly. If one conducted an auto-correlation analysis on a limited-volume survey samplesof galaxies, the local wavenumber ˜ k and the local cor-relation length ˜ l should be mapped to the global valuesas k = ˜ k ˜ aa (25)and l = ˜ l a ˜ a . (26)Equation (25) should be applied to the scaling law sug-gested for the k -space fluctuation and the power spec-trum by Goldberg & Vogeley (2004, Equs. 21 and 26),where they forgot to scale the wavenumber.Halo identification schemes should also be approachedcarefully, which is of our keen interest. Let us restrict thediscussion to one specific halo identification scheme: theFriends-of-Friends (FoF) algorithm. A halo is identifiedif a collection of particles are connected in lengths thatare smaller than the FoF linking length b , in units of themean particle separation. A given linking length b , is animplicit indicator of the mean density of resulting halos h ρ m i : ( h ρ m i ≃ b/ . − ¯ ρ m or δ lin > δ crit ∼ .
67, asin Lacey & Cole 1994). Because the relation h ρ m i ¯ ρ m ≃ (cid:18) b . (cid:19) − (27)roughly holds regardless of cosmology for any simulationbox (Lacey & Cole 1994), this can be interpreted in ourcase as h ˜ ρ m i ¯˜ ρ m ≃ ˜ b . ! − (28)for any local patch, where halos with a local linkinglength ˜ b in units of the local mean particle separationwill have the local overdensity 180 (cid:16) ˜ b/ . (cid:17) − . Now, anyhalo inside our universe is defined in terms of the globalmean density ¯ ρ m to make h ˜ ρ m i / ¯ ρ m = 180( b/ . − , andthus ˜ b = b (cid:16) a ˜ a (cid:17) . (29)Equation (29) is indeed consistent with a simple factthat a single proper linking length, b × (global mean par-ticle separation), should be applied universally if one imagines a very large, mean-density simulation box thatencloses many overdense and underdense patches.Once a local patch is treated as a separate universe,then FoF halos there should be generated using ˜ b givenby equation (29), which requires scaling ˜ b in a time-dependent way if a constant b , e.g. b = 0 .
2, is used.The minimum number of N-body CDM particles for haloidentification need not change, because a universal cri-terion should be used. For example, let us imagine acertain overdense patch with ˜ a = 0 .
05 but a = 0 . tH i = 10. If the “usual” linkinglength is 0.2 of the mean particle separation in the mean-density patch, then the linking length of the overdensepatch should be 0 . × (0 . / .
05) = 0 . δ th = 1 for triggering a certain astrophysical pro-cess, this needs to be scaled too. To have a universaldensity threshold, ρ th = (1 + ˜ δ th )˜ ρ b = (1 + δ th ) ρ b , weneed ˜ δ th = (˜ a/a ) (1 + δ th ) −
1. If a threshold is giveninstead in terms of n th or ρ th because a physical pro-cess of interest is dependent on the proper density, thenthere is no need for scaling. APPLICATIONAs an application, we use BCCOMICS for generat-ing initial conditions and perform a suite of cosmo-logical simulations of structure formation using Enzo(Bryan et al. 2014), sampling a few patches of varying∆ c and V bc . One can consider other large-scale vari-ables (Θ c , ∆ b , Θ b , and ∆ T ) as well. However, strongcorrelation between Θ c and ∆ c alleviates the need forconsidering Θ c , and { ∆ b , Θ b , ∆ T } are of less impor-tance than ∆ c and V bc . The rather quick transition of avery loose correlation between ∆ c and ∆ b at z = 1000to a very strong correlation at z = 200 is one of thereasons for ignoring ∆ b in this work. Nevertheless, oneshould be wary of the separate impact of ∆ b on theever-existing BAO feature in the matter density fluc-tuation ∆ + and a subsequent impact on the small-scalestructure formation, if e.g. some type of cosmology withfirst galaxies is considered (e.g. McQuinn & O’Leary2012). The set of large-scale variables are { ∆ c , V bc ( z =1000) / (km s − ) } = {
0, 0 } , {
0, 26.5 } , { σ ∆ c , 26.5 } , and { σ ∆ c , 38 } , which are sampled over 4-Mpc patches at z = 1000. V bc = 0 case becomes too rare to be realizedin our (604 Mpc) box when ∆ c = 2 σ ∆ c (see Fig 1f).These sampling parameters are listed in Table 1 withcase names.For fair comparison, in generating small-scale fluctu-ations, we apply a single universal Gaussian randomseed to all these cases. Grid quantities are generated irst Galaxies in Peaks and Voids D0V0 D0V2 D2V2 D2V3∆ c σ ∆ c σ ∆ c V bc (km/s) 0 26.5 26.5 38 Table 1.
All cases use initial conditions with a universalrandom seed, a grid resolution and the CDM-particle numberof 512 each, and the initial box size of 1 Mpc. High densitycases D2V2 and D2V3 gradually detaches from the Hubbleexpansion and thus the comoving size of these boxes shrinksin time, in contrast to the cases D0V0 and D0V2. on a uniform grid of 512 cells, and dark matter parti-cle displacements are based on uniform spacing. Eventhough the large-scale quantities are sampled over 4-Mpc patches, we let simulation boxes to be of 1 Mpcto resolve minihalos down to M ∼ M ⊙ . Chemistryof and cooling by primordial elements (H, He and theirions) are calculated, and formation of Pop III stars aretracked by sink particles when the number density cri-terion for baryons, n b > cm − , is met. This is still apilot study, and we do not calculate the radiation trans-fer and its effects on gas, and turn off the mesh refine-ment. The subsequent paper (Paper II, in preparation)will adopt a more aggressive configuration that is suit-able for the study of first star formation.We stress this fact again: it is not a good practice toassume V bc = 0 or impose a sudden V bc (as is done bymany simulations: Greif et al. 2011b; Maio et al. 2011;Stacy et al. 2011; Schauer et al. 2017; Hirano et al. 2017) in small-scale structure formation simulations, whichis one of a few reasons why one should use at least CI-CsASS for ∆ c = ∆ b = 0 case or BCCOMICS for moregeneric cases of non-zero ∆ c and other large-scale quan-tities. First, V bc = 0 lies at the very end of the Maxwell-Boltzmann distribution tail (Equ. 14 of Ahn16) andthus is too a rare event. Second, a sudden impositionof non-zero V bc on an initial condition based on CAMBtransfer functions will underestimate the large-scale en-vironmental effect ( V bc only for CICsASS, and all possi-ble variants for BCCOMICS) on small-scales, which hascontinued since the recombination epoch (see also thediscussion by O’Leary & McQuinn 2012). CAMB andother Boltzmann solvers do not solve equation (1), notto mention equation (2), and thus a sudden impositionof V bc onto an initial condition based on these Boltz-mann solvers suffers from this problem.We show three types of field maps at z=200 in Fig.4, which are CDM overdensity, baryon overdensity andthe ratio of kinetic-to-thermal energies of baryons, un-der varying overdensity and streaming-velocity environ-ments. In this figure, CDM overdensity maps look indis-tinguishable from one another, because the dynamics ofCDM is dominated by the self-gravity of CDM, while baryon overdensity maps show distinguishable smearthat gets stronger as V bc increases. The ratio of kinetic-to-thermal energies of baryons overall increases as V bc increases. If we imagine the practice of a sudden im-position of V bc , the baryon overdensity maps would notshow any mutual difference. The relative importanceof V bc , in terms of energetics, also gets stronger as V bc increases, as seen in figures of the energy ratio.After simulating structure formation based on initialconditions generated by BCCOMICS, we identified ha-los using the FoF scheme. As described in Section 3.3,a universal linking length b = 0 . b = 0 . a/ ˜ a ). We used yt anal-ysis tool for halo identification, and because this tool iskeen only to local values including the local mean par-ticle separation, we fed this time-varying ˜ b ( z ) into yt for any redshift z . Figure 5 is the halo mass function dn/dM for the net DM mass M of FoF halos.A few features are notable. First, overdense patchescontain more halos across the full mass range than mean-density patches, as expected. Second, when overdensityenvironment is the same, higher V bc yields a strongersuppression. Third, the impact of V bc weakens in time(see bottom panels of Fig. 5), which is expected fromthe fact that V bc ∝ /a . These features are consistentwith the positive effect of overdensity and the negativeeffect of V bc on the clumping of baryons inside DM halosand even DM clumping itself, as seen in the k -space vari-ance of matter (CDM+baryon), ∆ m ≡ P m ( k ) k / (2 π ),and that of CDM, ∆ ≡ P c ( k ) k / (2 π ). Fig. 6 showsthat ∆ c boosts both P m ( k ) and P c ( k ) at all values of k ,while V bc suppresses both P m ( k ) and P c ( k ) in a boundregion of k . Aside from very rare and massive halos( M & M ⊙ ) we could identify until z ≃
20, the num-ber density of halos are ∼ [2-4] times as high as that inmean-density patches. Because we allowed gas cooling,it is probable that the negative effect of V bc is some-what reduced compared to the case without cooling (e.g.O’Leary & McQuinn 2012). One subtle feature is thatsuppression of halo formation is not biased toward thelow-mass end. TH predicts that at z ∼
40, halos withmass M ∼ M ⊙ will be more strongly suppressed thanthose with e.g. M ∼ M ⊙ . In Fig. 5, this tendencycan be barely observed. For high-mass end, however,such a comparison is not meaningful due to poor statis-tics. We plan to obtain better statistics in Paper II byenlarging the simulation box size (including larger-scale http://yt-project.org Ahn & Smith
Figure 4.
Overdensity maps of CDM (top) and baryons (middle), and the kinetic-to-thermal energy ratio of baryons (bottom)at z = 200, generated by BCCOMICS. From left to right, the CDM overdensity and streaming-velocity environments are { ∆ c , V bc } = {
0, 0 } (D0V0), {
0, 26.5 km/s } (D0V2), { σ ∆ c , 26.5 km/s } (D2V2), and { σ ∆ c , 38 km/s } (D2V3), where V bc = V bc ( z =1000). Note that directions of V bc are not identical, and neither are directions of the baryon-density smear. modes of fluctuation, equivalently) and thus allowingmore frequent formation of massive halos.The amount of gas that is first gravitationally boundin halos or filaments and then undergoes cooling isclosely related to the star formation process. Becausewe did not allow adaptive mesh refinement (AMR), wedefer our analysis on the star formation to Paper IIbut instead analyzed the total amount of cooling gas.We define the cooling gas by the criterion t cool < t dyn ,where t cool and t dyn are the local cooling time and thelocal dynamical time, respectively. Figure 7 shows the evolution of the total amount of cooling gas ( M cool ) ineach simulation box. We can take this quantity as arough indicator of star formation activity. As nonlin-ear structures grow in time, M cool increases in all cases.The positive effect of overdensity and the negative ef-fect of V bc on gas clumping (or cooling) is also clearlyobserved. In the early phase, at z ∼
30, D0V0 case hasthe largest M cool because the density environments ofD2V2 and D2V3 have not deviated too much from thoseof D0V0 and D0V2, and therefore the hierarchy roughlyfollows that of V bc in descending order: M cool (D0V0) > irst Galaxies in Peaks and Voids M (M ⊙ )02 . . r a t i o − l og (cid:0) d n d M (cid:1) ( M ⊙− M p c − ) D0V0D0V2D2V2D2V3 z=30 M (M ⊙ )0246 r a t i o − l og (cid:0) d n d M (cid:1) ( M ⊙− M p c − ) D0V0D0V2D2V2D2V3 z=20
Figure 5.
Halo mass functions under different { ∆ c , V bc } environments (D0V0: black solid; D0V2: red dotted; D2V2: bluedotted-dashed; D2V3: orange dashed), obtained from N-body+hydro simulations based on initial conditions at z = 200 depictedin Fig. 4. Plotted are halo mass functions (upper) and the ratios with respect to D0V0 case (lower) at z = 30 (left) and z = 20(right).
10 100 1000 10000 k (Mpc − )0 . r a t i o ∆ m D0V0D0V2 ( µ = 0 . µ = 0 . µ = 0 . (a)
10 100 1000 10000 k (Mpc − )11 . r a t i o ∆ c D0V0D0V2 ( µ = 0 . µ = 0 . µ = 0 . (b) Figure 6. (a) Matter power spectrum at z = 200, expressed in terms of the k -space variance ∆ m , for different patches (D0V0:black solid; D0V2: red dotted; D2V2: blue dotted-dashed; D2V3: orage dashed). Except for the D0V0 patch, a constant µ ≡ cos( k , V bc ) = 0 . V bc . The ratio in thebottom is with respect to the D0V0 case. (b) Same as (a), but for the CDM power spectrum. Ahn & Smith
20 25 30z0 . . r a t i o M c oo l ( M ⊙ ) D0V0D0V2D2V2D2V3
Figure 7.
Evolution of the total cooling mass inside thesimulation box (top) and the ratio to the D0V0 case (bot-tom). The line-type convention is the same as Figure 5. Notethat this figure should be considered only for a qualitativecomparison, because the numerical resolution of these sim-ulations, 512 uniform-grid on a simulation box of volume(1 h − Mpc) , is currently not adequate for accurately pre-dicting this quantity. M cool (D0V2) ≃ M cool (D2V2) > M cool (D2V3). Afterthat, however, the positive effect of overdensity takesover as density peaks are more detached from the globalHubble flow, and D2V2 case becomes the site for themost efficient cooling. At z .
20, the positive effect fromdensity dominates over the early-phase negative effectfrom V bc , but still retains the memory of the negative ef-fect from V bc , such that M cool (D2V2) > M cool (D2V3) >M cool (D0V0) > M cool (D0V2). This is indeed an epochwhere minihalos provide a significant contribution tocosmic reionization (e.g. Ahn et al. 2012), and thus ourstudy is expected to improve upon the existing scenariosof reionization.We will extend this application with a more self-consistent treatment in Paper II, especially with theAMR capability on. The quantitative result of this sec-tion, therefore, is likely to be changed. Nevertheless, thequalitative result is expected to remain the same. SUMMARYFirst stars and first galaxies are created inside miniha-los, which are the first nonlinear structure in the historyof the universe. The evolution of small-scale fluctuationsof both CDM and baryons are found strongly affectedby the large-scale environment, leading to cosmic vari- ance in the formation and evolution of these first ob-jects. This effect is caused predominantly by the large-scale density and the baryon-CDM streaming velocity.Because this effect has not been fully incorporated inexisting initial condition generators, we have developedan initial condition generator, BCCOMICS, that fullyincorporates this effect for the first time. BCCOMICSfirst calculates the evolution equations, which were givenby TH for only the mean-density environment and byAhn16 for a fully generic, non-zero overdensity environ-ment, and then generates three-dimensional fields of gridand particle quantities.Study of this cosmic variance requires realizing the lo-cal environment in simulations. This can be realized aszoom-in patches in one big simulation box, or as an in-dividual patch with a periodic box condition. For thelatter, we have developed a systematic scheme to sim-ulate the growth of small-scale structure, inside den-sity peaks and voids, by treating the environment asa separate universe with local cosmological quantities.The affected quantities are the local Hubble parameter,the local cosmic abundance of various contents (CDM,baryon, radiation, cosmological constant), and the lo-cal scale factor. Analysis of simulation data requires ascaling law that maps the local quantity to the globalquantity. For example, a correlation length of galaxyclustering or the spatial BAO peak of the correlationfunction, which are found through a local galaxy sur-vey inside an overdense environment, will differ fromthe global values found through a unlocalized galaxysurvey. We provided a trivial but important scaling lawthat allows one to easily deduce the corresponding globalquantities.As a pilot study, we generated initial conditions byBCCOMICS and performed a suite of N-body+hydrosimulations of small-scale structure formation undervarying large-scale environments. As expected, the over-density environment yields positive feedback effects andthe streaming-velocity environment negative feedbackeffects. Compared to the mean-density environment, ha-los are generated in higher population and gas coolingbecomes more efficient in the overdense environment.The higher the overdensity is, the stronger this positivefeedback will become. In contrast, the streaming veloc-ity tends to suppress halo formation and also the gascooling.We need to improve upon the quantitative predic-tion of this pilot study, and will conduct a more self-consistent study by allowing AMR and possibly the ra-diation transfer as well. At the same time, a wider pa-rameter space of environmental variation will be con- irst Galaxies in Peaks and Voids15sidered, and its result will be described in Paper II (inpreparation).We thank the anonymous referee for a clear report. This work was supported by the NRF grant 2016R1D1A1B04935414and a research grant from Chosun University (2016).REFERENCES