A NEW APPROACH TO DETERMINE OPTICALLY THICK H₂ COOLING AND ITS EFFECT ON PRIMORDIAL STAR FORMATION

Since H₂ is a symmetric, diatomic molecule, it has no permanent dipole moment and can therefore only radiate via ro-vibrational transitions (Le Bourlot et al. 1999). In the optically thin regime (n ⩽ 10⁹ cm⁻³) we expect all radiation to escape the cloud freely and the cooling rate can be calculated fairly accurately as a function of the local density, temperature, and chemical composition of the gas. A number of different parameterizations of the H₂ cooling rate in terms of these quantities are available in the literature. In this work, we use the rates given in Glover & Abel (2008). In the optically thin case we assume that photons, which are emitted in the ro-vibrational transitions, can escape the cloud without being scattered or absorbed and therefore transport thermal energy outward efficiently. In the optically thick case, however, scattering and absorption events might capture the photons and hence decrease the overall cooling efficiency. One generally assumes complete redistribution of the frequency between these scattering events and uses the "escape probability method" in order to solve this problem. Most important for our further discussion are the applications of the escape probability method by Yoshida et al. (2006) and Clark et al. (2011a) to the case of optically thick cooling in primordial gas. The cooling rate in an optically thick medium is given by

$$\begin{equation} \Lambda _{\mathrm{H}_2\,,\mathrm{thick}} = \sum _{l,u} E_{lu}A_{lu} \beta _{\mathrm{esc},lu} n_u, \end{equation} \tag{ 1 }$$

where E_lu = hν_lu is the energy separation between the lower level l and upper level u, A_lu is the spontaneous radiative transition rate for transitions between u and l, β_{esc, lu} is the probability for an emitted line photon to escape without absorption, and n_u is the population density of hydrogen molecules in the upper level. Following Yoshida et al. (2006), we assume all energy levels to be populated according to local thermodynamic equilibrium and we consider rotational levels from J = 0 to 20 and vibrational levels ν = 0, 1, 2. At the temperatures of interest, the contribution of other levels will be negligibly small (<10⁻³). We take values for the level energies from the compilation made available by P. G. Martin on his Web site¹ and the radiative transitions rates from Wolniewicz et al. (1998). Since all other quantities are known for typical conditions in primordial gas, the remaining task is the determination of the photon escape probability.

Based on the work by Zanstra (1934), Sobolev (1947) derived the escape probability for constant velocity gradients in his study of expanding envelopes (for an English translation see Sobolev 1960). This approximation has been reviewed during the last few decades (Castor 1970; Hummer & Rybicki 1982; Yoshida et al. 2006) and is widely used in simulations of primordial star formation (Yoshida et al. 2006, 2008; Turk et al. 2011; Clark et al. 2011a, 2011b; Greif et al. 2011; Wolcott-Green et al. 2011; Greif et al. 2012; Hirano & Yoshida 2013; Stacy et al. 2013; Stacy & Bromm 2013; Greif et al. 2013; Hirano et al. 2014; Greif et al. 2014a; Stacy & Bromm 2014). Following the derivation by Emerson (1996), the photon escape probability is given by

$$\begin{equation} \beta _{\mathrm{esc},lu} = \frac{1-e^{-\tau _{lu}}}{\tau _{lu}}, \end{equation} \tag{ 2 }$$

where τ_lu is the opacity at the line center. The absorption coefficient for a transition from the lower to the upper level is given by

$$\begin{equation} \alpha _{lu} = \frac{E_{lu}}{4 \pi }n_l B_{lu} \left[ 1- \exp \left(\frac{-E_{lu}}{k_B T} \right) \right] \phi (\nu), \end{equation} \tag{ 3 }$$

where B_lu is the Einstein coefficient for absorption and ϕ(ν) is the line profile function. The opacity can be written as

$$\begin{equation} \tau _{lu} = \alpha _{lu} L_\mathrm{char}, \end{equation} \tag{ 4 }$$

where L_char is a characteristic length scale (Castor 1970; Goldreich & Kwan 1974; de Jong et al. 1975; Stahler & Palla 2005). According to Equation (3), $\alpha _{lu}/n_{\mathrm{H}_2\,}$ is only a function of temperature and, following Clark et al. (2011a), we can express the optical depth by

$$\begin{equation} \tau _{lu} = \left(\frac{\alpha _{lu}}{n_{\mathrm{H}_2\,}} \right) n_{\mathrm{H}_2\,} L_\mathrm{char}, \end{equation} \tag{ 5 }$$

where

$$\begin{equation} n_{\mathrm{H}_2\,} L_\mathrm{char} = N_{\mathrm{H}_2\,,\mathrm{eff}} \end{equation} \tag{ 6 }$$

defines an effective column density. Hence, the last task for the determination of optically thick cooling is the calculation of the characteristic length and the associated effective column density. Knowing these quantities, one can generally determine the angle-dependent escape probability and afterward average it over all lines of sight and all relevant ro-vibrational lines. Since H₂ cooling is due to a number of lines without a single dominant line (Haardt et al. 2002), we can directly average the escape probability over these lines. Furthermore, one often assumes spherical symmetry, whereas Yoshida et al. (2006) proposed an average over three orthogonal directions

$$\begin{equation} \beta = \frac{\beta _x + \beta _y + \beta _z}{3}. \end{equation} \tag{ 7 }$$

The averaged escape probability finally relates the optically thin and optically thick cooling rate by

$$\begin{equation} \Lambda _{\mathrm{H}_2\,,\mathrm{thick}} = \beta \Lambda _{\mathrm{H}_2\,,\mathrm{thin}} \end{equation} \tag{ 8 }$$

and is therefore also known as the "opacity correction." In order to determine this escape probability, we have to understand the dynamics of the cloud. If the photon is emitted in the center of the cloud and an envelope of gas is moving toward it with a constant radial velocity gradient dv_r/dr, then the photon observes the spectral lines of the envelope to be Doppler-shifted with respect to its rest frame. According to Sobolev (1947), a photon is not absorbed and can escape freely, if the spectral lines of a possibly absorbing H₂ molecule are shifted by more than one thermal line width. The line width of thermal line broadening is given by

$$\begin{equation} \Delta \nu _{\mathrm{th}} = \nu _0 \frac{v_{\mathrm{th}}}{c} = \frac{\nu _0}{c} \sqrt{\frac{2k_BT}{m_{\mathrm{H}_2}}}, \end{equation} \tag{ 9 }$$

where ν₀ is the central frequency of the line, v_th is the thermal velocity of molecular hydrogen, T is the temperature, and m_H is the mass of a hydrogen atom. Using this, we can determine the characteristic distance L_char beyond which the Sobolev criterion is fulfilled. This length scale is typically known as the Sobolev length

$$\begin{equation} L_{\mathrm{s}}= \frac{v_\mathrm{th}}{|\mathrm{d}v_\mathrm{r} / \mathrm{d}r|}. \end{equation} \tag{ 10 }$$

Phrased differently, all relevant matter that might reabsorb a photon is within its Sobolev length. Assuming a constant density within this Sobolev length, the effective column density can be determined by

$$\begin{equation} N_{\mathrm{H}_2\,,\mathrm{eff}} = n_{\mathrm{H}_2\,} L_{\mathrm{s}}. \end{equation} \tag{ 11 }$$

In order to capture the three-dimensional dynamics of the collapse, one normally uses

$$\begin{equation} L_{\mathrm{s}}= \frac{v_\mathrm{th}}{|\mathbf {\nabla }\cdot {\mathbf {v}}|} \end{equation} \tag{ 12 }$$

for the determination of the Sobolev length (Neufeld & Kaufmann 1993). However, a fundamental problem of the Sobolev approximation was already mentioned by several authors: both the velocity gradient and the number density have to be constant within one Sobolev length (Lucy 1971; Bujarrabal et al. 1980; Hummer & Rybicki 1992; Neufeld & Kaufmann 1993; Wolcott-Green et al. 2011). As we will see below, this assumption is generally not valid.

For simplicity, Zanstra (1934) assumed the absorption coefficient α_lu to be zero outside the interval [ν₀ − Δν_th, ν₀ + Δν_th] and Sobolev (1947) used the same simplification. Following this approach, a photon can escape freely from the optically thick gas after one Sobolev length L_s, because it will not be reabsorbed thereafter. The actual absorption probability however, is related to the true overlap of spectral lines. In the present context, the shape of the H₂ lines is dominated by thermal broadening and can therefore be described by a normalized Gaussian profile. If one line with central frequency ν₀ is Doppler-shifted to the frequency ν, the relative line displacement in units of the thermal gas velocity between these two cases is given by

$$\begin{equation} x=\frac{\nu - \nu _0}{\Delta \nu _{\mathrm{th}}}. \end{equation} \tag{ 13 }$$

Accordingly, the relative overlap of these two line profiles can be calculated by

$$\begin{equation} o(x)=\int _{- \infty } ^{x/2} \frac{2}{\sqrt{2 \pi }} (e^{- \nu ^2/2} - e^ {- (\nu -x)^2/2} ) \rm d \nu. \end{equation} \tag{ 14 }$$

The overlap of spectral lines for the special case x = 1 is illustrated in Figure 1. This displacement represents the case after one Sobolev length, though the relative overlap is still 62%. While the original Sobolev approximation ignores all possible absorption events beyond this point, there is a non-negligible absorption possibility beyond one Sobolev length and we therefore need to introduce a correction term in order to account for the additional matter. Even after three Sobolev lengths the relative overlap, for example, is still 13% and reabsorption might be possible.

Zoom In Zoom Out Reset image size

In order to find a proper correction for this overlap, we start from the basic definition of the column density. However, we are not interested in the total column density but rather in the effective column density, which includes only the gas that could be relevant for the reabsorption of escaping H₂ line photons. The Sobolev approximation gives a very simple answer to the question of which gas we have to include in the effective column density, namely all gas within one Sobolev length. Expressed in terms of line overlap we write this as

$$\begin{equation} N_{\mathrm{s}}= \underbrace{\int _0 ^{L_{\mathrm{s}}} 1 \cdot n_{\mathrm{H}_2\,} \mathrm{d} r}_{100\% \, \mathrm{overlap}} + \underbrace{\int _{L_{\mathrm{s}}} ^{\infty } 0 \cdot n_{\mathrm{H}_2\,} \mathrm{d} r}_{0\% \, \mathrm{overlap}} = n_{\mathrm{H}_2\,} L_{\mathrm{s}}, \end{equation} \tag{ 15 }$$

where a constant number density of molecular hydrogen is assumed in the last step. Thus, Sobolev (1947) overestimated the overlap and the matter contribution within one Sobolev length but neglected all matter contributions beyond this point. Since we want to account for the true overlap of spectral lines, we introduce the relative line overlap o(x) as a weighting function into the determination of the effective column density,

$$\begin{equation} N_{\mathrm{s,corr}} = \int _0 ^{\infty } o(x) \cdot n_{\mathrm{H}_2\,} \mathrm{d} r. \end{equation} \tag{ 16 }$$

The remaining problem is to relate the line displacement x to the distance r along the line of sight. Assuming a constant velocity gradient dv_r/dr, we can rewrite Equation (10) into

$$\begin{equation} \mathrm{d} r = \frac{L_{\mathrm{s}}}{v_{\mathrm{th}}} \mathrm{d} v_{\mathrm{r}} \end{equation} \tag{ 17 }$$

and transform the radial integration into an integration over the relative line displacement

$$\begin{equation} N_{\mathrm{s,corr}} = n_{\mathrm{H}_2\,} L_{\mathrm{s}}\int _0 ^{\infty } o(x) \mathrm{d} x, \end{equation} \tag{ 18 }$$

which finally yields

$$\begin{equation} N_{\mathrm{s,corr}} \simeq 1.694 N_{\mathrm{sob}}. \end{equation} \tag{ 19 }$$

In other words, the relevant column density for the determination of the escape probability of line photons is about 1.7 times higher than originally assumed by Sobolev. Although we correct for the line overlap, we should still keep in mind that this derivation implies several assumptions like constant velocity gradient and constant density within the Sobolev length. Also some gas further away may be absorbing, because of similar velocities. In particular in turbulent clouds this may become relevant. This is related to the question of how "real" clumps identified in position–position–velocity space are in position–position–position space (see, e.g., Ballesteros-Paredes & Mac Low 2002; Beaumont et al. 2013).

Molecular hydrogen has more than two hundred spectral lines in the range 1.0 μm ⩽ λ ⩽ 32 μm which could be relevant for cooling primordial gas (Ripamonti et al. 2002). Although these lines appear to be very close to each other, the possibility that a photon is emitted in one line with a certain frequency and absorbed by another line which is Doppler-shifted into the emitting frequency is negligibly small (∼10⁻⁶). Consequently, this effect can be ignored.

Besides the commonly used Sobolev approximation there are several other approaches, which we present in the following sections.

2.4.1. Gnedin Approximation

Gnedin et al. (2009) model molecular hydrogen and star formation in cosmological simulations. They determine H₂ column densities for the self-shielding of H₂ against Lyman–Werner photons by using a "Sobolev-like" approximation. They define a characteristic length scale

$$\begin{equation} L_{\mathrm{g}}= \frac{n_{\mathrm{H}_2\,}}{|\mathbf {\nabla } n_{\mathrm{H}_2\,}|} \end{equation} \tag{ 20 }$$

based on the number density and its gradient. The column density is then simply given by

$$\begin{equation} N_{\mathrm{g}} = n_{\mathrm{H}_2\,} L_{\mathrm{g}}. \end{equation} \tag{ 21 }$$

This approach accounts for the density gradient in the gas distribution. In the following, we will label this method as "Gnedin." They claim that this approximation provides a very good estimate for the column density in the range $3 \times 10^{20} \mathrm{cm}^{-2} < N_{\mathrm{HI}}+2N_{\mathrm{H}_2\,} < 3 \times 10^{23}\, \mathrm{cm}^{-2}$. Although we are dealing with star formation on much smaller scales than they do, their method seems to be a reasonable, scale-invariant approach for the determination of effective column densities. Furthermore, this method depends only on local quantities and is easy to implement in different codes. Nevertheless, this method includes no information about the velocity profile of the cloud and therefore neglects the enhanced photon escape probability due to the Doppler-shifting of lines.

2.4.2. Analytic Fit Functions

We are interested in the column densities in order to determine the photon escape probabilities for H₂ line cooling (Equation (8)). Besides the previously presented methods, there are two analytical fitting functions that directly relate a given number density of the gas to the photon escape probability. The first fitting function

$$\begin{equation} \beta = \min\, [ 1, (n/n_\mathrm{RA})^{-b_{\mathrm{RA}}} ] \end{equation} \tag{ 22 }$$

with n_RA = 8 × 10⁹ cm⁻³ and b_RA = 0.45 was proposed by Ripamonti & Abel (2004) and has been applied in several (mainly grid-based) simulations (O'Shea & Norman 2006; Turk et al. 2011; Hirano & Yoshida 2013; Greif et al. 2014a). This formula was obtained from the detailed one-dimensional calculation by Ripamonti et al. (2002). A second method was proposed by Greif et al. (2013) who study the chemo-thermal instability in primordial star-forming clouds. The idea follows the approach by Ripamonti & Abel (2004) but with a smooth transition and therefore a continuous derivative toward the optically thin regime. The formula is given by

$$\begin{equation} \beta =\left\lbrace \displaystyle\begin{array}{@{}l@{\quad }l@{}}\displaystyle\frac{(1+b_{\mathrm{G}})x}{x^{(1+b_{\mathrm{G}})}+b_{\mathrm{G}}} &\mathrm{for}\, x \ge 1 \\ 1 &\mathrm{for}\, x < 1 \end{array}\right. \end{equation} \tag{ 23 }$$

where x = n/n_G, n_G = 4 × 10⁹ cm⁻³, and b_G = 0.45. These fits are the easiest and most direct way to determine the photon escape probability, but their use does not account for any information about the temperature, velocity or density profiles. Recently, Greif et al. (2014a) investigated the collapse of primordial gas using a multi-line, multi-frequency ray-tracing scheme in order to accurately model the transfer of H₂ line emission. Although he provides a new fit formula, we will not include this in our analysis, since it is very similar to the original fit function by Ripamonti & Abel (2004).

2.4.3. Reciprocal Approach

In order to find a method for the determination of optically thick cooling that only depends on local quantities of the collapse and that is easy to implement, but nevertheless captures the dynamics of the cloud, we propose a combination of the Sobolev and Gnedin approximation. The corrected Sobolev approximation takes the line overlap into account but neglects the decreasing density. On the other hand, the Gnedin approximation takes the decreasing number density into account but neglects the Doppler-shifting of lines. Since each method on its own generally overestimates column densities, the general idea behind this new approach is to combine these two methods in order to overcome their individual shortcomings. The reciprocal sum

$$\begin{equation} \frac{1}{L_{\mathrm{rec}}} = \frac{1}{L_{\mathrm{g}}} + \frac{1}{L_{\mathrm{s,corr}}} \end{equation} \tag{ 24 }$$

of the two characteristic lengths provides the right behavior. If both lengths are of the same order (L₁  ≃  L₂), the result should be smaller than both lengths (L_rec  <  L₁  ≃  L₂), since each length individually overestimates the column density. Whereas if one length is significantly smaller than the other length L₁  ≪  L₂, the result should be equal to the smaller one (L_rec ≃ L₁), because beyond this smaller distance the photons can escape freely anyway. Following this method, the number density is simply given by

$$\begin{equation} N_{\mathrm{rec}} = n_{\mathrm{H}_2\,} L_{\mathrm{rec}}. \end{equation} \tag{ 25 }$$

The main problem of the Sobolev approximation is the dependence on local quantities and therefore the neglect of all information about density gradients, velocity profile and the actual shape of the cloud. Generally, the most exact way to determine effective column densities is to sum up all relevant mass along all possible lines of sight. In order to avoid this extremely high computational effort, Clark et al. (2012) designed the TreeCol algorithm, which determines column densities based on a tree structure, used by many gravitational N-body solvers. TreeCol uses a spherical pixelation with diamond-shaped pixels based on HEALPix (Górski et al. 2005). During the walk of the tree, all relevant data for the column density map are collected and projected onto a spherical grid. Since the data are already stored in a tree, TreeCol can use this information and therefore scales as Nlog N with the number of cells or particles N. However, in our simulations, the usage of TreeCol slows down the simulation by a factor of about five with respect to a run without TreeCol. This slowdown is mainly related to the evaluation of several inverse trigonometric functions. Although we use this method in an smoothed particle hydrodynamics (SPH)-based simulation, the only requirement for its implementation is the clustering of matter in a tree-like structure, as indeed, the TreeCol method has already been implemented in the Arepo moving mesh code (Smith et al. 2014) and the FLASH AMR code (R. Wünsch et al., in preparation).

TreeCol overcomes several shortcomings of the Sobolev approximation. First of all, we can use it for any density distribution because we directly sum up the individual mass contributions. Furthermore, we can use the velocity information of the tree nodes and do not have to assume a constant velocity gradient. Additionally, we account for the actual spatial matter distribution and do not have to stick to a rough one dimensional approximation. While the original TreeCol algorithm computes the H₂ column densities of the entire cloud, in what follows, we will present three improved versions that account for the Doppler-shifting of the line by only including mass that lies within the appropriate velocity range.

3.1.1. Sobolev-like

The Sobolev approximation assumes that all relevant mass for the column density is located within one Sobolev length. Translated into velocities, we should only include particles or tree nodes whose relative velocity is smaller than the thermal velocity. Thus, we modify TreeCol in order to include only the mass contributions of nodes that fulfill this criterion. A schematic illustration of this approach is given in Figure 2. From there we already see that the relevant volume around a particle does not necessarily have to be spherical, but usually follows the dynamic shape of the cloud.

Zoom In Zoom Out Reset image size

3.1.2. Corrected Sobolev

3.1.3. Lookup Method

In this section, we describe the general properties of our code, discuss the implementation of the TreeCol method, comment on its computational effort, and present our initial conditions.

3.2.1. Methods

We use a modified version of the cosmological simulation code Gadget-2. The original code was written by Springel (2005) in order to simulate structure formation by means of SPH. For the nonlinear collapse, we follow the chemical and thermal evolution of the gas cloud, solve the rate equations and couple the relevant heating and cooling terms to the hydrodynamic equations. For this purpose, we use a chemical network and cooling functions developed by Glover and collaborators (see Glover & Jappsen 2007; Glover & Abel 2008; Clark et al. 2011b and Clark et al. 2011a for full details). The primordial chemistry network includes H, D, He, H₂, H⁺, H⁻, D⁺, H$_2^+$, HD, He⁺, He⁺⁺ and e⁻. The rate equations between these species are solved self-consistently for every time step. For the three-body H₂ formation reaction

$$\begin{equation} \mathrm{H} + \mathrm{H} + \mathrm{H} \rightarrow \mathrm{H}_2\,+ \mathrm{H}, \end{equation} \tag{ 26 }$$

we adopt the rate coefficient of Glover (2008), derived by applying the principle of thermal balance to the collisional dissociation rate of Martin et al. (1996). The time scales associated with chemical heating and cooling can become extremely short when species come close to chemical equilibrium. Since the particle time step depends on the thermal timescale, this can result in some particles having very short timescales once the gas density exceeds ∼10¹⁷ cm⁻³. At this point it can become computationally prohibitive to run the simulation further.

Besides H₂ line cooling, which is the most relevant process for this analysis, we also track the contribution by several other heating and cooling processes, such as electronic excitation of H and He, recombination, photodissociation, HD line cooling, Compton cooling, and bremsstrahlung (see Ripamonti & Abel 2004; Glover & Jappsen 2007; Clark et al. 2011a for a more detailed description). For the collisional induced emission (CIE) cooling, we follow the implementation by Clark et al. (2011b), which is based on studies by Ripamonti et al. (2002), Ripamonti & Abel (2004), and Yoshida et al. (2008). Above densities of n_CIE ≈ 10¹⁴ cm⁻³, CIE cooling from H₂ becomes more efficient than H₂ line cooling, but the gas is most susceptible to fragmentation at densities below n_CIE (see Section 6.2). Consequently, the implementation of CIE cooling is unlikely to significantly affect the conclusions of our study.

Furthermore, we implement a simple feedback model for the accretion heating by assuming that the accreting protostars produce an accretion luminosity of $L_\mathrm{acc} = G M_{*} \dot{M}/R_{*}$, where $\dot{M}$ is the mass accretion rate and M_* and R_* are the mass and radius of the protostars, respectively. We use the models by Smith et al. (2011) for the protostellar radius and assume a constant accretion rate of $\dot{M}\,{=}\,10^{-2}\, M_\odot\, \mathrm{yr}^{-1}$, which is consistent with Clark et al. (2011b), Hirano & Yoshida (2013), and with the actual accretion rates in our model. The heating rate is then given by

$$\begin{equation} \Gamma _\mathrm{acc} = \rho \kappa _P \left(\frac{L_\mathrm{acc}}{4 \pi r^2} \right), \end{equation} \tag{ 27 }$$

where ρ is the gas density, κ_P is the Planck mean opacity and r is the distance from the source.

In order to minimize computational effort in the high-density regime, we use sink particles, based on the implementation by Jappsen et al. (2005), which was already used before (e.g., Boss & Black 1982; Boss 1987, 1989; Bate et al. 1995). Above a certain density threshold n_crit all SPH particles are merged into one single sink particle, which now contains all mass and momenta of the merged particles. This approach conserves mass and momentum, avoids small dynamical and chemical time steps, and we can identify the protostars by the newly formed sink particles. The critical density threshold is chosen so that the local Jeans mass is resolved by at least 100 SPH particles, following the resolution criterion by Bate & Burkert (1997).

Nevertheless, one should keep in mind that sink particles are not physical entities in their own right, but rather computationally motivated and consequently may cause problems by introducing a discontinuity in the mass and number of particles as well as a lack of pressure forces at the accretion radius. Strictly speaking, they also violate the hydrodynamic equations because the accretion onto the sink particle happens instantaneously (Greif et al. 2012). In order to guarantee that the formation of a sink particle actually represents the local collapse to a protostar, we introduce several formation criteria. The gas clump converted into a sink particle must exceed a certain density threshold, it must be gravitationally unstable, at the center of a locally convergent flow, and it must have a certain distance to preexisting sink particles (see, e.g., Federrath et al. 2010).

Another crucial value is the accretion radius r_acc, because a too large accretion radius might artificially influence the fragmentation behavior (Clark et al. 2011b; Greif et al. 2011; Machida & Doi 2013). On the other hand, Greif et al. (2012) and Machida & Doi (2013) state that fragmentation does not occur for densities above n ≃ 10¹⁷ cm⁻³, because there are no more efficient chemical or radiative cooling mechanisms. Smith et al. (2011), who additionally include the effect of heating by accretion feedback, find this value to be n ≃ 10¹⁵ cm⁻³. Phrased differently, the choice of the accretion radius might have an influence on the fragmentation behavior but as long as the critical density is n_crit  ≳  10¹⁶ cm⁻³, we should capture all fragmentation of the cloud. We set the critical density to n_crit  =  n_res  =  3.42 × 10¹⁵ cm⁻³ in order to resolve the Jeans mass throughout the simulation. The Jeans length is λ_J = 0.3AU under these conditions. Such a small accretion radius might lead to tiny dynamical time steps in the vicinity of sink particles and hence increase the computational effort. Therefore, we set the accretion radius to r_acc = 10 AU, which clearly fulfills the resolution criterion, but may suppress small-scale fragmentation.

3.2.2. Implementation of Cooling Approaches

The general idea of the different cooling approaches has already been discussed in Section 3.1. Here, we address the actual implementation of these approaches into the code. The optically thick cooling rate $\Lambda _{\mathrm{H}_2\,,\mathrm{thick}}$ is given by Equation (8), where the photon escape probability β can be expressed as a function of the column density divided by the thermal velocity of the H₂ molecules (N/v_th) and the temperature T. These values are stored in a lookup table for 31.6 K ⩽ T ⩽ 31, 600 K and 10¹⁷s  cm⁻³ ⩽ N/v_th ⩽ 10²⁷s  cm⁻³.

For the TreeCol-based determination of the effective column densities we create another lookup table, which relates the relative velocities in units of the thermal velocities (v_r/v_th) to an overlap of spectral lines. Since the relative velocities are distributed roughly equally throughout the simulation, we create this lookup table with linear steps in velocity space. For each node that might contribute to the effective column density we first check if

$$\begin{equation} \left(\frac{v_\mathrm{r}}{v_{\mathrm{th}}} \right) ^2 \le 43.3 \end{equation} \tag{ 28 }$$

because otherwise the line overlap is smaller than 10⁻³ and can be neglected anyway. The computational effort for the lookup of relative overlaps is negligibly small compared to the computational cost of TreeCol itself.

3.2.3. Initial Conditions

Before we apply our new method to cosmological halos, we test its accuracy in a simple test scenario. We start from a spherical, primordial cloud with a mass of 10³M_☉, an initial density of ∼10⁴ cm⁻³ and a temperature of ∼250 K. The gas is represented by 64³ SPH particles and we follow the simulation to densities above ∼10¹³ cm⁻³, where we insert sink particles. For each different method, we run the simulation independently and compare the effective H₂ column densities in Figure 3. The lookup approach explicitly accounts for the true line overlaps and therefore is the most exact method. While "TreeCol corrected Sobolev" accounts for the correction factor and reproduces "TreeCol lookup" remarkably well, "TreeCol Sobolev" corresponds to the uncorrected Sobolev method and generally underestimates the effective column density. Consequently, this correction factor is mandatory to account for the exact line overlaps.

Zoom In Zoom Out Reset image size

Furthermore, the commonly used Sobolev method and the Gnedin approach, which are only based on local gas quantities, overestimate the column density in the optically thick regime above ∼10¹⁰ cm⁻³. These two local methods yield almost the same results for the effective column density, which means that the relation $v_\mathrm{th}/|\mathbf {\nabla }\cdot {\mathbf {v}}| \approx n_{\mathrm{H}_2\,}/|\mathbf {\nabla } n_{\mathrm{H}_2\,}|$ holds for the spherical collapse scenario. This illustrates that the density is generally not constant within the Sobolev length, but rather varies on the same length scale as the velocity field. Hence, even in this simple test case, the local methods are not able to determine the column density properly.

In this section, we test our approaches in more realistic cosmological halos. We first use the TreeCol-based simulations and determine all relevant information by post-processing these output files. By doing so, we can focus on the intrinsic differences of the methods (determined under the same physical conditions) rather than comparing different simulations with presumably different dynamics. The plots in this section represent averages over all four cosmological halos.

We want to find an accurate cooling approach that reproduces the TreeCol method best and hence, we compare the commonly used Sobolev approximation with the corrected Sobolev approximation, the Gnedin approach and our new reciprocal method, which combines the corrected Sobolev and Gnedin approach. The column density for these different approaches can be seen in Figure 4. Generally, the local methods overestimate the effective column density in the optically thick regime (n ≳ 10¹⁰ cm⁻³). Especially at later stages of the collapse these differences increase, because the slope of the photon escape fraction as a function of density flattens with time in the optically thick regime for the TreeCol method. The (corrected) Sobolev method overestimates the column density all the time, whereas the Gnedin and reciprocal approaches overestimate the column density only for high densities. The latter yield accurate fits for the density regime 10⁹ cm⁻³ ≲ n ≲ 10¹² cm⁻³, whereby one should keep in mind that the relative importance of H₂ line cooling decreases above n ≃ 10¹³ cm⁻³. Realizing that the Sobolev approximation already overestimates column densities, one might ask why we need this additional correction factor, which makes the approximation even worse. A detailed answer to this question is given in Section 6.1, but we already want to emphasize the total neglect of any density gradient. The Sobolev approximation assumes a constant density, although the density of molecular hydrogen generally decreases when moving radially outward. Hence, this approximation is not valid, but leads to an overestimation of the column density (already for the uncorrected Sobolev method).

Zoom In Zoom Out Reset image size

The photon escape probability as a function of density can be seen in Figure 5 for different methods. Although there are slight differences between the individual approaches, all methods seem to agree well with the TreeCol approach prior to the formation of the first sink particles. However, at later times, the TreeCol approach yields higher values for the photon escape probability, corresponding to a smaller effective opacity of the cloud. The other methods, which depend only on local quantities, cannot reproduce this behavior and therefore underestimate the photon escape probability.

Zoom In Zoom Out Reset image size

We compare the analytical fitting formulae to the TreeCol method by using the parameterized formulae (Equations (22) and 23)) with the fit parameters n_RA, b_RA, n_G, and b_G. In Figure 6 we see the original functions and the best fits to the TreeCol data. The newly adjusted fits minimize the weighted scatter sum

$$\begin{equation} \sum _{i=1} ^{N_\mathrm{SPH}} \frac{|\beta _\mathrm{analytic} (i) - \beta _\mathrm{TreeCol} (i)|}{\beta _\mathrm{TreeCol} (i)} \frac{\Lambda _{\mathrm{H}_2} (i)}{\Lambda _{\mathrm{H}_2} (i) + \Lambda _\mathrm{CIE}(i)}, \end{equation} \tag{ 29 }$$

where N_SPH is the number of particles, β_analytic(i) is the analytic fit, β_TreeCol(i) is the photon escape probability based on the TreeCol approach, $\Lambda _{\mathrm{H}_2} (i)$ is the H₂ cooling rate, and Λ_CIE(i) is the CIE cooling rate for the ith particle, respectively. The weighting by the relative cooling rate accounts for the decreasing relevance of H₂ cooling at higher densities.

Zoom In Zoom Out Reset image size

For each single snapshot, the exact data can be fitted remarkable well by an analytic formula, but we note that both fit parameters vary strongly during the collapse. Although the original fits are a satisfying approximation at the beginning of the collapse, their accuracy decreases during the collapse and they totally miss the true photon escape probabilities at later stages.

The time evolution of the fit parameter b (Equations (22) and 23)) reveals an interesting insight into the structure of the collapse. In Figure 7 we compare the evolution of b for both fit formulae with time. It quantifies the slope and therefore represents the dependence of the opacity on density. The parameters b and n₀ are fitted simultaneously, but the plot only contains the slope because it reveals more information about the underlying physics. For large b, the cloud becomes opaque with increasing density very promptly, whereas for a shallow slope the opacity remains low although the density increases. The blue dash-dotted line illustrates the original slope for both fits of b = 0.45, but the newly fitted formulae to the TreeCol data reveal that the actual slope is shallower and most notably varies with time. The small value of b for t < −500 yr is a numerical artifact, because at these times there are almost no particles in the optically thick regime, which might define a distinct slope. On the other hand, the shallow slope of 0.2 ⩽ b ⩽ 0.4 during the later stages of the collapse can be related to the flattening of the cloud. Consequently, it is impossible to find one single analytic fit, which describes the dynamics of the collapse completely. However, the mean parameters that fit the TreeCol data best from 500 yr before formation of the first protostar until the end of our simulations are n_G = 5 × 10⁹ cm⁻³ and b_G = 0.32.

Zoom In Zoom Out Reset image size

In order to analyze the accuracy of the different cooling approaches quantitatively, we determine the relative error of the photon escape probability for each method. Therefore, we compare the photon escape probabilities of all particles above a certain density threshold to the ones determined with the TreeCol method. The density threshold is necessary, because below n = 10⁹ cm⁻³, the escape probability is very close to one anyway and consequently there are no significant deviations between the methods. The relative error |β_method − β_TreeCol|/β_TreeCol is weighted by the relative H₂ cooling rate, so that the mean relative error of the photon escape probability can be expressed as

$$\begin{eqnarray} \overline{\Delta \beta } &=& \frac{1}{N_\mathrm{thresh}} \sum _{i=1} ^{N_\mathrm{thresh}} \frac{|\beta _\mathrm{method} (i) - \beta _\mathrm{TreeCol} (i)|}{\beta _\mathrm{TreeCol} (i)} \nonumber \\ &&\times \frac{\Lambda _{\mathrm{H}_2} (i)}{\Lambda _{\mathrm{H}_2} (i) + \Lambda _\mathrm{CIE}(i)}, \end{eqnarray} \tag{ 30 }$$

where N_thresh is the number of particles above a certain density threshold. We chose the density threshold 10⁹ cm⁻³ to analyze how accurate the methods are in this density regime. The time evolution of the mean relative photon escape probability error for this threshold can be seen in Figure 8. The relative error is small at early times in the collapse because all photon escape probabilities are close to one, whereas at later times, more particles enter the optically thick regime and the relative errors rise to values between 5%  and 30%. The reciprocal approach is slightly better than the Gnedin approximation, although both yield errors between 5%  and 15% throughout the simulations. Interestingly, the corrected Sobolev approximation yields higher errors (15%–30%), whereas the uncorrected one yields errors between 10%  and 25%. This behavior is caused by the fact that the Sobolev approximation already overestimates column densities because it neglects the density gradient. Since the reciprocal method is mostly dominated by the Gnedin approximation, we can conclude that the effective column density is more strongly influenced by the density gradient than by the relative velocities. However, the velocity information still improves the fit and is therefore important for a proper treatment. The original analytic fits are comparatively inaccurate with errors between 20%  and 50% and are therefore not displayed in this figure. However, the analytical fitting function proposed by Greif et al. (2013) yields slightly better results than the one proposed by Ripamonti & Abel (2004). The same analysis for higher density thresholds (10¹⁰ cm⁻³, 10¹¹ cm⁻³) shows qualitatively similar results, but with a trend to higher mean relative errors.

Zoom In Zoom Out Reset image size

We are interested in how the treatment of the cooling affects the characteristic mass range for Pop III stars. Based on our simulations, we compare the number of sink particles in each halo at the end of the simulations. However, due to a prohibitively small chemical time step and a too simplified treatment of feedback, we end the simulations before accretion has terminated. In order to compare the number of sink particles, we choose for each halo the snapshots with the same amount of accreted mass onto the sink particles. This condition enables us to compare the clouds for the Sobolev-based and the TreeCol-based runs at the same stage of the accretion process. The resulting numbers of sink particles are given in Table 1.

The TreeCol approach yields more sink particles in all but one halo. However, the time evolution of the number of sink particles (see Figure 9) suggests that this effect might be caused by a delay of the fragmentation in the Sobolev-based runs. While the Sobolev approximation causes the disk to fragment into 1–6 protostars, the TreeCol approach causes slightly earlier fragmentation with at least two protostars in each halo. Since we had to end the simulations at this point, the final number of protostars cannot be fully constrained in this case. However, the TreeCol method accelerates the collapse and the gas fragments earlier compared to the Sobolev approximation. The accretion rates onto individuals protostars are in the range $3.3 \times 10^{-4}\, M_\odot\, \mathrm{yr}^{-1} \le \dot{M} \le 1.2 \times 10^{-2}\, M_\odot\, \mathrm{yr}^{-1}$ for the Sobolev approximation and $1.2 \times 10^{-3}\, M_\odot\, \mathrm{yr}^{-1} \le \dot{M} \le 1.7 \times 10^{-2}\, M_\odot\, \mathrm{yr}^{-1}$ for the TreeCol approach in all four halos.

Zoom In Zoom Out Reset image size

In order to study the fragmentation behavior independently of the number of sinks, we first analyze the dynamics of the collapse. This analysis is mainly based on the fact that the temperature profile of a gas cloud directly influence the collapse to protostars by affecting the local Jeans mass, which is the absolute minimum requirement for collapse and is a strong function of temperature (Clark et al. 2011b). Hence, we compare the basic quantities between the four runs using the TreeCol approach and the four runs using the Sobolev approximation. The first quantity to compare is naturally the photon escape probability. Figure 10 illustrates the time evolution of the photon escape probability for both methods. While the spatial distribution of the photon escape probability seems to be smoother for the TreeCol-based run, the distribution is comparatively structured for the Sobolev approximation. Additionally, we see another important difference. Using the TreeCol method yields values for the photon escape probability close to one at the end of the simulations, whereas the Sobolev approximation yields smaller photon escape probabilities. Obviously, the Sobolev approximation is not able to capture the flattening of the cloud (see Section 6.1). Furthermore, the Sobolev-based simulation develops only one central core, while an elongated core with two peaks is formed in the TreeCol-based run. Since the photon escape probability has a direct influence on the cooling rates, we are consequently interested in the H₂ cooling rate. This cooling rate is not significantly higher for the TreeCol-based run, although the photon escape probabilities are higher. This can be explained by the cooling implementation: due to the presumably higher cooling rate in the TreeCol-based run, thermal equilibrium is reached for lower temperatures. Thus, the cloud simply remains cooler instead of increasing its cooling rate significantly. Consequently, the gas in the Sobolev-based runs is generally hotter, which can also be seen in Figure 11.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The difference of temperatures in this inner regime can be up to ∼500 K, which influences the fragmentation behavior significantly.

Finally, we compare different fragmentation criteria of the clouds. The possibility of fragmentation in primordial clouds has been considered by many authors (e.g., Sabano & Yoshii 1977; Clark et al. 2008, 2011a, 2011b; Greif et al. 2011, 2012, 2013; Bromm 2013; Stacy & Bromm 2014). Fragmentation is a very chaotic, non-deterministic process and the actual outcome depends sensitively on the initial conditions (Girichidis et al. 2011, 2012a, 2012b). Nevertheless, there are three analytic expressions that help to quantify the possibility of a gas cloud fragmenting.

• 1.  
  In order to locally contract instead of globally collapse, a necessary criterion for fragmentation of a gas cloud is (Rees & Ostriker 1977)

  $$\begin{eqnarray} &&\frac{t_\mathrm{cool}}{t_\mathrm{ff}} < 1, \end{eqnarray} \tag{ 31 }$$

  where t_cool = nk_BT/(Λ[γ − 1]) is the cooling time with the net cooling rate Λ and $t_\mathrm{ff} = \sqrt{3 \pi /(32 G \rho)}$ is the free-fall time.
• 2.  
  Toomre (1964) analyzed the stability of rotating gas disks and derives the instability criteria

  $$\begin{eqnarray} &&Q=\frac{c_s \kappa }{\pi G \Sigma } <1, \end{eqnarray} \tag{ 32 }$$

  where κ is the epicyclic frequency and Σ is the surface density of the disk. Formally, this criterion is only valid for thin disks whereas Goldreich & Lynden-Bell (1965) extended the criterion by requiring Q < 0.676 for a finite-thickness isothermal disk to fragment.
• 3.  
  Gammie (2001) investigated the nonlinear outcome of a stability analysis of a Keplerian accretion disk. Based on numerical experiments, he derived the instability criterion

  $$\begin{eqnarray} &&\frac{t_\mathrm{cool}}{3 \omega ^{-1}} < 1, \end{eqnarray} \tag{ 33 }$$

  where ω is the orbital frequency. The Gammie criterion expresses the possibility that pieces of the disk cool and collapse before they have the opportunity to collide with one another in order to reheat the disk.

In their study on the formation and evolution of primordial protostellar systems, Greif et al. (2012) find the Toomre criterion insufficient for the quantification of gravitational instability and additionally use the Gammie criterion. None of these criteria guarantees fragmentation (Yoshida et al. 2006), but a combination of them yields a reliable quantification of instabilities in the gas disk.

We compare the radial and density profile of these fragmentation criteria for the Sobolev-based and the TreeCol-based runs. We do this directly after formation of the first protostar to ensure comparability between the different runs and methods (Figure 12). The requirement that all three fragmentation criteria have to be fulfilled simultaneously is very strict, because each criterion alone already quantifies stability. Applying this conservative criterion shows that all simulations are susceptible to fragmentations in the density regime 10¹⁰ cm⁻³  ⩽  n ⩽ 10¹¹ cm⁻³ at a radius of ∼200 AU from the densest region. However, the Gammie criterion seems to be a less restrictive version of the Rees–Ostriker criterion and consequently yields no additional constraint to this regime, because it is fulfilled anyway. The regime of fragmentation and the instability of the gas cloud is slightly higher for the TreeCol approach. Especially the cooling time criterion remains significantly longer under the stability threshold for the TreeCol-based run, which is mainly a consequence of the lower temperatures, as discussed before. Consequently, the high-density gas is more susceptible to fragmentation in the TreeCol-based runs. Since these are averaged profiles over four halos, we expect these results to be representative for primordial star formation.

Zoom In Zoom Out Reset image size

The local and isotropic approximations fail in determining the proper effective column densities for H₂ cooling. Although they yield acceptable results up to the formation of the first protostars, most approaches break down during later phases of the collapse, because they cannot capture the evolving dynamical structure, particularly once a flattened accretion disk has formed. The appropriate consideration of the overlap of spectral lines is very important, although even the uncorrected Sobolev approximation generally overestimates column densities. Further shortcomings are mostly related to the assumption that all relevant quantities are constant within the characteristic length scale. This simplification avoids the evaluation of integrals along the line of sight but is formally only valid for large velocity gradients. The analysis of the relevant quantities shows that the H₂ number density, thermal velocity and the gradient of the velocity vary by up to two orders of magnitude within one Sobolev length. While the thermal velocity seems to be rather constant, the H₂ number density varies most strongly. Since the distribution around the central value is not symmetric, these effects do not cancel out and generally lead to an overestimation of column densities.

Due to angular momentum conservation in the infalling material, a disk forms around the first protostar and the photon escape probability is enhanced perpendicular to the disk. An isotropic column density estimation that depends only on local quantities, cannot capture this feature and therefore generally underestimates the photon escape probability (see Figure 13).

Zoom In Zoom Out Reset image size

The Sobolev approximation is not able to capture this angular dependence and consequently underpredicts the angle-averaged photon escape probability. Comparing the escape probabilities in different directions for one and the same particle quantitatively yields high angle-dependent variations. Shortly before formation of the first sink particle, the photon escape probabilities in the inner ∼100 AU vary by factors of up to 10–100 in different HEALPix pixels. During the further collapse, this value decreases but is still between 2  and 20 on average in the center of the cloud. This clearly demonstrates the failure of the local approximation in modeling inhomogeneous, non-spherical density distributions. Hirano & Yoshida (2013) also mention the problem of direction-dependent escape probabilities during their comparison of different cooling implementations in simulations of primordial star formation. Moreover, the Sobolev approximation is highly sensitive to small-scale details of the velocity field, and hence improving the resolution and accuracy with which one models the sub-sonic turbulence in the gas actually makes the method worse, because locally you start to see more of the turbulent velocity gradients, but these persist only over short length scales, not globally (Greif et al. 2014a).

Although the TreeCol-based method is computationally expensive, it is able to capture the non-uniformity of the H₂ distribution in a way that local approximations cannot manage. It is also significantly cheaper than the more robust method employed by Greif et al. (2014a), although a comparison between the accuracy of the method proposed here and that by Greif et al. (2014a) still needs to be done.

Regarding the local methods, there is no general recommendation. Most of these methods assume spherical symmetry and hence their validity breaks down when a disk forms around the first protostar. Since the accretion process through the disk lasts for several thousand years or more, one should use an approach that yields a suitable long-time accuracy. The reciprocal method is the only approach that considers both gradients in density and the Doppler-shifting of spectral lines. Therefore, this method is the most accurate approach for the determination of the effective column density even after subsequent sink particle formation. The difference between whether one uses the corrected or uncorrected Sobolev length for the determination of the reciprocal sum is small, but the corrected Sobolev length provides a more accurate method for the long-time evolution of the disk. Analytic fits might be useful during the initial collapse of the primordial cloud, but also fail when the disk-like structure starts to form. Nevertheless, at any stage of the collapse, the escape probability as a function of mass can generally be fitted very accurately by an analytic formula, which could also be done on the fly during the simulation. Consequently, a hybrid version between the proper TreeCol-based determination of escape probabilities every several time steps and an analytic fit to these data for the steps in between should yield reliable results with less computational effort.

This new method for the determination of effective column densities can also be applied in other astrophysical scenarios. For example the photodissociation of H₂ in protogalaxies is a crucial process for the formation of supermassive black hole seeds (Latif et al. 2013). To model these processes properly, the effective column densities are needed to account for the effect of H₂ self-shielding against Lyman Werner photons (Wolcott-Green et al. 2011). It can generally be applied in most scenarios of line transfer, where the radially infalling velocities of gas are higher than the turbulent gas velocities.

Fragmentation is a highly chaotic process and slight changes in the initial conditions or in the implementation of the governing physics can completely change the outcome. Therefore, looking just at the number of sink particles is not a valid quantification, especially since we do not capture the entire duration of the accretion and fragmentation process. Nevertheless, it is instructive to determine and compare different fragmentation criteria for the individual cooling approaches. We find that the TreeCol method seems to promote fragmentation. In other words, commonly used cooling approximations generally underestimate the number of Pop III stars. Moreover, the disk has just formed, when we end our simulations. Therefore, the differences between the individual cooling approaches and their effect on fragmentation might even be more pronounced for later stages of the collapse. Since the overall mass accretion rates are roughly equal, regardless of the cooling implementation, we expect the Pop III stars in our TreeCol-based runs generally to have smaller masses than previous Sobolev-based studies have yielded.

There are several open questions, shortcomings and approximations, which one should keep in mind, when interpreting the previously presented results. Since we follow the fragmentation of the cloud only for the first few hundred years after formation of the first protostar, we have no information about the physical conditions during the late stages of disk accretion. It is, however, likely that the disk-like structure will proceed to grow and most of our statements remain valid.

Furthermore, we do not account for mergers, which influence the final number of Pop III stars (Greif et al. 2011, 2012; Stacy & Bromm 2014), nor do we include the effects of ionizing radiation. We determine radiative feedback under the assumption of a constant mass accretion rate and although the value of 10⁻² M_☉ yr⁻¹ seems to be justified by similar accretion rates in the simulations, a proper treatment is necessary. Besides, we should also keep in mind that the escape probability, as the theoretical basis for our cooling implementation, is an approximation by itself, e.g., we use one average escape probability for all photons, instead of determining the individual escape probabilities for each line separately. In this context, Greif et al. (2014a) find that a multi-line, multi-frequency ray-tracing scheme does not alter these results significantly. Magnetic fields might influence the collapse and act as a stabilizing force. We have not included these effects in our simulation, but for a detailed discussion of the effects of magnetic fields in primordial star formation, see Machida & Doi (2013 and references therein).