{'fraction_non_alphanumeric': 0.036483685960769195, 'fraction_numerical': 0.02578437515060967, 'mean_word_length': 4.2835212804656235, 'pattern_counts': {'":': 0, '<': 42, '<?xml version=': 0, '>': 42, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Is H3+ cooling ever important in primordial gas?

Studies of the formation of metal-free Population III stars usually focus primarily on the role played by H2 cooling, on account of its large chemical abundance relative to other possible molecular or ionic coolants. However, while H2 is generally the most important coolant at low gas densities, it is not an effective coolant at high gas densities, owing to the low critical density at which it reaches local thermodynamic equilibrium (LTE) and to the large opacities that develop in its emission lines. It is therefore possible that emission from other chemical species may play an important role in cooling high density primordial gas. A particularly interesting candidate is the H3+ molecular ion. This ion has an LTE cooling rate that is roughly a billion times larger than that of H2, and unlike other primordial molecular ions such as H2+ or HeH+, it is not easily removed from the gas by collisions with H or H2. It is already known to be an important coolant in at least one astrophysical context -- the upper atmospheres of gas giants -- but its role in the cooling of primordial gas has received little previous study. In this paper, we investigate the potential importance of H3+ cooling in primordial gas using a newly-developed H3+ cooling function and the most detailed model of primordial chemistry published to date. We show that although H3+ is, in most circumstances, the third most important coolant in dense primordial gas (after H2 and HD), it is nevertheless unimportant, as it contributes no more than a few percent of the total cooling. We also show that in gas irradiated by a sufficiently strong flux of cosmic rays or X-rays, H3+ can become the dominant coolant in the gas, although the size of the flux required renders this scenario unlikely to occur.


INTRODUCTION
Over the past decade, we have made substantial progress in understanding how the very first stars in the universe formed. We know that in cosmological models based on cold dark matter (CDM), the first stars will form in small protogalaxies, with total masses of the order of 10 5 -10 6 M⊙, and that by a redshift z ∼ 30 we expect to find at least one such star-forming system per comoving Mpc 3 (Yoshida et al. 2003). We also know that although molecular hydrogen formation is inefficient, owing to the absence of dust, it is nevertheless the most abundant molecule in primordial gas, and is the main source of cooling at low densities, for temperatures between ∼ 200 K and 8000 K.
Furthermore, simple semi-analytical estimates (e.g., Tegmark et al. 1997), later confirmed by detailed simulations (e.g., Yoshida et al. 2003), demonstrate that H2 provides enough cooling in these small protogalaxies to allow the gas to collapse under the influence of gravity on a timescale comparable to its gravitational free-fall timescale, thereby allowing star formation to occur. High-resolution, adaptive mesh simulations performed by Abel, Bryan, & Norman (2000 have taught us much about the dynamics of the gas in these first protogalaxies. They consider gas cooled only by H2, and find that gravitational fragmentation of the collapsing gas is inefficient and that therefore only a single, massive star will form during the collapse. This will then suppress further star formation through a variety of feedbacks (see the recent review by Ciardi & Ferrara 2005). Other simulations, making simplifying assumptions such as the adoption of spherical symmetry, have allowed us to examine the importance of various aspects of the physics of the gas which are currently difficult to treat in the high-resolution adaptive mesh simulations (for instance, the development of large optical depths in the rotational and vibrational lines of H2 at high densities; see e.g. Omukai & Nishi 1998;Ripamonti et al. 2002).
Much of the theoretical uncertainty that remains concerns the behaviour of dense gas during the later stages of gravitational collapse, and during the period of accretion that follows the formation of the first protostar. There is much that is still unknown here -for instance, a complete understanding of the mechanism by which the collapsing gas loses much of its initial angular momentum still eludes us -but in this paper we intend to focus on one relatively simple aspect: the identification of the dominant coolant(s) in the dense gas.
As previously noted, H2 cooling has long been known to dominate at low densities, and it is frequently assumed that it also dominates at high densities. However, it is not at all obvious that this is actually the case. Two factors dramatically reduce the effectiveness of H2 as a coolant in high density gas. The first is the fact that the excited rotational and vibrational levels of H2 have small radiative transition probabilities, and hence long radiative lifetimes (τ > ∼ 10 6 s; Wolniewicz, Simbotin, & Dalgarno 1998). This means that collisional de-excitation of excited H2 becomes competitive with radiative de-excitation at fairly low gas densities (n ∼ 10 4 cm −3 ), and so as the number density n increases, the cooling rate of H2 quickly reaches its local thermodynamic equilibrium (LTE) value, given by where Aji is the transition probability for a transition from j → i, Eji is the energy of this transition, nj is the number density of H2 molecules in level j, computed assuming LTE, and where we sum over all bound levels i and over all bound levels j with energies greater than i. In the LTE limit, the cooling rate per H2 molecule is independent of the gas density, and is largely determined by the size of the transition probabilities. Since these are small, the LTE cooling rate is also small. Consequently, at high gas densities, a molecular species whose excited states have much shorter radiative lifetimes than those of H2 will provide far more cooling per molecule than H2. The second factor reducing the effectiveness of H2 cooling at very large n is the fact that the gas eventually becomes optically thick in the cores of the main rovibrational lines of H2. The effects of this cannot currently be treated fully in threedimensional simulations, due to the high computational cost of solving the resulting radiative transfer problem, but it has been modeled accurately in simple spherically symmetric, one-dimensional simulations (Omukai & Nishi 1998;Ripamonti et al. 2002;Ripamonti & Abel 2004) and in an approximate fashion in three dimensions (Yoshida et al. 2006). These studies confirm that at densities n > ∼ 10 10 cm −3 , optical depth effects significantly suppress H2 cooling. Together, these factors combine to render H2 a fairly ineffective coolant in high density gas, despite the fact that at n > 10 10 cm −3 , several three-body H2 formation reactions H + H + H → H2 + H, H + H + He → H2 + He, H + H + H2 → H2 + H2, rapidly convert almost all of the hydrogen in the gas to H2 (Palla, Salpeter, & Stahler 1983). It is therefore reasonable to ask whether cooling from any of the other molecular species present in the gas will become competitive with H2 cooling at these densities. One obvious possibility is deuterated hydrogen, HD. Its excited rotational and vibrational levels have radiative lifetimes that are about a factor of 100 shorter than those of H2, and so the HD cooling rate does not reach its LTE limit until n ∼ 10 6 cm −3 . It is also a far more effective coolant than H2 at low temperatures (T < ∼ 200 K; see e.g., Flower et al. 2000). This is due primarily to the fact that radiative transitions can occur between rotational levels with odd and even values of J, allowing cooling to occur through the J = 1 → 0 transition. The corresponding odd ↔ even transitions in the case of H2 represent conversions from ortho-H2 to para-H2 or vice versa, and are highly forbidden. Furthermore, at low temperatures the ratio of HD to H2 can be significantly enhanced with respect to the cosmological D:H ratio by chemical fractionation (see e.g., Glover 2008).
The role of HD cooling in early protogalaxies has been investigated by a number of authors. In the case of the earliest generation of protogalaxies, which form from very cold neutral gas that is never heated to more than a few thousand K during the course of the galaxy formation process, HD cooling appears to be unimportant, as the collapsed gas does not become cold enough for sufficient fractionation to occur to make HD cooling dominant (Bromm, Coppi, & Larson 2002).
The situation is rather different, however, in primordial gas cooling from an initially ionized state. In that case more H2 is formed, allowing the gas to cool to a lower temperature, at which point fractionation becomes effective and HD cooling rapidly becomes dominant (see e.g., Nakamura Shchekinov & Vasiliev 2006). However, the initial ionization required is much larger than that expected to be present in the earliest protogalaxies.
Another molecule to have attracted considerable attention is lithium hydride, LiH. This molecule has a very large dipole moment, µ = 5.888 debyes (Zemke & Stwalley 1980), and consequently its excited levels have very short radiative lifetimes. Therefore, despite the very low lithium abundance in primordial gas (xLi = 4.3 × 10 −10 , by number; see Cyburt 2004), it was thought for a time that LiH would dominate the cooling at very high densities (see e.g., Lepp & Shull 1984). However, accurate quantal calculations of the rate of formation of LiH by radiative association Gianturco & Gori Giorgi 1996;Bennett et al. 2003) Li + H → LiH + γ, have shown that the rate is much smaller than was initially assumed, while recent work by Defazio et al. (2005) has shown that the reaction LiH + H → Li + H2, has no activation energy and so will be an efficient destruction mechanism for LiH for as long as some atomic hydrogen remains in the gas. Consequently, the amount of lithium hydride present in the gas is predicted to be very small, even at very high densities, and so LiH cooling is no longer believed to be important (Mizusawa, Omukai, & Nishi 2005). In contrast to HD or LiH, the various molecular ions present in the gas, such as H + 2 , H + 3 or HeH + , have attracted little attention. Some early work on H + 2 cooling in ionized primordial gas can be found in Suchkov & Shchekinov (1977, 1978, and its possible importance in hot, highly ionized conditions has recently been re-emphasized by Yoshida et al. (2007), but aside from this, there has been little exploration of the role that cooling from these species might play in the evolution of primordial gas, presumably because the abundances of these species are expected to be small. It is this absence that the current paper attempts to rectify.
We present here the results of a set of simulations of the chemical and thermal evolution of gravitationally collapsing primordial gas. These simulations use a very simple one-zone dynamical model for the gas, but couple this with a detailed chemical network and a comprehensive model of the various heating and cooling processes at work. Besides the coolants considered above, we include the effects of cooling from (in no particular order) H + 2 , HD + , D + 2 , H + 3 , H2D + , HD + 2 , D + 3 , HeH + , HeD + , He + 2 , LiH + , LiD + , LiD and LiH + 2 . We focus in particular on the possible role of H + 3 . This ion has a very large number of excited rotational and vibrational levels that are energetically accessible at the temperatures of interest in primordial gas, and its vibrational levels have much shorter radiative lifetimes than those of H2 or HD. In LTE, its cooling rate per molecule is roughly 10 9 times larger than that of H2. It is known to be an important coolant in planetary atmospheres (Miller et al. 2000) and may also be an effective coolant in high density primordial gas . Unlike ions such as H + 2 or HeH + , it does not react with H2, and is not easily destroyed by collisions with H, and so its abundance in high density gas should be large compared to the other molecular ions included in our model. The layout of this paper is as follows. In §2 we outline the numerical method used to simulate the thermal and chemical evolution of the protostellar gas. The chemical reactions included in the model are discussed in §3, and the thermal processes, in particular H + 3 cooling, are discussed in §4. We present the results of our simulations in §5 and close with a brief discussion in §6.

NUMERICAL METHOD
We treat the thermal and chemical evolution of the gas using a one-zone model, in which the density is assumed to evolve as where t ff = 3π/32Gρ is the free-fall timescale of the gas and where η is an adjustable constant. In most of our simulations, we set η = 1, corresponding to free-fall collapse, but in a few runs we examine the effect of slowing down the collapse by setting η < 1 (see §5.7 for details).
To follow the temperature evolution, we solve the energy equation where e is the internal energy density, p is the thermal pressure, Λ is the total cooling rate (which includes contributions from both radiative and chemical cooling, as outlined in §4) and Γ is the total heating rate. Since the temperature evolution is strongly coupled to the chemical evolution, we solve Equation 8 simultaneously with the chemical rate equations using the DVODE implicit ordinary differential equation solver (Brown, Byrne, & Hindmarsh 1989). To model the chemistry we use an extensive chemical network consisting of 392 reactions among 30 atomic and molecular species. Table 1 lists all 30 species considered. Tables A1-A14 in Appendix A list the reactions included in this network, broken down by the type of process involved. These tables also give details of the rate coefficient or rate adopted for each reaction and the source of the data. In these tables and elsewhere in the paper, T is the gas temperature in K, T3 = T /1000 K, and Te is the temperature in units The quoted values are fractional abundances relative to the number density of hydrogen nuclei. Chemical species listed without initial abundances are assumed to be in chemical equilibrium, as described in the text.
of eV. (Note that we assume, both here and throughout, that the kinetic temperature of the electrons is the same as that of the atoms and molecules). For many of these species, we followed the full time-dependent, non-equilibrium chemistry, but in some cases -the ions H − , D − , Li − , H + 2 , HD + , D + 2 , HeH + , HeD + , He + 2 , LiD + and LiH + 2 -chemical equilibrium is reached very rapidly, on a timescale of the order of teq = 1/(k dest n) ∼ 10 9 n −1 s, where k dest ∼ 10 −9 cm 3 s −1 is a characteristic destruction rate coefficient and n is the number density of hydrogen nuclei. At the gas densities considered in this work (1 < n < 3 × 10 13 cm −3 in the majority of our simulations), this chemical equilibrium timescale for these rapidly reacting ions is orders of magnitude smaller than the free-fall timescale of the gas, which is approximately 1.5 × 10 15 n −1/2 s −1 . It is therefore sufficient to use the equilibrium abundances for these ions. Some of the species remaining in our non-equilibrium model may also be close to chemical equilibrium during a large portion of the collapse, but are included in the non-equilibrium treatment because we cannot be sure that they always remain in equilibrium. Further details of our chemical network are given in §3.
We assume elemental abundances relative to hydrogen of 0.083, 2.6 × 10 −5 and 4.3 × 10 −10 for helium, deuterium and lithium, respectively (Cyburt 2004). In §5.5 we explore the effects of reducing the deuterium and/or the lithium abundance to zero, in order to asses the impact of the deuterium and lithium chemistry on the evolution of the gas. The initial abundances of the various molecular and ionic species in our standard model are summarized in Table 1. The values used for the initial H + , He + and H2 abundances, and the ratio of ionized to neutral lithium are taken from the calculations of  -specifically their model V. Deuterated species are assumed to have abundances that are a factor (2.6 × 10 −5 ) N D smaller than the abundances of the undeuterated equivalents, where ND is the number of deuterium nuclei in the species in question. The electron abundance is computed assuming charge conservation. To assess our sensitivity to these initial values, we have also run several simulations with different initial abundances; the results of these runs are discussed in §5.5.
For most simulations, we adopt an initial density ni = 1 cm −3 and an initial temperature Ti = 1000 K. The effects of altering ni and Ti are examined in §5.5. All of the simulations are run until the density exceeds n f = 3 × 10 13 cm −3 . At higher densities, collision-induced emission from H2 quickly comes to dominate the cooling, and the minor species considered here are unlikely to be important coolants in this very high density regime.
Finally, we ran all of our simulations starting at a redshift z = 20. However, the main influence of the redshift is to set a minimum temperature for the gas (since the gas cannot cool radiatively to below the CMB temperature, TCMB). As the gas temperature T ≫ TCMB at the densities of interest in this paper, we do not anticipate that changing z by a moderate amount will significantly affect our results.

CHEMISTRY
The chemical evolution of the gas is modeled using a chemical network consisting of 392 reactions amongst 30 neutral and ionized species. A list of all of the reactions included can be found in Tables A1-A14. This network, which to the best of our knowledge is the largest used to date for the study of primordial gas chemistry, is based in part on previous compilations by Abel et al. (1997), , Galli & Palla (1998, Lepp, Stancil, & Dalgarno (2002), Wang & Stancil (2002), and Walmsley, Flower, & Pineau des Forêts (2004), supplemented with additional reactions drawn directly from the chemical literature, as well as some whose rates have not (to our knowledge) been previously discussed. These latter are generally rates involving one or more deuterium nuclei in place of hydrogen nuclei; and in estimating rates for these reactions, we have generally followed the same procedure as in : for a non-deuterated reaction with a reaction rate coefficient that has a power-law temperature dependence k ∝ T m , we have estimated a rate coefficient for the deuterated reaction by multiplying this rate coefficient by a scaling factor (µH/µD) m , where µH and µD are the reduced masses of the reactants in the non-deuterated and deuterated reactions, respectively. Some notable exceptions to this strategy are discussed in more detail in section §3.1 below.
For reactions where the presence of a deuteron increases the number of distinguishable outcomes and where no good information exists on the branching ratio of the reaction, we have assumed that each outcome is equally likely. An example of this is the dissociative attachment of HD with e − (reactions AD6-AD7), which can produce either H and D − or H − and D, in contrast to the dissociative attachment of H2 with e − (reaction AD5) which can only produce H − and H. For this particular example, this assumption gives branching ratios of 50% for reactions AD6 and AD7, respectively.
In spite of the size of our chemical network, there remain a number of chemical processes that are not included. These are discussed in §3.4, along with our justifications for omitting them.  Table A2) We assume that the ionizing photons produced by the recombination of H + to ground state H are immediately consumed by the ionization of atomic hydrogen (the on-the-spot approximation), and so we use the case B rate coefficient for hydrogen recombination (Ferland et al. 1992). We note that although the fractional abundance of atomic hydrogen becomes small at densities greater than 10 10 cm −3 , the number density of atomic hydrogen remains considerable, and so case B remains a good approximation.
We also use the on-the-spot approximation to treat He + recombination, but in this case the net recombination rate is larger than the case B rate, as some of the photons produced by recombination directly into the n = 1 ground state are lost through photoionization of H rather than He. For He + recombinations directly into the ground state, occurring in primordial gas with a low fractional ionization and low molecular abundance, approximately 68% of the resulting photons are absorbed by H, with the remaining 32% being absorbed by He (Osterbrock 1989). Therefore, the effective He + recombination rate coefficient in these conditions is given by kPR3 = 0.68kPR3,rr,A + 0.32kPR3,rr,B + k PR3,di cm 3 s −1 , where kPR3,rr,A and kPR3,rr,B are the case A and case B rate coefficients, and k PR3,di is the dielectronic recombination rate coefficient. This formula becomes incorrect once the H2 fraction of the gas becomes large, but as this occurs only at densities n > ∼ 10 10 cm −3 at which the He + fractional abundance is negligible, the error that is introduced by using this prescription for kPR3 throughout the simulation is unimportant.
It is also necessary to take account of the photoionization of H caused by the He + recombination emission. In addition to the contribution coming from He + recombination into the n = 1 ground state, there is an additional contribution made by photons produced during transitions from n = 2 to n = 1 in He; in other words, even pure case B recombination of He + produces H-ionizing photons. The proportion of case B recombinations that yield photons capable of ionizing hydrogen depends upon the relative populations of the n = 2 singlet and triplet states, and hence upon the electron density, but in the low density limit, 96% of all recombinations to excited states produce photons that will ionize hydrogen (Osterbrock 1989). To model the effects of these photons, along with those produced by recombination direct to the n = 1 ground-state and by dielectronic recombination, we include in our chemical network a local H ionization rate per unit volume Rpi, with a value Rpi = [0.68(kPR3,rr,A − kPR3,rr,B) + 0.96kPR3,rr,B + 2k PR3,di ] nen He + cm −3 s −1 , where the three terms in brackets on the first line correspond to the contributions from recombination direct to the n = 1 ground state, pure case B recombination (i.e., recombination to all states n 2), and dielectronic recombination, respectively. (Note that every dielectric recombination produces two photons capable of ionizing hydrogen: one due to the radiative stabilization of the process, and one as the captured electron cascades to the 1s level). If the electron density is large (ne > ∼ 10 3 cm −3 ; Clegg 1987), then more helium recombinations will result in two-photon transitions from 2 1 S-1 1 S, reducing the number of photons produced that are capable of ionizing hydrogen (Osterbrock 1989). However, the effect on Rpi is relatively small, and in any case, we do not expect to encounter large abundances of He + in dense gas in the particular scenario that we are investigating. Therefore, adopting this simplified treatment at all n should be sufficient for our purposes. Table A3) For a long time, considerable disagreement has existed on the subject of the H + 3 dissociative recombination rate. Measurements of the rate in merged beam experiments (e.g., Sundström et al. 1994) typically give values of the order of 10 −7 cm 3 s −1 for the rate coefficient at temperatures near room temperature, while measurements made in flowing afterglow experiments (e.g., Smith &Španel 1993a,b) often give values of the order of 10 −8 cm 3 s −1 , an order of magnitude smaller. At the same time, most theoretical calculations have indicated a smaller rate still, of the order of 10 −11 cm 3 s −1 (see Orel, Schneider, & Suzor-Weiner 2000, and references therein), which is in complete disagreement with the experimental measurements.
However, it has recently become clear that three-body effects play a highly important role in the recombination of H + 3 in flowing afterglow experiments (see e.g., Glosík et al. 2005). When proper allowance is made for these effects, the inferred two-body recombination rate is in good agreement with the results of the merged beam experiments (see the discussion in Glosík et al. 2007). Moreover, recent theoretical calculations of the rate coefficient by Kokoouline & Greene (2003) that account for the effects of Jahn-Teller coupling between the electronic and vibrational degrees of freedom produce a result that is in good agreement with the experimental measurements, although disagreement at the level of a factor of two or so is still present at some energies. In our calculations, we therefore take our value for the total H + 3 dissociative recombination rate coefficient from the recent ion storage ring measurements of McCall et al. (2004). To convert this total rate coefficient -the sum of the rate coefficients for reactions DR4 and DR5 -into a rate coefficient for each individual reaction, we adopt a branching ratio of 0.25 for reaction DR4 and 0.75 for reaction DR5, based on the measurements of Datz et al. (1995). Strictly speaking, these values are only appropriate for temperatures T < 3000 K, but in practice the behaviour of the gas is not particularly sensitive to the values chosen.
3.1.3 H + 3 formation by radiative association (reaction RA18, Table A5) The rate coefficient we quote for reaction RA18, the radiative association of H2 and H + to form H + 3 , was taken from the study of Gerlich & Horning (1992), and was the rate coefficient quoted by Galli & Palla (1998) for this reaction. However,  prefer a much smaller rate coefficient of 10 −20 cm 3 s −1 for this reaction. In §5.6, we examine the effect of adopting this smaller rate coefficient.
3.1.4 H2 formation by associative detachment of H − (reaction AD1, Table A6) The rate of this reaction is quite uncertain, and we have shown in previous work that this uncertainty can lead in some cases to a substantial uncertainty in the amount of H2 formed in the gas (Glover, Savin, & Jappsen 2006). However, we do not expect this uncertainty to significantly affect the results in this paper. At the high densities at which H + 3 cooling is potentially important, the dominant H2 formation pathway is three-body formation. This can produce a much larger molecular fraction than is possible via two-body reactions, and so uncertainty in the amount of H2 produced via the H − ion has a negligible impact on the evolution of the gas at high densities. In our simulations, we adopt a default value for the rate coefficient for reaction AD1 of kAD1 = 1.5 × 10 −9 T −0.1 3 cm 3 s −1 , where T3 = T /300 K, taken from Launay et al. (1991). In §5.6 we demonstrate that our results are insensitive to this choice.
We note also that the other reaction discussed in Glover, Savin, & Jappsen (2006), the mutual neutralization of H − by H + (reaction MN1, Table A8), is no longer a source of significant uncertainty in chemical models of primordial gas. Recent measurements of the cross-section for this reaction at astrophysically-relevant energies made by X. Urbain (private communication) yield values in good agreement with those obtained by Fussen & Kubach (1986) and used as a basis for the rate coefficient of Croft et al. (1999). These measurements strongly suggest that the earlier measurements of the mutual neutralization cross-section made by Moseley et al. (1970) were somehow in error, and that rate coefficients based on them (see e.g. Galli & Palla 1998) are incorrect. The error in the rate coefficient for this reaction has thus been reduced from the order of magnitude discussed in Glover, Savin, & Jappsen (2006) to an uncertainty of about 50% (X. Urbain, private communication).
3.1.5 Collisional dissociation of H2 (reactions CD9-CD12, Table A7) The rate coefficients for the collisional dissociation of H2 by H (CD9), H2 (CD10), He (CD11) and e − (CD12) are represented by functions of the form log ki = n/ncr 1 + n/ncr log ki,LTE + 1 1 + n/ncr log ki,v=0, where n is the number density of hydrogen nuclei, ki is the collisional dissociation rate for collisions with species i, and kv=0,i and kLTE,i are the rate coefficients for this reaction in the limits in which all of the H2 molecules are in the vibrational ground-state (appropriate in low-density gas), or have their LTE level populations (appropriate for high density gas). The critical density, ncr, for H2 collisional dissociation in a gas containing a mix of H, H2, He and electrons is not well determined. For simplicity, we therefore assume that it is given by a weighted harmonic mean of the (better known) critical densities corresponding to reactions CD9, CD10, and CD11 considered individually, i.e., where xH, xH 2 , and xHe are the fractional abundances of H, H2 and He relative to the total number of hydrogen nuclei, we use the approximation that xH + 2xH 2 = 1, and where ncr,H 2 = dex 4.845 − 1.3 log T4 + 1.62 (log T4) 2 cm −3 , ncr,He = dex 5.0792 with T4 = T /10000 K. The expression for ncr,H is from Lepp & Shull (1983), but has been decreased by an order of magnitude, as recommended by Martin, Schwarz, & Mandy (1996). The expression for ncr,H 2 comes from Shapiro & Kang (1987), and the expression for ncr,He comes from Dove et al. (1987). Note that this expression for the critical density assumes that in high density gas, ne ≪ nH, so that electron excitation of H2 does not significantly affect the value of ncr. Other forms of averaging to obtain ncr are possible, of course, but we would not expect our results to be sensitive to our particular choice here, as any differences will only be seen for densities n ∼ ncr, and in our simulations, gas at these densities is always far too cold for collisional dissociation of H2 to be important.
To ensure that our adopted collisional dissociation rate coefficients and three-body H2 formation rate coefficients are consistent, we used the fact that in LTE, the equilibrium constant K obeys and varies with temperature as (Flower & Harris 2007) to determine values for the rate coefficients of reactions CD9, CD10, and CD11 in the LTE limit. However, we also ran some test simulations where we used rate coefficients from Lepp & Shull (1983) and Shapiro & Kang (1987) for reactions CD9 and CD10, regardless of the value of kTB1 or kTB2. These simulations produced almost identical results, demonstrating that our results here are insensitive to our treatment of H2 collisional dissociation.
3.1.6 Collisional dissociation of HD and D2 (reactions CD13-CD20, Table A7) For collisions with electrons, accurate rate coefficients are available from Trevisan & Tennyson (2002a) and Trevisan & Tennyson (2002b). For collisions with H, H2, or He, however, we are unaware of a treatment in the literature. We have therefore assumed that the rate coefficients of these reactions in the v = 0 and LTE limits are the same as for the corresponding H reactions (nos. CD9-CD11). For D2, we have also adopted the same value for the critical density, while for HD, we have increased it by a factor of 100 to account for the larger radiative transition probabilities. Note that although these rate coefficients are highly approximate, this probably does not introduce much uncertainty into the chemical model, as reactions IX18 and IX20 (Table A10) become effective at much lower temperatures and therefore dominate the destruction of HD and D2 in warm gas.
3.1.7 Three-body H2 formation (reactions TB1-TB3, Table A9) Although unimportant at low densities, three-body reactions are the dominant source of H2 in high density primordial gas, and so these reactions represent an important part of our chemical network. Unfortunately, rate coefficients for these reactions are, in general, not known to a high degree of accuracy. For three-body collisions in which atomic hydrogen is the third body (reaction TB1), the situation is particularly bad. One commonly adopted rate coefficient is that of Palla, Salpeter, & Stahler (1983), who quote a rate coefficient kTB1,PSS = 5.5 × 10 −29 T −1 cm 6 s −1 for this reaction, based on experimental work by Jacobs, Giedt, & Cohen (1967). Also in common usage is the rate coefficient adopted by Abel, Bryan, & Norman (2002), which is kTB1,ABN = 1.14 × 10 −31 T −0.38 cm 6 s −1 T 300 K = 3.90 × 10 −30 T −1.00 cm 6 s −1 T > 300 K.
The low temperature portion of this rate coefficient is based on Orel (1987), while the high temperature portion is an extrapolation by Abel, Bryan, & Norman (2002). This reaction is also discussed by Cohen & Westberg (1983) in their large compilation and review of chemical kinetic data. They summarize a large number of different experimental measurements and argue that the precision of the data does not justify anything more elaborate than a constant rate coefficient kTB1,CW = 8.8 × 10 −33 cm 6 s −1 .
Another possibility is found in Schwenke (1990), who gives calculated values at T = 3000K and T = 5000K of 1.4×10 −32 cm 6 s −1 and 8.2 × 10 −33 cm 6 s −1 , respectively. These values are roughly an order of magnitude larger than those given by the Abel, Bryan, & Norman (2002) rate coefficient, and about 30% lower than the values given by the Palla, Salpeter, & Stahler (1983) rate coefficient. More recently, Flower & Harris (2007) have argued in favor of a rate coefficient which they derived from the rate coefficient of the inverse process (H2 collisional dissociation by atomic hydrogen, reaction CD9) by using the principle of detailed balance. This rate coefficient is approximately six times larger than the Palla, Salpeter, & Stahler (1983) rate coefficient at T = 1000K, or approximately ninety times larger than the Abel, Bryan, & Norman (2002) rate coefficient. Unfortunately, the accuracy of a rate coefficient derived using detailed balance depends upon the accuracy with which the rate coefficient of the inverse process is known. In this case, that accuracy is poor, as the H2 collisional dissociation rate coefficient is not well constrained by experiment at low temperatures (i.e., T < 2000 K) owing to its small size at these temperatures. Flower & Harris (2007) used the Jacobs, Giedt, & Cohen (1967) fit to the collisional dissociation rate coefficient, but if we instead use the calculated rate coefficient from Martin, Schwarz, & Mandy (1996), then a much smaller three-body H2 formation rate coefficient is obtained, which can be fit to within ∼ 20% by (Glover 2008) kTB1,GL = 7.7 × 10 −31 T −0.464 cm 6 s −1 .
There is thus an uncertainty of almost two orders of magnitude in the rate coefficient for reaction TB1. In our simulations, we adopt the Abel, Bryan, & Norman (2002) rate coefficient as our default value, but in §5.6 we examine the effect of using a different rate coefficient. The rate coefficient for three-body H2 formation in collisions where H2 is the third body (reaction TB2) is known with greater precision, but nevertheless substantial uncertainty remains. Palla, Salpeter, & Stahler (1983) quote a rate coefficient kTB2,PSS = 6.9 × 10 −30 T −1.0 cm 6 s −1 for this reaction, again taken from Jacobs, Giedt, & Cohen (1967), while Cohen & Westberg (1983) recommend instead kTB2,CW = 2.8 × 10 −31 T −0.6 cm 6 s −1 .
Finally, calculations by Schwenke (1988) using orbital resonance theory find a rate coefficient that is about a factor of two lower than the Cohen & Westberg (1983) values, but Schwenke (1990) shows that one of the assumptions underlying his own orbital resonance calculations is invalid, and provides revised values, obtained from a master equation approach, that agree well with the Cohen & Westberg (1983) recommendation. These rate coefficients agree to within a factor of a few at T = 3000 K, consistent with the scatter in experimental determinations of the rate coefficient at this temperature (Cohen & Westberg 1983), but differ more by more than an order of magnitude at low temperatures. In our simulations, we use the Palla, Salpeter, & Stahler (1983) rate coefficient as our default value, but we examine in §5.6 the effect of altering kTB2.
We also included three-body formation of H2 via collisions with He (reaction TB3), using a rate coefficient from Walkauskas & Kaufman (1975). This reaction has not been included in previous treatments of the evolution of dense primordial gas and so we wished to assess its effects. We found that reaction TB3 could be responsible for anywhere between 0.1% and 10% of the total three-body H2 formation rate, depending on the temperature, the H2 abundance and the choice of rate coefficients for reactions TB1 and TB2. Moreover, this estimate does not take into account the uncertainty in the rate of reaction TB3, which we have been unable to quantify, but which should probably be assumed to be comparable to the uncertainty in the other three-body rates. Thus, although it probably never dominates, reaction TB3 should be included if accurate modeling of H2 formation in dense gas is desired.
3.1.8 Deuterated three-body reactions (TB4-TB9, TB11-TB13, TB17-TB31, & TB34-TB35, Table A9) In view of the large uncertainties present in the rate coefficients of many of the three-body reaction rates (particularly for reactions TB1 and TB2, as discussed above), we consider that the most prudent course of action when estimating rates for the deuterated forms of these reactions is simply to adopt the same values as for the non-deuterated reactions. Any uncertainty introduced by this assumption is likely dwarfed by the uncertainties arising from our poor knowledge of the non-deuterated reaction rates. We note that Flower & Harris (2007) follow a similar course of action in their study of three-body H2 and HD formation in primordial gas.
3.1.9 Destruction of D2 by collision with H (reaction IX20, Table A10) Our fit to the data collated by Mielke et al. (2003) for this reaction is accurate to within a few percent over the temperature range of the tabulated data, 200 T 2200 K. Outside of this range, our fit may be significantly inaccurate (although at low temperatures, the reaction rate is small enough that any inaccuracy is unlikely to be important).
3.1.10 Photodissociation of H2 and HD (reactions BP7 & BP8, Table A12) Table A12 lists the rates of these reactions in optically thin gas, given our assumed incident UV spectrum (see §3.2). In optically thick gas, however, self-shielding of H2 by H2 and HD by HD can significantly reduce both of these rates, by factors f sh,H 2 and f sh,HD respectively. In static, isothermal gas, these self-shielding factors can be calculated approximately using the prescription of Draine & Bertoldi (1996) together with an appropriate set of molecular data, provided that one knows the H2 and HD column densities. In gas which is in motion, with internal velocities comparable to or larger than the thermal velocity of the gas, the Draine & Bertoldi (1996) treatment breaks down, and one must use approaches that are either less accurate or more computationally expensive, as discussed in Glover & Jappsen (2007). However, in the one-zone calculations presented here, we know neither the H2 and HD column densities, nor the velocity structure of the gas, and so including even a highly approximate treatment of the effects of self-shielding is problematic. In our runs with a non-zero UV background, we therefore consider two limiting cases: one in which self-shielding is highly efficient and f sh,H 2 = f sh,HD = 0 throughout the run, and one in which it is ineffective, and we remain in the optically thin limit throughout (i.e., f sh,H 2 = f sh,HD = 1). These two limiting cases bracket the true behaviour of the gas.
3.1.11 Formation and destruction of LiH + 2 (reactions RA20, DR19, DR20, DR21 & CD26) To date, the LiH + 2 ion has attracted little attention in the astrochemical literature. Kirby & Dalgarno (1978) considered the reaction chain LiH + 2 + e − → LiH + H (reactions RA20 and DR20 in our chemical model) as a possible source of LiH in the interstellar medium, but showed that even if the rate coefficient for the radiative association were assumed to be very large (kRA20 ∼ 10 −16 cm 3 s −1 ), the resulting LiH abundance would be far too small to be observable. More recently,  considered the LiH + 2 ion in their comprehensive study of the lithium chemistry of the primordial intergalactic medium, but again reached the conclusion that its abundance would be very small, and so chose not to include it in their chemical model. However, to the best of our knowledge, there has been no previous investigation of the role that this ion may play in regulating the fractional ionization of very dense primordial gas.
Previous work focussing on modelling the fractional ionization at high densities (Maki & Susa 2004, 2007, in the context of the study of ambipolar diffusion in dense population III prestellar cores, has shown that once the free electron fraction falls below x ∼ 10 −10 , ionized lithium takes over from ionized hydrogen as the primary positive ion in the gas. It is therefore important to ensure that all of the major loss routes for Li + are represented in the chemical model. In addition to photorecombination (reaction PR4), Li + can be removed from the gas by a number of reactions with atomic or molecular hydrogen: Li Most of these processes are highly endothermic, and so are not competitive with photorecombination even when x is small. However, the two radiative association reactions are exothermic and deserve closer scrutiny. Radiative association with atomic hydrogen (reaction RA10) has been included in a number of previous models of primordial gas chemistry (e.g. Galli & Palla 1998), and accurate quantal calculations of the rate coefficient for this reaction are available Gianturco & Gori Giorgi 1996). However, at the densities of interest in the present case, efficient three-body formation of H2 ensures that the hydrogen is primarily molecular rather than atomic, and so renders radiative association with H2 (reaction RA20) the more important reaction. Unfortunately, we have been unable to locate any calculation of the rate coefficient of this reaction. Kirby & Dalgarno (1978) quote an upper limit of kRA20 = 10 −16 cm 3 s −1 , while  quote an upper limit of kRA20 = 10 −17 cm 3 s −1 , but the true rate coefficient could be orders of magnitude smaller. In our reference model, we make the conservative assumption that the rate coefficient is of the same order of magnitude as that for reaction RA10, and hence adopt a value kRA20 = 10 −22 cm 3 s −1 . We investigate the effects of adopting a larger value in section §5.6.4. The inclusion of reaction RA20 in our chemical network necessitates the inclusion of additional reactions: the dominant destruction mechanisms for LiH + 2 . Unfortunately, there has been very little theoretical or experimental study of any reactions involving LiH + 2 and so it is not even clear which processes dominate. The best studied destruction process is dissociative recombination (reactions DR19, DR20 and DR21): Thomas et al. (2006) have studied this process experimentally using the CRYRING heavy ion storage ring, and have reported preliminary results regarding the branching ratio of the reaction, but have not yet reported any value for the total rate. C. Greene and collaborators are currently involved in a theoretical calculation of the total rate coefficient, but again have yet to publish any results. However, their preliminary findings suggest a total rate coefficient that is about 2.5-3 times larger than that for the dissociative recombination of H + 3 . (C. Greene, private communication). We therefore adopt a total rate coefficient 2 × 10 −7 T 300 −1/2 cm 3 s −1 for the dissociative recombination of LiH + 2 , and use the values quoted by Thomas et al. (2006) for the branching ratios.
None of these reactions appears to have previously been studied in the astrochemical literature, and so their rates are unknown. For simplicity, therefore, we include only a single representative example from this set of reactions, namely the collisional dissociation of LiH + 2 by H2 (reaction CD26). As we are primarily interested in the role of LiH + 2 within the highly molecular dense core, it is likely that this reaction will dominate, unless its rate coefficient is unusually small. In our reference model, we adopt a rate coefficient of kCD26 = 1.0 × 10 −9 exp − 3250 T cm 3 s −1 for this reaction; however, in §5.6.4 we examine the effect of adopting a smaller value.

Photochemistry
To compute rates for the photochemical reactions listed in Table A12, we assume that the gas is illuminated by an external background radiation field with the spectral shape of a 10 5 K black-body at energies hν < 13.6 eV, and which is zero at higher energies. This choice of spectrum is motivated by the fact that the brightest population III stars are expected to have high effective temperatures, T eff ≃ 10 5 K (Cojazzi et al. 2000), while the cutoff at the Lyman limit is intended to account for the effects of absorption by neutral hydrogen in the intergalactic medium. We quantify the strength of this background radiation field in terms of the flux at the Lyman limit, J(να) = 10 −21 J21 erg s −1 cm −2 Hz −1 sr −1 . The rates listed in Table A12 are computed assuming that J21 = 1.0, but scale linearly with J21 and so can easily be rescaled for other values of the background radiation field strength.

Cosmic rays
If cosmic rays are present, then they will directly ionize some species and indirectly photoionize and photodissociate others. Direct ionization is simple to treat, and the appropriate rates are listed in Table A13, normalized by the cosmic ray ionization rate of atomic hydrogen, ζH, which we treat as a free parameter. However, the indirect effects of the cosmic rays are harder to model accurately.
The basic physics is straightforward, and was first discussed by Prasad & Tarafdar (1983). They noted that the secondary electrons produced by cosmic ray ionizations are energetic enough to excite the electronic states of H2, and that the subsequent radiative decay of these excited states would produce ultraviolet photons. In the Galactic context, the mean free path of these photons is small, and so the cosmic-ray induced photochemistry can be modelled as a purely local process (see e.g., Gredel et al. 1989).
In the population III star formation context in which we are interested, however, there are two main factors that complicate matters. First, the H2 fraction in the gas is small at densities n ≪ 10 10 cm −3 , and so most of the secondary electrons produced by the cosmic rays lose energy by exciting and ionizing atomic hydrogen, rather than molecular hydrogen. Second, an accurate treatment of the propagation of the photons produced by excited H and H2 is far more involved than in the Galactic case. The continuum opacity of metal-free gas is very small at most densities of interest (Lenzuni, Chernoff, & Salpeter 1991), owing to the absence of dust absorption, and so the majority of the photons produced by cosmic-ray induced excitation of H2 have large mean free paths. On the other hand, Lyman-α photons produced by the excitation of atomic hydrogen have small mean free paths, but scatter many, many times before escaping the gas (see e.g., Dijkstra, Haiman, & Spaans 2006). In neither case is it a good approximation to assume that all of the photons are absorbed locally in the gas, and so the simple local treatment developed for Galactic dark clouds no longer applies.
If we consider only the effects of emission from H2, then to compute RX, the photoionization (or photodissociation) rate per unit volume of species X, one must use an equation of the form where the volume V that we integrate over corresponds to the entirety of the protogalactic core, and where τ (ν, x ′ , x) is the optical depth of the gas between point x and point x ′ at a frequency ν. The photon emissivity ǫ(ν, x ′ ) is given in this case by where PH 2 (ν) dν is the probability that the cosmic ray ionization of H2 leads to the production of a photon with a frequency in the range ν → ν + dν, and ζH 2 is the cosmic ray ionization rate of H2. For gas at the center of a spherically symmetric protogalactic core, we can simplify Equation 42 to where τ (ν, r) is the optical depth between the center of the halo and gas at a radius r, ǫ(ν, r) is the emissivity at r, and R is the core radius. However, even after making this simplification, calculation of RX still requires more information than we have available in our one-zone calculation, namely the radial profiles of density, temperature and chemical abundances, which determine both τ (ν, r) and ǫ(ν, r). In their absence, we are forced to approximate. If we assume that the protogalactic core has a density structure with a 'core plus halo' form, i.e., n(r) = nc r < rc nc(rc/r) α , r > rc (45) then provided that α > 1, the integral in Equation 42 will be dominated by the contribution from the core of the density profile.
If we further assume that the core is chemically homogeneous, and that τ (ν, r) ∼ 0, then we can approximate Equation 42 as where nc,H 2 is the number density of H2 within the core. Note that even if the point we are considering is not directly at the center of the core, Equation 46 remains a reasonable approximation to RX, provided that we are considering a point within rc, and that rc ≪ R. Provided that our approximations hold, Equation 46 allows us to reduce what is formally a non-local problem into one that can be treated as if it were local. To properly include the effects of hydrogen excitation, one would have to solve for the radiative transfer of the Lyman alpha photons within the collapsing protostellar core. However, as the outcome would be highly sensitive to the assumed density and velocity profiles of the gas, which are not available from our one-zone calculation, the wisdom of performing such a detailed calculation for each of our simulations that include cosmic rays is questionable; we run the risk of getting an answer that is completely determined by our assumptions, and that therefore is not robust. Instead, we have chosen a more conservative course, and have attempted to put limits on the effects of the Lyman alpha photons by considering two limiting cases: one in which they do not propagate significantly into the core of the protogalaxy, and do not contribute to RX (which is thus given in this case by Equation 46 above), and another in which the optical depth of the gas to the Lyman alpha photons is negligible, and RX is given by a generalization of Equation 46: where PH(ν) dν is the probability that the cosmic ray ionization of H leads to the production of a photon with a frequency in the range ν → ν + dν, ζH is the cosmic ray ionization rate of H, and nc,H is the number density of atomic hydrogen in the core To evaluate RX, it remains necessary to specify rc. In our calculations, we assume, following Omukai (2000) and Omukai et al. (2005), that rc is given by the current Jeans length.
In Table A14, we list estimated values for σ X,eff,H 2 = ∞ 0 σX(ν)PH 2 (ν) dν and σ X,eff,H = ∞ 0 σX(ν)PH(ν) dν for both photoionization and photodissociation for a number of different chemical species. To compute these values, we assumed that all of the photons produced by excited hydrogen are emitted in the Lyman-α line, implying that PH = δ(να − ν), where να is the frequency of Lyman-α and δ is the Dirac delta function. For PH 2 , we used estimated values based on the emission spectra given in Sternberg, Dalgarno, & Lepp (1987); note that these are likely accurate only to within a factor of a few. Given these values, the rates for the cosmic-ray induced photoionization and photodissociation of these species in our model cores can be calculated using Equation 46 or 47, as appropriate.

Neglected processes
Although the chemical network presented in this paper is, to the best of our knowledge, the most comprehensive network used to date to simulate primordial gas chemistry, there remain a large number of possible reactions that we have not included. Below, we discuss which types of processes have been omitted, and why.
(i) We do not include the formation of H2 or D2 by the radiative association of ground state atomic hydrogen or deuterium, respectively, on the grounds that the rate coefficients for these processes are negligible.
(ii) We have not included reactions that involve electronically excited atomic hydrogen (as considered in Latter & Black 1991or Rawlings, Drew, & Barlow 1993, for instance). We justify this omission by noting that the population of the n = 2 electronic level of atomic hydrogen will be very small at the densities and temperatures considered in this work, on account of the large Einstein coefficient associated with the Lyman-α transition and the large energy separation of the n = 1 and n = 2 levels. For similar reasons, we have not included any reactions that require electronically excited deuterium, helium or lithium.
(iii) We have restricted the range of chemical species considered to those with three or fewer atoms. In principle, the formation of larger species is possible -for instance, H + 5 can be formed from H + 3 by the radiative association reaction (Paul et al. 1995) However, the chemical abundances of the three-atom species in our model are very small, and we expect the abundance of even larger species to be much smaller still. It therefore seems unlikely that they will play any significant role in the cooling or chemistry of the gas.
(iv) We have omitted any photochemical reactions that require photons with energies greater than 13.6 eV, under the assumption that any such photons emitted by external sources of radiation will be absorbed in the intergalactic medium, or in the interstellar medium of the protogalaxy, before reaching the particular collapsing core under study. Moreover, since we consider only the initial collapse of the core, internal sources of radiation, such as a central protostar, fall outside of the scope of this work (v) We do not include processes that have negligible reaction rates at all temperatures treated in this work. This includes the production of doubly-ionized helium, He 2+ , or doubly or triply ionized lithium, Li 2+ and Li 3+ by collisional ionization (see e.g., Lepp, Stancil, & Dalgarno 2002), which are therefore omitted from the chemical model.
(vi) We have ignored the effects of stimulated radiative association and stimulated radiative attachment, i.e., reactions of the form where γ b represents a background photon. The influence of stimulated radiative association or attachment on the production of various molecules in primordial gas (LiH, HD, H − , Li − and HeH + ) has been investigated (Stancil & Dalgarno 1997a,b, 1998Zygelman, Stancil, & Dalgarno 1998), generally for the case of a black-body radiation field. However, significant effects are found only for radiation temperatures T rad > 500 K, much larger than the CMB temperature at the redshifts of interest in this study. The background radiation fields considered in §5.4 have the same shape below hν = 13.6 eV as a 10 5 K black-body, but have intensities that are orders of magnitude weaker than a true black-body radiation field with this temperature, and so are also unimportant in this context. Therefore, it is clear that the influence of these simulated processes will be negligible.
(vii) We have omitted collisional processes such as or that have energy thresholds corresponding to temperatures significantly higher than those considered in this work, as we do not expect these processes to play an important role in low temperature gas.
(viii) We do not include dissociative charge transfer reactions involving H + 2 or its isotopologues, e.g.
because at the temperature of interest in this study, the cross-sections for these processes are far smaller than those for the equivalent non-dissociative charge transfer reactions (Krstić 2002;) (ix) We have not included charge transfer from He + to Li or its inverse The first of these reactions is unimportant in comparison to charge transfer from H + owing to the low He + abundances we find in our simulations. The second reaction is negligible at T < 10000 K due to its large endothermicity. (x) We have omitted all collisional dissociation reactions caused by minor molecules or ions, e.g. HD, Li, LiH, etc. Collisional dissociation reactions involving HD, such as H2 + HD → H + H + HD, will be unimportant compared to the analogous reactions involving H2, while reactions of the form will be unimportant due to the very small abundance of lithium relative to hydrogen.
(xi) We have ignored a large number of possible three-body processes: specifically, every process that involves any species other than H, H2 or He as the third body. At the densities at which three-body reactions become significant, the abundances of these three species are orders of magnitude higher than the abundances of any other species, and so it is easy to justify the omission of these minor contributions.
(xii) We do not include transfer reactions involving two ions of the same charge, e.g.
as the mutual Coulomb repulsion of the ions renders these reactions ineffective at the temperatures considered in this work.
(xiii) We do not include the double ionization of H2 or He by cosmic rays, i.e.
as the fraction of cosmic ray ionization events leading to these outcomes is expected to be very small (Glassgold & Langer 1973).
(xiv) We assume that in mutual neutralization reactions involving H + 3 or one of its isotopologues, complete breakup of the ion is unlikely to occur; this is in line with e.g. the detailed chemical network of Le Teuff et al. (2000), which includes the processes and but not (xv) A number of possible reactions involving LiH + 2 have been omitted, as have any reactions involving LiHD + or LiD + 2 . As we have already discussed in §3.1.11, we know very little regarding the values of the rate coefficients for many of the most important formation and destruction processes for LiH + 2 , rendering its abundance highly uncertain. In view of this large uncertainty, there is little to be gained by adding in additional, equally uncertain but less important processes involving LiH + 2 , or by considering the chemistry of the deuterated forms of the molecular ion.

THERMAL PROCESSES
Our model of the thermal behaviour of the gas includes the effects of heating and cooling from a large number of different radiative and chemical processes. A full list of the processes included is given in Table 2, while a more detailed discussion is given below.

H2 cooling
In our treatment of H2 rovibrational cooling at low densities, we include the effects of collisions between H2 and H, H2, He, H + and e − , using fitting formulae from Glover & Abel (2008). At high densities, we use the standard LTE cooling function (see e.g., Hollenbach & Mckee 1979). We assumed the usual value of 3:1 for the H2 ortho-para ratio, which Glover & Abel (2008) have demonstrated is a good approximation for all temperatures in the range in which H2 cooling is important. Although the revised treatment of H2 cooling presented by Glover & Abel (2008) can make a significant difference to the thermal evolution of the gas in some circumstances -notably, in gas with a substantial fractional ionization -we do not expect it to have a significant impact on our current results, as at the densities of greatest interest in this paper, H2 is well within the LTE cooling regime.
An important source of inaccuracy at high densities is the treatment of the opacity of the H2 emission lines. In our models, we follow Ripamonti & Abel (2004) and model optically thick H2 cooling with the expression where Λ H 2 ,thick and Λ H 2 ,thin are the optically thick and optically thin cooling rates, respectively, n is the number density of hydrogen nuclei, and where n0 = 8 × 10 9 cm −3 and β = 0.45. This approximation works well for modelling gas at the center of a collapsing core, but is less accurate when used to treat H2 cooling in the surrounding envelope (N. Yoshida, private communication). Yoshida et al. (2006) present a more accurate approach based on the computation of escape probabilities 4.0429912 × 10 3 −1.1130317 × 10 2 0.0 a 7 −8.2863818 × 10 2 1.4187579 × 10 1 0.0 a 8 6.5975582 × 10 1 −6.9969136 × 10 −1 0.0 for each individual H2 line using the Sobolev approximation. However, this treatment requires dynamical information, in the form of the local velocity gradient, that is not available in any meaningful form in our simple models, and so for our current study we must be content with the Ripamonti & Abel (2004) approach. At very high densities (n > 10 14 cm −3 ), cooling from H2 becomes dominated by collision-induced emission (CIE). When an H2 molecule collides with a hydrogen or helium atom, or another H2 molecule, the particles involved briefly act as a 'supermolecule' with a non-zero electric dipole, which has a high probability of emitting a photon. Because the collision time is very short, the resulting collision-induced emission lines are very broad, and typically merge into a continuum. In our model, we model the CIE cooling rate with a power-law approximation taken from Ripamonti & Abel (2004), valid for gas in which xH 2 > 0.5: Although cooling from collision-induced emission contributes only 12% of the total cooling at the density reached at the end of our simulation, n f = 3 × 10 13 cm −3 , we have verified in test runs that at higher densities it very rapidly becomes the dominant form of cooling, justifying our decision to end our simulations at this point.
4.2 H + 3 cooling In the LTE limit, it is straightforward to calculate the H + 3 cooling rate using the data presented by Neale, Miller, & Tennyson (1996). In our simulations, we use tabulated values computed directly from the Neale, Miller, & Tennyson line list. However, for the convenience of readers, we also provide an analytical fit of the form where Λ H + 3 ,LTE is the H + 3 cooling rate per unit volume in the LTE limit, and where the values of the ai coefficients are listed in Table 3. This fit is accurate to within 25% for temperatures in the range 20 < T < 400 K and to within a few percent for temperatures in the range 400 < T < 10000 K. At high temperatures, the Neale, Miller, & Tennyson (1996) line list is known to be incomplete, and so for T > 3000 K, we may systematically underestimate the cooling due to H + 3 . However, as we shall see in §5, we never find gas at these temperatures in the regime where H + 3 is potentially important, and so this incompleteness will not affect our results.
At densities where H + 3 is not in LTE, the calculation of the H + 3 cooling rate presents more of a problem. A commonly used approximation for dealing with the cooling from molecular species (see e.g., Hollenbach & Mckee 1979) is to compute the cooling at a density n using the expression where ΛLTE is the LTE cooling rate per unit volume, and where the critical density ncr is given by ncr/n = ΛLTE/Λn→0, where Λn→0 is the cooling rate per unit volume in the n → 0 limit. For n ≪ ncr and n ≫ ncr, this expression is highly accurate, while for n ∼ ncr it does a reasonable job of capturing the basic behaviour, at a far smaller computational cost than a full level population calculation would require. We adopt this approximation in our treatment of H + 3 cooling, reducing the problem of calculating Λ H + 3 to one of calculating Λ H + 3 ,n→0 . Here, however, we hit a problem. To compute Λ H + 3 ,n→0 , we must evaluate where C0j is the rate of collisional excitation from the H + 3 ground state 1 , here denoted as level 0, to an excited level j, and E0j is the energy difference between level 0 and level j. The collisional excitation rate for transitions from 0 → j is simply where q0j,H, q0j,H 2 , q0j,He and q 0j,e − are the collisional excitation rate coefficients for collisions with H, H2, He and e − , respectively; nH, nH 2 , nHe and n e − are the corresponding particle number densities; and where we have ignored the effect of collisions with protons (which are unimportant in the case of a positively charged ion such as H + 3 ) or with minor ionic or molecular species such as H − or HD. The difficulty in computing the collisional terms, and hence the low density limit of the H + 3 cooling rate, arises because most of the required collisional excitation rate coefficients are unknown. Faure & Tennyson (2003) give rate coefficients for the collisional excitation of a number of low-lying rotational states by collisions with electrons, but in the high density, low ionization conditions of interest in this study, collisions with electrons are unimportant, and analogous datasets for collisions with H, He or H2 are not available.
To deal with this problem, we have used an approach introduced by Oka & Epp (2004). They computed rate coefficients for rotational transitions in H + 3 caused by collisions with H2 by making use of the principle of detailed balance and by assuming that the collisional transitions are completely random (i.e., that they obey no selection rules). These assumptions led them to suggest rate coefficients of the form for transitions between an initial level i and final level j, where gi and gj, are the statistical weights of levels i and j, respectively, Ei and Ej are the corresponding level energies, and Kij is a normalizing factor given by where C is the total collision rate, which is independent of i and j, and where the summation does not include levels i or j. Although Oka & Epp (2004) consider only pure rotational transitions, the same scheme can be used to treat ro-vibrational transitions.
In our treatment, we assume that the Oka & Epp scheme can be used to treat collisions with atomic hydrogen and helium as well as H2, and hence are able to determine the temperature dependence of the set of collisional excitation rate coefficients for each collider (q0j,H, q0j,H 2 and q0j,He); we ignore collisions with electrons, on the grounds of the very small electron abundance that exists at the densities where H + 3 cooling is potentially important. Using these collisional excitation rate coefficients, we can then construct C0j via Equation 72, from which Λ H + 3 ,n→0 follows via Equation 71. The overall normalization of the cooling rate remains uncertain, as it depends on the H + 3 number density, and on the total collision rates with each of H, H2 and He, which we can write as CH, CH 2 and CHe. If we define the total collision rate C to be the sum of these three unknowns and write the low density H + 3 cooling rate as then it is easy to show that the combination L H + 3 ,n→0 /C is completely determined. We have computed an analytical fit to this quantity, using a fit of the form of Equation 69. The fitting coefficients are listed in Table 3. This fit is accurate to within 1% over the temperature range 20 < T < 10000 K. To convert from L H + 3 ,n→0 /C to L H + 3 ,n→0 , we must fix the size of our remaining free parameter, the total collision rate C. In most of our simulations, we assume that C is given by C = 2.2 × 10 −9 nH + 1.9 × 10 −9 nH 2 + 8.1 × 10 −10 nHe s −1 which is the sum of the Langevin rates for collisions between H + 3 and H, H2 and He, respectively. These Langevin rates were computed using polarizabilities for He and H2 taken from Huiszoon & Briels (1993); the exact value was used for H. The true value of C is unlikely to be very much larger than this, but could be considerably smaller, and so in §5.2 we examine the sensitivity of our results to our choice of value for C. Finally, we note that as both Λ H + 3 ,LTE and Λ H + 3 ,n→0 are directly proportional to n H + We assume that the H + 3 emission remains optically thin throughout our simulations. For a subsonic collapse in which the effect of large-scale velocity gradients are unimportant compared to the local Doppler broadening of the emission lines, the optical depth at line center corresponding to a given emission line can be written as (Mihalas & Mihalas 1984) τji = gj gi where gj and gi are the statistical weights of levels i and j, νij is the frequency of the transition from level j to level i, Aij is the corresponding spontaneous radiative transition rate, Ni is the column density of absorbers in level i, and ∆νD = (νij/c)(2kT /m) 1/2 is the Doppler width of the line, where m is the mass of the H + 3 ion. For simplicity, we have neglected the effects of stimulated emission. Illustrative values for νij and Aij for a strong vibrational transition are νij ≃ 8.46 × 10 13 Hz and Aij ≃ 94 s −1 (Neale, Miller, & Tennyson 1996), and so for this transition where we have assumed a gas temperature of 1000 K. Therefore, in this example, the line becomes optically thick only once Ni > ∼ 10 15 cm −2 . As many of the H + 3 emission lines are considerably weaker than this example, and as the column density of H + 3 ions in any particular level i can clearly be no larger than the total H + 3 column density, N H +

3
, it is safe to conclude from this analysis that optical depth effects are unlikely to significantly affect the H + 3 cooling rate until N H + 3 > 10 15 cm −2 .
We now investigate whether we expect our models to reach this column density in H + 3 . Our one-zone dynamical model does not contain any information about the overall structure of the collapsing core and so does not predict N H + 3 directly.
However, based on the results of more detailed numerical simulations (Abel, Bryan, & Norman 2002;Yoshida et al. 2006), we assume that the protostellar core has a density profile that is well approximated by a power law n(r) ∝ r −2.2 , and that it is collapsing subsonically. With this assumed density profile, the column density of hydrogen nuclei along a radial ray from a point r to the edge of the core is given by = 5 6 rn(r) 1 − r rcore where n(r) is the number density of hydrogen nuclei at r and rcore is the radius of the core. As we shall see in §5, the H + 3 abundance in the collapsing gas typically varies only slightly with density below some threshold density n thr and then declines sharply for n > n thr . The value of n thr depends on factors such as the cosmic ray ionization rate and the speed of the collapse, but even in the most extreme models (e.g., run CR5; see §5.3) n thr = 10 11 cm −3 , while in general it is much smaller. Therefore, almost all of the contribution to N H + 3 comes from gas at densities n < n thr , and hence at radii r > r thr , where r thr is the radius such that n(r thr ) = n thr . If we assume thar r thr ≪ rcore, or in other words that the density distribution of the collapsing core extends to densities n ≪ n thr , then Equation 81 tells us that the column density of hydrogen nuclei between r thr and the edge of the core is approximately Ntot(r thr ) ≃ 5 6 r thr n thr .
Denoting the fractional abundance of H + 3 at r thr as x H + 3 (r thr ), and assuming that it remains constant for r > r thr , we can therefore write the H + 3 column density between r thr and the edge of the core as 2005) rate remains reasonably accurate down to temperatures as low as 50 K, with errors no greater than 20%, and that even at T = 30K it remains accurate to within a factor of two. At temperatures T ≫ 100K, the Lipovka, Núñez-López, & Avila-Reese (2005) cooling rate slightly underestimates the effects of HD cooling compared to the newer calculations of Wrathmall, Gusdorf, & Flower (2007), presumably owing to the more accurate vibrational excitation rates used in the latter study, but the differences are never greater than about 50%, and in any case occur in the temperature regime in which H2 cooling dominates. The breakdown of the Lipovka, Núñez-López, & Avila-Reese (2005) fit at very high temperatures (T > 20000 K) is unimportant, as the gas in our models never reaches this temperature.
To correctly model the effects of HD cooling at low temperatures, it is necessary to take the effects of the cosmic microwave background into account. We do this approximately, by using a modified HD cooling rate, Λ ′ HD , defined as where ΛHD(T ) and ΛHD(TCMB) are the unmodified HD cooling rates at the gas temperature T and the CMB temperature TCMB, respectively. The quoted range of densities for which the Lipovka, Núñez-López, & Avila-Reese (2005) cooling function is valid is 1 < n < 10 8 cm −3 . To extend the range of the cooling function to densities n < 1 cm −3 , we assume that at these densities the HD cooling rate scales proportionately to the number density of HD times the number density of colliders (i.e. as n 2 ), and hence that ΛHD(n = n ′ ) = (n ′ ) 2 ΛHD(n = 1) for n ′ 1 cm −3 , where ΛHD(n) is the HD cooling rate per unit volume (with units erg cm −3 s −1 ) at gas number density n. To extend the cooling function to high densities, n > 10 8 cm −3 , we assume that the HD molecule is in LTE and hence that the HD cooling rate per unit volume is independent of the number density of colliders and scales linearly with n; or in other words, that ΛHD/nHD is independent of n. In view of the fact that 1 ≪ ncr,HD ≪ 10 8 cm −3 , where ncr,HD is the HD critical density, both of these assumptions appear well justified.
The Lipovka et al. (2005) cooling function only includes the effects of collisions between HD and H. However, Flower et al. (2000) have shown that the influence of the H2/H ratio on the HD cooling rate is very small, and so the Lipovka et al. (2005) cooling function should remain reasonably accurate even after the molecular fraction becomes large. Moreover, collisions between HD and other species (electrons, H + , etc.) can be neglected compared to collisions with H on account of the much larger abundance of the latter at the gas densities of interest.
Finally, we assume that the HD rovibrational lines remain optically thin throughout all of our runs. In practice, this is probably not the case: the strongest HD lines will become saturated once the HD column density exceeds NHD = 10 22 cm −2 , and an analysis similar to that performed for H + 3 in the previous section suggests that this will occur in our model cores once n > ∼ 3 × 10 13 cm −3 . However, HD is only a minor coolant at these high densities, and so we can safely neglect optical depth effects on HD cooling without significantly affecting the thermal evolution of the gas in our models.

LiH cooling
To treat cooling from LiH, we use the cooling function given in Galli & Palla (1998). This is the low density limit of the LiH cooling rate and so is strictly valid only for gas densities significantly below the LiH critical density, ncr,LiH. However, the transition probabilities for the rotational and vibrational transitions of LiH are very large, on account of the molecule's large dipole moment. This means that the LiH critical density is large, ncr,LiH ≃ 10 12 cm −3 (Lepp & Shull 1984), and so the Galli & Palla (1998) cooling function is a reasonable choice over most of the range of densities covered by our simulations. At very high densities, we would expect LiH to begin to reach LTE, and our approximation to break down; at these densities, our continued use of the Galli & Palla (1998) cooling function means that we will overestimate the effectiveness of LiH cooling. Despite this, we find LiH cooling to be ineffective at all densities (see §5.1 below), suggesting that if we were to use a more accurate treatment of LiH cooling at high densities it would not significantly alter our conclusions.
In principle, we should adjust the LiH cooling rate to account for the effects of the CMB, just as we do for the HD cooling rate. However, as LiH cooling proves to be unimportant at all densities, this correction is also unimportant and its omission does not significantly affect the thermal evolution of the gas.
We assume that the LiD cooling rate is the same as the LiH cooling rate. While this is a crude approximation, in practice the LiD abundance is so small that its contribution is always negligible and thus the choice of LiD cooling rate is unimportant.

H +
2 , HD + and D + 2 cooling To treat cooling from these molecular ions, we use the same approach as in Glover & Abel (2008). At low densities, most H + 2 cooling occurs due to collisions with free electrons and neutral hydrogen atoms; collisions with He and H2 excite H + 2 at comparable rates to collisions with H (Roberge & Dalgarno 1982), but are unimportant due to the low abundances of these species relative to atomic hydrogen. To model the cooling due to collisions with electrons, we use the vibrational excitation rate coefficients of Sarpal & Tennyson (1993), while for collisions with neutral hydrogen, we use a fit to the Suchkov & Shchekinov (1978) rate coefficient provided to us by D. Galli (private communication); note that this is a factor of ten smaller than the rate coefficient given in Galli & Palla (1998), owing to a normalization error in the latter paper.
At high densities, the vibrational levels of H + 2 will be in LTE. In this regime, the cooling rate is given approximately by This fit is from Glover & Abel (2008) and includes contributions from all vibrational states v 8 (higher vibrational states are not expected to contribute significantly at the temperatures of interest in this work). It was computed using level energies from Karr & Hilico (2006) and radiative transition rates from Posen, Dalgarno, & Peek (1983). The effects of rotational excitation were not included, but are unlikely to change this expression by a large amount, owing to the very small Einstein coefficients associated with pure rotational transitions in H + 2 . At intermediate densities, we assume that the H + 2 vibrational cooling rate is given approximately by the function where ncr/n = Λ H + 2 ,LTE /Λ H + 2 ,n→0 is the H + 2 critical density, and where Λ H + 2 ,n→0 is the H + 2 cooling rate in the low density limit. This is given by where L H + 2 ,e and L H + 2 ,H are the cooling rates per H + 2 ion per unit collider density for collisions with electrons and atomic hydrogen, respectively, taken from Sarpal & Tennyson (1993) and Suchkov & Shchekinov (1978) as noted above.
To model cooling from vibrational transitions in HD + , we assume, in the absence of better information, that the low density cooling rate is the same as that used for H + 2 . However, since HD + has much larger radiative transition rates than H + 2 , the LTE cooling rate for HD + is much larger than that for H + 2 . Accordingly, we use the following functional fit for the HD + LTE cooling rate: Λ HD + ,LTE = 1.09 × 10 −11 T 0.03 exp − 2750 T n HD + at temperatures T 1000 K and Λ HD + ,LTE = 5.07 × 10 −12 T 0.14 exp − 2750 T n HD + at T > 1000 K. These fits are from Glover & Abel (2008) and were calculated using HD + level energies from Karr & Hilico (2006) and radiative transition rates from Peek, Hashemi-Attar, & Beckel (1979). For densities between the low density and LTE limits, we again use a function of the form of Equation 87 to compute the cooling rate. Finally, to model D + 2 cooling, we simply assume that the same rates apply as for H + 2 cooling. In practice, the very small size of the typical D + 2 abundance renders this process completely unimportant.

Other radiative coolants
In addition to the coolants discussed above, we also include a treatment of cooling from a number of other minor molecular ions that are present in the gas. Specifically, we include the effects of cooling from H2D + , HD + 2 , D + 3 , HeH + , HeD + , He + 2 , LiH + , LiD + and LiH + 2 . Since appropriate collisional data are not available for most of these species, we treat their contribution to the cooling rate in an extremely simple fashion. We assume that the contribution to the cooling rate made by a species i with number density ni can be written as where nc is the number density of a collider c, Cic is the rate coefficient for inelastic collisions between i and c, and where we sum over all possible colliders. For collisions with H, H2 or He, we assume that Cic is given by the Langevin rate, while for collisions with electrons we conservatively assume that Cic = 10 −6 cm 3 s −1 , which is comparable to the total rate coefficients found for other molecules, such as H + 3 (Faure & Tennyson 2003). Collisions with all other species can be and are neglected. In constructing this approximation we have assumed that each collision with i transfers an amount of energy kT , all of which is subsequently radiated. In practice, this procedure is likely to significantly overestimate the cooling provided by i, for several reasons. For one thing, it is not clear that the mean amount of energy transferred in a collision will always be ∼ kT , since collisions that transfer ∆E ≪ kT are possible while collisions that transfer ∆E ≫ kT are highly unlikely. More importantly, this procedure neglects effects such as the collisional de-excitation of excited levels that will significantly limit cooling at high densities. However, these simplifications are unlikely to significantly affect our results, since even when we use these overestimates for the cooling rates, we find that the contribution of these minor species to the total cooling rate is negligible (see §5).
As well as these minor coolants, we also include two forms of cooling that are of great importance in hot, ionized gas. The first is cooling from electron impact excitation of atomic hydrogen (Lyman-α cooling). We treat this using a rate from Cen (1992), which is itself based on a rate in Black (1981). However, we note that for temperatures T < 8000 K, Lyman-α cooling is completely negligible, and so it does not significantly affect the outcome of the simulations presented in this paper.
The second process is cooling due to the Compton scattering of CMB photons by free electrons. This is also treated using a rate from Cen (1992), but again plays very little role in the thermal evolution of the gas, since it is important primarily at low densities (n < ∼ 1 cm −3 ), even when the gas is initially highly ionized. Finally, we note that we do not include the effects of cooling from D2. As D2 is a homonuclear molecule, it suffers from the same drawbacks that H2 does with regard to low temperature cooling. However, as the results in §5.1 demonstrate, it generally has a chemical abundance that is many orders of magnitude smaller than H2. Thus, in contrast to HD, it appears highly unlikely that D2 cooling is ever significant.

Radiative heating
We include two forms of radiative heating that can be significant if a strong ultraviolet background is present. The first is the photodissociation of H2. We calculate the H2 photodissociation rate as discussed in §3, and then, following Black & Dalgarno (1977), we assume that each photodissociation deposits 0.4 eV of heat into the gas. Note that although the photodissociation of other ionic and molecular species (e.g., H + 2 , HD) will also heat the gas, they are unimportant when compared to H2 photodissociation owing to the low abundances of the other species relative to H2.
The second form of radiative heating included in our thermal model arises due to the population of excited vibrational states of H2 produced by radiative pumping by the UV field. At high densities, this leads to heating of the gas, as most of the excited molecules undergo collisional de-excitation. We adopt a radiative pumping rate that is 8.5 times larger than the photodissociation rate (Draine & Bertoldi 1996), and assume that each excitation transfers an average of 2 (1 + ncr/n) −1 eV to the gas (Burton, Hollenbach, & Tielens 1990), where ncr is the H2 critical density, calculated as discussed in §3 above.

Cosmic ray heating
Following Goldsmith & Langer (1978), we assume that each primary ionization deposits 20 eV of energy into the gas, giving us a heating rate Γcr = 3.2 × 10 −28 ζ 10 −17 s −1 n ergs s −1 cm −3 , where ζ = i ζi and we sum over all species listed in Table A13.

Chemical heating and cooling
Any exothermic chemical reaction will potentially heat the gas, while any endothermic reaction will cool it. In practice, however, the effect of most reactions on the gas temperature is small, and only in a few cases do we need to take chemical heating or cooling into account. In highly ionized gas, cooling due to the recombination of hydrogen and helium can be a significant effect, particularly at temperatures which are too low for Lyman-α cooling to be effective. However, recombination cooling becomes ineffective once the fractional ionization of the gas falls below x ∼ 0.01, and it therefore plays no role at the high gas densities of interest in this study.
Other forms of chemical cooling included in our model occur due to the collisional dissociation of H2 (reactions CD9, CD10, CD11 and CD12), the destruction of H2 by charge transfer with H + (reaction CT2), and the collisional ionization of hydrogen and helium (reactions CI1 & CI2). Cooling from these processes may be of some importance at very early times in simulations starting at high temperatures (T > ∼ 10 4 K), but in general the gas temperature is too low for these sources of cooling to be significant.
As far as chemical heating is concerned, the most significant process is H2 formation heating. When H2 is formed by reaction AD1 or reaction CT1, it preferentially forms in an excited vibrational state, with an energy comparable to the exothermicity of the reaction (3.73 eV for reaction AD1, 1.83 eV for reaction CT1). In low density gas, this energy is simply radiated away, but for n > 10 4 cm −3 , most is instead converted into thermal energy by collisional de-excitation of the newly formed H2. H2 formation via H − or H + 2 therefore acts as a minor heat source in gas with n > 10 4 cm −3 .
Three-body formation of H2 also heats the gas, since the third body in the collision generally carries away additional energy equal to the binding energy of the new H2 molecule, 4.48 eV. In our reference simulation, this is the dominant heat source for densities 5 × 10 10 < n < 2 × 10 12 cm −3 .

RESULTS
Our simple one-zone dynamical model of gravitationally collapsing primordial gas contains a number of free parameters. To fully explore the role of H + 3 cooling and its sensitivity to these free parameters, it is necessary to perform a large number of calculations. However, discussion of the results of all of these calculations to the same level of detail would not only be extremely tedious, but would also run the risk of obscuring our main results. Therefore, we proceed by first discussing in detail in §5.1 the results of a single calculation -our reference model, hereafter denoted as computational run REF -before highlighting in the subsequent sections the differences in outcome (if any) that result from alterations in our free parameters. Full details of all of the runs discussed here can be found in Table 4.

The role of H +
3 cooling We begin our study by investigating the outcome of our reference calculation, run REF, whose parameters are indicated in Table 4. In Figure 1 we show how the fractional abundances of 28 of our 30 chemical species vary with density during the course of the collapse. For clarity, we have divided these species into four sets on the basis of the elements that they contain, and illustrate the evolution of each set separately in Figures 1a-1d. The two species that are not plotted -He + and He + 2have abundances that remain negligibly small throughout the calculation. Figure 1 demonstrates that the evolution of the H + 3 abundance passes through four distinct phases. In the first phase, at n < ∼ 10 8 cm −3 , the ratio of H + 3 to H2 remains approximately constant: x H + 3 /xH 2 ∼ 10 −9 . This is a consequence of the balance between two main processes: the formation of H + 3 by the radiative association of H2 and H + (reaction RA18) and its destruction by dissociative recombination (reactions DR4 and DR5). If we assume that these reactions dominate the formation and destruction of H + 3 , then the corresponding equilibrium abundance of H + 3 is given by which reduces to if xe = x H + , which is a very good approximation at these densities, as Figure 1a demonstrates. The second phase in the evolution of the H + 3 abundance begins at a density of around 10 8 cm −3 , when there is a sudden decrease in the H + 3 abundance. This decline is caused by the fact that at these densities, dissociative recombination is no longer the only significant destruction mechanism. Increasingly, H + 3 is also destroyed by reaction TR17: Although this reaction is endothermic, the temperature of the gas at these densities (T > ∼ 800 K; see Figure 2) is high enough to make this mechanism significant in comparison to reactions RA18, DR4 and DR5, owing to the very low fractional ionization of the gas.
The third phase occurs at n ∼ 10 10 cm −3 as the decline in the H + 3 abundance is briefly halted by an increase in the H + 3 formation rate. This is caused by the increase in the H2 abundance at these densities, which itself is driven by the onset of efficient three-body formation of H2. Finally, at a density n > ∼ 10 11 cm −3 , the H + 3 abundance decreases once more, owing to the rapid loss of the few remaining free H + ions from the gas, and the consequent disruption of the major H + 3 formation mechanisms. The reaction responsible for this loss of H + ions is the same as the reaction driving the formation of H + 3 , namely RA18. If this were the only process operating, then it would convert all of the H + in the gas into H + 3 in a time tconv, given approximately by Comparing this timescale with the free-fall timescale, t ff ≃ 1.4 × 10 15 n −1/2 s, and taking kRA18 = 10 −16 cm 3 s −1 , we find that tconv < t ff if n > 50/x 2 H 2 . For xH 2 ∼ 10 −3 , this gives a critical density nconv ∼ 5 × 10 7 cm −3 . Therefore, in the absence of any other effects, conversion of H + to H + 3 should occur rapidly once the number density exceeds nconv. In practice, however, a second effect intervenes. The steady increase in the gas temperature at these densities soon results in reaction TR17 becoming a major destruction mechanism for H + 3 , as noted above. Destruction of H + 3 by reaction TR17 produces H + 2 ions, most of which are then destroyed by reaction CT3 producing H + ions. Therefore, most of the H + ions that are removed from the gas by reaction RA18 are replaced by this chain of reactions. A net loss of H + from the gas occurs only if the H + 3 ion produced by reaction RA18 is destroyed by dissociative recombination, rather than by reaction TR17. The proportion of the H + 3 destroyed by dissociative recombination is given by where we have assumed that dissociative recombination (reactions DR4 and DR5) and reaction TR17 are the only significant processes destroying H + 3 . Consequently, the net rate at which H + ions are removed from the gas is a factor fDR slower than was assumed in our calculation of tconv above, and hence the actual timescale on which the majority of the H + ions are removed is given by Now, t loss depends on the electron density through fDR, and so as long as H + is the dominant source of free electrons, decreasing its abundance increases t loss , thereby preventing rapid removal of the H + ions from the gas. However, once the H + abundance falls below ∼ 10 −11 , it is singly ionized lithium, Li + , that becomes the dominant positive ion. At this point, further decreases in x H + have very little effect on t loss . Taking x e − = 10 −11 and T = 1000 K, and assuming that (kDR4 + kDR5)ne ≪ kTR17nH, we find that fDR ≃ 10 −4 x −1 H , and hence t loss ≃ 10 20 xH/nH 2 s. Thus, at the point at which Li + first becomes the dominant positive ion, which occurs around n ∼ 10 9 cm −3 in our reference simulation, t loss ≫ t ff . However, the wholesale conversion of H to H2 by three-body reactions that begins to set in at around this density rapidly decreases t loss , and at n ∼ 10 11 cm −3 it becomes shorter than the free-fall timescale of the gas. At this point, most of the remaining H + ions are lost from the gas, following which H + 3 formation largely ceases. Since the destruction of H + 3 by reactions DR4, DR5 and TR17 is unaffected by the fall-off in the H + abundance, the end result is a very rapid fall-off in the H + 3 abundance. To determine whether the small amount of H + 3 that forms in the gas is enough to significantly affect its thermal evolution, we have compared the H + 3 cooling rate per unit volume to the total cooling rate per unit volume. The results are plotted in Figure 3. We see that even at its moment of peak effectiveness, which occurs at n ∼ 10 8 cm −3 , H + 3 contributes no more  than about 3% of the total cooling rate. At lower densities, the H + 3 ions, which are not yet in LTE, contribute less of the cooling because they undergo fewer collisions. At higher densities, on the other hand, the effectiveness of H + 3 is reduced by the significant decrease in its chemical abundance, even though each individual H + 3 ion contributes more cooling than at n = 10 8 cm −3 .
If H + 3 is ineffective, then what about other potential coolants, such as LiH, H + 2 and its deuterated counterparts, or other ions such as HeH + or H2D + ? As far as LiH is concerned, Mizusawa, Omukai, & Nishi (2005) have already shown that far too little forms for it ever to be a significant coolant, a result which we confirm (c.f., our Fig. 1d with their Fig. 1b). The contributions made by the other possible coolants are assessed in Figures 4a and 4b, where for convenience we plot only the sum of the contributions from two sets of species. For the ions in one of these sets (H + 2 , HD + and D + 2 ; Fig. 4a), we have cooling functions that should be at least reasonably accurate; for those in the other set (H2D + , HD + 2 , D + 3 , HeH + , HeD + , He + 2 , LiH + , LiD + and LiH + 2 ; Figure 4b), we use the highly approximate treatment described in §4.6. It is clear from the figures that none of these species contribute significantly to the total cooling rate, which is unsurprising given their extremely small abundances throughout the range of densities examined here (see also Figure 1).

Sensitivity to the choice of H +
3 collision rate As we saw in the previous sub-section, one of the factors preventing H + 3 from becoming a dominant coolant in our reference calculation is the fact that its abundance begins to decrease, owing to the increasing importance of destruction by collisions with hydrogen atoms, before its cooling rate has reached its LTE limit. If the H + 3 ion were to reach LTE earlier than we have  assumed -in other words, if the low density limit of its cooling rate were to be larger -then H + 3 cooling would have more effect. We have therefore explored the effect of altering C, the total H + 3 collisional excitation rate coefficient that is the single free parameter in our treatment of H + 3 cooling. The value of C in our reference model, hereafter C ref , is given by Equation 77. Increasing it or decreasing it compared to this value has the effect of, respectively, increasing or decreasing the low-density H + 3 cooling rate; or, equivalently, decreasing or increasing the critical density at which H + 3 reaches LTE. In Figure 5, we show the effect that increasing or decreasing C by a factor of ten has on the contribution made by H + 3 to the total cooling rate. As one would expect, the result is rather dramatic. In particular, it is clear that if C = 10 C ref , then H + 3 cooling does contribute significantly to the total cooling rate around densities n ∼ 10 8 cm −3 . However, such a large value for C seems unrealistic, given that in our expression for C ref , we are already assuming that collisions occur at the Langevin rate. Collisional excitation by electrons could in ideal circumstances give one a large value for C, but it is clear from Figure 1 that the electron abundance is orders of magnitude too low in the present case for collisions with electrons to be important. Moreover, even if C were as large as 10 C ref , the extra cooling provided by the H + 3 ions would have only a small effect on the temperature evolution, as Figure 6 demonstrates.
On the other hand, if C is smaller than we have assumed, then H + 3 cooling has even less effect. Therefore, despite the uncertainties in our treatment of H + 3 cooling at low densities, our main result -that H + 3 cooling is, in general, unimportant -seems robust.

Influence of cosmic rays
In §5.1 we saw that a major reason for the dramatic decrease in the H + 3 abundance at gas densities n > 10 8 cm −3 is the loss of the remaining H + ions from the gas. In the absence of any H + ions, H + 3 can no longer be produced directly by reaction RA18 (formerly, the dominant production mechanism), while its production from reactions involving H + 2 or HeH + is also disrupted, as the main formation routes for these species also depend upon the availability of H + . Clearly, therefore, one way of significantly enhancing the abundance of H + 3 at high densities would be to provide a source of H + ions at high densities. Alternatively, the H + 3 abundance could also be enhanced if there were a suitable source of H + 2 ions in dense gas. One mechanism capable of producing both H + and H + 2 ions in dense gas is the partial ionization of the gas by an external flux of cosmic rays. In view of the many uncertainties and unknowns regarding the composition, energy spectrum and energy density of the cosmic rays produced by the earliest supernovae, summarized in Stacy & Bromm (2007) and Jasche, Ciardi, & Ensslin (2007), we use a highly simplified treatment of their effects. We assume that all of the uncertainties can be folded into a single free parameter, ζH, the cosmic ray ionization rate of atomic hydrogen, and that the cosmic ray ionization rates of other atoms and molecules have the same scaling with ζH as they are commonly assumed to have in the local ISM. Secondary effects resulting from the Prasad-Tarafdar mechanism are treated as outlined in §3.3; note that we assume in this first set of models that Lyman-α photons make a negligible contribution to the secondary photochemical rates.
In Figure 7a, we show how the temperature evolution of the gas changes as we increase ζH. We show in the figure results from five runs that included cosmic rays: CR1, CR2, CR3, CR4 and CR5, with ζH = 10 −20 , 10 −19 , 10 −18 , 10 −17 and 10 −16 s −1 , respectively. We also plot the temperature evolution in our reference run REF, for the purposes of comparison. We see that as we increase the cosmic ray ionization rate, the gas gets colder. In particular, for ζH 10 −18 s −1 , the gas is able to cool to T < 100 K, indicative of the fact that in these runs, HD becomes the dominant low-temperature coolant. This is a simple consequence of the ionization produced by the cosmic rays: the additional free electrons allow more H2 to be produced than in our reference run (see Figure 7b), and so the gas can cool to lower temperatures. Stacy & Bromm (2007) find a similar effect in their recent study of the effects of cosmic rays on primordial star formation, and also show that the effect of the cosmic rays on the temperature evolution becomes significant once ζH 10 −19 s −1 ; we find a slightly larger critical value here, possibly due to the differences in our treatment of H2 cooling.
In Figure 7c, we show how the increase in ζH affects the contribution that H + 3 cooling makes to the total cooling rate. We see that as the ionization rate increases, H + 3 cooling becomes steadily less effective at low densities, particularly for ζH 10 −18 s −1 , owing to the growing importance of HD cooling at these densities. Above n ∼ 10 8 cm −3 , however, the contribution from H + 3 cooling increases with increasing ζH. As Figure 7d illustrates, this is a consequence of a significant increase in the high density H + 3 abundance in these runs compared to run REF, which is an expected consequence of the greater availability of H + and H + 2 at high densities in the runs with non-zero ζH. Figure 7c also demonstrates that if ζH 10 −18 s −1 , then H + 3 is responsible for > 10% of the total cooling over several orders of magnitude in gas density. Moreover, if ζH 10 −17 s −1 , then there is a brief period in which it provides > 50% of the total cooling. Thus, for cosmic ray ionization rates of this order of magnitude, H + 3 cooling is clearly significant and H + How plausible is it that the cosmic ray ionization rate in primordial gas will be as large as 10 −17 s −1 ? This value is comparable with the standard estimates for the cosmic ray ionization rate in dense gas in the local ISM (see e.g., Bergin et al. 1999;van der Tak & van Dishoeck 2000), and lower than recent estimates of the rate in diffuse gas (see e.g., , and so this value is not prima facie unreasonable. However, our requirement that the cosmic rays penetrate to very high gas densities means that they must be highly energetic. If we use the same simple model for our collapsing protostellar core as in §4.2, then at n ∼ 10 10 cm −3 , the core radius is r ∼ 10 15 cm, and the column density of the core is N ∼ 10 25 cm −2 . To penetrate to this depth, the cosmic rays must have energies of at least 100 MeV (Stacy & Bromm 2007), which means that given reasonable assumptions regarding the shape of the cosmic ray energy spectrum, the main contribution to the cosmic ray ionization rate will come from cosmic rays with roughly this energy. Following Stacy & Bromm (2007), we can estimate the required energy density in 100 MeV cosmic rays as If we further assume that most high-redshift cosmic rays are produced in the supernova remnants left by pair-instability supernovae (Heger & Woosley 2002), then Stacy & Bromm (2007) show that the total cosmic ray energy density UCR is Figure 8. Temperature evolution as a function of gas density in runs CR5 (solid line) and CR6 (dashed line). Both runs share the same set of input parameters, including a cosmic ray ionization rate ζ H = 10 −16 s −1 , but in run CR6, H + 3 cooling has been artificially disabled.
related to the cosmological star-forming rate per unit comoving volume, Ψ * , by where pCR is the fraction of the supernova explosion energy, ESN, that is used to accelerate cosmic rays, and fPISN is the number of pair-instability supernovae per solar mass of stars formed. Given reasonable values for pCR, ESN and fPISN, this relationship implies that to produce a cosmic ray energy density of order 10 −13 erg cm −3 , we require a star formation rate per unit volume of order 1 M⊙ yr −1 Mpc −3 , two to three orders of magnitude larger than current estimates of the population III star formation rate (Yoshida et al. 2003;Bromm & Loeb 2006). Moreover, this estimate assumes that essentially all of the cosmic ray energy density is in 100 MeV cosmic rays; if we allow for a significant fraction of cosmic rays with smaller energies, as are required to produce the ionization in low density gas in our simplified model, then the required cosmic star formation rate increases still further. From this argument, we can conclude that any extragalactic cosmic ray background will be too small to produce the effect that we desire. How about local sources of cosmic rays? Stacy & Bromm (2007) show that much higher energy densities can be produced close to individual supernova remnants, but to get an energy density of 10 −13 erg cm −3 one would have to be within ∼ 10 pc of the remnant, near enough that the gas would have been strongly processed by the ultraviolet radiation of the supernova progenitor (Glover & Brand 2001;Whalen, Abel, & Norman 2004;Susa 2007). Consequently, this scenario for producing a high cosmic ray ionization rate also does not appear promising.
Furthermore, even if we assume that it is possible to maintain a large ζH at high densities, and that H + 3 cooling does briefly become dominant, it is possible to show that its effects on the temperature evolution of the gas remain small. In Figure 8, we compare the temperature evolution in runs CR5 and CR6. In both runs, we have set ζH = 10 −16 s −1 , but in run CR6 we have artificially disabled H + 3 cooling. We see that at densities 10 8 < ∼ n < ∼ 10 12 cm −3 the temperature in run CR5 is smaller than the temperature in run CR6, as expected. However, the difference is relatively small, and the temperature evolution is qualitatively similar in both cases.
Finally, we have examined the importance of cosmic-ray induced photoionization and photodissocation (the Prasad-Tarafdar mechanism, discussed in section 3.3) by peforming several additional simulations. In runs CR7, CR8 and CR9, we took ζH = 10 −20 , 10 −18 and 10 −16 s −1 , respectively, but neglected the effects of the Prasad-Tarafdar mechanism completely. In runs CR10, CR11 and CR12, we adopted the same cosmic ray ionization rates, but maximized the effects of the Prasad-Tarafdar mechanism by assuming that the Lyman-α photons produced by secondary excitations of atomic hydrogen could propagate freely into the core of the protogalaxy, and could contribute to the total secondary photoionization and photodissociation rates there (cf. our standard treatment, where we assume that the Lyman-α photons are unable to penetrate into the core).
In Figure 9a we compare the temperature evolution of the gas in these six runs with the evolution in runs CR1, CR3, and CR5, which have the same values of ζH, but which include the effects of cosmic-ray induced photoprocesses. In Figure 9b, we show a similar comparison of the ratio of H + 3 cooling to total cooling in these runs. Figure 9. (a) Sensitivity of the temperature evolution to our treatment of the Prasad-Tarafdar mechanism. Solid lines correspond to runs using our default treatment, dashed lines to runs that neglect its effect entirely, and dash-dotted lines to runs that maximize its effects by including the effects of Lyman-α photons as if the gas were optically thin to them. The lower, middle and upper sets of curves correspond to runs with ζ H = 10 −20 , 10 −18 and 10 −16 s −1 , respectively. Note that many of the lines in this figure are indistinguishable. (b) As (a), but showing the ratio of the H + 3 cooling rate to the total cooling rate in the same set of runs. The lower, middle and upper solid lines on the right-hand side of this plot correspond to runs CR1, CR3 and CR5, respectively, the lower, middle and upper dashed lines on the same side correspond to runs CR7, CR8 and CR9, and while the lower, middle and upper dash-dotted lines on that side correspond to runs CR10, CR11 and CR12. (Note that the results of runs CR5 and CR12 are indistinguishable in the plot, as are the results of runs CR1, CR7 and CR10).
Below n = 10 8 cm −3 , the Prasad-Tarafdar mechanism has no discernable effect on the evolution of the gas. At higher densities, however, its effect is to suppress H + 3 cooling. More specifically, cosmic-ray induced photodissociation of H + 2 (reaction CP3) reduces the amount of H + 3 formed via reaction TR3, leading to a reduction in the H + 3 abundance at these densities, and hence an overall reduction in the effectiveness of H + 3 cooling. Nevertheless, the effect of the enhanced cooling in runs including its effects is slight, as can be seen from Figure 9a. Figure 9 also allows us to assess the impact of the inaccuracy in our treatment of the Lyman-α photons photons produced by secondary excitations of atomic hydrogen. Our results indicate that when ζH is small, the Prasad-Tarafdar mechanism has little effect on the amount of H + 3 cooling that occurs, and so any inaccuracies in our treatment of it are unimportant. In high-ζH runs, the Prasad-Tarafdar mechanism is more important, but we find that we obtain very similar results with and without the inclusion of the Lyman-α photons, and so any inaccuracy in their treatment is again unimportant.

Influence of a radiation background
In most of our calculations, we have assumed that any external sources of radiation have a negligible effect on the evolution of our collapsing protogalactic gas. At the epoch corresponding to the formation of the very first stars, this assumption is well-justified: the cosmic microwave background does not significantly affect the gas chemistry at redshifts z < 100 and no other sources of radiation yet exist. Once population III star formation has begun, however, the situation changes. Radiation from massive population III stars or from their remnants can affect primordial gas through a variety of mechanisms, as discussed in detail in the recent reviews by Ciardi & Ferrara (2005) and Ciardi (2008). In this section, we explore how an external radiation field can influence the thermal evolution of primordial gas and affect the role of H + 3 cooling.

Ultraviolet radiation
One of the most important forms of radiative feedback in the high-redshift Universe is the build-up of a soft ultraviolet background at photon energies hν < 13.6 eV. Photons from this background that are absorbed in the Lyman and Werner band transitions of H2 can cause photodissociation, and since these photons can propagate to large cosmological distances through the intergalactic medium, the strength of the background and the size of the associated photodissociation rate can both become considerable. This soft UV background is therefore expected to have a significant effect on the evolution of primordial gas within protogalaxies (  and UV3 (dotted line). The strength of the ultraviolet background in these runs was J 21 = 0.0, 10 −4 , 10 −2 and 1.0, respectively, and the gas was assumed to remain optically thin throughout its evolution. (b) As (a), but for runs REF, UV4, UV5 and UV6. Runs UV4, UV5 and UV6 had the same UV background field strengths as runs UV1, UV2 and UV3, respectively, but in this case we assumed that the gas was optically thick in the Lyman-Werner lines of H 2 and HD. The results of the four runs are indistinguishable in the plot.
To investigate the impact of such an ultraviolet background on our results, we have run several models with non-zero backgrounds: runs UV1, UV2, UV3, UV4, UV5 and UV6. As previously noted in §3.2, we assume that the spectral shape of the background is that of a diluted 10 5 K black-body, with a sharp cutoff at 13.6 eV. The only free parameter is then the normalization of this spectrum. We choose to normalize it by specifying its strength at the Lyman limit. In runs UV1, UV2 and UV3, the field strength is J21 = 10 −4 , 10 −2 , and 1.0, respectively, where J21 is the flux at the Lyman limit in units of 10 −21 erg s −1 cm −2 Hz −1 sr −1 . In these three runs, we assume that self-shielding by H2 and HD is not effective, and that the gas remains optically thin to the external radiation field throughout the simulation. In runs UV4, UV5 and UV6, the field strength is the same as in runs UV1, UV2 and UV3, respectively, but we assume that H2 and HD self-shielding is so effective that the H2 and HD photodissociation rates are negligible. The true behaviour of the gas lies between these two limiting cases.
In Figure 10a, we show how the gas temperature evolves in optically thin runs UV1, UV2 and UV3, as well as in run REF for comparison. The corresponding behaviour in runs UV4, UV5 and UV6 is shown in Figure 10b. In the optically thin runs, the effect of the ultraviolet background is to increase the temperature of the gas at early times; quite dramatically so in the case of run UV3, where the minimum temperature reached by the gas is T ∼ 900 K, compared to only T ∼ 200 K in our reference calculation. This temperature increase is an obvious consequence of the photodissociation of H2 in low density gas, as can be seen clearly by comparing these results with those from the runs in which H2 photodissociation was assumed to be negligible.
As far as H + 3 cooling is concerned, Figure 11 demonstrates that it remains ineffective in both sets of runs. In the optically thin runs, H2 dissociation at early times reduces the effectiveness of H2 cooling, but also significantly reduces the H + 3 abundance. The net effect is to reduce the amount of cooling coming from H + 3 to below the level that it has in our reference run. At densities n > ∼ 10 9 cm −3 , on the other hand, the effect of the UV background is to enhance cooling from H + 3 . This occurs because the H + abundance does not decline so quickly in the runs with the UV background, owing to the higher gas temperature, and so more free protons are available for making H + 3 at very high densities. This effect boosts the contribution of H + 3 to the cooling rate in run UV3 by about a factor of four compared to our reference run. Despite this, however, the contribution from H + 3 remains unimportant. In runs UV4, UV5 and UV6, the behaviour of the contribution from H + 3 is somewhat different. The H2 abundance in these runs evolves in almost the same manner as in run REF, as does the temperature. Therefore at early times, there is no difference in the H + 3 contribution. At n > 10 7 cm −3 , however, a difference does become apparent between run REF and runs UV5 and UV6, with the H + 3 contribution falling off faster the more the strength of the UV background is increased. This behaviour is again a result of a change in the behaviour of the H + abundance at high densities. In this case, the H + abundance falls off faster at high density when the UV field strength is increased. This occurs because the ultraviolet background maintains a higher Li + fraction in the gas than in our reference calculation. Because the Li + abundance is larger, it contributes more free electrons to the gas, and so the fraction of H + 3 ions that are destroyed by dissociative recombination in the high density regime becomes larger. Consequently, fewer of the H + ions destroyed by reaction RA18 are recycled by reaction TR17, and so the H + and H + 3 abundances fall off more rapidly than in our reference model. A similar effect is not seen in the optically thin runs because it is more than offset by the effects of the higher gas temperature, which increases the rate of reaction TR17, thereby decreasing the fraction of H + 3 ions that are destroyed by dissociative recombination.

X-rays
X-rays are another form of radiation that can affect the evolution of primordial gas. By providing an additional source of ionization in dense gas, they can promote H2 formation (Haiman, Rees, & Loeb 1997;Haiman, Abel, & Rees 2000;Glover & Brand 2003;Machacek, Bryan, & Abel 2003), in much the same manner as the cosmic rays considered in section 5.3. However, just as in the case of cosmic rays, for X-rays to materially affect the importance of H + 3 cooling, they must be able to penetrate deeply into the collapsing protostellar core, to densities n ∼ 10 10 cm −3 or more, corresponding to column densities N ∼ 10 25 cm −2 or more. As a core with this column density is opaque even to hard X-ray photons, the photon flux required to produce a significant photoionization rate is considerable.
Given reasonable assumptions regarding the shape of any high-redshift hard X-ray background, the dominant contribution to the photoionization rate at n = 10 10 cm −3 comes from X-ray photons with energies E ∼ 3 keV (Yan, Sadeghpour, & Dalgarno 1998). The hydrogen ionization cross-section at this energy is approximately 10 −24 cm 2 , and so for a column density N = 10 25 cm −2 , the gas has an optical depth τ ∼ 10 for a 3 keV photon. Each of the photons that is absorbed is responsible for roughly 100 ionizations, once the effects of secondary ionization are taken into account (Dalgarno, Yan, & Liu 1999). Therefore, a flux FX of 3 keV photons incident on the cloud exterior produces an ionization rate RX given approximately by RX ∼ 100 10 −24 FX∆ν e −10 s −1 , where the second line follows if we assume that ∆ν ∼ 1 keV/h. Inverting this expression, we obtain FX ∼ 9 × 10 −9 RX 10 −17 s −1 photons s −1 cm −2 Hz −1 .
This corresponds to an X-ray background field strength at E = 3 keV of IX = 4.2 × 10 −17 RX 10 −17 s −1 erg s −1 cm −2 Hz −1 , which is orders of magnitude larger than any plausible range of values for the high-redshift X-ray background (see e.g., Glover & Brand 2003). We can therefore rule out an extragalactic X-ray background as the source of the required hard X-ray photons. Furthermore, although higher X-ray fluxes can be maintained close to strong X-ray sources such as miniquasars (Haiman, Abel, & Rees 2000;Kuhlen & Madau 2005), even in this case it is difficult to produce a significant ionization rate. For example, the luminosity at 3 keV of the model miniquasars considered by Kuhlen & Madau (2005), assuming their hardest spectral model, is LX ≃ 10 22 erg s −1 Hz −1 . To see a flux IX = 4.2 × 10 −17 erg s −1 cm −2 Hz −1 from this miniquasar, one must therefore be within a distance of it. Gas this close to the miniquasar would have been strongly processed by the ultraviolet radiation of its progenitor, and is not a promising place to expect to find further star formation. We therefore consider it likely that the hard X-ray flux seen by most collapsing protostellar cores will be far too small to significantly affect the production of H + 3 at high densities. Moreover, even if somehow a sufficiently large flux was produced, we would expect its effects to be very similar to those of the cosmic rays considered in §5.3. Accordingly, we do not consider it necessary or time-efficient to examine the effects of a hard X-ray flux in any greater detail.

Altering the initial density and temperature
In order to verify that our main results are not sensitive to the initial temperature or density assumed in our models, we have performed several calculations with different initial densities or temperatures. In runs N1 and N2, we set ni = 0.03 and 30 cm −3 , respectively, while keeping all of the other input parameters fixed. The effect that this has on the thermal state of the gas is illustrated in Figure 12a, where we compare the temperature evolution in runs N1 and N2 with the evolution in our reference calculation, run REF. It is clear from the Figure that the temperature evolution of the three runs is strongly convergent, in line with previous findings (see e.g., Palla, Salpeter, & Stahler 1983;Omukai 2000). Consequently, it comes as no surprise to find that the contribution that H + 3 makes to the cooling in the three runs is not greatly affected by the choice of ni, as shown in Figure 12b.
In runs T1 and T2, we set Ti = 100 and 10000 K, respectively, and performed a similar comparison, which is illustrated in Figure 13. Again we find that the results of runs T1, T2 and our reference run REF converge well, although the differences in this case are slightly larger than those that occur when ni is varied.

Altering the initial fractional ionization
We have also explored the effect of altering the initial fractional ionization of the gas. In runs X1, X2 and X3, we set the initial H + abundance to 10 −6 , 10 −2 or 1.0, respectively. We also rescaled the initial D + abundances by a similar amount. However, the initial He + abundance was not altered. As Figure 14a demonstrates, altering the initial fractional ionization in this way has a dramatic effect on the temperature evolution of the gas. In run X1, the low abundances of free electrons and of H + delay the formation of H2 and limit the amount that can form. The gas therefore undergoes a period of adiabatic heating that lasts for much longer than in our reference calculation. Furthermore, once enough H2 has formed to cool the gas, the gas temperature remains significantly higher than in run REF. However, the two runs eventually converge at n ∼ 10 10 cm −3 , as at this density three-body processes dominate the formation of H2, and so the H2 abundance, and hence the thermal evolution of the gas, are no longer sensitive to the fractional ionization.
In runs X2 and X3, on the other hand, the enhanced initial ionization has the effect of promoting the formation of H2, and cooling the gas more than in our reference run. Moreover, the cooling provided by this extra H2 is sufficient to lower the gas temperature to a level at which chemical fractionation between HD and H2 becomes highly effective. The resulting boost in the HD abundance allows it to dominate the gas cooling and to cool the gas to a lower temperature than could be reached by H2 cooling alone. Note, however, that in contrast to the results of previous studies (see e.g., Nakamura & Umemura 2002;Nagakura & Omukai 2005;Johnson & Bromm 2006;Yoshida et al. 2007), the gas does not reach the CMB temperature. This is a consequence of the relatively rapid rate of collapse assumed here (c.f. §5.7) and of the revised treatment of H2 cooling used in this work, which tends to render H2 cooling less effective, as explored in more detail in Glover & Abel (2008). In any case, the period of HD dominance lasts for only a short time, as illustrated in Figure 14b. As the HD level populations near their LTE values, HD cooling becomes less effective and the gas starts to warm. As it warms, chemical fractionation becomes less effective and the HD:H2 ratio declines. Once the gas has warmed to T ∼ 200 K, which occurs at a density between 10 5 and 10 6 cm −3 , the amount of HD remaining in the gas is no longer sufficient to maintain HD as the dominant coolant; H2 becomes dominant once more, with HD thereafter relegated to a minor role.
In Figure 14c, we show how the ratio of the H + 3 cooling rate to the total cooling rate varies in these runs. We see that if the initial fractional ionization is lowered, the contribution of H + 3 to the cooling is lowered at n < 10 9 cm −3 and increased at n > 10 9 cm −3 relative to our reference calculation. The lowered importance of H + 3 cooling at low densities is a result of the higher gas temperature: the greater temperature sensitivity of the H2 cooling rate compared to the H + 3 cooling rate makes the former more effective in comparison to the latter as the temperature is raised. At high densities, the gas temperatures converge, and the difference between the runs has a different cause: the lower H2 abundance (discussed at the beginning of this subsection) increases the time required to convert all of the H + to H + 3 . For this reason, the H + abundance in n ∼ 10 9 -10 10 cm −3 gas in the low ionization run is higher than in the reference run, with the result that the H + 3 abundance and the H + 3 contribution to the cooling are also marginally higher. If the initial fractional ionization is increased, an interesting effect occurs. The contribution of H + 3 to the total cooling rate is considerably suppressed at densities n < 10 6 cm −3 and n > 10 8.5 cm −3 compared to the contribution in our reference run, but is slightly enhanced at densities 10 6 < n < 10 8.5 cm −3 . The reduction in the effectiveness of H + 3 at low densities in these runs is a result of the previously noted strong enhancement of the HD cooling rate at n < 10 6 cm −3 , as can clearly be seen by comparing Figures 14b and 14c. On the other hand, the suppression of H + 3 cooling at n > 10 8.5 cm −3 results from the rapid loss of H + in the gas driven by reaction RA18, which occurs more rapidly than in the reference run owing to the greater H2 abundance in these runs. Between these two density regimes, there is a small range of densities in which the H + 3 :H2 ratio remains relatively large, and where the gas temperature is ∼ 300-500 K. At these temperatures, the ratio of the H + 3 to H2 cooling rates is larger than at T = 1000 K, owing to the greater temperature dependence of the H2 cooling rate, but the temperature is not low enough for the HD abundance to be significantly enhanced by chemical fractionation. These conditions are therefore close to ideal for H + 3 cooling, and the fact that even in this case the H + 3 contributes no more than a few percent of the total cooling helps to strengthen our conclusion that it is generally of little or no importance. , and X3 (dotted line), which had initial H + abundances of x H + = 2.2 × 10 −4 , 10 −6 , 10 −2 and 1.0, respectively. (b) As (a), but showing the ratio of the HD cooling rate to the total cooling rate in these four runs. (c) As (b), but for the ratio of the H + 3 cooling rate to the total cooling rate.

Changing the elemental composition
Finally, we have examined the effect of changing the elemental composition of the gas by removing all of the deuterium (run EL1), lithium (run EL2) or both (run EL3). Although not physically realistic, these runs do provide a convenient way to examine the roles that the deuterium and lithium play in the overall thermal evolution of the gas. In Figure 15, we show how the gas temperature evolves in run EL1 in comparison to our reference calculation, run REF. We see that the omission of deuterium has a noticeable effect on the temperature evolution at densities 10 3 < n < 10 6 cm −3 , and a very slight effect on the temperature at densities n > 10 12 cm −3 . Examination of the contribution of HD cooling to the total cooling rate in run REF (plotted as the solid line in Figure 14b) demonstrates that at these densities, HD contributes significantly to the total cooling rate; indeed, at its peak at n ∼ 10 4 cm −3 , it contributes almost a third of the total cooling. At first sight, this result is rather surprising, as it is often assumed that HD cooling is unimportant in primordial gas unless the gas has a large initial fractional ionization, as in runs X2 or X3 discussed above. However, it appears that the conventional wisdom is wrong on this point; our results here are consistent with those of previous studies that have included HD (see e.g., Bromm, Coppi, & Larson 2002;Mizusawa, Omukai, & Nishi 2005), and show that although HD is never the dominant coolant, it does contribute enough to the total cooling rate at densities n ∼ 10 4 -10 5 cm −3 to warrant inclusion in future models of population III star formation. In run EL2, the temperature evolution is essentially the same as in run REF, while in run EL3, the evolution is the same as in run EL1, indicating that lithium does not play a significant role in the cooling of the gas, in agreement with the results presented in §5.1; note that runs EL2 and EL3 are not plotted in Figure 15, as they would not be distinguishable from the existing lines. In Figure 16, we investigate the size of the contribution that H + 3 makes to the total cooling rate in runs EL1, EL2 and EL3; we also plot the result from run REF for purposes of comparison. The most obvious point to note is that in the runs without lithium, namely EL2 and EL3, the H + 3 contribution no longer falls off sharply at high densities, although it  remains too small to be significant. This is easy to understand, given our previous discussion of the chemistry of the gas at high densities and low fractional ionizations (see §5.1). As previously noted, in the absence of lithium, the net rate of removal of H + ions from the gas via reaction RA18 decreases as the fractional ionization decreases, since an increasing fraction of the H + 3 created by reaction RA18 is converted back to H2 and H + by reactions TR17 and CT3, instead of being destroyed by reactions DR4 and DR5. Therefore, the H + removal timescale, t loss , remains longer than the free-fall timescale throughout the simulation, and the H + abundance falls off gradually at high densities; the rapid fall-off that occurs in our reference run once t loss < t ff does not take place. Consequently, the corresponding rapid fall-off in the H + 3 abundance also does not occur, as the H + Figure 17. Ratio of the H + 3 cooling rate to the total cooling rate in runs REF (solid line) and RA (dashed line). These runs used values for k RA18 that differed by four orders of magnitude, and that were taken from Gerlich & Horning (1992) and Stancil et al. (1998), respectively. 5.6 Sensitivity to uncertainties in the chemical rate coefficients 5.6.1 Reaction RA18 As we have already discussed in §3, large uncertainties exist in the rate coefficients of a number of the processes included in our chemical model. Of particular relevance to this paper is the huge uncertainty that appears to exist in the value of the rate coefficient for reaction RA18, the formation of H + 3 by the radiative association of H2 with H + . In Figure 17, we compare the contribution of H + 3 cooling in our reference calculation (solid line), in which we adopt the large Gerlich & Horning (1992) rate coefficient for reaction RA18, with the contribution of H + 3 cooling in run RA, a similar calculation that adopts the smaller  rate coefficient (dashed line). We see that H + 3 is less effective in the latter case, and that its effectiveness also peaks at a later point in the simulation.
However, the reduction in the H + 3 contribution is less than one might expect given the very large difference in kRA18 between the two runs. The reason for this is that although the rate of H + 3 formation by radiative association is strongly suppressed, other H + 3 formation mechanisms remain unaffected. The H + 3 formed in run RA is produced primarily by the familiar reaction with the necessary H + 2 coming mainly from reaction RA3, namely The persistence of a significant H + 3 contribution at later times in run RA than in run REF is a clear consequence of the fact that the rate at which H + ions are removed from the gas by reaction RA18 is smaller in the former run than in the latter. As a result, the familiar rapid fall-off in the H + abundance that occurs once the H + removal timescale, t loss , becomes smaller than the free-fall time (see §5.1) takes place at a later time in run RA than in run REF, and hence the corresponding fall-off in the H + 3 abundance also occurs at a later point in the evolution of the gas.

Reactions TB1 and TB2
Large uncertainties also exist in the rate coefficients for three-body H2 formation (reactions TB1 and TB2). As illustrated in Figure 18, these uncertainties significantly affect the temperature evolution of the gas at densities n > 10 8 cm −3 , particularly the uncertainty in the rate of reaction TB1. The use of larger values for the three-body reaction rates leads to faster production of H2 at high densities, and hence a greater H2 cooling rate. This has the effect of slowing the rise in the gas temperature at these densities, which may affect the ability of the gas to fragment at late times (see Clark, Glover, & Klessen 2008). However, in the present context, the effect of the faster H2 formation rates is to make cooling by H + 3 even less effective at late times than in our reference calculation, as demonstrated in Figure 19, where we examine the effect that varying the rate of both reactions has on the contribution that H + 3 makes to the total cooling rate.

Reaction AD1
The uncertainty in the rate of reaction AD1 discussed in Glover, Savin, & Jappsen (2006) also affects the temperature evolution of the gas. We have examined two cases, runs AR1 and AR2. In run AR1, we use a rate coefficient kAD1 = 0.65 × 10 −9 cm 3 s −1 for reaction AD1, taken from Glover, Savin, & Jappsen (2006), which is a plausible lower limit on the rate coefficient. In run AR2, on the other hand, we use a rate coefficient kAD1 = 5.0 × 10 −9 cm 3 s −1 for reaction AD1, again taken from Glover, Savin, & Jappsen (2006), which is a plausible upper limit.
In Figure 20a, we show how the temperature of the gas evolves in these two runs, as well as in run REF for comparison. It is clear from the figure that the uncertainty in kAD1 has only a slight impact on the temperature evolution of the gas. In Figure 20b we show a similar plot of the ratio of the H + 3 cooling rate to the total cooling rate. Again, the rate coefficient uncertainty has only a small effect. This result is in line with previous work showing that these are unimportant when starting from cold, low ionization initial conditions (Glover, Savin, & Jappsen 2006). If, instead, we start with hot, ionized gas, then the effect of the uncertainties on the temperature evolution is much greater (Glover & Abel 2008). However, even in this case, the largest effects are seen at densities n < ∼ 10 4 cm −3 , far below the densities at which H + 3 could conceivably become important, and so our basic conclusion regarding the unimportance of H + 3 cooling remains unaffected. In run AR1, the value for the rate coefficient of reaction AD1 is chosen so as to minimize H 2 production, while in run AR2, H 2 production is maximized, as discussed in more detail in the text. (b) As (a), but showing how the ratio of the H + 3 cooling rate to the total cooling rate varies in these three models.

Reactions RA20 and CD26
Finally, we have investigated the effects of varying the rates of two of the key reactions involved in the LiH + 2 chemistry: LiH + 2 formation by radiative association of Li + and H2 (reaction RA20) and the collisional dissociation of LiH + 2 by H2 (reaction CD26). As we have already discussed in section 3.1.11, the rate coefficients for both of these reactions are unknown. In our reference model, we adopted values of kRA20 = 10 −22 cm 3 s −1 and kCD26 = 1.0 × 10 −9 exp − 3250 T cm 3 s −1 for the rate coefficients; i.e., a small value for reaction RA20 and a large value for reaction CD26. These choices serve to minimize the role played by LiH + 2 in the chemical evolution of the gas, and hence we considered them to be the most conservative options in the circumstances. In run LP1, we adopted instead a much larger value for the rate coefficient of reaction RA20, kRA20 = 10 −17 cm 3 s −1 , but kept the same value for kCD26 as in our reference model. In run LP2, we used our reference value for kRA20, but adopted a much smaller value for the rate coefficient of reaction CD26: kCD26 = 1.0×10 −13 exp − 3250 T cm 3 s −1 . Finally, in run LP3, we altered both rate coefficients, using the larger value for kRA20 and the smaller for kCD26.
In Figure 21, we show how the contribution of H + 3 cooling to the total cooling rate varies with density in runs LP1, LP2 and LP3, along with run REF for comparison. It is clear that in runs LP1 and LP2, the behaviour is essentially the same as in our reference run. This can be understood if we consider the timescale on which Li + ions are destroyed by the reaction sequence followed by This sequence of reactions removes Li + on a timescale where here fDR is given by fDR = (kDR19 + kDR20 + kDR21)n e − (kDR19 + kDR20 + kDR21)n e − + kCD26nH 2 , and represents the fraction of LiH + 2 ions destroyed by dissociative recombination, rather than by collisional dissociation. In our reference run, at a gas density n = 10 10 cm −3 , the temperature T ≃ 1000 K and the electron abundance x e − ≃ 5 × 10 −11 , and hence fDR ≃ 1.4 × 10 −7 and t loss ≃ 1.4 × 10 19 s, many orders of magnitude longer than the dynamical timescale. Thus, in our reference run, the LiH + 2 chemistry has almost no effect on the Li + abundance. In run LP1, kRA20 is a factor of 10 5 larger than in our reference run, and in run LP2, kCD26 is a factor of 10 4 smaller, and so in both runs, t loss is significantly reduced. However, it still remains far greater than the free-fall timescale, which is ∼ 10 10 s at this density. Thus, in these runs, the LiH + 2 chemistry still has almost no effect. In run LP3, however, where we both increase kRA20 and decrease kCD26, t loss is reduced by a factor of 10 9 , making it ; note that the first three of these lines are not distinguishable in the plot. In run LP1, the rate coefficient for reaction RA20 was increased by a large factor compared to our reference value, while in run LP2, the rate coefficient for reaction CD26 was decreased by a large factor. In run LP3, both changes were made. t loss ∼ 1.4 × 10 10 s at n = 10 10 cm −3 , of the same order of magnitude as the free-fall collapse time. Moreover, as t loss scales with density as t loss ∝ n −1 , while the free-fall time scales as t ff ∝ n −1/2 , t loss becomes smaller than t ff at densities not very much greater than 10 10 cm −3 . Therefore, in this run, the LiH + 2 chemistry does have a noticeable effect on the Li + abundance, reducing it by a factor of roughly fifty in comparison to the reference run by the end of the simulation. This reduction in the Li + abundance reduces the number of free electrons available for destroying H + 3 , and so limits the rate at which its abundance declines at very high densities, much as in runs performed without any lithium (c.f. section 5.5.3). Nevertheless, it is clear from Figure 21 that this change in the lithium chemistry does not change our basic results: H + 3 cooling remains ineffective, albeit somewhat less ineffective at high densities than in our reference run.

Sensitivity to the details of the dynamical model
A major limitation of our current study is the highly simplified dynamical treatment that we use in our one-zone model. In all of the calculations that we have presented so far, we have assumed that the gas is collapsing gravitationally at the free-fall rate. However, it is well known from more detailed three-dimensional hydrodynamical models (e.g., Abel, Bryan, & Norman 2002;Yoshida et al. 2006) that in realistic primordial clouds, the collapse speed is significantly slower than the free-fall rate, owing to the non-negligible gas pressure. As a result, the gas takes longer to evolve than assumed here, and the impact of compressional heating is also somewhat smaller. A crude way of taking this into account in a one-zone calculation is to artifically slow down the collapse of the gas. In other words, instead of assuming that the density evolves according to the standard free-fall relationship where t ff is the free-fall time, we can instead assume that with η < 1. The effect of this change is to lengthen the time taken for the gas to collapse to any given density by a factor (1/η). In Figure 22a we show the effect that slowing the collapse in this fashion has on the temperature evolution of the gas by comparing the results of three runs, DYN1, DYN2 and DYN3, with η = 0.6, η = 0.3 and η = 0.1, respectively, with the results of our reference run REF. We see that reducing η leads to a reduction in the temperature of the gas throughout the run, a simple consequence of the reduction in the compressional heating rate. Interestingly, in the η = 0.3 and η = 0.1 models, the reduced heating allows the gas to cool to temperatures low enough for chemical fractionation to strongly enhance the HD fraction, allowing HD cooling to further cool the gas down to temperatures close to TCMB. The fact that this effect is not seen in more realistic hydrodynamical models (e.g., Bromm, Coppi, & Larson 2002) suggests that in these models we are overestimating the extent to which gas pressure slows the collapse, at least at early times. In Figure 22b, we show how reducing η affects the contribution that H + 3 cooling makes to the total cooling rate. In the η = 0.6 run, H + 3 cooling is slightly more effective than in the reference run at densities 10 5 < ∼ n < ∼ 10 8 cm −3 . This is a consequence of the slightly lower temperature of the gas at these densities in run DYN1 compared to run REF, which decreases the H2 cooling rate more than the H + 3 cooling rate. At lower densities, H + 3 cooling in run DYN1 is slightly less effective than in run REF, as the lower gas temperature makes HD cooling more effective at these densities in the former run than in the latter. At higher densities, H + 3 cooling becomes far less effective in run DYN1 than in run REF owing to a more rapid falloff in the H + 3 abundance, as illustrated in Figure 22c. As in earlier runs, the reason for this rapid falloff is that the timescale for the removal of H + from the gas by conversion to H + 3 followed by destruction of the H + 3 by dissociative recombination, becomes shorter than the dynamical timescale of the gas at these densities. In most of the runs that we have studied, the high gas temperature at these densities allows reaction TR17 to interfere with this process by converting most of the H + 3 to H + 2 , following which reaction CT3 restores the original proton to the gas. As we saw in §5.1, the effect of this is to lengthen the time required to remove all of the H + , which delays the precipitous falloff until n ∼ 10 11 cm −3 . However, in run DYN1, the lower gas temperature means that reaction TR17 is less effective, and so the delay is much shorter. Moreover, the dynamical timescale itself is longer. Consequently, the rapid falloff occurs at a lower density.
A similar kind of behaviour is seen in runs DYN2 and DYN3. However, in these runs the much lower temperature at low densities renders HD cooling dominant for significantly longer, while the lower temperature of the gas at high densities, plus the longer dynamical timescale, allow the rapid falloff in the H + 3 abundance to occur sooner. Consequently, H + 3 cooling never becomes significant in these runs.

CONCLUSIONS
We have examined the contribution that H + 3 cooling makes to the total cooling rate of gravitationally collapsing primordial gas in a wide range of different models using a newly-developed H + 3 cooling function along with the most detailed model of primordial gas chemistry published to date. Our results demonstrate that in general H + 3 cooling is not important, although it comes close to being so at densities n = 10 7 -10 9 cm −3 , contributing at its peak a few percent of the total amount of cooling. We come to this conclusion despite making several assumptions (regarding the collapse rate of the gas, the formation rate of H + 3 by radiative association, and the collisional excitation rate of its excited vibrational states) that favour H + 3 cooling, and thus we have confidence that our conclusion is robust.
As H + 3 comes so close to being important, it is instructive to examine why it ultimately fails to dominate. This can be ascribed to a combination of two main effects. First, at high densities (n > ∼ 10 8 cm −3 ), the gas temperature becomes high enough to make the endothermic reaction TR17 the most important H + 3 destruction mechanism, which significantly suppresses the H + 3 abundance at these densities. Second, the formation rate of H + 3 is strongly suppressed by the rapid removal of H + from the gas at densities n > ∼ 10 11 cm −3 . At these densities, the fractional ionization of the gas is so low that the main loss route for the H + is conversion to H + 3 , followed by H + 3 dissociative recombination, and once the timescale for H + removal via this combination of reactions becomes short compared to the free-fall time, the H + abundance decreases by orders of magnitude within a short space of time, effectively switching off the formation of H + 3 . Moreover, destruction of H + 3 by dissociative recombination remains effective despite the fall-off in x H + thanks to the contribution of electrons from ionized lithium, Li + , which for n > 3 × 10 8 cm −3 is the most abundant positive ion in the gas.
In our study, the only situation in which we found H + 3 cooling to be important is if the gas is illuminated by a strong flux of cosmic rays or X-rays. If the incident flux is strong enough to produce an ionization rate > ∼ 10 −18 s −1 at densities n 10 10 cm −3 , then the high density H + 3 abundance can be significantly increased, and H + 3 can even become the dominant coolant. However, the necessary flux of cosmic rays or hard X-rays is difficult to produce in the high-redshift Universe. As the estimates in § § 5.3 and 5.4.2 demonstrate, the required flux is orders of magnitude greater than the size of any plausible extragalactic background, and will only be achieved within gas that is very close to a local source (i.e., within 5-10 pc). However, gas that is this close to a supernova remnant or miniquasar will have been strongly affected by radiative feedback from the progenitor star and so is not a promising place to expect to find ongoing population III star formation.
Furthermore, even if the ionization rate is high enough to make H + 3 an important or dominant coolant at high densities, the effect of H + 3 cooling on the thermal evolution of the gas remains relatively small; the difference it makes to the temperature evolution at n > 10 8 cm −3 is smaller than the error introduced by the uncertainty in the three-body H2 formation rate (reaction TB1). Finally, our model has also allowed us to explore the effects of cooling from the other minor ionic and molecular species present in the gas (e.g., H + 2 , H2D + , LiH, etc.). Despite making rather optimistic assumptions regarding the cooling from these species, we find that they are orders of magnitude less effective than H + 3 at cooling high density gas, and hence are never significant.