The nonlinear evolution of baryonic overdensities in the early universe: Initial conditions of numerical simulations
We run very large cosmological -body hydrodynamical simulations in order to study statistically the baryon fractions in early dark matter halos. We critically examine how differences in the initial conditions affect the gas fraction in the redshift range . We test three different linear power spectra for the initial conditions: (1) A complete heating model, which is our fiducial model; this model follows the evolution of overdensities correctly, according to Naoz & Barkana (2005), in particular including the spatial variation of the speed of sound of the gas due to Compton heating from the CMB. (2) An equal- model, which assumes that the initial baryon fluctuations are equal to those of the dark matter, while conserving of the total matter. (3) A mean model, which assumes a uniform speed of sound of the gas. The latter two models are often used in the literature. We calculate the baryon fractions for a large sample of halos in our simulations. Our fiducial model implies that before reionization and significant stellar heating took place, the minimum mass needed for a minihalo to keep most of its baryons throughout its formation was M. However, the alternative models yield a wrong (higher by about ) minimum mass, since the system retains a memory of the initial conditions. We also demonstrate this using the ”filtering mass” from linear theory, which accurately describes the evolution of the baryon fraction throughout the simulated redshift range.
Recent measurements of anisotropies of the cosmic microwave background (CMB) radiation have revealed the detailed distribution of matter in the Universe a few hundred thousand years after the Big Bang (Spergel et al., 2003, 2007; Komatsu et al., 2009, 2010). Observations utilizing large ground-based telescopes and space telescopes have discovered galaxies and black holes that were in place when the age of the Universe was less than a billion years. Moreover, many galaxies have been found at (Bouwens et al., 2010; McLure et al., 2010) in the Hubble Ultra Deep Field, whereas already a few gamma-ray bursts at have been detected by the Swift satellite (Tanvir et al., 2009; Salvaterra et al., 2009; Lin, Liang, & Zhang, 2010). These first objects are probably the building blocks of the present day galaxies, thus, solving the puzzle behind their formation will have a profound implication on our understanding of the Universe (see for recent reviews Bromm et al., 2009; Yoshida, 2009, and references therein).
The formation of the first generation of galaxies in the Universe has been studied for many years. High resolution cosmological simulations can follow complex astrophysical processes, while analytical calculations can provide an over-all understanding, and can be used to decouple different physical effects seen in simulations. Analytic models are also useful for estimating the limitations of numerical simulations such as insufficient resolution and small boxsizes (Yoshida et al. 2003; Barkana & Loeb 2004; Naoz & Barkana 2005). Combining the two approaches may offer many of the advantages of both.
The initial conditions (hereafter ICs) in a cosmological simulation can have a large effect on the formation of the first galaxies in simulations, i.e., both on the formation time (or on the halo abundance at a given time) and the halo properties at formation time (such as the average gas fraction). Yoshida et al. (2003) studied high-redshift structure formation and reionization while testing two different models for power spectra as their ICs. They found that different models have a profound effect on the abundance of primordial star-forming gas clouds and thus on when the reionization was initiated and its progress. In the analytical point of view, Naoz, Noter & Barkana (2006) and Naoz & Barkana (2007) showed that the ICs at high redshift have a significant effect on the halo abundance and the gas fraction at virialization. While these effects are largest at the highest redshift, e.g., for the first star in the universe (Naoz, Noter & Barkana, 2006), they are still significant for halos forming at . The first gas rich halos at these redshifts are expected to host the first stars ( Naoz, Noter & Barkana, 2006; Yoshida, 2006; Gao et al., 2007; Trenti, Stiavelli, & Michael Shull, 2009a) and even the first gamma-ray bursts (e.g., Bromm & Loeb, 2006; Naoz & Bromberg, 2007). Thus, investigating the formation properties of these halos is of prime importance .
Gas rich halos in the early Universe may very well be a nurturing ground for dwarf galaxies, which at high redshift can form stars (e.g., Bromm, Coppi, & Larson, 2002, 1999; Abel, Bryan, & Norman, 2002; Yoshida et al., 2006; Yoshida, Omukai & Hernquist, 2008, and references therein) perhaps even at a high star formation rate (Ricotti, Gnedin, & Shull, 2002; Greif et al., 2010; Clark et al., 2010). Their properties are very important as they are responsible for metal pollution and the ionizing radiation at these early times (e.g., Shapiro, Iliev, & Raga, 2004; Ciardi et al., 2006; Gnedin, Kravtsov, & Chen, 2008; Trenti & Stiavelli, 2009b). Moreover, halos that are too small for efficient cooling via atomic hydrogen, i.e., minihalos, are most susceptible to the effect of initial conditions. While they may not normally host astrophysical sources, minihalos may produce a 21-cm signature (Kuhlen, Madau, & Montgomery (2006); Shapiro et al. (2006); Naoz & Barkana (2008) but see Furlanetto & Oh (2006)), and they can block ionizing radiation and produce an overall delay in the initial progress of reionization (e.g., Barkana & Loeb, 2002; Iliev et al., 2003, 2005; McQuinn et al., 2007). The evolution of the halo gas fraction at various epochs of the universe is of prime importance, particularly in the early universe.
In this paper, we examine the effect of using different initial conditions in simulations on the resulting minimum gas-rich halo mass in the redshift regime . We perform Gadget-2 (Springel, Yoshida, & White, 2001; Springel, 2005) simulations using a total of particles. We compare the initial conditions presented in Naoz & Barkana (2005), which describe the linear evolution of overdensities in a fully consistent way, to two other alternative ICs, often used in the literature. We also compare to the prediction of the gas-rich mass from linear theory. We describe our different initial conditions and simulations in sections 2 and 3, respectively. Our simulation results are presented in section 4 where we divide our discussion to the evolution of the non-linear power spectra (section 4.1) and to the minimum gas-rich halo mass resulting from either linear theory or from the simulations (section 4.2). Finally, we discuss our conclusions (section 5).
Throughout this paper, we adopt the following cosmological parameters: (, , , n, , )= (0.72, 0.28, 0.046, 1, 0.82, 70 km s Mpc) (Komatsu et al., 2009).
2 Different Initial Condition models - Basic Equations
2.1 The fiducial ICs - ”fid”
We follow Naoz & Barkana (2005), who studied the linear evolution of both dark matter and baryon overdensities. The fluctuations of the temperature of the baryons () cannot be described as a simple function of a spatially uniform baryonic sound speed , as was previously assumed (e.g., Ma & Bertschinger, 1995). Furthermore, at high redshifts, the baryon density fluctuations () are not equal to those of dark matter () (contrary to a common assumption in simulations; four redshift examples are shown in figure 1 of Naoz & Barkana, 2005). We label the power spectrum model as the ”fid” (fiducial) ICs since it follows the evolution of linear overdensities in a complete and consistent way.
Following Naoz & Barkana (2005) we write the basic equations that describe the evolution of the dark matter, baryon density and temperature fluctuations:
where and are the mean cosmic dark matter and baryonic fraction respectively. Here we follow the standard notations for cosmological parameters such as . The baryons are also subject to a pressure term:
where is the mean molecular weight, is the Boltzmann constant and is the wavenumber. Using the first law of thermodynamics, Naoz & Barkana (2005) derived the equations for the evolution of the baryon average temperature and temperature fluctuations:
where is the mean CMB temperature, and the first-order equation for the perturbation:
with the second term on the right-hand side accounting for the Compton scattering of the CMB photons on the residual electrons from recombination, where is the electron fraction out of the total number density of gas particles at time , and
where is the Thomson scattering cross section and is the photon energy density. The first term on the right-hand-side of each of these two equations (3) and (4) accounts for adiabatic expansion of the gas, and the remaining terms capture the effect of the thermal exchange with the CMB. Following Naoz & Barkana (2005) we have numerically calculated the evolution of the perturbations by modifying the CMBFAST code (Seljak & Zaldarriaga, 1996) according to these equations. Note that similar physics was also explored by Yamamoto et al. (1997, 1998).
We solve the complete set of equations to obtain the power spectrum at different redshifts which can be used as initial conditions for our simulations. Figure 1 shows the ratio between this initial condition to the two alternative models tested in this paper.
2.2 Alternative model I - equal - ”E-”
In many cosmological ICs for N-body simulations and semi-analytical calculations, the fluctuations of the baryons are assumed to be equal to the fluctuations of the dark matter. We construct a model that includes this incorrect assumption while maintaining the correct overall (i.e., conserving at , see appendix A for more details). Thus, in our “E-” model we calculate the correct as a combination of and from the fiducial calculation in section (2.1), but then take the baryon perturbation to be the same as for the dark matter, namely:
where () is the resulting linear over-density from the fiducial calculation (E- model) for the baryons and dark matter, respectively. We then compare the equal model to our fiducial calculation. Figure 1 shows the ratio between the fid ICs and the E- model for both the baryons and dark matter. We find that the E- model overestimates the baryon fluctuations by on large scales ( kpc) while the overestimate grows to a much larger factor on small scales.
Before recombination the baryons were tightly coupled to the radiation, resulting in suppression of the growth of their overdensity. However, the dark matter component, which is not affected by the photons, could basically grow once the fluctuation wavelength entered the Hubble horizon (in the linear regime, before equality, the dark mater fluctuations grew logarithmically with the scale factor, where after equality they grew linearly with the scale factor). Therefore, this resulted in a suppression of the baryonic overdensity by about three orders of magnitude compare to the dark matter at recombination (e.g., fig. 1 in Naoz & Barkana, 2005). While the baryons subsequently fall into the potential wells of the dark matter, it takes them some time to catch up, and the baryon fluctuations are still suppressed even at lower redshifts. This point is often overlooked in simulations and analytical calculations.
2.3 Alternative model II - the mean sound speed approximation - ”mean cs”
Naoz & Barkana (2005) showed that the presence of spatial fluctuations in the sound speed modifies the calculation of perturbation growth significantly. Nevertheless, for completeness and as a case of comparison with previous results, we compare the simulation results with the results obtained using this approximation. Thus, we proceed by presenting the basic equations of the growth of density fluctuations, in this approximation of a uniform sound speed (hereafter “mean ”). The evolution of the density fluctuations is described by a different set of coupled second order differential equations:
where is assumed to be spatially uniform (i.e., independent of ) and is thus calculated from the thermal evolution of a uniform gas undergoing Hubble expansion. With this assumption, the temperature fluctuations (as a function of ) are simply proportional at any given time to the gas density fluctuations:
Naoz & Barkana (2005) showed that this approximation leads to an underestimation of the baryon density fluctuations by up to 30% at and 10% at for large wavenumbers. Figure 1 shows the ratio between the mean initial conditions and the fiducial ones for both the baryons and dark matter. It agrees with our previous results, showing that the underestimate by the mean model is greatest at kpc. The non-linear evolution resulting from these initial conditions will result in shallower potential wells compared to the fiducial calculation,
Even though it is clear that the precise baryon temperature fluctuations at high redshift are very significant, still many simulations use initial conditions that assume a uniform speed of sound in the Universe. As shown below this leads to significantly different estimates for the gas content of the early halos.
3 The simulation
3.1 Basic parameters
We run a Gadget 2 simulation (Springel, Yoshida, & White, 2001; Springel, 2005) starting from redshift , for a total of particles ( particles each for the Dark Matter and baryon components) and our box size is: Mpc. We choose this box size so that a halo mass of M would have particles. This way according to Naoz, Barkana & Mesinger (2009) we are able to estimate the gas fraction in M halos correctly (see below for the halo definition). Our softening length is 0.2 comoving kpc.
For all runs, glass-like cosmological ICs were generated using the Zel’dovich approximation. The transfer functions were generated using the various models described above. We have used a glass file which was randomly displaced thus removing the coupling between nearby DM and gas particles. Using this randomization procedure we achieve essentially the same effect to that shown in Yoshida, Sugiyama, & Hernquist (2003). In generating the ICs, a convolution between the glass file and the transfer function from the different models was done, thus taking into account the different velocities of the DM and baryons (for the fiducial and mean cs models). We note that we have used the same phases for the DM and baryons, in all of the simulations.
We set the initial temperature to be K (as derived from linear theory), and thus Gadget assumes neutral and monoatomic gas, and converts to thermal energy (i.e., adiabatic initial conditions). Although this work emphasizes the need for a precise calculation of the baryon overdensities resulting from temperature fluctuations, we actually neglect the temperature fluctuations in the initial conditions. This may not be a bad approximation since the halos we study are already somewhat non-linear at our initial redshift, and the Compton heating is quite small compared to the adiabatic heating during non-linear gravitational collapse (see Appendix C). A more complete treatment would be to include in the simulation the precise temperature fluctuations, which we leave for future work. Nevertheless, even with the current treatment our results show consistency with linear theory.
3.2 Halo definition
We locate dark matter halos by running a FOF group-finder algorithm with a linking parameter of . We then find the center of mass of each halo and calculate the density profile of the dark matter and baryons, separately. In order to derive the density profile we assume a spherical halo, and divide it to 2000 shells. Combining these density profiles, we find the virial radius at which the overdensity is times the background density, and the gas fraction of each halo.
Recently, Trenti et al. (2010) performed a resolution analysis in order to study the mass definition of halos in simulations. Their conclusion (their figure 2) is that using the FOF algorithm and assuming about 500 particles per spherical halo introduces an error of in the mass definition. In our gas fraction analysis we have chosen only halos with a number of particles larger or equal to , i.e., we limit our errors in halo mass definition to below . Also, according to Naoz, Barkana & Mesinger (2009), this way we can estimate the gas fraction inside a halo accurately.
4 Results and comparison among the models
4.1 Non-linear power spectrum evolution
One way to probe cosmic structure particularly on small scales is through the non-linear power spectrum. We begin our simulation at with linear initial conditions111This is, of course, an approximation, ,since as shown in Naoz, Noter & Barkana (2006) at overdensities are already slightly non-linear. The effect of starting the simulation at high redshifts is studied elsewhere (Naoz et. al, in prep.). The main disagreement between the three models lies in the baryonic component (although the E- calculation also underestimates the dark matter overdensities by ). This input difference is then modified by the non-linear evolution.
Following Yoshida, Sugiyama, & Hernquist (2003) we compared the linear power spectrum for the fid model, as computed from Naoz & Barkana (2005), for the dark matter and baryon components, with the non-linear power spectra from the simulation (see figure 2). The two power spectra agree well as expected in the linear regime. We note that the other two models approach the fid model at low redshifts (see appendix A figure 5).
Fig 3 shows the differences among the fid, E- and mean ICs, in terms of the non-linear power spectra at the later redshifts at which halos were formed in our simulation. The mean model maintains over time roughly the same level of discrepancy with the fid model, while in the E- model both the baryonic and dark matter differences decline slightly slower than with the inverse scale factor. As clearly can be seen from figure 3, the non-linear evolution of halos is still strongly affected by the choice of initial conditions even at redshift . The fid ICs (Naoz & Barkana, 2005) describe the linear evolution consistently and thus represent the best available prescription for the initial conditions.
4.2 The minimum gas rich mass
Studying the galaxy evolution and reionization either by using simulations (both AMR and SPH) or by using analytical calculations relies on knowing the amount of gas within the dark matter halos. The simplest assumption, often used in the literature, is that a dark matter halo has the mean cosmic fraction. This can lead to incorrect results, especially when one tries to study star formation, galaxy mergers, and related phenomena.
Consider the various scales involved in the formation of non-linear objects containing DM and gas. On large scales (small wavenumbers) gravity dominates halo formation and gas pressure can be neglected. On small scales, on the other hand, the pressure dominates gravity and prevents baryon density fluctuations from growing together with the dark matter fluctuations. The relative force balance at a given time can be characterized by the Jeans (1928) scale, which is the minimum scale on which a small gas perturbation will grow due to gravity overcoming the pressure gradient. As long as the Compton scattering of the CMB on the residual free electrons after cosmic recombination kept the gas temperature coupled to that of the CMB, the Jeans mass was constant in time. However, at the gas temperature decoupled from the CMB temperature and the Jeans mass began to decrease with time as the gas cooled adiabatically. Any overdensity on a scale more massive than the Jeans mass at a given time can begin to collapse, due to a lack of sufficient pressure. However, the Jeans mass is related only to the evolution of perturbations at a given time. When the Jeans mass itself varies with time, the overall suppression of the growth of perturbations depends on a time-averaged Jeans mass.
Gnedin & Hui (1998) defined a “filtering mass” that describes the highest mass scale on which the baryonic pressure still manages to suppress the linear baryonic fluctuations significantly. Gnedin (2000) suggested, based on a simulation, that the filtering mass also describes the largest halo mass whose gas content is significantly suppressed compared to the cosmic baryon fraction. The latter mass scale, in general termed the “characteristic mass”, is defined as the halo mass for which the enclosed baryon fraction equals half the mean cosmic fraction. Thus, the characteristic mass distinguishes between gas-rich and gas-poor halos. Many semi-analytical models of dwarfs galaxies often use the characteristic mass scale in order to estimate the gas fraction in halos (e.g., Bullock, Kravtsov, & Weinberg, 2000; Benson et al., 2002a, b; Somerville, 2002). Theoretically this sets an approximate minimum value on the mass that can still form stars.
4.2.1 Prediction from linear theory
In linear theory the filtering mass, first defined by Gnedin & Hui (1998), describes the highest mass scale on which the baryon density fluctuations are suppressed significantly compared to the dark matter fluctuations. Naoz & Barkana (2007) included the fact that the baryons have smoother ICs than the dark matter (see Naoz & Barkana, 2005) and found a lower value of the filtering mass (by a factor of , depending on the redshift). Following Naoz & Barkana (2007), the filtering scale (specifically, the filtering wavenumber ) is defined by expanding the ratio of baryonic to total density fluctuations to first order in :
where is the wavenumber, and and are the baryonic and total (i.e., including both baryons and dark matter) density fluctuations, respectively. The parameter (a negative quantity) describes the relative difference between and on large scales (Naoz & Barkana, 2007), i.e.,
where , (see also Barkana & Loeb, 2005). The ratio is independent of , and its magnitude decreases with time approximately , since is roughly constant and is dominated by the growing mode (see figure 1 top panel in Naoz & Barkana, 2007).
The filtering mass is defined from simply as:
where is the mean matter density today. This relation is one eighth of the definition in Gnedin (2000) (who also used a non-standard definition of the Jeans mass). In figure 4 (bottom panel) we show the filtering mass (solid curve) resulting from eq. (11), as calculated in Naoz & Barkana (2007) (see also their figure 3).
For each of the models we calculate the filtering mass as described here, assuming the model’s initial conditions. Since the simulation is limited in box size, all of the perturbations on large scales are effectively frozen in the simulation. Therefore, we do not extract directly from the simulations, but instead calculate it based on the initial conditions as , where the subscript ”in” refers to initial. Thus, for example, for the E- case, . Figure 4 (bottom panel) shows the analytical results of the filtering mass for the fid calculation, the mean approximation and E- (solid, dashed and dotted curves, respectively). Since the fid calculation is the most consistent calculation, we compare the two other models to it.
The filtering mass represents the competition between gravity and pressure, as it measures the largest scale at which pressure has had a significant overall effect on halo formation. Since it measures an integrated effect over the formation, this mass scale is also very sensitive to the evolution history and the initial conditions (as shown in Naoz & Barkana, 2007). In the mean model, the temperature fluctuations are greatly overestimated on all relevant scales (see Naoz & Barkana, 2005), while in reality the coupling to the CMB (in the fid model) keeps the temperature fluctuations highly suppressed for some time after recombination. Moreover, as mentioned in section 3.1 (and see also Appendix C), we do not include explicitly the effect of initial temperature fluctuations in the simulations. However, the temperature fluctuations from higher redshifts influence the baryon density at the initial redshift (see figure 1) and suppress the baryon density on small scales. As demonstrated in Naoz & Barkana (2007) the system remembers the initial conditions. In other words, the initially enhanced filtering mass (compared to the fid model) helps maintain a higher filtering mass even at moderately low redshift.
In the E- model, the baryon perturbations start out much higher than in the other models, so one might expect that the final baryon fraction in halos would tend to be higher as well; here, however, it is important to separate two issues. The high initial baryon perturbations in the E- model are present at all scales, so they affect even high-mass halos that are unaffected by pressure. This can explain why the simulation with the E- ICs produced the highest baryon fraction in high-mass halos (see the top panel of Figure 4). However, when we consider the differences between large and small scales, the high baryon perturbations produce a large pressure term, increasing the effect of pressure relative to gravity and producing a higher filtering mass in the E- model than in the fid model. Note that the filtering mass is particularly sensitive to the importance of pressure at the very highest redshifts (above 100), since at lower redshifts the gas cools and the Jeans mass decreases, reducing the contribution of these redshifts to the final filtering mass.
We note that in Naoz & Barkana (2007) the calculation of the filtering mass in the fiducial model was compared to the time integrated filtering mass in a model that assumes the mean speed of sound model, neglects the factor, and starts out with initial conditions as in the E- model. Here, we have separated our discussion into several different cases.
4.2.2 The non-linear characteristic mass
There is no apriori reason to think that the filtering mass can also accurately describe properties of highly non-linear, virialized objects. For halos, Gnedin (2000) defined a characteristic mass for which a halo contains half the mean cosmic baryon fraction . In his simulation he found the mean gas fraction in halos of a given total mass , and fitted the simulation results to the following formula:
where is the gas fraction in the high-mass limit. In this function, a higher causes a sharper transition between the high-mass (constant ) limit and the low-mass limit (assumed to be ). Gnedin (2000) found a good fit for , with a characteristic mass that in fact equaled the filtering mass by his definition. By our definition, the claim from Gnedin (2000) is that .
Naoz, Barkana & Mesinger (2009) found that, given their errors, the filtering mass from linear theory is consistent with the characteristic mass fitted from the simulations, for two (pre-reionization) scenarios that they tested: the NoUV case (i.e., no stellar heating) and the Flash case (i.e., after a sudden flash of stellar heating). For clarity, we emphasize that this statement () refers to our definition of in equation (11).
The characteristic mass is essentially a non-linear version of the filtering mass, and so it also measures the competition between gravity and pressure. At high masses, where pressure is unimportant, , while the low mass tail is determined by the suppression of gas accretion caused by high baryonic pressure.
4.2.3 Comparison between the simulation and the theoretical predictions
A major conclusion of the simulation results is that different ICs result in different gas fractions in the final halos. Specifically, we measure these differences through the characteristic mass at various redshifts. varies for different ICs. We determine for each simulation output the characteristic mass and the parameter using a two-dimensional fit to equation (12), with separately fixed to equal the average of the highest few mass bins (see Appendix B for a complete description of the fitting process, together with the errors). In figure 4 we show , and , for all the simulated cases. The characteristic mass clearly depends on the initial conditions, with the mean model and E- model both yielding gas suppression at systematically higher halo masses then for the fid model. The parameter shows a less clear pattern with redshift, but it is generally lowest for the fid model. Overall, the most important implication is that the gas fraction in halos is highly sensitive to the assumed initial conditions.
Comparing to linear theory allows us to understand some of these results. As noted in section 4.2.1, we calculated the filtering mass from linear theory for each of the ICs, and the linear calculation allows us to understand the relative importance of pressure in the various IC models, at least during the linear evolution. Although the simulation results come from non-linear, viralized halos, we find an approximate agreement (typically to within ) between the filtering mass, as defined here and in our previous work (Naoz, Barkana & Mesinger, 2009; Naoz & Barkana, 2007), and the characteristic mass as measured in the simulation, for all the models. In particular, the relative sizes of among the various models, and the slow decline of all the characteristic masses with time, are well matched by the corresponding values predicted from linear theory. This close match can be understood from the fact that while both gravity and pressure increase during the non-linear evolution, their relative strength only changes by a relatively small factor as a halo undergoes non-linear collapse and virialization. Halos in which pressure had a large effect during the early, linear evolution stage, keep sufficient pressure to maintain the suppressed baryon content all through the final collapse. On the other hand, in more massive halos in which gravity overcame pressure early on, the baryons keep up with the collapse of the dark matter and the pressure never has a major role.
For the E- alternative model, we find that the resulting characteristic mass is higher than the result in the fid model. Specifically, at we find M and ). This can be understood since setting the gas fluctuations to be equal to the dark matter’s means that the pressure of the gas is higher compared to the fid model. As can also be seen from comparison to linear theory, the system retains the memory of the pressure, due to the time integrated nature of the filtering mass. Therefore, the higher pressure translates to a higher filtering/characteristic mass.
The mean approximation starts with effectively smoother ICs than in the fid model ( underestimate of the small-scale baryon overdensity). Thus, the baryonic components lag behind the dark matter collapse, and the pressure is always overestimated for a given baryon overdensity (due to the overestimated temperature fluctuations), resulting in a lower gas fraction for any given halo mass, i.e., the characteristic mass is higher than in the fid model. Specifically, at we find M and . This can be compared with M and for the fid ICs.
Recently, Hoeft et al. (2006) and Okamoto, Gao, & Theuns (2008) showed that the characteristic mass scale does not agree with the Gnedin & Hui (1998) filtering mass in the low-redshift, post-reionization regime. However, it is important to note that at these low redshifts, the heating/cooling and other feedback mechanisms are complex and highly inhomogeneous, so that the “filtering mass” calculated from linear theory is not really precisely defined, and the comparison of the linear and non-linear results cannot really be considered a direct and precise test. In contrast, Naoz, Barkana & Mesinger (2009) found that the filtering mass gives a good approximation to the characteristic mass, even in the presence of a ”flash” heating event (see also Mesinger, Bryan & Haiman, 2006) that is physically somewhat contrived but allows for a clear comparison of the linear and non-linear results.
Summarizing our results, we find a good agreement between the characteristic mass and the filtering mass in all the models. Figure 4 shows the best fitted parameters at various redshifts for and , and our value for , for all models (the 1- () confidence regions are listed in table 1). It is important to emphasize that in this statement we are referring to our definition in equation (11), which is one eighth of the original definition which Gnedin (2000) claimed was a good fit to the characteristic mass. While we have been careful to select halos with at least 500 particles, based on the results of Naoz, Barkana & Mesinger (2009), we do not have the even higher mass resolution needed to perform a resolution convergence test as they did. Our main conclusion is that at least in the redshift range the filtering mass provides a fairly good estimate for the characteristic mass. This extends the redshift range of the agreement between the filtering mass and the characteristic mass found in Naoz, Barkana & Mesinger (2009) (). Another significant result from this agreement is that previous work (either analytical, semi-analytical, or using simulations) that used the filtering or characteristic mass without accounting for the correct initial conditions resulted in inaccurate results. This is due to the significant (factor of 2–3) variation among the predictions of the filtering/characteristic mass in the various models. Since this mass scale is of prime importance in early structure formation it is imperative to calculate it accurately.
We have used three-dimensional hydrodynamical simulations to investigate the effect of different initial conditions on the gas fraction in halos in the early universe. Specifically, we studied the minimum “gas-rich” mass defined to have half of the mean cosmic baryon fraction. We tested three different models for the initial conditions (see text for more details)
”fid” ICs; this model is based on the linear evolution from Naoz & Barkana (2005), which allows the baryonic speed of sound to spatially vary as a result of the Compton scattering with the CMB.
”E-” ICs; in this model, the linear evolution from Naoz & Barkana (2005) is modified to match a common assumption in the literature, where the linear initial overdensity of the baryons is taken to be equal to that of the DM, i.e., , while conserving from the fid model.
”mean ICs”; this model assumes that the baryonic speed of sound is spatially uniform. Although Naoz & Barkana (2005) showed that this assumption yields an inaccurate evolution of the baryon density and temperature perturbations, it is still often used in codes that generate initial conditions for simulations.
For all of the tests we used a total of particles of dark matter and baryons with a box size of Mpc, starting at .
There are two major findings from the analysis we present here. The first, shown throughout the paper, is the importance of assuming the correct initial conditions, both for analytical calculations and numerical simulations. Structure formation (both in the linear and non-linear regime) and halo gas fractions are very sensitive to the initial conditions even at relatively low redshifts (). The second major finding is the apparent agreement between the filtering mass and the characteristic mass (to within ). This suggests, as a broader implication, that one can use linear theory in order to predict the overall trend of highly non-linear behavior (at least in the case of determining the gas fraction of halos).
The the fiducial calculation, which was presented in Naoz & Barkana (2005), follows the time evolution of the linear overdensities correctly. However, the other ICs produce different results for the baryonic structure formation. For instance, the non-linear power spectrum (fig. 3) shows that the system still remembers its initial condition differences even at redshift 15. In particular, the model underestimates the non-linear baryonic fluctuations by about while the E- model overestimates them by on small scales.
The mean approximation and the E- model are often used to set the initial conditions in simulations, e.g., the CMBFAST code (Seljak & Zaldarriaga, 1996) assumes the mean approximation while Eisenstein & Hu (1999) is used with the E- assumption. We have shown that the non-linear evolution is very sensitive to the initial conditions (figure 3) and they affect the gas fraction in small halos down to redshift (figure 4). Our results emphasize the importance of the differences between the dark matter and baryons and of the spatial sound speed fluctuations, in both the linear calculation and the initial conditions of the simulations.
It is important to emphasize that although Compton heating is not included in the Gadget code that we used in this analysis (Gadget-2), the fiducial calculation still describes fairly well the non-linear behavior. Actually, the Compton heating contribution to the heating of the gas in non-linear objects is negligible compare to the adiabatic heating due to the gravitational collapse (see Appendix C). Also, as noted above, much of the contribution to the filtering mass comes from the highest redshifts, above our simulation starting redshift of 99, since the Jeans mass is highest then and so the pressure has the greatest impact at that early time.
In each simulation, we calculated the characteristic mass for which a halo keeps most of its baryons (eq. 12). We found that the fid calculation gives the lowest value, which suggests that with these correct ICs, pressure plays only a moderate role in galaxy formation. In particular, the characteristic mass of M is significantly below the minimum mass for molecular hydrogen cooling, so the gas content is not strongly suppressed even in the smallest star-forming halos. In other words we find that before significant heating took place the baryon fraction in halos is (eq. 12 with M and )
The other alternative models give incorrect higher value for the characteristic mass, closer to the minimum mass for forming stars. Even with the fid ICs, pressure does strongly limit the amount of gas in minihalos below the molecular hydrogen cooling mass. We note that this value of M assumes adiabatic evolution, in particular with no stellar heating. This value is consistent with the results of Naoz, Barkana & Mesinger (2009) for a somewhat higher redshift range.
We find that the theoretical linear filtering mass (as defined in section 4.2.1) is in fairly good agreement with the characteristic mass. This finding is true for all the models tested here, throughout a significant redshift range, so this may imply more generally a close relation between linear theory and non-linear halo formation. In addition, this is consistent with the findings by Naoz, Barkana & Mesinger (2009) from AMR simulations, where the filtering mass and the characteristic mass agreed in the ”E-” model, even when a sudden heating was introduced.
Finally, we emphasize that our results are valid only in the pre-reionization era. At the end of the reionization, Mesinger & Dijkstra (2008) concluded that the characteristic mass is likely to be close to the atomic-cooling threshold of , which is also close to the values found by Hoeft et al. (2006) and Okamoto, Gao, & Theuns (2008).
Recently Tseliakhovich & Hirata (2010) argued that the initial velocity difference between the baryons and dark matter after recombination has not been fully accounted for, because of a higher-order contribution that is not included in the linear theory approach. They estimated this higher-order effect within the mean approximation and found that it causes an additional suppression of the small-scale power spectrum, in turn affecting the formation of the first structures. This effect should be further investigated as in our detailed approach here, although this would be more difficult (analytically, it is a higher-order and anisotropic term, and to simulate it directly would require starting at quite high redshifts).
We thank the anonymous referee for useful and helpful comments. We thank Ikkoh Shimizo for supplying the code for calculating the non-linear power spectrum. We also thank Nick Gnedin, Andrey Kravtsov, Matt McQuinn and Michele Trenti for useful discussions. SN also expresses special thanks to Yoram Lithwick for interesting discussions. The authors acknowledge financial support by the Grants-in-Aid for Young Scientists (S) 20674003 by the Japan Society for the Promotion of Science. SN acknowledges NASA ATP grant NNX07AH22G and in part the German-Israeli Project Cooperation (DIP) grant STE1869/1-1.GE625/15-1 and the generous support of the National Post Doctoral Award Program for Advancing Women in Science (Weizmann Institute of Science). RB is grateful for the support of Israel Science Foundation grant 823/09.
- Abel, Bryan, & Norman (2002) Abel T., Bryan G. L., Norman M. L., 2002, Sci, 295, 93
- Barkana & Loeb (2002) Barkana R., Loeb A., 2002, ApJ, 578, 1
- Barkana & Loeb (2004) Barkana R., Loeb A., 2004, ApJ, 609, 474
- Barkana & Loeb (2005) Barkana R., Loeb A., 2005, MNRAS, 363, L36
- Benson et al. (2002a) Benson A. J., Frenk C. S., Lacey C. G., Baugh C. M., Cole S., 2002a, MNRAS, 333, 177
- Benson et al. (2002b) Benson A. J., Lacey C. G., Baugh C. M., Cole S., Frenk C. S., 2002b, MNRAS, 333, 156
- Bromm, Coppi, & Larson (2002) Bromm V., Coppi P. S., Larson R. B., 2002, ApJ, 564, 23
- Bromm, Coppi, & Larson (1999) Bromm V., Coppi P. S., Larson R. B., 1999, ApJ, 527, L5
- Bromm & Loeb (2006) Bromm V., Loeb A., 2006, ApJ, 642, 382
- Bromm et al. (2009) Bromm V., Yoshida N., Hernquist L., McKee C. F., 2009, Natur, 459, 49
- Bouwens et al. (2010) Bouwens R. J., et al., 2010, ApJ, 709, L133
- Bullock, Kravtsov, & Weinberg (2000) Bullock J. S., Kravtsov A. V., Weinberg D. H., 2000, ApJ, 539, 517
- Clark et al. (2010) Clark P. C., Glover S. C. O., Klessen R. S., Bromm V., 2010, arXiv, arXiv:1006.1508
- Ciardi et al. (2006) Ciardi B., Scannapieco E., Stoehr F., Ferrara A., Iliev I. T., Shapiro P. R., 2006, MNRAS, 366, 689
- Duffy et al. (2010) Duffy A. R., Schaye J., Kay S. T., Dalla Vecchia C., Battye R. A., Booth C. M., 2010, arXiv, arXiv:1001.3447
- Eisenstein & Hu (1999) Eisenstein D. J., Hu W., 1999, ApJ, 511, 5
- Fan et al. (2003) Fan X., et al., 2003, AJ, 125, 1649
- Furlanetto & Oh (2006) Furlanetto S. R., Oh S. P., 2006, ApJ, 652, 849
- Gao et al. (2007) Gao L., Yoshida N., Abel T., Frenk C. S., Jenkins A., Springel V., 2007, MNRAS, 378, 449
- Gnedin & Hui (1998) Gnedin, N. Y. & Hui, L. 2004, MNRAS, 296.
- Gnedin (2000) Gnedin N. Y., 2000, ApJ, 542, 535.
- Gnedin, Kravtsov, & Chen (2008) Gnedin N. Y., Kravtsov A. V., Chen H.-W., 2008, ApJ, 672, 765
- Greif et al. (2010) Greif T. H., Glover S. C. O., Bromm V., Klessen R. S., 2010, ApJ, 716, 510
- Hoeft et al. (2006) Hoeft M., Yepes G., Gottlöber S., Springel V., 2006, MNRAS, 371, 401
- Harford, Hamilton, & Gnedin (2008) Harford A. G., Hamilton A. J. S., Gnedin N. Y., 2008, MNRAS, 389, 880
- Iliev et al. (2003) Iliev, I. T., Scannapieco, E., Martel, H., & Shapiro, P. R. 2003, MNRAS, 341, 81
- Iliev et al. (2005) Iliev I. T., Scannapieco E., Shapiro P. R., 2005, ApJ, 624, 491
- Jeans (1928) Jeans, J. H. Cambridge University press, (1928).
- Komatsu et al. (2009) Komatsu E., et al., 2009, ApJS, 180, 330
- Komatsu et al. (2010) Komatsu E., et al., 2010, arXiv, arXiv:1001.4538
- Kuhlen, Madau, & Montgomery (2006) Kuhlen M., Madau P., Montgomery R., 2006, ApJ, 637, L1
- Iye et al. (2006) Iye M., et al., 2006, Natur, 443, 186
- Lin et al. (2006) Lin W. P., Jing Y. P., Mao S., Gao L., McCarthy I. G., 2006, ApJ, 651, 636
- Lin, Liang, & Zhang (2010) Lin L., Liang E., Zhang S. 2010, ScChG, 53, 64
- Ma & Bertschinger (1995) Ma, C. & Bertschinger, E. 1995, ApJ, 455, 7
- Mesinger, Bryan & Haiman (2006) Mesinger A., Bryan G., & Haiman Z. 2006, ApJ, 648, 835
- Mesinger & Dijkstra (2008) Mesinger A., Dijkstra M., 2008, MNRAS, 390, 1071
- McLure et al. (2010) McLure R. J., Dunlop J. S., Cirasuolo M., Koekemoer A. M., Sabbi E., Stark D. P., Targett T. A., Ellis R. S., 2010, MNRAS, 126
- McQuinn et al. (2007) McQuinn M., Lidz A., Zahn O., Dutta S., Hernquist L., Zaldarriaga M., 2007, MNRAS, 377, 1043
- Morales & Wyithe (2009) Morales M. F., Wyithe J. S. B., 2009, arXiv, arXiv:0910.3010
- Naoz & Barkana (2005) Naoz S., Barkana R., 2005, MNRAS, 362, 1047
- Naoz, Noter & Barkana (2006) Naoz S., Noter S., Barkana R., 2006, MNRAS, 373, L98
- Naoz & Barkana (2007) Naoz S., Barkana R., 2007, MNRAS, 377, 667
- Naoz & Barkana (2008) Naoz S., Barkana R., 2008, MNRAS, 385, L63
- Naoz, Barkana & Mesinger (2009) Naoz S., Barkana R., Mesinger A., 2009, MNRAS, 377, 667
- Naoz & Bromberg (2007) Naoz S., Bromberg O., 2007, MNRAS, 380, 757
- Okamoto, Gao, & Theuns (2008) Okamoto T., Gao L., Theuns T., 2008, arXiv, 806, arXiv:0806.0378
- Ricotti, Gnedin, & Shull (2002) Ricotti M., Gnedin N. Y., Shull J. M., 2002, ApJ, 575, 49
- Ricotti, Gnedin, & Shull (2008) Ricotti M., Gnedin N. Y., Shull J. M., 2008, ApJ, 685, 21
- Salvaterra et al. (2009) Salvaterra R., et al., 2009, Natur, 461, 1258
- Seljak & Zaldarriaga (1996) U., Seljak, & M., Zaldarriaga 1996, ApJ, 469, 437; see http://www.cmbfast.org
- Shapiro et al. (2006) Shapiro P. R., Ahn K., Alvarez M. A., Iliev I. T., Martel H., Ryu D., 2006, ApJ, 646, 681
- Shapiro, Iliev, & Raga (2004) Shapiro P. R., Iliev I. T., Raga A. C., 2004, MNRAS, 348, 753
- Somerville (2002) Somerville R. S., 2002, ApJ, 572, L23
- Springel (2005) Springel V., 2005, MNRAS, 364, 1105
- Spergel et al. (2003) Spergel D. N., et al., 2003, ApJS, 148, 175
- Spergel et al. (2007) Spergel D. N., et al., 2007, ApJS, 170, 377
- Springel, Yoshida, & White (2001) Springel V., Yoshida N., White S. D. M., 2001, NewA, 6, 79
- Tanvir et al. (2009) Tanvir N. R., et al., 2009, Natur, 461, 1254
- Trenti et al. (2010) Trenti M., Smith B. D., Hallman E. J., Skillman S. W., Shull J. M., 2010, arXiv, arXiv:1001.5037
- Trenti, Stiavelli, & Michael Shull (2009a) Trenti M., Stiavelli M., Michael Shull J., 2009, ApJ, 700, 1672
- Trenti & Stiavelli (2009b) Trenti M., Stiavelli M., 2009, ApJ, 694, 879
- Trenti & Shull (2010) Trenti M., Shull J. M., 2010, ApJ, 712, 435
- Tseliakhovich & Hirata (2010) Tseliakhovich D., Hirata C., 2010, arXiv, arXiv:1005.2416
- Trenti, Santos, & Stiavelli (2008) Trenti M., Santos M. R., Stiavelli M., 2008, ApJ, 687, 1
- Yamamoto et al. (1997) Yamamoto, K., Sugiyama, N., & Sato, H. 1997, PRD, 56, 7566
- Yamamoto et al. (1998) Yamamoto, K., Sugiyama, N., & Sato, H. 1998, ApJ, 501, 442
- Yoshida (2006) Yoshida N., 2006, NewAR, 50, 19
- Yoshida (2009) Yoshida N., 2009, arXiv, arXiv:0906.4372
- Yoshida et al. (2003) Yoshida N., Sokasian A., Hernquist L., Springel V., 2003, ApJ, 598, 73
- Yoshida et al. (2006) Yoshida N., Omukai K., Hernquist L., Abel T., 2006, ApJ, 652, 6
- Yoshida et al. (2003) Yoshida N., Abel T., Hernquist L., Sugiyama N., 2003, ApJ, 592, 645
- Yoshida et al. (2007) Yoshida N., Oh S. P., Kitayama T., Hernquist L., 2007, ApJ, 663, 687
- Yoshida, Omukai & Hernquist (2008) Yoshida N., Omukai K., Hernquist L., 2008, Sci, 321, 669
- Yoshida, Sugiyama, & Hernquist (2003) Yoshida N., Sugiyama N., Hernquist L., 2003, MNRAS, 344, 481
Appendix A conservation
We have defined the two different models such that they conserve . From linear theory we do not expect the evolution of the mean model to be significantly different from that of the fid model (in terms of halo abundance and total power spectrum). This is indeed the case for the evolution in time of the total fluctuations of the mean model compared to the fid model on large scales (small ), as shown in figure 5 (lower set of thin curves).
A more delicate treatment is needed for the E- model (see section 2.2). In this case, at high redshift (such as the initial ), the baryons are in the process of falling into the DM potential. This results in a faster growth of the total fluctuations compared to the case in which there is a relative velocity between the DM and the baryons (such as in the case of the mean cs and fid models, where the relative velocity for the E- model are negligible); see figure 5 dotted thick curve. At later times, the baryon fluctuations approach the dark matter fluctuations, and the large scale behavior (i.e., on linear scales) deviates from the fid model by less then (see the solid curve in fig. 5).
We also note that we have checked the overall effect of on the main results. We have performed two additional simulations for the E- model, where we increased or decreased by . We found that the calculated is within the fit errors (see appendix B and table 1) at . At , the difference in the best fitted value is below .”
Appendix B Fit properties
For each redshift snapshot for each run we find the characteristic mass and using a two dimensional fit. In figure 6 we consider two example redshifts (high, and low, ) for which we show the binned data points and the resulting fit. In table 1 we show our best fit parameters. We note that we have checked that the fits give consistent results if we lower the condition on the minimum number of particles per halo to 300 (instead of 500). We also note that our determination of relies on an extrapolation (via the fit) below our simulations’ resolution limit
The parameter in equation (12) is an average of the gas fraction values in the few highest mass bins. In our simulation the high-end tail of the masses has large scatter in the estimated gas fraction because of the low number of halos (each bin among the last 3 or 4 in figure 6 represents just 1 or 2 halos), thus we have to average over this scatter to get a reasonable result. This scatter is in part a result of assuming that the halos are spherical, and thus halos that are undergoing a major merger deviate greatly from a spherical shape and are treated inaccurately in our analysis. We have tested the resulting when taking a linking parameter of , which indeed resulted in more high-mass halos, but in any case was consistent with the value of we found with the linking parameter. Thus, in this paper, we use the standard value of .
As expected at high redshift, where we have fewer halos, the errors become quite large. We also tried, following Naoz, Barkana & Mesinger (2009), to bin the data and to perform the fit for the binned data with the weight for each bin. For the redshifts for which we had more than halos we got that the binned analysis gave results within the non-binned fit errors,and with comparable errors.
We also tried the approach of taking to be a free parameter, but this produced very problematic fits222Naoz, Barkana & Mesinger (2009) also found that treating as a free parameter was unproductive. . This is mainly because of the large scatter at the high mass end, so that a three-parameter fit could not strongly constrain the parameter values. We also note the fact that is lower than the mean cosmic fraction , by about - for the fid and mean models, and for the E- model (see figure 4 top panel). The result in the fid and mean models may reflect the real suppression of the large-scale baryon fluctuations in these models; the difference in linear theory is at (Naoz & Barkana, 2007), but the non-linear evolution may increase this effect. The discrepancy in the E- model may reflect a limitation of the simulation; we note that in Naoz, Barkana & Mesinger (2009) was also lower than and even lower by from our results at the overlapping redshifts (where we compare the E- model in both cases). This might be due to the fact that gas shocks in AMR are sharper than in Gadget simulations, and thus AMR may produce a more realistic gas profile, although the result is still below the universal cosmic baryon fraction (Lin et al., 2006). In our simulation, going to a larger radii can result in a more realistic value, but we used for consistency with the common definition.
Appendix C Heating of non-linear halos
The fiducial model follows correctly the baryon density and temperature perturbations due to Compton scattering on the residual free electrons after recombination. While this is fully incorporated in our fid ICs, our simulation does not take into account Compton heating. Below we show that for non-linear objects the heating is actually negligible compared to the adiabatic heating due to the gravitational collapse of baryons into the dark matter potential wells. Therefore, it is sufficient to include Compton heating in the linear stage only.
The heating of the gas due to Compton heating from the CMB (Naoz & Barkana, 2005) during the free-fall time of gravitational collapse is
where is the Thomson scattering cross section, is the photon energy density, and are the CMB and gas temperature and is the electron fraction out of the total number density of gas particles at time .
The virial theorem gives a relation in collapsed objects between the thermal energy and the gravitational energy , i.e., . Thus, for a halo mass with virial radius the thermal energy can be expressed as:
For all relevant redshifts and mass scales we find that . Therefore, neglecting the contribution of the Compton heating during the non-linear evolution is justified. However, as we have shown, neglecting the Compton heating in the linear evolution and in the initial conditions leads to inaccurate values for the gas fraction in halos.