# How Density Environment Changes the Influence of the Dark Matter-Baryon Streaming Velocity on the Cosmological Structure Formation

###### Abstract

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.

## 1 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. 2012 Planck 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 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 recombination 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 large-scale streaming velocity fields which are coherent over a few comoving Mpc scale. This has its own spatial fluctuation with km/s standard deviation at the epoch of recombination, and then decays in proportion to the inverse of the scale factor . Its impact on the small-scale 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 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; McQuinn & O’Leary 2012; 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; McQuinn & O’Leary 2012) 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; O’Leary & McQuinn 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 ) 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; McQuinn & O’Leary 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 (McQuinn & O’Leary 2012), 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 (“-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 (), TH used a trivial solution for the evolution of the streaming velocity () and the density ( and being overdensities of the CDM and baryons in large scale, respectively) environment: , . Treating this as a new 0th-order solution to the perturbation equations, they then examined how the perturbation in small scales grows. In doing so, they treated as a spatial quantity and perturbation variables in small scales as -space quantities. This is basically a high-order perturbation theory, coupling large-scale mode () and the small-scale modes ( and , small scale CDM and baryonic overdensities, respectively). However, is tightly linked to the fluctuating and 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 to and , and Schmidt (2016) finds that the divergence of is indeed an important term that one should not ignore. Blazek et al. (2015) also works on the generic basis of non-zero and .

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 is treated as the spatial quantity and small-scale overdensities and are treated as the -space quantity in the perturbation analysis. Most importantly, we explicitly include generically “non-zero” and 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.

## 2 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 modes of interest are in sub-horizon scale such that the Newtonian perturbation theory holds. Let us define the overdensity , where the subscript ={c, b} denotes either the CDM (c) or the baryonic (b) component, and where is the scale factor and is the proper peculiar velocity of component . 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)

(1) |

where is the scale factor, is the gradient in the comoving frame, , , , , is the present-day mean density of component in the unit of the critical density, is the sound speed, and is the Hubble constant at a given redshift. Even though a usual approximation for is a spatially uniform one given by

(2) |

which assumes a mean-density environment undergoing Hubble expansion
with the mean baryon temperature , for high modes a
more accurate treatment is required (Naoz & Barkana 2005). This requires
replacing the pressure term^{1}^{1}1We keep the cross term
in Equation (3)
in order to find the correct coupling of high- and low-
modes in Equations (10) and (11).,

(3) |

in Equation (1) and considering another rate equation

(4) |

where and are temperature fluctuations of the baryon and the photon, respectively, is the mean photon temperature, is the global electron fraction at time and . For , we use the fitting formula by TH.

Equation (1) allows a trivial solution: , and , where 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- and large- 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 is coherent over the length scale of a few comoving Mpc (contributed by modes with wave numbers in the range ), and averaged the “local” power spectra of density fluctuations over many such patches with varying .

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):

(5) |

where is the matter content with respect to the critical density at , and we ignored second-order terms and also the pressure term . This is indeed a valid approximation in the wave number range () relevant to the coherent (TH), where the second-order terms remain much smaller than the first-order terms and the baryonic sound speed keeps decreasing from km/s after recombination to make the pressure term negligible. Then, Equation (5) can be rewritten as

(6) |

where , , , and . In the matter-dominated () flat universe, allows both the growing mode (, , ) and the decaying mode (, , ). allows a slowly decaying (“streaming”) mode (, , ) and a compensated mode (=constant, , ). During , the non-negligible amount of the radiation component (CMB and neutrinos) makes , and most of the simple analytical forms above become no longer intact except for {, } of the streaming mode and {, , } of the compensated mode.

Using this mode decomposition, the large-scale perturbations evolve in the following form: {widetext}

(7) |

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. , , , and are the initial (at ) values of the growing, decaying, compensated, and streaming modes, respectively. , , and are growth factors of the growing, decaying, and streaming modes, respectively, and they are all normalized as at . 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 and , indicating that baryons tend to move together with CDMs in time. More interestingly, decays as throughout the evolution at any overdensity environment even though and grow roughly as individually (except in regions with where and decay as ). This fact may seem to make the analysis by TH valid in generic overdensity environments to some extent: TH relied on the trivial solution , in which all velocity components decay in time such that , , and most importantly . Because the suppression of the matter-density fluctuations () depends not on individual velocity components but on 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- and low- modes in general. We also require the evolution of , which is given by Equation (4):

(8) |

which is evolved in conjunction with Equation (1). Here we neglect term due to its smallness (see also the following discussion), even though we do not neglect when initializing (Section 2.2). We have found a useful fitting formula for for patches with volume :

(9) | |||||

which provides a good fit to at , and for higher redshift range we simply ignore altogether because of smallness of in general. Here , and , and an almost complete coupling of to at yields this simple, empirical relation to . We describe this fitting formula in more details in Appendix B. The decoupling of from the CMB is much earlier than the mean value experiences, which occur at , because the larger the is, the earlier the decoupling occurs (see e.g. Figure 1 of Naoz & Barkana 2005). Including explicitly may delay this decoupling to some extent, but we leave such an accurate calculation to future work. Using the evolution equation for 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 ), the net gravitational potential and the net baryon temperature as , , , and , respectively. Here and 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 (e.g. Bernardeau et al. 2002). Similarly, . We then have

(10) |

where we marked the terms that can be further ignored in square brackets. First, we can safely ignore any terms containing and , because the former is negligible at compared to (see Figure 1) and the latter is just too small ( at and decaying in time) to produce any appreciable impact on the baryon temperature of high modes. Secondly, the justification for ignoring is easily seen in viewpoint of the -space. With the Fourier expansion of a quantity in an actual, real space (“ space” henceforth), the background velocities have but with the condition . In contrast, the small-scale modes fluctuating against the background patches have intrinsically larger wave number , or . Then, at each , . Similarly, . 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. top-hat 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 (as in O’Leary & McQuinn 2012; TH chose the baryon rest frame), Equation (10), in the -space, finally becomes

(11) | |||||

where , and now denote fluctuations in the -space while , and are fluctuations of a given patch at in the space, given by Equation (7).

### 2.2 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 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 -space quantities whose distributions are all Gaussian, one needs to sample these values in the -space accordingly. For the small-scale modes, one just needs to track the evolution of the average value in the -space.

Let us first describe the statistics of background patches that we
expect. TH calculated the evolution of small-scale ()
fluctuations under different background patches but only of ,
and defined “local power spectrum”
averaged out over all possible opening angles between
and . In our case, there are extra dimensions to consider
which are , , ,
and , 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 completely compose the ensemble at any time through
Equation (7). At the minimal level^{2}^{2}2For a more accurate treatment or for a specific patch of interest,
one should also consider variations in other variables, such that
.
As the correlation between () and
() becomes tighter in time, the
initial variation gets gradually diluted, which roughly justifies
our restricting the parameter space only to and ., it would
suffice to just consider variation of in addition
to such that the local power spectrum is an explicit
function of the two background quantities, or .

For the initial condition for background patches, we generate 3D maps of , , , , , , and at on uniform grid cells inside a cubical volume of . We generate fluctuations of discrete modes that are randomized as

(12) |

for given , where and are random numbers drawn from mutually independent Gaussian distributions with mean 0 and standard deviation 1, stands for any kind of -space fluctuations, and TF is the transfer function of , whose sign should be multiplied because some ’s oscillate around zero in . The configuration is roughly equivalent to applying a smoothing filter of length Mpc. In practice, we use CAMB (Lewis et al. 2000) for , and use the continuity equations () for , with the help of two CAMB transfer-function outputs at mutually nearby redshifts for time differentiation. is obtained from the relation . is fixed by following the scheme by Naoz & Barkana (2005): we require at the initial redshift in Equation (4), which results in

(13) |

where all quantities are evaluated at , especially with the help of Equation (7) for and two adjacent CAMB transfer-function outputs for . Finally, all these -space fluctuations are Fourier-transformed to obtain -space fluctuations.

3D maps and 2D histograms of several initial quantities are presented in Figure 2. Fields of and on a part of a slice of the box at 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 dominates over , and thus the map of looks very similar to Figure 2(a). This occurs because baryons lag behind CDMs due to their coupling to CMB. (Figure 2b) is coupled to (Figure 2c) more strongly than to . and (similarly and ) are almost perfectly correlated (Figure 2d). and 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. and are not correlated (Figure 2f). Because of this fact, the probability distribution function (PDF) is simply a multiplication of PDFs and , at the minimal level. Due to Gaussianity at , we have

(14) |

where and are the standard deviations of and projected onto one Cartesian-coordinate axis, respectively. With our setup, we find that and km/s at (or the root-mean-square of is km/s). At the minimal level of only allowing the variance in and , the average power spectrum will then be given by the ensemble average

(15) | |||||

where the PDFs and integral arguments are the ones at . Of course, a more accurate and straightforward way is to just ensemble-average over the patches from a large-box realization, because for example and are too poorly correlated at .

We numerically integrate Equation (11) to examine the evolution of small-scale (high-) 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®(2015b, The MathWorks, Inc., Natick, Massachusetts, United States) and of GNU Octave, which use the 4th-order Runge-Kutta method, with the relative tolerance and the absolute tolerance . During the evolution, the number of integration steps is the highest for , because its amplitude changes from the initial, very small values around to final, much larger values close to . Therefore, taking sub-steps for while coarser steps for other ’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 and to obtain (Equation 15).

## 3 Result

### 3.1 Power spectrum of the matter density

We first examine how the evolution of depends on the density environment, and compare the result to the prediction by TH. Figure 3 shows the evolution of of three arbitrarily chosen wave numbers (={33, 150, 2000} /Mpc) when in different density environments (={-0.01, -0.005, 0, 0.005, 0.01}). Note again that at , and thus these samples correspond to and . First, as expected, the growth of small-scale fluctuations are biased when and anti-biased when , with respect to the mean-density case (prediction by TH). Secondly, when ’s are equal in amplitude but opposite in sign, the deviations of from reveal the same trend but only until when and when . 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 is, the earlier this unbalance starts. Thirdly, the fractional deviation from the mean-density case is almost universal regardless of the value of , and thus the timing of the unbalance is approximately a function only of . Finally, in some high- patches, our linear analysis based on Equation (11) without quadratic and higher-order terms in starts to break down at , because these modes enter the nonlinear regime (; 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, grows in time and . Then, in Equation (11), and work as sources terms in addition to to the growth of . This will boost the growth rate of of both overdense (, ) and underdense (, ) modes. In contrast, in underdense background patches, these terms suppress the growth rate of of both overdense and underdense modes. The distribution of is Gaussian, and therefore for each “bias” case with there exists an “anti-bias” case with . 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 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 boosts even in the range ( and ) where the power spectrum is almost unaffected by non-zero (Figure 4 shown in terms of the -space matter-density variance ). Note that the result for cannot be trusted, because our perturbation theory is based on the condition that large-scale modes () are well separated from small-scale modes in scale. Discrepancy of including non-zero ’s from the prediction by TH is negligible at , but later the discrepancy grows in time. Of course, individual patches may experience discrepancy in