text:
# Simulations of Pregalactic Structure Formation with Radiative Feedback

## 1. Introduction

A critical epoch in the evolution of the universe is the the post-recombination period during which the first stars form, affect their environment and help determine how the intergalactic medium (IGM) ultimately becomes reionized. Recent observations are beginning to constrain models of this epoch. The observation of quasars at redshifts as high as $`z=5.8`$ (Fan et al. 2000) with no accompanying Gunn-Peterson absorption trough in their spectra indicates that the onset of reionization occurs quite early ($`z7`$) and so the first generation of stars forms even earlier. Observations of quasar spectra (Ellison et al. 1999, Ellison et al. 2000) suggest that metallicities possibly as high as $`\mathrm{log}CIV/HI=2.6`$ may be present even in the low column density $`\mathrm{log}N_{\mathrm{HI}}<14.5`$ Lyman alpha forest. This too may be indicative of an active period of pre enrichment of the IGM by an early generation of stars (Ostriker & Gnedin 1996). Future multiwavelength observations may be sensitive probes of this star formation epoch and the subsequent onset of reionization. A partial list of proposed observations include direct imaging of quasars and star clusters at $`z>10`$ by the Next Generation Space Telescope (Haiman & Loeb 1999, Barkana & Loeb 2000), absorption measurements in the 21 cm line using the Square Kilometer Array to identify the first epoch of massive star formation (Tozzi et al. 1999), and searches for secondary anisotropies in the Cosmic Microwave Background by the next generation of millimeter telescopes such as ALMA to constrain properties of the reionization epoch (Benson et al. 2000). Thus it is important to understand what our best cosmological models predict for the evolution of these early star forming structures and what governs the amount of gas that can ultimately condense into stars.

It has long been understood that a high redshift population of stars could form naturally in hierarchical models of structure formation. In these models small masses collapse first and later merge into larger structures. The ability of these objects to form stars depends on the ability of gas within these objects to cool radiatively, reduce pressure support, and condense. However, the first structures to become self-gravitating in these models would have temperatures $`T<10^4\mathrm{K}`$, too cool for HI line cooling to be effective. The only other coolant available to gas with a primordial composition is molecular hydrogen. Therefore it is the complex chemistry of H2 in these objects that controls whether the gas can condense to sufficiently high densities for star formation to take place.

Early analytical models of the collapse of the first small scale structures studied the effect of H2 cooling within highly idealized cloud geometries, typically a homogeneous gas cloud radiatively cooling while undergoing a uniform free fall collapse (Saslow & Zipoy 1967, Peebles & Dicke 1968, Hirasawa 1969, Matsudo et al. 1969, Hutchins 1976, Silk 1983, Palla et al. 1983, Lepp & Shull 1984, Lahov 1986). 1-D simulations sought to improve upon these analytical estimates, but still had to impose the shape of the initial density profiles and mass of the virialized object as initial conditions within an assumed cloud geometry (Haiman, Thoul & Loeb 1996, Tegmark et al. 1997, Omukai & Nishi 1998, Nishi et al. 1998, Nakamura & Umemura 1999). Haiman, Thoul & Loeb (1996) followed the dynamics of gas and cold dark matter and the H2 chemistry in a single, spherical density peak and showed the importance of shell crossing between the gas and dark matter in determining the fraction of mass within an object that could collapse and the minimum mass for collapse. In the same spirit Tegmark et al. (1997) used top hat density profiles and a simplified collapse criterion, i.e. that the cooling time be less than the Hubble time for the object right after virialization, to obtain rough estimates of the critical H2 mass fraction needed for cooling and the minimum mass (as a function of virialization redshift) that could produce it. However, the minimum mass to collapse determined in this way is strongly dependent on the H2 cooling function used. Abel et al. (1998) and Fuller & Couchman (2000) recalculated this minimum mass for top hat density profiles using the Tegmark et al. (1997) collapse criterion with updated H2 cooling functions (Lepp & Shull 1983 and Galli & Palla 1998, respectively) and found that the minimum mass able to collapse (at virialization redshifts $`30`$) was reduced by as much as an order of magnitude.

Molecular hydrogen is fragile and easily photodissociated by photons whose energies lie below $`13.6\mathrm{eV}`$ (Stecker & Williams 1967). Thus even before reionization, radiation from the first stars could affect H2 abundances at large distances (Dekel & Rees 1987). Haiman, Rees, & Loeb (1997) quantified this negative radiative feedback in more detail. They studied the collapse of homogeneous, spherical gas clouds in the presence of a soft UV photodissociating flux assuming both ionization and chemical equilibrium within the cloud to determine which masses could still collapse. Haiman, Abel & Rees (1999) improved on that calculation by solving spherically symmetric radiative transfer coupled with the time dependent rate equations for the H2 chemistry across truncated isothermal spheres. They also modified the collapse criterion such that a structure was said to be able to collapse and form stars if at any time during the object’s history the H2 cooling time in the core became less than the time elapsed since the core’s formation. This allowed the fraction of H2 abundance required for cooling and collapse of the core to be lower in certain cases than previous estimates of the “critical” H2 fraction based on the Tegmark et al. (1997) collapse criterion. Analytical models that include a first population of OB-like stars have been used to probe the suppression of star formation elsewhere within the host cloud due to the destruction of H2 coolant by the direct radiation from the first stars (Omukai & Nishi 1999; Glover & Brand 2000).

The major drawback affecting all of the above studies whether with or without radiative feedback is their one dimensional nature. Hierarchical structure formation that predicts the growth of structure in and along filaments is inherently three dimensional. However, fully three dimensional simulations are computationally challenging due to the need for large dynamic range to model simultaneously the internal structure of the first star forming structures and the effects of gravitational tidal forces and merging on their growth. Such 3-D simulations have only recently become computationally feasible. Gnedin & Ostriker (1997) simulated the thermal history of the intergalactic medium and the population of first objects to collapse but had insufficient resolution to probe inside these first structures. Abel et al. (1998) used a static nested grid approach and constrained realizations of $`3`$$`4\sigma `$ primordial fluctuations in a cluster normalized SCDM model to study when these structures were able to cool, but were also limited by spatial resolution ($`1`$ kpc) from probing deeply within the peaks. Subsequent higher resolution 3-D simulations focused on the properties of a single collapsing cloud. Bromm, Coppi, & Larson (1999) studied the fragmentation of a $`2\times 10^6M_{}`$ rotating disk under highly idealized initial conditions using a TREESPH code. Fuller & Couchman (2000) used an adaptive smoothing-kernal SPH coupled to a P3M gravity solver in a simulation box of comoving length $`25h^1`$ kpc to investigate the properties of the first peak to collapse. However, because of the small box size they needed to impose an overdense cold dark matter cosmology to model the formation of the cloud within its host filament. The highest spatial and mass resolution self-consistent cosmological 3-D simulations to date of the collapse and fragmentation of primordial clouds (Abel, Bryan, & Norman 1998, 2000; Norman, Abel, & Bryan 2000) use Eulerian adaptive mesh refinement (AMR) techniques to achieve extremely large dynamic range (currently $`3\times 10^7`$). They start from cosmological initial conditions in a cluster normalized SCDM model to study the detailed internal evolution of a single, collapsing cloud in interaction with its environment in a box with comoving length $`L=128`$ kpc.

In this paper we present the results of fully 3-D Eulerian AMR simulations of the formation and collapse of primordial pregalactic structure in a realistic $`\mathrm{\Lambda }`$CDM cosmology. This study expands upon previous work in several important ways. First, we use a larger simulation volume ($`1`$ Mpc3 comoving) to model more accurately the gravitational tidal forces and effects of merging on the evolving structures. Second, we introduce a simple model for radiative feedback. We analyze the collapse properties of individual peaks with and without the presence of this soft UV background and compare our results with previous 1-D studies. Finally, we use the population of pregalactic objects evolving together in the most active structure forming region of our simulation to develop statistics on the fraction of gas that is able to cool due to the presence of H2 and the fraction of gas that is both cold and dense enough to be available for star formation as functions of cloud mass and the average soft UV radiation field it experiences. These statistics represent the key new result of this paper. The amount of gas available for star formation in pregalactic objects is a critical input parameter into recent semi-analytical attempts to model stellar feedback from the first stars on the intergalactic medium, subsequent structure formation, and the nature of reionization (Gnedin & Ostriker 1997, Ciardi et al. 2000, Ferrara et al. 2000, Ciardi, Ferrara, & Abel 2000).

While we focus on the negative effect of radiation in this paper, it should be noted that positive feedback is possible if the first objects also emit high energy X-ray photons (Haiman, Abel & Rees 1999; Oh 2000). We will address this point in future work.

This paper is organized in the following way. In Section §2 we describe the set-up of our simulations including simulation parameters and technique, cosmological model, H2 chemistry, and radiative feedback. In Section §3 we present the characteristics of the simulated data set of pregalactic objects identified in our simulations and used for further analysis. In Section §4 we use the data set as a whole to develop fitting formulaes for the fraction of gas that can cool, become dense and thus be available for star formation with and without the presence of a background H2 photodissociating flux. In Section §5 we use radial profiles to show the internal collapse characteristics of these structures. We consider both the time evolution of a single peak for a fixed background H2 photodissociating flux and also how the internal properties of a given cloud change at fixed redshift when the level of H2 photodissociating flux is varied. Armed with these examples of the internal evolution and structure of the collapsing clouds, we discuss in Section §6 the neglect of self-shielding and H- photodetachment in our simulations and then use a simple analytical argument to illuminate the essential physics of the collapse. We summarize our results in Section §7.

## 2. The Simulations

A self-consistent three dimensional cosmological simulation of the collapse of low mass ($`10^5`$$`10^7M_{}`$) protogalaxies at high redshift ($`30>z>19`$) requires both large enough simulation volume to model the gravitational tidal forces at work during the formation epoch and high spatial (and mass) resolution to probe within the small collapsing clouds. One technique that has recently been able to achieve the large dynamic range necessary to shed light on the collapse and fragmentation of these first star producing structures is adaptive mesh refinement (AMR). The three dimensional AMR algorithm used in our simulations is described more fully elsewhere (Bryan 1999; Bryan & Norman 1997, 1999; Norman & Bryan 1998). Briefly, the code uses an adaptive hierarchy of rectangular grid patches at various levels of resolution very similar to the algorithm described by Berger & Collela (1989). Each grid patch covers some region within its parent grid needing additional refinement and may itself become a parent grid to an even higher resolution child grid. We take the ratio of parent to child grid mesh spacing to be two. The dark matter is followed using methods similar to those presented by Couchman (1991) and the gas hydrodynamics uses the artificial viscosity method of Stone & Norman (1992).

We use a flat, low matter density $`\mathrm{\Lambda }`$CDM model for our simulations with $`\mathrm{\Omega }_0=0.4`$, $`\mathrm{\Omega }_b=0.05`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=0.6`$, $`h=0.65`$, $`\sigma _8=0.8`$, and $`n=1`$ where $`\mathrm{\Omega }_0`$, $`\mathrm{\Omega }_b`$, and $`\mathrm{\Omega }_\mathrm{\Lambda }`$ are the fraction of the critical energy density carried in nonrelativistic matter, baryons, and vacuum energy, respectively, $`h`$ is the dimensionless Hubble parameter in units of $`100`$ km s-1 Mpc-1 , $`\sigma _8`$ is the density fluctuation normalization in a sphere of radius $`8h^1`$ Mpc, and $`n`$ is the slope of the primordial density perturbation power spectrum. The parameters are chosen to provide good consistency with observation. The fluctuation normalization is consistent with the CMB quadrupole as measured by COBE (Bunn & White 1997) and also with observations of the number density of galaxy clusters (White, Efstathiou & Frenk 1993; Bond & Myers 1996). $`\mathrm{\Omega }_bh^2`$ is consistent with big bang nucleosynthesis and the measured abundance of primordial deuterium (Copi, Schramm, & Turner 1995; Burles & Tytler 1998), and $`\mathrm{\Omega }_\mathrm{\Lambda }`$ is consistent with the upper limit ($`\mathrm{\Omega }_\mathrm{\Lambda }<0.7`$) of Maoz & Rix (1993) and the best fit parameters of Ostriker & Steinhardt (1995).

The simulations are initialized at $`z=99`$ with density perturbations generated for the above $`\mathrm{\Lambda }`$CDM model using the Eisenstein & Hu (1998) transfer functions in a comoving simulation volume of $`1`$ Mpc3. This simulation volume is large enough to ensure that the fundamental mode in the box remains linear at the lowest redshifts considered here. A low resolution run is used first to identify a region of active structure formation. In order to achieve high mass resolution in the region of interest the simulation is reinitialized with multiple static refinement levels surrounding the chosen region such that each successive static parent grid contains the lagrangian volume of its child grid. This results in a mass resolution in the initial conditions for the region of interest of $`4.78(38.25)M_{}`$ for the gas (dark matter), respectively. In addition, that region is allowed to dynamically refine so that the local Jeans length is resolved at all times by at least four grid zones and no grid cell contains more than four times the initial gas mass element ($`4.78M_{}`$) or ten times the initial dark matter element ($`382.5M_{}`$). We limit the total refinement to $`14`$ levels within a $`64^3`$ top grid which results in a maximum dynamic range of $`1,048,576`$. The comoving spatial resolution of $`0.95`$ pc at maximum refinement translates into a physical spatial resolution of $`0.03`$ pc ($`0.05`$ pc) at $`z=30`$ ($`z=19`$), respectively. We call a peak “collapsed” when it reaches maximal refinement. Once this occurs, we still wish to follow the subsequent evolution of structure in the region, but are uninterested in the details of the collapsed peak. In order to prevent the peak from collapsing further (which could cause numerical instabilities), we introduce a form of artificial pressure support. This is done by defining, for each cell, an effective pressure which is the greater of the thermal pressure and $`KG\rho _b^2\mathrm{\Delta }x_f^2/\mu `$, where $`\rho _b`$ is the baryon density in the cell, $`\mathrm{\Delta }x_f`$ is the cell width on the finest level and $`\mu =1.22m_H`$ is the usual mean mass per particle. The dimensionless constant $`K`$ is set to 100, a value which spreads the mass over a spherical region with a radius of several cells. The form of this expression is chosen by matching the thermal energy and gravitational self-energy of the gas; notice also that it results in a polytropic equation of state. We also smooth the gravitational potential by $`4\mathrm{\Delta }x_f`$ which helps to reduce the amplitude of the artificial pressure support required. Although introducing artificial pressure support causes us to lose information about the detailed morphology of the inner $`1`$$`2`$ pc ($`1`$ % of the virial radius) of the peak after collapse, it should not affect the determination of cooled gas fractions in Section §4 since the cooling region lies outside of the collapsed H2 core.

Our data set consists of results from four simulations starting from identical cosmological initial conditions on the density fields. Three of the simulations follow the nonequilibrium, time dependent evolution of nine chemical species (H, H+, He, He+, He++, e-, H2 , H$`{}_{2}{}^{}{}_{}{}^{+}`$, H-) using the algorithm of Anninos, et al. (1997) initialized with post-recombination abundances (Anninos & Norman 1996). The reaction network for even this simple chemistry is complex. Numerous authors have worked to identify the dominant reactions and calculate their rate coefficients including recently Haiman, Thoul & Loeb (1996), Abel et al. (1997) and Galli & Palla (1998). We use primarily the rate coefficients of Abel et al. (1997) with the exception of rate coefficients for the processes $`H+e^{}H^{}+\gamma `$, $`H_2+e^{}2H+e^{}`$, and $`H^{}+H^+H_2^++e^{}`$ which are taken from Shapiro & Kang (1987). We have also included the three-body process $`2H+HH_2+H`$ (Orel 1987, Palli et al. 1983), although its contribution is negligible at the densities considered in this work. The small post-recombination abundance of H2 in the intergalactic medium ($`2\times 10^6`$) is too low for cooling to be effective. However, at the densities and temperatures in the interiors of collapsing proto-galactic clouds, H2 formation is catalyzed by H- ions through the process

$$H+e^{}H^{}+\gamma $$
(1)
$$H+H^{}H_2+e^{}.$$
(2)
and can reach abundances high enough to act as an efficient coolant allowing the first stars to form.

Once the first stars form, radiation above $`13.6`$ eV produced by these stars is quickly absorbed by the dense neutral hydrogen surrounding them. However, the intergalactic medium is optically thin to radiation below the HI ionization limit. This radiation escapes to form a soft ultraviolet background of radiation throughout neighboring regions with sufficient energy to photodissociate H2 by means of the Solomon process

$$H_2+\gamma H_2^{}H+H,$$
(3)
thus inhibiting continued star formation (Haiman, Abel & Rees 1999). $`H_2^{}`$ is any of $`76`$ Lyman-Werner resonances in the energy range $`11.18`$$`13.6\mathrm{eV}`$. We use a very simple model to investigate the negative feedback effects of this soft UV radiation. The radiation field is assumed to be constant throughout the simulation with mean photon energy of $`12.87\mathrm{eV}`$. In terms of the constant mean flux $`F_{LW}`$, the H2 photodissociation rate coefficient for the Solomon process denoted in Equation 3 is (Abel, et al. 1997)

$$k_{diss}=1.1\times 10^8F_{LW}\text{s}\text{-1}.$$
(4)
We distinguish three of our simulations by different levels of this background flux, i.e. $`F_{LW}=0`$ (first structure formation), $`10^{22}`$, and $`10^{21}`$ erg s-1 cm-2 Hz-1 , respectively. These flux levels correspond to $`J_{21}`$ normalizations of the spectral intensity (Haiman, Abel, & Rees 1999) of $`J_{21}=0`$, $`10^{2.1}`$ and $`10^{1.1}`$ erg s-1 cm-2 Hz-1 sr-1 , respectively. We also include results from a fourth simulation that starts from the same cosmological initial conditions, but evolves only six species (H, H+, He, He+, He++, e-) representing the extreme case when destruction of molecular hydrogen is complete.

We neglect the time dependent turn on and build-up of the soft radiative field. Thus the UV radiation experienced by the collapsing structures may be interpreted as due either to “direct”, nearby sources or a more diffuse soft background. We also neglect the effect of self-shielding by the collapsing cores. As will be discussed more fully in Section §6, this should have only a small affect on our results. In any case each of these neglected effects would tend to reduce the impact of the photodissociating flux on the collapsing clouds. Thus our results for the simulations with $`F_{LW}=10^{22}`$ and $`10^{21}`$ erg s-1 cm-2 Hz-1 should be interpreted as representing the maximal suppression of structure formation at these flux levels.

We also neglect the photo-detachment of H- by the UV radiation field. One might argue that the presence of the radiative background might photo-detach H- before it could catalyze H2 formation by the process given in Equation 2. This could hinder the collapse of clouds by preventing the build-up of coolant in the protogalactic cores. While the photons clearly have sufficient energy ($`E_\gamma >0.755\mathrm{eV}`$) to photo-detach H-, we will show in Section §6 that H- photo-detachment is suppressed by at least two orders of magnitude relative to H2 formation at the temperatures, densities and UV background flux levels we consider.

## 3. The Sample of Simulated Halos

The locations of high density pregalactic clouds are identified in the simulations using the HOP algorithm (Eisenstein & Hut 1998) acting on the dark matter density distribution. We find that the high density peaks form within and at the intersection of sheets and filaments as expected in any cold dark matter dominated cosmology (Abel, et al. 1998). We choose the $`8`$$`10`$ most massive peaks of the $`200`$ identified by HOP at each redshift for further analysis. In Figure 1 we present the mass of these clouds as a function of redshift for the simulation without H2 . Since we use identical initial conditions for each of the four simulations, this is essentially the same set of peaks used for analysis in the simulations with H2 present and varying levels of the radiative background flux. These peaks span a mass range from $`3\times 10^4`$ to $`7\times 10^6M_{}`$ and are modeled in our simulations by $`1000`$ to $`125,000`$ dark matter particles, respectively. The absence of low mass peaks in the bottom portion of Figure 1 is an artifact of our selection criterion. However, the increase in the maximum mass of the cloud with decreasing redshift is real and characteristic of the “bottom up” growth of structure in these models. Small clouds form first at high redshift, merging to form more massive structures later. There is clear evidence for mergers among the clouds we consider. See, for example, the merger of the two most massive peaks in Figure 1 into a single peak between redshifts $`z=21`$ and $`z=20.5`$. However since we do not follow all of the low mass clouds in the simulation, we do not construct complete merger trees for the sample. Thus we caution the reader that all the peaks are not statistically independent of each other. This is particularly noticeable for the high mass peaks which are rare in our simulation volume and may result in a decrease in the scatter of measured properties for the higher mass clouds.

We define the virial radius, $`r_{vir}`$, as the radius of a sphere within which the average density in the cloud exceeds $`200`$ times the mean matter density of the universe. The virial mass, $`M_{vir}`$ is the total mass contained within that radius. Before cooling becomes important the mean gas-mass-weighted temperature of the cloud obeys a virial relationship between temperature, mass and redshift similar to that observed in galaxy clusters (Bryan & Norman 1998). For a neutral mix of primordial gas

$$T_{vir}(\mathrm{K})=155f(h^2\mathrm{\Delta }_cE^2)^{1/3}(\frac{M_{vir}}{10^7M_{}})^{2/3}$$
(5)
where $`f=0.8`$ is a normalization factor, $`h=0.65`$ is the dimensionless Hubble parameter at $`z=0`$, $`\mathrm{\Delta }_c=200`$ is the collapse overdensity and

$$E^2=\mathrm{\Omega }_0(1+z)^3+\mathrm{\Omega }_\mathrm{\Lambda }\mathrm{\Omega }_0(1+z)^3$$
(6)
for $`\mathrm{\Omega }_0=0.4`$ and $`\mathrm{\Omega }_\mathrm{\Lambda }=0.6`$ at redshifts $`30>z>19`$. In Figure 2 we compare the mean gas-mass-weighted temperature $`T_{mean}`$, for the clouds in our sample to the virial temperature $`T_{vir}`$ defined by Equation 5 for the case when no molecular hydrogen is present and for the case $`F_{LW}=0`$ when molecular hydrogen cooling is maximal. All of the clouds in our sample have virial temperatures $`T_{vir}<8000`$$`10,000\mathrm{K}`$, the temperature necessary for HI line cooling to become important. Peaks analyzed from the simulation without H2 can not cool and their mean temperatures are well characterized by the virial relation given in Equation 5. Once cooling becomes important, the mean gas-mass-weighted temperature of the cloud will drop significantly below its virial value as demonstrated in Figure 2 by peaks from the simulation that includes the H2 chemistry.

In Table 1 we compare the properties of the first protogalaxies to “collapse”, i.e reach maximal refinement, for the simulations with H2 cooling and soft UV background flux levels $`F_{LW}=0`$, $`10^{22}`$, and $`10^{21}`$ erg s-1 cm-2 Hz-1 . Protogalactic clouds in the simulation without H2 have not collapsed by redshift $`z=19`$. The first structures ($`F_{LW}=0`$) in this $`\mathrm{\Lambda }`$CDM model collapse early ($`z30`$) and are small ($`r_{vir}50`$ pc), low mass ($`M_{vir}1.4\times 10^5M_{}`$) clouds. Simulations of the minimum collapse mass by Fuller & Couchman (2000) find the minimum mass to collapse at $`z=30`$ to be $`3\times 10^5M_{}`$ (expected to be robust within about a factor of two) thus consistent with our results. Abel, Bryan, & Norman (2000) find a $`7\times 10^5M_{}`$ peak collapses first (at $`z=19`$) in their high resolution SCDM simulation. This too is consistent with the minimum mass collapse curves of Fuller & Couchman. Thus despite differences in numerical method and cosmology there appears to be general agreement that the minimum mass of the first star forming structures is $`10^5M_{}`$. This is about an order of magnitude below the earlier estimates of Tegmark et al. (1997) due primarily to the use of improved H2 cooling functions.

As the amount of the H2 photodissociating flux increases, the mass threshold for collapse increases pushing the redshift of collapse to lower values. However, even with flux levels as high as $`F_{LW}=10^{21}`$ erg s-1 cm-2 Hz-1 , H2 can survive in the cores of the clouds such that objects ($`M_{vir}1.26\times 10^6M_{}`$, $`T_{vir}2262`$ K ) below the mass threshold for HI line cooling are still able to collapse. However, the collapse in this case is delayed until $`z=21.5`$. This is consistent with the 1-D simulations of Haiman, Abel & Rees (1999). See, in particular, their Figure 6 for $`J_{21}=10^{1.1}`$. Furthermore the ability of a cloud to collapse does not seem to be strongly dependent on the angular momentum $`L`$ of its dark matter halo, since the threshold objects we found to collapse spanned the full range of expected spin parameters $`0.02<\lambda _{dm}<0.11`$ where $`\lambda =(L|E|^{1/2})/(GM^{5/2})`$, $`E`$ is the total energy, $`M`$ is the mass, and $`G`$ is the gravitational constant.

In order to check the numerical robustness of our results, we ran several lower mass resolution simulations ($`\mathrm{\Delta }M=306(38)M_{}`$ for dark matter(gas)) using the same random seed for the initial conditions and region for dynamical refinement as the high mass resolution runs presented here. We also ran several different realizations of the initial conditions at this lower resolution for the $`F_{LW}=10^{21}`$ erg s-1 cm-2 Hz-1 case to ensure that there was nothing exceptional about the particular realization we present. We found that, although the properties of the first peak to collapse were very similar for the two mass resolutions, the higher mass resolution run delayed the redshift of first collapse to slightly lower values ($`\mathrm{\Delta }z1`$). This difference in the redshift at which maximal refinement is reached, however, is already less than the cosmic variance ($`\mathrm{\Delta }z2`$) in this quantity observed between simulations run with different realizations of the initial conditions. The mass threshold for collapse was also reduced by $`30`$ % in the higher mass resolution run compared to the lower resolution simulation. This is slightly larger than the $`\pm 20`$ % variation in mass found from simulations with fixed mass resolutions but initial conditions generated using different random seeds.

## 4. Cold Gas Fraction

A critical input parameter into semi-analytical models that seek to include stellar as well as radiative feedback is the amount of gas in these first pregalactic objects that is available for star formation. Abel et al. (1998) estimated the fraction of gas able to cool by adding up all the mass within the radius of a $`4\sigma `$ peak (in three different simulations) for which the temperature in radially averaged profiles began to decrease. However, cooling alone is not sufficient to determine whether the gas is available to form stars. It must also become dense. We improve upon these previous estimates by measuring both the fraction of gas that is able to cool via molecular hydrogen cooling, $`f_c`$, and the fraction of gas that becomes both cold and dense, $`f_{cd}`$, in the clouds comprising the data sample described in Section §3. It is this latter fraction $`f_{cd}`$ of cold, dense gas that is available to form stars in these objects.

In Figure 3 we present both these statistics as functions of the cloud mass $`M_{vir}`$ for various levels of the background radiative flux $`F_{LW}`$. We define the following:

* $`f_c`$ is the fraction of total gas within the virial radius with temperature $`T<0.5T_{vir}`$ and gas density $`\rho >1000\rho _{mean}`$ where $`\rho _{mean}`$ is the mean gas density of the universe. This is the amount of gas within the cloud that has been able to cool due to molecular hydrogen cooling.
* $`f_{cd}`$ is the fraction of total gas within the virial radius with temperature $`T<0.5T_{vir}`$ and gas density $`\rho >10^{19}M_{}\text{Mpc}\text{-3}330`$ cm-3 . This is the fraction of gas within the cloud that is available for star formation.

The above criteria to determine the amount of gas in each category (cooled or both cold and dense) are applied on a cell by cell basis within the virial radius of each peak in the data sample. The temperature threshold has been chosen to ensure that the gas is substantially cooler than the virial temperature of the peak as given in Equation 5. The density threshold for cold, dense gas corresponds to the gas density at which the baryons become important to the gravitational potential and thus to the subsequent evolution of the core. However, setting a density threshold for gas that has been able to cool through the H2 chemistry but may not be dense is more subtle. We must be careful to exclude gas that is cool only because it belongs to lower mass clouds that are infalling onto the more massive structure. With the high mass resolution in our simulations, this substructure due to merging is resolvable. Since there is no cooling in the peaks extracted from the simulation without H2 , any cold gas present must be due to infall and merging in the outer regions of the cloud. We confirm that this is indeed true in the data set without H2 cooling by looking at radial profiles for the cold gas fraction in several peaks. We then used this data to determine the density threshold for cooled gas given above that first minimizes the cold infalling component. As shown in the bottom panel of Figure 3, the resultant background due to merging is negligible for all but the most massive peaks in the sample. The bump in $`f_c`$ in the high mass end of the bottom panel is due to the merger of the two most massive peaks in the box at $`z20.5`$ seen in Figure 1. Even in this extreme case the above density threshold keeps the infalling cold component of $`f_c`$ limited to a few percent. Note also that $`f_{cd}`$, the fraction of gas that can become dense and form stars, is zero for all peaks in the simulation without H2 , as expected.

Figure 3 represents one of the key results of this paper. Both the fraction of gas $`f_c`$ that can cool and the fraction $`f_{cd}`$ of cold, dense gas available for star formation increase logarithmically with the pregalactic cloud mass over the mass range dominated by H2 cooling. This dependence can be summarized in a simple fitting formula:

$$f_i=B_i\mathrm{ln}(M/M_{TH})$$
(7)
for $`i=c`$ (cooled gas) or $`i=cd`$ (cold, dense gas) and $`M>M_{TH}`$. The best fit slopes $`B_c`$ for the fraction of gas $`f_c`$ that cools are $`0.138\pm 0.006`$, $`0.139\pm 0.008`$, and $`0.085\pm 0.009`$ for radiative background fluxes of $`F_{LW}=0`$, $`10^{22}`$, and $`10^{21}`$ erg s-1 cm-2 Hz-1 , respectively, where the quoted errors are standard errors of the fit. However, the fraction of cold, dense gas $`f_{cd}`$ increases much less rapidly with mass. Best fit slopes $`B_{cd}`$ for this cold, dense component are $`0.058\pm 0.006`$, $`0.066\pm 0.007`$, and $`0.030\pm 0.004`$ for radiative background fluxes $`F_{LW}=0`$, $`10^{22}`$, and $`10^{21}`$ erg s-1 cm-2 Hz-1 , respectively. The equality of the slopes for flux levels $`J_{LW}=0`$ and $`10^{22}`$ erg s-1 cm-2 Hz-1 is striking and probably reflects the fact that the cooling chemistry is the same in each case independent of radiative flux level. The apparent softening of the slopes for $`F_{LW}=10^{21}`$ erg s-1 cm-2 Hz-1 is most likely an artifact of the small number of independent high mass clouds in our simulation volume. The mass thresholds $`M_{TH}`$ for nonzero $`f_c`$ and $`f_{cd}`$ are about equal at fixed $`F_{LW}`$ and increase as the level of soft UV flux is increased. This is a characterization of the negative feedback expected due to the photodissociation of H2 in the lower mass clouds. Our results for the mass threshold can be parameterized as:

$$M_{TH}(M_{})=1.25\times 10^5+8.7\times 10^5(\frac{F_{LW}}{10^{21}})^{0.47}$$
(8)
The presence of a soft H2 photodissociating flux at the levels considered here delays, but does not prevent gas in these pregalactic objects from becoming dense due to H2 cooling in their cores. The fraction of gas available for star formation due to H2 cooling can be represented as:

$$f_{cd}0.06\mathrm{ln}(M/M_{TH})$$
(9)
with $`M_{TH}`$ given by Equation 8.

## 5. Understanding Peak Formation and Collapse

In this section we use spherically averaged radial profiles of the pregalactic cloud properties to illustrate the internal dynamics of its collapse. We first consider the evolution with redshift of a collapsing pregalactic object when exposed to a fixed level of photodissociative flux. We then consider how the internal properties of a given cloud change at fixed redshift when the level of the photodissociating flux changes.

### 5.1. Evolution with redshift

In three of the panels of Figure §4 we show the evolution of the gas density, temperature, and H2 mass fraction from $`z=30`$ to maximal refinement at $`z=21.5`$ for the first object to collapse at a radiative flux level $`F_{LW}=10^{21}`$ erg s-1 cm-2 Hz-1 . (See Table 1). In the fourth panel (lower right) we show the evolution of the cooling and dynamical timescales for $`25z21.5`$ after cooling begins. Starting at large radius the upper set of curves represent the H2 cooling time $`t_{cool}`$ given by

$$t_{cool}=\frac{1.5n_gk_BT}{\mathrm{\Lambda }n_Hn_{H2}}$$
(10)
where $`n_g`$, $`n_H`$, and $`n_{H2}`$ are number densities for the gas, neutral hydrogen, and H2 respectively, $`k_B`$ is the Boltzmann constant, $`T`$ is the temperature, and $`\mathrm{\Lambda }`$ is the H2 cooling function. The lower set of curves at large radii show the dynamical timescale $`t_{dyn}`$ where

$$t_{dyn}=\left(\frac{3\pi }{16G\rho }\right)^{1/2}$$
(11)
for the same redshifts. The solid (dot-dashed) horizontal line is the Hubble time

$$t_{Hubble}=\frac{2}{3H_0}(\mathrm{\Omega }_0(1+z)^3)^{1/2}$$
(12)
for $`z=21.5(25)`$, respectively, and the heavy dashed line is the difference between the two, i.e. the elapsed time between the onset of cooling and collapse. The vertical line locates the virial radius $`r_{vir}=141`$ pc at $`z=21.5`$ when the structure has grown to a mass of $`1.26\times 10^6M_{}`$. The increases in density, temperature and H2 mass fraction seen outside the virial radius at $`1`$ kpc are due to averaging over other massive structures forming nearby.

During the initial stages of evolution (from $`30z25`$) H2 abundances are too low to be an effective coolant and the cloud grows from a mass of $`6.5\times 10^4M_{}`$ to $`5.3\times 10^5M_{}`$ through accretion and merging from within its host filament. The radial profiles within the core are well described by constant densities and temperatures that increase as expected as the mass of the object increases, while the outer edges of the profiles steepen with decreasing $`z`$. The H2 mass fractions are also roughly constant in the core (mirroring the density profiles) and grow with decreasing $`z`$. By $`z=25`$ the core temperature has risen to about $`2900\mathrm{K}`$, the core gas density has reached $`3\times 10^{17}M_{}\text{Mpc}\text{-3}10`$ cm-3, the H2 mass fraction is $`2.3\times 10^5`$ and the core has broadened from about $`10`$ pc to $`20`$ pc. For $`F_{LW}=10^{21}`$ erg s-1 cm-2 Hz-1 the photodissociation timescale

$$t_{diss}(s)=k_{diss}^1=9\times 10^9/F_{LW},$$
(13)
is $`9\times 10^{12}`$ s, much shorter than other timescales in the problem. Thus we can understand the H2 mass fraction in the core by a simple equilibrium argument. If we assume that H2 formation is dominated by the H- process given in Equations 1 and 2, the H2 equilibrium number density is given by

$$n_{H2}^{equil}\frac{k_7n_Hn_e}{k_{diss}},$$
(14)
where $`k_71.8\times 10^{18}T^{0.88}`$ cm3 s-1 over this temperature range, and $`n_H`$($`n_e`$) is the number density of neutral hydrogen (electrons), respectively. At $`z=25`$ Equation 14 predicts an equilibrium H2 mass fraction $`\chi _{H2}^{equil}2.6\times 10^5`$ in the core that is nearly identical to the value above found in our simulation.

By $`z=24`$ the H2 mass fraction has reached $`4\times 10^5`$ and gas in the central region has begun to cool. As cooling continues the core density steepens, evolving into a $`r^2`$ profile at $`z=21.5`$ in the inner $`10`$ pc. The gas density reaches number densities $`n4\times 10^3`$ exceeding the dark matter density for $`r<1`$ pc. Notice, however, that most of the cooling occurs in the region $`2<r<20`$ pc, outside this very dense central core. The temperature profile at collapse is very similar to that found by Abel, Bryan & Norman (2000) showing with decreasing radius cosmic infall, shocking, and a cooling flow. We also find for this object that the cooling time becomes comparable to the dynamical time over most of the cooling flow region just prior to collapse, and that the outer radius of the cooling region corresponds closely to the radius at which the cooling time has dropped below the elapsed time (heavy dashed line in Figure 4) in agreement with the collapse criterion of Haiman, Abel & Rees (1999).

Although spherically averaged radial profiles are easy to picture, it is important to remind ourselves that the evolution of these objects is fully three-dimensional and they do not collapse in isolation from their surroundings. In Figure 5 we show images of the log of the projected dark matter density, the baryonic density, and the temperature at $`z=21.5`$ for a portion of the dynamically refined region in the $`F_{LW}=10^{22}`$ erg s-1 cm-2 Hz-1 simulation. This region contains the two most massive structures in the simulation. The image is $`680`$ pc (proper) on a side and is projected along the y-axis. The masses (virial radii) for these protogalaxies are $`2.4\times 10^6M_{}`$($`174`$ pc) for the upper cloud and $`2.2\times 10^6M_{}`$($`170`$ pc) for the lower cloud, while their proper center-to-center separation is only $`300`$ pc. Thus their outer regions are already beginning to merge. The temperature map shows an intricate pattern of asymmetrical shocks with maximum temperatures around $`5400\mathrm{K}`$ surrounding elongated knots of cooled gas at the clouds’ centers. These cooling knots have temperatures of a few hundred Kelvin and are the regions of highest H2 concentration and thus gas density. The dark matter halos, containing $`42,000(37,000)`$ particles for the upper (lower) structures, respectively, clearly show substructure near their centers due to the prior infall and merger of lower mass clouds that have not yet completely relaxed. This merger activity is responsible in part for the aspherical geometry of the cooling knots and corresponding gas density.

### 5.2. Changing the photodissociating flux

In Figure 6 we compare the spherically averaged radial profiles of dark matter density, gas density, temperature and H2 mass fraction at fixed redshift ($`z=21.5`$) but with varying levels of H2 photodissociating flux for the $`2.2\times 10^6M_{}`$ (lower) structure of Figure 5. The cloud has only been able to maximally refine in the simulation with $`F_{LW}=0`$, although it is close to collapse in the $`F_{LW}=10^{22}`$ case also. This is in spite of the fact that the mass of this object is above the masses of the first objects to collapse in both the $`F_{LW}=10^{22}`$ and $`10^{21}`$ erg s-1 cm-2 Hz-1 simulations (See Table 1). Probably this is due to the strong effect of tidal interactions and high level of merging activity in and surrounding this object. The features located at $`300`$ pc and $`1`$ kpc in the profiles reflect the presence of the upper, $`2.4\times 10^6M_{}`$ protogalaxy of Figure 5 and the $`1.26\times 10^6M_{}`$ structure we analyzed earlier in Figure 4, respectively. The dark matter density profiles are essentially the same for all of the simulations, independent of H2 fraction or radiative flux. This is to be expected as long as the gravitational potential is dominated by the dark matter. The dark matter density only steepens away from this common profile in the inner $`2`$ pc of the $`F_{LW}=0`$ case when the gas density equals and then exceeds the dark matter density. For the $`F_{LW}=0`$ case the gas density steepens in the inner several parsecs with gas number densities reaching $`10^4`$, characteristic of the formation of a dense collapsed core. The cooling region extends inward from about $`40`$ pc until the temperature drops to the minimum value possible under H2 cooling, $`200\mathrm{K}`$, at around $`3`$ pc. The H2 mass fraction has saturated at $`10^3`$. The sharp drops in temperature that correspond to enhancements in the H2 mass fraction at $`9`$ and $`25`$ pc are the result of averaging over infalling clouds that were able to form H2 in their cores and cool before merging with this more massive object. As the radiative flux level is increased, the negative feedback becomes apparent. The maximum gas density drops in the central region and forms a homogeneous core whose core radius also increases as the flux is increased. The radius of the cooling region moves inward and cooling is suppressed due to the decreased amount of H2 coolant. The cloud is clearly in photodissociation equilibrium for the high flux ($`F_{LW}=10^{21}`$) case with the H2 fraction given by Equation 14. For the lower flux case, the H2 fraction is somewhat smaller than the photodissociation equilibrium value in the inner $`10`$ pc due to the influence of recombination on the formation of H- ions and thus the H2 formation rate. We also see, as expected, that there is no cooling at all for the case in which H2 is absent.

We can gain a better understanding of the above radial profiles and their scaling properties by comparing the relevant physical timescales in the problem. These times are shown for both the $`F_{LW}=10^{22}`$ case (upper panel) and the $`F_{LW}=10^{21}`$ case (lower panel) in Figure 7. In addition to profiles for the Hubble, dynamical, photodissociation, and cooling times defined earlier, we include profiles for the sound crossing time $`t_{cross}`$, the recombination time $`t_{rec}`$ given by

$$t_{rec}=7.7\times 10^9T^{0.64}/n_e,$$
(15)
and the H2 formation time $`t_{H2}`$ given by

$$t_{H2}=\frac{n_{H2}}{k_7n_Hn_e}.$$
(16)
All times are in seconds. We also show that the two-body relaxation time $`t_{relax}`$ (Equation 8-71 in Binney & Tremaine 1987) for the dark matter particles is long compared to all other timescales confirming that two-body numerical effects are unimportant.

For both cases the photodissociation and H2 formation times are much shorter than the Hubble time everywhere within the virial radius. In the high flux case they are also much shorter than the recombination time so that photodissociation equilibrium ($`t_{diss}=t_{H2}`$) is established throughout the cloud (see equation 14). Even so, enough H2 forms for the cooling time to drop below the Hubble time in the inner $`40`$ pc — core collapse is not far off. The low flux case is similar in the outer parts of the cloud: the recombination time is long and the H2 fraction is governed by its equilibrium value. However, within $`10`$ pc recombination becomes faster than photodissociation. This reduces the amount of $`H^{}`$ allowed to form and suppresses H2 formation below the photodissociation equilibrium value. However, enough H2 has formed that cooling and collapse proceed faster than photodissociation ($`t_{dyn}<t_{diss}`$ and $`t_{cool}<t_{diss}`$) at these small radii. We also find here — as in Abel, Bryan, & Norman (2000) — that $`t_{cross}<t_{dyn}`$ in the inner regions of the cloud so that the collapse proceeds hydrostatically.

## 6. Discussion

Now that we have detailed examples of the collapse properties of clouds in our simulations, we are better able to assess the importance of neglected processes in our simulations. In this section we discuss in more detail why the neglect of self-shielding and H- photodetachment should not significantly change our results. We then use a simple analytical model, similar to that by Tegmark et al. (1997), to illustrate the essential physics of the collapse.

### 6.1. Neglected processes

We first return to the issue of $`H_2`$ self-shielding. Although this effect could, in principle, be important for our objects, we argue here that its actual impact is slight. The maximum column densities of H2 we encounter in our simulations are around $`10^{17}`$ cm-2, while a column density of $`5\times 10^{14}`$ cm-2 is sufficient to start blocking the soft-UV flux. So, at first glance, it appears that self-shielding should be important; however, this value has been derived for a static gas distribution. The Lyman-Werner bands consist of a large number of individual lines whose width is dominated (in our case) by Doppler broadening.

There are substantial velocity gradients within all of our simulated objects, even those which appear to be relaxed. These take the form not just of radial inflow, but also large-scale disordered velocity fields (i.e. turbulence). Typical rms line-of-sight velocities are around 3-4 km/s. These values are larger than the $``$2 km/s line widths due to thermal Doppler broadening. This means that even if a given parcel of gas is behind a large column of molecular hydrogen, it may still be photo-dissociated if it has a velocity shift relative to the blanketing material. There is very little continuum opacity.

It is a difficult radiative transfer calculation to determine the precise effect of these velocity shifts; however we may consider two limiting cases. The first is the previously mentioned static case in which velocities are negligible compared to the thermal Doppler widths ($`v/b1`$) and an $`H_2`$ column density of $`5\times 10^{14}`$ cm-2 is sufficient to start affecting the photo-dissociation rate. In the other extreme, in which $`v/b1`$, the rate is not significantly reduced until a substantial fraction of the 11-13 eV range is blocked by the Lyman-Werner lines. Since this does not occur until the damping wings of the lines become important, it requires a much larger molecular hydrogen column density, of order $`10^{22}`$ cm-2 (Draine & Bertoldi 1996; Glover & Brand 2000). In the simulations, $`v/b>1`$ so we are likely to be closer to the second regime than the first. Even if it were marginally important, there are three further reasons to think that self-shielding is unimportant. First, these large column densities occur at late times, once the core has already collapsed. For most of the evolution the objects are clearly in the optically thin regime. Second, the maximum column density quoted above is to the center of the object, but the cooling front lies outside of the $`H_2`$ core, and so beyond most of the covering material. Finally, once a sufficient amount of $`H_2`$ has formed to produce a substantial blanketing column, the cooling time has dropped below the photo-dissociation time, so that photo-dissociation is no longer important (and thus even if the soft-UV flux is attenuated, the subsequent evolution of the gas would be unchanged). So while the results described earlier are properly lower limits to the cooled fraction because of our neglect of self-shielding, a full radiative transfer calculation should produce nearly identical results.

Let us also reexamine whether H- photo-detachment (which we neglect) could significantly affect the results. In order to bracket the size of this effect we consider two spectral shapes for the radiation field, $`J\nu ^\alpha `$ with $`\alpha =1`$ and $`1`$. The softer spectrum, $`\alpha =1`$, is the spectral slope expected in this energy range from a typical OB star. The steeper slope, $`\alpha =1`$, is the type of spectral shape one might expect for a background radiation field built up from a number of sources. The cross section (in cm2) for H- photo-detachment is given by (Kang & Shapiro 1987, Abel et al. 1997),

$$\sigma _{detach}=7.928\times 10^5(\nu \nu _{th})^{3/2}\frac{1}{\nu ^3}\mathrm{for}\nu >\nu _{th}$$
(17)
where $`h_P`$ is Planck’s constant and $`h_P\nu _{th}=0.755`$ eV. The H- photo-detachment rate coefficient $`k_{detach}`$ can then be calculated

$$k_{detach}=_{\nu _{th}}^{\nu _H}\frac{4\pi J(\nu )}{h_P\nu }\sigma _{detach}(\nu )𝑑\nu $$
(18)
with $`h_P\nu _H=13.6`$ eV. For the two spectral shapes given above we find $`k_{detach}=1.9\times 10^{12}F_{LW}`$ s-1 for $`\alpha =1`$ and $`6.4\times 10^{11}F_{LW}`$ s-1 for $`\alpha =1`$ where $`F_{LW}`$ is the radiation flux in erg s-1 cm-2 Hz-1 .

The dominant processes other than photodetachment that determine the abundance of H- are shown in Equations 1 and 2 with rate coefficients $`k_7`$ and $`k_8`$, respectively. Furthermore the processes proceed sufficiently rapidly compared to the H2 chemistry that $`n_H`$ is well approximated by its equilibrium value (Abel et al. 1997) such that (including photo-detachment):

$$n_H\frac{k_7n_Hn_e}{k_8n_H+k_{detach}}.$$
(19)
We see that H- photo-detachment can be neglected as long as $`k_{detach}k_8n_H`$. This leads to a condition on the hydrogen number density $`n_Hk_{detach}/k_8`$ that can be checked against the conditions found within the collapsing structures. For gas temperatures found in our simulations, $`k_81.43\times 10^9`$ cm3s-1 (Abel et al. 1997). Thus photo-detachment is unimportant provided that $`n_H0.045(F_{LW}/10^{21})`$ for the steeper $`\alpha =1`$ spectrum or, a much weaker constraint, that $`n_H1.4\times 10^3(F_{LW}/10^{21})`$ for the softer $`\alpha =1`$ spectrum. Let us check these constraints against the hydrogen number densities found within the first peak to collapse at our highest flux level ($`F_{LW}=10^{21}`$ erg s-1 cm-2 Hz-1 ) shown in Figure 4. At $`z=25`$, before cooling commences, hydrogen number densities in the core region where H2 will form already exceed $`n_H7`$. This number density satisfies the constraint for the steepest spectrum ($`\alpha =1`$) by more than two orders of magnitude so that the neglect of H- photo-detachment does not significantly affect our results.

### 6.2. An Analytic Estimate of the Collapse Threshold

Analytical arguments often provide qualitative insight into the key physical principles at work in a complex process. Since we have seen that the H2 fraction is close to photodissociation equilibrium throughout our objects, we use these equilibrium abundances to develop analytical estimates of the critical H2 fraction and mass thresholds for collapse. We adopt the simplified collapse criterion, similar to that of Tegmark et al. (1997), that an object will collapse if the cooling time $`t_{cool}`$, given by Equation 10, is less than the Hubble time $`t_{Hubble}`$, defined in Equation 12. Thus

$$\frac{1.5n_gk_BT}{\mathrm{\Lambda }n_Hn_{H2}}\frac{2}{3H_0}(\mathrm{\Omega }_0(1+z)^3)^{1/2}$$
(20)
The equality gives the critical number density of H2 necessary for collapse. Solving for $`n_{H2}^{crit}`$ we obtain

$$n_{H2}^{crit}=2.25\frac{n_gk_BT}{n_H\mathrm{\Lambda }}H_0\mathrm{\Omega }_0^{1/2}(1+z)^{3/2}$$
(21)
In this equation, $`\mathrm{\Lambda }`$ is, in our case, the Lepp & Shull (1983) cooling function. In the low density limit ($`n_{H2}10^4`$ cm-3 ) for temperatures $`300\mathrm{K}T5000\mathrm{K}`$ of interest here it is well approximated by $`\mathrm{\Lambda }6.0\times 10^{11}Tk_Be^{730/T}`$ such that

$$n_{H2}^{crit}1.3\times 10^7he^{730/T}\mathrm{\Omega }_0^{1/2}(1+z)^{3/2}$$
(22)
If we now equate this critical abundance of H2 required to satisfy our collapse criterion to the photodissociation equilibrium abundance from Equation 14 and solve for $`F_{LW}`$, we obtain an expression for the critical flux that just allows a cloud with virial temperature $`T`$ to collapse.

$$\left(\frac{F_{LW}^{crit}}{10^{21}}\right)\frac{126n_H^2x_eT^{0.88}e^{730/T}}{h\mathrm{\Omega }_0^{1/2}(1+z)^{3/2}}.$$
(23)
where the electron fraction $`x_e=n_e/n_H`$.

To complete the above relationship we need to approximate the electron fraction and the gas number density in the cloud core. If the collapse occurs well within one recombination time, the electron fraction, $`x_e`$, is expected to be close to its residual value $`x_e=1.2\times 10^5\mathrm{\Omega }_0^{1/2}/(\mathrm{\Omega }_bh)`$. The average baryon number density in the universe is given by $`n_B9.3\times 10^6(1+z)^3\mathrm{\Omega }_bh^2`$ cm-3 . The gas temperature after thermal decoupling from the CMB at $`z200`$ but before reionization $`z>7`$ drops adiabatically with $`T_{IGM}0.0135(1+z)^2\mathrm{K}`$. If we assume that the baryons in the central regions of dark matter halos were never shocked, we find the maximum density that can be reached by adiabatic evolution is $`n_{max}=n_B(T_{vir}/T_{IGM})^{1/(\gamma 1)}`$. Here $`T_{vir}`$ denotes the virial temperature of the dark matter halo and $`\gamma =5/3`$ is the ratio of specific heats for an ideal gas. This gives a collapse redshift independent maximal density as a function of virial temperature for adiabatic evolution of the cloud as $`n_{max}187\mathrm{\Omega }_bh^2(T/1000)^{3/2}`$ cm-3 . We expect adiabatic evolution to be a good approximation for low mass clouds, but caution that for larger mass structures merging and shock heating become important. Inserting the expressions for $`n_{max}`$ and $`x_e`$ into Equation 23 we obtain:

$$\left(\frac{F_{LW}^{crit}}{10^{21}}\right)223\mathrm{\Omega }_bh^2\left(\frac{T}{1000}\right)^{3.88}e^{730/T}\left(\frac{1+z}{20}\right)^{3/2}$$
(24)
The structures that are first able to collapse in our simulations all have virial temperatures between $`700\mathrm{K}`$ and $`3000\mathrm{K}`$ where the exponential in Equation 24 can be approximated by a power law, $`e^{730/T}0.48(T/1000)^{0.7}`$. Inserting this approximation into Equation 24 and inverting, we obtain an analytic approximation for the critical virial temperature for collapse at a given soft-UV flux level.

$$\left(\frac{T_{vir}^{crit}}{1000}\right)0.36\left((\mathrm{\Omega }_bh^2)^1\left(\frac{F_{LW}}{10^{21}}\right)\left(\frac{1+z}{20}\right)^{3/2}\right)^{0.22}$$
(25)
For $`F_{LW}=10^{21}`$ erg s-1 cm-2 Hz-1 at $`z=25`$ in our cosmological model ($`h=0.65`$, $`\mathrm{\Omega }_b=0.05`$) this gives a $`T_{vir}^{crit}920\mathrm{K}`$, somewhat lower than but comparable to the $`T_{vir}1470\mathrm{K}`$ found at that redshift for the first peak to reach maximal refinement at that flux level in our simulation (See Figure 4). The disagreement in detail is not surprising since the approximation invokes a simplified collapse criterion and ignores the nonadiabatic evolution of the collapse. However, the analytical argument does serve to ellucidate the essential physics at work in the collapse.

## 7. Summary

In this paper, we have used high-resolution numerical simulations to investigate the effect of radiative feedback on the formation of $`10^510^7M_{}`$ clouds. This range of objects is important because they are large enough to form molecular hydrogen, but too small to cool by hydrogen line cooling. Therefore they are susceptible to the destruction of molecular hydrogen from soft-UV flux in the 11-13 eV range (the Lyman-Werner bands). This radiation comes from the very first generation of stars and so its amplitude is unclear; we test four cases, ranging from zero flux to complete H2 destruction. The initial conditions are drawn from a $`\mathrm{\Lambda }`$CDM cosmological model, and the simulations evolve the non-equilibrium rate equations for twelve species of hydrogen and helium. Our main conclusions may be summarized as follows.

* A soft-UV flux does delay cooling and star formation in the cores of clouds in this mass range. An increasing amount of radiation delays collapse until later redshifts when larger objects have collapsed. However, even without any radiative feedback, objects with masses less than $`1.4\times 10^5M_{}`$ (at $`z30`$) do not cool efficiently due to their low H2 fraction. We provide a fit to the mass threshold for collapse, given the mean flux in the Lyman-Werner bands (equation 8).
* The amount of cold, dense gas which is available for star formation in the cores of the collapsed objects depends primarily on two numbers, with relatively little scatter: the flux of soft UV radiation and the mass of the cloud. We show the cold, dense fraction increases logarithmically with the mass of the halo and we provide a fit for this relation as well (equation 9).
* The evolution of these objects is relatively straightforward. With a negligible soft-UV radiation field, the formation of H2 continues until it is suppressed by the removal of electrons from the core due to recombination. However, by the time this occurs, generally enough H2 has formed to allow core collapse. When radiative feedback is important ($`F_{LW}5\times 10^{24}`$ erg s-1 cm-2 Hz-1 ), the destruction and formation processes are sufficiently fast so that the H2 fraction is in equilibrium (equation 14) and recombination plays little role until after core collapse begins. This allows the derivation of an analytic estimate for the collapse threshold which comes close to matching the full simulation results.

Finally, there are a number of other results which deserve further investigation. One of these is that the collapse redshift does not seem to depend on the angular momentum of the halo (see section §3). Another is that we do not appear to witness fragmentation in these objects. The few times in which we clearly see more than one collapsed core in a halo appears to be due to the merging of a number of halos, each of which already contained a core.

We also examined a number of process which we omitted from the simulation, including H2 self-shielding and $`H^{}`$ photo-detachment, none of which appears to be important for this calculation. These results were derived within the context of a currently popular $`\mathrm{\Lambda }`$-dominated cosmological model. However, at these high redshifts, the effective $`\mathrm{\Omega }`$ should be close to one, so the most important parameter is the amplitude of the power spectrum on these scales. Relatively small changes to the amplitude (such as would arise from the current uncertainty in $`\sigma _8`$) are unlikely to have a major impact on the results presented here, although it should always be kept in mind that we are extrapolating the power spectrum to scales significantly smaller than constrained by observations.

This work is supported in part by NSF grant AST-9803137 under the auspices of the Grand Challenge Cosmology Consortium. NASA also supported this work through Hubble Fellowship grant HF-0110401-98A from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc. under NASA contract NAS5-26555. The computations used the SGI Origin2000 at the National Center for Supercomputing Applications. We acknowledge useful discussions with Zoltan Haiman.


format: text

text:
# Experimental realization of order-finding with a quantum computer

## Abstract

We report the realization of a nuclear magnetic resonance (NMR) quantum computer which combines the quantum Fourier transform (QFT) with exponentiated permutations, demonstrating a quantum algorithm for order-finding. This algorithm has the same structure as Shor’s algorithm and its speed-up over classical algorithms scales exponentially. The implementation uses a particularly well-suited five quantum bit molecule and was made possible by a new state initialization procedure and several quantum control techniques.

The quest for the experimental realization of quantum computers has culminated in the creation of specific entangled quantum states, most recently with four quantum bits (qubits) using trapped ions, and seven qubits using liquid state nuclear magnetic resonance (NMR), and in the successful implementation of Grover’s search algorithm  and the Deutsch-Jozsa algorithm  on two, three, and five qubit systems (see  for additional references).

However, a key step which remains yet to be taken is a computation with the structure of Shor’s factoring algorithm, which appears to be common to all quantum algorithms that achieve exponential speed-up . This structure involves two components: exponentiated unitary operations and the quantum Fourier transform (QFT). Implementing these components is challenging because they require not just the creation of static entangled states , but also precise dynamic quantum control over the evolution of multiple entangled qubits, over the course of tens to hundreds of quantum gates for the smallest meaningful instances of this class of algorithms. The evolution of the states is precisely where NMR quantum computers appear to have an exponential advantage over classical computers .

Here we report the experimental implementation of a quantum algorithm for finding the order of permutation; its structure is the same as for Shor’s factoring algorithm and it scales exponentially faster than any classical algorithm for the problem. The realization of this algorithm was made possible by the synthesis of an unusual molecule with five pairwise coupled, easily addressable 19F spins, and by the introduction of two new techniques: an efficient and effective temporal labeling scheme for initial state preparation, and a method for precise simultaneous rotations of multiple spins at nearby frequencies .

The order of a permutation $`\pi `$ can be understood via the following analogy: imagine $`2^n`$ rooms and $`2^n`$ one-way passages connecting the rooms, with exactly one entrance and one exit in each room (for some rooms, the passage going out may loop back to the room itself). These rules ensure that when making transitions from one room to the next going through the passages, you must eventually come back to the room you started from. Define the order $`r`$ as the minimum number of transitions needed to return to the starting room $`y`$, where $`r`$ may depend on $`y`$. The order-finding problem is to determine $`r`$ solely by trials of the type “make $`x`$ transitions using $`\pi `$ starting from room $`y`$ and check which room you are in”. Mathematically, we will describe such trials as queries of an oracle or black box which outputs $`\pi ^x(y)`$. The goal then is to find $`r`$ with the least possible number of queries.

Cleve showed that order-finding using no other resource or information than the oracle is hard both for deterministic and probabilistic classical computers , i.e. there exists a lower bound on the number of oracle queries needed for order-finding which is exponential in $`n`$. In contrast, this problem can be solved much faster on a quantum computer, because finding the order of $`\pi (y)`$ is equivalent to finding the period of the function $`f(x)=\pi ^x(y)`$. The latter can be done with a constant probability of success in a constant number of function evaluations using a generalization of Shor’s quantum algorithm . This is because, in some sense, the quantum computer can make transitions to many rooms at once. Thus, in terms of the number of oracle queries required, the gap between the quantum and classical case is exponential .

We experimentally implemented the order-finding quantum algorithm to determine the order of a representative subset of all $`4!=24`$ permutations on $`4`$ elements, including instances of each possible order. It can be proven that the best classical algorithm needs two queries of the oracle to determine $`r`$ with certainty, and that using only one query of the oracle, the probability of finding $`r`$ using a classical algorithm can be no more than $`1/2`$. One optimal classical strategy is to first ask the oracle for the value of $`\pi ^3(y)`$: when the result is $`y`$, $`r`$ must be $`1`$ or $`3`$; otherwise $`r`$ must be $`2`$ or $`4`$. In either case, the actual order can be guessed only with probability $`1/2`$. In contrast, the probability of success is $`0.55`$ with only one oracle-query using the quantum algorithm on a single quantum computer. In fact, since in our implementation an ensemble of $`10^{18}`$ quantum computers contribute to the signal, our output data enables the order to be deduced with virtual certainty.

The quantum algorithm is as follows (see Fig. 1 for the quantum circuit): (0) initialize the first register of three qubits in the state $`|0`$ and set the second register of two qubits to $`|y_1y_0`$ or for short $`|y`$, where $`y_1y_0`$ is the binary representation of the number $`y`$; (1) apply a Hadamard transform $`H`$ to qubits $`1`$, $`2`$ and $`3`$, which puts the first register in the state $`|x=(|0+|1+\mathrm{}+|7)/\sqrt{8}`$; (2) apply the unitary transformation $`|x|y|x|\pi ^x(y)`$, which is the oracle query; (3) perform the quantum Fourier transform (QFT) on the first three qubits ; (4) measure the first three qubits — for an ideal single quantum computer, the possible measurement outcomes $`m`$ and their probabilities are listed in Fig. 2 for each possible value of $`r`$ (the translation to ensemble averaged measurements is discussed later); (5) depending on the measurement outcome, make a probabilistic guess $`r^{}`$ as shown also in Fig. 2. It is easy to verify that Pr\[$`r^{}=r`$\] is $`0.55`$, regardless of the probability distribution of $`r`$ or $`\pi `$.

In order to implement this algorithm, we custom synthesized a molecule  containing five 19F spins which served as the qubits (Fig. 3). When placed in a static magnetic field, each spin has two discrete energy eigenstates, spin-up, $`|0`$, and spin-down, $`|1`$, described by the Hamiltonian $`\mathrm{}\omega _iI_{zi}`$, where $`\omega _i`$ is the transition frequency between the spin-up and spin-down states and $`I_z`$ is the $`\widehat{z}`$ component of the spin angular momentum operator. In this molecule, all five spins are remarkably well-separated in frequency, $`\omega _i`$, and are mutually coupled with a coupling Hamiltonian of the form $`2\pi \mathrm{}J_{ij}I_{zi}I_{zj}`$ (Fig. 3). The linewidths of the NMR transitions are $``$ 1 Hz, so the $`T_2`$ quantum coherence times of the spins were at least $``$ 0.3 s. The $`T_1`$ time constants were measured to be between $`3`$ and $`12`$ s.

This nuclear spin system was used at room temperature. The thermal equilibrium state is then highly mixed, i.e. the probabilities that each spin is $`|0`$ or $`|1`$ differ by only $`1`$ part in $`10^5`$, which is not a suitable initial state for a quantum computation . Instead, as has previously been shown, an “effective pure” initial state must be created, in which only the spins in the $`|00000`$ state produce a net output signal . We devised a new procedure to prepare effective pure states which is best understood in terms of product operators . The equilibrium density matrix for a homonuclear spin system is a sum of $`n=5`$ terms: $`IIIIZ+IIIZI+IIZII+IZIII+ZIIII`$. The desired effective pure state density matrix is $`IIIIZ+\mathrm{}+ZIIII+IIIZZ+\mathrm{}+ZZIII+IIZZZ+\mathrm{}+ZZZII+IZZZZ+\mathrm{}+ZZZZI+ZZZZZ`$, a sum of $`2^n1=31`$ terms. Using short sequences of controlled-NOT operations ($`\mathrm{C}_{ij}`$ flips spin $`j`$ if and only if $`i`$ is $`|1`$), the five terms obtained in equilibrium can be transformed into different sets of five terms. For homonuclear spin systems, the summation of only $`(2^n1)/n=7`$ different experiments thus suffices to create all $`31`$ terms, although it may be advantageous to use slightly more experiments in order to keep the preparation sequences short. In contrast, both conventional temporal averaging  and later improvements  require up to $`2^n1`$ experiments, and furthermore suffer from higher complexity and/or a lower signal-to-noise ratio. The overhead is still exponential though, so even this improved technique is not scalable. But importantly, a scalable approach to NMR quantum computation exists  and may become practical if large polarization enhancements can be achieved. We used $`9`$ experiments, giving a total of $`45`$ product operator terms. The $`14`$ extra terms were canceled out pairwise, using NOT ($`\mathrm{N}_i`$) operations to flip the sign of selected terms. The $`9`$ state preparation sequences were

$`\mathrm{C}_{51}\mathrm{C}_{45}\mathrm{C}_{24}\mathrm{N}_3,\mathrm{C}_{14}\mathrm{C}_{31}\mathrm{C}_{53}\mathrm{N}_2,\mathrm{C}_{54}\mathrm{C}_{51}\mathrm{N}_2,\mathrm{C}_{31}\mathrm{C}_{43}\mathrm{C}_{23}\mathrm{N}_5,\mathrm{C}_{21}\mathrm{C}_{52}\mathrm{C}_{45}\mathrm{C}_{34},\mathrm{C}_{53}\mathrm{C}_{25}\mathrm{C}_{12}\mathrm{N}_4,\mathrm{C}_{12}\mathrm{C}_{15}\mathrm{C}_{13}\mathrm{C}_{41},\mathrm{C}_{32}\mathrm{C}_{13}\mathrm{C}_{25}\mathrm{N}_4,\mathrm{C}_{35}\mathrm{C}_{23}\mathrm{N}_1.`$

The resulting data are remarkably clean, as illustrated in Fig. 4. After the state preparation, only the $`0000`$ line should remain visible, reflecting that only molecules with all spins in the ground state contribute to the signal. This is clearly observed in the measured spectra.

This effective pure state served as the initial state for the order-finding algorithm. The actual computation was realized via a sequence of $`50`$ to $`200`$ radio-frequency (RF) pulses, separated by time intervals of free evolution under the Hamiltonian, for a total duration of $`50`$ to $`500`$ ms, depending on $`\pi `$. The pulse sequences for the order-finding algorithm were designed by translating the quantum circuits of Fig. 1 into one- and two-qubit operations, employing several simplification methods . These pulse sequences were implemented on a custom modified four-channel Varian Unity INOVA spectrometer, and a Nalorac HFX probe. The frequency of one channel was set at $`(\omega _2+\omega _3)/2`$, and the other three channels were set on the resonance of spins $`1,4`$ and $`5`$. The chemical shift evolutions of spins $`2`$ and $`3`$ were calculated with the help of a time-counter, which kept track of the time elapsed from the start of the pulse sequence. On-resonance excitation of spins $`2`$ and $`3`$ was achieved using phase-ramping techniques . All pulses were spin-selective and Hermite shaped . Rotations about the $`\widehat{z}`$-axis were implemented by adjusting the phases of the subsequent pulses . Unintended phase shifts  of spins $`i`$ during a pulse on spin $`ji`$ were calculated and accounted for by adjusting the phase of subsequent pulses. During simultaneous pulses, the effect of these phase shifts was largely removed by shifting the frequency of the pulses via phase-ramping. The pulse frequency shifts are designed such that they track the shifting spin frequencies and thereby greatly improve the accuracy of the simultaneous rotations of two or more spins . This technique circumvents the need to avoid simultaneous pulses at nearby frequencies , and thus permits more efficient pulse sequences.

Upon completion of the pulse sequence, the states of the three spins in the first register were measured and the order $`r`$ was determined from the read-out. Since an ensemble of quantum computers rather than a single quantum computer was used, the measurement gives the bitwise average values of $`m_i(i=1,2,3)`$, instead of a sample of $`m=m_1m_2m_3`$ with probabilities given in Fig. 2 . Formally, measurement of spin $`i`$ returns $`O_i=12m_i=2\text{Tr}(\rho I_{zi})`$, where $`\rho `$ is the final density operator of the system. The $`O_i`$ are obtained experimentally from integrating the peak areas in the spectrum of the magnetic signal of spin $`i`$ after a $`90^{}`$ read-out pulse on spin $`i`$, phased such that positive spectral lines correspond to positive $`O_i`$. The theoretically predicted values of $`O_i`$ ($`i=1,2,3`$) for each value of $`r`$ follow directly from the probabilities for $`m`$ in Fig. 2. For reference, we also include the values of $`O_4`$ and $`O_5`$ (for $`y=0`$; if $`y0`$, $`O_4`$ and $`O_5`$ can be negative): for the case $`r=1`$ the $`O_i`$ are $`1,1,1,1,1`$; for $`r=2`$ they are $`1,1,0,1,0`$; and for $`r=4`$ they are $`1,0,0,0,0`$. For $`r=3`$, the $`O_i`$ ($`i=1,2,3`$) are $`0,1/4,5/16`$, and $`O_4`$ and $`O_5`$ can be $`0`$, $`\pm 1/4`$ or $`\pm 1/2`$, depending on $`y`$. The value of $`r`$ can thus be unambiguously determined from the spectra of the three spins in the first register. This was confirmed experimentally by taking spectra for these three spins, which were in excellent agreement with the theoretical expectations.

In fact, the complete spectrum of any one of the first three spins uniquely characterizes $`r`$ by virtue of extra information contained in the splitting of the lines. For the spectrum of spin $`1`$ the values of $`O_i`$ given above indicate that for $`r=1`$, only the $`0000`$ line (see Fig. 4) will be visible since spins $`25`$ are all in $`|0`$. Furthermore, this line should be absorptive and positive since spin $`1`$ is also in $`|0`$. Similarly, for $`r=2`$ the $`0000,0001,0100`$ and $`0101`$ lines are expected to be positive, and for $`r=4`$ all $`16`$ lines should be positive. Finally, for $`r=3`$, the net area under the lines of spin $`1`$ should be zero since $`O_1=0`$, although most individual lines are expected to be non-zero and partly dispersive. These unambiguous characteristics are reflected in the data. Results for four representative permutations are presented in Fig. 5. In all cases, the spectrum is in good agreement with the predictions, both in terms of the number of lines, and their position, sign and amplitude. Slight deviations from the ideally expected spectra are attributed mostly to incomplete refocusing of undesired coupled evolutions and to decoherence.

The success of the order-finding experiment required the synthesis of a molecule with unusual NMR properties and the development of several new methods to meet the increasing demands for control over the spin dynamics. The major difficulty was to address and control the qubits sufficiently well to remove undesired couplings while leaving select couplings active. Furthermore, the pulse sequence had to be executed within the coherence time. Clearly, the same challenges will be faced in moving beyond liquid state NMR, and we anticipate that solutions such as those reported here will be useful in future quantum computer implementations, in particular in those involving spins, such as solid state NMR , electron spins in quantum dots  and ion traps .

We thank X. Zhou, A. Verhulst, M. Sherwood, S. Smallcombe, A. Brooke and R. Laflamme for help and discussions, J. Harris and W. Risk for support, and the Aspen Center for Physics for its hospitality. L.V. gratefully acknowledges support by a Yansouni Family Stanford Graduate Fellowship. This work was supported by DARPA under the NMRQC initiative.

Email address: lieven@snow.stanford.edu.