text:
# A Dynamical Model for Hot Gas in the Galactic Center

## 1 Introduction

Chandra observations of the center of the Milky Way reveal diffuse gas within several parsecs of the central massive black hole, Sgr A\* (Baganoff et al. 2003). This gas undoubtedly originates from the interaction of the strong stellar winds produced by the several dozen massive stars in the central parsec star cluster (e.g., Krabbe et al. 1991; Najarro et al. 1997). These stellar winds, and the associated hot X-ray emitting gas, are believed to be the primary reservoir of material for accretion onto the central black hole (e.g., Melia 1992).

In this paper I present a model for the dynamics of the observed hot gas on scales of $`0.011`$ pc in the Galactic Center. I am motivated by several considerations. First, the rate at which Sgr A\* accretes surrounding gas is usually estimated using the Bondi accretion formula. The resulting accretion rate, $`10^510^6\mathrm{M}_{}\mathrm{yr}^1`$ (e.g., Melia 1992; Baganoff et al. 2003), is much less than the total mass loss rate by stars in the central parsec, $`10^3\mathrm{M}_{}\mathrm{yr}^1`$ (Najarro et al. 1997). Nearly all of the mass lost by stars must therefore be driven out of the Galactic Center in a wind (e.g., Chevalier 1992). To accurately model the gas gravitationally captured by Sgr A\*, one should also self-consistently account for the dynamics of the unbound wind. Moreover, during most epochs, Chandra observes an extended X-ray source coincident with Sgr A\*, which has a size $`R_B1^{\prime \prime }`$ (Baganoff et al. 2003), where $`R_BGM/c_s^2`$ is the Bondi accretion radius for gas of sound speed $`c_s`$ around a black hole of mass $`M`$. In a previous paper (Quataert 2002; see also Yuan et al. 2002), I argued that this extended source is due to thermal emission from hot gas at $`R_B`$, but I did not present a detailed model for the dynamics of this gas, nor did I quantitatively model the observed surface brightness profile. If the thermal emission interpretation of the extended source is correct, it would be direct evidence for gas being gravitationally captured by the central black hole, confirming that $`10^5\mathrm{M}_{}\mathrm{yr}^1`$ is flowing in on scales of $`R_B`$. The significance of this inference motivates a better model for the dynamics of the hot gas observed by Chandra.

In the next section I incorporate stellar mass loss as a source term in the hydrodynamic equations and calculate the dynamics of both accreting and outflowing gas assuming spherical symmetry. This extends previous work on a Galactic Center wind (e.g., Chevalier 1992) to incorporate the effects of the central black hole. The spherical assumption is relatively simplistic since a small number of stars dominate the mass supply in the central parsec star cluster. I show, however, that this model reproduces the Chandra observations well. Melia and collaborators (Coker & Melia 1997; Rockefeller et al. 2003) have presented 3D simulations that address some of the issues considered here.

## 2 The Fate of Stellar Winds

The interaction of multiple stellar winds leads to shocks which heats the gas to X-ray emitting temperatures. The dynamics of the hot gas can be modeled by incorporating stellar winds as a source of mass and energy in the equations of hydrodynamics. Assuming spherical symmetry, the resulting equations are given by (e.g., Holzer & Axford 1970)

$$\frac{\rho }{t}+\frac{1}{r^2}\frac{}{r}r^2\rho v=q(r),$$
(1)
$$\rho \frac{dv}{dt}=\frac{p}{r}\rho \frac{GM}{r^2}q(r)v,$$
(2)
and

$$\rho T\frac{ds}{dt}=q(r)\left[\frac{v^2}{2}+\frac{v_w^2}{2}\frac{\gamma }{\gamma 1}c_s^2\right],$$
(3)
where $`\rho `$, $`v`$, $`c_s`$, and $`s`$ are the mass density, radial velocity, isothermal sound speed, and entropy per unit mass, respectively; $`M=3.6\times 10^6M_{}`$ is the mass of the black hole which dominates the gravity on the scales of interest. For the densities and temperatures appropriate to the Galactic Center, radiative cooling is negligible and so has been dropped in eq. (3). In eqs. (1)-(3), $`q(r)`$ is the stellar mass loss rate per unit volume and $`v_w^2/2`$ is the rate of energy injection per unit mass from stellar winds with velocity $`v_w`$. The total rate of mass injection is given by $`\dot{M}_w=4\pi r^2q(r)`$. Incorporating mass loss as a source term eliminates the need to specify boundary conditions on the density and temperature of gas at an arbitrary ’fiducial’ radius, as is required in Bondi accretion and Parker wind models. It should be noted that equations analogous to those above have been used to study gas flow in several other environments, such as ’galactic winds’ from elliptical galaxies (e.g., Mathews & Baker 1971) and winds from star clusters without black holes (e.g., Canto et al. 2000).

The source terms in the above equations are well constrained by IR spectroscopy of the Galactic Center, which reveals a cluster of massive stars within $`10`$” of the black hole (e.g., Genzel et al. 2003). These include blue supergiants with mass loss rates $`10^4\mathrm{M}_{}\mathrm{yr}^1`$ and wind speeds $`v_w6001000`$ km s-1 (e.g., Najarro et al. 1997). The total stellar mass loss is $`10^3\mathrm{M}_{}\mathrm{yr}^1`$ and is dominated by IRS 13E which is $`3.5^{\prime \prime }`$ from Sgr A\* on the sky. There is an additional cluster of massive stars much closer to the black hole, namely those whose orbits have recently been measured (e.g., Schödel et al. 2002; Ghez et al. 2003a). Spectroscopy of one such star, S0-2, suggests that it is a main sequence O/B star (Ghez et al 2003b; Eisenhauer et al. 2003), in which case its mass loss rate is probably much smaller than that of the evolved stars further from the black hole. In what follows I neglect mass loss from the closer-in star cluster; further observations are required to check this assumption.

I solve equations (1)-(3) to determine the fate of hot gas in the Galactic Center. I choose model parameters based on the observations described above; for standard parameters I take $`v_w1000`$ km s-1 and $`\dot{M}_w10^3\mathrm{M}_{}\mathrm{yr}^1`$. The biggest uncertainty is how to model the spatial distribution of mass loss in this simplified one-dimensional calculation. Since the observed mass losing stars are located several arcsec from the black hole, I set $`q(r)r^\eta `$ for $`rϵ[2^{\prime \prime },10^{\prime \prime }]`$, and $`q(r)=0`$ otherwise. The local mass injection rate is given by $`d\dot{M}_w/d\mathrm{ln}rr^{\eta +3}`$ so that $`\eta =0`$ corresponds to mass injection that is concentrated at large radii while $`\eta =3`$ corresponds to equal mass injection per decade in radius. Modest variations about this choice of $`q(r)`$ yield similar results to those described below.

After several sound crossing times, the solution of equations (1)-(3) settles into a steady state in which gas in the inner region is captured and flows in towards the black hole, while gas further away is blown out of the system in a wind. Figures 1 & 2 show the steady state radial velocity, temperature, and density as a function of distance from the black hole for the above parameters with $`\eta =0,2,3`$. The total accretion rate through the inner boundary in the three solutions ranges from $`0.00150.03\dot{M}_w10^5\mathrm{M}_{}\mathrm{yr}^1`$; the majority of the gas, $`\dot{M}_w10^3\mathrm{M}_{}\mathrm{yr}^1`$, is driven out of the central star cluster. As Fig. 1 shows, the division between inflowing and outflowing gas occurs quite close to the black hole, at $`23^{\prime \prime }`$;111Asymptotically the solutions in Figs. 1 and 2 approach that of Bondi accretion at small radii (Bondi & Hoyle 1944) and that of a thermally driven wind at large radii (e.g., Parker 1960). in 1D this separation occurs at a stagnation point where $`v=0`$. The temperature near the stagnation point, at $`110^{\prime \prime }`$, is set by the stellar wind velocity; for $`vc_s`$, the steady state solution to equation (3) is $`c_s^2v_w^2/5`$, which yields a characteristic temperature of $`1`$ keV, as is observed by Chandra (Baganoff et al. 2003).

The density profile of the gas is shown in Fig. 2. Also shown are two Chandra measurements (Baganoff et al. 2003); the first is an inferred density of $`27`$ cm-3 in the central 10” from the diffuse thermal emission in the Galactic Center. The second is an inferred density of $`130`$ cm-3 in the central 1.5” from the extended X-ray source coincident with Sgr A\* (discussed more below). The agreement between the models and the Chandra observations shows that, although the spherically symmetric approximation has its limitations, it captures the overall dynamics of hot gas in the Galactic Center reasonably well. It is worth stressing that in our model, the two Chandra observations probe gas with very different dynamics: the $`10^{\prime \prime }`$ observation probes the majority of the gas that is being driven out of the central star cluster away from the black hole, while the $`1.5^{\prime \prime }`$ observation probes the small fraction of the gas that is gravitationally captured by the black hole.

Figure 3 shows the surface brightness profile within 3” of Sgr A\* derived from the density and temperature profiles shown in Fig. 1-2. Also shown is the surface brightness profile observed by Chandra coincident with Sgr A\*, and that of a nearby point source (as an indication of the Chandra point spread function; see Baganoff et al. 2003). Aside from epochs when the source flares (Baganoff et al. 2001; Porquet et al. 2003), the X-ray source coincident with Sgr A\* is significantly extended with respect to nearby point sources, although there are uncertainties in the surface brightness profile of Sgr A\* because X-ray scattering by dust grains along the line of sight can make the profile somewhat more extended (see Tan & Draine 2003, who estimate that dust scattering of an unresolved source can, at most, account for $`50\%`$ of the extended emission coincident with Sgr A\*).

The models presented here of thermal emission produced by gas gravitationally captured by Sgr A\* reproduce the extended X-ray source reasonably well. Naturally, models with steeper density profiles (e.g., $`\eta =3`$) produce steeper surface brightness profiles more in accord with the Chandra observations. I have found it difficult to fully reproduce the rapid drop in surface brightness that is observed. Models with such a steep surface brightness profile either overpredict the density by a factor of few or underpredict the temperature by a comparable factor. There are several possible explanations for this discrepancy. For example, it may be due to limitations of the spherically symmetric dynamical model. In addition, the timescale for electrons and protons to come into thermal equilibrium with each other is $`200n_{100}^1T_1^{3/2}`$ years, where $`n_{100}=n/100\mathrm{cm}^3`$ and $`T_1=T_e/1\mathrm{keV}`$; this is longer than the characteristic flow time at 1”, $`R/v_w40`$ years. Thus it is possible that the electrons and protons do not have the same temperature on the scales observed by Chandra (if, e.g., they are shock heated to different temperatures, as is typically observed in supernova shocks).222The good agreement between the observed electron temperature and the temperature expected from shocked stellar winds suggests that this is not likely to be a very large effect; even a factor of $`1.52`$ difference would, however, modify the comparisons made here. In this case the observed electron temperature would not be a good proxy for the total pressure, as is assumed in the models considered here. In principle, this could be tested by observing lines from ions such as iron, oxygen, or nitrogen, as is done to probe the temperature structure of supernova remnants (e.g., Vink et al. 2003). In practice this may not be possible because many of the relevant lines are in the soft X-rays which are heavily absorbed towards the Galactic Center, and because the lines will be too narrow to be resolved by Chandra (Astro-E2 may help with the latter problem).

## 3 Discussion

The models presented here describe the dynamics of the hot gas produced by shocked stellar winds in the Galactic Center, assuming for simplicity spherical symmetry. They quantitatively account for the observed level of diffuse X-ray emission in the central parsec, predicting an electron density on $`10^{\prime \prime }`$ scales of $`2030`$ cm-3, in good agreement with Chandra observations (Fig. 2). This emission is produced by gas that is not bound to the black hole and is being thermally driven out of the central star cluster in a wind. This wind can have important dynamical effects on the surrounding interstellar medium (e.g., Yusef-Zadeh & Wardle 1993); it may also be an important source of mass for the thermal X-ray emitting ’lobes’ observed symmetrically around Sgr A\* by Chandra (Morris et al. 2002).

In our models, a few percent of the mass supplied by stellar winds to the central parsec is gravitationally captured by Sgr A\*, implying an accretion rate (at large radii) of $`10^5\mathrm{M}_{}\mathrm{yr}^1`$. As this gas moves in towards the black hole it is compressed, resulting in an increase in the gas density and X-ray surface brightness close to Sgr A\* that are in reasonable agreement with Chandra observations (Fig. 2 & 3). I suggest that this agreement provides strong evidence that Chandra has directly observed gas being gravitationally captured by Sgr A\*, confirming one of the long-standing predictions of theoretical accretion models (e.g., Melia 1992; Narayan et al. 1995). There is also evidence from the linear polarization of mm emission from Sgr A\* that the density close to the black hole is much less than a straightforward extrapolation of the Bondi accretion rate to small radii (e.g., Bower et al. 2003). This is in accord with theoretical predictions that very little of the gas captured at large radii actually accretes onto the black hole (e.g., Blandford & Begelman 1999).

###### Acknowledgements.

I thank Fred Baganoff, Bruce Draine, and Jonathan Tan for useful discussions; conversations with Anatoly Spitkovsky on numerical methods for solving PDEs were particularly useful. This work was supported in part by NSF grant AST 0206006, NASA Grant NAG5-12043, and an Alfred P. Sloan Fellowship.


format: text

text:
# Unusually Large Fluctuations in the Statistics of Galaxy Formation at High Redshift

## 1. Introduction

Recent observations of the cosmic microwave background (Spergel et al., 2003) have confirmed the notion that the present large-scale structure in the universe originated from small-amplitude density fluctuations at early cosmic times. Due to the natural instability of gravity, regions that were denser than average collapsed and formed bound halos, first on small spatial scales and later on larger and larger scales. At each snapshot of this cosmic evolution, the abundance of collapsed halos, whose masses are dominated by cold dark matter, can be computed from the initial conditions using numerical simulations and can be understood using approximate analytic models (Press & Schechter, 1974; Bond et al., 1991). The common understanding of galaxy formation is based on the notion that the constituent stars formed out of the gas that cooled and subsequently condensed to high densities in the cores of some of these halos (White & Rees, 1978).

The standard analytic model for the abundance of halos (Press & Schechter, 1974; Bond et al., 1991) considers the small density fluctuations at some early, initial time, and attempts to predict the number of halos that will form at some later time corresponding to a redshift $`z`$. First, the fluctuations are extrapolated to the present time using the growth rate of linear fluctuations, and then the average density is computed in spheres of various sizes. Whenever the overdensity (i.e., the density perturbation in units of the cosmic mean density) in a sphere rises above a critical threshold $`\delta _c(z)`$, the corresponding region is assumed to have collapsed by redshift $`z`$, forming a halo out of all the mass that had been included in the initial spherical region. In analyzing the statistics of such regions, the model separates the contribution of large-scale modes from that of small-scale density fluctuations. It predicts that galactic halos will form earlier in regions that are overdense on large scales (Kaiser, 1984; Bardeen et al., 1986; Cole & Kaiser, 1989; Mo & White, 1996), since these regions already start out from an enhanced level of density, and small-scale modes need only supply the remaining perturbation necessary to reach $`\delta _c(z)`$. On the other hand, large-scale voids should contain a reduced number of halos at high redshift. In this way, the analytic model describes the clustering of massive halos.

As gas falls into a dark matter halo, it can fragment into stars only if its virial temperature is above $`10^4`$K for cooling mediated by atomic transitions \[or $`500`$ K for molecular $`\mathrm{H}_2`$ cooling; see, e.g., Figure 12 in Barkana & Loeb (2001)\]. The abundance of dark matter halos with a virial temperature above this cooling threshold declines sharply with increasing redshift due to the exponential cutoff in the abundance of massive halos at early cosmic times. Consequently, a small change in the collapse threshold of these rare halos, due to mild inhomogeneities on much larger spatial scales, can change the abundance of such halos dramatically. The modulation of galaxy formation by long wavelength modes of density fluctuations is therefore amplified considerably at high redshift. In this paper we show that this results in major new predictions for high-redshift observations. The implications are particularly significant for cosmic reionization and all observational probes of this epoch.

This paper is organized as follows. In § 2 we quantify the scatter in the statistics of galaxy formation produced by this amplification effect. We first explain in § 2.1 the basic physical ideas and implications using the well-established extended Press-Schechter model. We then present in § 2.2 a simple idea that yields a much more accurate model that fits an array of previous simulations at low redshift. We demonstrate the qualitative correctness of our basic assumptions as well as the quantitative accuracy of our model by matching results from recent simulations at high redshift. Since high-redshift galaxies provide the UV photons that lead to the reionization of the intergalactic medium (hereafter IGM), a large scatter is also expected in the reionization redshift within different regions in the universe. We consider this scatter and the modified character of reionization in § 3.1, and show in § 3.2 that existing numerical simulations do not include fluctuations on sufficiently large scales at high redshift. In § 3.3 we discuss the observational implications of the large cosmic scatter expected at high redshift. Finally, we summarize our main results in § 4.

## 2. Halo Mass Function in Different Environments

### 2.1. Basic Model: Amplification of Density Fluctuations

Galaxies at high redshift are believed to form in all halos above some minimum mass $`M_{\mathrm{min}}`$ that depends on the efficiency of atomic and molecular transitions that cool the gas within each halo. This makes useful the standard quantity of the collapse fraction $`F_{\mathrm{col}}(M_{\mathrm{min}})`$, which is the fraction of mass in a given volume that is contained in halos of individual mass $`M_{\mathrm{min}}`$ or greater. If we set $`M_{\mathrm{min}}`$ to be the minimum halo mass in which efficient cooling processes are triggered, then $`F_{\mathrm{col}}(M_{\mathrm{min}})`$ is the fraction of all the baryons in the universe that lie in galaxies. In a large-scale region of comoving radius $`R`$ with a mean overdensity $`\overline{\delta }_R`$, the standard result is (Bond et al., 1991)

$$F_{\mathrm{col}}(M_{\mathrm{min}})=\mathrm{erfc}\left[\frac{\delta _c(z)\overline{\delta }_R}{\sqrt{2\left[S(R_{\mathrm{min}})S(R)\right]}}\right],$$
(1)
where $`S(R)`$ is the variance of fluctuations in spheres of radius $`R`$, and $`S(R_{\mathrm{min}})`$ is the variance in spheres of radius $`R_{\mathrm{min}}`$ corresponding to the region at the initial time that contained a mass $`M_{\mathrm{min}}`$. In particular, the cosmic mean value of the collapse fraction is obtained in the limit of $`R\mathrm{}`$ by setting $`\overline{\delta }_R`$ and $`S(R)`$ to zero in this expression. Throughout this section our results assume this standard model, known as the extended Press-Schechter model, which we apply to a universe with cosmological parameters matching the latest observations \[specifically, the running index model of Spergel et al. (2003)\]. Whenever we consider a cubic region, we estimate its halo abundance by applying the model to a spherical region of equal volume. Note also that we consistently quote values of comoving distance, which equals physical distance times a factor of $`(1+z)`$.

Our results are based on a simple idea. At high redshift, galactic halos are rare and correspond to high peaks in the Gaussian probability distribution of initial fluctuations. A modest change in the overall density of a large region modulates the threshold for high peaks in the Gaussian density field, so that the number of galaxies is exponentially sensitive to this modulation. This amplification of large-scale modes is responsible for the large statistical fluctuations that we find.

In numerical simulations, periodic boundary conditions are usually assumed, and this forces the mean density of the box to equal the cosmic mean density. The abundance of halos as a function of mass is then biased in such a box (see Figure 1), since a similar region in the real universe will have a distribution of different overdensities $`\overline{\delta }_R`$. At high redshift, when galaxies correspond to high peaks, they are mostly found in regions with an enhanced large-scale density. In a periodic box, therefore, the total number of galaxies is artificially reduced, and the relative abundance of galactic halos with different masses is artificially tilted in favor of lower-mass halos. We illustrate our results for two sets of parameters, one corresponding to the first galaxies and early reionization ($`z=20`$) and the other to the current horizon in observations of galaxies and late reionization ($`z=7`$). We consider a resolution equal to that of state-of-the-art cosmological simulations that include gravity and gas hydrodynamics. Specifically, we assume that the total number of dark matter particles in the simulation is $`N=324^3`$, and that the smallest halo that can form a galaxy must be resolved into 500 particles; Springel & Hernquist (2003) showed that this resolution is necessary in order to determine the star formation rate in an individual halo reliably to within a factor of two. Therefore, if we assume that halos that cool via molecular hydrogen must be resolved at $`z=20`$ (so that $`M_{\mathrm{min}}=7\times 10^5M_{}`$), and only those that cool via atomic transitions must be resolved at $`z=7`$ (so that $`M_{\mathrm{min}}=10^8M_{}`$), then the maximum box sizes that can currently be simulated are $`l_{\mathrm{box}}=1`$ Mpc and $`l_{\mathrm{box}}=6`$ Mpc at these two redshifts, respectively.

At each redshift we only consider cubic boxes large enough so that the probability of forming a halo on the scale of the entire box is negligible. In this case, $`\overline{\delta }_R`$ is Gaussian distributed with zero mean and variance $`S(R)`$, since the no-halo condition $`\sqrt{S(R)}\delta _c(z)`$ implies that at redshift $`z`$ the perturbation on the scale $`R`$ is still in the linear regime. We can then calculate the probability distribution of collapse fractions in a box of a given size (see Figure 2). This distribution corresponds to a real variation in the fraction of gas in galaxies within different regions of the universe at a given time. In a numerical simulation, the assumption of periodic boundary conditions eliminates the large-scale modes that cause this cosmic scatter. Note that Poisson fluctuations in the number of halos within the box would only add to the scatter, although the variations we have calculated are typically the dominant factor. For instance, in our two standard examples, the mean expected number of halos in the box is 3 at $`z=20`$ and 900 at $`z=7`$, resulting in Poisson fluctuations of a factor of about 2 and 1.03, respectively, compared to the clustering-induced scatter of a factor of about 16 and 2 in these two cases.

Within the extended Press-Schechter model, both the numerical bias and the cosmic scatter can be simply described in terms of a shift in the redshift (see Figure 3). In general, a region of radius $`R`$ with a mean overdensity $`\overline{\delta }_R`$ will contain a different collapse fraction than the cosmic mean value at a given redshift $`z`$. However, at some wrong redshift $`z+\mathrm{\Delta }z`$ this small region will contain the cosmic mean collapse fraction at $`z`$. At high redshifts ($`z>3`$), this shift in redshift can be easily derived from eq. 1 to be

$$\mathrm{\Delta }z=\frac{\overline{\delta }_R}{\delta _0}(1+z)\times \left[1\sqrt{1\frac{S(R)}{S(R_{\mathrm{min}})}}\right],$$
(2)
where $`\delta _0\delta _c(z)/(1+z)`$ is approximately constant at high redshifts (Peebles, 1980), and equals 1.28 for our assumed cosmological parameters. Thus, in our two standard examples, the bias is -2.6 at $`z=20`$ and -0.4 at $`z=7`$, and the one-sided $`1\sigma `$ scatter is 2.4 at $`z=20`$ and 1.2 at $`z=7`$.

### 2.2. Improved Model: Matching Numerical Simulations

In this subsection we develop an improved model that fits the results of numerical simulations more accurately. The model constructs the halo mass distribution (or mass function); cumulative quantities such as the collapse fraction or the total number of galaxies can then be determined from it via integration. We first define $`f(\delta _c(z),S)dS`$ to be the mass fraction contained at $`z`$ within halos with mass in the range corresponding to $`S`$ to $`S+dS`$. The halo abundance is then

$$\frac{dn}{dM}=\frac{\overline{\rho }_0}{M}\left|\frac{dS}{dM}\right|f(\delta _c(z),S),$$
(3)
where $`dn`$ is the comoving number density of halos with masses in the range $`M`$ to $`M+dM`$. In the model of Press & Schechter (1974),

$$f_{\mathrm{PS}}(\delta _c(z),S)=\frac{1}{\sqrt{2\pi }}\frac{\nu }{S}\mathrm{exp}\left[\frac{\nu ^2}{2}\right],$$
(4)
where $`\nu =\delta _c(z)/\sqrt{S}`$ is the number of standard deviations that the critical collapse overdensity represents on the mass scale $`M`$ corresponding to the variance $`S`$.

However, the Press-Schechter mass function fits numerical simulations only roughly, and in particular it substantially underestimates the abundance of the rare halos that host galaxies at high redshift. The halo mass function of Sheth & Tormen (1999, see also ) adds two free parameters that allow it to fit numerical simulations much more accurately (Jenkins et al., 2001). We note that these simulations followed very large volumes at low redshift, so that cosmic scatter did not compromise their accuracy. The matching mass function is given by

$$f_{\mathrm{ST}}(\delta _c(z),S)=A^{}\frac{\nu }{S}\sqrt{\frac{a^{}}{2\pi }}\left[1+\frac{1}{(a^{}\nu ^2)^q^{}}\right]\mathrm{exp}\left[\frac{a^{}\nu ^2}{2}\right],$$
(5)
with best-fit parameters (Sheth & Tormen, 2002) $`a^{}=0.75`$ and $`q^{}=0.3`$, and where normalization to unity is ensured by taking $`A^{}=0.322`$.

In order to calculate cosmic scatter we must determine the biased halo mass function in a given volume at a given mean density. Within the extended Press-Schechter model (Bond et al., 1991), the halo mass distribution in a region of comoving radius $`R`$ with a mean overdensity $`\overline{\delta }_R`$ is given by

$$f_{\mathrm{bias}\mathrm{PS}}(\delta _c(z),\overline{\delta }_R,R,S)=f_{\mathrm{PS}}(\delta _c(z)\overline{\delta }_R,SS(R)).$$
(6)
The corresponding collapse fraction in this case is given simply by eq. (1). Despite the relatively low accuracy of the Press-Schechter mass function, the relative change is predicted rather accurately by the extended Press-Schechter model. In other words, the prediction for the halo mass function in a given volume compared to the cosmic mean mass function provides a good fit to numerical simulations over a wide range of parameters (Mo & White, 1996; Casas-Miranda et al., 2002; Sheth & Tormen, 2002).

For our improved model we adopt a hybrid approach that combines various previous models with each applied where it has been found to closely match numerical simulations. We obtain the halo mass function within a restricted volume by starting with the Sheth-Tormen formula for the cosmic mean mass function, and then adjusting it with a relative correction based on the extended Press-Schechter model. In other words, we set

$`f_{\mathrm{bias}}(\delta _c(z),\overline{\delta }_R,R,S)=`$
$`f_{\mathrm{ST}}(\delta _c(z),S)\times \left[{\displaystyle \frac{f_{\mathrm{PS}}(\delta _c(z)\overline{\delta }_R,SS(R))}{f_{\mathrm{PS}}(\delta _c(z),S)}}\right].`$ (7)
As noted, this model is based on fits to simulations at low redshifts, but we can check it at high redshifts as well. Figure 4 shows the number of galactic halos at $`z1530`$ in two numerical simulations run by Yoshida et al. (2003), and our predictions given the cosmological input parameters assumed by each simulation. The close fit to the simulated data (with no additional free parameters) suggests that our hybrid model (solid lines) improves on the extended Press-Schechter model (dashed lines), and can be used to calculate accurately the cosmic scatter in the number of galaxies at both high and low redshifts. The simulated data significantly deviate from the expected cosmic mean \[eq. (5), shown by the dotted line\], due to the artificial suppression of large-scale modes outside the simulated box. We note that Yoshida et al. (2003) mentioned that the lack of large-scale modes might produce a systematically low halo abundance, particularly in the RSI model, but they did not quantify this effect.

As an additional example, we consider the highest-resolution first star simulation (Abel et al., 2002), which used $`l_{\mathrm{box}}=128`$ kpc and $`M_{\mathrm{min}}=7\times 10^5M_{}`$. The first star forms within the simulated volume when the first halo of mass $`M_{\mathrm{min}}`$ or larger collapses within the box. To compare with the simulation, we predict the redshift at which the probability of finding at least one halo within the box equals $`50\%`$, accounting for Poisson fluctuations. We find that if the simulation formed a population of halos corresponding to the correct cosmic average \[as given by eq. (5)\], then the first star should have formed already at $`z=24.0`$. The first star actually formed in the simulation box only at $`z=18.2`$ (Abel et al., 2002). Using eq. (7) we can account for the loss of large-scale modes beyond the periodic box, and predict a first star at $`z=17.8`$, a close match given the large Poisson fluctuations introduced by considering a single galaxy within the box.

The artificial bias in periodic simulation boxes can also be seen in the results of extensive numerical convergence tests carried out by Springel & Hernquist (2003). They presented a large array of numerical simulations of galaxy formation run in periodic boxes over a wide range of box size, mass resolution, and redshift. In particular, we can identify several pairs of simulations where the simulations in each pair have the same mass resolution but different box sizes; this allows us to separate the effect of large-scale numerical bias from the effect of having poorly-resolved individual halos. Specifically, their simulations Z1 and R4 \[see Table 1 in Springel & Hernquist (2003)\] used the same particle mass but R4 had a box length larger by a factor of 3.375 . The simulations Q1 and D4 are similarly related, as are Q2 and D5. In each case, the smaller simulation substantially underestimated the star formation rate at high redshift \[see Figure 10 in Springel & Hernquist (2003)\], Z1 by a factor of 5 at $`z=15`$ compared to R4, Q1 by a factor of 3 at $`z=8`$ compared to D4, and Q2 by a factor of 3 at $`z=10`$ compared to D5.

We note that there have been previous attempts to develop a model for the halo mass function in different environments, so that the model would be consistent with the Sheth-Tormen mass function of eq. (5) which accurately fits the cosmic mean mass function measured in numerical simulations. In order to identify the specific requirements for such a consistency, we first consider the analogous case of the extended Press-Schechter model and its relation to the Press-Schechter formula for the cosmic mean mass function. The extended Press-Schechter model is consistent with the mean mass function in the sense that $`f_{\mathrm{bias}\mathrm{PS}}`$ evaluated for an infinite box (i.e., in the limit where $`\overline{\delta }_R`$ and $`S(R)`$ both vanish) yields the Press-Schechter mass function:

$$f_{\mathrm{bias}\mathrm{PS}}(\delta _c(z),0,\mathrm{},S)=f_{\mathrm{PS}}(\delta _c(z),S).$$
(8)
This condition does not suffice, however, since there is an additional self-consistency test that any viable model must satisfy. Consider any fixed scale $`R`$. Suppose we consider a very large number $`N`$ of spheres of radius $`R`$ within the universe. The mean density $`\overline{\delta }_R`$ in each sphere is determined according to a probability distribution $`p(\overline{\delta }_R)`$; we assume that $`R`$ is large enough so that the probability of forming a halo out of all the mass on the scale $`R`$ is negligible, and so the distribution is a Gaussian with zero mean and variance $`S(R)`$ (see also §2.1). The number of galaxies in each sphere is given in the extended Press-Schechter model by eq. (6). As $`N\mathrm{}`$, the halo mass function averaged over all these spheres must approach the cosmic mean value, and it must also approach the ensemble-averaged mass function, where the averaging is performed over the probability distribution of $`\overline{\delta }_R`$. This yields the following self-consistency requirement, which is indeed satisfied by the extended Press-Schechter model:

$`{\displaystyle f_{\mathrm{bias}\mathrm{PS}}(\delta _c(z),\overline{\delta }_R,R,S)p(\overline{\delta }_R)𝑑\overline{\delta }_R}=`$
$`f_{\mathrm{bias}\mathrm{PS}}(\delta _c(z),0,\mathrm{},S).`$ (9)
Now we again consider attempts to construct an improved model that is consistent with the Sheth-Tormen mass function. Such a model must satisfy eq. (8) (except with $`f_{\mathrm{ST}}`$ on the right-hand side), and it must also satisfy eq. (9) in order to be self-consistent. The latter equation must be satisfied separately for every scale $`R`$ large enough to avoid collapsing \[i.e., that satisfies $`\sqrt{S(R)}\delta _c(z)`$\]. Previous proposed models (Sheth & Tormen, 2002; Gottlöber et al., 2003) satisfied simple consistency but not the self-consistency test. Our hybrid model of eq. (7) satisfies both requirements (with respect to the Sheth-Tormen mass function), a result that follows immediately from the fact that the extended Press-Schechter model also satisfies both requirements (with respect to the Press-Schechter mass function). Thus, our hybrid model is the first self-consistent model that is also consistent with the Sheth-Tormen mass function, at least when fluctuations are considered on large scales $`R`$ for which $`\overline{\delta }_R`$ is Gaussian distributed. As demonstrated in this section, our model also matches results from a wide array of numerical simulations (Abel et al., 2002; Yoshida et al., 2003; Mo & White, 1996; Casas-Miranda et al., 2002; Jenkins et al., 2001).

## 3. Implications

### 3.1. The nature of reionization

The photons of the cosmic microwave background have traveled to us mostly undisturbed after neutral atoms first formed in the universe at the cosmic recombination epoch. Radiation from the first generation of stars is thought to have reionized the hydrogen throughout the universe, transforming the IGM back into a hot and highly-ionized plasma.

The popular view developed in the literature (Arons & Wingert, 1972; Fukugita & Kawasaki, 1994; Shapiro et al., 1994; Haiman & Loeb, 1997; Gnedin, 2000; Barkana & Loeb, 2001) maintains that reionization ended with a fast, simultaneous, overlap stage throughout the universe. This view has been based on simple arguments and has been supported by numerical simulations with small box sizes. The underlying idea was that the ionized hydrogen (H II) regions of individual sources began to overlap when the typical size of each H II bubble became comparable to the distance between nearby sources. Since these two length scales were comparable at the critical moment, there is only a single timescale in the problem – given by the growth rate of each bubble – and it determines the transition time between the initial overlap of two or three nearby bubbles, to the final stage where dozens or hundreds of individual sources overlap and produce large ionized regions. Whenever two ionized bubbles were joined, each point inside their common boundary became exposed to ionizing photons from both sources, reducing the neutral hydrogen fraction and allowing ionizing photons to travel farther before being absorbed. Thus, the ionizing intensity inside H II regions rose rapidly, allowing those regions to expand into high-density gas that had previously recombined fast enough to remain neutral when the ionizing intensity had been low. Since each bubble coalescence accelerates the process, it has been thought that the overlap phase has the character of a phase transition and occurs rapidly. Indeed, the best simulations of reionization to date (Gnedin, 2000) found that the average mean free path of ionizing photons in the simulated volume rises by an order of magnitude over a redshift interval $`\mathrm{\Delta }z=0.05`$ at $`z=7`$.

Our results substantially modify this commonly accepted picture for the development of reionization. Overlap is still expected to occur rapidly, but only in localized high-density regions, where the ionizing intensity and the mean free path rise rapidly even while other distant regions are still mostly neutral. In other words, the size of the bubble of an individual source is about the same in different regions (since most halos have masses just above $`M_{\mathrm{min}}`$), but the typical distance between nearby sources varies widely across the universe. The strong clustering of ionizing sources on length scales as large as 30–100 Mpc introduces long timescales into the reionization phase transition. The sharpness of overlap is determined not by the growth rate of bubbles around individual sources, but by the ability of large groups of sources within overdense regions to deliver ionizing photons into large underdense regions. Simply put, the common view assumes that reionization occurred in patches a few Mpc in size, and that it ended nearly simultaneously in all of them. In reality, however, these two statements are contradictory. If the patches are a few Mpc in size, then there is a very large spread in their reionization redshifts. Conversely, if the spread is small, this implies (from Figure 3) that the patches must be far larger than is commonly assumed.

Note that the recombination rate is higher in overdense regions because of their higher gas density. These regions still reionize first, though, despite the need to overcome the higher recombination rate, since the number of ionizing sources in these regions is increased even more strongly as a result of the dramatic amplification of large-scale modes discussed earlier.

### 3.2. Limitations of current simulations

The shortcomings of current simulations do not amount simply to a shift of $`10\%`$ in redshift and the elimination of scatter, for several reasons. First, the effect that we have identified can be expressed in terms of a shift in redshift only within the context of the extended Press-Schechter model, and only if the total mass fraction in galaxies is considered and not its distribution as a function of galaxy mass. The halo mass distribution should still have the wrong shape, resulting from the fact that $`\mathrm{\Delta }z`$ in eq. 2 depends on $`M_{\mathrm{min}}`$. Furthermore, in our more accurate hybrid model (§ 2.2), the effect on the collapse fraction is no longer exactly equivalent to a shift in redshift. In any case, a self-contained numerical simulation cannot rely on approximate models and must directly evolve a very large volume<sup>1</sup><sup>1</sup>1Ciardi, Stoehr, & White (2003) simulated a 30 Mpc box, larger than previous simulations, but only gravity was directly simulated, and the mass resolution was three orders of magnitude lower than the minimum necessary to resolve most galactic halos at high redshift..

The second reason that current simulations are limited is that at high redshift, when galaxies are still rare, the abundance of galaxies grows rapidly towards lower redshift. Therefore, a $`10\%`$ relative error in redshift implies that at any given redshift around $`z10`$–20, the simulation predicts a halo mass function that can be off by an order of magnitude for halos that host galaxies (see Figures 1 and 4). This large underestimate suggests that the first generation of galaxies formed significantly earlier than indicated by recent simulations. This makes it easier to explain recent observations of the cosmic microwave background (Spergel et al., 2003) that suggest an early reionization at $`z15`$–20.

The third reason for the failure of simulations arises from the large cosmic scatter. This scatter can fundamentally change the character of any observable process or feedback mechanism that depends on a radiation background. Simulations in periodic boxes eliminate any large-scale scatter by assuming that the simulated volume is surrounded by identical periodic copies of itself. In the case of reionization, for instance, current simulations neglect the collective effects described above, whereby groups of sources in overdense regions may influence large surrounding underdense regions. In the case of the formation of the first stars due to molecular hydrogen cooling, the effect of the soft ultraviolet radiation from these stars, which tends to dissociate the molecular hydrogen around them (Haiman et al., 1997; Ricotti et al., 2002; Oh & Haiman, 2003), must be reassessed with cosmic scatter included.

### 3.3. Observational consequences

The spatial fluctuations that we have calculated fundamentally affect current and future observations that probe reionization or the galaxy population at high redshift. For example, there are a large number of programs searching for galaxies at the highest accessible redshifts (6.5 and beyond) using their strong Ly$`\alpha `$ emission (Hu et al., 2002; Rhoads et al., 2003; Maier et al., 2003; Kodaira et al., 2003). These programs have previously been justified as a search for the reionization redshift, since the intrinsic emission should be absorbed more strongly by the surrounding IGM if this medium is neutral. For any particular source, it will be hard to clearly recognize this enhanced absorption because of uncertainties regarding the properties of the source and its radiative and gravitational effects on its surroundings (Barkana & Loeb, 2003a, b; Santos, 2003). However, if the luminosity function of galaxies that emit Ly$`\alpha `$ can be observed, then faint sources, which do not significantly affect their environment, should be very strongly absorbed in the era before reionization. Reionization can then be detected statistically through the sudden jump in the number of faint sources (Haiman & Spaans, 1999; Haiman, 2002). Our results alter the expectation for such observations. Indeed, no sharp “reionization redshift” is expected. Instead, a Ly$`\alpha `$ luminosity function assembled from a large area of the sky will average over the cosmic scatter of $`\mathrm{\Delta }z1`$–2 between different regions, resulting in a smooth evolution of the luminosity function over this redshift range. In addition, such a survey may be biased to give a relatively high redshift, since only the most massive galaxies can be detected, and as we have shown, these galaxies will be concentrated in overdense regions that will also get reionized relatively early.

The distribution of ionized patches during reionization will likely be probed by future observations, including small-scale anisotropies of the cosmic microwave background photons that are rescattered by the ionized patches (Aghanim et al., 1996; Gruzinov & Hu, 1998; Santos et al., 2003), and observations of 21 cm emission by the spin-flip transition of the hydrogen in neutral regions (Tozzi et al., 2000; Carilli et al., 2002; Furlanetto et al., 2003). Previous analytical and numerical estimates of these signals have not included the collective effects discussed above, in which rare groups of massive galaxies may reionize large surrounding areas. These photon transfers will likely smooth out the signal even on scales significantly larger than the typical size of an H II bubble due to an individual galaxy. Therefore, even the characteristic angular scales that are expected to show correlations in such observations must be reassessed.

The cosmic scatter also affects observations in the present-day universe that depend on the history of reionization. For instance, photoionization heating suppresses the formation of dwarf galaxies after reionization, suggesting that the smallest galaxies seen today may have formed prior to reionization (Bullock et al., 2001; Somerville, 2002; Benson et al., 2002). Under the popular view that assumed a sharp end to reionization, it was expected that denser regions would have formed more galaxies by the time of reionization, possibly explaining the larger relative abundance of dwarf galaxies observed in galaxy clusters compared to lower-density regions such as our Local Group of galaxies (Tully et al., 2002; Benson et al., 2003a). Our results undercut the basic assumption of this argument and suggest a different explanation altogether. Reionization occurs roughly when the number of ionizing photons produced starts to exceed the number of hydrogen atoms in the surrounding IGM. If the processes of star formation and the production of ionizing photons are equally efficient within galaxies that lie in different regions, then reionization in each region will occur when the collapse fraction reaches the same critical value, even though this will occur at different times in different regions. Since the galaxies responsible for reionization have the same masses as present-day dwarf galaxies, this estimate argues for a roughly equal abundance of dwarf galaxies in all environments today. This simple picture is, however, modified by several additional effects. First, the recombination rate is higher in overdense regions at any given time, as discussed above. Furthermore, reionization in such regions is accomplished at an earlier time when the recombination rate was higher even at the mean cosmic density; therefore, more ionizing photons must be produced in order to compensate for the enhanced recombination rate. These two effects combine to make overdense regions reionize at a higher value of $`F_{\mathrm{col}}`$ than underdense regions. In addition, the overdense regions, which reionize first, subsequently send their extra ionizing photons into the surrounding underdense regions, causing the latter to reionize at an even lower $`F_{\mathrm{col}}`$. Thus, a higher abundance of dwarf galaxies today is indeed expected in the overdense regions.

The same basic effect may be even more critical for understanding the properties of large-scale voids, 10–30 Mpc regions in the present-day universe with an average mass density that is well below the cosmic mean. In order to predict their properties, the first step is to consider the abundance of dark matter halos within them. Numerical simulations show that voids contain a lower relative abundance of rare halos (Cen & Ostriker, 2000; Somerville et al., 2001; Mathis & White, 2002; Benson et al., 2003b), as expected from the raising of the collapse threshold for halos within a void. On the other hand, simulations show that voids actually place a larger fraction of their dark matter content in dwarf halos of mass below $`10^{10}M_{}`$ (Gottlöber et al., 2003). This can be understood within the extended Press-Schechter model. At the present time, a typical region in the universe fills halos of mass $`10^{12}M_{}`$ and higher with most of the dark matter, and very little is left over for isolated dwarf halos. Although a large number of dwarf halos may have formed at early times in such a region, the vast majority later merged with other halos, and by the present time they survive only as substructure inside much larger halos. In a void, on the other hand, large halos are rare even today, implying that most of the dwarf halos that formed early within a void can remain as isolated dwarf halos till the present. Thus, most isolated dwarf dark matter halos in the present universe should be found within large-scale voids (Barkana, 2003).

However, voids are observed to be rather deficient in dwarf galaxies as well as in larger galaxies on the scale of the Milky Way (e.g., Kirshner et al., 1981; Eder et al., 1989; Grogin & Geller, 1999, 2000; El-Ad & Piran, 2000; Peebles, 2001). A deficit of large galaxies is naturally expected, since the total mass density in the void is unusually low, and the fraction of this already low density that assembles in large halos is further reduced relative to higher-density regions. The absence of dwarf galaxies is harder to understand, given the higher relative abundance expected for their host dark matter halos. The standard model for galaxy formation may be consistent with the observations if some of the dwarf halos are dark and do not host stars. Large numbers of dark dwarf halos may be produced by the effect of reionization in suppressing the infall of gas into these halos. Indeed, exactly the same factors considered above, in the discussion of dwarf galaxies in clusters compared to those in small groups, apply also to voids. Thus, the voids should reionize last, but since they are most strongly affected by ionizing photons from their surroundings (which have a higher density than the voids themselves), the voids should reionize when the abundance of galaxies within them is relatively low. A quantitative analysis of how the reionization redshift varies with environment may help establish a common framework for explaining the observed properties of dwarf galaxies in environments ranging from clusters to voids.

## 4. Conclusions

We have shown that the important milestones of high-redshift galaxy formation, such as the formation of the first stars and the completion of reionization, occurred at significantly different times in different regions of the universe. This conclusion results from the fact that the temperature threshold, above which cooling and fragmentation of gas are possible, selects out dark matter halos that become exceptionally rare at high redshifts. Consequently, density fluctuations on large scales modulate the threshold for the collapse of high density peaks on small scales in the exponential tail of the Gaussian random field of density fluctuations, and introduce a remarkably large scatter in the abundance of star-forming galaxies at early cosmic times.

We have developed an improved method to calculate the cosmic scatter (see §2.2). This yields the first self-consistent analytic model that matches the halo mass function measured in various regions in numerical simulations that covered a wide range of the parameter space of region size, mean density, and redshift.

Since the characteristic distance between nearby sources of ionizing radiation varies widely across the universe, the overlap of the H II regions produced by these individual sources in the IGM occurs at significantly different times in different cosmic environments. Quantitatively, we find that the spread in the redshift of reionization should be at least an order of magnitude larger than previous expectations that argued for a sharp end to reionization (see § 3.1).

Current numerical simulations that treat gravity and hydrodynamics (Gnedin, 2000; Abel et al., 2002; Yoshida et al., 2003) largely eliminate this real cosmic scatter, and are artificially biased toward late galaxy formation since they exclude large-scale modes (see Figures 3 and 4). We find that galaxy formation within state-of-the-art simulations with $`324^3`$ particles is artificially biased to occur too late by a redshift interval $`\mathrm{\Delta }z0.5`$ at $`z=7`$ and $`\mathrm{\Delta }z2.5`$ at $`z=20`$. The box length used in state-of-the-art simulations of reionization (Gnedin, 2000; Yoshida et al., 2003) is 1.5–2 orders of magnitude below the minimum size necessary to treat the scatter reliably, and so alternative computational schemes (Barkana & Loeb, 2003c) must be implemented in order to quantify the implications of the large cosmic scatter on the reionization history. This scatter should affect the statistical fluctuations in the number and clustering properties of sources in surveys with a narrow field of view (such as the Hubble Deep Field), the luminosity function of Ly$`\alpha `$-emitting galaxies around the reionization redshift, the fluctuations in the 21 cm flux produced by the neutral IGM, the power spectrum of the secondary anisotropies in the cosmic microwave background, and the present abundance of dwarf galaxies in various environments (see § 3.3). Simulations limited to a small box may be able to study the scatter in the number density of galaxies by varying the mean density of the box, but such simulations cannot probe the global structure of reionization since this would involve the radiative transfer of ionizing photons over distances larger than the box size.

We thank Paul Steinhardt for suggesting to apply our work to galaxy formation in voids. We acknowledge support by NSF grant AST-0204514 and NATO grant PST.CLG.979414. R.B. is grateful for the kind hospitality of the Harvard-Smithsonian CfA and the Institute for Advanced Study, and the support of an Alon Fellowship at Tel Aviv University and of Israel Science Foundation grant 28/02/01. A.L. acknowledges sabbatical support from the John Simon Guggenheim Memorial Fellowship. This work was also supported in part by NSF grant AST-0071019 and NASA grant NAG 5-13292 (for A.L.).


text:
# HD 199143 and HD 358623: Two Recently Identified Members of the 𝛽 Pictoris Moving Group

## 1 Introduction

The discovery of an increasingly large number of young stars near Earth is providing fertile ground for studies of Galactic star formation, including the initial kinematics of the local stellar population, the evolution of the stars themselves, and the formation of their primordial disks. During the past two years, a number of authors have presented data on the rapidly rotating F8V star HD 199143 and the K7-M0 dwarf HD 358623. Based on their common proper motion, van den Ancker et al. (2000) (hereafter vdA00) concluded that these stars are associated and have an age of $``$ 20 Myr. To explain EUVE and ROSAT (Voges et al., 1999, 2000) detections, rapid rotation, variability in the continuum and emission lines of HD 199143, vdA00 postulated a protoplanetary disk or unseen stellar companion. Abundant lithium (EW = 400 mÅ; Mathioudakis et al. (1995)) in HD 358623 and chromospheric activity of both stars also indicate young ages. van den Ancker et al. (2001) (vdA01) drew on the similar proper motions of these two stars to classify them as a new young cluster, which they named the Capricornus association.

Adaptive optics (AO) imaging by Jayawardhana & Brandeker (2001) (hereafter, JB01) revealed companion candidates to HD 199143 and HD 358623 at separations of 1 and 2″ respectively. Astrometry by Neuhäuser et al. (2002) firmly established HD 358623 B as a true companion of HD 358623 A. vdA01 measured N (11 $`\mu `$m) and Q (19 $`\mu `$m) band excesses around both HD 199143 and HD 358623 implying circumstellar disks. Near-IR photometry by JB01 resulted in an extremely red color (J-K = 1.4) for HD 199143B which was interpreted as evidence of a circumsecondary disk. However, Chauvin et al. (2002) later showed that this photometry was highly dependent upon the deconvolution parameters and suggested that this influenced the assignment of an infrared excess to one of the stars.

Despite these discoveries, the reported ages of these stars have been inconsistent. Since age is the key factor in constraining models of lithium depletion, dust disk dissipation, and planetary formation, it is desirable to have an accurate estimate of the ages for these stars. In this paper, we present near- and mid-infrared photometry as well as optical spectroscopy to show that the space motions and spectral features of HD 199143 and HD 358623 are consistent with those of the $`\beta `$ Pictoris Moving Group; from this, and from high lithium abundances, we infer young ages for the two systems.

## 2 Observations

### 2.1 Near-IR Photometry with Adaptive Optics

We observed HD 199143 and HD 358623 at four epochs (Table 1). The AO system at the Shane 3-m telescope at Lick observatory was used with IRCAL, a 256 $`\times `$ 256 Rockwell HgCdTe array with a plate scale of 0.0756 arcsec/pixel (Lloyd et al., 2000). Strehl ratios (at 2.1$`\mu `$m) were in the range 0.3–0.45. Absolute photometry for HD 199143 and HD 358623 was obtained on 30 June 2002 by comparing open-loop images with the standard star HD 201941 (Elias et al., 1982). Measurements were done with IRAF using an aperture radius of 40 pixels (3″). Although the K’and Ks filters used for this study have slightly shorter bandpasses than the Elias et al. (1982) K-band filter, Persson et al. (1998) find that the color-correction between K and Ks is $`<`$ 0.02 mag even for stars as red as HD 199143 and HD 358623. Therefore, we have made no color correction here, but have added in quadrature an additional 0.02 mag error to the quoted K magnitude (Table 2).

At the Keck II telescope, we used the AO system with NIRSPEC’s slit-viewing camera, SCAM (McLean et al., 1998). The plate scale for SCAM’s 256 $`\times `$ 256 HgCdTe array is 0$`\stackrel{}{.}`$$`0167\pm `$0$`\stackrel{}{.}`$001/pixel. Strehl ratios for this run (K-band) were 0.10–0.46.

Photometric measurements of HD 199143 AB were performed using Starfinder routine (Diolaiti et al., 2000) and are reported in Table 2. Our estimates of the magnitudes of these stars agree well with those of other studies except in the case of the JB01 J-band measurement of HD 199143 B. However, Chauvin et al. (2002) note that the JB01 value may be brought into agreement with the others by the choice of a more robust myopic deconvolution algorithm.

On 14 June 2001, we measured the binary separations and position angles as 1$`\stackrel{}{.}`$02 $`\pm `$ 0$`\stackrel{}{.}`$03 and -36$`\stackrel{}{.}`$0 $`\pm `$ 0$`\stackrel{}{.}`$5 for HD 199143 and 2$`\stackrel{}{.}`$12 $`\pm `$ 0$`\stackrel{}{.}`$02 and 137$`\stackrel{}{.}`$9 $`\pm `$ 0$`\stackrel{}{.}`$5 for HD 358623. Our position angles differ slightly from those published by JB01, however due to the superior weather conditions under which our observations were made (Chauvin et al., 2002), we believe that our measurements are more accurate.

### 2.2 Optical Spectroscopy

We obtained high-resolution ($`0.15`$ Å) spectra of HD 199143 and HD 358623 (Fig. 1) with the Hamilton Echelle spectrometer (Vogt, 1987) at the Shane Telescope. HD 358623 was observed on UT 18 June 2001 and both HD stars were observed on 29 July 2001. Seeing was $``$ 1″ on both nights. Spectra were reduced using IRAF following a standard procedure (flat-field correction, etc.). Equivalent widths of H$`\alpha `$ and Li6708 features were measured using the IRAF task, splot. Terrestrial atmospheric features were removed using an A-type rapid rotator, HR 7235.

The heliocentric radial velocity of HD 358623 (Table 3) was measured by cross-correlating its spectrum with radial velocity standards HR 911, HR 6349, HR 8232, and HR 8969. All standards gave consistent results. The weighted mean of the measured radial velocities of HD 358623 is -6.8$`\pm `$1.3 km/s.

To verify that our H$`\alpha `$ spectrum of HD 358623 B was uncontaminated by light from HD 358623, we placed the slit first on B and then on a blank field at equal distance from and opposite position angle with respect to A. Figure 2 confirms that there is no scattered light from the primary star at this distance. The projected dimensions of the slit used for this observation are 1.5 arcsec wide in the dispersion direction and 2.5 arcsec long in the spatial direction.

### 2.3 Mid-infrared Photometry

Mid-infrared observations of HD 199143 and HD 368623 were taken on 7 Dec and 21 Dec 2001 at the W. M. Keck Observatory with the facility Long Wavelength Spectrograph (LWS) (Jones & Puetter, 1993). LWS’s 128 $`\times `$ 128 pixel array has a plate scale of 0.08 arcsec pixel-1 and a focal-plane field of view of 10$`\stackrel{}{.}`$2 square. The conditions were photometric during both sets of observations. The image quality was limited by seeing to FWHM of 0$`\stackrel{}{.}`$5 at 12 $`\mu `$m and 0$`\stackrel{}{.}`$6 at 18 $`\mu `$m on 7 Dec and 0$`\stackrel{}{.}`$8 at 8–12 $`\mu `$m and 0$`\stackrel{}{.}`$6 at 18 $`\mu `$m on 21 Dec.

The images were taken in a chop-nod mode with a chop amplitude of 10″, a nod of 10″ in the same direction, and a chopper frequency of 5 Hz. This scheme produces four images, two of which contain the object and two of which contain background only. The four images were double differenced for optimal background subtraction of the sky and telescope thermal emission. Filter characteristics are tied to the IRAS calibration.

The photometric standard star HR 8728 (Fomalhaut) was observed near in time and close in airmass to the observations of the HD stars. From its color-corrected IRAS flux densities, we calculate Fomalhaut’s mid-infrared magnitude as 0.88 mag from 8–18 $`\mu `$m with an uncertainty of 10%. On 21 Dec, observations of other standards at smaller airmasses are consistent with this magnitude and uncertainty for HR 8728 given typical Mauna Kea airmass corrections.

For HD 199143, the flux ratio of the binary was found by PSF fitting at 12 $`\mu `$m to be 5.6 $`\pm `$ 0.1 and the total flux density from both components was measured in a 4$`\stackrel{}{.}`$8 synthetic aperture. The 18 $`\mu `$m image had insufficient signal-to-noise ratio for the same PSF fitting procedure to work, so the flux densities were measured using small, diameter 1″, aperture photometry around each star. The total flux density in a very large aperture agreed with the sum in the small ones. For HD 358253, the photometry of each component was found in a synthetic aperture of diameter 1″ with residual background measured in an annulus 3$`\stackrel{}{.}`$5 from the source and subtracted. Statistical uncertainties calculated from the background were added in quadrature. The total uncertanties are reported in Table 4. Column two of this table represents the total integration time.

### 2.4 ROSAT and Near-infrared Imaging

Voges et al. (2000) report a bright (0.28 cts/s) X-ray source about one arc minute north of HD 199143 which itself is listed in the ROSAT All-Sky Survey. We examined ROSAT images of the field around HD 199143. With 30″ pixelation, we see no indication at all of an X-ray source one minute north of the star. We also investigated the area north of the star with UCLA’s twin-channel Gemini camera (McLean et al., 1994) at the Shane Telescope on 5 August 2001. A 60 s J-band exposure with limiting magnitude (3$`\sigma `$) of 19.1 mag reveals no source 1′ to the north of HD 199143. There is also no indication on the corresponding ROSAT image of an X-ray source at the listed position.

## 3 Results and Discussion

We use our photometry to place these stars in a color-magnitude diagram along isochrones by Baraffe et al. (2002) (Figure 3). All four stars lie on or to the right of the 20 Myr isochrone.

To investigate the reliability of Baraffe et al. isochrones at low masses, we selected 29 K and M-type dwarfs with less than 10% errors in parallax from the Gliese (1969) and Gliese & Jahreiss (1991) catalogues. When added to the color-magnitude diagram (Fig. 3), these field stars, which we expect to be of main sequence age, parallel the Baraffe 100 Myr isochrone. This is consistent with a pre-main-sequence phase for the HD 199143 and HD 358623 pairs.

Further evidence for the youth of HD 358623 comes from the 400 mÅ equivalent width of its Li $`\lambda `$6708 absorption line. Figure 1 of Song et al. (2002) – a plot of (B-V) color vs. Li $`\lambda `$6708 EW – shows that all stars with equivalent widths $``$400 mÅ belong to very young groups such as the TW Hydrae association, or the open cluster NGC 2264.

The ages of HD 199143 and HD 358623 are consistent with the 12$`{}_{4}{}^{}{}_{}{}^{+8}`$ age derived by Zuckerman et al. (2001) for members of the $`\beta `$ Pictoris Moving Group (BPMG) and confirmed by Ortega et al. (2002), who obtained a kinematic age of 11.5 Myr by tracing the orbits of all proposed BPMG members backwards in time through a 3-D Galactic potential. Furthermore, the Galactic space motion of HD 358623, relative to the Sun: (-9.8$`\pm `$0.6, -15.5$`\pm `$0.3, -11.1$`\pm `$0.5), derived using Hipparcos distances and proper motions, is consistent with the mean motions (-10.5, -15.8, -9.4) (Ortega et al., 2002) and velocity dispersions ($`\pm `$2km/sec in each direction) of the published BPMG stars (Zuckerman et al., 2001).

Figures 4 (HD 199143) and 5 (HD 358623) show synthetic spectra fit to the near-infrared photometry. In all four stars, the mid-infrared flux densities are consistent with those expected from the photospheres. These results contradict the measurements of vdA00 and vdA01 by showing no substantial infrared excess. The small excess that appears around HD 199143B at 18$`\mu `$m is only significant at the 2$`\sigma `$ level.

## 4 Conclusions

Data presented in this report show that HD 199143 AB and HD 358623 AB are members of the $`\beta `$ Pictoris Moving Group, rather than a separate Capricornus group. The high lithium abundance of HD 358623 suggests an age of less than roughly 10 Myr for this star. Near-IR adaptive optics photometry and models by Baraffe et al. (2002) also suggests a young age (less than 20 Myr). Our near- and mid-IR photometric results show no firm evidence for the presence of warm dust around either star. Investigation of a listed ROSAT All-Sky Survey source one arcminute north of HD 199143 revealed no near-IR counterpart to a limiting magnitude of J = 19.1 mag.

###### Acknowledgements.

We would like to thank Randy Campbell, Jay Farihi, and Patrick Lowrance for their useful input to this study, as well as the referees, Drs. A. Brandeker and R. Liseau. This research was supported in part by a NASA grant to UCLA and by NASA’s Astrobiology Institute. Some of the data presented herein were obtained at the W.M. Keck Observatory, which is operated as a scientific partnership among the California Institute of Technology, the University of California, and the National Aeronautics and Space administration. The Observatory was made possible by the generous financial support of the W.M. Keck foundation.


text:
# Damped Lyman-𝛼 Absorbers in Cosmological SPH Simulations: the “metallicity problem”

## 1. Introduction

DLAs are considered to be important reservoirs of neutral hydrogen in the universe at $`z3`$ (Storrie-Lombardi & Wolfe 2000). Studying the physical properties of DLAs, such as the distribution of their star formation rates (SFRs) and metallicity, will therefore provide us with important information on the history of star formation, galaxy formation, and chemical enrichment of the universe. As such, this information complements that provided by the emitted light from stars in high-redshift galaxies in a powerful and independent way.

Here we use cosmological SPH simulations of the $`\mathrm{\Lambda }`$CDM model with varying resolution and feedback strength to study the physical properties of DLAs. In our first study (Nagamine, Springel, & Hernquist 2003a), we showed that the $`\mathrm{\Lambda }`$CDM model is able to account for the observed abundance of DLAs at redshift $`z=34.5`$ quite well. Another important conclusion of our study was that earlier numerical work overestimated the DLA abundance significantly due to insufficient numerical resolution, lack of efficient feedback processes, and inaccuracies introduced in cooling processes when conventional formulations of SPH are used, giving rise to an overcooling problem.

In this conference proceedings, we focus our attention on the metallicity and the SFRs of DLAs. In particular, we discuss the “metallicity problem” that we currently face in both simulations and observations. We refer the readers to Nagamine, Springel, & Hernquist (2003b) for the details of our work.

## 2. Simulations

In this section, we briefly describe the SPH simulations that we use for our study. Our simulations include radiative cooling and heating with a uniform UV background, star formation, supernova (SN) feedback, as well as a phenomenological model for feedback by galactic winds. The latter allows us to examine, in particular, the effect of galactic outflows on the distribution of the SFR and metallicity of DLAs. For the details of these models, we refer readers to Springel & Hernquist (2003a) and a concise summary in Nagamine et al. (2003b).

We employ a “conservative entropy” formulation (Springel & Hernquist 2002) of SPH which alleviates numerical overcooling problems that affected earlier simulations. In addition, we utilize a series of simulations of varying boxsize and particle number to investigate the impact of numerical resolution on our results. The simulation parameters are summarized in Table 1. The adopted cosmological parameters of all runs are $`(\mathrm{\Omega }_m,\mathrm{\Omega }_\mathrm{\Lambda },\mathrm{\Omega }_b,\sigma _8,h)=(0.3,0.7,0.04,0.9,0.7)`$.

## 3. DLA metallicity in the simulations

In Figure 1, we show the projected gas metallicity vs. neutral hydrogen column density $`N_{\mathrm{HI}}`$ for random sight-lines in the ‘Q5’-run at $`z=3`$. We show the results from the ‘Q5’-run because it has the highest resolution at $`z=3`$. Each point in the figure represents one sight-line, and the contours are equally spaced on a logarithmic scale. The solid square symbols give the median value in each $`\mathrm{log}N_{\mathrm{HI}}`$ bin, with error bars indicating the quartiles on both sides. The median metallicity increases as $`N_{\mathrm{HI}}`$ increases, reaching solar metallicity at $`\mathrm{log}N_{\mathrm{HI}}23`$. This trend is expected because star formation is more vigorous in high $`N_{\mathrm{HI}}`$ systems.

An important result is that the median metallicity we find for DLAs in our simulations is much higher than the observed value. Note that observers typically find DLAs with $`20<\mathrm{log}N_{\mathrm{HI}}<22`$ to have a metallicity of $`\mathrm{log}(Z/Z_{})1.5`$ (e.g. Boissé et al. 1998; Pettini et al. 1999; Prochaska & Wolfe 2000), as indicated by the shaded region.

Another point to notice is the existence of systems that have both high $`N_{\mathrm{HI}}`$ and high metallicity; these systems tend to be absent in observations. Dust obscuration is sometimes invoked to reconcile this result with observations, however, recent observational tests suggest that the dust extinction effect is not so strong (Ellison et al. 2001; Prochaska & Wolfe 2002), and the solution could rather lie in a more adequate treatment of star formation and supernova feedback. For example, Schaye (2001) argues that the conversion of neutral hydrogen atoms into a molecular form, which we have not yet implemented in our simulations, would introduce a physical limit to the highest $`N_{\mathrm{HI}}`$ that can be attained. This process would eliminate the highest $`N_{\mathrm{HI}}`$ systems, but would not reduce the metallicity of low $`N_{\mathrm{HI}}`$ systems because it would make the star formation even more efficient than in our current simulations. Other possibilities include the existence of metal-loaded winds, which we will discuss in more detail in Section 6. It hence remains to be seen whether future cosmological simulations with a more sophisticated modeling of star formation and SN feedback processes confirm the existence of high-metallicity, high-$`N_{\mathrm{HI}}`$ systems, which could then turn into an interesting challenge for the CDM model.

## 4. Star formation rates in DLAs

In Figure 2, we show the distribution of sight-lines in the ‘Q5’-run at $`z=3`$ on the plane of metallicity vs. projected star formation rate surface density $`\mathrm{\Sigma }_{\mathrm{SFR}}`$ (in proper units of $`M_{}\mathrm{yr}^1\mathrm{kpc}^2`$). The shaded area roughly indicates the region of recent observational estimate by Wolfe et al. (2003b).

The shape of the simulated distribution is easy to understand: Because $`\mathrm{\Sigma }_{\mathrm{SFR}}`$ is tightly correlated with $`N_{\mathrm{HI}}`$ in our simulation (as shown in Figure 3; the Kennicutt  law), the distribution seen in Figure 2 is simply a reflection of Figure 1 around a diagonal line. The star formation model adopted in our simulation depends only on local physical quantities; e.g. the local gas density. The free parameter of the model was chosen to reproduce the Kennicutt law in isolated disk galaxies at low redshift, but was kept fixed as a function of time. We hence implicitly assumed that the Kennicutt law holds at all redshifts, and our simulation results reflect this assumption.

The observational estimates of SFR by Wolfe et al. (2003b) agree well with the simulation, although they do not follow the Kennicutt law tightly. Using the DLA abundance information, Wolfe et al. have estimated the volume-averaged SFR density of DLAs and found that it is comparable to that of the Lyman-break galaxies’ (LBGs). This leads to the so-called global “missing metals”-problem (e.g. Pagel 2002; Pettini 2003), where the total amount of metals seen in DLAs are not sufficient to account for all the metals expected from the observed SFR density at $`z=3`$. However, even before taking the volume-average, both our simulations and the observations by Wolfe et al. (2003a,b) seem to be pointing to the same problem; i.e. the observed DLA metallicities are much lower than those expected from the (either observed or simulated) DLA star formation rates.

## 5. Discussion

While the projected SFR density at $`z=3`$ in our simulations is plausible and agrees with current observations well, the median metallicity of DLAs appears to be too high compared to the values typically observed for DLAs. There are a number of possible explanations for this problem, and we will briefly discuss some of the most prominent possibilities.

One potential reason for high metallicity in DLAs is that the feedback by galactic winds is not efficient enough in blowing out metals from DLAs. Clearly, if the feedback by winds were stronger, then star formation and hence metal creation in DLAs would be more strongly suppressed. However, simply making the winds stronger and blowing out more gas will not necessarily decrease the metallicity of DLAs much, because in our current simulation model, the winds transport away metals and gas at the same time; i.e. the wind’s initial metallicity is assumed to be equal to that of the gas of the DLA, leaving the ratio of metal and gas mass in the DLA unchanged.

It is however quite plausible that the wind is metal-loaded compared to the gas in the DLA, as is for example suggested by simulations of SN explosions (e.g. MacLow & Ferrara 1999, Bromm et al. 2003). After all, the ejecta of SN are heavily enriched and inject large parts of the energy that is assumed to ultimately drive the outflow. If the mixing with other DLA-gas is not extremely efficient before the outflow occurs, it can then be expected that the wind material has potentially much higher metallicity than the DLA, thereby selectively removing metals.

A related possibility concerns the metallicities of the cold and diffuse phases of the DLA. In the present study, we assumed that metals are always efficiently and rapidly mixed between the gas of the cold clouds and the ambient medium, such that there is a homogeneous metal distribution in the DLA (operationally, we used only a single metallicity variable for each gas particle, reflecting this assumption). However, this assumption may not be fully correct. If the metals were preferentially kept in the hot phase of the ISM after they are released by SNe, then they would not be observed in the cold gas that is responsible for the DLAs. Since we did not track the metal distribution in cold and hot phases separately in the current simulations, we may then have overestimated the amount of metals in DLAs by counting those in the hot phase as well as those in the cold phase. Note that a more detailed tracking of metals in the simulation, separately for hot and cold phases of the gas, could in principle be done easily on a technical level. The difficulty however lies in obtaining a reasonable description of the physics that governs the exchange of metals between the different phases of the ISM, something that is presently not attainable from either observation or theory.

The viability of the feedback model in the simulations can also be tested by comparing with observations of the Lyman-$`\alpha `$ forest, which is generated by systems of much lower column density than the DLAs studied in this work. Curiously, an analysis using the current simulation series (Springel et al. 2003, in preparation), as well as a study by Theuns et al. (2002), suggest that the spectral features of the Ly-$`\alpha `$ forest are not significantly affected by the feedback from galactic outflows, despite the fact that the wind strength is taken to be on the ‘strong side’ in these studies, and despite the fact that a non-negligible fraction of the IGM volume is heated by the winds. Note that in our simulations with a strong wind model, the $`\mathrm{H}_\mathrm{I}`$ mass density in the entire simulation box is somewhat lower than suggested by observational estimates (Nagamine et al. 2003a), therefore it appears problematic to increase the wind strength even beyond the present value.

The high metallicity of DLAs may also be related to the steep luminosity function of galaxies in our SPH simulations (Nagamine et al. 2003c). An analysis using a population synthesis model shows that the luminosity function still has a very steep slope at the faint end even at low redshift, similar to that of the dark matter halo mass function. (It is not obvious whether the steep faint end in the simulation at $`z>1`$ is a problem, because it is not well constrained observationally at $`z>1`$ yet.) This means that the formation of low-mass galaxies in our simulation was not suppressed enough, or equivalently, that star formation was too efficient in low-mass haloes. Stronger winds may help to alleviate this problem, but it appears unlikely that our present feedback model can solve it satisfactorily simply by adopting a higher efficiency parameter for feedback. It is more plausible that additional physical processes need to be considered in a more faithful way. One simple possibility for this is related to the UV background field, which is turned on by hand at $`z=6`$ in our present simulations to mimic reionization of the Universe at a time when the first Gunn-Peterson troughs in spectra to distant quasars are observed (Becker et al. 2001). However, it is possible (in fact suggested by the WMAP satellite) that the Universe was reionized at much higher redshift. The associated photoheating may have then much more efficiently impaired the formation of low-mass galaxies than in our present simulations. We plan to explore this possibility in future work by adopting different treatments of the UV background radiation field.

In conclusion, we have shown that the DLAs found in our simulation series have many plausible properties. In particular, they are in good agreement with recent observations of the total neutral hydrogen mass density, the $`N_{\mathrm{HI}}`$ distribution function, the abundance of DLAs, and the distribution of SFR in DLAs. However, our simulated DLAs show typically considerably higher metallicity than what is presently observed for the bulk of these systems. This likely indicates that metal transport and mixing processes have not been efficient enough in our simulations. It will be interesting to study more sophisticated metal enrichment models in future simulations in order to further improve our understanding of the nature of DLAs in hierarchical CDM models.

## References

Becker, R. H., et al. 2001, AJ, 122, 2850

Boissé P., Le Brun V., Bergeron J., Deharveng J. M. 1998, A&A, 333, 841

Bromm, V., Yoshida, N., Hernquist, L. 2003, ApJ, 596, 135

Ellison S. L., Yan L., Hook I. M., Pettini M., Wall J. V., Shaver P. 2001, A&A, 379, 393

Kennicutt, R. C. Jr., 1998, ARA&A, 36, 189

Mac Low, M. M. & Ferrara, A. 1999, ApJ, 513, 142

Nagamine, K., Springel, V., Hernquist, L. 2003a, MNRAS, submitted (astro-ph/0302187)

Nagamine, K., Springel, V., Hernquist, L. 2003b, MNRAS, accepted (astro-ph/0305409)

Nagamine, K., Springel, V., Hernquist, L., & Machacek, M. 2003c, in preparation.

Pagel 2002, in ASP Conf. Series 253, ‘Chemical Enrichment of Intracluster and Intergalactic Medium’, eds. R. Fusco-Femiano & F. Matteucci, San Francisco: ASP, 489

Pettini 2003, Lectures given at the XIII Canary Islands Winter School of Astrophysics, ‘Cosmochemistry: The Melting Pot of Elements’ (astro-ph/0303272)

Pettini M., Ellison S. L., Steidel C. C., Bowen D. V. 1999, ApJ, 510, 576

Prochaska, J. X. & Wolfe, A. M., 2002, ApJ, 566, 68

Prochaska, J. X. & Wolfe, A. M. 2000, ApJ, 533, L5

Schaye, J. 2001, ApJ, 562, L95

Springel, V., Hernquist, L. 2002, MNRAS, 333, 649

Springel, V., Hernquist, L. 2003a, MNRAS, 339, 289

Springel, V., Hernquist, L. 2003b, MNRAS, 339, 312

Springel, et al., 2003, in preparation

Storrie-Lombardi L. J., Wolfe A. M. 2000, ApJ, 543, 552

Theuns, T., Viel, M., Kay, S., Schaye, J., Carswell, R. R., & Tzanavaris, P. 2002, ApJ, 578, L5

Wolfe A. M., Prochaska, J. X., Gawiser, E. 2003a, ApJ, 593, 215

Wolfe A. M., Gawiser, E., Prochaska, J. X. 2003b, ApJ, 593, 235