{'fraction_non_alphanumeric': 0.04803562280805454, 'fraction_numerical': 0.019308822949915278, 'mean_word_length': 4.205107168495539, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 10, 'https://': 0, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 11, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
How Density Environment Changes the Influence of the Dark Matter-Baryon Streaming Velocity on the Cosmological Structure Formation

We study the dynamical effect of relative velocities between dark matter and baryonic fluids, which remained supersonic after the epoch of recombination. The impact of this supersonic motion on the formation of cosmological structures was first formulated by Tseliakhovich&Hirata (2010), in terms of 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 those of the global mean density. We improve on their formulation by allowing variation in the density environment 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: faster in the overdense environment and slower in the underdense environment. We also find that the net effect on the global power spectrum of the density fluctuation is to boost its overall amplitude from the prediction by Tseliakhovich&Hirata (2010). Correspondingly, the conditional mass function of cosmological halos and the halo bias parameter are both affected in a similar way. The discrepancy between our prediction and that by Tseliakhovich&Hirata (2010) is significant, and therefore the related cosmology and high-redshift astrophysics should be revisited. The mathematical formalism of this study can be used for generating cosmological initial conditions of small-scale perturbations in generic, overdense (underdense) background patches.


INTRODUCTION
The Λ-cold dark matter (ΛCDM) scenario, combined with the theory of cosmic inflation, is the successful, concurrent model describing the past and the present of our universe, consistent with a wide range of observations. In this scenario, cosmological structures grow out of an extremely uniform density field but with tiny fluctuations that are seeded by the cosmic inflation. The growth of the CDM density fluctuations and the growth of baryon density fluctuations are not in perfect synchronization, because baryons were tightly coupled to photons before the epoch of recombination and thus their motion was different from the motion of CDM which only reacts to gravity. Only after recombination baryons gradually decoupled from photons, and followed the motion of the CDM under gravity.
Cosmological observations have verified the ΛCDM scenario in scales large enough to make the baryonic physics almost irrelevant (e.g. Komatsu et al. 2011;Reichardt et al. 2012Planck Collaboration et al. 2015. However, once in the regime where the baryonic physics becomes important, the growth of baryon fluctuations is affected by hydrodynamics and the growth of the CDM is affected by the gravitational feedback from the baryon fluctuations. Some of the usual assumptions that are made for treating very large scales, therefore, should be taken carefully or modified when treating relatively small scales. For example, Naoz & Barkana (2005) improved on the previous estimation of the linear density power spectrum in small scales, by replacing the usual assumption made in cosmology that the sound speed of baryons is uniform in space with the fact that the sound speed fluctuates in space in small scales. They showed that more than ∼ 10 % change occurs in the baryon density power spectrum and even more change in the baryon temperature power spectrum. A sheer inclusion of the sub-dominant, yet non-negligible baryonic component in the analysis changes the prediction on the matter density power spectrum at a few percent level even in large scales, as was shown in the framework of the high-order perturbation theory (Shoji & Komatsu 2009;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 the first time, properly calculated the growth of small-scale density fluctuations under the influence of the streaming velocity. Small-scale fluctuations are coupled to largescale streaming velocity fields which are coherent over a few comoving Mpc scale. This has its own spatial fluctuation with ∼ 30 km/s standard deviation at the epoch of recombination, and then decays in proportion to the inverse of the scale factor a. Its impact on the smallscale structure formation is non-negligible because the streaming velocity remains supersonic (until the intergalactic medium is strongly heated). This then leads to the suppression of small-scale matter density fluctuations. The wave-modes that are the most strongly affected are around k ∼ 200/Mpc, and the impact is on the matter power spectrum, the conditional mass function and the halo bias parameter, to name a few (TH).
Subsequent studies have considered the impact of the streaming velocity in the perspective of both cosmology and astrophysics. The boost of amplitude and the shift of the peak of the baryonic acoustic oscillation (BAO) feature due to the streaming velocity, in the gas intensity mapping or the galaxy survey, were intensively investigated (Dalal et al. 2010;Yoo et al. 2011;Slepian & Eisenstein 2015;Lewandowski et al. 2015;Blazek et al. 2015;Schmidt 2016). They find that both high-redshift (e.g. Dalal et al. 2010; and low-redshift (e.g. Yoo et al. 2011) surveys will be affected. The impact of the streaming velocity on BAO may be separated out from the impact of the matter density itself (Slepian & Eisenstein 2015), which is important because otherwise it will become another nuisance parameter in cosmology. The astrophysical impact of the streaming velocity has been investigated with focus on the formation of the nonlinear structure such as cosmological halos and stellar 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;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 that the formation of minihalos (roughly in the mass range M = [10 4 − 10 8 ] M ⊙ ) and the formation of stellar objects in them are suppressed. It may induce baryon-dominated objects such as globular clusters (Naoz & Narayan 2014;Popa et al. 2015), but the actual star formation process leading to globular clusters has yet to be simulated. It may be responsible even for the generation of the primordial magnetic field (Naoz & Narayan 2013). Because minihalos are the most strongly affected among cosmological halos and they are responsible for the early phase of the cosmic reionization process, how they change the high-redshift 21-cm background is also of a prime interest (Visbal et al. 2012;Fialkov et al. 2013). It is noteworthy that some of these numerical simulation results (Maio et al. 2011;Stacy et al. 2011;Greif et al. 2011), which are based on the initial condition generated by the usual Boltzmann solver such as the Code for Anisotropies in the Microwave Background (CAMB: Lewis et al. 2000), need to be re-examined. This is because the impact of the streaming velocity is cumulative and inherent even at z ∼ 200 , at which or later these simulations start with a streaming velocity implemented by hand.
The original formalism by TH has been re-investigated in terms of the high-order perturbation theory in the wave-number space ("k-space" henceforth), and some "missing terms" previously neglected were found important (Blazek et al. 2015;Schmidt 2016). Basically, for the large-scale modes responsible for the streaming velocity (k ∼ [0.01 − 1]/Mpc), TH used a trivial solution for the evolution of the streaming velocity (V bc ) and the density (∆ c and ∆ b being overdensities of the CDM and baryons in large scale, respectively) environment: V bc ∝ a −1 , ∆ c = ∆ b = 0. Treating this as a new 0thorder solution to the perturbation equations, they then examined how the perturbation in small scales grows. In doing so, they treated V bc as a spatial quantity and perturbation variables in small scales as k-space quantities. This is basically a high-order perturbation theory, coupling large-scale mode (V bc ) and the small-scale modes (δ c and δ b , small scale CDM and baryonic overdensities, respectively). However, V bc is tightly linked to the fluctuating ∆ c and ∆ b through the density continuity equation, and thus the trivial solution adopted by TH cannot be used for generically overdense and underdense regions. The continuity equation connects the divergence of V bc to ∆ c and ∆ b , and Schmidt (2016) finds that the divergence of V bc is indeed an important term 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 considering the non-zero overdensities. Toward this end, different from Blazek et al. (2015) and Schmidt (2016), we inherit the original method by TH and focus on the impact on small-scale modes: large-scale V bc is treated as the spatial quantity and small-scale overdensities δ c and δ b are treated as the k-space quantity in the perturbation analysis. Most importantly, we explicitly include generically "non-zero" ∆ c and ∆ b as a new set of spatial quantities, and consequently the divergence of CDM and baryon velocities as well. We find that this leads to a set of mode-mode coupling terms, including the velocity divergence-density coupling. We carefully include all the coupling terms to the leading order in our perturbation analysis. We also include the baryonic physics, namely fluctuations in the sound speed (Naoz & Barkana 2005), the gas temperature and the photon temperature. Our formalism is also suitable for generating initial conditions for N-body+hydro numerical simulations. Because the new set of mode-mode couplings is imprinted in the initial condition, the initial condition generator considering the streaming-velocity effect by O' Leary & McQuinn (2012), CICsASS, should also be improved on if one were to numerically simulate the structure formation inside overdense or underdense regions.
This paper is organized as follows. After the introduction, we lay out the basic formalism and describe the statistics of large-scale fluctuations in Section 2. In Section 3, we show results on the matter density power spectrum, the conditional halo abundance and the halo bias as applications of the formalism. We conclude this work in Section 4 with a summary, discussion and future prospects. Some details left out in the main body are described in Appendices.

FORMALISM AND NUMERICAL METHOD
2.1. Fluctuation under non-zero overdensity and relative velocity: perturbation formalism We start from a set of equations for perturbations of relevant physical variables. All the k modes of interest are in sub-horizon scale such that the Newtonian perturbation theory holds. Let us define the overdensity δ j ≡ (ρ j −ρ j ) /ρ j , where the subscript j={c, b} denotes either the CDM (c) or the baryonic (b) component, and θ j ≡ (1/a) ∇ · v j where a is the scale factor and v j is the proper peculiar velocity of component j. When the universe is in the regime where we can ignore fluctuations of photons and neutrinos due to their rapid diffusion after recombination, we have (e.g. Bernardeau et al. 2002;TH) where a is the scale factor, ∇ is the gradient in the comoving frame, , Ω j,0 ≡ρ j,0 /ρ crit,0 is the present-day mean density of component j in the unit of the critical density, c s is the sound speed, and H is the Hubble constant at a given redshift. Even though a usual approximation for c s is a spatially uniform one given by which assumes a mean-density environment undergoing Hubble expansion with the mean baryon temperaturē T , for high k modes a more accurate treatment is required (Naoz & Barkana 2005). This requires replacing the pressure term 1 , in Equation (1) and considering another rate equation (4) where δ T and δ Tγ are temperature fluctuations of the baryon and the photon, respectively,T γ ≡ 2.725 K (1 + z) is the mean photon temperature, x e (t) is the global electron fraction at time t and t γ ≡ 1.17 × 10 12 yr. For T , we use the fitting formula by TH.
Equation (1) allows a trivial solution: where a i is the initial scale factor. TH took this as the zeroth-order solution of a spatial patch and developed a linear perturbation theory of small scale modes inside the patch. The physical process is in principle a coupling of small-k and large-k modes (mode-mode coupling), which is beyond the linear theory where all modes are assumed to be mutually independent. TH used the fact that the relative velocity V bc ≡ v b − v c is coherent over the length scale of a few comoving Mpc (contributed by modes with wave numbers in the range 0.01 (k/Mpc −1 ) 1), and averaged the "local" power spectra of density fluctuations over many such patches with varying V bc .
Even though a trivial solution exists, there also exists a nontrivial solution to Equation (1), exact to the first order. This nontrivial solution is suited to describe the physics inside patches with non-zero overdensity (see also Blazek et al. 2015). In order to obtain the nontrivial solution, we first linearize Equation (1): where Ω m ≡ρ m (a)/ρ crit (a) is the matter content with respect to the critical density ρ crit at a, and we ignored second-order terms and also the pressure term a −1 c 2 s ∇δ b . This is indeed a valid approximation in the wave number range (0.01 (k/Mpc −1 ) 1) relevant to the coherent V bc (TH), where the second-order terms remain much smaller than the first-order terms and the baryonic sound speed keeps decreasing from c s ∼ 6 km/s after recombination to make the pressure term negligible. Then, Equation (5) can be rewritten as where In the matter-dominated (Ω m = 1) flat universe, δ + allows both the growing mode (δ + ∝ a, θ . During 1000 z 50, the non-negligible amount of the radiation component (CMB and neutrinos) makes Ω m = 1, and most of the simple analytical forms above become no longer intact except for {θ − , v − } of the streaming mode and {δ − , θ − , v − } of the compensated mode. Using this mode decomposition, the large-scale perturbations evolve in the following form: where we used upper-case letters to denote the "background" fluctuations for each patch of a few comoving Mpc, over which we will develop the small-scale perturbation. ∆ 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, and streaming modes, respectively, and they are all normalized as D g = D d = D s = 1 at z = 1000. We describe the details of these modes in Appendix A. It is noteworthy, as is well known already, that both CDM and baryonic components tend to approach the same asymptotes ∆ = ∆ gro D g (a) and Θ = −aH∆ gro dD g (a)/da, indicating that baryons tend to move together with CDMs in time. More interestingly, V bc decays as a −1 throughout the evolution at any overdensity environment even though V c and V b grow roughly as a 1/2 individually (except in regions with ∆ c = ∆ b = 0 where V c and V b decay as a −1 ). This fact may seem to make the analysis by TH valid in generic overdensity environments to some extent: TH relied on the trivial solution ∆ c = ∆ b = 0, in which all velocity components decay in time such that Because the suppression of the matter-density fluctuations (f c δ c + f b δ b ) depends not on individual velocity components but on V bc only, different temporal behavior of individual velocity components among the trivial and generic solutions do not matter. Nevertheless, quantitative prediction by TH will be questioned in Section 3.1, because we use nontrivial solutions (Equation 7) which result in the new type of coupling of high-k and low-k modes in general. We also require the evolution of ∆ T , which is given by Equation (4): which is evolved in conjunction with Equation (1). Here we neglect ∆ Tγ term due to its smallness (see also the following discussion), even though we do not neglect ∆ Tγ (a i ) when initializing ∆ T (a i ) (Section 2.2). We have found a useful fitting formula for ∆ T (a) for patches with volume (4 Mpc) 3 : which provides a good fit to ∆ T at 300 z 5, and for higher redshift range we simply ignore ∆ T altogether because of smallness of ∆ T in general. Here , and an almost complete coupling of ∆ T (a) to ∆ b (a) at z 300 yields this simple, empirical relation to ∆ b (a). We describe this fitting formula in more details in Appendix B. The decoupling of ∆ T (a) from the CMB is much earlier than the mean value experiences, which occur at z ≃ 150, because the larger the k is, the earlier the decoupling occurs (see e.g. Figure 1 of Naoz & Barkana 2005). Including ∆ Tγ explicitly may delay this decoupling to some extent, but we 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 the 4 modes. Now we expand Equations (1), (3) and (4) to the linear order, taking Equation (7) as the zeroth-order solution. We first define the net density, the net velocity (of the fluid component i), the net gravitational potential and the net baryon temperature as Here X and x denote the comoving-coordinate position of the center of a background patch and that of a small-scale fluid component, respectively, and we use lower-case letters for small-scale fluctuations. Φ is the gravitational potential sourced only by the background fluctuations such that Bernardeau et al. 2002). Similarly, where we marked the terms that can be further ignored in square brackets. First, we can safely ignore any terms containing δ Tγ and ∆ Tγ , because the former is negligible at z 1000 compared to δ T (see Figure 1) and the latter is just too small (∆ Tγ 10 −5 at z = 1000 and decaying in time) to produce any appreciable impact on the baryon temperature of high k modes. Secondly, the justification for ignoring a −1 (v i · ∇) V i is easily seen in viewpoint of the k-space. With the Fourier expansion A(x) = k A(k) exp(ik · x) of a quantity A(x) in an actual, real space ("r space" henceforth), the background velocities have V j (x) = K V j (K) exp(iK · x) but with the condition K 1 Mpc −1 . In contrast, the small-scale modes fluctuating against the background patches have intrinsically larger wave number k, or We finally note that we do not include the coupling terms between similar wavenumbers in Equation (10), which will involve quadratic and higher-order polynomicals of ∆ and δ. Therefore, the validity of Equation (10) will break down in the non-linear regime. Nevertheless, our approach is more suitable for a crude estimattion of the conditional halo mass function (Section 3.2) in terms of the extended Press-Schechter formalism, which is based on the mapping of the linear density growth to the nonlinear growth of the e.g. tophat density perturbation.
The evolution of small-scale perturbations, therefore, is coupled to large-scale perturbations on which they are sitting. Ignoring the terms in square brackets, and shifting the viewpoint to the CDM rest frame in which where δ j , θ j and δ T now denote fluctuations in the kspace while ∆ j , Θ j and V j are fluctuations of a given patch at (a, X) in the r space, given by Equation (7).

Evolution of perturbation inside patches: Numerical Scheme
Evolution of small-scale perturbations can be calculated by integrating the rate equation (Equation 11) from some initial redshift, preferentially not too long after the recombination epoch when the relative motion has not yet influenced the evolution. We take z i ≡ 1000 as the initial redshift. The initial condition should be generated for both the background quantities and the small-scale modes. For the background, as perturbations in Equation (11) are r-space quantities whose distributions are all Gaussian, one needs to sample these values in the r-space accordingly. For the small-scale modes, one just needs to track the evolution of the average value in the k-space.
Let us first describe the statistics of background patches that we expect. TH calculated the evolution of small-scale (k 10) fluctuations under different background patches but only of ∆ c = ∆ b = 0, and defined "local power spectrum" P loc,m (k; V bc ) averaged out over all possible opening angles between V bc and k. In our case, there are extra dimensions to consider which are ∆ c , ∆ b , Θ c , and Θ b , resulting in a much higher computational demand. Fortunately, some of these quantities are in perfect correlation with one another with linear proportionality. In addition, their initial values at a = a i completely compose the ensemble at any time through Equation (7). At the minimal level 2 , it would suffice to just consider variation of ∆ c in addition to V bc such that the local power spectrum is an explicit 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, we generate 3D maps of ∆ c , ∆ b , Θ c , Θ b , V c , V b , and ∆ T at z = 1000 on 151 3 uniform grid cells inside a cubical volume of V box = (604 Mpc) 3 . We generate fluctuations of discrete modes that are randomized as for given k, where G 1 and G 2 are random numbers drawn from mutually independent Gaussian distributions with mean 0 and standard deviation 1, ∆ k stands for any kind of k-space fluctuations, and TF is the transfer function of ∆ k , whose sign should be multiplied because some ∆ k 's oscillate around zero in k.
The configuration is roughly equivalent to applying a smoothing filter of length (604/151) = 4 Mpc. In practice, we use CAMB (Lewis et al. 2000) for ∆ j (k, a i ), and use the continuity equations (∂∆ j /∂t = −Θ j ) for Θ j (k, a i ), with the help of two CAMB transferfunction outputs at mutually nearby redshifts for time is fixed by following the scheme by Naoz & Barkana (2005): we require ∂∆ T /∂t = ∂∆ Tγ /∂t at the initial redshift in Equation (4), which results in where all quantities are evaluated at a i , especially with 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/4)∂∆ γ (k, a i )/∂t. Finally, all these k-space fluctuations are Fourier-transformed to obtain r-space fluctuations.
3D maps and 2D histograms of several initial quantities are presented in Figure 2. Fields of ∆ c and V bc on a part of a slice of the box at z = z i = 1000 are shown in Figure 2(a). As expected, the velocity field converges on overdense regions and diverges on underdense regions. We find that in most patches V c dominates over V b , and thus the map of V c looks very similar to Figure 2(a). This occurs because baryons lag behind CDMs due to their 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 due to the tight coupling of baryons to photons at the redshift (Figure 2e), but the correlation becomes tighter in time. ∆ c and V bc are not correlated (Figure 2f). Because 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 Gaussianity at z i , we have 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.0042 and σ V bc = 17.8 km/s at z i (or the root-mean-square of V bc is √ 3σ V bc = 30.9 km/s). At the minimal level of only allowing the variance in ∆ c and V bc , the average power spectrum P m (k) will then be given by the ensemble average where the PDFs and integral arguments are the ones at z i . Of course, a more accurate and straightforward way is to just ensemble-average P loc,m = P loc,m k; V bc ; ∆ c ; ∆ b ; Θ c ; Θ b ; ∆ T ; ∆ Tγ over the 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 examine the evolution of small-scale (high-k) fluctuations at any overdense (underdense) patch to the linear order, with the help of Equations (7) and (9) for the evolution of background quantities. In practice, we used the ODE45 modules of MATLAB R (2015b, The MathWorks, Inc., Natick, Massachusetts, United States) and of GNU Octave, which use the 4th-order Runge-Kutta method, with the relative tolerance 10 −4 and the absolute tolerance 10 −2 δ k (a i ). During the evolution, the number of integration steps is the highest for δ T , because its amplitude 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 's is expected to boost the computational efficiency, even though we did not yet implement the method in our computation. The end result is then ensemble-averaged over varying ∆ c and V bc to obtain P m (k) (Equation 15).

Power spectrum of the matter density
We first examine how the evolution of P loc,m (k; V bc ; ∆ c ) depends on the density environment, and compare the result to the prediction by TH. Figure 3 shows the evolution of ∆ 2 loc,m ≡ k 3 P loc,m /2π 3 of three arbitrarily chosen wave numbers (k={33, 150, 2000} /Mpc) when V bc = 22 km/s (a/a i ) −1 in different density environments (∆ c (a i )={-0.01, -0.005, 0, 0.005, 0.01}). Note again that σ ∆c = 0.0042 at z = 1000, and thus these samples correspond to ±2.4σ ∆c and ±1.2σ ∆c . First, as expected, the growth of small-scale fluctuations are biased when ∆ c > 0 and anti-biased when ∆ c < 0, with respect to the mean-density case (prediction by TH). Secondly, when ∆ c (a i )'s are equal in amplitude but opposite 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 but only until z ≃ 850 when |∆ c (a i )| = 0.01 and z ≃ 730 when |∆ c (a i )| = 0.005. Afterwards, the bias and the anti-bias are not balanced by more than 1% and such off-balance keeps growing in time. The higher the |∆ c ( a i )| is, the earlier this unbalance starts. Thirdly, the fractional deviation from the mean-density case is almost universal regardless of the value of k, and thus the timing of the unbalance is approximately a 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 δ starts to break down at z ∼ 20, because these modes enter the nonlinear regime (∆ 2 loc,m (k) 0.1; see Figure 3) at this epoch. A similar breakdown of the formalism will occur for perturbations inside very low-∆ patches as well. Therefore, a higher-order scheme than our work is required when one is to predict the low-redshift evolution of δ's. On the other hand, if one were to generate initial conditions for numerical simulations in the linear regime, our formalism would provide the sufficient accuracy. How large-scale overdensity impacts the evolution of small-scale inhomogeneities is reflected in the density continuity 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 rate of δ j of both overdense (δ j > 0, θ j < 0) and underdense (δ j < 0, θ j > 0) modes. In contrast, in underdense background patches, these terms suppress the growth rate of δ j of both overdense and underdense modes. The distribution of ∆ j (a i ) is Gaussian, and therefore for each "bias" case with ∆ j > 0 there exists an "anti-bias" case with ∆ j < 0. Nevertheless, the unbalance described above is expected to boost the average power spectrum from that by TH.
The overall effect of including the Gaussian distribution of ∆ c is thus to mitigate the negative impact by the relative 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 100/Mpc and k 1000/Mpc) where the power spectrum is almost unaffected by non-zero V bc (Figure 4 shown in terms of the k-space matter-density variance ∆ 2 m ≡ k 3 P m (k)/2π 2 ). Note that the result for k 10/Mpc cannot be trusted, because our perturbation theory is based on the condition that large-scale modes (0.01 K Mpc 1) are well separated from small-scale modes in scale. Discrepancy of P m (k) including non-zero ∆'s from the prediction by TH is negligible at z 45, but later the discrepancy grows in time. Of course, individual patches may experience discrepancy in P loc,m (k; V bc ; ∆ c ) much earlier than this epoch (Figure 3).

Halo abundance
Understanding the abundance and the spatial distribution of cosmological halos is crucial in modern astrophysics and cosmology. In this section, we examine the halo 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. They adopted the extended Press-Schechter formalism and calculated the local halo abundance, in terms of the con-  ditional mass function, using the peak-background split scheme. A patch with ∆ (let us use this notation for the matter overdensity, to avoid confusion with the k-space matter-density variance ∆ 2 m ) and V bc will have the number of halos per unit Eulerian comoving volume per M given by where δ crit is the critical overdensity of spherical collapse, and σ 2 is the variance of density field smoothed with the window function W M corresponding to mass M , One should note that σ 2 (M, V bc ) should be that of high k modes only, or more accurately a reduced value σ 2 (M, V bc ) − σ 2 patch where σ 2 patch is the variance of density field smoothed with the window function corresponding to the mass of the patch 3 (e.g. Bond et al. 1991;Mo & White 1996;Ahn et al. 2015). One can instead put a lower bound k min (= 6π 2ρ m /M patch 1/3 ) in the integral of Equation (17), which would be identical to the reduced variance if a sharp k-space window function is used. The global mass function (dn/dM ) g is simply an average of the local mass function, dn(M |∆, V bc )/dM , over the ensemble of patches. When dn(M |∆, V bc )/dM is averaged only over V bc for a given ∆, which is equivalent to visiting only those patches with the same ∆ and taking the average, it leads to the conditional mass function (dn/dM ) ∆ ≡ dn(M |∆)/dM . This calculation should be modified, because σ 2 depends also on ∆ through the dependence of P loc,m on ∆: which will enter Equation (16). Here we used the sharp k-space filter, and thus k max = 6π 2ρ m /M 1/3 . An overdense patch will then have a boost in dn(M |∆, V bc )/dM from the value by TH because P m (k, V bc , ∆ > 0) > P m (M, V bc , ∆ = 0) and thus σ 2 (M, V bc , ∆ > 0) > σ 2 (M, V bc , ∆ = 0), and vice versa (a decrease from the 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 usual peak-background split scheme and the one by TH. In the "standard" scheme, if the density field is purely Gaussian, all the wave modes are assumed mutually independent in the linear regime. Therefore, the local, highk modes have a universal 4 variance σ 2 k whether or not they are placed inside a patch with non-zero ∆. The way how the halo formation is biased in an overdense region is simply through the shift in the density (+∆). TH then realized the fact that the variance is not universal but should depend on V bc . Because non-zero V bc tends to suppress σ 2 k , the odds to cross δ crit decrease relative to the standard picture. We find that there is another dependency of the variance, which is ∆. In other words, we find that there are two biasing effects in an overdense region compared to a mean-density region: getting closer to δ crit because of the shift in the density (+∆, also in the standard scheme), and having a larger degree of fluctuation in δ k due to the mode-mode coupling (e.g. source terms −∆ c θ c and −Θ c δ c in ∂δ c /∂t in Equation 11, which is a new finding). Therefore, by not fully implementing the effect of non-zero ∆, TH in effect underestimates and overestimates the halo mass functions in overdense and underdense regions, respectively.
We note that this additional bias effect should be present even in the standard picture with V bc = 0, because this is due to the natural coupling between the large-scale and small-scale density perturbations. In this case, however, we are not sure about the quantitative validity of the extended Press-Schechter formalism on the conditional mass function (Equation 16), which is based on the linear theory guaranteeing Gaussianity at any filtering scales without the mode-mode couplings. Qualitatively, we believe that the boost of local δ k and σ 2 k under ∆ should boost the conditional mass function to the level estimated by Equations (16) and (18) as described above anyways. We defer a further investigation of this issue, which can be clarified with numerical simulations of the halo formation under different ∆'s.
Discrepancy between the conditional mass functions by this work and by TH is significant if we focus on individual patches. Figure 5 illustrates how our prediction differs from that by TH. For example, at z ∼ 44 − 19, under ∆ c (a i ) = 0.005 we predict [100 -2000] % boost in (dn/dM ) ∆ compared to the values by TH (let us denote them by (dn/dM ) ∆, TH ). For ∆ c (a i ) = −0.005, we predict 90 % or more decrease in (dn/dM ) ∆ compared to the values by TH. This is the obvious result of the mode-mode coupling of ∆ and δ described above. It is also noteworthy that the discrepancy is the largest for the rarest halos: first, at any redshift, the discrepancy increases as the halo mass increases and secondly, for any given halo mass, the discrepancy decreases in time.
The discrepancy among the conditional mass functions by this work, by TH and by the standard picture also influences (dn/dM ) g . Let us just take the example of z = 19 ( Figure 6). Compared to the standard prediction (V bc = 0), conditional mass functions by TH stay lower regardless of ∆. This is due to the suppression of structure formation by the relative velocity. In contrast, (dn/dM ) ∆ in this work is either higher or lower than the standard prediction (≡ (dn/dM ) V bc =0 ) depending on the halo mass and ∆. At z = 19, for ∆ c (a i ) = −0.005 (dn/dM ) ∆ < (dn/dM ) V bc =0 and for ∆ c (a i ) = 0.005 (dn/dM ) ∆ > (dn/dM ) V bc =0 for any halo mass. The tendency for (dn/dM ) ∆ to overshoot (dn/dM ) V bc =0 when ∆ > 0 is not generic, because the influence of the positive ∆'s on the structure formation appears only late in its evolution (see Figs. 3 and 4). At any rate, the contrast in (dn/dM ) ∆ among overdense and underdense regions is increased compared to TH 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 indistinguishable from that by TH, which undershoots (dn/dM ) g, V bc =0 .
We note that (dn/dM ) ∆ and (dn/dM ) g have the usual problem of not correctly predicting the actual mass function, if one sticks to the original extended Press-Schechter formalism. Minihalos at high redshifts are usually underestimated by the extended Press-Schechter formalism. The usual peak-background split method suffers from large discrepancies between its prediction and the N-body simulation results for rare halos in general. In this case, a hybrid method to connect the peakbackground-split halo bias parameter to the betterfitting mean mass function types (e.g. Barkana & Loeb 2004;Ahn et al. 2015) is much more appropriate. We will apply this method in the future for a better estimation of the conditional mass function and the k-space halo bias parameter (Section 3.3).

Halo bias and stochasticity
The conditional mass function we examined in Section 3.2 is an indicator of how halo formation is biased toward overdense regions. The halo bias can be viewed also in the the k-space. This is a crucial parameter in cosmology when trying to probe the fluctuation of the matter density from surveys of galaxies through, for example, the power spectrum analysis.
The halo bias parameter in k-space is defined as where P h (k) is the power spectrum of halo overabundance .
δ n (M, x) can then be Fourier-transformed in order to calculate P h (k). We expect b(k) to be larger than that predicted by TH, because the "contrast" in (dn/dM ) ∆ between overdense and underdense regions has increased from that by TH (Section 3.2; Figure 5). This is clearly seen in Figure 7, where we show b(k)'s computed in this work and by TH. At z = 19 and for M = 10 6 M ⊙ , we find that b(k) is about [1.5 − 2] times as large as the one by TH. This is again caused by the mode-mode coupling. As was pointed out by TH, b(k) oscillates in k due to the baryonic acoustic oscillations (BAO) 5 , which is tied to the modulation of the streaming velocity, the density fluctuation, and the baryon fraction in the existence of the compensated mode (Barkana & Loeb 2011). This makes it difficult to deduce P m (k) from P h (k), not to mention from the galaxy surveys where galaxy formation mechanism, strongly influenced by baryonic physics, is another nuisance parameter (but see Slepian & Eisenstein 2015 for how to separate out the streaming velocity effect). For the halo mass range treated in this paper, b(k) should modulate the distribution of the first stars which grow predominantly inside minihalos. As was also noted by TH, the difference in b(k) should also influence the formation of much larger-mass halos, which are used for galaxy surveys. It is also possible that the nonlinear effect described by O' Leary & McQuinn (2012), or the heating of the intergalactic medium (IGM) due to the velocity difference between CDM and baryons, is modulated in space depending on the overdensity. Then the power spectrum in the 21-cm background, which may be dominated by the velocity fields if the heating is efficient , is likely to be boosted. Because such a power spectrum shows a very clear BAO feature and the 21-cm observation usually suffers from the low sensitivity, the signal boosted even more from the prediction by McQuinn & O'Leary (2012) will be a very promising target for the high-redshift 21-cm cosmology.
Let us briefly discuss the halo stochasticity. The halo stochasticity is defined as  Figure 5. Conditional halo mass functions, at varying redshifts (in title), for patches with different initial overdensities: ∆c(ai) = -0.005 (dotted), 0 (solid) and 0.005 (dashed). Thick (red) curves represent our predictions, and thin curves (blue) represent predictions by TH. There is no distinction between the predictions when ∆c(ai) = 0, because then V bc is the only variable for 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.
where P hm (k) is the cross power spectrum between the halo density and the matter density. Both in TH and in this work, stochasticity is caused by the fluctuation in (dn/dM ) ∆ because patches with the same ∆ can have different V bc 's which affect (dn/dM ) ∆ . One should note that the fluctuation should be caused also by the sampling variance. The conditional mass functions usually 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 function just as TH did, in both works χ does not reflect the sampling variance. Including this effect requires the calculation of the sub-cell correlation function (Ahn et al. 2015), which we delay to future work.

DISCUSSION
We investigated the impact of the relative velocity (streaming velocity) between CDM and baryons on the small-scale structure formation. TH first studied this effect by adopting a trivial solution to the large-scale velocity and density fields. Because velocity fields are correlated with density fields, however, such a trivial solution cannot accurately describe the physics in regions with non-zero overdensity. We thus improved on the work by TH by implementing a non-trivial solution to the large-scale velocity and density fields, and we find that this causes a new type of coupling between largescale and small-scale modes. This results in boosting the small-scale structure formation in overdense regions and suppressing that in underdense regions, aside from the suppression originating from the streaming velocity. The net effect on the structure formation is to boost the overall fluctuation, in terms of P m (k), and thus the "negative" effect by TH is mitigated to some extent. Depending on the wave mode (k) and the observing redshift, P m (k) can even be larger than that in the standard picture with V bc = 0. The conditional halo mass function and the halo bias are also affected in similar ways.
The results of this work show that the formation and evolution of small-scale structures depend strongly on not only the streaming velocity but also the density environment. The most important aspect of our work is that in contrast to TH, who predict that regardless of the underlying density the local matter power spectrum of 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 on the formalism by TH. We already showed that P loc,m (k), P m (k), (dn/dM ) ∆ , (dn/dM ) g , and b(k) are affected. If one were to simulate the nonlinear evolution of density perturbations and the structure formation in ∼ 4 Mpc patches, he should generate initial conditions based on this work. As is seen in Figure 3, we cannot neglect the impact of overdensity even when the simulation starts at e.g. z = 200, because at ∼ 1σ ∆c the discrepancy between our prediction and that by TH is already a few percent at that redshift.
Both the previous semi-analytical work (Tseliakhovich et al. 2011;Fialkov et al. 2012;Bovy & Dvorkin 2013;Naoz & Narayan 2014;Asaba et al. 2016) and the seminumerical work (e.g. Fialkov et al. 2013;Visbal et al. 2012) should be re-examined. Attempts to numerically simulate the nonlinear evolution of small-scale structures have been mostly limited to the physics inside mean-density regions Maio et al. 2011;Stacy et al. 2011;Greif et al. 2011) or special, isolated regions (Tanaka & Li 2014). These numerical simulations thus need to be extended to incorporate ∆ which varies in space. In doing so, a reasonable method would be to use adaptive mesh refinement (AMR) codes with nested grids, so that one or a few interesting regions (∼ 4 Mpc patches with ∆ = ±1σ ∆ , for example) are treated with fine meshes and other regions with coarse meshes for computational efficiency.
We also showed that cosmology through galaxy surveys should carefully consider the impact of the modemode couping, because the halo bias (and galaxy bias as well) would be boosted from not only the standard prediction with V bc = 0 but also the prediction by TH. Cosmology with the intensity mapping may also be affected. The post-reionization intensity mapping targets the large-angle, diffuse 21-cm background 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, small or large, contribute to this cumulative 21-cm background, such observations will be affected by the streaming velocity through b(k). The pre-reionization intensity mapping targets the large-angle, diffuse 21cm background from the intergalactic neutral hydrogen atoms (Scott & Rees 1990;Bharadwaj & Ali 2004;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 mechanism and boost the power spectrum of the velocity field , and therefore the new findings of our work should be incorporated.
Application to the study of the cosmic reionization process is of a prime interest in terms of the high-redshift astrophysics.
The complex nature of the process usually requires numerical simulations, and they are performed through either 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 cosmic reionization must have been driven by the first stars, possibly forming first in minihalos, as these are the first luminous objects in the universe. A very important factor that modulates the formation of the first stars inside minihalos is the Lyman-Werner intensity, which have been properly treated in simulations in a box large enough for statistical reliability but implementing subgrid physics for the first star formation inside minihalos (Ahn et al. 2012;Fialkov et al. 2013). To predict reionization scenarios to our best knowledge, especially on its early phase, this work should be properly incorporated because the star formation inside minihalos is strongly modulated by the streaming velocity as well.
The results of this paper have rooms for further improvements. This paper is based on the presumption that the fluctuations at a few Mpc scale remain linear even at later epochs. However, high density regions will reach the nonlinear regime earlier than the rest, and then our formalism will break down in such regions. The halo bias is more strongly pronounced in the nonlinear patches (e.g. Ahn et al. 2015) than in the linear theory, and thus one should use the actual values of the overdensity in such circumstances. One could achieve this goal by adopting the quasi nonlinear calculation (e.g. 2LPT by Crocce et al. 2006), adopting the top-hat collapse model as in Mo & White (1996) and Ahn et al. (2015), or for the best accuracy running N-body+hydro simulations which resolve the density fluctuation at ∼ Mpc scale. Then, in each patch of a few comoving Mpc, Equation (11) can be integrated with the newly computed values of ∆'s. Wave modes in the range k ≃ [1 − 10]/Mpc are not accurately treated, because we based our formalism on the separability of the largescale modes (k 1/Mpc) and the small-scale modes (k 10/Mpc). The code we used will be released for the public use, but it requires technical improvements such as allowing parallel computation and porting to more generic computation languages. We will maintain and control its development through the website http://www.chosun.ac.kr/kjahn. We thank P. R. Shapiro, F. Schmidt and R. Barkana for helpful discussions. We also thank the anonymous referee for the clear report which led to a significant, quantitative improvement of the paper. This work was supported by a research grant from Chosun University  (5) and (6). Their evolution can then be described by the 4 normal modes believe that this is partly due to our neglect of the fluctuations of radiation components, because they can affect the evolution of density fluctuations. We find that our current choice of z 1 and z 2 is optimal for k 0.01 Mpc: using these modes and evolving them with the growth factors, we find a good match between {∆ + , ∆ − } "constructed" by using these modes and those calculated by CAMB, at any redshifts with at most a several percent error. We show their comparison 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 the present. In our case, accuracy in ∆ j and Θ j is required mostly at a high redshift range, 1000 z 50, because the small-scale modes are continuously affected by large-scale modes from the epoch of recombination and the linear perturbation analysis on small-scales modes breaks down later when they become nonlinear.

B. FITTING FORMULA FOR THE LARGE-SCALE TEMPERATURE FLUCTUATION
For the evolution of ∆ T , we integrate Equation (8) on each spatial patch. Early on, its value is strongly affected by the initial fluctuation of the CMB temperature. Later, it decouples from the CMB fluctuation and is determined mostly 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.1. Regardless of the variance in ∆ T , therefore, one can find a fitting formula for ∆ T (a) after its decoupling from the CMB temperature.
The fitting formula is given by Equation (9). We note that Equation (9) is valid only when the patch size is 4 comoving Mpc. For patches in different size, we believe that a generic form of a two-parameter fit,  Figure A3. Comparison of the density fluctuations constructed by the normal modes (with superscript "con") and those calculated by CAMB (with superscript "CAMB"), at different redshifts. It is difficult to see the difference in the upper panels because the match is good, but depending on the value of k, some error is inherent as seen in the lower panels. The difference is with respect to the values calculated by CAMB. We allow this degree of error in this paper, but for higher accuracy one is advised to use the CAMB outputs for the large-scale fluctuations. Y = (β − 2) γ log 10 |∆ T, B | − (β − 1) γ log 10 |∆ T, A | (β − 1) γ − (β − 2) γ , will serve as a good fit regardless of the patch size. The linear relation between ∆ T and ∆ b at a = 0.01 and 0.1 can be found by comparing the two quantities. In case of the 4 comoving Mpc patch, we find that β = 2.8 and γ = 0.33 (Equation 9) provides an excellent fit to ∆ T when ∆ T 10 −4 , with the linear relations ∆ T /∆ b = 0.279 at a = 0.01 and ∆ T /∆ b = 0.599 at a = 0.1 This is demonstrated in Figure B4, where the evolution of different ∆ T 's are shown depending on the initial ∆ c , together with the corresponding fits.