The numerical frontier of the high-redshift Universe

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 ~10^6 M_sun. 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 ~10^8 M_sun, 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.


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 20 (Bromm & Larson 2004;Glover 2005;Bromm 2013;Glover 2013). They formed at the center of dark matter (DM) 'minihalos' with virial masses M vir ∼ 10 6 M , 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 (H 2 ). The importance of H 2 for the cooling of low-mass gas clumps that condense out of the expanding Universe was recognized in the late 1960's (Saslaw & Zipoy 1967;Peebles & 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;Carlberg 1981;Kashlinsky & Rees 1983;Palla et al. 1983;Silk 1983;Carr et al. 1984;Izotov & Kolesnik 1984;Couchman & Rees 1986;Susa et al. 1996;Uehara et al. 1996;Tegmark et al. 1997;Nishi & Susa 1999;Vasiliev & Shchekinov 2003). 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;Nakamura & Umemura 1999), while three-dimensional simulations had to await improvements in numerical simulation techniques in the late 1990's (Abel et al. 1998;Bromm et al. 1999Abel et al. 2000Abel et al. , 2002.
One of the main results of these studies is that the minimum temperature to which the gas can cool via H 2 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 * ∼ 100 M . Following the introduction of this 'standard model' of primordial star for-mation, 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 & Barkana 2001;Ciardi & Ferrara 2005;Fan et al. 2006;Stiavelli 2009). Their radiation also dissociates H 2 on cosmological scales and may suppress star formation until halos massive enough to activate atomic hydrogen cooling assemble (Ciardi & 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 & 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 & Yoshida 2011;Loeb & 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 & Abel 2007a;Greif et al. 2008). Radiative and chemical feedback from stars in progenitor halos influence their formation as well (e.g. Wise & 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 & Loeb (2001), Loeb (2008), and Wiklind et al. (2013). The high-redshift IGM is discussed in Barkana & 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 & Larson (2004), Glover (2005), Bromm (2013), and Glover (2013), while the properties of the first galaxies are described in Bromm & Yoshida (2011), Johnson (2013), and Loeb & 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 & 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 δ = (ρ −ρ) /ρ, where ρ is the local mass density andρ 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.

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 of the order of the Earth mass if DM consists of a weakly interacting massive particle (WIMP). Following the collapse of a DM halo, it achieves virial equilibrium within a few dynamical times, with the virial density given by ρ vir = ∆ virρ , where ∆ vir 200 (e.g. Barkana & Loeb 2001). The balance between gravitational potential energy and kinetic energy allows a virial velocity to be defined, and is given by v 2 vir = GM vir /R vir , where G is the gravitational constant, M vir the virial mass, and R vir the virial radius of the halo. The latter is related to the mass of the halo via R 3 vir = 3M vir / (4π∆ virρ ). Halos with virial masses up to ∼ 10 6 M are commonly referred to as minihalos, since they are massive enough to cool and host star formation. A 2-3σ peak typically forms at z ∼ 20, and their virial radius may be conveniently written as: 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 & Hernquist 2001;Heitmann et al. 2006;Lukić et al. 2007;Reed et al. 2007;Cohn & White 2008;Davis & Natarajan 2009Davis 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 vir = m H v 2 vir /k B , where m H is the mass of the hydrogen atom, and k B Boltzmann's constant. In terms of a typical minihalo, the virial temperature may be written as: T vir 2 × 10 3 K M vir 10 6 M 2/3 1 + z 20 . (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, λ J = c s t ff , where c s is the sound speed of the gas, and t 2 ff = 3π/ (32Gρ) the free-fall time. The corresponding Jeans mass is given by M J = πρλ 3 J /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 100, the cosmological Jeans mass initially remains constant at M J 10 5 M (Bromm & Larson 2004). At lower redshifts, the gas expands adiabatically, and the Jeans mass is given by M J 10 4 M [(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 & Hui 1998;. These calculations find that the filtering mass remains nearly constant at a few times 10 4 M at z 100. 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.

H 2 formation and cooling
For the gas 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 10 4 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 vir = 5 × 10 7 M , which corresponds to a virial temperature of 10 4 K at z = 10, must rely on another coolant. Saslaw & Zipoy (1967) and Peebles & Dicke (1968) realized that molecular hydrogen (H 2 ) may be such a coolant.
The rotational and vibrational states of H 2 in the electronic ground state are excited by collisions with other particles, which decay and allow the gas to cool. These transitions operate in the temperature range 200 K T 5000 K. Despite its importance at high redshifts, the H 2 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 → 0, where J is the rotational quantum number. The least energetic transition with a non-negligible probability of occurring is the J = 2 → 0 transition, which corresponds to a temperature of 500 K. In practice, the Maxwellian tail of the velocity distribution allows the gas to cool to 200 K. The absence of a dipole moment in the H 2 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 & Salpeter 1963).
For this reason, alternative formation channels have been considered (Haiman et al. 1996b;Abel et al. 1997;Galli & 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 & Dicke 1968): Associative detachment of the H − atoms with neutral hydrogen atoms then results in the formation of H 2 : One of the most important parameters that governs the amount of H 2 that can be formed is the electron abundance. It decreases from its relic abundance by recombining with ionized hydrogen: Most of the H 2 therefore forms within a recombination time. A straightforward calculation shows that the asymptotic H 2 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 H 2 abundance scales approximately as y H 2 ∝ T 1.5 , and for T = 1000 K lies in the range of y H 2 10 −4 − 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 H 2 , which allows the gas to cool more efficiently.
The H 2 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 on which the virial temperature of a halo changes significantly (Rees & Ostriker 1977;Silk 1977). A straightforward calculation shows that for a virial temperature of about 1000 K, enough H 2 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: 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 ∼ 20, they are not the very first stars in the observable Universe. Instead, these are believed to have formed at z 50 − 65 Naoz et al. 2006).

Jeans instability
During the initial stages of the collapse, the level populations of H 2 are not yet populated according to local thermal equilibrium (LTE), and the cooling rate scales as Λ H 2 ∝ n 2 H , where n H is the number density of hydrogen nuclei. Once enough H 2 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 H 2 molecules and other species become so frequent that the ro-vibrational levels of H 2 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 as Λ H 2 ∝ n 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 H n H,crit = 10 4 cm −3 , where the gas has cooled to 200 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 . The Jeans mass at this characteristic density and temperature is 100 M . Once this mass has been accreted, the cloud collapses under its own gravity. Since the H 2 cooling rate is a strong function of temperature, the cooling time in the LTE regime adjusts to the free-fall time, and the effective equation of state is related to the power-law exponent of the cooling function α via γ eff = 1 + 1/ (2α − 2). Typically, α 4, which results in γ eff 1.1 − 1.2. The temperature therefore increases very slowly with increasing density for n H n H,crit . In a near-isothermal, Jeans-unstable cloud the accretion rate is approximately given byṀ = M J /t ff ∝ 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.

Three-body H 2 formation
Following the nearly isothermal collapse of the cloud to densities n H 10 8 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 with a smaller contribution from reactions involving H 2 molecules and helium atoms. The rapidly increasing H 2 fraction results in a similarly rapid increase in the H 2 line cooling rate. However, since each formation process is accompanied by the release of the binding energy of the H 2 molecule of 4.48 eV, the cooling is offset by chemical heating. The net effect is a mild drop in temperature as the H 2 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 & Yoshii 1977;Ripamonti & Abel 2004;Nakamura & 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 H 2 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 1000 K, the published rates differ by up to two orders of magnitude (Abel et al. 2002;Palla et al. 1983;Flower & 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.

Optically thick H 2 line cooling
Another important transition occurs at densities n H 10 10 cm −3 , where the gas becomes optically thick to H 2 line emission. For each of the 200 rovibrational 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 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 & 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 esc = min (n H /n H,line ) −0.45 , 1 , where n H,line = 6 × 10 9 cm −3 . This fit has also been used in three-dimensional simulations (Turk et al. , 2010(Turk et al. , 2011, 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 & 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 esc = [1 − exp (−τ )] /τ , where τ = αL Sob , L Sob is the Sobolev length, and α the absorption coefficient. The Sobolev length is given by L Sob = v therm / |dv r /dr|, where v therm is the thermal velocity, and dv r /dr 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,b;Greif et al. 2011aGreif et al. , 2012Greif et al. , 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 & 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-consistently treatment of H 2 line transfer in a three dimensional simulation was presented in Greif (2014). This study used a multi-line, multifrequency ray-tracing scheme to follow the propagation of H 2 line radiation through the computational domain (see Figure 1). The spherically averaged escape fraction agreed relatively well with the fit of Ripamonti & 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 H 2 fraction at densities n H 10 15 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.
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 50 100 150 x [au] 1200 1600 2000 2400 50 100 150 x [au] 2.4 1.8 1.2 0.6 log f esc Figure 1: Hydrogen number density, temperature, and escape fraction of H 2 line photons in the central 200 au of a minihalo, shown for various treatments of the optically thick H 2 line cooling rate. The top row shows the outcome for the density-dependent fit of Ripamonti & 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).
number of frequency bins (the extra factor of N 1/3 cells 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 cells 2 × 10 7 , N rays = 48, N lines = 32, and N ν = 8, nearly 10 14 opacity calculations per time step were carried out. A single ray-tracing step therefore took about 2000 s on 1024 modern computing cores. To evolve the gas from n H 10 10 cm −3 to n H 10 15 cm −3 , where the gas becomes optically thick to H 2 line emission, about 2000 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 cell N rays N lines N ν bytes must be reserved, which corresponds to 1 GB per core for 1024 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. (2014). This study used the treecol algorithm of Clark et al. (2012) to determine the total column density of H 2 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. (2014) to follow the build-up and fragmentation of the disk around the primary pro-tostar. 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.

Collision-induced emission and absorption
For densities n H 10 14 cm −3 , collision-induced emission and absorption become important Ripamonti et al. 2002;Ripamonti & 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 H 2 molecules. As opposed to H 2 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 H 2 line emission. The various stages in the collapse of the gas up to this point are shown in Figure 2.
In the one-dimensional dynamical model of Ripamonti & Abel (2004), the gas becomes optically thick to collision-induced emission at densities n H 10 16 cm −3 . This rapid transition compared to H 2 line emission is due to the absence of absorption wings, which increases the optical depth of the gas. Ripamonti & Abel (2004) also found that a chemothermal instability similar in nature to the one induced by three-body H 2 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 with 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 & 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 The H 2 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 H 2 fraction remains nearly constant, the cloud becomes fully molecular once three-body reactions set in. Adapted from Yoshida et al. (2006). with one-dimensional calculations for clouds that are approximately spherically symmetric (Hirano & Yoshida 2013).

Collisional dissociation
The last process that is able to cool the gas is collisional dissociation of H 2 . It acts as a thermostat by converting the compressional energy of the gas into energy that is used for the dissociation of H 2 , which has a binding energy of 4.48 eV. The temperature and density at which H 2 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 2000 K at a density of n H 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 H 2 line cooling rate (Turk et al. 2011;Greif 2014).
From a computational standpoint, obtaining a self-consistent H 2 fraction at these densities is computationally expensive, since the three-body formation rate scales with the cube of the density. However, the H 2 abundance approaches chemical equilibrium at n H 10 15 cm −3 , and the electron abundance at n H 10 19 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 Ripamonti et al. 2002;Yoshida et al. 2006Yoshida et al. , 2008Greif et al. 2012).

Protostar formation and evolution
Following the dissociation of H 2 , the temperature rises steeply with increasing density. At densities n H 10 20 cm −3 , the collapse stalls and an accretion shock forms that heats the gas to 10 4 K Yoshida et al. 2008). This marks the formation of a protostar at the center of the cloud with an initial mass of 0.01 M and a size of 0.1 au.
The first calculations of the evolution of primordial protostars with continuous accretion were carried out by Stahler et al. (1986a,b) and Omukai & Palla (2001. 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 × 10 −3 M 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 acc = M * /Ṁ * , where M * is the mass of the protostar, is smaller than the Kelvin-Helmholtz time t 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 100 R . Once t acc > t KH , this trend is reversed and the protostar undergoes Kelvin-Helmholtz contraction, which continues until M * 50 M . Following a phase of deuterium burning, hydrogen burning begins after 10 4 yr, and the protostar settles onto the main sequence after 10 5 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  and Ohkubo et al. (2009).
One of the major uncertainties of these models is the assumption of spherical symmetry. This was relaxed in Tan & 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. Omukai & 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Ṁ ∼ M J /t ff ∼ c 3 s /G ∝ T 3/2 , where the Jeans mass is evaluated at the density and temperature when the gas first becomes gravitationally unstable, i.e. n H 10 4 cm −3 and T 200 K. Since the temperature in present-day star-forming regions is 10 K, the accretion rate in primordial gas clouds is expected to be nearly 100 times higher, withṀ ∼ 10 −3 M 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 10 yr (Ripamonti et al. 2002;Greif et al. 2012). Since primordial stars are expected to accrete for more than ∼ 10 5 yr before radiation feedback terminates accretion (Tan & McKee 2004;McKee & Tan 2008), various methods have been used to estimate the final masses of Population III stars. In , the self-similar form of the spherically symmetric infall solution for γ eff 1.1 was extended to the accretion phase, resulting iṅ M = 8.3 × 10 −2 M yr −1 (t/1 yr) −0.27 , where t denotes the time (Larson 1969;Penston 1969;Shu 1977;Yahil 1983;Suto & Silk 1988). The mass of the star after 10 5 yr is thus M * 500 M . A similar power-law for the accretion rate was found in the calculations of Ripamonti et al. (2002) and Tan & McKee (2004).

Accretion
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;O'Shea & 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 & Loeb 2004;Clark et al. 2008;Stacy et al. 2010Stacy et al. , 2011bStacy et al. , 2012aStacy & Bromm 2013;Clark et al. 2011a,b;Greif et al. 2011a;Smith et al. 2011Smith et al. , 2012a. At sufficiently late times, they agree with those found in previous studies, but show substantial variability.

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;Kudritzki 2002;Omukai & Inutsuka 2002;Krtička et al. 2003;Marigo et al. 2003;Krtička & Kubát 2006;McKee & Tan 2008;Krtička & Kubát 2009;Krtička et al. 2010;Sonoi & Shibahashi 2011;Muijres et al. 2012;Sonoi & Shibahashi 2012;Sonoi & Umeda 2012), the most important feedback mechanisms are the dissociation of H 2 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 H 2 and removes the most important coolant of the gas (Omukai & Nishi 1999;Glover & 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).
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 13.6 eV thermalizes and heats the gas to 10 4 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 & 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 & 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 10-1000 M , 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. Figure 3: Two-dimensional simulation of primordial star formation including radiation from the central protostar. The panels shows the temperature and 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 40 M . Adapted from Hosokawa et al. (2011).
They found that most Population III stars had final masses in the range 10 M * 100 M .
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.

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 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 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). They are typically inserted at a threshold density, and accrete the gas within a predefined accretion radius. This prevents the gas from collapsing to increasingly high densities, and allows the 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.
Some of the first studies that employed sink particles inserted them on scales significantly larger than the expected sizes of the protostars Bromm & Loeb 2004;Stacy et al. 2010;Stacy & 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 H = 10 16 cm −3 , which corresponds to a radius of 1 au and is the approximate size of a metal-free protostar that has a realistic accretion rate (Stahler et al. 1986a,b). The initial conditions were designed to reproduce the end state of the primordial gas clouds found in , 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 10 M . 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 & White 1994;Tormen et al. 1997;Gao et al. 2005). The resolution was increased several times with a particle-splitting technique (Kitsionas & Whitworth 2002;Bromm & Loeb 2003a). This allowed the density at which sink particles were created to be increased to n H = 10 17 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 H 2 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 H 2 cooling is inefficient at the densities and temperature at which the disk forms ( Tan (2014) found that primordial disks remain marginally stable despite the H 2 cooling. 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 H 2 line cooling rate (see Greif 2014). However, the more accurate method used in Hartwig et al. (2014) showed that this does not substantially affect the fragmentation of the gas.
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 10 3 yr, with masses ranging from 0.1 to 10 M . 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 case, 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) andSmith 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.

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 10 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 s κ/ (πGΣ). 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 cool < 3/Ω, where t cool is the cooling time, to be satisfied. Greif et al. (2012) showed that both of these conditions are fulfilled on a scale of 1 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 H 2 , which acts as a thermostat by extracting compressional energy from the gas and using it to unbind H 2 . In this respect the results of Greif et al. (2012) differ somewhat from Clark et al. (2011a), where H 2 line emission was found to be the dominant coolant. This is likely due to the higher resolution employed in Greif et al. (2012).
The dependence of the stability of the disk on the abundance of H 2 implies that differences in the three-body formation rates and the treatment of the optically thick H 2 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 H 2 line cooling rate used in Turk et al. (2012). Indeed, this trend has been confirmed by the simulations of Hartwig et al. (2014). 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.
In a qualitatively similar study, Latif et al. (2013c) artificially truncated the collapse of the gas at a density of n H 10 13 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 5 × 10 4 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 Figure 6: From left to right: the first 10 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).
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.

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;. 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 thus endowed with significant rotation, and since it accretes material from a Keplerian disk, the angular momentum of the protostar will likely remain high as it evolves towards the main sequence. 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 5000 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.

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 & Shull 1984;Puy et al. 1993;Stancil et al. 1998;Galli & Palla 1998Flower & Roueff 1999;Flower et al. 2000;Flower 2000;Flower & Pineau des Forêts 2001;Lipovka et al. 2005;Johnson & Bromm 2006;Glover & Abel 2008). As opposed to H 2 , 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 H,crit ∼ 10 6 cm −3 . Furthermore, since ortho and para states do not exist, the transition J = 1 → 0 is accessible, corresponding to an excitation temperature of 130 K. Despite the very low cosmological abundance of deuterium of a few times 10 −5 , chemical fractionation can allow HD to become a more important coolant than H 2 at very low temperatures. The primary formation reaction is: showing that the HD formation rate depends on the ionization state of the gas and the abundance of H 2 .
In minihalos, the temperature usually does not become low enough for HD cooling to become important . However, more recent studies have found that HD cooling may in fact dominate over H 2 cooling in certain halos. Ripamonti (2007) and McGreer & Bryan (2008) showed that HD cooling is activated in high-redshift minihalos with masses below 3 × 10 5 M . This allows the gas to cool down to the temperature of the CMB, which lies between 50 and 100 K at z 20. The corresponding Jeans mass is up to an order of magnitude lower than in the pure H 2 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 & 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 10 M . Clark et al. (2011a) came to a somewhat different conclusion. They found that the gas in these clouds heated more rapidly at densities n H n H,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, that the total mass in stars is likely lower in the HD cooling case.
If the electron abundance is significantly enhanced with respect to the postrecombination 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 & Bromm 2006;Greif et al. 2008;Shchekinov & Vasiliev 2006;Vasiliev & Shchekinov 2008;Prieto et al. 2012;Bovino et al. 2014b;Prieto et al. 2014), as well as in SN remnants (Mac Low & Shull 1986;Uehara & Inutsuka 2000;Mackey et al. 2003;Machida et al. 2005;Vasiliev & Shchekinov 2005a;Vasiliev et al. 2008). The H 2 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 & Omukai 2005;Yoshida et al. 2007a,b;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.

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 & Silk 1989;Beck et al. 1994;Lesch & Chiba 1995;Kulsrud et al. 1997). The magnetic field strength at the time of the formation of the first stars is constrained to B 1 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 & 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 small compared to the evolutionary time of a system, the magnetic field is 'frozen' into the gas and moves with it (but see Maki & Susa 2004. In a spherically symmetric, contracting gas cloud the magnetic field grows ∝ ρ 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 ∝ ρ −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 & 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 showed 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. 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 10 5 . 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).
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 hydro- Figure 7: Volume-averaged magnetic energy scaled by ρ 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). magnetic jet that removes material along the poles of the disk and aids the transport of angular momentum (Maki & Susa 2004;Tan & Blackman 2004;Silk & 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 & 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. However, it 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 that decrease the strength of the magnetic field.

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 & Vasiliev (2004) and Vasiliev & Shchekinov (2006) found that cosmic rays pre-ionize the gas, which in turn facilitates the formation of H 2 and HD. This may decrease the minimum halo mass required to cools the gas by an order of magnitude at z ∼ 20. Next to the enhanced cooling, cosmic rays also directly heat the gas. In the calculations of Stacy & Bromm (2007), the additional cooling provided by the enhanced H 2 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 H 2 , a larger LW flux was necessary to suppress cooling and star formation prior to the onset of Ly-α cooling.

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 & Hirata 2010). This relative velocity has also been termed the 'streaming velocity'. At z 1000, 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 sound speed drops to that of a gas with a temperature of ∼ 10 4 K, such that the streaming velocity becomes supersonic. A 1σ-peak at z 1000 has a relative velocity of v rel 30 km s −1 and a Mach number of M 5. Following recombination, the relative velocity decreases ∝ a −1 , such that at z 20 the streaming velocity is close to 1 km 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 5 − 100 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 & McQuinn 2012). They found that streaming velocities reduce the baryon overdensity in minihalos and possibly delay the onset of cooling (see . 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 & 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 21cm signal Visbal et al. 2012;Fialkov et al. 2012Fialkov et al. , 2013Fialkov et al. , 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.

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 ∼ 100 GeV. These are expected to have a very small cross section for selfannihilation of σv 3 × 10 −26 cm 3 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 electronpositron 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 is boosted to a level where it may become 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. Enclosed Star Forming Mass (M ) 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. 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 10−100 Mpc scales. Adapted from O' Leary & McQuinn (2012).
the H 2 line cooling rate at a density of n H 10 13 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 & 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 ∼ 10 6 M . With a luminosity of up to 10 10 L in the infrared, they may be observable by the James Webb Space Telescope ( 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 for the standard WIMP scenario this is not the case. 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 H 2 (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.

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 1 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 WDM = 10 keV on the formation of the first stars. In this case, the matter power spectrum has an exponential cut-off at 0.05 Mpc, which corresponds to a mass of 5 × 10 6 M . As a result, they found that star formation in minihalos is nearly entirely suppressed, and must await the virialization of larger halos. O' Shea & Norman (2006) used somewhat more detailed simulations to investigate the influence of WDM particles with masses in the range 10 keV m WDM 50 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 H 10 5 cm −3 , the thermodynamic evolution of the cloud became nearly indistinguishable from that in the CDM paradigm. Gao & 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 & Umemura 1999Bessho & Tsuribe 2012. Since the particle noise exceeds the WDM fluctuation power on small scales, these simulations could not derive the resulting fragment mass. Maio & 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 H 2 molecules (Biermann & 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 NL , Maio & Iannuzzi (2011) and Maio & 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

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 & 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 10 4 K. In these halos, cooling and star formation 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 10 4 K. The corresponding relation between virial mass and virial temperature may be obtained via equation 2: The characteristic virial mass of an atomic cooling halo is therefore 10 7 M vir 10 8 M , with a typical formation redshift of z 10-15 for 2-3σ peaks.

Turbulence
The virialization of primordial gas in atomic cooling halos without prior star formation in minihalos has been investigated by Wise & 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, H 2 line cooling takes over and the turbulence becomes supersonic with Mach numbers M 5. 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 & 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 & 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 & Birnboim 2006). In the simulations of Wise & 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).

Radiative feedback
Analytic considerations as well as simulations have indicated that the first stars had masses M * ∼ 100 M , 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 & Siess 2001;Siess et al. 2002;Gil-Pons et al. 2005Suda 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 & Shull 2000;Bromm et al. 2001b;Schaerer 2002Schaerer , 2003. These studies concluded that massive Population III stars radiate approximately as blackbodies with an effective temperature of 10 5 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 11.2-13.6 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).

Photodissociating radiation
One possible way to prevent the formation of H 2 is to photodetach H − , which is one of the main reaction partners in forming H 2 . 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 & Haiman 2012). Furthermore, the direct dissociation of H 2 via radiative excitation to the vibrational continuum is highly forbidden. Most of the H 2 is therefore dissociated by the two-step Solomon process (Field et al. 1966;Stecher & Williams 1967). Radiation in the LW bands excites a higher electronic state of H 2 , 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 diss 10 8 F LW s −1 , where F LW is the average flux density in the LW bands in erg s −1 cm −2 Hz −1 (Draine & 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 & Nishi 1999;Nishi & Tashiro 2000;Glover & Brand 2001). Subsequent work concentrated on the effects of LW radiation on cosmological scales, finding that the H 2 in the IGM with a fractional abundance of only y H 2 ∼ 10 −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 10 4 K well before reionization (Haiman et al. 1997;Kepner et al. 1997;Ciardi et al. 2000a;Haiman et al. 2000;Kitayama et al. 2001). However, only a modest intensity of the order of J LW = 10 −23 erg s −1 cm −2 Hz −1 sr −1 is built up by z 20, which is not sufficient to quench star formation in minihalos Kitayama et al. 2001;Ricotti et al. 2002;Yoshida et al. 2003c;Greif & Bromm 2006;MacIntyre et al. 2006;Johnson et al. 2008;Trenti & Stiavelli 2009). Instead, the minimum virial mass required for efficient H 2 cooling increases relatively slowly (Machacek et al. 2001;Yoshida et al. 2003c;Mesinger et al. 2006;Johnson et al. 2007;Wise & Abel 2007b;O'Shea & Norman 2008;Latif et al. 2014d;Visbal et al. 2014c). Unless the LW flux is extremely high, H 2 cooling may even become important in atomic cooling halos Omukai 2001;Omukai & Yoshii 2003;Oh & 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 & 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 H 2 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 & Shapiro 2007;Susa & Umemura 2006;Susa 2007;Susa et al. 2009;Whalen et al. 2008aWhalen et al. , 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 H 2 formation in the thin shell ahead of the H ii region due to the dissociation of H − and H + 2 .

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 30 km s −1 that begins to blow out the gas. After 10 5 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 threedimensional simulations qualitatively confirmed these results, and found that 100 M 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 × 10 4 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 10 4 K. The elevated electron fraction facilitates the formation of H 2 to a level of y H 2 10 −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 et al. 2007;Yoshida et al. 2007a,b). 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 et al. , 2001, 2004Susa et al. 2009;O'Shea et al. 2005;Mesinger et al. 2006Mesinger et al. , 2009Whalen et al. 2008aWhalen et al. , 2010Hasegawa et al. 2009). Whalen & 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.

X-rays
Although the direct emission of X-rays from Population III stars is likely cosmologically unimportant (e.g. Venkatesan & 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 ≥ 1 keV (Haiman & Loeb 1999;Venkatesan et al. 2001;Glover & Brand 2003;Machacek et al. 2003;Madau et al. 2004;Ricotti & Ostriker 2004b;Salvaterra et al. 2005;Tanaka et al. 2012). They may also increase the intergalactic H 2 abundance by more than an order of magnitude and reduce the clumping of the IGM (e.g. Kuhlen & Madau 2005). Similar to the evolution of relic H ii regions, the enhanced H 2 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 & Bromm 2006;Alvarez et al. 2009;Milosavljević et al. 2009a,b;Park & Ricotti 2011Aykutalp et al. 2013). It is therefore unlikely that miniquasars will grow fast enough to explain the presence of super-massive BHs at z 6 (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;Machida 2008;Machida et al. 2009b;Turk et al. 2009;Stacy et al. 2010;Greif et al. 2013;Stacy & 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;Jeon et al. 2012;Power et al. 2013;Jeon et al. 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 & Brand 2003).

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 140 − 260 M , a so-called pair-instability SN disrupts the entire star Heger & Woosley 2002;Heger et al. 2003;Joggerst & Whalen 2011;Baranov et al. 2013;Chen et al. 2014d,b). 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 10 53 ergs. In the range 100-140 M , pulsational instabilities drive episodic outbursts, while in the range 40-100 M the entire star collapse to a BH. Below 40 M , a SN with ∼ 10 51 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 & Limongi 2002, 2004Umeda & Nomoto 2002Iwamoto et al. 2005;Tominaga et al. 2007b;Zhang et al. 2008;Joggerst et al. 2009;Heger & Woosley 2010;Joggerst et al. 2010b;Limongi & Chieffi 2012). For masses 260 M , 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 10 5 -10 6 M 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 10 4 -10 5 M (Heger et al. 2003;Begelman et al. 2006Begelman et al. , 2008Begelman 2010;Volonteri & Begelman 2010;Montero et al. 2012;Hosokawa et al. 2012a;Volonteri 2012;Inayoshi et al. 2013;Hosokawa 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 10 55 ergs of kinetic energy (Whalen et al. 2013d,c,h). The various fates of Population III stars are illustrated in Figure 9.
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 & Maeder 2002b,a;Meynet et al. 2006;Heger et al. 2005;Chiappini et al. 2006;Hirschi 2007;Chiappini et al. 2008;Ekström et al. 2008;Takahashi et al. 2014), and their final fates (Suwa et al. 2007b;Joggerst et al. 2010a; Chat- Above 30 M , 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 & Woosley (2002).

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 & 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 Greif et al. 2007;Wise & 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 H 2 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 10 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 Sakuma & 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 supermassive stars with up to 10 5 M in atomic cooling halos that remained metalfree (Johnson et al. 2013b;Whalen et al. 2013d,c). 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 10 55 ergs, these SNe were able to completely disrupt their host halos.

Chemical Enrichment
In addition to their mechanical feedback, SN explosions from the first stars enrich the IGM with metals Scannapieco et al. 2001;Mori et al. 2002;Scannapieco et al. 2002;Mackey et al. 2003;Scannapieco et al. 2003;Wada & Venkatesan 2003;Ricotti & Ostriker 2004a;. 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 ∼ 10 −3 Z (Wise & 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-100 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 ∼ 20 (Greif & Bromm 2006;Trenti & 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 Tornatore et al. 2007;Ricotti et al. 2008;Maio et al. 2010Maio et al. , 2011aMuratov et al. 2013a,b). 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 & Truran 2003;Fang & Cen 2004;Ricotti & Ostriker 2004a;Yoshida et al. 2004;Matteucci & Calura 2005;Greif & Bromm 2006). Despite the progress made, the degree to which the metals mix with primordial gas is not yet fully understood , and depends on the employed sub-grid model (Greif et al. 2009a;Ritter et al. 2012Ritter et al. , 2014.

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 ∼ 10 −3 Z ,  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 10 6 and 10 9 M , and gradually enriches the IGM with metals. Adapted from Wise et al. (2012).
which reduces the characteristic Jeans mass of the cloud and facilitates fragmentation (see also Bromm & Loeb 2003b;Santoro & Shull 2006;Safranek-Shrader et al. 2010). Later simulations included H 2 cooling and confirmed that metal line cooling reduces the fragment mass (see Figure 11; Smith & Sigurdsson 2007;Smith et al. 2008Smith et al. , 2009Safranek-Shrader et al. 2014b).
Other studies argued that H 2 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(Jappsen et al. , 2009a. 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 H 2 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 10 −3 Z cooled to the temperature of the CMB in the presence of a strong UV background that dissociated the H 2 .
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 ∼ 10 −6 Z . At high redshifts, dust is thought to be produced in the SN remnants of Population III stars (Todini & Ferrara 2001;Nozawa et al. 2003;Salvaterra et al. 2004;Nozawa et al. 2006Nozawa et al. , 2007Gall et al. 2011a,b;Nozawa et al. 2012Nozawa et al. , 2014. A characteristic mass of ∼ 1 M 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;Omukai et al. 2005;Schneider et al. 2006;Tsuribe & Omukai 2008;Omukai et al. 2010;Schneider & Omukai 2010;Schneider et al. 2012;Chiaki et al. 2013aChiaki et al. , 2014. Tsuribe & Omukai (2006) modeled dust cooling in threedimensional simulations using idealized initial conditions and found that the central clump becomes elongated and fragments for Z 10 −6 Z . 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 10 −5 Z . 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 Figure 11: Thermodynamic evolution of metal-enriched gas in an atomic cooling halo at z 16. 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). 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 H 2 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 0.1 M .
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 & 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.

Direct collapse black holes
If the local LW flux is high enough, the formation of H 2 in minihalos may be suppressed until the halo has become massive enough for Ly-α cooling to become important (Omukai 2001;Bromm & Loeb 2003a;Volonteri & Rees 2005;Spaans & Silk 2006;Schleicher et al. 2010b;Johnson et al. 2013c). The gas may then contract isothermally at a temperature of 10 4 K to n H 10 6 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 H 10 16 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 Omukai et al. 2008).
Initial calculations showed that the critical LW flux required to suppress H 2 cooling for a photospheric temperature of 10 5 K corresponding to Pop-ulation III stars is of the order of 10 3 in units of J 21 , where J 21 = 10 −21 erg s −1 cm −2 Hz −1 sr −1 (Omukai 2001). A flux of this level dissociates H 2 up to a density of n H 10 4 cm −3 , where the level populations of H 2 transition to LTE and H 2 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 H 2 is used. The inclusion of radiation from Population I/II stars and the treatment of H − dissociation further reduces the critical flux (Wolcott-Green & Haiman 2012). Shang et al. (2010) find a value of only J 21,crit 30 − 300 for a photospheric temperature of 10 4 K corresponding to Population I/II stars. On the other hand, Latif et al. (2014b) employ higher resolution, a more accurate H 2 self-shielding formula, and investigate halos that virialize at slightly higher redshifts. They find a critical flux of J 21,crit ∼ 10 3 .  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,crit ∼ 10 3 , while Van Borm & Spaans (2013) argue that magnetic fields and turbulence potentially reduce J 21,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 10 5 than 10 4 K (Sugimura et al. 2014). Latif et al. (2014a) accounted for this fact and found that hydrodynamical effects such as shockheating and inhomogeneous collapse, which are only captured in realistic three-dimensional simulations, increase the critical LW flux to a few times 10 4 instead of 10 3 . 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 & Bromm 2006). However, the local LW flux can be boosted by many orders of magnitude near star formation sites (Dijkstra et al. 2008;Agarwal et al. 2012;Dijkstra et al. 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,crit = 30 and found ∼ 10 −3 direct collapse black holes per comoving Mpc at z = 10, while Dijkstra et al. (2014) used J 21,crit = 300 and included metal enrichment, finding a number density of 10 −10 − 10 −5 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,crit = 10 4 (e.g. Latif et al. 2014a;Sugimura et al. 2014;Yue et al. 2014).
Other mechanisms for suppressing H 2 formation and cooling have been suggested as well. Inayoshi & Omukai (2012) proposed that cold accretion flows may penetrate deep into atomic cooling halos and shock-heat the gas at densities n H 10 4 cm −3 , where H 2 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 H 2 cooling becomes unimportant. Similarly, Tanaka & 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 H 2 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.
If H 2 cooling is indeed suppressed until the halo reaches a virial temperature of 10 4 K, the gas contracts nearly isothermally and becomes gravitationally unstable at a Jeans mass of 10 5 − 10 6 M . 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 & Haiman 2002;Koushiappas et al. 2004;Begelman et al. 2006;Lodato & Natarajan 2006;Wise et al. 2008;Begelman & 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 M 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 .
The evolution of the central gas beyond the initial collapse of the gas was investigated in a number of studies. Regan & Haehnelt (2009) used adaptive mesh refinement simulations with a maximum resolution of 0.01 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 H 10 9 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 H 10 6 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 & Haehnelt (2009). They found that even though fragmentation occurs, the central object continues to grow via turbulent accretion and mergers at a rate of 1 M yr −1 . At this rate, a super-massive star with 10 6 M would form after only 1 Myr (see also Inayoshi & Haiman 2014). Similar fragmentation and accretion was found in the high-resolution simulations of Regan et al. (2014a). In a complementary approach,  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 H 10 17 cm −3 , they find a minimum Jeans mass of 0.2 M , 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 , 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 1 M 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 15 M . 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: 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, 1000 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 2010). There may also be a connection between the first stars and globular clusters (Padoan et al. 1997;Bromm & Clarke 2002;Beasley et al. 2003;West et al. 2004;Kravtsov & 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 Gnedin & Kravtsov 2006;Moore et al. 2006;Read et al. 2006;Ricotti et al. 2008;Salvadori et al. 2008;Bovill & Ricotti 2009;Muñoz et al. 2009;Ricotti 2009;Salvadori & Ferrara 2009;Frebel et al. 2010;Bovill & Ricotti 2011a,b;Frebel & Bromm 2012;Karlsson et al. 2012;Milosavljević & Bromm 2014). There is even a possibility of finding true Population III stars (Greif et al. 2011a;, 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;Fujimoto et al. 2000;Weiss et al. 2000;Schlattl et al. 2002;Picardi et al. 2004;Suda et al. 2004;Weiss et al. 2004;Lucatello et al. 2005;Lau et al. 2007;Tumlinson 2007a;Lau et al. 2009;Suda & Fujimoto 2010;Starkenburg et al. 2014). Low-mass Population III stars may also accrete metal-enriched gas from the IGM as they move through the Galaxy Frebel et al. 2009;.

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 threedimensional simulations of primordial star formation that started from cosmological initial conditions and included primordial chemistry networks (e.g. Abel et al. 2002;. They established the 'standard model' of primordial star formation, in which the physics of H 2 cooling leads to a characteristic stellar mass of ∼ 100 M . 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. , 2008. A key insight gained from these simulations is that primordial gas clouds are prone to fragmentation, but that the sec-ondary 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 & 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 & 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 & Abel 2008b;Greif et al. 2010). The UV radiation of massive Population III stars dissociates molecules and photoheats the gas to 10 4 K (e.g. Alvarez et al. 2006;Johnson et al. 2007;O'Shea & 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 & Omukai 2006;Dopcke et al. 2013;Safranek-Shrader et al. 2014b). The first galaxies started the process of reionization (e.g. Wyithe & 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 & Loeb 2003a;Dijkstra et al. 2008;Latif et al. 2013a;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. Toma, K., Sakamoto, T., & Mészáros, P. 2011, ApJ, 731, 127 Tominaga, N. 2009 Tominaga, N., Maeda, K., Umeda, H., et al. 2007a, ApJ, 657, L77 Tominaga, N., Umeda, H., & Nomoto, K. 2007b, ApJ, 660, 516 Toomre, A. 1964, ApJ, 139, 1217