Skip to main content

The numerical frontier of the high-redshift Universe

Abstract

The first stars are believed to have formed a few hundred million years after the big bang in so-called dark matter minihalos with masses \(\sim 10^{6}\mbox{ M}_{\odot}\). Their radiation lit up the Universe for the first time, and the supernova explosions that ended their brief lives enriched the intergalactic medium with the first heavy elements. Influenced by their feedback, the first galaxies assembled in halos with masses \(\sim10^{8}\mbox{ M}_{\odot}\), and hosted the first metal-enriched stellar populations. In this review, I summarize the theoretical progress made in the field of high-redshift star and galaxy formation since the turn of the millennium, with an emphasis on numerical simulations. These have become the method of choice to understand the multi-scale, multi-physics problem posed by structure formation in the early Universe. In the first part of the review, I focus on the formation of the first stars in minihalos - in particular the post-collapse phase, where disk fragmentation, protostellar evolution, and radiative feedback become important. I also discuss the influence of additional physical processes, such as magnetic fields and streaming velocities. In the second part of the review, I summarize the various feedback mechanisms exerted by the first stars, followed by a discussion of the first galaxies and the various physical processes that operate in them.

1 Introduction

The formation of the first stars marked a fundamental transition in the history of the Universe. They initiated the transformation of the homogeneous intergalactic medium (IGM) to one filled with the rich structure we observe today. Ending the so-called ‘cosmic dark ages’, when the Universe contained no visible light, they lit up the Universe at redshifts \(z\gtrsim20\) (Bromm and Larson 2004; Glover 2005, 2013; Bromm 2013). They formed at the center of dark matter (DM) ‘minihalos’ with virial masses \(M_{\mathrm{vir}}\sim10^{6}\mbox { M}_{\odot}\), which are the smallest building blocks in the hierarchy of galaxy formation. These objects accrete the pure hydrogen and helium gas forged in the Big Bang, and after continued cooling and contraction form a stellar embryo that begins to accrete from the surrounding gas cloud.

Since the virial temperature of minihalos is not high enough to activate atomic hydrogen cooling, the gas can only cool via molecular hydrogen (H2). The importance of H2 for the cooling of low-mass gas clumps that condense out of the expanding Universe was recognized in the late 1960’s (Saslaw and Zipoy 1967; Peebles and Dicke 1968; Hirasawa 1969; Matsuda et al. 1969; Takeda et al. 1969). Later on, simplified one-zone models accounted for the dynamics of collapsing gas clouds next to the radiative cooling of the gas (Yoneyama 1972; Hutchins 1976; Silk 1977, 1983; Carlberg 1981; Kashlinsky and Rees 1983; Palla et al. 1983; Carr et al. 1984; Izotov and Kolesnik 1984; Couchman and Rees 1986; Susa et al. 1996; Uehara et al. 1996; Tegmark et al. 1997; Nishi and Susa 1999; Vasiliev and Shchekinov 2003, 2005b). The first one-dimensional calculations of the simultaneous collapse of DM and gas were carried out in the context of the cold dark matter (CDM) paradigm (Haiman et al. 1996b; Omukai and Nishi 1998; Nakamura and Umemura 1999), while three-dimensional simulations had to await improvements in numerical simulation techniques in the late 1990’s (Abel et al. 1998, 2000, 2002; Bromm et al. 1999, 2002).

One of the main results of these studies is that the minimum temperature to which the gas can cool via H2 lines is more than an order of magnitude higher than in present-day star formation regions, resulting in a greatly increased Jeans mass. Since the accretion rate is directly related to the Jeans mass, Population III stars are expected to have masses of the order of \(M_{*}\sim100\mbox { M}_{\odot}\). Following the introduction of this ‘standard model’ of primordial star formation, recent work has focused on refining this picture. In particular, the influence of fragmentation, protostellar interactions, radiation, and magnetic fields were initially neglected. These may have a substantial effect on the collapse of the gas. In this review, I summarize the progress made on these topics since the advent of the first three-dimensional simulations of primordial star formation.

Despite the complications brought about by these processes, it is likely that Population III stars were much more massive than present-day stars. They are therefore strong emitters of ultraviolet (UV) radiation, which heats the IGM and begins the process of reionization (Loeb and Barkana 2001; Ciardi and Ferrara 2005; Fan et al. 2006; Stiavelli 2009). Their radiation also dissociates H2 on cosmological scales and may suppress star formation until halos massive enough to activate atomic hydrogen cooling assemble (Ciardi and Ferrara 2005). Next to radiative feedback, the supernova (SN) explosions of massive Population III stars exert mechanical as well as chemical feedback (e.g. Wise and Abel 2008b; Greif et al. 2010). The ejected metals facilitate the transition to Population I/II star formation by allowing the gas to cool to much lower temperatures than would otherwise be possible (e.g. Omukai 2000; Bromm et al. 2001a).

Due to the intricate nature of feedback, the formation of the first galaxies in so-called ‘atomic cooling halos’ is much more complicated than the formation of the first stars in minihalos (Bromm and Yoshida 2011; Loeb and Furlanetto 2013). Fully self-consistent simulations that predict the properties of the first galaxies are not yet available. Nevertheless, recent work has shown that supersonic turbulence is one of the key factors governing the properties of the first galaxies (e.g. Wise and Abel 2007a; Greif et al. 2008). Radiative and chemical feedback from stars in progenitor halos influence their formation as well (e.g. Wise and Abel 2008b; Greif et al. 2010). The degree of sophistication of first galaxy simulations has continued to increase, with the most recent studies including star formation and feedback recipes that rival those of present-day galaxy formation simulations (e.g. Wise et al. 2012).

The organization of this review is as follows. In Section 2, I give a brief introduction to structure formation in the high-redshift Universe, followed by a description of the collapse of gas in minihalos. I then review the influence of fragmentation on the accretion phase, the evolution of the nascent protostellar system, and the effects of protostellar radiation on the initial mass function (IMF) of the first stars (Section 3). In Section 4, I discuss additional physics that may affect the formation of the first stars, including hydrogen deuteride (HD) cooling, magnetic fields, cosmic rays, streaming velocities, DM annihilation, and alternative cosmologies. Finally, in Section 5 I discuss the radiative, mechanical and chemical feedback exerted by Population III stars, next to the progress made with respect to the formation of the first galaxies. I focus on theoretical work, with a particular emphasis on numerical simulations, but briefly mention neglected processes and empirical signatures in Section 6. The summary is presented in Section 7. All distances are quoted in proper units, unless noted otherwise.

Finally, a reference to related reviews is in order. A general description of structure formation in the early Universe may be found in Barkana and Loeb (2001), Loeb (2008), and Wiklind et al. (2013). The high-redshift IGM is discussed in Barkana and Loeb (2007) and Meiksin (2009), and the effects of the relative velocity offset between DM and baryons may be found in Fialkov (2014). The formation of the first stars in minihalos is reviewed in Bromm and Larson (2004), Glover (2005, 2013), and Bromm (2013), while the properties of the first galaxies are described in Bromm and Yoshida (2011), Johnson (2013), and Loeb and Furlanetto (2013). Less focused reviews of star formation at high redshifts may be found in Bromm et al. (2009) and Loeb (2010). Finally, feedback by Population III is summarized in Ciardi and Ferrara (2005).

2 First stars: the initial collapse

2.1 Structure formation

On the largest scales, the Universe appears nearly uniform and isotropic. However, the presence of stars and galaxies indicates that below a certain scale the Universe must have begun to deviate from its uniformity. It is believed that these structures grew from infinitesimally small perturbations seeded by quantum fluctuations in the very early Universe. Within variants of the CDM model, where the mass density of the Universe is dominated by DM, the matter overdensity \(\delta= (\rho-{\bar{\rho}} )/{\bar{\rho}}\), where ρ is the local mass density and \({\bar{\rho}}\) the mean density of the Universe, grows in proportion to the scale factor \(a=1/ (1+z )\), where z denotes redshift. Once the overdensity becomes of order unity, the associated region decouples from the expanding background Universe and collapses under its own gravity. The collapse may occurs simultaneously in one, two, or three dimensions. Structures that collapse in one dimension are termed ‘sheets’, the collapse of two sheets results in a ‘filament’, and DM ‘halos’ form at the intersection of filaments. The likelihood of formation decreases with increasing dimensionality, such that there are many more sheets than halos. However, halos provide the deepest potential wells and are therefore the sites of star and galaxy formation.

2.2 Dark matter halos

One of the most important characteristics of the CDM paradigm is the hierarchical, ‘bottom-up’ nature of collapse. Increasingly massive halos form via accretion and merging of low-mass halos. The smallest collapse scale is set by the free-streaming length, which is about 1 au if DM consists of a weakly interacting massive particle (WIMP; Boyanovsky et al. 2008). Following the collapse of a DM halo, it achieves virial equilibrium within a few dynamical times, with the virial density given by \(\rho_{\mathrm{vir}}=\Delta_{\mathrm{vir}}{\bar{\rho}}\), where \(\Delta_{\mathrm{vir}}\simeq200\) (e.g. Barkana and Loeb 2001). The balance between gravitational potential energy and kinetic energy allows a virial velocity to be defined, and is given by \(v_{\mathrm{vir}}^{2}=GM_{\mathrm{vir}}/R_{\mathrm{vir}}\), where G is the gravitational constant, \(M_{\mathrm{vir}}\) the virial mass, and \(R_{\mathrm{vir}}\) the virial radius of the halo. The latter is related to the mass of the halo via \(R_{\mathrm{vir}}^{3}=3M_{\mathrm{vir}}/ (4\pi\Delta_{\mathrm{vir}}{\bar{\rho}} )\). Halos with virial masses up to \({\sim}10^{6}\mbox { M}_{\odot}\) are commonly referred to as minihalos, since they are massive enough to cool and host star formation. A \(2\mbox{-}3\sigma\) peak typically forms at \(z\sim20\), and their virial radius may be conveniently written as:

$$ R_{\mathrm{vir}}\simeq100\mbox { pc} \biggl(\frac{M_{\mathrm{vir}}}{10^{6}\mbox { M}_{\odot}} \biggr)^{1/3} \biggl(\frac{1+z}{20} \biggr)^{-1}. $$
(1)

A number of studies have systematically investigated the properties of minihalos, finding that they are usually denser and more clustered than their high-mass counterparts (Jang-Condell and Hernquist 2001; Heitmann et al. 2006; Lukić et al. 2007; Reed et al. 2007; Cohn and White 2008; Davis and Natarajan 2009, 2010; Davis et al. 2011; de Souza et al. 2013a; Sasaki et al. 2014).

The DM sets the gravitational potential for the gas, which virializes within the halo and heats to the virial temperature \(T_{\mathrm{vir}}=\mu m_{\mathrm{H}}v_{\mathrm{vir}}^{2}/k_{\mathrm{B}}\), where μ is the mean molecular weight, \(m_{\mathrm{H}}\) the mass of the hydrogen atom, and \(k_{\mathrm{B}}\) Boltzmann’s constant. In terms of a typical minihalo, the virial temperature may be written as:

$$ T_{\mathrm{vir}}\simeq2\times10^{3}\mbox { K} \biggl( \frac{M_{\mathrm{vir}}}{10^{6}\mbox { M}_{\odot}} \biggr)^{2/3} \biggl(\frac{1+z}{20} \biggr). $$
(2)

In contrast to the pressure-less DM, a minimum scale exists below which the gas cannot collapse. This scale is approximately given by the Jeans length, \(\lambda_{\mathrm{J}}=c_{\mathrm{s}}t_{\mathrm{ff}}\), where \(c_{\mathrm{s}}\) is the sound speed of the gas, and \(t_{\mathrm{ff}}=\sqrt{3\pi/ (32G\rho )}\) the free-fall time. The corresponding Jeans mass is given by \(M_{\mathrm{J}}=\pi\rho\lambda_{\mathrm{J}}^{3}/6\), and may be derived from the density and temperature of the IGM. Since the gas temperature is coupled to the cosmic microwave background (CMB) by Compton scattering at \(z\gtrsim100\), the cosmological Jeans mass initially remains constant at \(M_{\mathrm{J}}\simeq10^{5}\mbox { M}_{\odot}\) (Bromm and Larson 2004). At lower redshifts, the gas expands adiabatically, and the Jeans mass is given by \(M_{\mathrm{J}}\simeq10^{4}\mbox { M}_{\odot}[ (1+z )/20 ]^{3/2}\). This is many orders of magnitude larger than the minimum DM halo mass. More careful calculations that use perturbation theory to compute the growth of density fluctuations define a filtering mass below which the gas content in a halo is significantly suppressed with respect to the cosmic mean (e.g. Gnedin and Hui 1998; Naoz and Barkana 2007). These calculations find that the filtering mass remains nearly constant at a few times 104 M at \(z\lesssim100\). Even though this is the minimum mass of a halo at which the gas can contract by a substantial amount, it is not yet a sufficient criterion for star formation.

2.3 H2 formation and cooling

For the collapse of the gas to continue beyond the formation of a hydrostatic object, it must be able to radiate away its thermal energy. One of the most efficient coolants in primordial gas is collisional excitation cooling of atomic hydrogen, also termed Ly-α cooling. This coolant is most effective at temperatures \(\simeq10^{4}\mbox { K}\), where the first excited state of atomic hydrogen begins to be populated. However, halos with masses greater than the filtering mass, but below approximately \(M_{\mathrm{vir}}=5\times10^{7}\mbox { M}_{\odot}\), which corresponds to a virial temperature of \({\simeq}10^{4}\mbox { K}\) at \(z=10\), must rely on another coolant. Saslaw and Zipoy (1967) and Peebles and Dicke (1968) realized that molecular hydrogen (H2) may be such a coolant.

The rotational and vibrational states of H2 in the electronic ground state are excited by collisions with other particles, which decay and allow the gas to cool. The most important transitions operate in the temperature range \(200\mbox { K}\lesssim T\lesssim5\text{,}000\mbox { K}\). At higher temperatures, Ly-α cooling becomes dominant. Despite its importance at high redshifts, the H2 molecule is not a particularly efficient coolant. Since the hydrogen molecule is symmetric, it does not have a permanent electric dipole moment, resulting in correspondingly lower transition probabilities. In addition, the existence of ortho and para states for hydrogen rules out the lowest energy transition \(J=1\rightarrow0\), where J is the rotational quantum number. The least energetic transition with a non-negligible probability of occurring is the \(J=2\rightarrow0\) transition, which corresponds to a temperature of \({\simeq}500\mbox { K}\). In practice, the Maxwellian tail of the velocity distribution allows the gas to cool to \({\simeq}200\mbox { K}\). The absence of a dipole moment in the H2 molecule also rules out its simplest formation channel: the direct association of two hydrogen atoms. The excess energy of the collision cannot be radiated away quickly enough, and so the intermediate system decays again (Gould and Salpeter 1963).

For this reason, alternative formation channels have been considered (Haiman et al. 1996b; Abel et al. 1997; Galli and Palla 1998; Stancil et al. 1998). For the purpose of primordial star formation, the most important formation channel is via the intermediary reaction of radiative association of free electrons with neutral hydrogen atoms (McDowell 1961; Peebles and Dicke 1968):

$$ \mathrm{H}+\mathrm{e}^{-}\rightarrow\mathrm{H}^{-}+\gamma. $$
(3)

Associative detachment of the H atoms with neutral hydrogen atoms then results in the formation of H2:

$$ \mathrm{H}^{-}+\mathrm{H}\rightarrow\mathrm{H}_{2}+\mathrm{e}^{-}. $$
(4)

One of the most important parameters that governs the amount of H2 that can be formed is the electron abundance. It decreases from its relic abundance by recombining with ionized hydrogen:

$$ \mathrm{H}^{+} +\mathrm{e}^{-}\rightarrow\mathrm{H}+\gamma. $$
(5)

Most of the H2 therefore forms within a recombination time. A straightforward calculation shows that the asymptotic H2 abundance depends primarily on the ratio of the H formation to the recombination rate constant (Glover 2013). If the recombination rate is significantly smaller than the radiative association rate, more electrons remain to form H. For typical values of the rate equations, the asymptotic H2 abundance scales approximately as \(y_{ \mathrm{H}_{2}}\propto T^{1.5}\), and for \(T=1\text{,}000\mbox { K}\) lies in the range of \(y_{ \mathrm{H}_{2}}\simeq10^{-4}\mbox{-}10^{-3}\) (Tegmark et al. 1997). This leads to the somewhat counterintuitive result that a more massive halo with a higher virial temperature forms more H2, which allows the gas to cool more efficiently.

The H2 abundance necessary to cool the gas and facilitate its collapse to high densities may be obtained by comparing the cooling time with the Hubble time. The latter estimates the time scale on which the virial temperature of a halo changes significantly (Rees and Ostriker 1977; Silk 1977). A straightforward calculation shows that for a virial temperature of about \(1\text{,}000\mbox { K}\), enough H2 is produced to cool the gas (Glover 2013). This value does not depend strongly on redshift, such that the minimum halo mass required for efficient cooling may be obtained from equation (2):

$$ M_{\mathrm{cool}}\simeq3.5\times10^{5}\mbox { M}_{\odot}\biggl(\frac {1+z}{20} \biggr)^{-3/2}. $$
(6)

This is typically an order of magnitude higher the filtering mass, showing that not all halos that are massive enough to allow the gas to collapse are also massive enough to facilitate cooling and star formation.

It is important to note that the above derivation of the cooling mass is only approximately valid, since in reality the gas is constantly heated by minor mergers (Yoshida et al. 2003a). Furthermore, even though the first stars typically formed at \(z\sim20\), they are not the very first stars in the observable Universe. Instead, these are believed to have formed at \(z\simeq50\mbox{-}65\) (Gao et al. 2005; Naoz et al. 2006).

2.4 Jeans instability

During the initial stages of the collapse, the level populations of H2 are not yet populated according to local thermal equilibrium (LTE), and the cooling rate scales as \(\Lambda_{ \mathrm{H}_{2}}\propto n_{\mathrm{H}}^{2}\), where \(n_{\mathrm{H}}\) is the number density of hydrogen nuclei. Once enough H2 has formed, a runaway process ensues in which cooling allows the central gas cloud to contract, which results in even more cooling. This collapse phase ends when collisions between H2 molecules and other species become so frequent that the ro-vibrational levels of H2 are populated according to LTE. In this case, the collisional excitation rate is approximately balanced by the collisional de-excitation rate, such that the cooling rate scales approximately as \(\Lambda_{ \mathrm{H}_{2}}\propto n_{\mathrm{H}}\). The cooling time thus remains slightly larger than the free-fall time, and the collapse rate decreases.

The transition from non-LTE to LTE occurs approximately at a density of \(n_{\mathrm{H}}\simeq n_{ \mathrm{H}, \mathrm{crit}}=10^{4}\mbox { cm}^{-3}\), where the gas has cooled to \({\simeq}200\mbox { K}\). This phase in the collapse of the gas has been termed the ‘loitering phase’, since the collapse momentarily slows while the central gas cloud continues to accrete mass (Bromm et al. 2002). The Jeans mass at this characteristic density and temperature is \({\gtrsim}100\mbox { M}_{\odot}\). Once this mass has been accreted, the cloud collapses under its own gravity. Since the H2 cooling rate is a strong function of temperature, the cooling time \(t_{\mathrm{cool}}\) in the LTE regime adjusts to the free-fall time. Equating these two time scales implies \(\rho^{-1/2}\propto T^{1-\alpha }\), where \(\alpha\simeq4\) is the power-law exponent of the cooling function. Writing this in the form \(T\propto\rho^{\gamma_{\mathrm{eff}}}\), where \(\gamma_{\mathrm{eff}}\) is the effective equation of state, one finds \(\gamma_{\mathrm{eff}}=7/6\). The temperature therefore increases very slowly with increasing density for \(n_{\mathrm{H}}\gtrsim n_{ \mathrm{H}, \mathrm{crit}}\). In a near-isothermal, Jeans-unstable cloud the accretion rate is approximately given by \({\dot{M}}=M_{\mathrm{J}}/t_{\mathrm{ff}}\propto T^{3/2}\). Since the temperature of present-day star formation regions is about 10 K, this simple physical argument implies that Population III stars were typically 100 times more massive than Population I/II stars.

2.5 Three-body H2 formation

Following the nearly isothermal collapse of the cloud to densities \(n_{\mathrm{H}}\simeq10^{8}\mbox { cm}^{-3}\), three-body reactions begin to convert the mostly atomic gas to molecular hydrogen (Palla et al. 1983). The most important formation reaction is

$$ \mathrm{H}+\mathrm{H}+\mathrm{H}\rightarrow\mathrm{H}_{2}+\mathrm{H}, $$
(7)

with a smaller contribution from reactions involving H2 molecules and helium atoms. The rapidly increasing H2 fraction results in a similarly rapid increase in the H2 line cooling rate. However, since each formation process is accompanied by the release of the binding energy of the H2 molecule of \({\simeq}4.48\mbox { eV}\), the cooling is offset by chemical heating. The net effect is a mild drop in temperature as the H2 fraction approaches unity. The resulting thermal instability of the gas caused by a chemical transition corresponds to the chemothermal instability envisaged by Yoneyama (1973). Previous studies have shown that it may trigger gravitational instability and fragmentation (Sabano and Yoshii 1977; Ripamonti and Abel 2004; Nakamura and Umemura 1999; Yoshida et al. 2006). Indeed, in the simulations of Turk et al. (2009) and Greif et al. (2013), a subset of the clouds fragmented into two distinct clumps, although the strength of the instability in Greif et al. (2013) was likely overestimated due to the approximate treatment of the optically thick H2 line cooling rate (see Section 2.6). The further evolution of the clumps cannot be reliably extrapolated from the state of the cloud at the instant they fragment, but it appears that in some halos a wide binary system may form.

Until recently, one of the major caveats in the chemical evolution of primordial gas clouds was the large uncertainty in the three-body formation rate coefficient. At \(1\text{,}000\mbox { K}\), the published rates differed by up to two orders of magnitude (Abel et al. 2002; Palla et al. 1983; Flower and Harris 2007; Glover 2008). The influence of the three-body formation rate on the properties of primordial gas clouds was investigated by Turk et al. (2011) and Bovino et al. (2014c). They found that radially averaged quantities as well as the morphologies of the clouds varied significantly between simulations with different rates, and are comparable to the differences for different initial conditions. This uncertainty may have recently been alleviated by the quantum-mechanical calculations of Forrey (2013). They found that the three-body formation rate coefficient at the relevant temperatures is approximately two times lower than that of the commonly employed Glover (2008) rate.

2.6 Optically thick H2 line cooling

Another important transition occurs at densities \(n_{\mathrm{H}}\gtrsim 10^{10}\mbox { cm}^{-3}\), where the gas becomes optically thick to H2 line emission. For each of the 200 ro-vibrational transitions that are important at the relevant densities and temperatures, thermal Doppler broadening dominates the frequency dependence of the emitted radiation. The cross section for each line shows a similar frequency dependence, but in addition depends on the relative velocity along the line of sight. Until recently, this complicated radiative transfer problem has only been solved accurately in one-dimensional calculations (Omukai et al. 1998; Omukai and Nishi 1998; Ripamonti et al. 2002). These studies found that the wings of the lines allow the radiation to escape more efficiently than in the case of a grey opacity. In a spherically symmetric calculation, Ripamonti et al. (2002) and Ripamonti and Abel (2004) found that the escape fraction, which is defined as the probability of a photon to escape from the cloud, may be fit by \(f_{\mathrm{esc}}=\min [ (n_{\mathrm{H}}/n_{ \mathrm{H}, \mathrm{line}} )^{-0.45}, 1 ]\), where \(n_{ \mathrm{H},\mathrm{line}}=6\times10^{9}\mbox { cm}^{-3}\).

This fit has also been used in three-dimensional simulations (Turk et al. 2009, 2010, 2011, 2012), due to the high computational cost of multi-frequency line transfer calculations in three dimensions. The accuracy of this treatment depends on how spherically symmetric the cloud is, and on how well the radially averaged profiles agree with those found in Ripamonti et al. (2002). In particular, the functional form of the fit depends on the detailed chemical, thermal, and dynamical evolution of the cloud, which in turns depends on the various rate equations employed (Turk et al. 2011). According to the simulations of Turk et al. (2009), the requirement of spherical symmetry is relatively well satisfied, while simulations that use the smoothed particle hydrodynamics method tend to find a more pronounced disk component (e.g. Yoshida et al. 2006; Clark et al. 2011b). The moving-mesh simulations of Greif et al. (2013) found clouds with a morphology similar to those of Turk et al. (2009). Only a few showed substantial deviations from spherical symmetry. In cases where these deviations exist, Hirano and Yoshida (2013) found that a simple density-dependent fitting function does not yield accurate results.

Another method to estimate the escape fraction is the Sobolev method (Sobolev 1960). In a cloud with a uniform velocity gradient, the escape fraction is given by \(f_{\mathrm{esc}}= [1-\exp (-\tau ) ]/\tau\), where \(\tau=\alpha L_{\mathrm{Sob}}\), \(L_{\mathrm{Sob}}\) is the Sobolev length, and α the absorption coefficient. The Sobolev length is given by \(L_{\mathrm{Sob}}=v_{\mathrm{therm}}/\vert \mathrm{d}v_{\mathrm{r}}/\mathrm{d}r\vert \), where \(v_{\mathrm{therm}}\) is the thermal velocity, and \(\mathrm{d}v_{\mathrm{r}}/\mathrm{d}r\) the velocity gradient. The Sobolev length estimates the scale on which a line is Doppler-shifted out of resonance. This method was first used in the simulations of Yoshida et al. (2006), where the average escape fraction along the three principal axes of the computational domain was computed. In subsequent studies, this was replaced by the average velocity gradient, and the Sobolev length was limited to be smaller than the Jeans length for very small velocity gradients (Clark et al. 2011a, 2011b; Greif et al. 2011a, 2012, 2013). In general, the Sobolev method is only valid if the Sobolev length is small compared to the scales on which the properties of the gas change significantly. This is not the case in primordial gas clouds, where transonic turbulence is present and the variation in the radial velocity is comparable to the sound speed of the gas. Nevertheless, Turk et al. (2011) found that the escape fractions obtained with the different methods disagree only by a factor of a few, while Yoshida et al. (2006) and Hirano and Yoshida (2013) found agreement within a factor of two. However, the results of Turk et al. (2011) were obtained using different hydrodynamical schemes for the two methods, such that it is unclear to what extent the discrepancy is related to the treatment of the optically thick cooling rate.

The first self-consistent treatment of H2 line transfer in a three-dimensional simulation was presented in Greif (2014). This study used a multi-line, multi-frequency ray-tracing scheme to follow the propagation of H2 line radiation through the computational domain (see Figure 1). The spherically averaged escape fraction agreed relatively well with the fit of Ripamonti and Abel (2004), while the Sobolev method yielded escape fractions that were up to an order of magnitude higher. This discrepancy is due to the fact that the Sobolev approximation is not valid in primordial gas clouds. As a result, the gas temperature in Greif (2014) is systematically higher than in previous studies that used the Sobolev method, and the cooling instability found in Greif et al. (2013) also disappears. The reduced temperature results in a somewhat reduced H2 fraction at densities \(n_{\mathrm{H}}\simeq10^{15}\mbox { cm}^{-3}\), which is in agreement with Turk et al. (2012). An effect that is not captured in the previous approximate methods is the diffusion of radiation through the cloud, which tends to smooth out density perturbations. As a result, the central gas cloud becomes increasingly spherically symmetric as it collapses to higher densities. The fitting function also does not account for changes in the geometry of the cloud, which are expected once the collapse stalls and a disk forms. It neglects the complex dependence of the escape fraction on the properties of the gas, such as the density, velocity, and temperature profiles, and thus an implicit dependence on the various chemical and thermal rate equations employed.

Figure 1
figure 1

Number density of hydrogen nuclei, temperature, and escape fraction of \(\pmb{\mbox{H}_{2}}\) line photons in the central 200 au of a minihalo, shown for various treatments of the optically thick \(\pmb{\mbox{H}_{2}}\) line cooling rate. The top row shows the outcome for the density-dependent fit of Ripamonti and Abel (2004), the middle row for the Sobolov method, and the bottom row for the self-consistent radiative transfer calculation of Greif (2014). The Sobolev method greatly overestimates the escape fraction at high densities, whereas the fitting function shows better agreement with the ray-tracing solution, even though the gas fragments in this particular simulation. Adapted from Greif (2014).

Unfortunately, it is not yet computationally feasible to continue the radiation hydrodynamics simulations of Greif (2014) beyond the initial collapse of the gas (e.g. Greif et al. 2012). This is in part due to the large number of opacity calculations necessary for an accurate integration, which is given by \(N_{\tau}=N_{\mathrm{cells}}^{4/3}N_{\mathrm{rays}}N_{\mathrm{lines}}N_{\nu}\), where the individual factors denote the number of cells, the number of rays sent from each cell, the number of lines, and the number of frequency bins (the extra factor of \(N_{\mathrm{cells}}^{1/3}\) accounts for the average number of cells an individual ray must traverse). In the simulation of Greif (2014), these parameters were chosen to achieve an overall accuracy of 5%. For \(N_{\mathrm{cells}}\simeq2\times10^{7}\), \(N_{\mathrm{rays}}=48\), \(N_{\mathrm{lines}}=32\), and \(N_{\nu}=8\), nearly 1014 opacity calculations per time step were carried out. A single ray-tracing step therefore took about \(2\text{,}000\mbox { s}\) on 1,024 modern computing cores. To evolve the gas from \(n_{\mathrm{H}}\simeq10^{10}\mbox { cm}^{-3}\) to \(n_{\mathrm{H}}\simeq10^{15}\mbox { cm}^{-3}\), where the gas becomes optically thick to H2 line emission, about 2,000 ray-tracing steps had to be performed. This corresponds to a total runtime of about one month and a computational cost of nearly one million CPU hours. The simulation was run using a hybrid MPI/Open-MP parallelization scheme on the Sandy-bridge cores of the supercomputer Stampede at the Texas Advanced Supercomputer (TACC). Among many other optimizations, the opacity calculation exploited the Intel AVX instruction set, which allows up to eight single-precision operations to be performed simultaneously.

Another challenge of the radiative transfer calculation is the amount of memory required to store the energy associated with the individual rays. In single precision, \(4N_{\mathrm{cell}}N_{\mathrm{rays}}N_{\mathrm{lines}}N_{\nu}\) bytes must be reserved, which corresponds to \({\simeq}1\mbox { GB}\) per core for 1,024 cores. Each global communication step therefore requires nearly the entire memory to be exchanged between MPI tasks, which accounts for approximately 50% of the total computational cost. Another problem is the significant imbalance in the number of rays stored on the tasks. Since the domain decomposition is optimized for hydrodynamical simulations instead of radiative transfer simulations, a significant imbalance accumulates as the rays are traversed. Typically, this amounts to another factor of two, and leads to a similar reduction in performance. It is therefore not yet beneficial to tailor the scheme for graphics processing units (GPU’s) or co-processors.

An alternative method was recently introduced by Hartwig et al. (2014b). This study used the treecol algorithm of Clark et al. (2012) to determine the total column density of H2 molecules in various directions around each cell. The method accounts for relative velocities in a simplified manner, and computes the average photon escape fraction from the optical depth of the individual transitions. The computational cost is only about five times higher than in a simulation without radiative transfer, and allowed Hartwig et al. (2014b) to follow the build-up and fragmentation of the disk around the primary protostar. During the initial collapse, the escape fraction agrees relatively well with that of Greif (2014), while in the later stages of the collapse the disk allows the radiation to escape more easily than previous approximate methods would imply. This reduces the thermal stability of the gas, and promotes fragmentation.

2.7 Collision-induced emission and absorption

For densities \(n_{\mathrm{H}}\gtrsim10^{14}\mbox { cm}^{-3}\), collision-induced emission and absorption become important (Omukai and Nishi 1998; Ripamonti et al. 2002; Ripamonti and Abel 2004). These processes are characterized by two hydrogen molecules approaching each other and forming a ‘super-molecule’ with a dipole moment induced by van der Waals forces. The resulting super-molecule can emit or absorb radiation via dipole transitions, and has higher transition probabilities than the quadrupole transitions of isolated H2 molecules. As opposed to H2 line emission, where the line width is small compared to the separation of the individual lines, the short lifetime of the super-molecule results in line widths large enough that they merge into a continuum. Yoshida et al. (2006) found that collision-induced emission continues to cool the gas even after it has become optically thick to H2 line emission. The various stages in the collapse of the gas up to this point are shown in Figure 2.

Figure 2
figure 2

Top panel : thermal evolution of primordial gas as it collapses over many orders of magnitude in density. The labels denotes various milestones in the collapse: heating to the virial temperature of the halo (A), runaway cooling via H2 line emission (B), cooling to the minimum temperature of \({\simeq}200\mbox { K}\) (C), onset of three-body H2 formation (D), gas becomes optically thick to H2 line emission (E), onset of collision-induced emission (F), and collisional dissociation of H2 (G). Bottom panel: H2 fraction versus density. The H2 fraction increases from its cosmological abundance of 10−6 to 10−3 via associative detachment of H and H. Following an extended plateau where the H2 fraction remains nearly constant, the cloud becomes fully molecular once three-body reactions set in. Adapted from Yoshida et al. (2006).

In the one-dimensional dynamical model of Ripamonti and Abel (2004), the gas becomes optically thick to collision-induced emission at densities \(n_{\mathrm{H}}\simeq10^{16}\mbox { cm}^{-3}\). This rapid transition compared to H2 line emission is due to the absence of absorption wings, which increases the optical depth of the gas. Ripamonti and Abel (2004) also found that a chemothermal instability similar in nature to the one induced by three-body H2 formation may be triggered. However, the growth rate is longer than the free-fall time, such that it most likely does not result in sub-fragmentation of the cloud (see also Yoshida et al. 2006). In three-dimensional simulations, the continuum opacity of the gas has been taken into account by using the escape fraction method. The escape fraction is computed by evaluating the optical depth along the principal axes of the computational domain (Yoshida et al. 2008; Hirano and Yoshida 2013), or by using a density-dependent fitting function (e.g. Clark et al. 2011a; Greif et al. 2011a). The results achieved with these treatments can differ substantially for asymmetric cloud configurations, but agree relatively well with one-dimensional calculations for clouds that are approximately spherically symmetric (Hirano and Yoshida 2013).

2.8 Collisional dissociation

The last process that is able to cool the gas is collisional dissociation of H2. It acts as a thermostat by converting the compressional energy of the gas into energy that is used for the dissociation of H2, which has a binding energy of 4.48 eV. The temperature and density at which H2 begins to dissociate depends on the ratio of the three-body formation rate to the collisional dissociation rate. In Turk et al. (2009), this occurs already at temperatures \(\gtrsim2\text{,}000\mbox { K}\) at a density of \(n_{\mathrm{H}}\simeq 10^{15}\), while other studies find this to occur at significantly higher densities (Yoshida et al. 2008; Clark et al. 2011b; Greif et al. 2012). This discrepancy is likely related to differences in the reaction rates and the treatment of the optically thick H2 line cooling rate (Turk et al. 2011; Greif 2014).

From a computational standpoint, obtaining a self-consistent H2 fraction at these densities is computationally expensive, since the three-body formation rate scales with the cube of the density. However, the H2 abundance approaches chemical equilibrium at \(n_{\mathrm{H}}\gtrsim10^{15}\mbox { cm}^{-3}\), and the electron abundance at \(n_{\mathrm{H}}\gtrsim10^{19}\mbox { cm}^{-3}\). Once the density exceeds these critical values, an equilibrium solver instead of a non-equilibrium solver may be used, which greatly reduces the computational effort. Instead of evolving a system of stiff differential equations on a time scale that may be much smaller than the Courant time, the chemistry simplifies to a root-finding problem (Omukai and Nishi 1998; Ripamonti et al. 2002; Yoshida et al. 2006, 2008; Greif et al. 2012).

3 First stars: the accretion phase

3.1 Protostar formation and evolution

Following the dissociation of H2, the temperature rises steeply with increasing density. At densities \(n_{\mathrm{H}}\simeq 10^{20}\mbox { cm}^{-3}\), the collapse stalls and an accretion shock forms that heats the gas to \({\gtrsim}10^{4}\mbox { K}\) (Omukai and Nishi 1998; Yoshida et al. 2008). This marks the formation of a protostar at the center of the cloud with an initial mass of \({\simeq}0.01\mbox { M}_{\odot}\) and a size of \({\simeq}0.1\mbox { au}\).

The first calculations of the evolution of primordial protostars with continuous accretion were carried out by Stahler et al. (1986a, 1986b) and Omukai and Palla (2001, 2003). They solved the stellar structure equations for a hydrostatic core of radius \(R_{*}\) bounded by an accretion shock with a constant accretion rate of about \(5\times10^{-3}\mbox { M}_{\odot} \mbox { yr}^{-1}\). A radiative precursor forms ahead of the accretion shock due to the high opacity of the gas, which may be many times larger than the protostar. Up to a mass of about 10 M, the accretion time \(t_{\mathrm{acc}}=M_{*}/\dot{M}_{*}\), where \(M_{*}\) is the mass of the protostar, is smaller than the Kelvin-Helmholtz time \(t_{\mathrm{KH}}=GM_{*}^{2}/ (R_{*}L_{*} )\), where \(L_{*}\) is the luminosity of the protostar. As a result, the protostar evolves nearly adiabatically and the radius increases to \({\gtrsim}100\mbox { R}_{\odot}\). Once \(t_{\mathrm{acc}}>t_{\mathrm{KH}}\), this trend is reversed and the protostar undergoes Kelvin-Helmholtz contraction, which continues until \(M_{*}\simeq50\mbox { M}_{\odot}\). Following a phase of deuterium burning, hydrogen burning begins after \({\simeq}10^{4}\mbox { yr}\), and the protostar settles onto the main sequence after \({\simeq}10^{5}\mbox { yr}\). Accretion is finally terminated when the luminosity becomes comparable to the Eddington luminosity, and the star has grown to a few hundred solar masses. These results agree with the more recent calculations of Hosokawa and Omukai (2009) and Ohkubo et al. (2009).

One of the major uncertainties of these models is the assumption of spherical symmetry. This was relaxed in Tan and McKee (2004), where a semi-analytic prescription for an accretion disk was employed. They found that the reduced density in the polar regions around the protostar results in a reduced optical depth, such that the photosphere of the spherically symmetric models disappears. This may have a substantial effect on the propagation of radiation from the protostar, and the maximum mass that can be reached.

3.2 Accretion

Omukai and Palla (2003) showed that one of the main parameters that governs the evolution of protostars is the accretion rate. Neglecting radiative transfer effects and assuming that the gas does not fragment, the time-averaged accretion rate may be estimated as \(\dot{M}\sim M_{\mathrm{J}}/t_{\mathrm{ff}}\sim c_{\mathrm{s}}^{3}/G\propto T^{3/2}\), where the Jeans mass is evaluated at the density and temperature when the gas first becomes gravitationally unstable, i.e. \(n_{\mathrm{H}}\simeq10^{4}\mbox { cm}^{-3}\) and \(T\simeq200\mbox { K}\). Since the temperature in present-day star-forming regions is \({\simeq}10\mbox { K}\), the accretion rate in primordial gas clouds is expected to be nearly 100 times higher, with \(\dot{M}\sim10^{-3}\mbox { M}_{\odot} \mbox { yr}^{-1}\). Due to the Courant constraint, it is not possible to model the accretion phase accurately for more than a few free-fall times, or \({\simeq}10\mbox { yr}\) (Ripamonti et al. 2002; Greif et al. 2012). Since primordial stars are expected to accrete for more than \({\sim}10^{5}\mbox { yr}\) before radiation feedback terminates accretion (Tan and McKee 2004; McKee and Tan 2008), various methods have been used to estimate the final masses of Population III stars. In Omukai and Nishi (1998), the self-similar form of the spherically symmetric infall solution for \(\gamma_{\mathrm{eff}}\simeq1.1\) was extended to the accretion phase, resulting in \(\dot{M}=8.3\times10^{-2}\mbox { M}_{\odot} \mbox { yr}^{-1}\ (t/1\mbox { yr} )^{-0.27}\), where t denotes the time (Larson 1969; Penston 1969; Shu 1977; Yahil 1983; Suto and Silk 1988). The mass of the star after 105 yr is thus \(M_{*}\simeq 500\mbox { M}_{\odot}\). A similar power-law for the accretion rate was found in the calculations of Ripamonti et al. (2002) and Tan and McKee (2004).

In three-dimensional simulations, the accretion rate may be obtained by measuring the spherically averaged density and velocity profiles (Abel et al. 2002; Yoshida et al. 2006; Gao et al. 2007; O’Shea and Norman 2007; Turk et al. 2011). These instantaneous accretion rates generally agree with the power laws found in one-dimensional calculations. Studies that employed ‘sink particles’ to represent growing protostars determined the accretion rate from the mass accreted by the sink particles (Bromm and Loeb 2004; Clark et al. 2008; Stacy et al. 2010, 2011b, 2012a; Stacy and Bromm 2013, 2014; Clark et al. 2011a, 2011b; Greif et al. 2011a; Smith et al. 2011, 2012a). At sufficiently late times, they agree with those found in previous studies, but show substantial variability due to the instability and fragmentation of the accretion disk around the central protostar.

3.3 Protostellar radiation

Stellar radiation may shut down the accretion flow from the surrounding gas reservoir before the star reaches masses well in excess of 100 M. Since the influence of both stellar winds and radiation pressure are expected to be small in the absence of metals and dust grains (Baraffe et al. 2001; Bromm et al. 2001b; Marigo et al. 2001, 2003; Kudritzki 2002; Omukai and Inutsuka 2002; Krtička et al. 2003, 2010; Krtička and Kubát 2006, 2009; McKee and Tan 2008; Sonoi and Shibahashi 2011, 2012; Muijres et al. 2012; Sonoi and Umeda 2012), the most important feedback mechanisms are the dissociation of H2 and the heating that accompanies the ionization of neutral hydrogen atoms. The former occurs due to radiation in the so-called Lyman-Werner (LW) bands, which indirectly dissociates H2 and removes the most important coolant of the gas (Omukai and Nishi 1999; Glover and Brand 2001). Hosokawa et al. (2011) argued that this process is not particularly important, since the disk remains optically thick to LW radiation, while the molecules in the bipolar regions perpendicular to the disk are collisionally dissociated by photoheated gas (see Figure 3).

Figure 3
figure 3

Two-dimensional simulation of primordial star formation including radiation from the central protostar. The panels shows the temperature and number density of the cloud at various output times. The ionizing radiation clears out the gas along the poles, while the accretion continues nearly unabated in the plane of the disk. As the photo-heated region expands, it shuts down the mass flow from the envelope onto the disk, which eventually terminates accretion. Here, the star reaches a final mass of \({\simeq}40\mbox { M}_{\odot}\). Adapted from Hosokawa et al. (2011).

Once the star undergoes Kelvin-Helmholtz contraction, the effective temperature becomes high enough for significant amounts of ionizing radiation to be produced. The energy above \({\simeq}13.6\mbox { eV}\) thermalizes and heats the gas to \({\gtrsim}10^{4}\mbox { K}\), which is significantly higher than the virial temperature of a minihalo. The high pressure therefore begins to drive gas from the halo. However, in the one-dimensional calculations of Omukai and Inutsuka (2002), the nascent H ii region remains compact due to the high inflow velocity of the gas, and does not impede the accretion flow. The accretion disk model of McKee and Tan (2008), on the other hand, predicts a reduced density and inflow velocity along the poles, which allows the ionizing radiation to escape and begin to photo-evaporate the disk. This is confirmed by the two-dimensional simulations of Hosokawa et al. (2011) and the three-dimensional simulations of Stacy et al. (2012a), which used ray-tracing methods to compute the propagation of the radiation from the protostars. In the case of Hosokawa et al. (2011), flux-limited diffusion was used to model the diffuse component of the radiation, and the effects of radiation pressure were included. These studies found an upper mass limit of a few tens of solar masses.

Recently, Hirano et al. (2014) used the methodology of Hosokawa et al. (2011) to investigate the influence of radiation in a sample of 100 different minihalos that were extracted from three-dimensional cosmological simulations. They found final stellar masses in the range \(\simeq10\mbox{-}1\text{,}000\mbox { M}_{\odot}\), which are correlated with the thermal evolution of the gas during the initial collapse phase (see Figure 4). In contrast to these studies, Susa (2013) and Susa et al. (2014) found that the molecule-dissociating radiation from the central protostar is more important and indirectly shuts down accretion. However, their resolution was not high enough to resolve the ionization front along the poles. They found that most Population III stars had final masses in the range \(10\lesssim M_{*}\lesssim100\mbox { M}_{\odot}\).

Figure 4
figure 4

IMF of Population III according to the two-dimensional radiation hydrodynamics simulations of Hosokawa et al. ( 2011 ), but applied to 100 different minihalos. The various colors and hatchings represent different paths of protostellar evolution: P1 denotes a Kelvin-Helmholtz contracting protostar, P1hd a halo in which HD cooling was activated, P2 an oscillating protostar, and P3 and P3p super-giant protostars. Adapted from Hirano et al. (2014).

Despite the significant progress made by these studies, some caveats remain. The simulations of Hosokawa et al. (2011) and Hirano et al. (2014) included the most detailed treatment of radiative transfer, but were limited to two dimensions. The three-dimensional simulations of Stacy et al. (2012a), Susa (2013) and Susa et al. (2014), on the other hand, suffered from limited resolution and a simplified treatment of the radiative transfer. Significant progress could be made if the physical detail of the simulations of Hosokawa et al. (2011) were applied to three dimensions.

3.4 Fragmentation

Fragmentation during the accretion phase may further reduce the mass of the central protostar. If the cloud fragments shortly after the formation of the first protostar, much of the material in the surrounding envelope may be accreted by secondary protostars before they reach the center of the cloud. This process has been termed fragmentation-induced starvation (Peters et al. 2010). Although fragmentation may occur already during the initial collapse phase (Turk et al. 2009; Greif et al. 2013), most of the fragmentation is expected to occur after the first protostar has formed. Due to the Courant constraint, numerical simulations are not able to self-consistently evolve the collapse of the gas for a significant period of time beyond the initial collapse (Omukai and Nishi 1998; Ripamonti et al. 2002; Greif et al. 2012). For this reason, sink particles have been used to represent growing protostars and circumvent the need to model their interior structure and evolution (Bate et al. 1995; Krumholz et al. 2004; Jappsen et al. 2005; Federrath et al. 2010).

Next to criteria involving the thermal, kinetic, and gravitational binding energy of the gas, as well as the properties of the local velocity field, sink particles are typically inserted at a threshold density and accrete gas within a predefined accretion radius. This prevents the gas from collapsing to increasingly high densities, and allows the simulations to be continued for a much longer period of time than would otherwise be possible. By nature, the sink particle method has its own limitations. The boundary conditions imposed on the gas near the accretion radius are necessarily artificial and may lead to unphysical results. The complicated torques on the scale of the accretion radius are not captured accurately, which might artificially enhance or prevent fragmentation. The loss of resolution on the scale of the accretion radius also sets a minimum scale on which fragmentation can occur. Close encounters between sink particles may result in dynamical ejections, although the orbital energy may in fact be dissipated through unresolved torques. Nevertheless, the sink-particle technique is currently the only method to investigate the evolution of primordial gas clouds over significant periods of time, and efforts are underway to increase their physical realism (Hubber et al. 2013).

Some of the first studies that employed sink particles inserted them on scales significantly larger than the expected sizes of the protostars (Bromm et al. 2002; Bromm and Loeb 2004; Stacy et al. 2010; Stacy and Bromm 2013). It was therefore not clear whether the gas represented by the sink particles would further sub-fragment. Clark et al. (2008) addressed this problem by using a threshold density of \(n_{\mathrm{H}}=10^{16}\mbox { cm}^{-3}\), which corresponds to a radius of \({\simeq}1\mbox { au}\) and is the approximate size of a metal-free protostar that has a realistic accretion rate (Stahler et al. 1986a, 1986b). The initial conditions were designed to reproduce the end state of the primordial gas clouds found in Bromm et al. (2002), and they used the tabulated equations-of-state of Omukai et al. (2005) for various metallicities instead of a chemical network. In the primordial case, the central gas cloud fragmented into 25 protostars after only a few hundred years, with masses ranging from 0.02 to \({\simeq}10\mbox { M}_{\odot}\). As opposed to the simulations with a non-negligible metallicity, the mass function in the primordial case was relatively flat, implying that most of the mass is locked up in high-mass protostars. In a subsequent study, Clark et al. (2011a) improved upon their previous work by using a detailed chemistry and cooling network. They investigated how the fragmentation depends on the mean Mach number of the turbulence, and found that a higher degree of turbulence generally leads to more fragmentation.

The remaining caveat of idealized initial conditions was addressed in Clark et al. (2011b). They employed cosmological initial conditions and a hierarchical zoom-in procedure to extract the central, Jeans-unstable cloud (Navarro and White 1994; Tormen et al. 1997; Gao et al. 2005). The resolution was increased several times with a particle-splitting technique (Kitsionas and Whitworth 2002; Bromm and Loeb 2003a). This allowed the density at which sink particles were created to be increased to \(n_{\mathrm{H}}=10^{17}\mbox { cm}^{-3}\). Clark et al. (2011b) also included the effects of radiative heating due to accretion using a Planck mean opacity, based on the assumption that the gas was optically thin. They found that an accretion disk forms around the central protostar, which develops pronounced spiral arms. The first fragments appear after about 100 yr, and by the end of the simulation three additional fragments have formed. The fragmentation is attributed to the limited rate at which mass can be processed through the disk, which is significantly smaller than the accretion rate from the surrounding envelope onto the disk. The surface density of the disk increases, while the compressional energy is radiated away by H2 line emission. This quickly drives the disk towards gravitational instability (see Figure 5). Previous semi-analytic studies found that the disk remains stable, since they assumed that H2 cooling is inefficient at the densities and temperature at which the disk forms (Tan and McKee 2004; Tan and Blackman 2004; Mayer and Duschl 2005). Tanaka and Omukai (2014) found that primordial disks remain marginally stable despite the H2 cooling (but see Latif and Schleicher 2014a, 2014b). It is possible that the ability of the gas to cool may have been overestimated in previous studies that employed the Sobolev method to compute the optically thick H2 line cooling rate (see Greif 2014). However, the more accurate method used in Hartwig et al. (2014b) showed that this does not substantially affect the fragmentation of the gas.

Figure 5
figure 5

From top left to bottom right : Surface density, temperature, \(\pmb{\mbox{H}_{2}}\) fraction, and Toomre Q parameter in the Keplerian disk that forms around the central protostar. For \(Q\lesssim1\), the disk becomes gravitationally unstable and fragments. In this simulation, the first fragment forms after \({\simeq}90\mbox { yr}\) at a distance of \({\simeq}20\mbox { au}\) from the central protostar. Adapted from Clark et al. (2011b).

The simulations of Greif et al. (2011a) were similar in nature to those of Clark et al. (2011a), but employed a moving-mesh (Springel 2010) instead of an SPH approach (Springel et al. 2001; Springel 2005). Using the same chemical model and a similar sink-particle scheme, Greif et al. (2011a) investigated fragmentation in five different minihalos over a ten times longer period of time than in Clark et al. (2011a). They found that on average ten protostars per halo had formed after \({\simeq}10^{3}\mbox { yr}\), with masses ranging from 0.1 to \({\simeq}10\mbox { M}_{\odot}\). The mass function was similarly flat as in Clark et al. (2008). Greif et al. (2011a) also found that a number of protostars with masses below 0.8 M formed, which would allow them to survive to the present day. These protostars were ejected from the central gas cloud by N-body interactions with more massive protostars and stopped accreting. In some cases, their velocities exceeded the escape velocity of the halo. However, the fraction of ejected protostars also depended strongly on the implementation of the sink-particle method, since the gas surrounding the protostars may absorb a significant fraction of their orbital energy.

The effects of radiative heating by the protostars was investigated by Smith et al. (2011) and Smith et al. (2012a), using the halos of Greif et al. (2011a) and the radiative transfer method employed in Clark et al. (2011b). Due to the low effective temperature of the protostars during the adiabatic expansion phase, most of the radiation is in the infrared and therefore serves to only slightly heat the gas. Fragmentation occurs somewhat later than in Greif et al. (2011a), and the number of fragments that form is somewhat reduced.

3.5 Protostellar system

In an attempt to avoid the uncertainties introduced by sink particles, Greif et al. (2012) self-consistently evolved the collapse of the gas beyond the formation of the first protostar. Due to the substantial computational effort involved in resolving the Jeans length with 32 cells at all times, only the first \({\simeq}10\mbox { yr}\) in the evolution of the protostellar disk could be modeled. Five different halos with the same initial conditions as in Greif et al. (2011a) were investigated. In analogy to previous studies, a Keplerian disk formed that developed spiral arms and fragmented. The gravitational stability of the disk was evaluated using the so-called Toomre (1964) and Gammie (2001) criteria. The Toomre Q parameter, which quantifies the stability of the disk to perturbations, is given by \(Q=c_{\mathrm{s}}\kappa/ (\pi G\Sigma )\). Here, κ is the epicyclic frequency of the perturbation, which is equal to the angular velocity Ω in a Keplerian disk, and Σ is the surface density. For \(Q<1\), perturbations in the disk can grow, but these perturbations do not necessarily result in fragmentation. The latter requires the Gammie criterion, \(t_{\mathrm{cool}}<3/\Omega\), where \(t_{\mathrm{cool}}\) is the cooling time, to be satisfied. Greif et al. (2012) showed that both of these conditions are fulfilled on a scale of \({\simeq}1\mbox { au}\), which agrees with the location of the fragmentation in the simulations. The efficient cooling of the gas necessary to fulfill Gammie’s criterion stems mainly from the dissociation of H2, which acts as a thermostat by extracting compressional energy from the gas and using it to unbind H2. In this respect the results of Greif et al. (2012) differ somewhat from Clark et al. (2011a), where H2 line emission was found to be the most important coolant. This is most likely due to the higher resolution employed in Greif et al. (2012), which allowed the collapse of the gas to be followed to significantly higher densities, where collisional dissociation cooling begins to dominate over H2 line emission.

The dependence of the stability of the disk on the abundance of H2 implies that differences in the three-body formation rates and the treatment of the optically thick H2 line cooling rate may significantly affect the fragmentation of the gas (e.g. Turk et al. 2012). The recent simulations of Greif (2014) show that the ability of the gas to cool may indeed have been overestimated. However, the strong asymmetry of the cloud that develops over time should allow the radiation to escape more easily than in the spherically symmetric model for the optically thick H2 line cooling rate used in Turk et al. (2012). Indeed, this trend has been confirmed by the simulations of Hartwig et al. (2014b). In the future, it will be useful to evolve the self-consistent simulations of Greif (2014) to higher densities, but with modifications to the radiative transfer scheme that make them computationally feasible.

The explicit resolution of the interface between the protostars and the disk in Greif et al. (2012) allowed the complex interactions of the protostars with other protostars and the surrounding gas to be modeled self-consistently (see Figure 6). This study found that the protostars are subject to strong gravitational torques that lead to the migration of about half of the secondary protostars formed in the disk to the center of the cloud, where they merge with the primary protostar. The aggressive migration and merging of the protostars occurs on a free-fall time scale. In analogy to Greif et al. (2011a), some low-mass protostars migrate to higher orbits due to N-body interactions, but do not become nearly as unbound as in Greif et al. (2011a). This is mainly because the protostars have extremely large radii of the order of 100 R, such that protostellar interactions do not resemble those of the point masses in Greif et al. (2011a). Much of the orbital energy is transferred to the surrounding gas instead of the protostars. At the end of the simulation, about five protostars per halo are present, with a tendency to further increase. The mass budget is dominated by the primary protostar, which grows to about five times the mass of all other protostars.

Figure 6
figure 6

From left to right : the first \(\pmb{{\simeq}10\mbox { yr}}\) of the evolution of the protostellar system that forms at the center of four different minihalos. In all cases, a disk forms that becomes gravitationally unstable and fragments into a small system of protostars. Complicated gravitational torques result in a fraction of the protostars migrating to the center of the cloud, where they merge with the primary protostar, while others migrate to higher orbits. Adapted from Greif et al. (2012).

In a qualitatively similar study, Latif et al. (2013c) artificially truncated the collapse of the gas at a density of \(n_{\mathrm{H}}\simeq 10^{13}\mbox { cm}^{-3}\). Even though the size of the disk was much larger than in Greif et al. (2012), they found a similar susceptibility to fragmentation. Vorobyov et al. (2013) employed two-dimensional simulations with a predefined effective equation-of-state to investigate the evolution of a metal-free protostellar system for \({\simeq}5\times10^{4}\mbox { yr}\). In this study, the secondary protostars quickly merged with the primary protostar, but the disk continued producing new fragments, with of order 10 protostars present at any given time.

Despite the limitations of these simulations, they demonstrate the advantages and disadvantages of the sink-particle technique. Sink particles allow the simulations to be continued for much longer than would otherwise be possible, and probe the chemical and thermal evolution of the gas on scales larger than the accretion radius. On the other hand, they are not well suited to predict the properties of the protostars themselves, which depend sensitively on the properties of the gas on scales comparable to or smaller than the accretion radius. Since self-consistent simulations are extremely expensive, sink particles remain the only method capable of probing sufficiently late times in the Population III star formation process.

3.6 Stellar rotation

Early simulations showed that the cloud out of which the first protostellar seed forms tends to also have a strong rotational component (e.g. Abel et al. 2002; Bromm et al. 2002). Due to turbulence and other instabilities, the angular momentum is constantly redistributed as the gas continues to collapse, maintaining rotational velocities that are a significant fraction of the Keplerian velocity (e.g. Yoshida 2006; Greif et al. 2012). The protostar that forms at the center of the cloud is therefore endowed with significant rotation. Since it accretes material from a Keplerian disk, the angular momentum of the protostar may remain high as it evolves towards the main sequence. In simulations of present-day star formation, magnetic braking may reduce the rotation rate of the protostar, but their influence in the primordial case is not yet clear. Rapidly rotating Population III stars can have very different properties than their non-rotating counterparts, since rotation affects their evolution on the main sequence and beyond, as well as the degree to which the nucleosynthetic products mix with each other. The amount of rotation may also affect the types and properties of their SN explosions (see Section 5.7).

Stacy et al. (2011a) employed sink particles that traced the angular momentum of the accreted gas to investigate the rotation rate of primordial protostars. They found that the protostar represented by the sink particle rotated at near break-up speeds throughout the \({\simeq}5\text{,}000\mbox { yr}\) that they simulated. In a follow-up study, Stacy et al. (2013) analyzed the rotation of the protostars formed in the simulations of Greif et al. (2012), which self-consistently modeled the interface between the protostars and the surrounding gas. Similar to Stacy et al. (2011a), they found that the central protostar rotated at a significant fraction of the Keplerian velocity, despite strong gravitational torques due to frequent mergers with secondary protostars. Stellar evolution models should therefore account for the rotation of the star.

4 First stars: additional physics

4.1 HD cooling

The second low-temperature coolant in primordial gas that may become important is hydrogen deuteride (HD; Lepp and Shull 1984; Puy et al. 1993; Stancil et al. 1998; Galli and Palla 1998, 2002; Flower and Roueff 1999; Flower et al. 2000; Flower 2000; Flower and Pineau des Forêts 2001; Lipovka et al. 2005; Johnson and Bromm 2006; Glover and Abel 2008). As opposed to H2, HD possesses a permanent dipole moment, with correspondingly higher radiative transition probabilities. This shifts the critical density at which the level populations transition to LTE to \(n_{ \mathrm{H}, \mathrm{crit}}\sim10^{6}\mbox { cm}^{-3}\). Furthermore, since ortho and para states do not exist, the transition \(J=1\rightarrow0\) is accessible, corresponding to an excitation temperature of \({\simeq}130\mbox { K}\). The primary formation reaction is:

$$ \mathrm{H}_{2} +\mathrm{D}^{+}\rightarrow\mathrm{HD}+ \mathrm{H}^{+}, $$
(8)

showing that the HD formation rate depends on the ionization state of the gas and the abundance of H2. Since HD is primarily destroyed by the inverse reaction of equation (8), the difference in the zero-point energies of H2 and HD can lead to a significant degree of chemical fractionation at low temperatures, which boosts the HD/H2 ratio well above the cosmological D/H ratio of about 10−5.

In minihalos, the temperature usually does not become low enough for HD cooling to become important (Bromm et al. 2002). However, more recent studies have found that HD cooling may in fact dominate over H2 cooling in certain halos. Ripamonti (2007) and McGreer and Bryan (2008) showed that HD cooling is activated in high-redshift minihalos with masses below \({\simeq}3\times10^{5}\mbox { M}_{\odot}\). This allows the gas to cool down to the temperature of the CMB, which lies between 50 and \({\simeq}100\mbox { K}\) at \(z\gtrsim20\). The corresponding Jeans mass is up to an order of magnitude lower than in the pure H2 cooling case, and may lead to the formation of a distinct population of metal-free stars with a characteristic mass of 10 M (see also Nakamura and Umemura 2002).

A similar decrease of the initial Jeans mass was found in some of the halos investigated in Greif et al. (2011a) and Greif et al. (2013). In the radiation hydrodynamics simulations of Hosokawa et al. (2012b) and Hirano et al. (2014), the lower accretion rates in halos in which HD cooling was activated led to final stellar masses \({\gtrsim}10\mbox { M}_{\odot}\). Clark et al. (2011a) came to a somewhat different conclusion. They found that the gas in these clouds heated more rapidly at densities \(n_{\mathrm{H}}\gtrsim n_{ \mathrm{H}, \mathrm{crit}}\), which suppresses turbulence and results in less fragmentation. The accreted mass is therefore distributed to fewer protostars, which allows them to become more massive. The influence of fragmentation on the growth of the protostars could not be addressed in the two-dimensional simulations of Hosokawa et al. (2012b) and Hirano et al. (2014), while Clark et al. (2011a) evolved the central cloud for only a few thousand years, and did not model the radiation from the protostars. It is therefore not yet clear whether HD cooling acts to increase or decrease the typical mass of Population III stars. However, the total mass in stars is likely lower in the HD cooling case. The number of minihalos in which HD cooling is activated may be significantly reduced by a global LW background, which dissociates H2 and thus reduces the rate at which HD forms (Wolcott-Green and Haiman 2011).

If the electron abundance is significantly enhanced with respect to the post-recombination value, HD cooling may become important in other circumstances as well. An elevated electron abundance is produced by the virialization shocks of atomic cooling halos or during mergers (Greif and Bromm 2006; Greif et al. 2008; Shchekinov and Vasiliev 2006; Vasiliev and Shchekinov 2008; Prieto et al. 2012, 2014; Bovino et al. 2014b), as well as in SN remnants (Mac Low and Shull 1986; Uehara and Inutsuka 2000; Mackey et al. 2003; Machida et al. 2005; Vasiliev and Shchekinov 2005a; Vasiliev et al. 2008). The H2 abundance in the post-shock gas increases to values well above those found in minihalos, which allows the gas to cool to temperatures low enough that chemical fractionation occurs and HD cooling takes over. A similar process occurs in relic H ii regions, where the recombination time is long enough to facilitate the formation of HD (Nagakura and Omukai 2005; Johnson and Bromm 2007; Yoshida et al. 2007a, 2007b; Machida et al. 2009a). A distinct population of metal-free stars may thus arise under circumstances where the gas has been affected by radiation, but is not yet enriched with metals.

4.2 Magnetic fields

Magnetic fields have been neglected in most studies of primordial star formation, even though it has been shown that they are important in protogalactic halos (Pudritz and Silk 1989; Beck et al. 1994; Lesch and Chiba 1995; Kulsrud et al. 1997). The magnetic field strength at the time of the formation of the first stars is constrained to \(B\lesssim1\mbox { nG}\) (Schleicher et al. 2008), although they are likely much weaker. Potential seed fields are generated during inflation, the electroweak phase transition, and the quark-hadron phase transition (for a recent review, see Widrow et al. 2012). A more robust formation mechanism is the Biermann battery (Biermann 1950), which requires the density gradient in an ionized gas to be misaligned with the pressure gradient, and is thus coupled to the vorticity of the gas. Xu et al. (2008) investigated the growth of the magnetic field in a three-dimensional simulation of a minihalo via the Biermann battery, and found that it generates a seed field of the order of 10−18 G (see also Doi and Susa 2011). Another possible formation mechanism at a later stage in the collapse of the gas is via radiative forces (e.g. Shiromoto et al. 2014).

If the time scale for ambipolar or Ohmic diffusion is large compared to the evolutionary time of a system, the magnetic field is ‘frozen’ into the gas and moves with it (but see Maki and Susa 2004, 2007). In a spherically symmetric, contracting gas cloud the magnetic field grows \(\propto\rho^{2/3}\), while more asymmetric configurations lead to a shallower power-law. For a gas in which the temperature does not evolve substantially, the ratio of the thermal pressure to the magnetic pressure thus decreases \(\propto\rho^{-1/3}\), such that a seed field generated by the Biermann battery and amplified by flux-frozen collapse does not become dynamically important. A much more potent amplification mechanism is the turbulent dynamo (Parker 1963; Kraichnan and Nagarajan 1967). Small-scale turbulent motions in the gas repeatedly fold the magnetic field, which can grow by many e-foldings in a free-fall time. Simple calculations have shown that this mechanism can amplify the magnetic field in primordial gas clouds to appreciable levels (Schleicher et al. 2010a; Schober et al. 2012). This was confirmed by three-dimensional simulations that included the effects of magnetic fields (Sur et al. 2010; Peters et al. 2012; Latif et al. 2013d, 2014e). These studies also found that the Jeans length must be resolved by at least 32 cells for significant amplification to occur (see also Federrath et al. 2011), and that the amplification rate does not converge with increasing resolution. This is expected, since saturation requires the viscous scale to be resolved, while the Reynolds number in primordial gas clouds is of order 105. Using cosmological initial conditions and a detailed chemical model, Turk et al. (2012) confirmed the results of Sur et al. (2010), but found that 64 instead of 32 cells per Jeans length are necessary (see Figure 7).

Figure 7
figure 7

Volume-averaged magnetic energy scaled by \(\pmb{\rho^{4/3}}\) as a function of density in a simulation of primordial star formation that includes magnetic fields. The various line styles denote the resolution of the Jeans length. If the magnetic flux is frozen into the flow, the magnetic energy evolves along the lower black dashed line. Depending on resolution, the turbulent dynamo amplifies the magnetic field strength well above this value. Here, at least 64 cells per Jeans length are required for significant amplification to occur. The shaded regions show the distribution of the gas in the \(J=64\) case. There is no indication that the magnetic energy converges with increasing resolution. Magnetic fields will likely become dynamically important during the initial collapse, and affect the susceptibility of the gas to fragmentation. Adapted from Turk et al. (2012).

The turbulent dynamo rapidly amplifies the magnetic field to a level where it becomes dynamically important. In an accretion disk, a strong field can trigger the magneto-rotational instability, or lead to the formation of a hydromagnetic jet that removes material along the poles of the disk and aids the transport of angular momentum (Maki and Susa 2004; Tan and Blackman 2004; Silk and Langer 2006). The idealized, three-dimensional simulations of Machida et al. (2006) demonstrated that a jet indeed forms, and is capable of removing a significant amount of mass from the disk. Later simulations showed that strong magnetic fields may also suppress fragmentation (Machida et al. 2008a; Machida and Doi 2013). Peters et al. (2014) included sink particles and used a polytropic equation of state to model the thermal evolution of the gas. They found that dynamically important magnetic fields delay the onset of fragmentation by an order of magnitude compared to the case where no magnetic field is present.

Irrespective of the strength of the seed field, it appears that the turbulent dynamo amplifies the magnetic field rapidly enough that it becomes dynamically important already during the initial collapse phase. Simple one-dimensional calculations have shown that the effects of ambipolar diffusion are negligible at this point (e.g. Schober et al. 2012). The magnetic field is expected to have the largest effect following the formation of the circumstellar disk. Magnetic braking and the launch of a hydromagnetic jet increase the rate of angular momentum transport through the disk and affects its ability to fragment. Exactly how effective these processes are depends on the initial strength of the magnetic field and how well the turbulent dynamo is resolved. In the near future, it may become possible to also include dissipative processes such as ambipolar diffusion in multi-dimensional simulations, which decrease the strength of the magnetic field.

4.3 Cosmic rays

A background of cosmic rays at high redshift is primarily produced by SNe, and may affect the formation of the first stars. Due to their long mean free paths, they can affect the chemical evolution of primordial gas on cosmological scales. Shchekinov and Vasiliev (2004) and Vasiliev and Shchekinov (2006) found that cosmic rays pre-ionize the gas, which in turn facilitates the formation of H2 and HD. This may decrease the minimum halo mass required to cools the gas by an order of magnitude at \(z\sim20\). Next to the enhanced cooling, cosmic rays also directly heat the gas. In the calculations of Stacy and Bromm (2007), the additional cooling provided by the enhanced H2 and HD formation rate dominates over the heating, allowing the gas in a minihalo to cool to the temperature of the CMB. As discussed in Section 4.1, this may change the characteristic mass of the star. Similar results were found in Jasche et al. (2007), where a larger range of parameters was explored. Nakauchi et al. (2014) investigated the combined effects of cosmic rays and LW radiation on atomic cooling halos. Since the cosmic rays enhanced the formation of H2, a larger LW flux was necessary to suppress cooling and star formation prior to the onset of Ly-α cooling.

4.4 Streaming velocities

An important effect that has been neglected in most studies of primordial star formation is the relative velocity between the DM and gas sourced by baryon acoustic oscillations before recombination (Tseliakhovich and Hirata 2010). This relative velocity has also been termed the ‘streaming velocity’. At \(z\gtrsim1\text{,}000\), sound waves propagate through the IGM and affect the DM and gas differently. While the gas pressure decelerates the gas, the DM only interacts gravitationally and maintains its velocity. The magnitude of the relative velocity is related to the effective sound speed of the gas, which is of the order of the speed of light before recombination due to Compton scattering. After recombination, the effective sound speed drops to that of a gas with a temperature of \({\sim}10^{4}\mbox { K}\), such that the streaming velocity becomes supersonic. A 1σ-peak at \(z\simeq1\text{,}000\) has a relative velocity of \(v_{\mathrm{rel}}\simeq30\mbox { km}\mbox { s}^{-1}\) and a Mach number of \(\mathcal{M}\simeq5\). Following recombination, the relative velocity decreases \(\propto a^{-1}\), such that at \(z\gtrsim 20\) the streaming velocity is close to \(1\mbox { km}\mbox { s}^{-1}\). This is comparable to the virial velocity of minihalos, where the effect is expected to be strongest.

The streaming velocity damps density perturbations on the acoustic oscillation scale at recombination, i.e. in the range \(\simeq5\mbox{-}100\mbox { Mpc}\). However, this effect is relatively small and only leads to a 10% suppression of the number density of minihalos (Naoz et al. 2012). The effect on the virialization of the gas in minihalos is more pronounced. Numerical simulations included streaming velocities by either employing a fixed velocity offset between the DM and gas (Greif et al. 2011b; Maio et al. 2011b; Stacy et al. 2011a; Naoz et al. 2013; Richardson et al. 2013; Latif et al. 2014c), or by self-consistently evolving the linear equations in the presence of streaming velocities (O’Leary and McQuinn 2012). They found that streaming velocities reduce the baryon overdensity in minihalos and possibly delay the onset of cooling (see Figure 8). In addition, the center of the gas cloud may be shifted with respect to the center of the DM halo by a separation that is comparable to the virial radius. The moving DM halo induces a bow shock in the gas, which decelerates the halo via dynamical friction. Some of these effects were also found in more simplified semi-analytic calculations (Tseliakhovich and Hirata 2010; Tseliakhovich et al. 2011; Fialkov et al. 2012). The resulting modulation of star formation on large scales may also have a substantial effect on the 21-cm signal (McQuinn and O’Leary 2012; Visbal et al. 2012; Fialkov et al. 2012, 2013, 2014). Tanaka et al. (2013) investigated the influence of streaming velocities on the formation and growth of stellar seed black holes (BHs), finding that they have only a minor effect on the abundance of BHs at late times.

Figure 8
figure 8

Baryonic mass available for star formation in minihalos in a cosmological simulation that includes streaming velocities. The circles, stars, and triangles denote the outcome for Mach numbers 0, 1.9, and 3.8 at \(z=20\), respectively, and the black dashed line shows the gas mass assuming a cosmological baryon to dark matter ratio. As the streaming velocity is increased, the halos retain less gas, the central gas density is reduced, and less gas is able to cool and form stars. The degree to which star formation is suppressed increases with decreasing halo mass. Since the streaming velocities are sourced by acoustic oscillations before recombination, this leads to a modulation of Population III star formation on \({\simeq}10\mbox{-}100\mbox { Mpc}\) scales. Adapted from O’Leary and McQuinn (2012).

4.5 Dark matter annihilation

Despite the fact that the cosmological mass density of DM is well constrained, its nature remains unknown. The most popular model invokes a WIMP, such as the lightest particle predicted by supersymmetry, which has a mass of \({\sim}100\mbox { GeV}\). These are expected to have a very small rate coefficient for self-annihilation of \(\langle \sigma v\rangle \simeq3\times10^{-26}\mbox { cm}^{3}\mbox { s}^{-1}\). However, this may be high enough for DM annihilations to have an effect on the DM and gas. Following a complex chain of reactions, the end products of DM annihilations are electron-positron pairs, neutrinos, and gamma rays, which can ionize and heat the gas. The average DM density is too low for this to have a significant effect on the IGM, but within halos the density-squared dependence of the annihilation rate may render DM annihilations important. In particular, minihalos are expected to have significantly higher central DM densities than halos in the present-day Universe. Their average density is significantly higher, since the virial density scales with the cube of the redshift. They are also more centrally concentrated, since the concentration parameter increases with decreasing halo mass (Navarro et al. 1997). Finally, their formation is expected to be nearly monolithic, since the variance of matter fluctuations is nearly constant towards the low-mass end. Following the collapse of the DM alongside the gas via ‘adiabatic contraction’ (Young 1980; Blumenthal et al. 1986; Gnedin et al. 2004), DM annihilations begin to affect the gas.

Spolyar et al. (2008) and Freese et al. (2009) used a simple dynamical model to show that in the standard WIMP scenario the DM heating rate matches the H2 line cooling rate at a density of \(n_{\mathrm{H}}\simeq10^{13}\mbox { cm}^{-3}\). They conclude that the collapse stalls at this point and a ‘dark star’ powered by DM annihilations forms (see also Natarajan et al. 2009). Stellar structure calculations that included DM-baryon scattering indicate that these stars are much larger than normal Population III stars, have lower surface temperatures, and are more luminous (Freese et al. 2008; Iocco 2008; Iocco et al. 2008a; Taoso et al. 2008; Yoon et al. 2008; Spolyar et al. 2009; Hirano et al. 2011; Sivertsson and Gondolo 2011). Due to their low surface temperatures, they would appear dark at optical wavelengths, and radiative feedback would not be able to impede the accretion flow, allowing them to grow to a final mass of \({\sim}10^{6}\mbox { M}_{\odot}\). With a luminosity of up to 1010 L in the infrared, they may be observable by the James Webb Space Telescope (JWST; Freese et al. 2010; Zackrisson et al. 2010a, 2010b; Ilie et al. 2012), next to other unique signatures (Schleicher et al. 2009; Maurer et al. 2012; Sandick et al. 2012).

One of the main problems of the dark star formation scenario is the assumption that the collapse of the gas stalls once the DM heating rate becomes comparable to the cooling rate. Ripamonti et al. (2010) showed that this is not the case in the standard WIMP scenario. Instead, the temperature increases only marginally until the cooling rate exceeds the DM heating rate, and the collapse continues nearly unhindered. Another problem concerns the inherent assumption of spherical symmetry. If the DM remains aligned with the gas, the annihilation rate is maximized. However, previous studies have shown that primordial gas clouds are permeated by transonic turbulence, develop a disk, and fragment in the later stages of the collapse (e.g. Turk et al. 2009; Clark et al. 2011a; Greif et al. 2012). The influence of these perturbations on the DM profile was investigated by Stacy et al. (2012b). They found that they reduced the annihilation luminosity to a degree where it no longer had a substantial effect on the gas. While Smith et al. (2012b) found that the circumstellar disk was stabilized by annihilation heating, Stacy et al. (2014) attributed this to the use of a spherically symmetric DM density profile. When they allowed the DM potential to vary, they recovered the results of Stacy et al. (2012b).

In summary, while DM annihilation heating may affect the thermodynamic evolution of the gas at moderate densities by ionizing the gas and facilitating the formation of H2 (e.g. Ripamonti et al. 2010), it appears unlikely that the high-density evolution is changed to a degree at which the formation of dark stars instead of ‘normal’ Population III stars is favored.

4.6 Alternative cosmologies

Next to the standard ΛCDM paradigm, a number of alternative cosmologies may be viable. These could alter the formation of the first stars by changing the matter fluctuation power spectrum on small scales. In common gravitino warm dark matter (WDM) models, the particle mass is assumed to be \({\gtrsim}1\mbox { keV}\), since lower values are ruled out observationally (e.g. Narayanan et al. 2000). Yoshida et al. (2003c) investigated the effects of a WDM particle with \(m_{\mathrm{WDM}}=10\mbox { keV}\) on the formation of the first stars. In this case, the matter power spectrum has an exponential cut-off at \({\simeq}0.05\mbox { Mpc}\), which corresponds to a mass of \({\simeq}5\times10^{6}\mbox { M}_{\odot}\). As a result, they found that star formation in minihalos is nearly entirely suppressed, and must await the virialization of larger halos. O’Shea and Norman (2006) used somewhat more detailed simulations to investigate the influence of WDM particles with masses in the range \(10\mbox { keV}\gtrsim m_{\mathrm{WDM}}\gtrsim50\mbox { keV}\). At the lower mass end, they find a similar increase in the halo mass required for gas collapse as Yoshida et al. (2003c). In addition, they found a delay in the onset of runaway cooling. However, once the gas collapsed to a density of \(n_{\mathrm{H}}\simeq10^{5}\mbox { cm}^{-3}\), the thermodynamic evolution of the cloud became nearly indistinguishable from that in the CDM paradigm. Gao and Theuns (2007) used a WDM particle mass of 3 keV and found that the gas collapsed along a filament prior to the onset of runaway cooling (see also Nakamura and Umemura 1999, 2001, 2002; Bessho and Tsuribe 2012, 2013). Since the particle noise exceeds the WDM fluctuation power on small scales, these simulations could not derive the resulting fragment mass. Maio and Viel (2014) improved upon these aspects and found that Population III star formation was substantially suppressed at high redshifts.

WDM particles consisting of sterile neutrinos also suppress small-scale power, but have the additional effect that they decay into X-rays, which ionize the gas and enhance the abundance of H2 molecules (Biermann and Kusenko 2006; Stasielak et al. 2007). This may facilitate Population III star formation in minihalos. The influence of a running spectral index on the number density of minihalos was investigated by Somerville et al. (2003) and Yoshida et al. (2003b). In these models, the power-law index of the primordial power spectrum decreases with increasing wavenumber. For a plausible ’running’ of the spectral index, these studies found that the number density of minihalos was suppressed by about two orders of magnitude at \(z=20\). Another possible deviation from the standard CDM power spectrum may come from non-Gaussianities. However, for realistic values of the dimensionless coupling constant \(f_{\mathrm{NL}}\), Maio and Iannuzzi (2011) and Maio and Khochfar (2012) found that the properties of minihalos change only by a few percent. In quintessence models, the equation-of-state parameter w is as a function of redshift, and is considered to decrease from a value \(w>-1\) at \(z>0\) to \(w=-1\) at \(z=0\), which results in enhanced small-scale power at high redshifts. For quintessence models that do not violate other cosmological constraints, the number density of minihalos may increase by up to an order of magnitude at \(z=20\) (Maio et al. 2006).

5 From the first stars to the first galaxies

5.1 Definition

The first stars are unambiguously associated with the first DM halos in which primordial gas is able to collapse and cool. The ‘first galaxies’, however, are more difficult to define (see Bromm and Yoshida 2011). Since the term ‘galaxy’ refers to an association of stars in a gravitationally bound system, a minihalo hosting a binary stellar system may already be considered a first galaxy. This definition may also be favorable from an observational standpoint, since the stellar radiation from star-forming minihalos may be spectroscopically indistinguishable from that from more massive halos. An alternative definition is based on the transition induced by the onset of atomic hydrogen cooling in halos with virial temperatures \(\gtrsim10^{4}\mbox { K}\). In these halos, the onset of cooling does not depend on the presence on molecular hydrogen. In addition, gas photoionized and heated by stellar radiation remains gravitationally bound to the halo, which is generally not the case in minihalos. The self-sustaining cycle of star formation and feedback that is associated with galaxies can therefore operate in these halos.

Since this review is primarily concerned with the theory of primordial star and galaxy formation, I will use the definition involving the threshold for atomic hydrogen cooling. In this sense, the terms ‘first galaxy’ and ‘atomic cooling halo’ both refer to halos with virial temperatures \(\gtrsim10^{4}\mbox { K}\). The corresponding relation between virial mass and virial temperature may be obtained via equation (2):

$$ M_{\mathrm{vir}}\simeq5\times10^{7}\mbox { M}_{\odot}\biggl( \frac{T_{\mathrm{vir}}}{10^{4}\mbox { K}} \biggr)^{3/2} \biggl(\frac{1+z}{10} \biggr)^{-3/2}. $$
(9)

The characteristic virial mass of an atomic cooling halo is therefore \(10^{7}\lesssim M_{\mathrm{vir}}\lesssim10^{8}\mbox { M}_{\odot}\), with a typical formation redshift of \(z\simeq10\mbox{-}15\) for \(2\mbox{-}3\sigma\) peaks.

5.2 Turbulence

The virialization of primordial gas in atomic cooling halos without prior star formation in minihalos has been investigated by Wise and Abel (2007a) and Greif et al. (2008). They found that the gas accreted onto the halo first shock-heats to the virial temperature, after which Ly-α cooling is activated and virial equilibrium is attained via turbulence. At higher densities, H2 line cooling takes over and the turbulence becomes supersonic with Mach numbers \(M\simeq5\). This is an important difference to minihalos, where the gas is at most mildly supersonic. Prieto et al. (2011) showed that this turbulence leads to the formation of a number of gravitationally bound clumps within the halo. Wise and Abel (2007a) and Greif et al. (2008) also found that two distinct modes of accretion exist. In the standard ‘hot accretion’ mode, the gas accretes nearly radially onto the halo and shock-heats to the virial temperature near the virial radius (Birnboim and Dekel 2003). This is accompanied by a ‘cold accretion’ mode, in which cold intergalactic gas accumulates in filaments before streaming onto the halo (Kereš et al. 2005; Dekel and Birnboim 2006). In the simulations of Wise and Abel (2007a), these streams reach down to about 25% of the virial radius, while they penetrate the center of the halo in Greif et al. (2008). In the latter study, the prominence of cold streams may have been overestimated due to the aggressive accretion of the gas by a massive BH at the center of the halo, and the artificial suppression of mixing (Nelson et al. 2013; Fernandez et al. 2014).

5.3 Radiative feedback

Analytic considerations as well as simulations have indicated that the first stars had masses \(M_{*}\sim100\mbox { M}_{\odot}\), possibly with a large scatter around this value. For the purpose of stellar evolution, I therefore only consider massive Population III stars, even though calculations for their low-mass counterparts exist (Chieffi et al. 2001; Goriely and Siess 2001; Siess et al. 2002; Gil-Pons et al. 2005, 2007; Suda et al. 2007; Lawlor et al. 2008; Lau et al. 2008; Mocák et al. 2010). Massive Population III stars have low opacities due to the absence of metals, and ignite nuclear burning at very high temperatures as a result of inefficient proton-proton and CNO burning (e.g. El Eid et al. 1983; Bond et al. 1984; Marigo et al. 2001). They are therefore expected to be smaller and hotter than Population I/II stars of the same mass. The spectral shape of the radiation emitted by Population III stars on the main sequence may be derived by combining stellar structure calculations with detailed LTE and non-LTE model atmospheres (Cojazzi et al. 2000; Tumlinson and Shull 2000; Bromm et al. 2001b; Schaerer 2002, 2003). These studies concluded that massive Population III stars radiate approximately as blackbodies with an effective temperature of \({\simeq}10^{5}\mbox { K}\), and produce up to an order of magnitude more UV photons per stellar baryon than normal stars. Depending on mass, they emit most of their radiation in the LW bands in the range of \({\simeq}11.2\mbox{-}13.6\mbox { eV}\), or above the H  i, He i, or He ii ionizing thresholds of approximately 13.6, 24.6, and 54.4 eV, respectively. They are also strong emitter of X-ray radiation. Although important for the 21-cm signal and the opacity of the IGM, I will here not discuss Ly-α radiation (for a review, see Dijkstra 2014).

5.4 Photodissociating radiation

One possible way to prevent the formation of H2 is to photodetach H, which is one of the main reaction partners in forming H2. However, since the reaction rate for associative detachment is so large, this process is usually not important (but see Chuzhoy et al. 2007; Glover 2007; Wolcott-Green and Haiman 2012). Furthermore, the direct dissociation of H2 via radiative excitation to the vibrational continuum is highly forbidden. Most of the H2 is therefore dissociated by the two-step Solomon process (Field et al. 1966; Stecher and Williams 1967). Radiation in the LW bands excites a higher electronic state of H2, which is followed by decay to the vibrational continuum in approximately 15% of all cases. Detailed radiative transfer calculations have shown that the optically thin dissociation rate can be written as \(k_{\mathrm{diss}}\simeq10^{8}F_{\mathrm{LW}}\mbox { s}^{-1}\), where \(F_{\mathrm{LW}}\) is the average flux density in the LW bands in \(\mathrm{erg}\mbox { s}^{-1}\mbox { cm}^{-2}\mbox { Hz}^{-1}\) (Draine and Bertoldi 1996; Abel et al. 1997).

Initial studies showed that the LW radiation from a single massive Population III star is sufficient to prevent further cooling and star formation in the halo in which the star formed (Omukai and Nishi 1999; Nishi and Tashiro 2000; Glover and Brand 2001). Subsequent work concentrated on the effects of LW radiation on cosmological scales, finding that the H2 in the IGM with a fractional abundance of only \(y_{ \mathrm{H}_{2}}\sim10^{-6}\) is quickly dissociated. The optically thin IGM then becomes permeated by a global LW background, which may suppress star formation in halos with virial temperatures below \({\simeq}10^{4}\mbox { K}\) well before reionization (Haiman et al. 1997, 2000; Kepner et al. 1997; Ciardi et al. 2000a; Kitayama et al. 2001). However, only a modest intensity of the order of \(J_{\mathrm{LW}}=10^{-23}\mbox { erg}\mbox { s}^{-1}\mbox { cm}^{-2}\mbox { Hz}^{-1}\mbox { sr}^{-1}\) is built up by \(z\simeq20\), which is not sufficient to quench star formation in minihalos (Susa and Umemura 2000; Kitayama et al. 2001; Ricotti et al. 2002; Yoshida et al. 2003c; Greif and Bromm 2006; MacIntyre et al. 2006; Johnson et al. 2008; Trenti and Stiavelli 2009). Instead, the minimum virial mass required for efficient H2 cooling increases relatively slowly (Machacek et al. 2001; Yoshida et al. 2003c; Mesinger et al. 2006; Johnson et al. 2007; Wise and Abel 2007b; O’Shea and Norman 2008; Latif et al. 2014d; Visbal et al. 2014c). Unless the LW flux is extremely high, H2 cooling may even become important in atomic cooling halos (Susa and Kitayama 2000; Omukai 2001; Omukai and Yoshii 2003; Oh and Haiman 2002; Wolcott-Green et al. 2011; Safranek-Shrader et al. 2012).

Even though the LW flux is nearly uniform when averaged over cosmological distances, the flux seen by individual halos may fluctuate by many orders of magnitude due to their large formation bias (Dijkstra et al. 2008; Ahn et al. 2009; Holzbauer and Furlanetto 2012). This may allow the flux to become high enough to completely suppress the formation of molecules in atomic cooling halos, possibly resulting in the formation of direct collapse BHs (see Section 5.11). The combination of photodissociating and photoionizing radiation may enhance the formation of H2 in a thin shell ahead of an H ii region, which may trigger cooling and collapse in nearby halos (Haiman et al. 1996a; Ricotti et al. 2001; Kitayama et al. 2004; Ahn and Shapiro 2007; Susa and Umemura 2006; Susa 2007; Susa et al. 2009; Whalen et al. 2008a, 2010). Since this positive feedback effect requires fine-tuning of various relevant parameters, such as the distance of the star and the central density of the halo, this scenario is likely not of cosmological significance. Another argument against a positive feedback effect was provided by Glover (2007), who showed that recombination radiation from the H ii region may suppress H2 formation in the thin shell ahead of the H ii region due to the dissociation of H and \(\mbox{H}_{2}^{+}\).

5.5 Ionizing radiation

Photoionizing radiation from Population III stars has a strong effect on the gas in the halos in which they form. The first calculations that investigated the propagation of ionizing radiation from Population III stars in minihalos used a simple dynamical model, but solved the radiative transfer accurately (Kitayama et al. 2004; Whalen et al. 2004). They found that the ionizing radiation is initially trapped well within the halo by a D-type ionization front, which drives a hydrodynamic shock with a speed of \({\simeq}30\mbox { km}\mbox { s}^{-1}\) that begins to blow out the gas. After \({\simeq}10^{5}\mbox { yr}\), the shock has nearly evacuated the halo, and the ionization front becomes R-type and propagates into the IGM. The resulting density and velocity profiles are comparable to the self-similar solutions for the champagne flows discussed in Shu et al. (2002). Later three-dimensional simulations qualitatively confirmed these results, and found that \({\simeq}100\mbox { M}_{\odot}\) Population III stars create H ii regions with sizes of a few kpc, and photon escape fractions that approach unity (O’Shea et al. 2005; Alvarez et al. 2006; Abel et al. 2007). In the simulation of Greif et al. (2009b), the H ii region breaks out anisotropically due to the presence of a disk. Yoshida et al. (2007a) included helium-ionizing radiation and found that a significant He ii region with a temperature in excess of \(3\times10^{4}\mbox { K}\) develops.

The relic H ii region left behind after the star fades away cools faster than it recombines. The electron fraction therefore remains high even after the gas has cooled to below \({\simeq}10^{4}\mbox { K}\). The elevated electron fraction facilitates the formation of H2 to a level of \(y_{ \mathrm{H}_{2}}\simeq10^{-3}\), which allows the gas to cool to temperatures where chemical fractionation occurs. The enhanced HD abundance then facilitates cooling to the temperature of the CMB (Johnson and Bromm 2007; Johnson et al. 2007; Yoshida et al. 2007a, 2007b). As discussed in Section 4.1, metal-free stars that form in relic H ii region gas may have a lower characteristic mass than stars that form in minihalos unaffected by radiation. The time scale for relic H ii region gas to re-collapse is of order the Hubble time, such that continued star formation in a minihalo must await the virialization of larger halos.

The influence of ionizing radiation on neighboring halos has been investigated as well. The sign of the feedback depends on various parameters, such as the state of the collapse and the distance to the source. If a halo is close to an ionizing source or has not yet collapsed to high densities, the halo may be photoevaporated, while in other cases the halo may survive and continue to collapse (Kitayama and Ikeuchi 2000; Kitayama et al. 2000, 2001; Susa and Kitayama 2000; Susa and Umemura 2000, 2004, 2006; Susa et al. 2009; O’Shea et al. 2005; Mesinger et al. 2006, 2009; Whalen et al. 2008a, 2010; Hasegawa et al. 2009). Whalen and Norman (2008) and Vasiliev et al. (2012b) also showed that shadow and thin-shell instabilities may develop in the ionization fronts. Once a pervasive UV background has been established, star formation in minihalos will ultimately be shut down.

The UV radiation of massive Population III stars begins the process of reionization (e.g. Gnedin and Ostriker 1997; Haiman and Loeb 1997; Ciardi et al. 2000b, 2001; Wyithe and Loeb 2003; Shapiro et al. 2004; Sokasian et al. 2004; Wyithe and Cen 2007; Alvarez et al. 2012; Ahn et al. 2012). While minihalos likely did not contribute significantly to the total ionizing photon budget at high redshifts, the first galaxies are expected to have been much more important (Barkana and Loeb 2001; Cen 2003a, 2003b; Ciardi et al. 2003; Oh and Haiman 2003; Ricotti and Ostriker 2004a; Furlanetto and Loeb 2005; Greif and Bromm 2006; Haiman and Bryan 2006; Johnson et al. 2009; Trenti and Stiavelli 2009; Venkatesan and Benson 2011). Simulations that model the propagation of ionizing radiation from individual Population III star clusters typically find very high escape fractions that are of the order of unity (Wise and Abel 2008a, 2008b; Wise and Cen 2009; Wise et al. 2014; Razoumov and Sommer-Larsen 2010; Paardekooper et al. 2013). The ionizing radiation may even suppress star formation in atomic cooling halos (Barkana and Loeb 1999; Thoul and Weinberg 1996; Gnedin 2000; Kitayama and Ikeuchi 2000; Kitayama et al. 2001; Dijkstra et al. 2004). Another interesting effect stems from minihalos that are not massive enough to host star formation. These become increasingly common as reionization proceeds, and may act as sinks of ionizing radiation (Haiman et al. 2001; Barkana and Loeb 2002; Iliev et al. 2005a, 2005b; Ciardi et al. 2006).

5.6 X-rays

Although the direct emission of X-rays from Population III stars is likely cosmologically unimportant (e.g. Venkatesan and Benson 2011), they may indirectly contribute to the build-up a pervasive X-ray background. One of these sources is the accretion of gas onto the compact remnants of Population III stars. These ‘miniquasars’ may pre-heat and ionize the IGM due to the long mean free path of photons with energies \(\ge1\mbox { keV}\) (Haiman and Loeb 1999; Venkatesan et al. 2001; Glover and Brand 2003; Machacek et al. 2003; Madau et al. 2004; Ricotti and Ostriker 2004b; Ricotti et al. 2005; Salvaterra et al. 2005; Tanaka et al. 2012). They may also increase the intergalactic H2 abundance by more than an order of magnitude and reduce the clumping of the IGM (e.g. Kuhlen and Madau 2005). Similar to the evolution of relic H ii regions, the enhanced H2 abundance may lead to chemical fractionation of HD (Hummel et al. 2014). On smaller scales, the radiation pressure and heating from the quasar reduces the surrounding density and impedes the accretion flow onto the BH (Johnson and Bromm 2006; Alvarez et al. 2009; Milosavljević et al. 2009a, 2009b; Park and Ricotti 2011, 2012; Aykutalp et al. 2013). It is therefore unlikely that miniquasars will grow fast enough to explain the presence of super-massive BHs at \(z\gtrsim6\) (Fan et al. 2006). However, their radiation may delay star formation in neighboring minihalos until the atomic cooling threshold is surpassed (Alvarez et al. 2009; Jeon et al. 2012).

Another source of X-rays is Roche lobe overflow in a binary system. For massive Population III stars, the collapse of the star into a BH is more likely than in the Population I/II case (e.g. Heger et al. 2003). In addition, a number of studies have shown that a significant fraction of Population III stars in minihalos may have formed in binaries (Saigo et al. 2004; Machida et al. 2008b, 2009b; Machida 2008; Turk et al. 2009; Stacy et al. 2010; Greif et al. 2013; Stacy and Bromm 2013). Although the spectrum is different, the effects of the radiation from X-ray binaries are similar to those of miniquasars (Power et al. 2009, 2013; Jeon et al. 2012, 2014a; Xu et al. 2014). X-rays may also be produced by Population III SNe as a result of thermal bremsstrahlung or inverse Compton scattering of relativistic electrons off the CMB (Oh 2001; Glover and Brand 2003).

5.7 Final fates of Population III stars

Massive Population III stars burn their nuclear fuel very quickly and live only a few million years (e.g. Bond et al. 1984). Models of non-rotating Population III stars have shown that in the mass range \(\simeq140\mbox{-}260\mbox { M}_{\odot}\), a so-called pair-instability SN disrupts the entire star (Fryer et al. 2001; Heger and Woosley 2002; Heger et al. 2003; Joggerst and Whalen 2011; Baranov et al. 2013; Chen et al. 2014d, 2014b). In this case, the center of the star loses pressure support due to the creation of electron-positron pairs. This leads to explosive nucleosynthesis, which produces a metal yield of 50% and a kinetic explosion energy of up to 1053 ergs. In the range \({\simeq}100\mbox{-}140\mbox { M}_{\odot}\), pulsational instabilities drive episodic outbursts, while in the range \({\simeq}40\mbox{-}100\mbox { M}_{\odot}\) the entire star collapse to a BH. Below \({\simeq}40\mbox { M}_{\odot}\), a SN with \({\sim}10^{51}\mbox { ergs}\) partially disrupts the star. Depending on mass, a fraction of the star collapses into a BH, which leads to an elemental segregation of the nucleosynthetic products (Chieffi and Limongi 2002, 2004; Umeda and Nomoto 2002, 2003, 2005; Iwamoto et al. 2005; Tominaga et al. 2007b; Zhang et al. 2008; Joggerst et al. 2009, 2010b; Heger and Woosley 2010; Limongi and Chieffi 2012). For masses \(\gtrsim260\mbox { M}_{\odot}\), the entire star collapses directly to a BH without any significant SN explosion (but see Ohkubo et al. 2006; Inayoshi et al. 2013). Finally, super-massive stars with \({\simeq}10^{5}\mbox{-}10^{6}\mbox { M}_{\odot}\) may form in atomic cooling halos in which previous star formation was suppressed. A general relativistic instability develops and a fraction of the star may collapse into a BH with \({\simeq}10^{4}\mbox{-}10^{5}\mbox { M}_{\odot}\) (Heger et al. 2003; Begelman et al. 2006, 2008; Begelman 2010; Volonteri and Begelman 2010; Montero et al. 2012; Hosokawa et al. 2012a, 2013; Volonteri 2012; Inayoshi et al. 2013; Schleicher et al. 2013; Chen et al. 2014c). Recent studies have found that a super-massive star may also trigger an extremely energetic SN explosion with up to 1055 ergs of kinetic energy (Whalen et al. 2013d, 2013c, 2013h). The various fates of Population III stars are illustrated in Figure 9.

Figure 9
figure 9

Final fates of Population III stars in the absence of stellar rotation. Above \({\simeq}30\mbox { M}_{\odot}\), a fraction of the star may collapse to a BH, while at higher masses it collapses directly to a BH or explodes as a pair-instability SN. The chemical yields of the ejecta depend sensitively on the fate of the star. The mass limits change if rotation is taken into account. Adapted from Heger and Woosley (2002).

Most of the above studies have neglected the effects of rotation. Models including rotation show that metals may be mixed between nuclear burning layers and even to the surface of the star. This may have an effect on the evolution of the stars, the degree to which their elements are mixed (Meynet and Maeder 2002b, 2002a; Meynet et al. 2006; Heger et al. 2005; Chiappini et al. 2006, 2008; Hirschi 2007; Ekström et al. 2008; Takahashi et al. 2014), and their final fates (Suwa et al. 2007b; Joggerst et al. 2010a; Chatzopoulos and Wheeler 2012a, 2012b; Yoon et al. 2012; Chatzopoulos et al. 2013). For example, Chatzopoulos and Wheeler (2012a) found that strongly rotating Population III stars may already explode as pair-instability SNe if their mass exceeds \({\simeq} 65\mbox { M}_{\odot}\). If a significant amount of rotation persists until the collapse phase, the SN explosion may be accompanied by a gamma-ray burst (Fryer et al. 2001; Maeda and Nomoto 2003; Yoon and Langer 2005; Yoon et al. 2006; Tominaga et al. 2007a; Tominaga 2009; Komissarov and Barkov 2010; Suwa and Ioka 2011; Nagakura et al. 2012; Smidt et al. 2014).

5.8 Mechanical feedback

The kinetic energy released by the SN of a massive Population III star can have a substantial effect on the halo in which the progenitor star formed. Kitayama and Yoshida (2005) and Whalen et al. (2008b) used one-dimensional calculations to investigate the evolution of SN remnants in minihalos with various explosion energies. They found that more conventional core-collapse SNe fail to remove the gas from more massive halos, while pair-instability SNe are able to completely evacuate even the most massive halos. The expansion of the remnant also depends sensitively on the presence of an H ii region, which reduces the central density prior to the explosion. Numerical simulations showed that the remnant of an energetic pair-instability SN expands to a maximum radius of a few kpc, which is comparable to the size of the H ii region created by the progenitor star (Bromm et al. 2003; Greif et al. 2007; Wise and Abel 2008b; Seifried et al. 2014). In nearly all cases, the expansion can be divided into three distinct phases. In the free expansion phase, the momentum of the swept-up gas is negligible compared to that of the remnant. Once the inertia of the ambient medium becomes important, the remnant enters the energy-conserving Sedov-Taylor phase (Taylor 1950; Sedov 1959). Finally, radiative losses facilitate the transition to the momentum-conserving snowplow phase, where the expansion is driven solely by the inertia of remnant. One of the most important coolants in the final phase is inverse Compton scattering. This process is particularly important in the high-redshift IGM, due to the strong dependence of the cooling rate on the temperature of the CMB.

The chemical and thermal evolution of the gas in the SN remnant is similar to that in relic H ii regions. The elevated electron abundance facilitates the formation of H2 and HD, which may trigger secondary Population III star formation once the gas recollapses. For highly energetic explosions, the time required for the gas to recollapse is of the order of the Hubble time (Greif et al. 2010), while for less energetic explosions the collapse time is \({\simeq}10\mbox { Myr}\) (Ritter et al. 2012). Idealized simulations have also found that fragmentation in the dense shell swept up by the SN remnant may occur (Salvaterra et al. 2004; Machida et al. 2005; Vasiliev et al. 2008; Nagakura et al. 2009; Chiaki et al. 2013b). The influence of the remnant on neighboring halos depends primarily on the distance of the halo from the progenitor star and their density (Greif et al. 2007; Sakuma and Susa 2009; Whalen et al. 2010). If they are close enough and have not yet collapsed, star formation will be delayed or suppressed, while in rare cases the shock wave may compress the halo and trigger collapse. Recent studies also investigated the explosion of super-massive stars with up to 105 M in atomic cooling halos that remained metal-free (Johnson et al. 2013b; Whalen et al. 2013c, 2013d). These studies employed various methods to evolve the SN remnant through the distinct stages, and investigated their impact on the progenitor halo. With a kinetic energy of up to 1055 ergs, these SNe were able to completely disrupt their host halos.

5.9 Chemical enrichment

In addition to their mechanical feedback, SN explosions from the first stars enrich the IGM with metals (Madau et al. 2001; Scannapieco et al. 2001, 2002, 2003; Mori et al. 2002; Mackey et al. 2003; Wada and Venkatesan 2003; Ricotti and Ostriker 2004a; Scannapieco 2005). The chemical yield depends sensitively on the type of the SN. For example, a pair-instability SN mainly produces elements with an even nuclear charge, and almost no neutron-capture elements, while more conventional core-collapse SNe display a characteristic enhancement of α-elements. The enrichment pattern of the IGM also depends on the SN explosion energy and on how many SN remnants overlap prior to the re-collapse of the enriched gas. Greif et al. (2008) found that of order 10 star-forming minihalos merge to form a first galaxy, which implies that a similar number of SN ejecta mix with each other prior to second-generation star formation. Based on numerical simulations that modeled the formation and explosion of isolated Population III stars in minihalos, a number of studies found that the gas is quickly enriched to a metallicity of \(Z\sim10^{-3}\mbox { Z}_{\odot}\) (Wise and Abel 2008b; Greif et al. 2010; Ritter et al. 2012; Wise et al. 2012; Vasiliev et al. 2012a; Chen et al. 2014a). In a recent study, Ritter et al. (2014) employed tracer particles at very high resolution to follow the evolution of metal-enriched gas. They found that mixing in the IGM is suppressed due to the long eddy turnover time, while the turbulence associated with the virialization of the underlying DM halo facilitates complete mixing within the halo. In addition, the differing yields of the various SN mass shells are reflected in the enrichment pattern of the re-collapsing gas, which prevents a one-to-one mapping between the nucleosynthetic yield of the SN and the stars that form from its remnants.

The metal-enriched gas tends to re-collapse on a time scale of \(10\mbox{-}100\mbox { Myr}\) as the underlying atomic cooling halo virializes, such that the transition to Population I/II star formation occurs very rapidly (Jeon et al. 2014b). In fact, studies that investigated metal enrichment on cosmological scales found that the global star formation rate is dominated by normal stars already at \(z\sim20\) (Greif and Bromm 2006; Trenti and Stiavelli 2009; Crosby et al. 2013; Johnson et al. 2013a; Muratov et al. 2013b; Xu et al. 2013). However, due to the high spatial bias of minihalos, the enrichment of the IGM proceeds very anisotropically, with pockets of Population III star formation surviving to very low redshifts (Norman et al. 2004; Tornatore et al. 2007; Ricotti et al. 2008; Maio et al. 2010, 2011a; Muratov et al. 2013a, 2013b). Figure 10 shows the patchy enrichment of the IGM in the simulation of Wise et al. (2012). Earlier models assumed that the ejected metals promptly enriched the IGM and established a bedrock metallicity of the order of 10−5 Z (Oh et al. 2001b; Schneider et al. 2002; Mackey et al. 2003; Venkatesan and Truran 2003; Fang and Cen 2004; Ricotti and Ostriker 2004a; Yoshida et al. 2004; Matteucci and Calura 2005; Greif and Bromm 2006). Despite the progress made, the degree to which the metals mix with primordial gas is not yet fully understood (Pan et al. 2013), and depends on the employed sub-grid model (Greif et al. 2009a; Ritter et al. 2012, 2014).

Figure 10
figure 10

Three-dimensional simulation that models the transition from Population III to Population I/II star formation on cosmological scales. The physical model includes ionizing radiation from various stellar populations, and the mechanical energy input and chemical enrichment from SNe. The columns denote the redshift, and the rows the mass density, temperature, and metallicity originating from Population III and Population I/II stars, respectively. The box size is 1 Mpc (comoving). Star formation takes place in halos with masses between 106 and \({\simeq}10^{9}\mbox { M}_{\odot}\), and gradually enriches the IGM with metals. Adapted from Wise et al. (2012).

5.10 Critical metallicity

The nature of the transition from Population III to Population I/II star formation remains a matter of debate. Bromm et al. (2001a) argued that the fine-structure cooling provided by carbon and oxygen allows the gas to cool to the temperature of the CMB once the gas metallicity exceeds \(Z\sim10^{-3}\mbox { Z}_{\odot}\), which reduces the characteristic Jeans mass of the cloud and facilitates fragmentation (see also Bromm and Loeb 2003b; Santoro and Shull 2006; Safranek-Shrader et al. 2010). Later simulations included H2 cooling and confirmed that metal line cooling reduces the fragment mass (see Figure 11; Smith and Sigurdsson 2007; Smith et al. 2008, 2009; Safranek-Shrader et al. 2014b). Other studies argued that H2 cooling is just as important as fine-structure cooling at the relevant densities and temperatures, and that there is no clear distinction in the resulting fragment masses (Jappsen et al. 2007, 2009a, 2009b). However, these studies were carried out at very high redshifts where the CMB temperature was close to the minimum temperature that may be reached via H2 line cooling. Large differences in the cooling rates thus did not have a substantial effect on the thermodynamic evolution of the clouds. In addition, metal fine-structure cooling is relatively insensitive to the strength of the UV background, while molecules may be easily destroyed at the redshifts at which the first metal-enriched stars form. This effect was confirmed by Bovino et al. (2014a), who found that halos with a metallicity greater than \({\simeq}10^{-3}\mbox { Z}_{\odot}\) cooled to the temperature of the CMB in the presence of a strong UV background that dissociated the H2.

Figure 11
figure 11

Thermodynamic evolution of metal-enriched gas in an atomic cooling halo at \(\pmb{z\simeq16}\) . The panels show the results for metallicities \(Z=10^{-4}\), 10−3, and 10−2 Z, respectively (top to bottom). The gas mass per bin over the entire mass in the box is color-coded from blue (lowest) to red (highest). The dashed red lines show the temperature of the CMB, and the solid black lines on the right-hand side denote the threshold for sink particle formation. In the bottom panel, lines of constant Jeans mass are indicated as well. As the metallicity increases, metal fine-structure cooling allows the gas to cool to lower temperatures. This decreases the characteristic Jeans mass of the fragments, and facilitates the transition from Population III to Population I/II star formation. Adapted from Safranek-Shrader et al. (2014b).

Another possible trigger for a transition to normal star formation is dust cooling. In this case, the critical metallicity is expected to be as low as \(Z\sim10^{-6}\mbox { Z}_{\odot}\). At high redshifts, dust is thought to be produced in the SN remnants of Population III stars (Todini and Ferrara 2001; Nozawa et al. 2003, 2006, 2007, 2012, 2014; Salvaterra et al. 2004; Gall et al. 2011a, 2011b). A characteristic mass of \({\sim}1\mbox { M}_{\odot}\) may only be obtained for dust cooling, since the dip in the effective equation of state lies at much higher densities than for fine-structure cooling (Omukai 2000; Schneider et al. 2003, 2006, 2012; Omukai et al. 2005, 2010; Tsuribe and Omukai 2008; Schneider and Omukai 2010; Chiaki et al. 2013a, 2014). Tsuribe and Omukai (2006) modeled dust cooling in three-dimensional simulations using idealized initial conditions and found that the central clump becomes elongated and fragments for \(Z\gtrsim10^{-6}\mbox { Z}_{\odot}\). Clark et al. (2008) started from more realistic initial conditions, employed sink particles, and used a tabulated equation of state to model the effects of dust cooling. They found that the number of fragments greatly increases for \(Z\gtrsim 10^{-5}\mbox { Z}_{\odot}\). Dopcke et al. (2011) and Dopcke et al. (2013) improved upon the simulations of Clark et al. (2008) by explicitly modeling the dust temperature, and found that fragmentation occurs at all metallicities, but is much more prominent if dust cooling becomes effective. Meece et al. (2014) investigated the combined effects of metal-line cooling and grain-catalyzed molecule formation, finding that the gas temperature approaches the CMB temperature at densities that decrease as the metallicity is increased, due to metal fine-structure cooling and the formation of H2 on dust grains (see also Latif et al. 2012). They also found that the amount of fragmentation increases with increasing metallicity. In a simulation that started from cosmological conditions and included metal fine-structure cooling as well as dust cooling, Safranek-Shrader et al. (2014a) found that realistic turbulent velocities in the first galaxies enhance the density contrast well above that expected from monolithic collapse models, reducing the characteristic fragment mass to \({\simeq}0.1\mbox { M}_{\odot}\).

The above studies show that the transition from Population III to Population I/II star formation is not only governed by the critical metallicity. Among many other factors, the mass function of the first metal-enriched stars depends on the initial conditions, the temperature of the CMB, the radiation background (Aykutalp and Spaans 2011), the metallicity, the elemental composition of the metals, and the dust depletion factor. Despite this complexity, metal fine-structure cooling typically results in fragment masses of the order of 10 M, while dust cooling can lead to the formation of sub-solar or solar-mass fragments, since it operates at significantly higher densities.

5.11 Direct collapse black holes

If the local LW flux is high enough, the formation of H2 in minihalos may be suppressed until the halo has become massive enough for Ly-α cooling to become important (Omukai 2001; Bromm and Loeb 2003a; Volonteri and Rees 2005; Spaans and Silk 2006; Schleicher et al. 2010b; Johnson et al. 2013c). The gas may then contract isothermally at a temperature of \({\simeq}10^{4}\mbox { K}\) to \(n_{\mathrm{H}}\simeq10^{6}\mbox { cm}^{-3}\), when the Ly-α radiation becomes trapped (e.g. Latif et al. 2011). At this point, two-photon emission and H continuum cooling take over and extend the near-isothermal collapse phase to \(n_{\mathrm{H}}\simeq10^{16}\mbox { cm}^{-3}\) (Omukai 2001). Star formation in the progenitor halos of the atomic cooling halo must be suppressed, since the gas would otherwise become enriched with metals and have a low-temperature coolant (Johnson et al. 2008; Omukai et al. 2008).

Initial calculations showed that the critical LW flux required to suppress H2 cooling for a photospheric temperature of 105 K corresponding to Population III stars is of the order of 103 in units of \(J_{21}\), where \(J_{21}=10^{-21}\mbox { erg}\mbox { s}^{-1}\mbox { cm}^{-2} \mbox{Hz}^{-1}\mbox { sr}^{-1}\) (Omukai 2001). A flux of this level dissociates H2 up to a density of \(n_{\mathrm{H}}\simeq 10^{4}\mbox { cm}^{-3}\), where the level populations of H2 transition to LTE and H2 cooling becomes comparatively inefficient (see Section 2.3). Wolcott-Green et al. (2011) found that the critical LW flux may be reduced by an order of magnitude if a more accurate self-shielding formula for H2 is used. The inclusion of radiation from Population I/II stars and the treatment of H dissociation further reduces the critical flux (Wolcott-Green and Haiman 2012). Shang et al. (2010) find a value of only \(J_{ 21,\mathrm{crit}}\simeq30\mbox{-}300\) for a photospheric temperature of 104 K corresponding to Population I/II stars. On the other hand, Latif et al. (2014b) employ higher resolution, a more accurate H2 self-shielding formula, and investigate halos that virialize at slightly higher redshifts. They find a critical flux of \(J_{21,\mathrm{crit}}\sim10^{3}\). Agarwal and Khochfar (2015) and Sugimura et al. (2014) use a realistic stellar spectrum, which differs significantly from a blackbody spectrum, and find that this increases the critical flux to \(J_{21,\mathrm{crit}}\sim10^{3}\), while Van Borm and Spaans (2013) argue that magnetic fields and turbulence potentially reduce \(J_{21,\mathrm{crit}}\) by an order of magnitude. Recent work has shown that the critical flux from a realistic metal-enriched population at high redshifts has a radiation temperature closer to 105 than 104 K (Sugimura et al. 2014). Latif et al. (2014a) accounted for this fact and found that hydrodynamical effects such as shock-heating and inhomogeneous collapse, which are only captured in realistic three-dimensional simulations, increase the critical LW flux to a few times 104 instead of 103. Finally, Regan et al. (2014b) have argued that the radiation from a point source instead of a uniform background may further increase the critical LW flux.

The global LW background at high redshifts is much too low to reach these values (e.g. Greif and Bromm 2006). However, the local LW flux can be boosted by many orders of magnitude near star formation sites (Dijkstra et al. 2008, 2014; Agarwal et al. 2012, 2014). The timing of the incident LW flux with respect to the state of the collapse of the atomic cooling halo is an important factor (Visbal et al. 2014b). The expected number density of halos exposed to a super-critical flux is still highly uncertain. Agarwal et al. (2012) used \(J_{ 21,\mathrm{crit}}=30\) and found 10−3 direct collapse black holes per comoving Mpc at \(z=10\), while Dijkstra et al. (2014) used \(J_{21,\mathrm{crit}}=300\) and included metal enrichment, finding a number density of \(10^{-10}\mbox{-}10^{-5}\mbox { Mpc}^{-3}\) (comoving) at \(z=10\), respectively. These numbers seem more realistic in light of recent work showing that the critical flux is closer to \(J_{21,\mathrm{crit}}=10^{4}\) (e.g. Latif et al. 2014a; Sugimura et al. 2014; Yue et al. 2014).

Other mechanisms for suppressing H2 formation and cooling have been suggested as well. Inayoshi and Omukai (2012) proposed that cold accretion flows may penetrate deep into atomic cooling halos and shock-heat the gas at densities \(n_{\mathrm{H}}\gtrsim10^{4}\mbox { cm}^{-3}\), where H2 cooling is no longer efficient. However, Fernandez et al. (2014) demonstrate that the cold flows dissipate their energy at too low densities for the gas to enter the ‘zone of no return’, where H2 cooling becomes unimportant. Similarly, Tanaka and Li (2014) suggest that streaming velocities may suppress halo collapse and molecule formation until the atomic cooling threshold is surpassed. However, Visbal et al. (2014a) demonstrated that even in this case the density does not become high enough to suppress H2 cooling. Johnson et al. (2014) found that the ionizing flux that likely accompanies the LW radiation may boost the electron fraction and enhance molecule formation, which increases the critical LW flux. However, they noted that this should only occur in rare cases, since the IGM is more transparent to LW radiation, and the halo is likely able to shield itself from ionizing radiation. In a similar mechanism, X-rays may increase the critical LW flux by increasing the electron fraction and facilitating the formation of H2 (Inayoshi and Tanaka 2014).

If H2 cooling is indeed suppressed until the halo reaches a virial temperature of \({\simeq}10^{4}\mbox { K}\), the gas contracts nearly isothermally and becomes gravitationally unstable at a Jeans mass of \({\simeq}10^{5}\mbox{-}10^{6}\mbox { M}_{\odot}\). The angular momentum of the contracting cloud is redistributed by turbulence and bar-like instabilities, such that the collapse proceeds nearly unhindered until the cloud becomes optically thick to continuum emission and a protostar forms (Oh and Haiman 2002; Koushiappas et al. 2004; Begelman et al. 2006; Lodato and Natarajan 2006; Wise et al. 2008; Begelman and Shlosman 2009; Choi et al. 2013; Latif et al. 2013b; Prieto et al. 2013). Due to the high temperature of the gas, the time-averaged accretion rate onto the protostar is of order \(1\mbox { M}_{\odot} \mbox { yr}^{-1}\), eventually resulting in the formation of a super-massive star. Once enough mass has been accreted, a general relativistic instability develops and a fraction of the star collapses into a BH (see Section 5.7). As opposed to minihalos, photoheating is not able to significantly suppress accretion, since the virial temperature of the halo is of order the temperature to which the gas is heated, and momentum transfer by photons only mildly reduces the collapse rate (Johnson et al. 2011, 2012).

The evolution of the central cloud beyond the initial collapse of the gas was investigated in a number of studies. Regan and Haehnelt (2009) used adaptive mesh refinement simulations with a maximum resolution of \({\simeq}0.01\mbox { pc}\), 16 cells per Jeans length, and employed a pressure floor to avoid artificial fragmentation. Due to the limited resolution, the collapse stalled at a density of \(n_{\mathrm{H}}\simeq10^{9}\mbox { cm}^{-3}\), followed by the formation of a massive disk around the central hydrostatic core. In one of the three halos investigated, the cloud fragmented into three distinct clumps. Latif et al. (2013e) employed sink particles above a density of \(n_{\mathrm{H}}\simeq10^{6}\mbox { cm}^{-3}\) and found that some of the clouds were prone to fragmentation. In a follow-up study, Latif et al. (2013a) investigated nine different halos using 64 cells per Jeans length, and employed a pressure floor at a similar density as Regan and Haehnelt (2009). They found that even though fragmentation occurs, the central object continues to grow via turbulent accretion and mergers at a rate of \({\simeq}1\mbox { M}_{\odot} \mbox { yr}^{-1}\). At this rate, a super-massive star with \({\simeq}10^{6}\mbox { M}_{\odot}\) would form after only \({\simeq}1\mbox { Myr}\) (see also Inayoshi and Haiman 2014). Similar fragmentation and accretion was found in the high-resolution simulations of Regan et al. (2014a). In a complementary approach, Inayoshi et al. (2014) included the most comprehensive chemical and thermal model to date. Although they do not start from cosmological initial conditions and terminate the simulation once the density reaches \(n_{\mathrm{H}}\simeq10^{17}\mbox { cm}^{-3}\), they find a minimum Jeans mass of \({\simeq}0.2\mbox { M}_{\odot}\), which is more than an order of magnitude higher than in other star formation environments.

One of the highest resolution simulations of the collapse of gas in atomic cooling halo was presented by Becerra et al. (2014). This simulation employed a somewhat simpler chemical model than Inayoshi et al. (2014), but followed the evolution of the central gas cloud well beyond its initial collapse (see Figure 12). In analogy to previous studies, the disk fragmented into a protostellar system with 5-10 members, with the central protostar accreting at a rate of \({\simeq}1\mbox { M}_{\odot} \mbox { yr}^{-1}\). Due to the high computational cost of the simulation, the system could only be evolved for 18 yr. Nevertheless, the central protostar had already grown to \({\simeq}15\mbox { M}_{\odot}\). Near the end of the simulation, a second clump collapsed at about 150 au from the primary clump, potentially resulting in the formation of a wide binary system. However, the accretion rate onto the central protostars in both clumps remains extremely high, showing that fragmentation does not prevent the rapid growth of the central objects.

Figure 12
figure 12

Zoom-in on the gas cloud that forms at the center of an atomic cooling halo, showing the number density of hydrogen nuclei. Clockwise from the top left, the width of the individual cubes are 10 pc, 1 pc, 0.1 pc, \(1\text{,}000\mbox { au}\), 100 au, and 10 au. The cloud has an irregular morphology that continues to change shape and orientation throughout the collapse. The filamentary structure indicates that turbulence is present on all scales. Adapted from Becerra et al. (2014).

6 Empirical signatures

A number of studies have predicted that the first stars and galaxies have distinct observational signatures. One of the most promising signatures stems from Population III SNe (e.g. de Souza et al. 2014). The resulting lightcurves may show the characteristics of core-collapse SNe (Meiksin and Whalen 2013; Whalen et al. 2013b, 2013g), pair-instability SNe (Scannapieco et al. 2005; Weinmann and Lilly 2005; Hummel et al. 2012; Pan et al. 2012; de Souza et al. 2013b; Whalen et al. 2013a, 2013f, 2014), gamma-ray bursts (Bromm and Loeb 2002, 2006; Choudhury and Srianand 2002; Islam et al. 2004; Natarajan et al. 2005; Hirose et al. 2006; Belczynski et al. 2007; Inoue et al. 2007; Naoz and Barkana 2007; Salvaterra et al. 2007, 2008; Salvaterra and Chincarini 2007; Komissarov and Barkov 2010; Mészáros and Rees 2010; Campisi et al. 2011; de Souza et al. 2011, 2012; Toma et al. 2011; Nakauchi et al. 2012; Wang et al. 2012; Macpherson et al. 2013; Elliott et al. 2014; Maio and Barkov 2014; Mesler et al. 2014), or even super-massive stars (Whalen et al. 2013c, 2013d, 2013h). More exotic signatures include the neutrino emission accompanying a SN explosion (Fryer et al. 2001; Iocco et al. 2005, 2008b; Nakazato et al. 2006; Suwa et al. 2009), and gravitational waves from BH binaries (Fryer et al. 2001; Belczynski et al. 2004; Kulczycki et al. 2006; Suwa et al. 2007a; Kinugawa et al. 2014).

The properties of the first stars may also be probed by their nucleosynthetic signature, which is likely imprinted in second-generation stars that form from metal-enriched gas. If they survive to the present day, their surface abundances may reflect the metal enrichment pattern of the cloud out of which they formed. This avenue of probing the first stars has been termed ‘stellar archeology’, or ‘near-field cosmology’ (Beers and Christlieb 2005; Sneden et al. 2008; Frebel 2010; Karlsson et al. 2013). A number of stars with extremely low metallicities have been found in the Galactic halo, and their abundance patterns may reflect certain types of Population III SNe (Tsujimoto et al. 1999; Karlsson and Gustafsson 2001, 2005; Daigne et al. 2004, 2006; Frebel et al. 2005, 2007, 2009; Karlsson 2005, 2006; Venkatesan 2006; Salvadori et al. 2007; Tumlinson 2007a, 2007b; Karlsson et al. 2008; Klessen et al. 2012; Cooke and Madau 2014; Hartwig et al. 2014a; Marassi et al. 2014). Most second-generation stars are in fact expected to reside within the bulge of the Galaxy instead of the halo (Diemand et al. 2005; Tumlinson 2006, 2010; Trenti et al. 2008; Gao et al. 2010). There may also be a connection between the first stars and globular clusters (Padoan et al. 1997; Bromm and Clarke 2002; Beasley et al. 2003; West et al. 2004; Kravtsov and Gnedin 2005; Bekki 2006; Moore et al. 2006; Bekki et al. 2007; Boley et al. 2009), as well as extremely metal-poor stellar populations found in local dwarf galaxies (Ricotti and Gnedin 2005; Gnedin and Kravtsov 2006; Moore et al. 2006; Read et al. 2006; Ricotti et al. 2008; Salvadori et al. 2008; Bovill and Ricotti 2009; Muñoz et al. 2009; Ricotti 2009; Salvadori and Ferrara 2009; Frebel et al. 2010; Bovill and Ricotti 2011a, 2011b; Frebel and Bromm 2012; Karlsson et al. 2012; Milosavljević and Bromm 2014). There is even a possibility of finding true Population III stars (Greif et al. 2011a; Stacy and Bromm 2014), even though it may be impossible to distinguish them from other metal-poor stars due to self-enrichment or mass transfer in a binary system (Schlattl et al. 2001, 2002; Fujimoto et al. 2000; Weiss et al. 2000, 2004; Picardi et al. 2004; Suda et al. 2004; Lucatello et al. 2005; Lau et al. 2007, 2009; Tumlinson 2007a; Suda and Fujimoto 2010; Starkenburg et al. 2014; Johnson 2014). Low-mass Population III stars may also accrete metal-enriched gas from the IGM as they move through the Galaxy (Shigeyama and Tsujimoto 2003; Shigeyama et al. 2003; Frebel et al. 2009; Johnson and Khochfar 2011).

The radiation emitted from individual Population III stars is too faint to be detected directly (but see Zackrisson et al. 2010a, 2010b). However, the stellar emission from the first galaxies may be detected by existing and upcoming telescopes such as the Hubble Space Telescope (HST), the JWST (Gardner et al. 2006, 2009), the Atacama Large Millimeter Array (ALMA), the Square Kilometer Array (SKA), and the planned extremely large \(30\mbox{-}40\mbox { m}\) class telescopes. A number of studies have attempted to predict the characteristic signature of the first galaxies (Greif et al. 2009b; Johnson et al. 2009; Johnson 2010; Wise and Cen 2009; Raiter et al. 2010; Pawlik et al. 2011, 2013; Salvaterra et al. 2011; Zackrisson et al. 2011a, 2011b, 2013; Wise et al. 2012, 2014), with a particular emphasis on the strong He ii recombination lines characteristic of metal-free stellar populations (Tumlinson and Shull 2000; Oh et al. 2001a; Tumlinson et al. 2001, 2003; Kudritzki 2002; Schaerer 2002, 2003; Venkatesan et al. 2003). A significant boost to the received luminosity may be provided by gravitational lensing (Maizy et al. 2010; Wyithe et al. 2011; Zackrisson et al. 2012; Mashian and Loeb 2013; Pan and Loeb 2013; Rydberg et al. 2013; Whalen et al. 2013e). Finally, the first galaxies may also have contributed to the near-infrared background (Santos et al. 2002; Magliocchetti et al. 2003; Salvaterra and Ferrara 2003; Cooray et al. 2004; Cooray and Yoshida 2004; Kashlinsky et al. 2004; Dwek et al. 2005; Kashlinsky et al. 2005; Kashlinsky 2005; Fernandez and Komatsu 2006; Salvaterra et al. 2006; Salvaterra and Ferrara 2006; Cooray et al. 2012; Fernandez and Zaroubi 2013; Yue et al. 2013a, 2013b).

A number of reviews summarize the observational signatures of the first stars and galaxies. The connection between Population III stars and extremely metal-poor halo stars is reviewed by Beers and Christlieb (2005), Sneden et al. (2008), Frebel (2010), and Karlsson et al. (2013), while the properties of dwarf galaxies in the Local Group are reviewed by Tolstoy et al. (2009) and Ricotti (2010). Various reviews also focus on observations of high-redshift galaxies (Bland-Hawthorn 2006; Bouwens and Illingworth 2006; Stark and Ellis 2006; Ellis 2008; Robertson et al. 2010), while the near-infrared background is summarized in Hauser and Dwek (2001) and Kashlinsky (2005, 2009). The process of reionization is described in detail in Loeb and Barkana (2001), Fan et al. (2006), and Stiavelli (2009), and the connection between primordial star formation and the 21-cm signal is discussed in Furlanetto et al. (2006) and Morales and Wyithe (2010). Finally, the properties and observational signature of BH remnants from the first stars are discussed in Haiman (2006, 2009), Greene (2012), Volonteri (2012), and Volonteri and Bellovary (2012).

7 Summary

The numerical frontier of the high-redshift Universe has advanced considerably over the last 10-15 years. The turn of the millennium saw the first three-dimensional simulations of primordial star formation that started from cosmological initial conditions and included primordial chemistry networks (e.g. Abel et al. 2002; Bromm et al. 2002). They established the ‘standard model’ of primordial star formation, in which the physics of H2 cooling leads to a characteristic stellar mass of \({\sim}100\mbox { M}_{\odot}\). However, increasingly sophisticated simulations have begun to refine this picture. The inclusion of additional chemical and thermal processes as well as numerical and technological headway have allowed the entire collapse process to be modeled self-consistently (Yoshida et al. 2006, 2008). A key insight gained from these simulations is that primordial gas clouds are prone to fragmentation, but that the secondary protostars rapidly migrate to the center of the cloud and merge with the primary (Clark et al. 2011b; Greif et al. 2012). The most likely scenario therefore appears to be the formation of a massive central star or a binary system, surrounded by a number of significantly less massive stars.

Studies that investigated the influence of the radiation from the central protostar found that photoheating terminates accretion onto the central star and leads to the formation of Population III stars with a wide range of masses. This is a result of their varied formation environments (e.g. Hosokawa et al. 2011; Stacy et al. 2012a; Hirano et al. 2014). Recent simulations have also begun to include magnetic fields, and found that they are amplified to dynamically important levels via the turbulent dynamo during the initial collapse (e.g. Sur et al. 2010; Peters et al. 2012; Turk et al. 2012). A strong magnetic field enhances the rate of angular momentum transport and reduces the susceptibility of the gas to fragmentation (e.g. Machida and Doi 2013; Peters et al. 2014). A number of other physical processes may play a role as well. These include HD cooling, cosmic rays, streaming velocities, DM annihilation, and alternative cosmologies.

The second step in the hierarchy of structure formation is the formation of the first galaxies in atomic cooling halos. The first high-resolution simulations focused on the virialization of the gas in the host halo (e.g. Greif et al. 2008; Wise and Abel 2007a). Later on, they included star formation recipes that modeled the radiative, mechanical and chemical feedback from Population III stars in their progenitor halos (e.g. Wise and Abel 2008b; Greif et al. 2010). The UV radiation of massive Population III stars dissociates molecules and photoheats the gas to \({\simeq}10^{4}\mbox { K}\) (e.g. Alvarez et al. 2006; Johnson et al. 2007; O’Shea and Norman 2008; Whalen et al. 2010). The metals dispersed by their SNe mix with primordial gas as they recollapse into the halo of the nascent galaxy (e.g. Wise et al. 2012; Ritter et al. 2014). This leads to the formation of the first metal-enriched stellar clusters, with an IMF that resembles that of our Galaxy (e.g. Tsuribe and Omukai 2006; Dopcke et al. 2013; Safranek-Shrader et al. 2014b). The first galaxies started the process of reionization (e.g. Wyithe and Cen 2007; Ahn et al. 2012; Salvadori et al. 2014), and some may have seeded the first quasars by forming massive BHs in halos subjected to a strong LW background (e.g. Omukai 2001; Bromm and Loeb 2003a; Dijkstra et al. 2008; Latif et al. 2013a; Inayoshi et al. 2014; Regan et al. 2014a).

Thanks to advances in computer technology and simulation methods, our understanding of primordial star and galaxy formation has rapidly increased. The well-known initial conditions provided by observations of the CMB make this field particularly attractive. The underlying equations are well known, such that obtaining an accurate solution is merely a matter of complexity. Based on the current rate of progress, the first simulations of primordial star formation that include radiative transfer as well as magnetic fields will become possible within the next five years. This nicely coincides with the commissioning of the next generation of ground- and space-based telescopes, such as the upcoming 30-40 m class telescopes or the JWST.

References

Download references

Acknowledgements

The author would like to thank Volker Bromm, Anastasia Fialkov, Simon Glover, Chalence Safranek-Shrader, Dominik Schleicher, and Naoki Yoshida for comments and suggestions that helped improve the manuscript.

Author information

Authors and Affiliations

Authors

Corresponding author

Correspondence to Thomas H Greif.

Additional information

Competing interests

The author declares that he has no competing interests.

Rights and permissions

Open Access This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly credited.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Greif, T.H. The numerical frontier of the high-redshift Universe. Comput. Astrophys. 2, 3 (2015). https://doi.org/10.1186/s40668-014-0006-2

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s40668-014-0006-2

Keywords