---
abstract: 'We report results from numerical simulations of star formation in the early universe that focus on the dynamical behavior of metal-free gas under different initial and environmental conditions. In particular we investigate the role of turbulence, which is thought to ubiquitously accompany the collapse of high-redshift halos. We distinguish between two main cases: the birth of Population III.1 stars – those which form in the pristine halos unaffected by prior star formation – and the formation of Population III.2 stars – those forming in halos where the gas has an increased ionization fraction. We find that turbulent primordial gas is highly susceptible to fragmentation in both cases, even for turbulence in the subsonic regime, i.e. for rms velocity dispersions as low as 20 % of the sound speed. Fragmentation is more vigorous and more widespread in pristine halos compared to pre-ionized ones. If such levels of turbulent motions were indeed present in star-forming minihalos, Pop III.1 stars would be on average of somewhat lower mass, and form in larger groups, than Pop III.2 stars. We find that fragment masses cover over two orders of magnitude, suggesting that the Population III initial mass function may have been much broader than previously thought. This prompts the need for a large, high-resolution study of the formation of dark matter minihalos that is capable of resolving the turbulent flows in the gas at the moment when the baryons become self-gravitating. This would help to determine the applicability of our results to primordial star formation'
author:
- 'Paul C. Clark$^1$, Simon C.O. Glover$^1$, Ralf S. Klessen$^{1, 2}$ & Volker Bromm$^3$'
title: Gravitational fragmentation in turbulent primordial gas and the initial mass function of Population III stars
---

Introduction {#sec:intro}
============

Over the course of the last decade, work by a number of groups has led to the development of a widely accepted picture for the formation of the first stars, the so-called Population III (or Pop. III) stars. In this picture, the very first stars (Population III.1 in the nomenclature of @tm08 and @omha2008) form within small dark matter halos that have total masses $M \sim 10^{6} \: {\rm M_{\odot}}$, virial temperatures of around a thousand K, and that are assembled at redshifts $z \sim 25$–30 or above [@bl04; @g05; @byhm09]. Gas falling into one of these small dark matter halos is shock-heated to a temperature close to the virial temperature of the halo, and thereafter cools via H$_{2}$ rotational and vibrational line emission. The H$_{2}$ that enables the gas to cool is primarily formed via the gas-phase reactions $$\begin{aligned}
{\rm H} + {\rm e^{-}} & \rightarrow & {\rm H^{-}} + \gamma, \\
{\rm H^{-}} + {\rm H} & \rightarrow & {\rm H_{2}} + {\rm e^{-}},\end{aligned}$$ where the required free electrons are those that remain in the gas after cosmological recombination at $z \sim 1100$. The low abundance of these electrons, plus the fact that the gas starts to recombine further once it collapses into the dark matter halo, limits the amount of H$_{2}$ that can form in this manner; the typical fractional abundances of H$_{2}$ that result are a few times $10^{-3}$. The limited H$_{2}$ abundance, the low critical density of the H$_{2}$ molecule, and the large energy separation of its lowest accessible rotational levels combine to significantly limit the extent to which the gas can cool. The minimum temperature reached depends upon the dynamical history of the collapse, but is typically around $200 \: {\rm K}$. This minimum temperature is reached at a density of around $10^{4}
\: {\rm cm^{-3}}$, comparable to the critical density of H$_{2}$ (the density at which its level populations reach their local thermodynamical equilibrium (LTE) values), and at higher densities, the gas begins to reheat. Gas falling into the dark matter halo therefore tends to fragment once it reaches this temperature and density, and the resulting fragments have masses of the order of the local Jeans mass, $M_{\rm J} \sim 1000 \: {\rm M_{\odot}}$.[^1] Following this period of fragmentation, the gas in the fragments undergoes a further period of gravitational collapse, eventually reaching a density of $n \sim 10^{8}$–$10^{10} \: {\rm cm^{-3}}$, at which point three-body H$_{2}$ formation becomes effective. This rapidly converts almost all of the available atomic hydrogen into H$_{2}$, but the significant energy release that accompanies this process heats the gas, which typically attains a temperature of 1000–2000 K during this phase of the collapse.

Continued collapse next leads to the gas becoming optically thick in the main H$_{2}$ cooling lines (at $n \sim 10^{10} \: {\rm cm^{-3}}$; @ra04, @yoha06), the onset of collision-induced emission cooling (at $n
\sim 10^{14} \: {\rm cm^{-3}}$; @ra04), and finally to the gas becoming optically thick in the continuum as well as in the H$_{2}$ lines ($n \sim 10^{16}
\: {\rm cm^{-3}}$; @yoh08). H$_{2}$ dissociation cooling allows for further collapse, but once all of the H$_{2}$ is gone, the collapse becomes fully adiabatic. Detailed simulations of the collapse, using adaptive mesh refinement (AMR; see e.g. @abn02 [@on07]) or smoothed particle hydrodynamics (SPH; see e.g. @brlb04 [@yoha06; @yoh08]) show little or no evidence for fragmentation between $n \sim 10^{4} \: {\rm cm^{-3}}$ and the onset of this adiabatic evolution, and it is natural to associate the latter with the formation of a primordial protostar. The amount of gas incorporated into the protostar at this point is small, approximately $0.01 \: {\rm M_{\odot}}$, but it is surrounded by a massive, dense envelope and hence begins to accrete rapidly. If the envelope does not fragment, and if accretion onto the protostar is unhindered by radiative feedback, then the final mass of the star can be very large. For example, @on07 estimate that its mass could lie anywhere in the range $M \sim 20$–$2000 \: {\rm M_{\odot}}$, depending on the environment and dynamical history of the gas. The degree to which feedback will suppress accretion remains uncertain, but the most effective potential feedback mechanisms, Lyman-$\alpha$ scattering and photoionization [@mt08] become important only once the mass of the star exceeds 20 M$_{\odot}$. The expectation is therefore that all Population III.1 stars were massive, with masses typically of the order of $100 \: {\rm M_{\odot}}$ or more, thereby also explaining why none of these stars appear to have survived until the present day. If a Population III.1 star of this mass does form, then it will rapidly dissociate H$_{2}$ throughout the halo, thereby suppressing further star formation [@on99; @gb01]. This fact, together with the lack of fragmentation seen in the simulations, is often taken to imply that only a single Population III.1 star will form in each dark matter halo.

This widely accepted picture for Population III star formation also provides for a second mode of Pop. III star formation. This occurs in metal-free gas which has been ionized by radiation from a previous generation of Population III stars. Following the death of these stars, the gas recombines, and the elevated fractional ionization in the recombining gas allows more H$_{2}$ to form. The enhanced H$_{2}$ fraction enables the gas to cool to a lower temperature, which in turn increases the effectiveness of cooling by the singly-deuterated hydrogen molecule, HD. In the standard Pop. III.1 scenario, HD cooling is of limited importance, but in this second scenario, termed Population III.2 star formation by @tm08, it becomes dominant, cooling the gas down to the temperature of the cosmic microwave background [@no05; @jb06; @yoh07]. The lower gas temperature, plus the higher critical density of HD ($n_{\rm crit} \sim 10^{6} \: {\rm cm^{-3}}$, compared with $n_{\rm crit} \sim 10^{4} \: 
{\rm cm^{-3}}$ for H$_{2}$), means that fragmentation occurs somewhat later, and produces significantly lower mass fragments, with characteristic masses $M \sim 100 \: {\rm M_{\odot}}$. Subsequently, the evolution of these fragments is believed to proceed in the same fashion as described above, resulting ultimately in the formation of a primordial protostar of similar mass. However, since this protostar is embedded in a less massive envelope, with a lower accretion rate, the final stellar mass is thought to be an order of magnitude or so smaller, $M \sim 10 \: {\rm M_{\odot}}$. Nevertheless, this still corresponds to what we would call, by Galactic standards, a massive star, and moreover one which is more than capable of dissociating H$_{2}$ throughout a large volume of the halo, thereby suppressing further star formation.

In the past few years, it has become possible to partially test this picture using numerical simulations that begin with the proper cosmological initial conditions and follow the collapse of the gas all the way to protostellar densities. These simulations confirm that fragmentation during the initial collapse phase is ineffective: in general, only a single protostar forms, although in roughly 20% of cases, the gas fragments into two clumps, forming a wide binary (@tao09; M. Turk, private communication). However, most of these studies were prevented by technical limitations from following the evolution of the infalling gas after the formation of the first protostar. The technical problem involves the hydrodynamical timestep: as the collapse is followed down to protostellar densities, this becomes extremely short ($\Delta t \sim 10^{-3} \: {\rm yr}$), making it computationally infeasible to follow the evolution of the surrounding dense gas over any significant time period.

In studies of contemporary star formation, a similar problem is avoided by the use of sink particles [@bbp95]. When gravitationally bound regions of gas collapse below the scale on which they can be spatially resolved, they are replaced in the simulation by a sink particle of the same mass. This can accrete gas from its surroundings, but otherwise interacts only via gravity. Recent simulations of Population III star formation that have used sink particles to allow the evolution of the gas to be followed beyond the formation of the first protostar consistently find evidence for more extensive fragmentation than is assumed in the conventional picture of Population III star formation (@cgk08 [@sgb10]; Clark et al., in prep.). If fragmentation is as common as these results suggest, and most Population III stars form in multiple systems rather than as single stars, then this has profound implications for the final masses of the stars, their production of ionizing photons and metals, the rate of high redshift gamma ray bursts, and many other issues.

It is therefore important to better understand the physical basis of fragmentation in these systems. To do this, we would ideally like to have a large sample of simulations, as basing our arguments on only one or two realizations of Population III star formation leaves one open to the possibility that these realizations may not be typical. Unfortunately, simulations starting from cosmological initial conditions and following collapse all the way to protostellar densities are computationally costly, making it difficult to explore a large region of parameter space. In these simulations, it is also difficult to be certain about which aspects of the included physics are the most important for driving fragmentation. Simulations that start from simpler initial conditions and that focus on exploring the importance of a single free parameter therefore play an important role, which is complementary to that of fully cosmological simulations. A good example is the recent work by Machida and collaborators [@mach08; @mach09] which examined the influence of the initial rotational energy of the gas, and found that zero metallicity clouds with sufficient initial rotational energy could fragment into tight binaries. Another example, and one which explores a much lower density regime, is the study by @jap09a [@jap09b], who showed that the fragmentation behavior depends sensitively on the adopted density profile of the primordial halo.

In this paper, we perform a similar study into the effects of the initial turbulent energy of the gas. Our motivation comes from the fact that recent high-resolution cosmological simulations have shown that the self-gravitating regions in which Pop. III star formation occurs contain significant turbulent motions [see e.g. @2008MNRAS.387.1021G]. From studies of present-day star formation, we know that such turbulent motions can lead to the fragmentation of initially marginally unstable gas clouds , and so it is instructive to investigate how such conditions could affect the ‘standard’ picture of primordial star formation. In a Pop. III analogue of the studies on present-day star formation, we will look at how subsonic turbulence affects the collapse of marginally super-critical Bonnor-Ebert spheres, in an attempt to quantify the effects of turbulently induced fragmentation on the mass function of Pop. III stars. Central to this study will be the use of the sink particles to allow us to capture the evolution of the gas cloud beyond the collapse of the first region.

The paper is laid out in the following manner. Our modifications to the cosmological SPH code Gadget 2 [@springel05] are outlined in §\[sec:gadget2\] (with further details of the chemical networks and heating and cooling described in the Appendix). The initial conditions for this numerical experiment are described in §\[sec:ics\], including the definitions of Pop. III.1 and III.2 that are used in our study. The details of the fragmentation process seen in the simulations are discussed in §\[sec:frag\] and the long-term chemical and thermodynamical evolution of the infalling envelope are described in §\[sec:evol\]. We discuss the implications of this study for the initial mass function (IMF) of primordial stars in §\[sec:chat\], and summarize the main points of this paper in §\[sec:sum\].

Numerical method {#sec:gadget2}
================

We model the evolution of the gas in our simulations using a modified version of the Gadget 2 smoothed particle hydrodynamics code [@springel05]. We have modified the publicly available code in several respects. First, we have added a sink particle implementation, based on the prescription in @bbp95, to allow us to follow the evolution of the gas beyond the point at which the first protostar forms. Our particular implementation is derived from the one first described in @jap05. We briefly describe the ideas behind the algorithm here. The actual numerical values used for the parameters discussed here are given later, in § \[sec:ics\].

Sink particles are not really added to the code in the sense that a new particle is introduced, but instead a normal SPH particle is turned into a sink particle once certain criteria have been met. The particle undergoes a series of tests once it has reached a threshold density. The first is to see whether the candidate SPH particle is sufficiently far away from any other sinks, measured in terms of the sink particle’s accretion radius, $r_{\rm acc}$. We adopt a conservative value of $2 r_{\rm acc}$. Next, we check to see whether the smoothing length of the particle is less than the ‘accretion radius’ of the sink particle that it will become. This ensures that when the sink particle forms it can instantly accrete at least $\sim 50$ neighboring particles (we adopt 50 neighbors in these simulations). The third test is to make sure that the candidate sink particle and its neighbors are on the same integration time-step.

Once these three preliminary criteria are met, the dynamical state of the possible sink particle and its neighbors are assessed to ensure that the particles are indeed undergoing gravitational collapse and are not about to re-expand from their dense state. This takes the form of a further four tests. Firstly, we require that

$$\alpha \le \frac{1}{2}$$

where $\alpha$ is the ratio of the thermal energy to the magnitude of the gravitational energy of the particles. Secondly, we ensure that

$$\alpha + \beta \le 1$$

where $\beta$ is the ratio of rotational energy to the magnitude of the gravitational energy. The third condition is that the total energy of the particles must be negative (which actually renders the above checks redundant, but can help to improve the computational efficiency of the code). Finally, the forth test is that the divergence of the accelerations must be less than zero. This final check ensures that the group of particles is not in the process of being tidally disrupted or bouncing. If all of these tests are passed, the particle becomes a sink and the mass and linear momentum of the neighbors are added to the sink particle.

As the simulation progresses, the sink particles are then allowed to accrete other gas particles that fall within the accretion radius. As in the Bate et al. (1995) prescription, several tests must be passed before any SPH particles can be accreted by the sink. First, it must obviously be bound and moving toward the candidate sink. Second, in the case where there is more than one sink present, it must be more bound to the candidate sink than any other sink in the simulation. Finally, the SPH and sink particles need to be on the same integration time-step (to ensure temporal momentum conservation). Once these conditions are passed, the mass and linear momentum of the SPH particle are added to the sink particle. Note that accretion onto existing sink particles is done [*before*]{} any new candidate sinks are considered.

In addition to the sink particles, we have also implemented an external pressure term [e.g. @benz90], that enables us to model a constant pressure boundary, as opposed to the vacuum or periodic boundary conditions that are the only choices available in the standard version of Gadget. We modify the standard gas pressure contribution to the Gadget 2 momentum equation,

$$\frac{d v_{i}}{d t} = - \sum_{j} m_{j} \left[
f_{i}\frac{P_{i}}{\rho_{i}^{2}} \nabla_{i} W_{ij}(h_{i})
+ f_{j}\frac{P_{j}}{\rho_{j}^{2}}\nabla_{i} W_{ij}(h_{j}) \right] ,$$

by replacing $P_{i}$ and $P_{j}$ with $P_{i} -  P_{\rm ext}$ and $P_{j} -  P_{\rm ext}$ respectively, where $P_{\rm ext}$ is the external pressure, and all quantities have the usual meaning, consistent with those used by @springel05. The pair-wise nature of the force summation over the SPH neighbors ensures that $P_{\rm ext}$ cancels for particles that are surrounded by other particles. At the edge, where the term does not disappear, it mimics the pressure contribution from a surrounding medium. The values of $P_{\rm ext}$ used in this study are $3 \times 10^{7} k_{\rm B} \: {\rm K \: cm^{-3}}$ for the Pop III.1 clouds, and $7.5 \times 10^{6} k_{\rm B} \: {\rm K \: cm^{-3}}$ for the Pop III.2 clouds, where $k_{\rm B}$ is the Boltzmann constant. These confining pressures are similar to the internal pressure of the clouds in each case (see the temperatures and densities in §\[sec:ics\] below).

Finally, we have added a treatment of primordial gas chemistry and cooling to the code. To model the thermal evolution of the gas, we use an operator-split formalism, which treats the effects of radiative and chemical heating and cooling separately from compressional heating. The influence of radiative and chemical heating and cooling on the thermal energy of each particle can thus be formulated as an ordinary differential equation (ODE), which can then be solved simultaneously with the chemical rate equations (also ODEs) that describe the chemical evolution of the gas. As @tao09 have previously discussed, the strong coupling between the chemical and thermal evolution of high density primordial gas that results from the importance of three-body H$_{2}$ formation heating and H$_{2}$ collisional dissociation cooling renders it vital to treat cooling and chemistry simultaneously; treatments that do not do so will produce numerically stable results only with great difficulty. Full details of our chemical model and cooling function are given in the Appendix.

![image](fig01a){height="3.0in"} ![image](fig01b){height="3.0in"}

![image](fig01c){height="3.0in"} ![image](fig01d){height="3.0in"}

![\[fig:sinkmf\] The top panel shows the mass functions from those Pop. III.1 simulations in which fragmentation occurs. In all cases the mass function is plotted at the point at which the total mass of gas converted to sink particles is 100. Note that as accretion is ongoing, and the system is still young ($t \sim 1000$ yr), these will often not be the final masses of the sinks. The mass functions in the individual simulations differ substantially, although the combined mass function, shown in the bottom panel, exhibits a broad and flat distribution between masses of 0.4 and 20 .](fig02a){height="2.5in"}

![\[fig:sinkmf\] The top panel shows the mass functions from those Pop. III.1 simulations in which fragmentation occurs. In all cases the mass function is plotted at the point at which the total mass of gas converted to sink particles is 100. Note that as accretion is ongoing, and the system is still young ($t \sim 1000$ yr), these will often not be the final masses of the sinks. The mass functions in the individual simulations differ substantially, although the combined mass function, shown in the bottom panel, exhibits a broad and flat distribution between masses of 0.4 and 20 .](fig02b){height="2.5in"}

![image](fig03a){height="3.4in"} ![image](fig03b){height="3.4in"}

![image](fig04a){height="3.in"} ![image](fig04b){height="3.in"}

![image](fig04c){height="3.in"} ![image](fig04d){height="3.in"}

![\[fig:rhotevol\] Temperature as a function of number density for the Pop. III.1 (dark blue) and Pop. III.2 (light blue) $\Delta v_{\rm turb} = 0.1\,c_{\rm s}$ simulations. In both cases, the curves denote the state of the cloud at the point just before the formation of the sink particle. ](fig05a){height="3.0in"}

![image](fig06a){height="3.in"} ![image](fig06b){height="3.in"}

![\[fig:mdotcloud\] Accretion rates as a function of enclosed gas mass in the Pop. III.1 (upper lines; blue) and Pop. III.2 (lower lines; magenta) simulations, estimated as described in Section \[pop31\]. Note that the sharp decline in the accretion rates for enclosed masses close to the initial cloud mass is an artifact of our problem setup; we would not expect to see this in a realistic Pop. III halo.](fig07a){height="3.0in"}

![\[fig:tsf\] Time taken for the first sink particle to form in each simulation. For reference, the horizontal dashed line denotes the free-fall time of the clouds at their initial density.](fig08a){height="3.in"}

![image](fig09a){height="3.in"} ![image](fig09b){height="3.in"}

![\[fig:sinkmf32big\] (top) Mass functions for our two 1000 Pop.  III.2 clouds. As with the Pop.  III.1 data in Fig. \[fig:sinkmf\], we show the mass function after 100 of gas has been accreted by the sink particles. Due to the reduced fragmentation in these calculations, the sink particles are on average more massive than their Pop.  III.1 counterparts since they have less competition for the available mass. (bottom) The combined mass function for the two simulations.](fig10a){height="2.5in"}

![\[fig:sinkmf32big\] (top) Mass functions for our two 1000 Pop.  III.2 clouds. As with the Pop.  III.1 data in Fig. \[fig:sinkmf\], we show the mass function after 100 of gas has been accreted by the sink particles. Due to the reduced fragmentation in these calculations, the sink particles are on average more massive than their Pop.  III.1 counterparts since they have less competition for the available mass. (bottom) The combined mass function for the two simulations.](fig10b){height="2.5in"}

![image](fig11a){width="2.5in"} ![image](fig11b){width="2.5in"}

![image](fig11c){width="2.5in"} ![image](fig11d){width="2.5in"}

![image](fig12a){width="3.5in"} ![image](fig12b){width="3.5in"}

Initial conditions {#sec:ics}
==================

Since the aim of this study is to investigate the fragmentation properties of primordial gas, we use initial conditions that allow us to perform a controlled numerical experiment. The clouds start as unstable Bonnor-Ebert (BE) spheres into which we inject a subsonic turbulent velocity field. The BE sphere is made by allowing a gravitationally stable cloud of gas to evolve under its own self-gravity until the system has settled into a stable, centrally condensed configuration. During this initial settling phase, the gas temperature is held constant and the chemical evolution is not followed. On the scales of interest in this study, gravitational forces from the dark matter are negligible compared to the self-gravity of the gas, and so for simplicity, we do not include dark matter in our models.

By re-scaling the mass and temperature of the cloud, we are then able to choose initial conditions that are gravitationally unstable, and which are similar to the initial conditions for the Pop. III.1 and Pop. III.2 star formation channels that occur in cosmological simulations. All of the simulations presented in this paper start with a maximum central cloud number density of $10^{5}$ cm$^{-3}$. The Pop.  III.1 simulations start with an initial temperature of 300K and contain 1000 of gas. For the Pop. III.2 simulations, we adopt an initial temperature of 75 K and examine two different initial masses. The first set of Pop. III.2 simulations contain 150 of gas, chosen such that they have the same ratio of thermal to gravitational energy as the Pop. III.1 simulations, roughly 0.15. As such, these Pop. III.2 simulations contain the same initial number of Jeans masses of gas as the Pop. III.1 simulations. For a uniform sphere, the number of Jeans masses is given by ($E_{\rm therm}/|E_{\rm grav}|)^{-3/2}$, so these clouds have roughly 3 Jeans masses in the initial configuration. All else being equal, if these clouds were to evolve isothermally from this point on, they would have the same propensity to fragment as the Pop. III.1 clouds, since this is directly related to the initial ratio of gravitational to thermal energy. Note that we have chosen a slightly subvirial configuration for our initial setup, to ensure that the clouds are still able to collapse when the turbulent motions are included. In addition to these simulations, we also performed a second set of Pop. III.2 simulations that start with the same gas mass as in the Pop. III.1 case (1000 M$_{\odot}$). In this case, the Pop.  III.2 simulations are initially more Jeans unstable than their Pop.  III.1 counterparts, and therefore might be expected to fragment significantly more.

In our Pop. III.1 simulations, we set the initial fractional abundances of H$_{2}$, H$^{+}$, HD and D$^{+}$ to $x_{\rm H_{2}} = 10^{-3}$, $x_{\rm H^{+}} = 10^{-7}$, $x_{\rm HD} = 3 \times 10^{-7}$ and $x_{\rm D^{+}} = 2.6 \times 10^{-12}$, respectively. Our values for $x_{\rm H_{2}}$, $x_{\rm H^{+}}$ and $x_{\rm HD}$ are typical of the values found at these densities in cosmological simulations of Pop. III.1 star formation [see e.g. @2008MNRAS.387.1021G], and account for the fact that the HD/H$_{2}$ ratio is elevated over the cosmological D/H ratio of $2.6 \times 10^{-5}$ [@mol08] owing to the effects of chemical fractionation [@g08]. In the case of D$^{+}$, fractionation is unimportant at our starting temperature, and so we simply set ${\rm D^{+}} / {\rm H^{+}} 
= 2.6 \times 10^{-5}$. In our Pop. III.2 simulations, we adopt the same initial H$^{+}$ and D$^{+}$ abundances, but set $x_{\rm H_{2}} = 3 \times 10^{-3}$ and $x_{\rm HD} = 3 \times 10^{-6}$, following @2008MNRAS.387.1021G. In both the Pop. III.1 and Pop. III.2 simulations, we assumed that all of the helium remained neutral, and set the initial abundances of all of our other tracked species to zero.

Within the BE spheres, we impose a turbulent velocity field that has a power spectrum of $P(k) \propto k^{-4}$. We assume that the turbulence considered here has its origin in gravitationally driven flows that arise as the gas and dark matter virialize in mini-halos (Wise & Abel 2007; Greif et al. 2008; Klessen & Hennebelle 2010). As the gas is compressible in nature the turbulent velocity field will have a power spectrum that is somewhat steeper than the standard Kolmogorov (1941) description for incompressible flows. However, this deviation is small and we note that the the ability of a cloud to fragment is only weakly dependent on the power spectrum of the turbulence [@ddetal04]. The three-dimensional root-mean-squared velocity in the turbulent field – which we will refer to as $\Delta v_{\rm turb}$ – is then scaled to some fraction of the sound speed in the initial conditions. For the simulations presented here we use four different rms velocities: 0.1, 0.2, 0.4 and 0.8 $c_{\rm s}$. For an isothermal sound speed and an adiabatic index of $\gamma = 5/3$, the corresponding ratios of the turbulent to thermal energy are given by $1/3(\Delta v_{\rm turb}/c_{\rm s})^{2}$, yielding 0.0033, 0.0133, 0.0533 and 0.2133 respectively, for our chosen values of $c_{\rm s}$. In order to focus on the effects of the turbulence, we do not include any ordered rotation of the initial gas cloud. Note, however, that this does not imply that the initial angular momentum of the cloud is zero, since the imposed turbulent velocity field gives the cloud a small amount of angular momentum. Note that we only consider subsonic turbulence in this study since our clouds have only a few Jeans masses, and supersonic turbulence would unbind them. To study the effects of supersonic turbulence, one would have to look at clouds that are initially more Jeans unstable than those we study here.

The clouds in our study are all modeled using 2,000,000 SPH particles. Although this means that the mass resolution is higher in the 150 Pop.  III.2 clouds than in the other simulations, the Jeans mass at the point where sink particles form is well resolved in every case (see below). For this study, it is more important that the turbulence in all simulations is evolved with the same resolution, hence our choice of a constant particle number throughout.

Sink particles are created once the number density of the gas reaches $10^{13}$ cm$^{-3}$, at which point the gas has a temperature of around 1200K. The corresponding Jeans mass at this density and temperature is 0.08 . Our mass resolution in the Pop.  III.1 and the 1000 Pop.  III.2 clouds is $2 N_{\rm neigh} m_{\rm part} = 0.05$ , where $N_{\rm neigh}$ is the number of neighbors employed for force evaluations (in our case 50), and $m_{\rm part}$ is the mass of an SPH particle. The 150 Pop. III.2 simulations have a mass resolution of 0.008 . Once the candidate particle has passed the criteria described in Bate et al. (1995), it is replaced by a sink particle that can accrete gas particles that fall within its accretion radius $r_{\rm acc}$, which we fix at 20 AU. Note that this radius is significantly larger than the Jeans radius at the point that the gas reaches the density threshold for sink creation, which is around 6 AU. As such, all the fragmentation that we capture in this study is well resolved and not lying close to the limits of our resolution; if any ‘artificial’ fragmentation were to occur, the sink particle would immediately swallow the offending region and replace it with a single accreting point. We also prevent sink particles from forming within $2\, r_{\rm acc}$ of one another. This prevents the formation of sinks out of gas that in reality would by accreted by a neighboring sink particle before it could go into direct collapse by itself. Lastly, gravitational interactions between sinks, and between the sinks and the gas, are softened in the standard way in Gadget 2, using a fixed softening parameter of 5 AU for the sinks, and a variable softening parameter for the gas particles that is equivalent to their smoothing length.

We note that the sink algorithm employed here does not allow sink particles to merge, in contrast to the ‘sticky’ sinks introduced in Bromm et al. (2002). As such, our simulations may be biased towards low masses. Small, secondary sinks could be driven towards more massive ones via dynamical friction and coalesce with them. This effect, which acts to reduce the resulting number of fragments, has been seen in AMR simulations of massive star forming regions in the present-day universe (Krumholz et al. 2007, 2009). However, given the size of our sink particles, it is unclear whether sink-merging is really applicable: the protostellar radius is at most expected to be 0.5 au (and then only for a short time, e.g. @ho2009), which is significantly smaller than both our adopted accretion radius and the gravitational softening. In fact, in the simulations from present-day star formation in which extremely small sink particles are used – along with little, or no softening — merging events are found to be fairly rare (e.g. @Bate2009). They do however occur, and this should be borne in mind when interpreting the results of this study. Also, given the size of our sink particles, we are unable to resolve any of the disks that would invariably form around the young protostars and it has been suggested that such structures would be unstable to fragmentation on these scales (e.g. @momi2008). Finally, our calculations do not include any model for the feedback processes that accompany the formation of a young star (e.g. @mt08). Given these uncertainties, the sink particles are unable to say exactly what the shape of the final IMF will be, but rather measure how the gas can fragment at a given scale (our resolution), and how these fragments are likely to evolve, assuming feedback processes play only a minor role over the timescales investigated in this study. These caveats should be borne in mind when interpreting the results in the following sections.

The fragmentation of primordial gas {#sec:frag}
===================================

In our comparison of the different simulations performed in this study, we will compare the properties of the clouds after roughly 10 percent of their mass has been accreted (although in some cases we will mention in passing what happens as the simulations are advanced further). In terms of looking at the ability of the clouds to fragment, comparing the different simulations at this point in their evolution ensures that they have turned the same faction of their initial number of Jeans masses into the sink particles. The exception in our analysis is the comparison between the 1000 Pop. III.1 and III.2 simulations, since in this case the Pop. III.2 clouds are initially more Jeans unstable. We now go on to describe the fragmentation in the Pop. III.1 and III.2 channels in some detail.

Pop. III.1 clouds {#pop31}
-----------------

The panels in Fig. \[fig:image3.1\] show the column density distribution in the Pop. III.1 simulations after 10 percent of their mass has been accreted onto the sink particles. The images show the inner 1300 AU of the cloud, centred on the first sink particle to form in each simulation. We can see that the clouds with $\Delta v_{\rm turb} \geq 0.2 c_{\rm s}$ fragment into small clusters of sink particles; the runs with $\Delta v_{\rm turb} = 0.2 c_{\rm s}$, $0.4 c_{\rm s}$ and $0.8 c_{\rm s}$ form 5, 31 and 15 sink particles respectively. These calculations clearly demonstrate that turbulent subsonic motions are able to promote fragmentation in primordial gas clouds. Only in the case with $\Delta v_{\rm turb} = 0.1 c_{\rm s}$ does the cloud form just a single sink particle. In this case an extended disk builds up around the star, since the seed turbulence gives rise to some low level of rotation in the collapsing core, but otherwise the gas contains little structure.

Interestingly, there is no clear trend linking the number of fragments that form and the initial turbulent energy: the 0.4 run fragments more than the 0.8 run, despite containing only one quarter of the initial turbulent energy. The effects at play here are somewhat complex. First, the turbulent velocity field contained in the collapsing region will differ with each value of $\Delta v_{\rm turb}/c_{\rm s}$, since the different strengths of flow will push the gas around to different degrees. In addition, the nature of the turbulence that survives in the collapsing core will also affect the fragmentation. As we see from the images in Fig. \[fig:image3.1\], the cloud with $\Delta v_{\rm turb} = 0.4 c_{\rm s}$ – the most successful in terms of fragmentation – forms a large disk-like structure. As we will discuss in §\[sec:evol\], this configuration appears to aid the fragmentation of the infalling envelope. Therefore, much of the cloud’s ability to fragment depends on the level of rotation that happens to become locked-up in the collapsing region. Further, and to a lesser degree, the extra delay in the collapse caused by the increased turbulent support also gives the cloud more time to wash out anisotropies in the gas (see Fig. \[fig:tsf\] for the collapse times). Thus the ability of the cloud to fragment is a competition between these conflicting processes. For the randomly generated velocity field that we used in this study, the ability to fragment is better for $\Delta v_{\rm turb}/c_{\rm s} = 0.4$, than $\Delta v_{\rm turb}/c_{\rm s} = 0.8$, but we note that this may not always be the case. We stress that to make a quantitative statistical statement about the number of fragments that form as a function of the turbulent Mach number, we would need to run a series of different realizations of the turbulent velocity field in each case. Such a comparison lies outside the scope of our current study.

The mass functions of the sink particles from the simulations that undergo fragmentation are shown in Fig. \[fig:sinkmf\]. For clarity, we have omitted the single 100 sink particle that forms in the 0.1 cloud. For the 0.2 and 0.8 clouds, we see that the sink masses cluster around some central value – roughly 12 and 4 respectively – while the 0.4 cloud has a mass function that is skewed to lower masses, with a peak at around 1 , a sharp fall-off below this, and a broad distribution towards higher masses, extending up to around 13 . We emphasize that the mass functions presented here (as well as in Fig. 10) do not yet represent the final IMFs, as they correspond to an intermediate time in the overall accretion process, where only 10% of the cloud has been accreted. They also might be affected by the numerical details of our sink technique, in particular the absence of any sink-sink mergers, as discussed in Section 3.

The reason for the spread of masses becomes apparent when we look at the accretion properties of the sink particles, which are shown in Fig. \[fig:maccp1\]. Focusing on the 0.4 and 0.8 runs, and looking at the evolution of the individual sinks in more detail, we see that some appear to accrete rapidly and then suddenly stop. This behaviour is typical of what is seen in simulations of bound, fragmenting cores in the context of present-day star formation [@kb00; @Klessen01; @kb01; @bvb04; @sk04], and is a result of velocity kicks from dynamical three (or more) body interactions. The mass accretion rate in such a system [@bonnell01a] is given by $$\label{equ:bh}
\dot{m_{*}} \propto \rho \frac{m_*^2}{v_{\rm rel}^3},$$

where $m_{*}$ is the mass of the sink particle, $\rho$ is the gas density, and $v_{\rm rel}$ is the velocity of the sink particle relative to the gas. The ability of a sink to accrete more mass from the available reservoir is significantly reduced once its velocity increases. The effect is then exacerbated by the fact that an increase in velocity results in the sink particle moving to a more distant orbit (or even being kicked out of the system entirely) and hence into a position where the gas density is lower. Since this sink particle is now in a position where further accretion is difficult, its siblings are able to accrete its ‘share’ of the mass reservoir, with the majority going to those few sinks that sit right in the middle of the cloud’s potential well. In general, as the sinks accrete from the background gas, they tend to move towards the centre, due to mass-loading. Further, their increased mass makes them more likely to survive dynamical encounters with their less massive siblings. As such, the ‘rich get richer’, with a few sinks ending up significantly more massive than the rest. The process is typically termed ‘competitive accretion’, [@bonnell97; @bonnell01a; @bonnell01b; @bb06] and normally leads to the type of distribution of masses seen in our 0.4run (Fig. \[fig:sinkmf\]).

Looking at the mass evolution of the individual sink particles, we also see that the formation of new sink particles occurs in bursts. The turbulence generates structure in the gas which is enhanced by the gravitational collapse as the gaseous envelope falls in towards the central system of sink particles. The bursts in sink formation reflect the moments when these structures detach from the flow and become self-gravitating in their own right.

Figure \[fig:maccp1\] also shows the accretion rate of the cluster as a whole, and how that compares to an estimate made from the radial infall profile. To construct this estimate, we first compute the mass infall rate as a function of the radial distance $r$ from the densest SPH particle: $$\label{equ:mdotinfall}
\dot{m}(r) = 4 \pi \,r^{2}\,\rho(r)\, v_{r}(r),$$

where $\rho(r)$ is the gas density in a spherical shell with radius $r$ and width d$r$, and $v_{r}(r)$ is the radial velocity of the gas in this shell and all quantities are volume-averaged. Given the enclosed mass as a function of radial distance, $m_{\rm enc}(r)$, it is straightforward to convert from $\dot{m}(r)$ to an infall rate as a function of enclosed mass, $\dot{m}(m_{\rm enc})$, which we identify with the total mass accretion rate of the system of sink particles. Figure \[fig:maccp1\] demonstrates that this estimate provides a fairly close match to the actual accretion rate onto the sink particle population, except at very early times, which is a numerical artifact: the sink particles instantly accrete all gas within their accretion radius when they form. At later times, the main difference between the estimated accretion rate and the true accretion rate is that the latter is significantly noisier, owing to the bursts of sink formation and the clumpiness of the infalling gas. Finally, we note that the accretion rates are also fairly insensitive to the level of turbulence in the cloud, suggesting that random turbulent motions with the magnitudes considered here should not lead to a variation in the overall accretion rates of primordial stars or star clusters from minihalo to minihalo. Other sources of support against gravity, such as rotation and magnetic fields (assuming the latter can be efficiently generated, as in @tb04 and @schl10), are likely to play a greater role in regulating the accretion rate. Indeed, the important role played by rotation in regulating the infall and accretion of gas can be appreciated if we compare the accretion rates measured in our simulations with those estimated or measured in previous studies of Pop. III star formation starting from more realistic cosmological initial conditions [e.g. @abn02; @brlb04]. The absence of initial rotational support in our simulations leads to higher infall velocities, and hence to an accretion rate that is a factor of a few larger than these previous values.

Pop.  III.2 clouds
------------------

### Small clouds with 150 

The first clouds we will examine in the Pop. III.2 case are those with an initial energy balance similar to those studied in the Pop. III.1 case: that is, clouds that have only a few Jeans masses initially. Since the Pop. III.2 channel is cooler at number densities around $10^{5}$ cm$^{-3}$, with a typical temperature of around 75K, the same initial number of Jeans masses requires that the clouds have a lower mass of 150 .

The column density distribution in these clouds after 10 percent of the gas has been accreted by the sink particles is shown in the column density images in Fig. \[fig:image3.2\]. We see from the images that these clouds undergo significantly less fragmentation than their Pop. III.1 counterparts, forming at most 3 sink particles after 10 percent of the cloud mass has been accreted. While these calculations demonstrate that this primordial star formation channel is susceptible to fragmentation if turbulence is present in the collapsing gas, it appears to be significantly more stable than the Pop. III.1 channel that gives rise to the first stars in the universe.

The main reason why this mode of star formation is less susceptible to fragmentation has to do with the thermal evolution of the gas as it collapses. In Fig. \[fig:rhotevol\] we show the temperature of the gas as a function of its number density, for the Pop. III.1 and Pop. III.2 cases, taken from the simulations with $\Delta v_{\rm turb} = 0.1\: {c_{\rm s}}$. The elevated ${{\rm H_{2}}}$ and HD fractions in the Pop. III.2 simulations do not provide enough cooling to keep the gas close to the CMB temperature at these densities, owing to the increasing inefficiency of the HD as a coolant as it nears the critical density at which the populations of its rotational level reach their local thermodynamic equilibrium (LTE) values. The gas therefore heats up significantly as it collapses, increasing its temperature from 75 K at $n = 10^{5} \: {\rm cm^{-3}}$ to roughly 700 K at $n = 10^{7} \: {\rm cm^{-3}}$, corresponding to evolution that is almost adiabatic. This sharp increase in temperature temporarily increases the local Jeans mass at the center of the cloud, and significantly slows the collapse, giving the turbulence in the cloud time to decay. In the absence of additional physical processes able to replenish the turbulence, there is nothing to sustain the density inhomogeneities in the cloud that act as the seeds for later fragmentation (since they are not yet self-gravitating), and so the end result is a collapse with a much lower level of fragmentation than in most of the Pop. III.1 simulations. A similar effect has previously been noted by @to08 in their study of fragmentation in very metal poor gas clouds. They find that for metallicities of around ${\rm Z} \sim 10^{-4.5} \:
{\rm Z_{\odot}}$, heat input due to three-body H$_{2}$ formation at $n \sim 10^{8} \: {\rm cm^{-3}}$ leads to a sharp jump in the gas temperature, which delays the collapse, reduces the elongation of the collapsing core, and suppresses any fragmentation. Furthermore, @yoh07 also briefly addressed this issue in their study of Pop. III.2 star formation, and showed that their simulated Pop. III.2 prestellar core would be stable against gravitational deformation at similar densities, owing to its hard effective equation of state.

Although the evolution of the temperature with density is significantly different between the Pop. III.1 and Pop. III.2 channels, we see that the accretion rates are similar when we consider the inner 0.01 to 10 solar masses of the collapsing envelope (Fig. \[fig:maccp2\]), consistent with the results from @yoh07. However, if we consider the evolution of the accretion rates over the whole cloud, we see that at later times the two channels depart significantly from one another, with accretion occurring at a significantly slower rate in the Pop. III.2 case (Fig. \[fig:mdotcloud\]). Also, at very early times, the accretion rates in the Pop. III.2 clouds are more sensitive to the level of turbulence than the Pop. III.1 clouds, and the turbulence has a somewhat stronger effect in delaying the onset of star formation in these calculations (Fig. \[fig:tsf\]).

### Large clouds with 1000 

Given that the Pop. III.2 clouds seem to be much more stable against fragmentation, it is worthwhile investigating whether they can be made to fragment when the gas is initially more Jeans unstable. Recent simulations of the formation of the first galaxies show that regions where Pop. III.2 star formation occurs are fed by cold, supersonic turbulent streams of gas [e.g. @2008MNRAS.387.1021G; @2008ApJ...685...40W]. As such, the initial condition for the Pop. III.2 channel in this picture may have significantly more than one Jeans mass, due to the rapid assembly of the self-gravitating core. In this section we consider Pop. III.2 clouds that contain 1000 , but otherwise have the same properties as the 150 clouds (i.e., same initial temperature, density, and initial chemical composition). Since these clouds are colder than the Pop. III.1 clouds of the same mass and density, they are initially more Jeans unstable, having around 26 Jeans masses at the start of the simulation.

The results of two such calculations are shown in Fig. \[fig:image3.2big\], in which the initial levels of turbulence have been set to 0.4 and 0.8 $c_{\rm s}$. Again we show the evolution at the point where 10 percent of the cloud’s mass has been converted into (or accreted onto) sink particles. Although these clouds are now initially more unstable than the Pop. III.1 clouds, they form significantly fewer sink particles, only 8 in the 0.4 cloud and 7 in the 0.8 cloud (see Fig. \[fig:sinkmf32big\]). The fact that the mass functions of the sink particles appear to have two peaks is due to the cloud undergoing two distinct bursts of sink formation.

The interesting result here is that these clouds exhibit [*less*]{} fragmentation than their Pop. III.1 counterparts, despite containing [*more*]{} Jeans masses in the initial configuration. The main reason for this is the gas in the Pop. III.2 clouds follows an extremely ‘stiff’ effective equation of state (essentially adiabatic evolution), as can be seen from the temperature-density relationship shown in Fig. \[fig:rhotevol\]. Such a rapid increase in the temperature during the initial collapse makes generating structure in the gas difficult, helping to remove any anisotropies introduced by the turbulent flows. An additional effect is that for equations of state with effective adiabatic index of $\gamma \ga 4/3$ (as is the case with these clouds), the Jeans mass increases with increasing density. As such, the collapse halts until sufficient mass has been assembled, and the new Jeans mass has been reached. Since the cooling time is significantly longer than the free-fall time, the gas has a chance to remove structure that could potentially assist the fragmentation at higher densities, where the effective $\gamma$ is more similar to that found in the Pop. III.1 clouds. The combination of these effects results in a gas which is much more stable to fragmentation during collapse than in the Pop. III.1 simulations.

Long term evolution of the infalling envelope {#sec:evol}
=============================================

More insight into the stability of the clouds against fragmentation can be gained from Fig. \[fig:evol\], where we show the temperature and ${{\rm H_{2}}}$ fraction as a function of density, for the 1000 Pop. III.1 and III.2 clouds in which the initial level of turbulence was 0.4 $c_{\rm s}$. The different colors correspond to different points in the evolution of the clouds, with red, green and blue corresponding to the formation of the first sink particle, 50 accreted, and 100 accreted, respectively. The temperature-density plots show the same behavior as reported in @sgb10. The particle evolution appears to diverge at a density of around $10^{10}$cm$^{-3}$, with one group of SPH particles heating up to a maximum temperature of around 7000K and the other staying around 1500K and cooling significantly at higher densities. As discussed in @sgb10, the hot part of the diagram corresponds to gas that falls in at a later stage in the evolution of the cloud. At that point, the enclosed gas mass is larger than at earlier times, and so the free-fall velocity is correspondingly larger. The gas therefore shocks more strongly than at early times, causing it to become hot enough to collisionally dissociate its H$_{2}$ rapidly. With the H$_{2}$ gone, there is nothing to cool the gas until its temperature reaches $T \sim 7000 \: {\rm K}$. At this point, Lyman-$\alpha$ cooling becomes effective, allowing it to resist further heating. Although we see the same trends in both the Pop. III.1 and Pop. III.2 cases, in the latter the amount of gas departing from the ‘standard’ temperature-density evolution is significantly reduced, since the elevated ${{\rm H_{2}}}$ fractions allow the post-shock gas to cool more effectively, limiting the temperature rise and allowing more of the gas to retain its H$_{2}$.

From the accompanying graphs we see that the ${{\rm H_{2}}}$ is rapidly dissociated as the temperature rises, but given the relatively low amounts of ${{\rm H_{2}}}$ present in the gas at these densities, the cooling provided by the dissociation is clearly unable to offset the compression. However we see that at higher densities the ${{\rm H_{2}}}$ suddenly reforms, coinciding with the region in the temperature-density diagrams where the gas is cold. This prompts the question of whether the cold gas is the result of an isochoric cooling instability, brought on by rapid ${{\rm H_{2}}}$ formation? Does the turbulence trigger an instability that naturally exists in pure primordial gas, but that is less effective in gas that has been influenced by the presence of previous star formation?

By harnessing the Lagrangian nature of SPH, the plots in Fig. \[fig:ptrace\] shed some light on this issue. They show the temporal evolution of the density, temperature and ${{\rm H_{2}}}$ fraction for several SPH particles that eventually become sink particles, and as such trace the conditions in the gas in the run-up to gravitational fragmentation and collapse. Again the simulations are those shown in Fig. \[fig:evol\]. The quantities are calculated by averaging over each particle’s 50 nearest neighbors at each instant in time. In the case of the Pop. III.1 cloud, these sinks are the last to form in the simulation, while in the Pop. III.2 cloud, the five lines represent all of the sink formation that occurs between the formation of the first and last sinks.

The figure shows a number of interesting features. First, we see that none of those SPH particles destined to become sinks undergo the rapid rise in temperature – and accompanying loss of ${{\rm H_{2}}}$ – that is shown in Fig. \[fig:evol\]. In contrast, their temperatures remain close to the temperatures found within the first collapsing core (see Fig. \[fig:rhotevol\]). This demonstrates that the cool particles seen in Fig. \[fig:evol\] do not come from regions that undergo shock heating and subsequent loss of ${{\rm H_{2}}}$. Instead, we see that they come from regions of gas that first undergo a relatively quiescent collapse, before being involved in several expansions and contractions.

The fact that this first stage of the density evolution is fairly slow needs to be stressed: the free-fall times at densities of $n = 10^{9} \: {\rm cm^{-3}}$ and $n = 10^{10} \: {\rm cm^{-3}}$ are approximately 1600 yr and 500 yr, respectively. As such, these particles are not experiencing as much compression as those that end up losing their ${{\rm H_{2}}}$. In fact we see that they actually have a higher than average ${{\rm H_{2}}}$ fraction, when we compare them with the lower right-hand plot in Fig. \[fig:evol\], a property that helps them remain fairly cool as they collapse. The reason why they experience less compression is that they are collapsing in a rotating structure that has been formed during the collapse of the turbulent gas.

The evolution of these particles and their immediate surroundings alters abruptly once they enter the dynamically complicated swirling regions that we can see in the column density figures. First, they collide with other material, which results in a sharp increase in their temperature and density. This increase in the density in turn increases the rate of ${{\rm H_{2}}}$ formation, which very rapidly turns them fully molecular. As they re-expand, the adiabatic cooling and the now significantly enhanced ${{\rm H_{2}}}$ line-cooling act together to reduce the temperature. In some cases this happens several times, but for all of these particles, the end result is the same and they find themselves in a clump of gas that is now Jeans unstable and fully molecular: the perfect conditions for forming a new protostar. Interestingly, we see similar behaviour leading up to the formation of the sink particles in both the Pop. III.1 and III.2 clouds. The relative lack of fragmentation in the Pop. III.2 case simply results from the lack of structure in the collapsing envelope, rather than any thermal properties of the gas at these high densities.

Discussion {#sec:chat}
==========

The calculations presented in this study suggest that the first stars in the universe may form in small dense clusters, provided that the turbulent initial conditions we adopt are close to those found in minihalos. This suggests that the Population III IMF covered a broad range in masses, possibly exhibiting a scale-free, power-law extension similar to the present-day case. Previous arguments in favor of a peaked IMF, in the shape of a narrow Gaussian or even a delta function, would then need to be revisited (e.g., Bromm & Larson 2004). Although a similar prediction can be drawn from the study of Clark et al. (2008, hereafter CGK08), there are some subtle differences. Our present study employs a fully self-consistent treatment of the thermodynamics, rather than the piece-wise polytropic equation of state (EOS) approach used by CGK08, which was taken as a fit to the detailed one-zone calculations of @om05. Although both start with similar (Pop. III.1) initial conditions, the fragmentation seen in the current simulations has a different origin. In the self-consistent treatment employed here, the fragmentation is driven by the complicated thermodynamics of high density clumpy gas as it enters the disk-like regions that surround the first protostar. In contrast, the fragmentation in CGK08 was evident at much lower densities, with the turbulent flows forming structures that were enhanced during the collapse, and rotation providing a ‘window of opportunity’ for those structures to become gravitationally unstable in their own right, rather than being simply accreted onto the central protostellar core. Note however that the simulations in CGK08 also contained systematic, solid-body rotation in the initial conditions, which we do not study in this paper. This may help some of the structure to survive.

In general, the calculations presented in this paper suggest that the structure formed early on in the collapse is less likely to survive when the full thermodynamical behavior of the gas is taken into account. This is evident in the fact that the current calculations yield significantly fewer fragments than the Pop III simulations in CGK08 – 25 sink particles for 19 of accreted gas in CGK08, compared to 31 sink particles for 100 of accreted gas in the current study. One feature that should be stressed (and which was originally pointed out in @sgb10) is that the long-term thermodynamic evolution of the envelope differs significantly from the one-zone models. This brings into question the practice of using a piece-wise barotropic EOS as a proxy for the full thermodynamics in simulations that intend to study the evolution of the gas beyond the collapse of the first protostellar core.

Assuming that the turbulent initial conditions used in this paper are indeed representative of the gas in minihalos, it is worth considering the implications of the fragmentation that we see. The mass spectrum of the fragments in our calculations ranges from a few 0.1 to a few 10 . We would not expect stellar feedback to change this result significantly, since previous studies [e.g. @mt08] have shown that feedback becomes effective at limiting accretion onto Pop. III protostars only for protostellar masses greater than about 20 . However, such feedback effects are expected to become important, and in particular to suppress further fragmentation (e.g., Krumholz et al. 2009) during the further evolution of the cloud. The fragments will then likely grow in mass, and possibly even merge with each other. It is very difficult to extrapolate to the final situation where accretion and merging stops. But our main result that the Population III IMF was likely broad seems robust.

As we have discussed above, it is difficult to reach any definitive conclusions regarding the final Population III masses and the resulting IMF. It appears at least possible, however, that in rare cases truly metal-free stars with masses less than $\sim 0.8$ could have formed, and would still be present in the Milky Way today, thus providing a unique opportunity to directly probe the physical conditions at the end of the dark ages with Galactic observations. Current models of hierarchical galaxy formation predict that the first and most metal-poor halos to merge will become part of the bulge component of the resulting spiral galaxy [e.g. @tum10]. It is therefore highly interesting to survey the bulge of the Milky Way for extremely metal-poor and metal-free stars. However, this is also very challenging as the bulge is far away, has high stellar density, and contains a wide range of stellar populations of all ages and metallicities. The odds of finding a few truly metal-free stars amongst millions of other stars are very low, or even zero due to pollution from the interstellar medium (Frebel et al. 2009).

With the same caveat as above, it is interesting to consider the implications of a Pop. III IMF with a characteristic mass significantly smaller than the canonical $100 \: {\rm M_{\odot}}$. For instance, the flux of ionizing photons produced by a population of Pop. III stars has a significant dependence on the form of the Pop. III IMF, as does the metal-enrichment pattern produced by a collection of Pop. III supernovae [@tvs04]. Indeed, the abundance patterns observed to date in extremely metal-deficient stars in the Galactic halo [see e.g. the review by @bc05] are far more consistent with an IMF that produces primarily core-collapse supernovae, with progenitor masses of 10–$40 \: {\rm M_{\odot}}$, rather than with an IMF that produces only very massive pair-instability supernovae, or PISNe [@jog10]. On the other hand, there may be subtle selection effects at work that bias current surveys against finding PISN-enriched stars (Karlsson et al. 2008). The basic argument here is that PISNe have such high metal yields that abundances in stars that form out of this material are already quite high (see Greif et al. 2010), and would therefore be missed in searches that target the lowest metallicities. The interpretation of the large carbon enhancements seen in the population of carbon-enhanced metal-poor stars as the result of the enrichment of these stars by winds from binary companions that have passed through the AGB phase also implies the existence of a large number of intermediate-mass Pop. III stars, with masses $M = 1$–$8 \: {\rm M_{\odot}}$ [@tum07a; @tum07b]. Again, there are alternative models to explain the carbon enhancement in metal-poor stars in terms of nucleosynthesis in faint supernovae (e.g., Iwamoto et al. 2005), in line with a higher characteristic Pop. III mass.

It is important to note that one factor that may limit the impact of turbulent fragmentation on the primordial IMF is if most Pop. III stars form in conditions resembling our Pop. III.2 clouds. It is relatively straightforward to show that most Pop. III stars will form in halos that have been affected in some fashion by a previous episode of Pop. III star formation [@ts09; @ggbk10]. Thus, using the @tm08 terminology, most Pop. III stars will be Pop. III.2 stars. However, it is less obvious how many of these stars will have formed out of gas that has been cooled to temperatures $T \ll 200 \: {\rm K}$ by HD. The effectiveness of HD cooling in these systems depends on the balance between the enhanced formation of H$_{2}$ and HD, owing to the enhanced initial fractional ionization in this gas, and the destruction of H$_{2}$ and HD due to Lyman-Werner band absorption of ultraviolet photons from an extragalactic background (Haiman, Abel & Rees 2000; Johnson, Greif & Bromm 2008), or from local sources [@on99; @gb01]. The relative importance of these effects has not been explored in great detail, and in any case, the outcome is likely to be sensitive to uncertainties in the microphysics that have only recently been resolved [@gsj06; @ga08; @kr10].

To sum up, our calculations suggest that turbulent fragmentation may play an important role in the formation of Pop. III stars, and may strongly influence the form of the Population III IMF. However, the approximate nature of our calculation – specifically, our simplified choice of initial conditions – means that we should regard this at present as no more than a plausible hypothesis. To better establish the role of turbulent fragmentation in real Pop. III minihalos, we will need to be able to simulate a representative sample of both Pop. III.1 and Pop. III.2 minihalos with sufficiently high resolution, such that they can resolve the turbulent flows in the gas at the moment when it becomes self-gravitating. This is a challenging prospect, lying far beyond of the scope of this preliminary study, but our results suggest that it would prove extremely worthwhile.

Summary {#sec:sum}
=======

We have explored the effects of subsonic turbulence on the gravitational collapse of primordial gas clouds. The study employed sink particles to model the run-away collapse of protostellar cores, which allowed us to follow the evolution of the collapsing clouds beyond the formation of the first protostar. The calculations also used a full time-dependent chemical network that accounts for the thermodynamic behaviour of the gas. The current calculations contain neither magnetic fields nor feedback from the protostars. Our main findings can be summarized as follows:

-   Turbulent primordial gas is unstable to fragmentation when one considers the evolution [*beyond*]{} the formation of the first protostellar core.

-   Gas starting from conditions appropriate to Pop. III.1 collapse, rather than Pop. III.2 star formation, is more susceptible to fragmentation. As suggested by @yoh07, the thermal evolution of Pop. III.2 gas as it collapses helps to suppress further gravitational instability over and above the main collapse mode. However some fragmentation in the Pop. III.2 case is seen in this study, caused by the inhomogeneities introduced by the turbulence.

-   In the cases where fragmentation is efficient (in particular the Pop III.1 cloud with turbulent rms velocity $\Delta v_{\rm turb} = 0.4 c_{\rm s}$), the masses of the fragments extend over a large range, which results in a distribution of stellar masses exhibiting a power-law extension towards high mass, as seen in the present-day IMF [@kroupa02; @chabrier03]. However it should be noted that the exact form of the sink particle mass function may depend on our assumed initial conditions (see Appendix \[ICdiscuss\]) and the simplifications employed in the sink particle implementation itself.

-   Due to the relative lack of fragmentation in the Pop. III.2 clouds compared to the Pop. III.1 clouds, the stars formed in the Pop. III.2 calculations are of higher mass on average (for the same accreted mass) than their Pop. III.1 counterparts, even though the mass accretion rate of the star cluster as a whole is higher in the Pop. III.1 clouds.

-   Fragmentation tends to occur in gas which has temperatures of around 200 – 400 K at densities above $10^{11} \: {\rm cm}^{-3}$, significantly lower than the 1000 – 1500 K associated with the collapse of the first core. Such cold temperatures are a result of the expansion that occurs as gas enters the rotating, disk-like regions around the central core, coupled with relatively high fractions of ${{\rm H_{2}}}$, which can provide efficient line-cooling. These rotating structures are themselves a consequence of the angular momentum that is present in the turbulence.

In summary, we propose that if even small levels of turbulence, with velocity dispersions of order 20% of the sound speed or more, are present in the baryonic component of dark matter minihalos, primordial stars are likely to be born in small stellar groups, rather than in isolation, and to have a wide range of stellar masses. Further, the very first stars (Pop. III.1) may have lower masses on average than the second generation of stars (Pop. III.2), contrary to what has previously been assumed.

The authors would like to thank Tom Abel and Thomas Greif for stimulating discussions that helped shape this paper. The work presented in this paper was assisted by the European Commission FP6 Marie Curie RTN CONSTELLATION (MRTN-CT-2006-035890). P.C.C. acknowledges support by the [*Deutsche Forschungsgemeinschaft*]{} (DFG) under grant KL 1358/5. R.S.K. acknowledges financial support from the [*Landesstiftung Baden-Würrtemberg*]{} via their program International Collaboration II (grant P-LS-SPII/18) and from the German [*Bundesministerium für Bildung und Forschung*]{} via the ASTRONET project STAR FORMAT (grant 05A09VHA). R.S.K. furthermore acknowledges subsidies from the DFG under grants no. KL1358/1, KL1358/4, KL1358/5, KL1358/10, and KL1358/11, as well as from a Frontier grant of Heidelberg University sponsored by the German Excellence Initiative. R.S.K. also thanks the Kavli Institute for Particle Astrophysics and Cosmology at Stanford University and the Department of Astronomy and Astrophysics at the University of California at Santa Cruz for their warm hospitality during a sabbatical stay in spring 2010. V.B. acknowledges support from NSF grant AST-0708795 and NASA ATFP grant NNX08AL43G. Part of the simulations were carried out at the Texas Advanced Computing Center (TACC), under TeraGrid allocation TG-AST090003.

Abel, T., Bryan, G. L., & Norman, M. L. 2002, Science, 295, 93

Aldrovandi, S. M. V. & Pequignot, D. 1973, A&A, 25, 137

Attwood, R. E., Goodwin, S. P., Stamatellos, D., & Whitworth, A. P. 2009, A&A, 495, 201

Barlow, S. G. 1984, PhD thesis, Univ. Colorado

Bate, M. R. 2009, MNRAS, 392, 590

Bate, M. R., Bonnell, I. A., & Price, N. M. 1995, MNRAS, 277, 362

Beers, T. C., & Christlieb, N. 2005, ARA&A, 43, 531

Benz, W., 1990, in ‘Proceedings of the NATO Advanced Research Workshop on The Numerical Modelling of Nonlinear Stellar Pulsations Problems and Prospects’, ed. J. R. Buchler, (Dordrecht: Kluwer), 269

Bonnell, I. A., Bate, M. R., Clarke, C. J., & Pringle, J. E. 1997, MNRAS, 285, 201

Bonnell, I. A., Bate, M. R., Clarke, C. J., & Pringle, J. E. 2001a, MNRAS, 323, 785

Bonnell, I. A., Clarke, C. J., Bate, M. R., & Pringle, J. E. 2001b, MNRAS, 324, 573

Bonnell, I. A., Vine, S. G., & Bate, M. R. 2004, MNRAS, 349, 735

Bonnell, I. A., & Bate, M. R. 2006, MNRAS, 370, 488

Bromm, V., Coppi, P. S., & Larson, R. B. 2002, ApJ, 564, 23

Bromm, V., & Larson, R. B. 2004, ARA&A, 42, 79

Bromm, V., & Loeb, A. 2004, New Astron. 9, 353

Bromm, V., Yoshida, N., Hernquist, L. & McKee, C. F. 2009, Nature, 459, 49

Chabrier, G. 2003, PASP, 115, 763

Clark, P. C., Glover, S. C. O., & Klessen, R. S. 2008, ApJ, 672, 757

Dalgarno, A., & Lepp, S. 1987, in Astrochemistry, ed. M. S. Vardya & S. P. Tarafdar (Dordrecht: Reidel), 109

Delgado-Donate, E. J., Clarke, C. J., & Bate, M. R. 2004, MNRAS, 347, 759

Dove, J. E., Rusk, A. C. M., Cribb, P. H., & Martin, P. G. 1987, ApJ, 318, 379

Ferland, G. J., Peterson, B. M., Horne, K., Welsh, W. F., & Nahar, S. N. 1992, ApJ, 387, 95

Flower, D. R., & Harris, G. J. 2007, MNRAS, 377, 705

Frebel, A., Johnson, J. L., & Bromm, V. 2009, MNRAS, 392, L50

Frommhold, L. 1993, Collision-induced Absorption in Gases, (Cambridge University Press: Cambridge)

Galli, D., & Palla, F. 1998, A&A, 335, 403

Gerlich, D. 1982, in ‘Symposium on Atomic and Surface Physics’, eds.  Lindinger, W., Howorka, F., Märk, T. D., (Dordrecht: Kluwer), p. 304

Glover, S. 2005, Space Sci. Rev., 117, 445

Glover, S. C. O. 2008, in First Stars III, eds.  B. O’Shea, A. Heger & T. Abel, (New York: AIP), 25; arXiv:0708.3086

Glover, S. C. O., & Abel, T. 2008, MNRAS, 388, 1627

Glover, S. C. O., & Brand, P. W. J. L. 2001, MNRAS, 321, 385

Glover, S. C. O., & Jappsen, A.-K. 2007, ApJ, 666, 1

Glover, S. C. O., Savin, D. W., & Jappsen, A.-K. 2006, ApJ, 640, 553

Goodwin, S. P., Whitworth, A. P., & Ward-Thompson, D. 2004a, A&A, 414, 633

Goodwin, S. P., Whitworth, A. P., & Ward-Thompson, D. 2004b, A&A, 423, 169

Greif, T. H., Johnson, J. L., Klessen, R. S., & Bromm, V. 2008, MNRAS, 387, 1021

Greif, T. H., Glover, S. C. O., Bromm, V., & Klessen, R. S. 2010, ApJ, 716, 510

Haiman, Z., Abel, T., & Rees, M. J. 2000, ApJ, 534, 11

Hosokawa, T., & Omukai, K. 2009, ApJ, 691, 823

Hummer, D. G., & Storey, P.J. 1998, MNRAS, 297, 1073

Huq, M. S., Doverspike, L. D., Champion, R. L., & Esaulov, V. A. 1982, J. Phys. B, 15, 951

Iwamoto, N., Umeda, H., Tominaga, N., Nomoto, K., & Maeda, K. 2005, Science, 309, 451

Janev, R. K., Langer, W. D., Evans, K., & Post, D. E. 1987, Elementary Processes in Hydrogen-Helium Plasmas, Springer

Jappsen, A.-K., Klessen, R. S., Larson, R. B., Li, Y., & [Mac Low]{}, M.-M. 2005, A&A, 435, 611

Jappsen, A.-K., Low, M.-M. M., Glover, S. C. O., Klessen, R. S., & Kitsionas, S. 2009, ApJ, 694, 1161

Jappsen, A.-K., Klessen, R. S., Glover, S. C. O., & Mac Low, M.-M. 2009, ApJ, 696, 1065

Joggerst, C. C., Almgren, A., Bell, J., Heger, A., Whalen, D. J., & Woosley, S. E. 2010, ApJ, 709, 11

Johnson, J. L., & Bromm, V. 2006, MNRAS, 366, 247

Johnson, J. L., Greif, T. H., & Bromm, V. 2008, MNRAS, 388, 26

Karlsson, T., Johnson, J. L., & Bromm, V. 2008, ApJ, 679, 6

Karpas, Z., Anicich, V., & Huntress, W. T. 1979, J. Chem. Phys, 70, 2877

Kimura, M., Lane, N. F., Dalgarno, A., & Dixson, R. G. 1993, ApJ, 405, 801

Klessen, R. S., & Burkert, A. 2000, ApJS, 128, 287

Klessen, R. S., & Burkert, A. 2001, ApJ, 549, 386

Klessen, R. S. 2001, ApJ, 550, L77

Klessen, R. S., & Hennebelle, P. 2009, arXiv:0912.0288

Kreckel, H., Bruhns, H., Čížek, M., Glover, S. C. O., Miller, K. A., Urbain, X., & Savin, D. W., 2010, Science, 329, 69

Kroupa, P. 2002, Science, 295, 82

Krumholz, M. R., Klein, R. I., & McKee, C. F. 2007, ApJ, 656, 959

Krumholz, M. R., Klein, R. I., McKee, C. F., Offner, S. S. R., & Cunningham, A. J. 2009, Science, 323, 754

Larson, R. B. 1981, MNRAS, 194, 809

Machida, M. N. 2008, ApJ, 682, L1

Machida, M. N., Omukai, K., Matsumoto, T., & Inutsuka, S.-i. 2008, ApJ, 677, 813

Machida, M. N., Omukai, K., Matsumoto, T., & Inutsuka, S.-I. 2009, MNRAS, 399, 1255

, M.-M., & Klessen, R. S. 2004, Rev. Mod. Phys., 76, 125

Martin, P. G., Schwarz, D. H., & Mandy, M. E. 1996, ApJ, 461, 265

Martin, P. G., Keogh W. J., & Mandy, M. E. 1998, ApJ, 499, 793

McKee, C. F., & Tan, J. C. 2008, ApJ, 681, 771

Mielke, S. L., Peterson, K. A., Schwenke, D. W., Garrett, B. C., Truhlar, D. G., Michael, J. V., Su, M.-C., & Sutherland, J. W. 2003, Phys. Rev. Lett., 91, 063201

Molaro, P. 2008, in ‘Pathways Through and Eclectic Universe’, ASP Conf. Ser. 390, eds. J. H. Knapen, T. J. Mahoney & A. Vazdekis, (San Francisco:ASP), 472

Nagakura, T., & Omukai, K. 2005, MNRAS, 364, 1378

Neufeld, D. A., & Kaufman, M. J. 1993, ApJ, 418, 263

Omukai, K., & Nishi, R. 1999, ApJ, 518, 64

Omukai, K., Tsuribe, T., Schneider, R., & Ferrara, A. 2005, ApJ, 626, 627

O’Shea, B. W., & Norman, M. L. 2007, ApJ, 654, 66

O’Shea, B. W., McKee, C. F., Heger, A., & Abel, T. 2008, First Stars III, 990, 13

Palla, F., Salpeter, E. E., & Stahler, S. W. 1983, ApJ, 271, 632

Peart, B., & Hayton, D. A. 1994, J. Phys. B, 27, 2551

Poulaert, G., Brouillard, F., Claeys, W., McGowan, J. W., & [Van Wassenhove]{}, G. 1978, J. Phys. B, 11, L671

Ramaker, D. E., & Peek, J. M. 1976, Phys. Rev. A, 13, 58

Ripamonti, E., & Abel, T. 2004, MNRAS, 348, 1019

Savin, D. W. 2002, ApJ, 566, 599

Savin, D. W., Krstic, P. S., Haiman, Z., & Stancil, P. C. 2004, ApJ, 606, L167; erratum ApJ, 607, L147

Schleicher, D. R. G., Banerjee, R., Sur, S., Arshakian, T. G., Klessen, R. S., Beck, R., & Spaans, M. 2010, MNRAS, submitted; arXiv:1003.1135

Schmeja, S., & Klessen, 2004, A&A, 419, 405

Schneider, I. F., Dulieu, O., Giusti-Suzor, A., & Roueff, E. 1994, ApJ, 424, 983; erratum ApJ, 486, 580

Schulz, G. J., & Asundi, R. K. 1967, Phys. Rev., 158, 25

Shapiro, P. R., & Kang, H. 1987, ApJ, 318, 32

Shavitt, I. 1959, J. Chem. Phys., 31, 1359

Springel, V. 2005, MNRAS, 364, 1105

Stacy, A., Greif, T. H., & Bromm, V. 2010, MNRAS, 403, 45

Tan, J. C., & Blackman, E. G. 2004, ApJ, 603, 401

Tan, J. C., & McKee, C. F. 2008, in First Stars III, eds.  B. O’Shea, A. Heger & T. Abel, (New York: AIP), 47

Trenti, M., & Stiavelli, M. 2009, ApJ, 694, 879

Trevisan, C. S., & Tennyson, J. 2002a, Plasma Phys. Controlled Fusion, 44, 1263

Tsuribe, T., & Omukai, K. 2008, ApJ, 676, L45

Tumlinson, J., Venkatesan, A., & Shull, J. M. 2004, ApJ, 612, 602

Tumlinson, J. 2007a, ApJ, 664, L63

Tumlinson, J. 2007b, ApJ, 665, 1361

Tumlinson, J. 2010, ApJ, 708, 1398

Turk, M. J., Abel, T., & O’Shea, B. W. 2009, Science, 325, 601

Turk, M. J., Clark, P. C., Glover, S. C. O., Greif, T. H., Abel, T., Klessen, R. S., & Bromm, V. 2010, ApJ, [*accepted*]{}, arXiv:1011.1908

Walkauskas, L. P., & Kaufman, F. 1975, Symp. Int. Combust. Proc., 15, 691

Wise, J. H., & Abel, T. 2007, ApJ, 665, 899

Wise, J. H., & Abel, T. 2008, ApJ, 685, 40

Wishart, A. W. 1979, MNRAS, 187, 59P

Wolniewicz, L., Simbotin, I., & Dalgarno, A. 1998, ApJS, 115, 293

Yoshida, N., Omukai, K., Hernquist, L., & Abel, T. 2006, ApJ, 652, 6

Yoshida, N., Omukai, K. & Hernquist, L. 2007, ApJ, 667, L117

Yoshida, N., Omukai, K., & Hernquist, L. 2008, Science, 321, 669

Zygelman, B., Dalgarno, A., Kimura, M., & Lane, N. F. 1989, Phys. Rev. A, 40, 2340

Chemistry {#num:chem}
=========

To model the chemical evolution of the metal-free gas, we use the network for primordial hydrogen, helium and deuterium chemistry detailed in Table \[tab:chem\]. Our treatment of the hydrogen and helium chemistry largely follows @ga08, but our treatment of the deuterium chemistry is significantly simplified. This simplification arises from our neglect of ${\rm D^{-}}$, ${\rm HD^{+}}$ and ${\rm D_{2}}$, none of which play a significant role in controlling the HD abundance in the physical conditions relevant to our present study. For the most part, our choice of rate coefficients also follows @ga08.[^2] We adopt case B rate coefficients for the recombination of H$^{+}$ and He$^{++}$, and treat He$^{+}$ recombination as decribed in Section 2.1.4 of @ga08. We adopt rate coefficients from @gp98 for the associative detachment of H$^{-}$ ions by atomic hydrogen (reaction 2), and for the mutual neutralization of H$^{-}$ by H$^{+}$ (reaction 5). We also note that the uncertainties in these reaction rate coefficients discussed in @gsj06 are unlikely to be significant in the high density, low ionization gas modelled in our present simulations. For the rate coefficient of reaction 30, the three-body formation of H$_{2}$ with atomic hydrogen as the third body $${\rm H} + {\rm H} + {\rm H} \rightarrow {\rm H_{2}} + {\rm H},
\label{tbh2a}$$ we adopt the rate coefficient proposed by @g08. The value of the rate coefficient for this reaction is highly uncertain [@g08; @tcg10], but our chosen value lies intermediate between the fastest and slowest rates to be found in the literature (@fh07 and @abn02, respectively), and is based on the best currently available data for the inverse reaction (i.e. collisional dissociation of H$_{2}$ by H, reaction 9). The influence of the uncertainty in this rate coefficient has been studied in detail elsewhere [@tcg10]. To fix the rate coefficient for the analogous reaction $${\rm H} + {\rm H} + {\rm H_{2}} \rightarrow {\rm H_{2}} + {\rm H_{2}},$$ we follow @pss83 and assume that the rate coefficient is one-eighth of the size of the rate coefficient adopted for reaction \[tbh2a\]. For three-body formation in which neutral helium is the third body, we follow @ga08 and use a rate coefficient originally taken from @wk75. We assume, following @fh07, that the rate coefficients for the three-body formation of HD (reactions 43–45) are the same as those for the corresponding H$_{2}$ formation reactions (nos. 30–32).

We adopt values for the rate coefficients for collisional dissociation of H$_{2}$ by H, H$_{2}$ and He that are consistent with our choices for the three-body association rates in the sense that each pair of rate coefficients individually satisfies $$\frac{k_{\rm form}}{k_{\rm dest}} = K,$$ where $k_{\rm form}$ is the rate coefficient for three-body association, $k_{\rm dest}$ is the rate coefficient for collisional dissociation, and where the equilibrium constant $K$ is given in all three cases by [@fh07] $$K = 1.05 \times 10^{-22} T^{-0.515} \exp \left(\frac{52000}{T} \right).$$ This procedure is necessary in order to ensure the correct chemical behaviour of the gas at high densities and temperatures, where the H$_{2}$ formation and destruction timescales are both short, and most of the gas is close to chemical equilibrium.

[llc]{} 1 & ${{\rm H}}+ {{\rm e^{-}}}\rightarrow  {{\rm H}^{-}}+ \gamma $ & 1\
2 & ${{\rm H}^{-}}+ {{\rm H}}\rightarrow {{\rm H_{2}}}+ {{\rm e^{-}}}$ & 2\
3 & ${{\rm H}}+ {{\rm H}^{+}}\rightarrow  {{\rm H_{2}^{+}}}+ \gamma $ & 3\
4 & ${{\rm H}}+ {{\rm H_{2}^{+}}}\rightarrow {{\rm H_{2}}}+ {{\rm H}^{+}}$ & 4\
5 & ${{\rm H}^{-}}+ {{\rm H}^{+}}\rightarrow {{\rm H}}+ {{\rm H}}$ & 2\
6 & ${{\rm H_{2}^{+}}}+ {{\rm e^{-}}}\rightarrow  {{\rm H}}+ {{\rm H}}$ & 5\
7 & ${{\rm H_{2}}}+ {{\rm H}^{+}}\rightarrow  {{\rm H_{2}^{+}}}+ {{\rm H}}$ & 6\
8 & ${{\rm H_{2}}}+ {{\rm e^{-}}}\rightarrow  {{\rm H}}+ {{\rm H}}+  {{\rm e^{-}}}$ & 7\
9 & ${{\rm H_{2}}}+ {{\rm H}}\rightarrow  {{\rm H}}+ {{\rm H}}+ {{\rm H}}$ & 8\
10 & ${{\rm H_{2}}}+ {{\rm H_{2}}}\rightarrow  {{\rm H_{2}}}+ {{\rm H}}+ {{\rm H}}$ & 9, 10\
11 & ${{\rm H_{2}}}+ {{\rm He}}\rightarrow {{\rm H}}+ {{\rm H}}+ {{\rm He}}$ & 11\
12 & ${{\rm H}}+ {{\rm e^{-}}}\rightarrow  {{\rm H}^{+}}+ {{\rm e^{-}}}+ {{\rm e^{-}}}$ & 12\
13 & ${{\rm H}^{+}}+ {{\rm e^{-}}}\rightarrow {{\rm H}}+  \gamma$ & 13\
14 & ${{\rm H}^{-}}+ {{\rm e^{-}}}\rightarrow  {{\rm H}}+ {{\rm e^{-}}}+ {{\rm e^{-}}}$ & 12\
15 & ${{\rm H}^{-}}+ {{\rm H}}\rightarrow  {{\rm H}}+ {{\rm H}}+ {{\rm e^{-}}}$ & 12\
16 & ${{\rm H}^{+}}+ {{\rm H}^{-}}\rightarrow  {{\rm H_{2}^{+}}}+ {{\rm e^{-}}}$ & 14\
17 & ${{\rm He}}+ {{\rm e^{-}}}\rightarrow  {{\rm He^{+}}}+ {{\rm e^{-}}}+ {{\rm e^{-}}}$ & 12\
18 & ${{\rm He^{+}}}+ {{\rm e^{-}}}\rightarrow {{\rm He^{++}}}+ {{\rm e^{-}}}+ {{\rm e^{-}}}$ & 12\
19 & ${{\rm He^{+}}}+ {{\rm e^{-}}}\rightarrow {{\rm He}}+ \gamma$ & 15, 16\
20 & ${{\rm He^{++}}}+ {{\rm e^{-}}}\rightarrow {{\rm He^{+}}}+ \gamma$ & 13\
21 & ${{\rm H}^{-}}+ {{\rm H_{2}^{+}}}\rightarrow {{\rm H_{2}}}+ {{\rm H}}$ & 17\
22 & ${{\rm H}^{-}}+ {{\rm H_{2}^{+}}}\rightarrow {{\rm H}}+ {{\rm H}}+ {{\rm H}}$ & 17\
23 & ${{\rm H_{2}}}+ {{\rm e^{-}}}\rightarrow {{\rm H}^{-}}+ {{\rm H}}$ & 18\
24 & ${{\rm H_{2}}}+ {{\rm He^{+}}}\rightarrow {{\rm He}}+ {{\rm H}}+ {{\rm H}^{+}}$ & 19\
25 & ${{\rm H_{2}}}+ {{\rm He^{+}}}\rightarrow {{\rm H_{2}^{+}}}+ {{\rm He}}$ & 19\
26 & ${{\rm He^{+}}}+ {{\rm H}}\rightarrow {{\rm He}}+ {{\rm H}^{+}}$ & 20\
27 & ${{\rm He}}+ {{\rm H}^{+}}\rightarrow {{\rm He^{+}}}+ {{\rm H}}$ & 21\
28 & ${{\rm He^{+}}}+ {{\rm H}^{-}}\rightarrow {{\rm He}}+ {{\rm H}}$ & 22\
29 & ${{\rm He}}+ {{\rm H}^{-}}\rightarrow {{\rm He}}+ {{\rm H}}+ {{\rm e^{-}}}$ & 23\
30 & ${{\rm H}}+ {{\rm H}}+ {{\rm H}}\rightarrow {{\rm H_{2}}}+ {{\rm H}}$ & 24\
31 & ${{\rm H}}+ {{\rm H}}+ {{\rm H_{2}}}\rightarrow {{\rm H_{2}}}+ {{\rm H_{2}}}$ & 25\
32 & ${{\rm H}}+ {{\rm H}}+ {{\rm He}}\rightarrow {{\rm H_{2}}}+ {{\rm He}}$ & 26\
33 & ${{\rm D}^{+}}+ {{\rm e^{-}}}\rightarrow {\rm{D}}+ \gamma$ & 27\
34 & ${\rm{D}}+ {{\rm H}^{+}}\rightarrow  {{\rm H}}+ {{\rm D}^{+}}$ & 28\
35 & ${{\rm H}}+ {{\rm D}^{+}}\rightarrow  {\rm{D}}+ {{\rm H}^{+}}$ & 28\
36 & ${{\rm H_{2}}}+ {\rm{D}}\rightarrow {\rm{HD}}+ {{\rm H}}$ & 29\
37 & ${{\rm H_{2}}}+ {{\rm D}^{+}}\rightarrow {\rm{HD}}+ {{\rm H}^{+}}$ & 30\
38 & ${\rm{HD}}+ {{\rm H}}\rightarrow  {{\rm H_{2}}}+ {\rm{D}}$ & 31\
39 & ${\rm{HD}}+ {{\rm H}^{+}}\rightarrow {{\rm H_{2}}}+ {{\rm D}^{+}}$ & 30\
40 & ${\rm{D}}+ {{\rm e^{-}}}\rightarrow {{\rm D}^{+}}+ {{\rm e^{-}}}+ {{\rm e^{-}}}$ & 27\
41 & ${{\rm He^{+}}}+ {\rm{D}}\rightarrow {{\rm D}^{+}}+ {{\rm He}}$ & 32\
42 & ${{\rm He}}+ {{\rm D}^{+}}\rightarrow {\rm{D}}+ {{\rm He^{+}}}$ & 32\
43 & ${\rm{D}}+ {{\rm H}}+ {{\rm H}}\rightarrow {\rm{HD}}+ {{\rm H}}$ & See text\
44 & ${\rm{D}}+ {{\rm H}}+ {{\rm H_{2}}}\rightarrow {\rm{HD}}+ {{\rm H_{2}}}$ & See text\
45 & ${\rm{D}}+ {{\rm H}}+ {{\rm He}}\rightarrow {\rm{HD}}+ {{\rm He}}$ & See text\

Cooling function and thermodynamics {#sec:chem}
===================================

The dominant coolant in our simulations is molecular hydrogen. To model rotational and vibrational line emission from H$_{2}$, we use the detailed cooling function described in @ga08, that includes contributions from collisions of H$_{2}$ with H, He, H$_{2}$, protons and electrons. At densities $n \simgreat 10^{9} \: {\rm cm^{-3}}$, the strongest of the H$_{2}$ lines become optically thick, reducing its effectiveness as a coolant. To model H$_{2}$ cooling in this regime, we use an approach based on the Sobolev approximation that was first used in models of primordial star formation by @yoha06. We write the H$_{2}$ cooling rate in optically thick gas as $$\Lambda_{\rm H_{2}} = \sum_{\rm u,l} \Delta E_{\rm ul} A_{\rm ul} \beta_{\rm esc, ul} 
n_{\rm u},$$ where $n_{\rm u}$ is the number density of hydrogen molecules in upper energy level $u$, $\Delta E_{\rm ul}$ is the energy difference between this upper level and a lower level $l$, $A_{\rm ul}$ is the spontaneous radiative transition rate for transitions between $u$ and $l$, and $\beta_{\rm esc, ul}$ is the escape probability associated with this transition, i.e. the probability that the emitted photon can escape from the region of interest. We take values for the level energies from the compilation made available by P. G. Martin on his website[^3] and for the radiative transition rates from @wsd98, and we fix $n_{\rm u}$ by assuming that the H$_{2}$ level populations are in LTE. The problem of modelling optically thick H$_{2}$ cooling thereby reduces to one of computing the escape probability for each transition. We follow @yoha06 and write the escape probability for the transition ${\rm u \rightarrow l}$ as $$\beta_{\rm esc, ul} = \frac{1 - \exp(-\tau_{\rm ul})}{\tau_{\rm ul}},$$ where we approximate $\tau_{\rm ul}$ as $$\tau_{\rm ul} = \alpha_{\rm ul} L_{\rm s} \label{tau}$$ where $\alpha_{\rm ul}$ is the line absorption coefficient and $L_{\rm s}$ is the Soloblev length. In the classical, one-dimensional spherically symmetric case, the Sobolev length is given by $$L_{\rm s} = \frac{v_{\rm th}}{|{\rm d}v_{\rm r} / {\rm d}r|},$$ where $v_{\rm th}$ is the thermal velocity, and ${\rm d}v_{\rm r} / {\rm d}r$ is the radial velocity gradient. In our inherently three-dimensional flow, we generalize this as $$L_{\rm s} = \frac{v_{\rm th}}{|\nabla \cdot {\mathbf v}|},$$ following @nk93. If the velocity dispersion of the gas is very small, then $L_{\rm s}$ can become very large, much larger than the size of the collapsing core. To ensure that we do not reduce the H$_{2}$ cooling rate in this case to an artificially low value, we take as our actual length scale in Equation \[tau\] the smallest of the Soloblev length and the local Jeans length, $L_{\rm J}$.

Since the line absorption coefficient $\alpha_{\rm ul}$ is linearly proportional to the number density of H$_{2}$, we can write $\tau_{\rm ul}$ as $$\tau_{\rm ul} = \left(\frac{\alpha_{\rm ul}}{n_{\rm H_{2}}} \right) N_{\rm H_{2}, eff},$$ where $N_{\rm H_{2}, eff} \equiv n_{\rm H_{2}} L_{\rm s}$ is an effective H$_{2}$ column density, and where $\alpha_{\rm ul} / n_{\rm H_{2}}$ is a function only of temperature. We therefore tabulate the cooling rate per H$_{2}$ molecule in the optically thick limit as a function of two parameters: the gas temperature $T$ and the effective H$_{2}$ column density $N_{\rm H_{2}, eff}$, and compute cooling rates during the simulations by interpolation from a pre-generated look-up table.

At densities $n > 10^{14} \: {\rm cm^{-3}}$, a second form of H$_{2}$ cooling becomes important, called collision-induced emission. Although H$_{2}$ molecules have no electric dipole, the interacting pair in a collision of H$_{2}$ with H, He or H$_{2}$ briefly acts as a ‘supermolecule’ with an non-zero electric dipole, and a hence a non-zero probability of emitting or absorbing a photon through a dipole transition. Because the collision time is very short, the resulting collision-induced transition lines are very broad, effectively merging into a continuum [for more details, see e.g. @fromm93]. We include this process in our cooling function, using a rate taken from @ra04. We crudely account for the reduction of the CIE cooling rate by continuum absorption at very high number densities using the following prescription (M. Turk, private communication) $$\Lambda_{\rm CIE, thick} = \Lambda_{\rm CIE, thin} \times {\rm min}\left(\frac{1 - 
e^{-\tau_{\rm CIE}}}{\tau_{\rm CIE}}, 1 \right)$$ where $$\tau_{\rm CIE} = \left(\frac{n}{1.4 \times 10^{16} \: {\rm cm^{-3}}} \right)^{2.8}.$$ However, we note that optical depth effects do not strongly affect the CIE cooling rate for densities below our threshold for sink particle creation (see §\[sec:ics\]), and so the approximate nature of this opacity cutoff is unlikely to significantly affect our results.

In addition to H$_{2}$ ro-vibrational line emission and CIE cooling, our cooling function contains a number of other radiative processes: electronic excitation of H, He and He$^{+}$, cooling from the recombination of H$^{+}$ and He$^{+}$, Compton cooling and bremsstrahlung. Details of our treatment of these processes can be found in @gj07.

We also account for changes in the thermal energy of the gas due to changes in its chemical makeup. Specifically, we include the effects of cooling due to the collisional ionization of H, He and He$^{+}$, and due to the destruction of H$_{2}$ by charge transfer and by collisional dissociation, as well as heating due to the three-body H$_{2}$ formation. The balance between cooling due to H$_{2}$ collisional dissociation and heating due to three-body H$_{2}$ formation plays a very important role in regulating the thermal evolution of the gas at densities $n \simgreat 10^{8} \: {\rm cm^{-3}}$.

Evolution of the Bonnor-Ebert spheres {#ICdiscuss}
=====================================

Given our somewhat arbitrary choice of initial conditions, it is prudent to ask whether they are applicable to primordial star formation, and in particular, whether the choice of the initial density profile affects the details of the collapse. In Fig. \[fig:profiles\], we show, for the Pop III.1 clouds with 0.1 and 0.4 $c_{\rm s}$ turbulence, the number density and enclosed mass as a function of radius and how the radial velocity and temperature varies as a function of the enclosed mass. In the case of the 0.1$c_{\rm s}$ cloud, the radial profiles (mass and number density) are much like those presented for the standard Pop III studies [@abn02; @yoha06], in which fully cosmological initial conditions were used. In particular, by the time that the first sink particle forms in this run, the gas has developed the same $n \propto r^{-2.2}$ radial density profile as found in the fully cosmological runs, demonstrating that the density profile at this stage is insensitive to the shape of the initial density profile.

We do however see differences between the radial velocity and temperature profiles, compared to those published in the literature. The infall velocities at this stage are somewhat higher than those seen in fully cosmological studies, but it should be noted that we do not include systemic rotation in our study, and so a major source of support is missing. However, aside from the higher values of the radial velocity, the shape of the profile is also different. In the runs presented in this study, the infall velocities of the shells enclosing 1 and greater are systematically higher than those seen in the cosmological simulations published to date, with a peak in the radial velocity profile at about 10-100 (depending on the level of turbulent support) instead of around 0.3 . Some of this difference in the radial velocity profiles is reflected in the somewhat steeper than reported temperature profiles, since the gas will reacts to the increased $p {\rm d} V$ heating. However the temperature profiles are also affected by the differences in the chemistry, and in particular, our choice of the three-body H$_{2}$ formation rate.

So the obvious question is, do our high infall velocities make fragmentation more likely than it would be reality? To address this, we performed another Pop III.1 simulation with 0.4 $c_{\rm s}$ turbulence, but this time we set the confining external pressure to a factor of 4 less than was used in the original simulations. The radial and enclosed mass profiles at the point just before the first sink particle is formed are shown as the dashed lines in Fig. \[fig:profiles\] (right-hand panel). We see that radial velocity profile is affected by the change in boundary pressure, resulting in a factor of 4 reduction in the infall speed. Compared to the original run, which formed 31 sink particles, this simulation forms 22 sinks when 10 percent of the gas has been accreted. While this suggests that the fragmentation we see is sensitive to the details of the collapse, and we are likely overestimating the number of fragments formed (at least at the scale of our sink particles) when one compares to the standard Pop III collapse profiles, it seems that this is not the main cause of the fragmentation. The gas still fragments, and the processes leading up to the fragmentation are the same as those discussed in §\[sec:evol\].

Finally, we note that there is also new evidence to suggest that the radial collapse profiles seen in the literature so far, may not be representative of all Pop III star formation. For example, while the radial density profiles from the studies by @abn02 and @yoha06 are similar to the newer, higher-resolution, cosmological simulations of @tao09 and @tcg10, the radial (or enclosed mass) profiles of the other quantities in the newer calculations show considerable variation. In particular, the velocity and temperature profiles in @tcg10 lie somewhere in between those found by @yoha06 and our 0.4 $c_{\rm s}$ turbulent cloud. A very high-resolution study of the onset of the initial gravitational instability in the baryonic component of dark matter minihalos, may help to establish whether these differences between the studies of, say, @abn02 and @tao09 are due to the increased resolution and improved physical treatment used in the latter study (e.g. the inclusion of the effects of three-body H$_2$ formation heating), or just reflect the effects of cosmic variance.

![\[fig:profiles\] Radial profiles from the Pop III.1 clouds with $\Delta v_{\rm turb} = 0.1 c_{\rm s}$ (left panel) and $\Delta v_{\rm turb} = 0.4 c_{\rm s}$ (right panel). The lines show the initial conditions (solid line), the state of the gas just before the creation of the first sink (dotted line). In the case of the $\Delta v_{\rm turb} = 0.4 c_{\rm s}$ cloud, we also show the results of a collapse that has a factor of 4 less external pressure (dashed line). The dot-dashed line shows a slope of $n \propto r^{-2.2}$. The crosses and circular points denote approximate values taken, respectively, from @yoha06 and @abn02](fig13a "fig:"){height="3.4in"} ![\[fig:profiles\] Radial profiles from the Pop III.1 clouds with $\Delta v_{\rm turb} = 0.1 c_{\rm s}$ (left panel) and $\Delta v_{\rm turb} = 0.4 c_{\rm s}$ (right panel). The lines show the initial conditions (solid line), the state of the gas just before the creation of the first sink (dotted line). In the case of the $\Delta v_{\rm turb} = 0.4 c_{\rm s}$ cloud, we also show the results of a collapse that has a factor of 4 less external pressure (dashed line). The dot-dashed line shows a slope of $n \propto r^{-2.2}$. The crosses and circular points denote approximate values taken, respectively, from @yoha06 and @abn02](fig13b "fig:"){height="3.4in"}

[^1]: Note that the term ‘fragmentation’ is often used in the literature even in the case where only one object forms at the center of the clump, even though, strictly speaking, the term refers to a system separating into several distinct parts with separate evolutionary paths.

[^2]: Note that there are two typographical errors that we are aware of in the set of reaction rate coefficients listed in Table A1 in @ga08. First, the fitting function used to describe the rate coefficient for the reaction ${\rm H_{2}} + {\rm H^{+}} \rightarrow {\rm H_{2}^{+}} + {\rm H}$ (reaction number 7 in their table) should use natural logarithms, and not base 10 logarithms, as listed. However, the listed fitting coefficients are correct. Second, the temperature dependence of the reaction ${\rm H_{2}} + {\rm He} \rightarrow
    {\rm He} + {\rm H} + {\rm H^{+}}$ (their reaction 24) should be $\exp(-35/T)$, and not $\exp(+35/T)$.

[^3]: http://www.cita.utoronto.ca/$\sim$pgmartin/h2.html
