How Density Environment Changes the Influence of the Dark Matter-Baryon Streaming Velocity on the Cosmological Structure Formation
aa r X i v : . [ a s t r o - ph . C O ] A ug Draft version May 9, 2018
Preprint typeset using L A TEX style AASTeX6 v. 1.0
HOW DENSITY ENVIRONMENT CHANGES THE INFLUENCE OF THE DARK MATTER-BARYONSTREAMING VELOCITY ON THE COSMOLOGICAL STRUCTURE FORMATION
Kyungjin Ahn Department of Earth Sciences, Chosun University, Gwangju 61452, Korea [email protected] ABSTRACTWe study the dynamical effect of relative velocities between dark matter and baryonic fluids, whichremained supersonic after the epoch of recombination. The impact of this supersonic motion on theformation of cosmological structures was first formulated by Tseliakhovich & Hirata (2010), in termsof the linear theory of small-scale fluctuations coupled to large-scale, relative velocities in mean-density regions. In their formalism, they limited the large-scale density environment to be thoseof the global mean density. We improve on their formulation by allowing variation in the densityenvironment as well as the relative velocities. This leads to a new type of coupling between large-scale and small-scale modes. We find that the small-scale fluctuation grows in a biased way: fasterin the overdense environment and slower in the underdense environment. We also find that the neteffect on the global power spectrum of the density fluctuation is to boost its overall amplitude fromthe prediction by Tseliakhovich & Hirata (2010). Correspondingly, the conditional mass function ofcosmological halos and the halo bias parameter are both affected in a similar way. The discrepancybetween our prediction and that by Tseliakhovich & Hirata (2010) is significant, and therefore therelated cosmology and high-redshift astrophysics should be revisited. The mathematical formalism ofthis study can be used for generating cosmological initial conditions of small-scale perturbations ingeneric, overdense (underdense) background patches.
Keywords: cosmology: theory — dark ages, reionization, first stars — surveys INTRODUCTIONThe Λ-cold dark matter (ΛCDM) scenario, combinedwith the theory of cosmic inflation, is the successful,concurrent model describing the past and the present ofour universe, consistent with a wide range of observa-tions. In this scenario, cosmological structures grow outof an extremely uniform density field but with tiny fluc-tuations that are seeded by the cosmic inflation. Thegrowth of the CDM density fluctuations and the growthof baryon density fluctuations are not in perfect synchro-nization, because baryons were tightly coupled to pho-tons before the epoch of recombination and thus theirmotion was different from the motion of CDM whichonly reacts to gravity. Only after recombination baryonsgradually decoupled from photons, and followed the mo-tion of the CDM under gravity.Cosmological observations have verified the ΛCDMscenario in scales large enough to make the baryonicphysics almost irrelevant (e.g. Komatsu et al. 2011;Reichardt et al. 2012 Planck Collaboration et al. 2015).However, once in the regime where the baryonic physicsbecomes important, the growth of baryon fluctuations is affected by hydrodynamics and the growth of the CDMis affected by the gravitational feedback from the baryonfluctuations. Some of the usual assumptions that aremade for treating very large scales, therefore, shouldbe taken carefully or modified when treating relativelysmall scales. For example, Naoz & Barkana (2005) im-proved on the previous estimation of the linear densitypower spectrum in small scales, by replacing the usualassumption made in cosmology that the sound speedof baryons is uniform in space with the fact that thesound speed fluctuates in space in small scales. Theyshowed that more than ∼
10 % change occurs in thebaryon density power spectrum and even more changein the baryon temperature power spectrum. A sheerinclusion of the sub-dominant, yet non-negligible bary-onic component in the analysis changes the predictionon the matter density power spectrum at a few percentlevel even in large scales, as was shown in the frameworkof the high-order perturbation theory (Shoji & Komatsu2009; Somogyi & Smith 2010).Similarly, the relative velocity (“streaming velocity”)between baryons and the CDM after the recombina-tion should also be considered carefully in cosmology.Tseliakhovich & Hirata (2010, TH hereafter), for thefirst time, properly calculated the growth of small-scaledensity fluctuations under the influence of the streamingvelocity. Small-scale fluctuations are coupled to large-scale streaming velocity fields which are coherent over afew comoving Mpc scale. This has its own spatial fluc-tuation with ∼
30 km/s standard deviation at the epochof recombination, and then decays in proportion to theinverse of the scale factor a . Its impact on the small-scale structure formation is non-negligible because thestreaming velocity remains supersonic (until the inter-galactic medium is strongly heated). This then leadsto the suppression of small-scale matter density fluc-tuations. The wave-modes that are the most stronglyaffected are around k ∼ / Mpc, and the impact is onthe matter power spectrum, the conditional mass func-tion and the halo bias parameter, to name a few (TH).Subsequent studies have considered the impact of thestreaming velocity in the perspective of both cosmologyand astrophysics. The boost of amplitude and theshift of the peak of the baryonic acoustic oscillation(BAO) feature due to the streaming velocity, in thegas intensity mapping or the galaxy survey, wereintensively investigated (Dalal et al. 2010; Yoo et al.2011; McQuinn & O’Leary 2012; Slepian & Eisenstein2015; Lewandowski et al. 2015; Blazek et al. 2015;Schmidt 2016). They find that both high-redshift(e.g. Dalal et al. 2010; McQuinn & O’Leary 2012)and low-redshift (e.g. Yoo et al. 2011) surveys willbe affected. The impact of the streaming velocityon BAO may be separated out from the impactof the matter density itself (Slepian & Eisenstein2015), which is important because otherwise it willbecome another nuisance parameter in cosmology. Theastrophysical impact of the streaming velocity hasbeen investigated with focus on the formation of thenonlinear structure such as cosmological halos andstellar objects (Maio et al. 2011; Stacy et al. 2011;Greif et al. 2011; Tseliakhovich et al. 2011; Naoz et al.2012; Fialkov et al. 2012; O’Leary & McQuinn2012; Bovy & Dvorkin 2013; Richardson et al. 2013;Tanaka et al. 2013; Naoz & Narayan 2014; Popa et al.2015; Asaba et al. 2016). Most studies indicate thatthe formation of minihalos (roughly in the massrange M = [10 − ] M ⊙ ) and the formation ofstellar objects in them are suppressed. It may inducebaryon-dominated objects such as globular clusters(Naoz & Narayan 2014; Popa et al. 2015), but theactual star formation process leading to globularclusters has yet to be simulated. It may be responsibleeven for the generation of the primordial magneticfield (Naoz & Narayan 2013). Because minihalos arethe most strongly affected among cosmological halos and they are responsible for the early phase of thecosmic reionization process, how they change thehigh-redshift 21-cm background is also of a primeinterest (Visbal et al. 2012; McQuinn & O’Leary 2012;Fialkov et al. 2013). It is noteworthy that some ofthese numerical simulation results (Maio et al. 2011;Stacy et al. 2011; Greif et al. 2011), which are basedon the initial condition generated by the usual Boltz-mann solver such as the Code for Anisotropies in theMicrowave Background (CAMB: Lewis et al. 2000),need to be re-examined. This is because the impact ofthe streaming velocity is cumulative and inherent evenat z ∼
200 (McQuinn & O’Leary 2012), at which orlater these simulations start with a streaming velocityimplemented by hand.The original formalism by TH has been re-investigatedin terms of the high-order perturbation theory in thewave-number space (“ k -space” henceforth), and some“missing terms” previously neglected were found impor-tant (Blazek et al. 2015; Schmidt 2016). Basically, forthe large-scale modes responsible for the streaming ve-locity ( k ∼ [0 . − / Mpc), TH used a trivial solutionfor the evolution of the streaming velocity ( V bc ) andthe density (∆ c and ∆ b being overdensities of the CDMand baryons in large scale, respectively) environment: V bc ∝ a − , ∆ c = ∆ b = 0. Treating this as a new 0th-order solution to the perturbation equations, they thenexamined how the perturbation in small scales grows.In doing so, they treated V bc as a spatial quantity andperturbation variables in small scales as k -space quan-tities. This is basically a high-order perturbation the-ory, coupling large-scale mode ( V bc ) and the small-scalemodes ( δ c and δ b , small scale CDM and baryonic over-densities, respectively). However, V bc is tightly linkedto the fluctuating ∆ c and ∆ b through the density con-tinuity equation, and thus the trivial solution adoptedby TH cannot be used for generically overdense and un-derdense regions. The continuity equation connects thedivergence of V bc to ∆ c and ∆ b , and Schmidt (2016)finds that the divergence of V bc is indeed an importantterm that one should not ignore. Blazek et al. (2015)also works on the generic basis of non-zero ∆ c and ∆ b .We improve on the formalism of TH by also consider-ing the non-zero overdensities. Toward this end, differ-ent from Blazek et al. (2015) and Schmidt (2016), weinherit the original method by TH and focus on theimpact on small-scale modes: large-scale V bc is treatedas the spatial quantity and small-scale overdensities δ c and δ b are treated as the k -space quantity in the per-turbation analysis. Most importantly, we explicitly in-clude generically “non-zero” ∆ c and ∆ b as a new setof spatial quantities, and consequently the divergence ofCDM and baryon velocities as well. We find that thisleads to a set of mode-mode coupling terms, includingthe velocity divergence-density coupling. We carefullyinclude all the coupling terms to the leading order inour perturbation analysis. We also include the bary-onic physics, namely fluctuations in the sound speed(Naoz & Barkana 2005), the gas temperature and thephoton temperature. Our formalism is also suitable forgenerating initial conditions for N-body+hydro numeri-cal simulations. Because the new set of mode-mode cou-plings is imprinted in the initial condition, the initialcondition generator considering the streaming-velocityeffect by O’Leary & McQuinn (2012), CICsASS, shouldalso be improved on if one were to numerically simulatethe structure formation inside overdense or underdenseregions.This paper is organized as follows. After the intro-duction, we lay out the basic formalism and describethe statistics of large-scale fluctuations in Section 2. InSection 3, we show results on the matter density powerspectrum, the conditional halo abundance and the halobias as applications of the formalism. We conclude thiswork in Section 4 with a summary, discussion and futureprospects. Some details left out in the main body aredescribed in Appendices. FORMALISM AND NUMERICAL METHOD2.1.
Fluctuation under non-zero overdensity andrelative velocity: perturbation formalism
We start from a set of equations for perturbations ofrelevant physical variables. All the k modes of interestare in sub-horizon scale such that the Newtonian per-turbation theory holds. Let us define the overdensity δ j ≡ ( ρ j − ¯ ρ j ) / ¯ ρ j , where the subscript j = { c, b } denoteseither the CDM (c) or the baryonic (b) component, and θ j ≡ (1 /a ) ∇ · v j where a is the scale factor and v j is theproper peculiar velocity of component j . When the uni-verse is in the regime where we can ignore fluctuationsof photons and neutrinos due to their rapid diffusion af-ter recombination, we have (e.g. Bernardeau et al. 2002;TH) ∂δ c ∂t = − a − v c · ∇ δ c − a − (1 + δ c ) ∇ · v c ,∂ v c ∂t = − a − ( v c · ∇ ) v c − a − ∇ φ − H v c ,∂δ b ∂t = − a − v b · ∇ δ b − a − (1 + δ b ) ∇ · v b ,∂ v b ∂t = − a − ( v b · ∇ ) v b − a − ∇ φ − H v b − a − c s ∇ δ b , ∇ φ = 4 πGa ¯ ρ m δ m , (1)where a is the scale factor, ∇ is the gradient in the co-moving frame, ¯ ρ m = f c ¯ ρ c + f b ¯ ρ b , δ m = f c δ c + f b δ b , f c = Ω c , / (Ω c , + Ω b , ), f b = Ω b , / (Ω c , + Ω b , ),Ω j, ≡ ¯ ρ j, /ρ crit , is the present-day mean density of component j in the unit of the critical density, c s is thesound speed, and H is the Hubble constant at a givenredshift. Even though a usual approximation for c s is aspatially uniform one given by c s = k B ¯ Tµm H (cid:18) − ∂ log ¯ T∂ log a (cid:19) , (2)which assumes a mean-density environment undergoingHubble expansion with the mean baryon temperature¯ T , for high k modes a more accurate treatment is re-quired (Naoz & Barkana 2005). This requires replacingthe pressure term , c s ∇ δ b k B ¯ Tµm H ∇ ( δ b + δ T + δ b δ T ) , (3)in Equation (1) and considering another rate equation ∂δ T ∂t = 23 ∂δ b ∂t + x e ( t ) t γ a − (cid:26) δ T γ (cid:18) T γ ¯ T − (cid:19) − δ T ¯ T γ ¯ T (cid:27) (4)where δ T and δ T γ are temperature fluctuations of thebaryon and the photon, respectively, ¯ T γ ≡ .
725 K (1 + z ) is the mean photon temperature, x e ( t ) is the globalelectron fraction at time t and t γ ≡ . × yr. For¯ T , we use the fitting formula by TH.Equation (1) allows a trivial solution: δ c = δ b = φ =0, v c = v c ,i ( a/a i ) − and v b = v b ,i ( a/a i ) − , where a i isthe initial scale factor. TH took this as the zeroth-ordersolution of a spatial patch and developed a linear pertur-bation theory of small scale modes inside the patch. Thephysical process is in principle a coupling of small- k andlarge- k modes (mode-mode coupling), which is beyondthe linear theory where all modes are assumed to be mu-tually independent. TH used the fact that the relativevelocity V bc ≡ v b − v c is coherent over the length scaleof a few comoving Mpc (contributed by modes with wavenumbers in the range 0 . . ( k/ Mpc − ) . V bc .Even though a trivial solution exists, there also ex-ists a nontrivial solution to Equation (1), exact to thefirst order. This nontrivial solution is suited to describethe physics inside patches with non-zero overdensity (seealso Blazek et al. 2015). In order to obtain the nontriv-ial solution, we first linearize Equation (1): ∂δ c ∂t = − θ c , We keep the cross term δ b δ T in Equation (3) in order to findthe correct coupling of high- k and low- k modes in Equations (10)and (11). ∂θ c ∂t = − H Ω m ( f c δ c + f b δ b ) − Hθ c ,∂δ b ∂t = − θ b ,∂θ b ∂t = − H Ω m ( f c δ c + f b δ b ) − Hθ b , (5)where Ω m ≡ ¯ ρ m ( a ) /ρ crit ( a ) is the matter content withrespect to the critical density ρ crit at a , and we ig-nored second-order terms and also the pressure term a − c ∇ δ b . This is indeed a valid approximation in thewave number range (0 . . ( k/ Mpc − ) .
1) relevant tothe coherent V bc (TH), where the second-order termsremain much smaller than the first-order terms and thebaryonic sound speed keeps decreasing from c s ∼ ∂δ + ∂t = − θ + ,∂θ + ∂t = − H Ω m δ + − Hθ + , ∂δ − ∂t = − θ − ,∂θ − ∂t = − Hθ − , (6)where δ + ≡ f c δ c + f b δ b , θ + ≡ f c θ c + f b θ b , δ − ≡ δ c − δ b ,and θ − ≡ θ c − θ b . In the matter-dominated (Ω m = 1)flat universe, δ + allows both the growing mode ( δ + ∝ a , θ + ∝ a − / , v + ≡ f c v c + f b v b ∝ a / ) and the decayingmode ( δ + ∝ a − / , θ + ∝ a − , v + ∝ a − ). δ − allows aslowly decaying (“streaming”) mode ( δ − ∝ a − / , θ − ∝ a − , v − ≡ v c − v b ∝ a − ) and a compensated mode( δ − =constant, θ − = 0, v − = 0). During 1000 & z & m = 1, and most of thesimple analytical forms above become no longer intactexcept for { θ − , v − } of the streaming mode and { δ − , θ − , v − } of the compensated mode.Using this mode decomposition, the large-scale per-turbations evolve in the following form:∆ c ( a ) = (cid:8) ∆ gro D g ( a ) + ∆ dec D d ( a ) (cid:9) + f b { ∆ com + ∆ str D s ( a ) } , ∆ b ( a ) = (cid:8) ∆ gro D g ( a ) + ∆ dec D d ( a ) (cid:9) − f c { ∆ com + ∆ str D s ( a ) } , Θ c ( a ) = − aH (cid:26) ∆ gro dD g ( a ) da + ∆ dec dD d ( a ) da (cid:27) + f b (Θ c ,i − Θ b ,i ) (cid:18) aa i (cid:19) − , Θ b ( a ) = − aH (cid:26) ∆ gro dD g ( a ) da + ∆ dec dD d ( a ) da (cid:27) − f c (Θ c ,i − Θ b ,i ) (cid:18) aa i (cid:19) − , V bc ( a ) = ( V c ,i − V b ,i ) (cid:18) aa i (cid:19) − , (7)where we used upper-case letters to denote the “back-ground” fluctuations for each patch of a few comovingMpc, over which we will develop the small-scale pertur-bation. ∆ gro , ∆ dec , ∆ com , and ∆ str are the initial (at z = 1000) values of the growing, decaying, compensated,and streaming modes, respectively. D g ( a ), D d ( a ), and D s ( a ) are growth factors of the growing, decaying, andstreaming modes, respectively, and they are all normal-ized as D g = D d = D s = 1 at z = 1000. We describethe details of these modes in Appendix A. It is notewor-thy, as is well known already, that both CDM and bary-onic components tend to approach the same asymptotes∆ = ∆ gro D g ( a ) and Θ = − aH ∆ gro dD g ( a ) /da , indicat-ing that baryons tend to move together with CDMs intime. More interestingly, V bc decays as a − through-out the evolution at any overdensity environment eventhough V c and V b grow roughly as a / individually(except in regions with ∆ c = ∆ b = 0 where V c and V b decay as a − ). This fact may seem to make the analysis by TH valid in generic overdensity environments to someextent: TH relied on the trivial solution ∆ c = ∆ b = 0,in which all velocity components decay in time such that V c ∝ a − , V b ∝ a − , and most importantly V bc ∝ a − .Because the suppression of the matter-density fluctua-tions ( f c δ c + f b δ b ) depends not on individual velocitycomponents but on V bc only, different temporal behav-ior of individual velocity components among the trivialand generic solutions do not matter. Nevertheless, quan-titative prediction by TH will be questioned in Section3.1, because we use nontrivial solutions (Equation 7)which result in the new type of coupling of high- k andlow- k modes in general. We also require the evolutionof ∆ T , which is given by Equation (4): ∂ ∆ T ∂t = 23 ∂ ∆ b ∂t − x e ( t ) t γ a − ¯ T γ ¯ T ∆ T , (8)which is evolved in conjunction with Equation (1). Herewe neglect ∆ T γ term due to its smallness (see also thefollowing discussion), even though we do not neglect∆ T γ ( a i ) when initializing ∆ T ( a i ) (Section 2.2). We havefound a useful fitting formula for ∆ T ( a ) for patches withvolume (4 Mpc) :∆ T ( a ) = sign(∆ T, A ) dex h α (log ( a ) + 2 . . i | ∆ T, A | . | ∆ T, B | − . , (9)which provides a good fit to ∆ T at 300 & z & T al-together because of smallness of ∆ T in general. Here∆ T, A ≡ ∆ T ( a = 0 .
01) = 0 . b ( a = 0 . α ≡ log (∆ T, B / ∆ T, A ) / . T, B ≡ ∆ T ( a =0 .
1) = 0 . b ( a = 0 . T ( a ) to ∆ b ( a ) at z .
300 yields this simple,empirical relation to ∆ b ( a ). We describe this fitting for-mula in more details in Appendix B. The decoupling of∆ T ( a ) from the CMB is much earlier than the meanvalue experiences, which occur at z ≃ k is, the earlier the decoupling occurs (seee.g. Figure 1 of Naoz & Barkana 2005). Including ∆ T γ explicitly may delay this decoupling to some extent, butwe leave such an accurate calculation to future work.Using the evolution equation for ∆ b ( a ) in Equation (7),Equation (9) is determined solely by local values of the4 modes.Now we expand Equations (1), (3) and (4) to the lin-ear order, taking Equation (7) as the zeroth-order solu-tion. We first define the net density, the net velocity(of the fluid component i ), the net gravitational po-tential and the net baryon temperature as ρ i ( a, x ) =¯ ρ i ( a ) { i ( a, X ) + δ i ( a, x ) } , v i, net ( a, x ) = Ha x + V i ( a, X ) + v i ( a, x ), φ net ( a, x ) = − ( a/ ∂ ( Ha ) /∂t +Φ( a, X ) + φ ( a, x ), and T b , net = ¯ T (1 + ∆ T + δ T ), respec-tively. Here X and x denote the comoving-coordinateposition of the center of a background patch and thatof a small-scale fluid component, respectively, and weuse lower-case letters for small-scale fluctuations. Φ isthe gravitational potential sourced only by the back-ground fluctuations such that ∇ Φ = 4 πGa Ω m ( f c ∆ c + f b ∆ b ) (e.g. Bernardeau et al. 2002). Similarly, ∇ φ =4 πGa Ω m ( f c δ c + f b δ b ). We then have ∂δ c ∂t = − a − V c · ∇ δ c − (1 + ∆ c ) θ c − Θ c δ c ,∂ v c ∂t = − a − ( V c · ∇ ) v c − a − ∇ φ − H v c − (cid:2) a − ( v c · ∇ ) V c (cid:3) ,∂δ b ∂t = − a − V b · ∇ δ b − (1 + ∆ b ) θ b − Θ b δ b ,∂ v b ∂t = − a − ( V b · ∇ ) v b − a − ∇ φ − H v b − a − k B ¯ Tµm H ∇ { (1 + ∆ b ) δ T + (1 + ∆ T ) δ b } − (cid:2) a − ( v b · ∇ ) V b (cid:3) ,∂δ T ∂t = 23 (cid:26) ∂δ b ∂t + ∂ ∆ b ∂t ( δ T − δ b ) + ∂δ b ∂t (∆ T − ∆ b ) (cid:27) + x e ( t ) t γ a − (cid:26)(cid:20)(cid:18) T γ ¯ T − (cid:19) δ T γ (cid:21) − ¯ T γ ¯ T δ T + 4 (cid:2) ∆ T γ δ T + ∆ T δ T γ (cid:3)(cid:27) ∇ φ = 4 πGa Ω m ( f c δ c + f b δ b ) , (10)where we marked the terms that can be further ignoredin square brackets. First, we can safely ignore any termscontaining δ T γ and ∆ T γ , because the former is negligi-ble at z . δ T (see Figure 1) and thelatter is just too small (∆ T γ . − at z = 1000 anddecaying in time) to produce any appreciable impact onthe baryon temperature of high k modes. Secondly, thejustification for ignoring a − ( v i · ∇ ) V i is easily seenin viewpoint of the k -space. With the Fourier expan-sion A ( x ) = P k A ( k ) exp( i k · x ) of a quantity A ( x ) inan actual, real space (“ r space” henceforth), the back-ground velocities have V j ( x ) = P K V j ( K ) exp( i K · x )but with the condition K . − . In contrast, thesmall-scale modes fluctuating against the backgroundpatches have intrinsically larger wave number k , or K ≪ k . Then, at each K , | ( v j · ∇ ) V j ( K ) | ∼ Kv j V j ≪ | ( V j ( K ) · ∇ ) v j | ∼ kv j V j . Similarly, |∇ ∆ b | ∼ K ∆ b ≪|∇ δ b | ∼ kδ b . We finally note that we do not include thecoupling terms between similar wavenumbers in Equa-tion (10), which will involve quadratic and higher-orderpolynomicals of ∆ and δ . Therefore, the validity ofEquation (10) will break down in the non-linear regime.Nevertheless, our approach is more suitable for a crudeestimattion of the conditional halo mass function (Sec-tion 3.2) in terms of the extended Press-Schechter for-malism, which is based on the mapping of the lineardensity growth to the nonlinear growth of the e.g. top-hat density perturbation.The evolution of small-scale perturbations, therefore,is coupled to large-scale perturbations on which theyare sitting. Ignoring the terms in square brackets, andshifting the viewpoint to the CDM rest frame in which Figure 1 . Fluctuations of the CDM density (short-dashed,black), the baryon density (solid, blue), the baryon tempera-ture (dotted, red) and the photon temperature (long-dashed,cyan), represented by the k -space variance at z = 1000. V c = 0 (as in O’Leary & McQuinn 2012; TH chose thebaryon rest frame), Equation (10), in the k -space, finallybecomes ∂δ 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) − x e ( t ) t γ a − ¯ T γ ¯ T δ T , (11)where δ j , θ j and δ T now denote fluctuations in the k -space while ∆ j , Θ j and V j are fluctuations of a givenpatch at ( a, X ) in the r space, given by Equation (7).2.2. Evolution of perturbation inside patches:Numerical Scheme
Evolution of small-scale perturbations can be calcu-lated by integrating the rate equation (Equation 11)from some initial redshift, preferentially not too longafter the recombination epoch when the relative motionhas not yet influenced the evolution. We take z i ≡ r -space quantities whose dis-tributions are all Gaussian, one needs to sample these values in the r -space accordingly. For the small-scalemodes, one just needs to track the evolution of the av-erage value in the k -space.Let us first describe the statistics of backgroundpatches that we expect. TH calculated the evolution ofsmall-scale ( k &
10) fluctuations under different back-ground patches but only of ∆ c = ∆ b = 0, and de-fined “local power spectrum” P loc , m ( k ; V bc ) averagedout over all possible opening angles between V bc and k .In our case, there are extra dimensions to consider whichare ∆ c , ∆ b , Θ c , and Θ b , resulting in a much higher com-putational demand. Fortunately, some of these quan-tities are in perfect correlation with one another withlinear proportionality. In addition, their initial valuesat a = a i completely compose the ensemble at anytime through Equation (7). At the minimal level , itwould suffice to just consider variation of ∆ c in addi-tion to V bc such that the local power spectrum is anexplicit function of the two background quantities, or P loc , m = P loc , m ( k ; V bc ( a i ); ∆ c ( a i )).For the initial condition for background patches, wegenerate 3D maps of ∆ c , ∆ b , Θ c , Θ b , V c , V b , and ∆ T at z = 1000 on 151 uniform grid cells inside a cubicalvolume of V box = (604 Mpc) . We generate fluctuationsof discrete modes that are randomized asRe(∆ k ) = G N (cid:18) P ( k )2 V box (cid:19) / sign [ TF (∆ k )] , Im(∆ k ) = G N (cid:18) P ( k )2 V box (cid:19) / sign [ TF (∆ k )] , (12)for given k , where G and G are random numbersdrawn from mutually independent Gaussian distribu-tions with mean 0 and standard deviation 1, ∆ k standsfor any kind of k -space fluctuations, and TF is thetransfer function of ∆ k , whose sign should be multi-plied because some ∆ k ’s oscillate around zero in k .The configuration is roughly equivalent to applying asmoothing filter of length (604 / j ( k , a i ),and use the continuity equations ( ∂ ∆ j /∂t = − Θ j )for Θ j ( k , a i ), with the help of two CAMB transfer-function outputs at mutually nearby redshifts for timedifferentiation. V j ( k , a i ) is obtained from the rela-tion V j = − ( ia k /k )Θ j . ∆ T ( k , a i ) is fixed by follow-ing the scheme by Naoz & Barkana (2005): we require ∂ ∆ T /∂t = ∂ ∆ T γ /∂t at the initial redshift in Equation For a more accurate treatment or for a specific patch of in-terest, one should also consider variations in other variables, suchthat P loc , m = P loc , m (cid:0) k ; V bc ; ∆ c ; ∆ b ; Θ c ; Θ b ; ∆ T ; ∆ T γ (cid:1) . As thecorrelation between V b (∆ b ) and V c (∆ c ) becomes tighter in time,the initial variation gets gradually diluted, which roughly justifiesour restricting the parameter space only to V bc and ∆ c . (4), which results in∆ T = ∆ T γ (cid:18) − T ¯ T γ (cid:19) + t γ x e ( t i ) a i (cid:18) ∂ ∆ b ∂t − ∂ ∆ T γ ∂t (cid:19) (13)where all quantities are evaluated at a i , especiallywith the help of Equation (7) for ∂ ∆ b ( k , a i ) /∂t and two adjacent CAMB transfer-function outputs for ∂ ∆ T γ ( k , a i ) /∂t = (1 / ∂ ∆ γ ( k , a i ) /∂t . Finally, allthese k -space fluctuations are Fourier-transformed toobtain r -space fluctuations.3D maps and 2D histograms of several initial quanti-ties are presented in Figure 2. Fields of ∆ c and V bc ona part of a slice of the box at z = z i = 1000 are shown inFigure 2(a). As expected, the velocity field converges onoverdense regions and diverges on underdense regions.We find that in most patches V c dominates over V b , andthus the map of V c looks very similar to Figure 2(a).This occurs because baryons lag behind CDMs due totheir coupling to CMB. ∆ b (Figure 2b) is coupled to∆ T (Figure 2c) more strongly than to ∆ c . ∆ c and Θ c (similarly ∆ b and Θ b ) are almost perfectly correlated(Figure 2d). ∆ c and ∆ b are very loosely correlated dueto the tight coupling of baryons to photons at the red-shift (Figure 2e), but the correlation becomes tighter intime. ∆ c and V bc are not correlated (Figure 2f). Be-cause of this fact, the probability distribution function(PDF) P is simply a multiplication of PDFs P ( k ; V bc )and P ( k ; ∆ c ), at the minimal level. Due to Gaussianityat z i , we have P ( k ; ∆ c ) = 1 √ πσ ∆ c exp (cid:20) − ∆ σ c (cid:21) , P ( k ; V bc ) = r π V σ V bc exp " − V σ V bc , (14)where σ ∆ c and σ V bc are the standard deviations of ∆ c and V bc projected onto one Cartesian-coordinate axis,respectively. With our setup, we find that σ ∆ c = 0 . σ V bc = 17 . z i (or the root-mean-square of V bc is √ σ V bc = 30 . c and V bc , the averagepower spectrum P m ( k ) will then be given by the ensem-ble average P m ( k ) = Z ∞ dV bc Z ∞−∞ d ∆ c P ( k ; V bc ) P ( k ; ∆ c ) × P loc , m ( k ; V bc ; ∆ c ) , (15)where the PDFs and integral arguments are theones at z i . Of course, a more accurate andstraightforward way is to just ensemble-average P loc , m = P loc , m (cid:0) k ; V bc ; ∆ c ; ∆ b ; Θ c ; Θ b ; ∆ T ; ∆ T γ (cid:1) overthe patches from a large-box realization, because for ex- ample ∆ c and ∆ b are too poorly correlated at z i .We numerically integrate Equation (11) to examinethe evolution of small-scale (high- k ) fluctuations at anyoverdense (underdense) patch to the linear order, withthe help of Equations (7) and (9) for the evolution ofbackground quantities. In practice, we used the ODE45modules of MATLAB R (cid:13) (2015b, The MathWorks, Inc.,Natick, Massachusetts, United States) and of GNU Oc-tave, which use the 4th-order Runge-Kutta method,with the relative tolerance 10 − and the absolute tol-erance 10 − δ k ( a i ). During the evolution, the number ofintegration steps is the highest for δ T , because its ampli-tude changes from the initial, very small values around δ T γ to final, much larger values close to δ b . Therefore,taking sub-steps for δ T γ while coarser steps for other δ k ’sis expected to boost the computational efficiency, eventhough we did not yet implement the method in ourcomputation. The end result is then ensemble-averagedover varying ∆ c and V bc to obtain P m ( k ) (Equation 15). RESULT3.1.
Power spectrum of the matter density
We first examine how the evolution of P loc , m ( k ; V bc ; ∆ c ) depends on the density en-vironment, and compare the result to the pre-diction by TH. Figure 3 shows the evolution of∆ , m ≡ k P loc , m / π of three arbitrarily chosenwave numbers ( k = {
33, 150, 2000 } /Mpc) when V bc = 22 km / s ( a/a i ) − in different density environ-ments (∆ c ( a i )= { -0.01, -0.005, 0, 0.005, 0.01 } ). Noteagain that σ ∆ c = 0 . z = 1000, and thus thesesamples correspond to ± . σ ∆ c and ± . σ ∆ c . First,as expected, the growth of small-scale fluctuations arebiased when ∆ c > c <
0, withrespect to the mean-density case (prediction by TH).Secondly, when ∆ c ( a i )’s are equal in amplitude butopposite in sign, the deviations of P loc , m ( k ; V bc ; ∆ c )from P loc , m ( k ; V bc ; ∆ c = 0) reveal the same trend butonly until z ≃
850 when | ∆ c ( a i ) | = 0 .
01 and z ≃ | ∆ c ( a i ) | = 0 . | ∆ c ( a i ) | is, the earlier this unbalance starts. Thirdly,the fractional deviation from the mean-density caseis almost universal regardless of the value of k , andthus the timing of the unbalance is approximatelya function only of | ∆ c ( a i ) | . Finally, in some high-∆patches, our linear analysis based on Equation (11)without quadratic and higher-order terms in δ startsto break down at z ∼
20, because these modes enterthe nonlinear regime (∆ , m ( k ) & .
1; see Figure 3)at this epoch. A similar breakdown of the formalismwill occur for perturbations inside very low-∆ patches -0.01-0.00500.0050.01 (a) × -4 -3-2-1012 (b) × -5 -6-4-2024 (c)(d) (e) (f) Figure 2 . (a) 2D map of CDM overdensity ∆ c (colored cells) and relative velocity V cb ≡ − V bc = V c − V b (arrows; projectedon the plane) fields on a slice of 200 Mpc containing 50 cells. The plotted slice is an arbitrarily chosen part of the actualvolume of 604 Mpc we used, containing 151 cells in total. (b) 2D map of baryon overdensity ∆ b on the same slice. (c) 2Dmap of baryon temperature overdensity ∆ T on the same slice. (d) Distribution of ∆ c and CDM velocity divergence Θ c (in unitsof Myr − ). The color bar represents the number of cells in sampling bins. (e) Distribution of ∆ c (x-axis) and ∆ b (y-axis). (f)Distribution of ∆ c and the relative velocity V bc = | V bc | (in units of km/s). All figures use quantities at z = 1000. as well. Therefore, a higher-order scheme than ourwork is required when one is to predict the low-redshiftevolution of δ ’s. On the other hand, if one were togenerate initial conditions for numerical simulationsin the linear regime, our formalism would provide thesufficient accuracy.How large-scale overdensity impacts the evolution ofsmall-scale inhomogeneities is reflected in the densitycontinuity equation. In overdense background patches,∆ j grows in time and Θ j <
0. Then, in Equation (11), − ∆ j θ j and − Θ j δ j work as sources terms in addition to − θ j to the growth of δ j . This will boost the growth rateof δ j of both overdense ( δ j > θ j <
0) and underdense( δ j < θ j >
0) modes. In contrast, in underdensebackground patches, these terms suppress the growthrate of δ j of both overdense and underdense modes. Thedistribution of ∆ j ( a i ) is Gaussian, and therefore for each“bias” case with ∆ j > j <
0. Nevertheless, the unbalance describedabove is expected to boost the average power spectrumfrom that by TH.The overall effect of including the Gaussian distribu-tion of ∆ c is thus to mitigate the negative impact by therelative velocity, predicted by TH, to some extent. In addition, the universality of the unbalance in k boosts P m ( k ) even in the k range (10 . k . / Mpc and k & / Mpc) where the power spectrum is almost un-affected by non-zero V bc (Figure 4 shown in terms of the k -space matter-density variance ∆ m ≡ k P m ( k ) / π ).Note that the result for k . / Mpc cannot be trusted,because our perturbation theory is based on the condi-tion that large-scale modes (0 . . K Mpc .
1) are wellseparated from small-scale modes in scale. Discrepancyof P m ( k ) including non-zero ∆’s from the prediction byTH is negligible at z &
45, but later the discrepancygrows in time. Of course, individual patches may ex-perience discrepancy in P loc , m ( k ; V bc ; ∆ c ) much earlierthan this epoch (Figure 3).3.2. Halo abundance
Understanding the abundance and the spatial distri-bution of cosmological halos is crucial in modern astro-physics and cosmology. In this section, we examine thehalo abundance both in the local and the global sense,just as we did for the matter power spectrum.Let us first revisit the calculation by TH. Theyadopted the extended Press-Schechter formalism andcalculated the local halo abundance, in terms of the con-
Figure 3 . Growth of ∆ , m ≡ k P loc , m ( k ; V bc ; ∆ c ) / (2 π ) with wave-numbers k = {
33, 150, 2000 } Mpc − and V bc ( z =1000) = 22 km/s in initially overdense (red, long-dashed and orange, dot-dashed), mean-density (black, solid), and under-dense (cyan, short-dashed and blue, dotted) regions. Initial CDM overdensities are chosen to be ∆ c ( a i ) = { -0.01, -0.005, 0,0.005, 0.01 } . The mean-density case, or P loc , m ( k ; V bc ; ∆ c = 0), corresponds to P loc , m ( k ; V bc ) by TH. Fractional differences[ P loc , m ( k ; V bc ; ∆ c ) − P loc , m ( k ; V bc ; ∆ c = 0)] /P loc , m ( k ; V bc ; ∆ c = 0) in % are plotted in the bottom sub-panels, with the line-type convention same as in the top sub-panels.
100 1000 10000 k (Mpc − )10 − . . ∆ m V bc = 0T&Hthis workz=199
100 1000 10000 k (Mpc − )0 . . . ∆ m V bc = 0T&Hthis workz=44
100 1000 10000 k (Mpc − )0 . . . . ∆ m V bc = 0T&Hthis workz=19
100 1000 10000 k (Mpc − )0 . . . ∆ m V bc = 0T&Hthis workz=10 Figure 4 . Mean power spectrum P m ( k ) of the matter-density fluctuation, expressed in terms of the k -space variance ∆ m ( k )(defined in the text). Comparison is made for the case without the relative-velocity effect (black, solid), the case investigatedby TH (blue, dotted) and the case investigated in this work (red, short-dashed). k -spacematter-density variance ∆ m ) and V bc will have the num-ber of halos per unit Eulerian comoving volume per M given by dndM ( M | ∆ , V bc ) = r π ¯ ρ m M δ crit − ∆ σ (cid:12)(cid:12)(cid:12)(cid:12) dσdM (cid:12)(cid:12)(cid:12)(cid:12) (1 + ∆) × exp " − ( δ crit − ∆) σ , (16)where δ crit is the critical overdensity of spherical col-lapse, and σ is the variance of density field smoothedwith the window function W M corresponding to mass M , σ ( M, V bc ) = Z ∆ m ( k, V bc ) W M d ln k. (17)One should note that σ ( M, V bc ) should be that ofhigh k modes only, or more accurately a reduced value σ ( M, V bc ) − σ where σ is the variance of den-sity field smoothed with the window function corre-sponding to the mass of the patch (e.g. Bond et al.1991; Mo & White 1996; Ahn et al. 2015). One can in-stead put a lower bound k min (= (cid:2) π ¯ ρ m /M patch (cid:3) / )in the integral of Equation (17), which would beidentical to the reduced variance if a sharp k -spacewindow function is used. The global mass function( dn/dM ) g is simply an average of the local mass func-tion, dn ( M | ∆ , V bc ) /dM , over the ensemble of patches.When dn ( M | ∆ , V bc ) /dM is averaged only over V bc fora given ∆, which is equivalent to visiting only thosepatches with the same ∆ and taking the average, itleads to the conditional mass function ( dn/dM ) ∆ ≡ dn ( M | ∆) /dM .This calculation should be modified, because σ de-pends also on ∆ through the dependence of P loc , m on∆: σ ( M, V bc , ∆) = Z k max k min ∆ m ( k, V bc , ∆) d ln k, (18)which will enter Equation (16). Here we used the sharp k -space filter, and thus k max = (cid:2) π ¯ ρ m /M (cid:3) / . An over-dense patch will then have a boost in dn ( M | ∆ , V bc ) /dM from the value by TH because P m ( k, V bc , ∆ > >P m ( M, V bc , ∆ = 0) and thus σ ( M, V bc , ∆ > >σ ( M, V bc , ∆ = 0), and vice versa (a decrease from the It is not clear whether TH used this reduced variance. Inaddition, a factor of 2 should be multiplied to Equation (18) ofTH. value by TH) for an underdense patch. Obviously, both( dn/dM ) ∆ and ( dn/dM ) g will also be affected.It is important to compare our findings to the usualpeak-background split scheme and the one by TH. In the“standard” scheme, if the density field is purely Gaus-sian, all the wave modes are assumed mutually indepen-dent in the linear regime. Therefore, the local, high- k modes have a universal variance σ k whether or notthey are placed inside a patch with non-zero ∆. Theway how the halo formation is biased in an overdenseregion is simply through the shift in the density (+∆).TH then realized the fact that the variance is not uni-versal but should depend on V bc . Because non-zero V bc tends to suppress σ k , the odds to cross δ crit decreaserelative to the standard picture. We find that thereis another dependency of the variance, which is ∆. Inother words, we find that there are two biasing effectsin an overdense region compared to a mean-density re-gion: getting closer to δ crit because of the shift in thedensity (+∆, also in the standard scheme), and having alarger degree of fluctuation in δ k due to the mode-modecoupling (e.g. source terms − ∆ c θ c and − Θ c δ c in ∂δ c /∂t in Equation 11, which is a new finding). Therefore, bynot fully implementing the effect of non-zero ∆, TH ineffect underestimates and overestimates the halo massfunctions in overdense and underdense regions, respec-tively.We note that this additional bias effect should bepresent even in the standard picture with V bc = 0, be-cause this is due to the natural coupling between thelarge-scale and small-scale density perturbations. In thiscase, however, we are not sure about the quantitativevalidity of the extended Press-Schechter formalism onthe conditional mass function (Equation 16), which isbased on the linear theory guaranteeing Gaussianity atany filtering scales without the mode-mode couplings.Qualitatively, we believe that the boost of local δ k and σ k under ∆ should boost the conditional mass functionto the level estimated by Equations (16) and (18) asdescribed above anyways. We defer a further investiga-tion of this issue, which can be clarified with numericalsimulations of the halo formation under different ∆’s.Discrepancy between the conditional mass functionsby this work and by TH is significant if we focus on in-dividual patches. Figure 5 illustrates how our predictiondiffers from that by TH. For example, at z ∼ − c ( a i ) = 0 .
005 we predict [100 - 2000] % boostin ( dn/dM ) ∆ compared to the values by TH (let us de-note them by ( dn/dM ) ∆ , TH ). For ∆ c ( a i ) = − . dn/dM ) ∆ compared Rigorously speaking, it is not perfectly universal because thelower bound changes slightly in ∆ as M patch = ¯ ρ m (1 + ∆). δ described above. It isalso noteworthy that the discrepancy is the largest forthe rarest halos: first, at any redshift, the discrepancyincreases as the halo mass increases and secondly, forany given halo mass, the discrepancy decreases in time.The discrepancy among the conditional mass func-tions by this work, by TH and by the standard picturealso influences ( dn/dM ) g . Let us just take the exampleof z = 19 (Figure 6). Compared to the standard pre-diction ( V bc = 0), conditional mass functions by THstay lower regardless of ∆. This is due to the sup-pression of structure formation by the relative velocity.In contrast, ( dn/dM ) ∆ in this work is either higher orlower than the standard prediction ( ≡ ( dn/dM ) V bc =0 )depending on the halo mass and ∆. At z = 19, for∆ c ( a i ) = − .
005 ( dn/dM ) ∆ < ( dn/dM ) V bc =0 and for∆ c ( a i ) = 0 .
005 ( dn/dM ) ∆ > ( dn/dM ) V bc =0 for anyhalo mass. The tendency for ( dn/dM ) ∆ to overshoot( dn/dM ) V bc =0 when ∆ > dn/dM ) ∆ among over-dense and underdense regions is increased compared toTH and the standard picture, and the net effect e.g.at z = 19 is to boost ( dn/dM ) g from ( dn/dM ) g , V bc =0 .At much higher redshifts, our ( dn/dM ) g is almost in-distinguishable from that by TH, which undershoots( dn/dM ) g , V bc =0 .We note that ( dn/dM ) ∆ and ( dn/dM ) g have the usualproblem of not correctly predicting the actual massfunction, if one sticks to the original extended Press-Schechter formalism. Minihalos at high redshifts areusually underestimated by the extended Press-Schechterformalism. The usual peak-background split methodsuffers from large discrepancies between its predictionand the N-body simulation results for rare halos in gen-eral. In this case, a hybrid method to connect the peak-background-split halo bias parameter to the better-fitting mean mass function types (e.g. Barkana & Loeb2004; Ahn et al. 2015) is much more appropriate. Wewill apply this method in the future for a better estima-tion of the conditional mass function and the k -spacehalo bias parameter (Section 3.3).3.3. Halo bias and stochasticity
The conditional mass function we examined in Sec-tion 3.2 is an indicator of how halo formation is biasedtoward overdense regions. The halo bias can be viewedalso in the the k -space. This is a crucial parameter incosmology when trying to probe the fluctuation of thematter density from surveys of galaxies through, for ex-ample, the power spectrum analysis.The halo bias parameter in k -space is defined as b ( k ) = (cid:18) P h ( k ) P m ( k ) (cid:19) / , (19)where P h ( k ) is the power spectrum of halo over-abundance δ n ( M, x ) = (cid:0) dndM (cid:1) ∆ ( M, x ) − (cid:0) dndM (cid:1) g ( M ) (cid:0) dndM (cid:1) g ( M ) . (20) δ n ( M, x ) can then be Fourier-transformed in order tocalculate P h ( k ). We expect b ( k ) to be larger than thatpredicted by TH, because the “contrast” in ( dn/dM ) ∆ between overdense and underdense regions has increasedfrom that by TH (Section 3.2; Figure 5). This is clearlyseen in Figure 7, where we show b ( k )’s computed in thiswork and by TH. At z = 19 and for M = 10 M ⊙ , wefind that b ( k ) is about [1 . −
2] times as large as the oneby TH. This is again caused by the mode-mode coupling.As was pointed out by TH, b ( k ) oscillates in k dueto the baryonic acoustic oscillations (BAO) , whichis tied to the modulation of the streaming velocity,the density fluctuation, and the baryon fraction in theexistence of the compensated mode (Barkana & Loeb2011). This makes it difficult to deduce P m ( k ) from P h ( k ), not to mention from the galaxy surveys wheregalaxy formation mechanism, strongly influenced bybaryonic physics, is another nuisance parameter (butsee Slepian & Eisenstein 2015 for how to separate outthe streaming velocity effect). For the halo mass rangetreated in this paper, b ( k ) should modulate the dis-tribution of the first stars which grow predominantlyinside minihalos. As was also noted by TH, the dif-ference in b ( k ) should also influence the formation ofmuch larger-mass halos, which are used for galaxy sur-veys. It is also possible that the nonlinear effect de-scribed by O’Leary & McQuinn (2012), or the heatingof the intergalactic medium (IGM) due to the velocitydifference between CDM and baryons, is modulated inspace depending on the overdensity. Then the powerspectrum in the 21-cm background, which may be dom-inated by the velocity fields if the heating is efficient(McQuinn & O’Leary 2012), is likely to be boosted. Be-cause such a power spectrum shows a very clear BAOfeature and the 21-cm observation usually suffers fromthe low sensitivity, the signal boosted even more fromthe prediction by McQuinn & O’Leary (2012) will be avery promising target for the high-redshift 21-cm cos-mology.Let us briefly discuss the halo stochasticity. The halostochasticity is defined as Our box size is barely larger than the BAO scale, and is thustoo small to accurately estimate b ( k ) at low k ’s. We will increasethe size of the box in future work for a better estimation. M (M ⊙ )0 . r a t i o − − − l og (cid:0) d n d M (cid:1) ∆ ( M ⊙− M p c − ) z=44 M (M ⊙ )0 . . r a t i o − − l og (cid:0) d n d M (cid:1) ∆ ( M ⊙− M p c − ) z=19 M (M ⊙ )0 . . r a t i o − − l og (cid:0) d n d M (cid:1) ∆ ( M ⊙− M p c − ) z=10 M (M ⊙ )0 . . . r a t i o − l og (cid:0) d n d M (cid:1) ∆ ( M ⊙− M p c − ) z=6 Figure 5 . Conditional halo mass functions, at varying redshifts (in title), for patches with different initial overdensities: ∆ c ( a i )= -0.005 (dotted), 0 (solid) and 0.005 (dashed). Thick (red) curves represent our predictions, and thin curves (blue) representpredictions by TH. There is no distinction between the predictions when ∆ c ( a i ) = 0, because then V bc is the only variablefor the ensemble in both cases. In each bottom sub-panel, we show the ratio of our prediction to the prediction by TH, or( dn/dM ) ∆ / ( dn/dM ) ∆ , TH . χ = P hm ( k ) P h ( k ) P m ( k ) , (21)where P hm ( k ) is the cross power spectrum between thehalo density and the matter density. Both in TH andin this work, stochasticity is caused by the fluctuationin ( dn/dM ) ∆ because patches with the same ∆ canhave different V bc ’s which affect ( dn/dM ) ∆ . One shouldnote that the fluctuation should be caused also by thesampling variance. The conditional mass functions usu-ally show super-Poissonian distributions in ( dn/dM ) ∆ even in the standard picture with V bc = 0 (e.g.Saslaw & Hamilton 1984; Sheth 1995; Neyrinck et al.2014; Ahn et al. 2015), and obviously this should cause the stochasticity in addition to that by the varying V bc .Because we use the “mean” conditional mass functionjust as TH did, in both works χ does not reflect thesampling variance. Including this effect requires the cal-culation of the sub-cell correlation function (Ahn et al.2015), which we delay to future work. DISCUSSIONWe investigated the impact of the relative velocity(streaming velocity) between CDM and baryons on thesmall-scale structure formation. TH first studied thiseffect by adopting a trivial solution to the large-scalevelocity and density fields. Because velocity fields arecorrelated with density fields, however, such a trivial so-lution cannot accurately describe the physics in regions3 M (M ⊙ )0 . . r a t i o this workTH − − l og (cid:0) d n d M (cid:1) ∆ ( M ⊙− M p c − ) this workTH V bc = 0 z=19, ∆ c ( a i ) = − . M (M ⊙ )1236 r a t i o this workTH − − l og (cid:0) d n d M (cid:1) ∆ ( M ⊙− M p c − ) this workTH V bc = 0 z=19, ∆ c ( a i ) = 0 . M (M ⊙ )12510 r a t i o this workTH − − l og (cid:0) d n d M (cid:1) g ( M ⊙− M p c − ) this workTH V bc = 0 z=19, global Figure 6 . Conditional mass functions at z = 19 with ∆ c ( a i ) = -0.005 (left) and 0.005 (middle). Comparison is made amongthis work (dashed; red), TH (dotted; blue) and the standard picture with V bc = 0 (solid; black). Global mass functions areplotted on the right panel. In each bottom sub-panel, we show the ratio of mass functions to that of the standard picture, or( dn/dM ) / ( dn/dM ) V bc =0 . .
03 0 . . k (Mpc − )20406080 d i ff e r e n ce ( % ) b ( k ) this workTH z=19 Figure 7 . Halo bias parameter for halos with M = 10 M ⊙ at z = 19. Comparison is made between this work(black, x) and TH (red, square) in the top sub-panel. Inthe bottom sub-panel, we show the fractional difference[ b ( k ) − b ( k ) TH ] /b ( k ) TH . with non-zero overdensity. We thus improved on thework by TH by implementing a non-trivial solution tothe large-scale velocity and density fields, and we findthat this causes a new type of coupling between large-scale and small-scale modes. This results in boostingthe small-scale structure formation in overdense regionsand suppressing that in underdense regions, aside from the suppression originating from the streaming velocity.The net effect on the structure formation is to boostthe overall fluctuation, in terms of P m ( k ), and thus the“negative” effect by TH is mitigated to some extent. De-pending on the wave mode ( k ) and the observing red-shift, P m ( k ) can even be larger than that in the standardpicture with V bc = 0. The conditional halo mass func-tion and the halo bias are also affected in similar ways.The results of this work show that the formation andevolution of small-scale structures depend strongly onnot only the streaming velocity but also the density en-vironment. The most important aspect of our work isthat in contrast to TH, who predict that regardless ofthe underlying density the local matter power spectrumof small-scale structures will be identical as long as V bc is the same, the underlying large-scale ( ∼ a few Mpc)overdensity is another key parameter in addition to V bc .This then requires re-examining previous work based onthe formalism by TH. We already showed that P loc , m ( k ), P m ( k ), ( dn/dM ) ∆ , ( dn/dM ) g , and b ( k ) are affected. Ifone were to simulate the nonlinear evolution of densityperturbations and the structure formation in ∼ z = 200, because at ∼ σ ∆ c the discrepancy be-tween our prediction and that by TH is already a fewpercent at that redshift.Both the previous semi-analytical work(Tseliakhovich et al. 2011; Fialkov et al. 2012;McQuinn & O’Leary 2012; Bovy & Dvorkin 2013;Naoz & Narayan 2014; Asaba et al. 2016) and the semi-numerical work (e.g. Fialkov et al. 2013; Visbal et al.42012) should be re-examined. Attempts to numeri-cally simulate the nonlinear evolution of small-scalestructures have been mostly limited to the physicsinside mean-density regions (McQuinn & O’Leary 2012;Maio et al. 2011; Stacy et al. 2011; Greif et al. 2011)or special, isolated regions (Tanaka & Li 2014). Thesenumerical simulations thus need to be extended toincorporate ∆ which varies in space. In doing so, areasonable method would be to use adaptive meshrefinement (AMR) codes with nested grids, so thatone or a few interesting regions ( ∼ ± σ ∆ , for example) are treated with fine meshesand other regions with coarse meshes for computationalefficiency.We also showed that cosmology through galaxy sur-veys should carefully consider the impact of the mode-mode couping, because the halo bias (and galaxy biasas well) would be boosted from not only the stan-dard prediction with V bc = 0 but also the predic-tion by TH. Cosmology with the intensity mappingmay also be affected. The post-reionization intensitymapping targets the large-angle, diffuse 21-cm back-ground from neutral hydrogen atoms inside galaxies(Chang et al. 2008; Abdalla et al. 2010; Bandura et al.2014; Xu et al. 2015). Because any galaxies, smallor large, contribute to this cumulative 21-cm back-ground, such observations will be affected by the stream-ing velocity through b ( k ). The pre-reionization in-tensity mapping targets the large-angle, diffuse 21-cm background from the intergalactic neutral hy-drogen atoms (Scott & Rees 1990; Bharadwaj & Ali2004; Loeb & Zaldarriaga 2004; Barkana & Loeb 2005;McQuinn et al. 2006; ; Mao et al. 2012; Shapiro et al.2013). Even though this is free from the galaxy bias,the streaming velocity may act as a heating mecha-nism and boost the power spectrum of the velocity field(McQuinn & O’Leary 2012), and therefore the new find-ings of our work should be incorporated.Application to the study of the cosmic reion-ization process is of a prime interest in termsof the high-redshift astrophysics. The complexnature of the process usually requires numericalsimulations, and they are performed through ei-ther efficient semi-numerical methods (Furlanetto et al.2004; McQuinn et al. 2007; Mesinger et al. 2011;Alvarez & Abel 2012) or fully numerical methods(e.g. Gnedin & Abel 2001; Razoumov et al. 2002;Maselli et al. 2003; Mellema et al. 2006; Baek et al. 2009; Wise & Abel 2011). The early phase of cosmicreionization must have been driven by the first stars,possibly forming first in minihalos, as these are the firstluminous objects in the universe. A very importantfactor that modulates the formation of the first starsinside minihalos is the Lyman-Werner intensity, whichhave been properly treated in simulations in a box largeenough for statistical reliability but implementing sub-grid physics for the first star formation inside minihalos(Ahn et al. 2012; Fialkov et al. 2013). To predict reion-ization scenarios to our best knowledge, especially on itsearly phase, this work should be properly incorporatedbecause the star formation inside minihalos is stronglymodulated by the streaming velocity as well.The results of this paper have rooms for further im-provements. This paper is based on the presumptionthat the fluctuations at a few Mpc scale remain lin-ear even at later epochs. However, high density re-gions will reach the nonlinear regime earlier than therest, and then our formalism will break down in suchregions. The halo bias is more strongly pronounced inthe nonlinear patches (e.g. Ahn et al. 2015) than in thelinear theory, and thus one should use the actual val-ues of the overdensity in such circumstances. One couldachieve this goal by adopting the quasi nonlinear calcu-lation (e.g. 2LPT by Crocce et al. 2006), adopting thetop-hat collapse model as in Mo & White (1996) andAhn et al. (2015), or for the best accuracy running N-body+hydro simulations which resolve the density fluc-tuation at ∼ Mpc scale. Then, in each patch of a fewcomoving Mpc, Equation (11) can be integrated with thenewly computed values of ∆’s. Wave modes in the range k ≃ [1 − / Mpc are not accurately treated, because webased our formalism on the separability of the large-scale modes ( k . / Mpc) and the small-scale modes( k & / Mpc). The code we used will be released forthe public use, but it requires technical improvementssuch as allowing parallel computation and porting tomore generic computation languages. We will main-tain and control its development through the website .We thank P. R. Shapiro, F. Schmidt and R. Barkanafor helpful discussions. We also thank the anonymousreferee for the clear report which led to a significant,quantitative improvement of the paper. This work wassupported by a research grant from Chosun University(2016) and by NRF-2014R1A1A2059811.APPENDIX A. NORMAL MODES FOR THE LARGE-SCALE FLUCTUATIONSWhen the fluctuation of radiation components are neglected, the growth of large-scale density and velocity fluctua-tions are well approximated by Equations (5) and (6). Their evolution can then be described by the 4 normal modes5described in Section 2.1. The growing and decaying modes are the two solutions to the second-order equation d ∆ + dt + 2 H d ∆ + dt − H Ω m ∆ + = 0 , (A1)where ∆ + = f c ∆ c + f b ∆ b , and can be written also as d ∆ + da + (cid:18) a + d ln Hda (cid:19) d ∆ + da − a Ω m ∆ + = 0 . (A2)Similarly, the compensated and streaming modes are the two solutions to d ∆ − dt + 2 H d ∆ − dt = 0 , (A3)where ∆ − = ∆ c − ∆ b , and can be written also as d ∆ − da + (cid:18) a + d ln Hda (cid:19) d ∆ − da = 0 . (A4)We now take a convention of writing each mode as the product of its initial value at z = 1000 and its growthfactor: ∆ g+ ( a ) = ∆ gro D g ( a ), ∆ d+ ( a ) = ∆ dec D d ( a ), ∆ c − ( a ) = ∆ com = constant, and ∆ s − ( a ) = ∆ str D s ( a ), denotingthe growing, decaying, compensated, and streaming modes, respectively. These modes comprise ∆ + and ∆ − , as∆ + ( a ) = ∆ gro D g ( a ) + ∆ dec D d ( a ) and ∆ − ( a ) = ∆ com + ∆ str D s ( a ), with the normalization D g = D d = D s = 1 at z = 1000.The first task in finding these modes is to calculate the growth factors. Because the Hubble constant H ( a ) (= H p Ω r, a − + Ω m, a − + Ω Λ , ) and Ω m ( a ) have non-negligible radiation components during the period of interest,1000 & z &
50, growth factors should be calculated numerically. We factor out deviations from the analyticalform valid during the matter-dominated (Ω m = 1) era, as D g ( a ) = ( a/a i ) F g ( a ), D d ( a ) = ( a/a i ) − / F d ( a ), and D s ( a ) = ( a/a i ) − / F s ( a ), and then solve for the order-of-unity values of F g , F d , and F s . They are determined by d F g da + (cid:18) a + d ln Hda (cid:19) dF g da − a (cid:26) a (cid:18) Ω m − (cid:19) − d ln Hda (cid:27) F g = 0 , (A5) d F d da + d ln Hda dF d da − a (cid:26) a (cid:18) Ω m (cid:19) + 32 d ln Hda (cid:27) F d = 0 , (A6)and d F s da + (cid:18) a + d ln Hda (cid:19) dF s da − a (cid:26) a + 12 d ln Hda (cid:27) F s = 0 . (A7)In practice, the numerical integration of Equations (A5)-(A7) is started at z = 10 (backward for z ≥
10 and forward for z ≤ dF g , d , s /da = 0 at z = 10 because it is a matter-dominated epoch, and F g , d , s = 1 at z = 1000 because of our normalization convention. One can instead start the integration from the radiation dominatedepoch (just numerically assuming that Equations (A5)-(A7) are all valid at z ≃ dF g , d , s /da and F g , d , s at that time), but we find that Equations (A6) and (A7) become stiff if integrated forward inincreasing a at high z . We show the growth factors found this way in Figure A1.Now we can find the initial values of these modes by using the growth factors found above on the transfer functionoutputs from CAMB. We algebraically relate two redshift outputs, at z ≡ z i = 1000 and z ≡
800 in practice, to findthese modes: ∆ gro = ∆ + ( a ) − D d ( a )∆ + ( a ) D g ( a ) − D d ( a ) , ∆ dec = ∆ + ( a ) − ∆ gro , ∆ str = ∆ − ( a ) − ∆ − ( a ) D s ( a ) − , ∆ com = ∆ − ( a ) − ∆ str , (A8)where ∆ ± ( a , ) are those from CAMB. The modes found this way are shown in Figure A2.A few things are notable. An unperturbed Hubble flow requires ∆ gro / ∆ dec = 3 /
2, while we find that at z = 1000∆ gro / | ∆ dec | ≃
27 and is not constant over k . Even though one can choose different redshifts to extract these modesand they should not change in principle, the resulting modes at z i vary depeding on the choice of the redshifts. We6 − .
01 0 . a . . . F F g F d F s − .
01 0 . a − . D D g D d D s Figure A1 . Evolution of growth factors, for the growing (black, solid), decaying (blue, dotted), and streaming (red, dashed)modes. When F (left panel) is multiplied to each analytic power-law evolution valid for the matter-dominated era, it results inthe actual growth factor D (right panel). .
01 0 . k (Mpc − )10 − . ∆ ( M p c / ) ∆ gro ∆ dec ∆ com ∆ str Figure A2 . Growing, decaying, compensated, and streaming modes at z = 1000, plotted in varying line widths. The negativevalues are plotted by a dashed line for the decaying mode and a dotted line for the streaming mode, after flipping the sign. believe that this is partly due to our neglect of the fluctuations of radiation components, because they can affect theevolution of density fluctuations. We find that our current choice of z and z is optimal for k & .
01 Mpc: using thesemodes and evolving them with the growth factors, we find a good match between { ∆ + , ∆ − } “constructed” by usingthese modes and those calculated by CAMB, at any redshifts with at most a several percent error. We show theircomparison in Figure A3.Schmidt (2016) follows a similar approach for the mode extraction but uses CAMB transfer function outputs at z ≃
0. Their focus is on the low-redshift galaxy surveys, and thus the accuracy is required mostly at and near thepresent. In our case, accuracy in ∆ j and Θ j is required mostly at a high redshift range, 1000 . z .
50, becausethe small-scale modes are continuously affected by large-scale modes from the epoch of recombination and the linearperturbation analysis on small-scales modes breaks down later when they become nonlinear. B. FITTING FORMULA FOR THE LARGE-SCALE TEMPERATURE FLUCTUATIONFor the evolution of ∆ T , we integrate Equation (8) on each spatial patch. Early on, its value is strongly affectedby the initial fluctuation of the CMB temperature. Later, it decouples from the CMB fluctuation and is determinedmostly by ∆ b . Therefore, one expects a strong correlation between ∆ T and ∆ b long after recombination. For example,at a = 0 .
01 and 0 .
1, they follow the linear relation: ∆ T / ∆ b = 0 .
239 at a = 0 .
01, and ∆ T / ∆ b = 0 .
586 at a = 0 . T , therefore, one can find a fitting formula for ∆ T ( a ) after its decoupling from theCMB temperature.The fitting formula is given by Equation (9). We note that Equation (9) is valid only when the patch size is 4comoving Mpc. For patches in different size, we believe that a generic form of a two-parameter fit,∆ T ( a ) = sign(∆ T, A ) dex [ α (log a + β ) γ − Y ] ,α = log (∆ T, B / ∆ T, A )( β − γ − ( β − γ , .
01 0 . k (Mpc − )00 . d i ff e r e n ce ( % ) ∆ + ∆ − − . . ∆ ( M p c / ) ∆ con+ ∆ CAMB+ ∆ con − ∆ CAMB − z=900 .
01 0 . k (Mpc − ) − − . d i ff e r e n ce ( % ) ∆ + ∆ − . . ∆ ( M p c / ) ∆ con+ ∆ CAMB+ ∆ con − ∆ CAMB − z=500 .
01 0 . k (Mpc − ) − − d i ff e r e n ce ( % ) ∆ + ∆ − . . ∆ ( M p c / ) ∆ con+ ∆ CAMB+ ∆ con − ∆ CAMB − z=200 Figure A3 . Comparison of the density fluctuations constructed by the normal modes (with superscript “con”) and thosecalculated by CAMB (with superscript “CAMB”), at different redshifts. It is difficult to see the difference in the upper panelsbecause the match is good, but depending on the value of k , some error is inherent as seen in the lower panels. The differenceis with respect to the values calculated by CAMB. We allow this degree of error in this paper, but for higher accuracy one isadvised to use the CAMB outputs for the large-scale fluctuations. − .
01 0 . a − − − . . ∆ T ∆ c ( a i ) = 4 . × − ∆ c ( a i ) = 1 . × − ∆ c ( a i ) = 4 . × − Figure B4 . Actual evolution of ∆ T of three arbitrarily chosen spatial patches (4 comoving Mpc) with varying ∆ c ( a i ) (thick;solid, dashed, dot-dashed). Overlaid are the corresponding fitting functions given by Equation 9, in thin dotted lines. Y = ( β − γ log | ∆ T, B | − ( β − γ log | ∆ T, A | ( β − γ − ( β − γ , (B9)will serve as a good fit regardless of the patch size. The linear relation between ∆ T and ∆ b at a = 0 .
01 and 0 . β = 2 . γ = 0 . T when ∆ T & − , with the linear relations ∆ T / ∆ b = 0 .
279 at a = 0 . T / ∆ b = 0 .
599 at a = 0 . T ’s are showndepending on the initial ∆ c , together with the corresponding fits.REFERENCES Abdalla, F. B., Blake, C., & Rawlings, S. 2010, MNRAS, 401, 743Ahn, K., Iliev, I. T., Shapiro, P. R., et al. 2012, ApJL, 756, L16Ahn, K., Iliev, I. T., Shapiro, P. R., & Srisawat, C. 2015,MNRAS, 450, 1486Alvarez, M. A., & Abel, T. 2012, ApJ, 747, 126Asaba, S., Ichiki, K., & Tashiro, H. 2016, PhRvD, 93, 023518Baek, S., Di Matteo, P., Semelin, B., Combes, F., & Revaz, Y.2009, A&A, 495, 389 Bandura, K., Addison, G. E., Amiri, M., et al. 2014, inProc. SPIE, Vol. 9145, Ground-based and Airborne TelescopesV, 914522Barkana, R., & Loeb, A. 2004, ApJ, 609, 474—. 2005, ApJL, 624, L65—. 2011, MNRAS, 415, 3113Bernardeau, F., Colombi, S., Gazta˜naga, E., & Scoccimarro, R.2002, PhR, 367, 1Bharadwaj, S., & Ali, S. S. 2004, MNRAS, 352, 142 Blazek, J., McEwen, J. E., & Hirata, C. M. 2015, ArXiv e-prints,arXiv:1510.03554Bond, J. R., Cole, S., Efstathiou, G., & Kaiser, N. 1991, ApJ,379, 440Bovy, J., & Dvorkin, C. 2013, ApJ, 768, 70Chang, T.-C., Pen, U.-L., Peterson, J. B., & McDonald, P. 2008,Physical Review Letters, 100, 091303Crocce, M., Pueblas, S., & Scoccimarro, R. 2006, MNRAS, 373,369Dalal, N., Pen, U.-L., & Seljak, U. 2010, JCAP, 11, 007Fialkov, A., Barkana, R., Tseliakhovich, D., & Hirata, C. M.2012, MNRAS, 424, 1335Fialkov, A., Barkana, R., Visbal, E., Tseliakhovich, D., & Hirata,C. M. 2013, MNRAS, 432, 2909Furlanetto, S. R., Zaldarriaga, M., & Hernquist, L. 2004, ApJ,613, 1Gnedin, N. Y., & Abel, T. 2001, New Astronomy, 6, 437Greif, T. H., White, S. D. M., Klessen, R. S., & Springel, V.2011, ApJ, 736, 147Komatsu, E., Smith, K. M., Dunkley, J., et al. 2011, ApJS, 192,18Lewandowski, M., Perko, A., & Senatore, L. 2015, JCAP, 5, 019Lewis, A., Challinor, A., & Lasenby, A. 2000, Astrophys. J., 538,473Loeb, A., & Zaldarriaga, M. 2004, Physical Review Letters, 92,211301Maio, U., Koopmans, L. V. E., & Ciardi, B. 2011, MNRAS, 412,L40Mao, Y., Shapiro, P. R., Mellema, G., et al. 2012, MNRAS, 422,926Maselli, A., Ferrara, A., & Ciardi, B. 2003, MNRAS, 345, 379McQuinn, M., Lidz, A., Zahn, O., et al. 2007, MNRAS, 377, 1043McQuinn, M., & O’Leary, R. M. 2012, ApJ, 760, 3McQuinn, M., Zahn, O., Zaldarriaga, M., Hernquist, L., &Furlanetto, S. R. 2006, ApJ, 653, 815Mellema, G., Iliev, I. T., Alvarez, M. A., & Shapiro, P. R. 2006,NewA, 11, 374Mesinger, A., Furlanetto, S., & Cen, R. 2011, MNRAS, 411, 955 Mo, H. J., & White, S. D. M. 1996, MNRAS, 282, 347Naoz, S., & Barkana, R. 2005, MNRAS, 362, 1047Naoz, S., & Narayan, R. 2013, Physical Review Letters, 111,051303—. 2014, ApJL, 791, L8Naoz, S., Yoshida, N., & Gnedin, N. Y. 2012, ApJ, 747, 128Neyrinck, M. C., Arag´on-Calvo, M. A., Jeong, D., & Wang, X.2014, MNRAS, 441, 646O’Leary, R. M., & McQuinn, M. 2012, ApJ, 760, 4Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2015,ArXiv e-prints, arXiv:1502.01589Popa, C., Naoz, S., Marinacci, F., & Vogelsberger, M. 2015,ArXiv e-prints, arXiv:1512.06862Razoumov, A. O., Norman, M. L., Abel, T., & Scott, D. 2002,ApJ, 572, 695Reichardt, C. L., Shaw, L., Zahn, O., et al. 2012, ApJ, 755, 70Richardson, M. L. A., Scannapieco, E., & Thacker, R. J. 2013,ApJ, 771, 81Saslaw, W. C., & Hamilton, A. J. S. 1984, ApJ, 276, 13Schmidt, F. 2016, ArXiv e-prints, arXiv:1602.09059Scott, D., & Rees, M. J. 1990, MNRAS, 247, 510Shapiro, P. R., Mao, Y., Iliev, I. T., et al. 2013, Phys. Rev. Lett.,110, 151301Sheth, R. K. 1995, MNRAS, 274, 213Shoji, M., & Komatsu, E. 2009, ApJ, 700, 705Slepian, Z., & Eisenstein, D. J. 2015, MNRAS, 448, 9Somogyi, G., & Smith, R. E. 2010, PhRvD, 81, 023524Stacy, A., Bromm, V., & Loeb, A. 2011, ApJL, 730, L1+Tanaka, T. L., & Li, M. 2014, MNRAS, 439, 1092Tanaka, T. L., Li, M., & Haiman, Z. 2013, MNRAS, 435, 3559Tseliakhovich, D., Barkana, R., & Hirata, C. M. 2011, MNRAS,418, 906Tseliakhovich, D., & Hirata, C. 2010, PhRvD, 82, 083520Visbal, E., Barkana, R., Fialkov, A., Tseliakhovich, D., & Hirata,C. M. 2012, Nature, 487, 70Wise, J. H., & Abel, T. 2011, MNRAS, 414, 3458Xu, Y., Wang, X., & Chen, X. 2015, ApJ, 798, 40Yoo, J., Dalal, N., & Seljak, U. 2011, JCAP, 7, 018 his workTHhis workTH