EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 572 (December 2014) A&A, 572 (2014) A22 Full HTML
Free Access
Issue
A&A
Volume 572, December 2014
Article Number A22
Number of page(s) 18
Section Cosmology (including clusters of galaxies
DOI https://doi.org/10.1051/0004-6361/201424658
Published online 20 November 2014

© ESO, 2014

1. Introduction

Several very bright quasars have been detected at z> 6, which suggests that supermassive black holes (SMBHs) with masses of ~109 M⊙ already existed when the Universe was less than 1 Gyr old (Fan 2006; Willott et al. 2010; Mortlock et al. 2011; Venemans et al. 2013). It is challenging to explain how such SMBHs could have assembled so soon after the Big Bang, in particular how and when the “seeds” of these SMBHs formed and how their subsequent growth proceeded. Various scenarios for the formation of seed black holes in the early Universe have been proposed and are briefly discussed below (for a detailed review, see Volonteri 2010; Haiman 2013).

Perhaps the most obvious scenario assumes that SMBHs grow from the first stellar remnants. The first stars are thought to form at redshifts of ~20−50 in minihaloes of ~106 M⊙, cooled by molecular hydrogen (Tegmark et al. 1997; Abel et al. 2002; Bromm et al. 2002; Latif et al. 2013c; Bovino et al. 2013; Hirano et al. 2014). The first generation of stars was expected to have a more top-heavy initial mass function (typical stellar masses ~100 M⊙) than the current mode of star formation, resulting from the inefficient cooling in these minihaloes. However, more recent simulations that follow the collapse beyond the formation of the first core find that fragmentation may be effective after all, and thus the first stars may tend to form in clusters with much lower individual masses than initially expected, ≲10 M⊙ (Turk et al. 2009; Stacy et al. 2010, 2012; Clark et al. 2011b; Greif et al. 2011; Bovino et al. 2014b; Susa et al. 2014). Accretion luminosity does not seem to have much influence on the fragmentation behaviour (Smith et al. 2011, 2012). On the other hand, stellar UV feedback appears to inhibit accretion onto the protostar, which would result in an upper limit on the stellar mass of ~50−100 M⊙ (Hosokawa et al. 2011, 2012b; Susa 2013). Even if very massive stars were able to form and collapse into seed black holes, it would be difficult for them to accrete sufficient mass in the available time. It has been suggested that super-Eddington accretion may be necessary to accomplish this (e.g. Madau et al. 2014). However, the HII region formed around the seed black hole significantly reduces the accretion rate onto the seed (Milosavljević et al. 2009b,a; Alvarez et al. 2009; Johnson et al. 2011). In addition, it has been found that stars with masses ~100−10 000 M⊙ may form in massive primordial haloes irradiated by a moderate UV background, with a strong correlation between the strength of the UV flux and the mass of a protostar (Latif et al. 2014c).

Another scenario predicts the formation of seed black holes from very compact nuclear star clusters, which may form at redshifts of ~10−15, after some metal enrichment has occurred so that metal line-cooling becomes effective, and in the presence of trace amounts of dust (Schneider et al. 2003; Portegies Zwart et al. 2004; Schneider et al. 2012a,b; Omukai et al. 2005; Omukai 2012; Clark et al. 2008; Klessen et al. 2012; Dopcke et al. 2011, 2013; Safranek-Shrader et al. 2014; Bovino et al. 2014a). In such a cluster, stellar collisions can occur in a runaway fashion and lead to the formation of a very massive star, finally resulting in a seed black hole with a mass up to ~3000 M⊙ (Begelman & Rees 1978; Devecchi & Volonteri 2009; Devecchi et al. 2010, 2012; Lupi et al. 2014). Although this is significantly more massive than what is expected for the first generation of stars, it will still be difficult for such seeds to grow into the observed SMBHs in the available time.

In this work, we focus on a third pathway: the so-called direct collapse scenario. In this case, the primordial gas in a halo would collapse directly into a single central object, without fragmenting (e.g. Bromm & Loeb 2003; Koushiappas et al. 2004; Begelman et al. 2006; Lodato & Natarajan 2006; Spaans & Silk 2006; Schleicher et al. 2010; Latif et al. 2013a; Regan et al. 2014). The most likely host candidates are haloes with virial temperatures ≳104 K at redshifts ~10−15. For a direct collapse to occur, it is important that fragmentation is suppressed, which is possible if the gas in the halo is kept hot (and thus the Jeans mass high). Hence, the formation of H2 must be inhibited so cooling can occur only through atomic hydrogen, because otherwise molecular hydrogen cooling will lower temperatures to ~200 K and fragmentation may occur. In the absence of H2 cooling, self-gravitating gas will collapse nearly isothermally until it becomes optically thick and the adiabatic phase sets in.

One plausible mechanism for suppressing the formation of sufficient H2 is the presence of a UV radiation background. If massive, the first generation of stars (Pop III) is expected to have a stellar spectrum with a characteristic temperature of ~105 K (T5 spectrum; Tumlinson & Shull 2000; Bromm et al. 2001; Schaerer 2002), while the lower mass second generation of stars (Pop II) has a softer spectrum with several 104 K (T4 spectrum). These two spectral types have been used in several studies (e.g. Omukai 2001; Bromm & Loeb 2003; Shang et al. 2010). Lyman-Werner radiation (11.2−13.6) eV with an intensity above a certain threshold is able to photo-dissociate H2 and H− (important for the formation of H2) and keep their abundance very low. A T5 spectrum will mainly directly photo-dissociate H2, while a T4 spectrum will be better at photo-detaching H−. The critical intensity required to suppress H2 formation in massive haloes where direct gas collapse can occur has been estimated at J21crit≳102−103\hbox{$J_{21}^{\rm crit} \gtrsim 10^2{-}10^3$}, where J21 denotes the specific intensity just below the Lyman limit (13.6 eV), in units of 10-21 erg.cm-2.sr-1.s-1.Hz-1 (e.g. Omukai 2001; Bromm & Loeb 2003; Schleicher et al. 2010; Shang et al. 2010; Van Borm & Spaans 2013; Latif et al. 2014b,c). This is relatively high compared to the expected cosmic UV background at the relevant redshifts (J21bg~10\hbox{$J_{21}^{\mathrm{bg}} \sim 10$} at z ~ 10) (e.g. Ahn et al. 2009; Holzbauer & Furlanetto 2012; Dijkstra et al. 2014). However, the UV background distribution has a long bright-end tail, owing to the presence of close (about 10 kpc) luminous neighbours, which means that there is a small but significant subset of haloes that is exposed to supercritical intensities (Dijkstra et al. 2008, 2014; Agarwal et al. 2014; Visbal et al. 2014). Recently, though, it has been shown that it is important to consider spectra generated from realistic stellar populations, taking the mode of star formation (continuous or bursty) and the age, metallicity, and mass of the stars into account (Agarwal & Khochfar 2014; Sugimura et al. 2014). This has implications for the H2 photo-dissociation rate and the H− photo-detachment rate, and thus affects the value of the critical intensity, J21crit\hbox{$J_{21}^{\rm crit}$}. (Agarwal & Khochfar 2014) computed the reaction rate coefficients for H2 photo-dissociation and H− photo-detachment using realistic spectra resulting from a stellar synthesis code, and found that these depend on the age and metallicity of the stars, in contrast to the findings of (Sugimura et al. 2014). The latter used a one-zone model and realistic stellar spectra to also calculate J21crit\hbox{$J_{21}^{\rm crit}$}, finding values in the range between 1000 and 1400. Latif et al. (2014a) have studied the impact of varying the temperature of a blackbody spectrum in 3D cosmological simulations, to more closely resemble a realistic spectrum generated by Pop II stars. They found an even higher value for J21crit\hbox{$J_{21}^{\rm crit}$}, a few times 104, due to additional 3D effects. This value depends only weakly on the adopted radiation spectra in the range between Trad = 2 × 104 K and 105 K.

Alternative mechanisms for inhibiting H2 cooling comprise dissipation of a sufficiently strong magnetic field (Schleicher et al. 2009; Sethi et al. 2010; Van Borm & Spaans 2013) or the presence of strong shocks (Inayoshi & Omukai 2012), both of which result in collisional dissociation of H2.

Numerical 3D simulations have found fragmentation to be inhibited and thus show the feasibility of the direct collapse scenario. In some simulations, bar-like instabilities (Wise et al. 2008) or self-gravitating disks on parsec scales (Regan & Haehnelt 2009a) were found, though these employed a resolution of 16 cells per Jeans length. More recently, it was demonstrated that at least 32 cells, and preferably more, are required to properly resolve turbulence (Sur et al. 2010; Federrath et al. 2011; Turk et al. 2012; Latif et al. 2013b). New simulations employing a higher resolution find that it is likely that ~105 M⊙ objects will form (Latif et al. 2013a,b,e), though the peak density in these studies is not much higher than 1015 cm-3. In the simulations pursued here, we aim to complement these studies exploring collapse and fragmentation on smaller scales.

While simulating the formation of the protostar in 3D was not yet possible, various one-zone models employing detailed chemical models show the expected thermal pathway (Omukai 2000, 2001; Omukai et al. 2005, 2008). For a strong UV background, Omukai (2001) showed that clouds collapse nearly isothermally, cooled successively by Lyman-alpha emission of atomic hydrogen, two-photon emission of atomic hydrogen from the 2s state, and H− free-bound emission. Afterwards, the adiabatic phase sets in at ~1020 cm-3, at which point the minimum Jeans mass, and thus the characteristic mass of the protostar, has been reduced to 0.03 M⊙.

Once the protostar has formed, it will accrete and evolve into either a supermassive star or a quasi-star, depending on the accretion rate. The work by Schleicher et al. (2013) suggests that for accretion rates >0.14 M⊙/yr, a quasi-star will be the result, while lower accretion rates lead to the formation of a supermassive star. A supermassive star (SMS, with a mass in the range 103−106 M⊙) of fixed mass, supported by radiation pressure, is thought to evolve as an n = 3 polytrope and finally collapse into a black hole containing most of the stellar mass (Johnson et al. 2011; Whalen et al. 2013; Hosokawa et al. 2012, 2013). However, if the mass accretion rate is high enough, the outer layers of the SMS cannot thermally relax. In this case, it is not well-described by an n = 3 polytrope, but will have a more complex structure with a convective core surrounded by a convectively stable envelope that contains most of the mass. The core will burn up its hydrogen, and subsequently collapse into a black hole with a mass of a few M⊙. The resulting structure, where the black hole accretes material from the massive, radiation-pressure-supported envelope, is termed a “quasi-star” (Begelman et al. 2006, 2008; Begelman 2010; Volonteri & Begelman 2010; Ball et al. 2011, 2012).

As is known from present-day star formation, turbulence plays an important role in angular momentum transport and determining the fragmentation properties of collapsing gas clouds, since it can both locally compress the gas as well as provide additional support against collapse on larger scales (e.g. Larson 1981; Mac Low & Klessen 2004; McKee & Ostriker 2007; Federrath & Klessen 2012). Similar effects have been found at high redshifts in simulations of minihaloes, where turbulence plays a role in distributing angular momentum (Abel et al. 2002), and affects the fragmentation behaviour (Clark et al. 2011a; Turk et al. 2012; Latif et al. 2013c). Also in simulations of more massive, atomic cooling haloes, the importance of turbulence has been recognized (Greif et al. 2008; Wise et al. 2008). However, many of these older studies do not employ a sufficient Jeans resolution, as its impact was only recognized later. Latif et al. (2013b) found that the amount of turbulent structure increases significantly with increasing resolution, and in the study by Latif et al. (2013a) it was found that fragmentation occurs occasionally, but that this does not prevent the growth of a central massive object resulting from turbulent accretion and mergers.

Numerical simulations of collapsing gas in minihaloes show that fragmentation also depends on the amount of rotation, with stronger rotation inducing more fragmentation (Bromm et al. 2002; Machida 2008; Hocuk & Spaans 2010). The study by Clark et al. (2008) shows that massive disk-like structures are assembled, fragmenting to form protostars. In atomic cooling haloes the effects of rotation have not yet been studied in detail, though Bromm & Loeb (2003) found that a single black hole is formed in low-spin galaxies, while higher spin galaxies tend to form binary black holes. In their simulations of atomic cooling haloes, Regan & Haehnelt (2009b) observed the formation of massive compact self-gravitating disks, and found mild fragmentation in one of the three simulated haloes.

In this paper we present the first study in which the formation of a massive protostar is simulated in 3D up to unprecedented high central densities (1021 cm-3), owing to improved modelling of the chemistry. A high spatial resolution is obtained as well; starting from pc scales, we are able to resolve scales down to a few solar radii. In addition, we investigate how the fragmentation behaviour of collapsing primordial gas in the presence of a strong Lyman-Werner radiation background is affected by varying amounts of turbulence and rotation. For each case the formation of clumps and their accretion rates are studied.

In Sect. 2 some details are given on the methods and setup of the numerical simulations that have been performed. In Sect. 3 the results for both the one-zone calculations and the 3D simulations are presented and discussed, and we conclude with a summary of the results in Sect. 4.

2. Numerical methodology and simulation setup

ENZO is an open-source adaptive mesh refinement (AMR) simulation code, which provides high spatial and temporal resolution for the modelling of astrophysical fluid flows (Bryan et al. 2014). It contains a wide variety of physics modules, making it suitable for many different astrophysical applications. We use a modified version of ENZO 2.3, replacing the chemistry implementation by a customized build of the KROME chemistry package (Grassi et al. 2014), as discussed in the following subsections. The hydrodynamical equations are solved using the MUSCL scheme, which is a second-order accurate extension of Godunov’s method. The implementation in Enzo uses second-order Runge-Kutta time integration, and the Riemann solver employed is the HLLC solver (Harten-Lax-van Leer with Contact), with a fallback to the more diffusive HLL solver (Harten-Lax-van Leer) in case negative energies or densities are computed. The choice of this solver is due to the strong shocks which occur once the central core becomes adiabatic and the central protostars forms. Self-gravity is computed by solving the Poisson equation using a multigrid technique.

2.1. Initial conditions

We follow the gravitational collapse of an isolated spherical gas cloud with a radius of 15 pc and a top-hat density profile, embedded in a 100 pc simulation box. The Jeans length is resolved by at least 64 cells at all times. Additionally, a refinement criterion based on overdensity is used. These combined criteria result in the simulations using 29 refinement levels, at which point an adiabatic core is formed and no further refinement is necessary. The collapse is followed for another 1.67 × 10-2 years, corresponding to ~5 freefall times, after the highest refinement level is reached. To ensure pressure equilibrium between the sphere and its surroundings, we set the initial sphere density to 1000 cm-3 and its temperature to 500 K, while the surrounding gas is initialized with a density of 100 cm-3 and a temperature of 5000 K. The above combination of parameters also ensures that the mass of the cloud (~3.5 × 105 M⊙) is greater than the local Jeans mass (~3 × 104 M⊙), and thus the cloud will collapse. The total mass contained in the box is ~2.8 × 106 M⊙. This setup has been chosen in order to be able to explore the formation of protostars up to very high densities. The cloud is irradiated by a UV background with a T5 spectrum (see Sect. 2.2.7) of intensity 105 in units of J21, so that the abundance of H2 is kept low and cooling will occur mainly through atomic hydrogen.

Furthermore, we add a certain amount of initial turbulence to the gas, as well as some rotation of the cloud. These parameters are varied to study and quantify their effects on the collapse dynamics and fragmentation properties. An overview of the different simulations can be found in Table 1. The turbulent velocities are sampled from a Maxwellian distribution with a temperature equal to the initial temperature of the gas sphere, and subsequently multiplied by the percentage given in the table. Since the maximum of the Maxwell distribution function is of the order of the sound speed cs, the turbulent velocities are of the order of a given percentage of cs. The amount of rotation is given in percentage of the Keplerian velocity; i.e. 100% rotation means the cloud is rotationally supported.

Table 1

Overview of the different simulations and their initial turbulent and rotational velocities.

2.2. Chemistry, heating, and cooling

We employ the KROME1 chemistry package, which has been developed in order to simplify the embedding of the chemistry and the microphysics in numerical simulations. It builds the corresponding rate equations, the solver parameters, and includes a series of thermal processes which are coupled to the chemical evolution. A patch to embed KROME in ENZO is available with the package and has been used within this work. KROME solves the non-equilibrium chemistry together with the temperature equation using the adaptive high-order solver DLSODES, which was shown to be both accurate and efficient for networks that present a corresponding ordinary differential equation system with a sparse Jacobian, and that are typical in astrophysical applications (Bovino et al. 2013; Grassi et al. 2013). We have modified and extended the available package, mainly to obtain the desired cooling processes. The main improvements of our modified version of KROME are the addition of H− cooling, Rayleigh scattering, and a different evaluation of the critical density used for the chemical heating, following Glover & Abel (2008).

2.2.1. Chemical network

Our chemical network consists of 31 reactions, including 9 species: H, e, H+, H−, H2, H+2\hbox{$_2^+$}, He, He+, and He++. All the reactions and their associated rates can be found in Appendix A.

2.2.2. Molecular hydrogen cooling

The molecular hydrogen cooling rates were taken from Glover & Abel (2008), with an opacity correction from Ripamonti & Abel (2004), as implemented in KROME. However, we modified the opacity correction to use the molecular hydrogen density instead of the total density, rendering it usable for cases with a non-zero UV background. Hence, the H2 cooling rate is multiplied by a factor min[1,nH2(/(8×109cm-3))-0.45]\hbox{$\min{\left[1, \left(n_{\rm H2} / \left(8 \times 10^{9}~{\rm cm^{-3}}\right)\right)^{-0.45}\right]}$}, where nH2 is the H2 number density. Recent studies by Greif (2014) and Hartwig et al. (2014) calculate the escape fraction of cooling photons using a multi-line, multi-frequency ray-tracing scheme, and an approach based on the TreeCol algorithm, respectively. Greif (2014) find that the radially averaged escape fraction agrees well with the analytical fit from Ripamonti & Abel (2004), while the results of Hartwig et al. (2014) suggest that this fit underestimates the escape fraction after the initial stage of collapse. Presently, it has not yet been investigated which of these two methods yields the most accurate results. However, additional one-zone calculations suggest that even a significantly larger escape fraction does not influence our results, as the ineffectiveness of the cooling is mainly the result of the low H2 abundance. Of course, opacity effects would become more important when considering a case where H2 is the dominant coolant.

2.2.3. Collision-induced emission cooling

When a collision takes place between an H2 molecule and another H2 molecule, a He molecule, or a H atom, the interacting pair briefly acts as a “supermolecule” with a non-zero electric dipole, and there is a high probability of a photon being emitted. Collision-induced emission (CIE) may become important at high densities, depending on the gas temperature. We use the fit provided in KROME for the optically thin rate, but modified to ensure it is 0 if fH2< 0.1 and does not become important before fH2 ~ 0.5, where fH2 is the H2 mass fraction relative to H, as it is uncertain whether the fit is still valid for extremely dissociated media. The opacity correction at high densities has been adopted from Ripamonti & Abel (2004), τCIE=max[10-5,(n2×1016cm-3)2.8],\begin{eqnarray} \tau_{\rm CIE} = \max{\left[10^{-5}, \left(\frac{n}{2 \times 10^{16}~{\rm cm^{-3}}}\right)^{2.8}\right]}, \end{eqnarray}(1)where n is the total number density. The CIE cooling rate is then multiplied by min[1,(1−exp(−τCIE)) /τCIE], where (1−exp(−τ)) /τ is the usual spherical escape probability.

2.2.4. Atomic cooling

KROME employs the atomic cooling rates from Cen (1992). These include the collisional ionization of H, He, He+, and He(2s) by electrons, the recombination of H+, He+, and He++, the dielectronic recombination of He+, the collisional excitation of H (all n), He (n = 2,3,4 triplets), and He+ (n = 2), and bremsstrahlung for all ions. The main cooling channel relevant here is the collisional excitation of H. We have added an optical depth approximation for the Rayleigh scattering by H atoms, which will suppress this main channel, as τRl=σH,RlnHIλJ2,\begin{equation} \tau_{\rm Rl} = \sigma_{\rm H,Rl} n_{\rm HI} \frac{\lambda_{\rm J}}{2}, \end{equation}(2)where λJ is the Jeans length, nHI is the number density of atomic hydrogen, and σH,Rl=5.799×10-29λ-4+1.422×10-30λ-6+2.784×10-32λ-8cm2\begin{eqnarray} \sigma_{\rm H,Rl} &=& 5.799 \times 10^{-29} \lambda^{-4} + 1.422 \times 10^{-30} \lambda^{-6} \notag\\ &&+ 2.784 \times 10^{-32} \lambda^{-8}~{\rm cm}^2 \end{eqnarray}(3)is the Rayleigh scattering cross section of H for radiation with wavelength λ (in μm; Kurucz 1970). The cooling rate is then multiplied by exp(−τRl). Additionally, we have added two fudge factors to mimic optical depth effects and thus reduce cooling at high densities (n ≳ 1017 cm-3), in accordance with the findings of Omukai (2001). The first factor, f1, represents that the gas should be optically thick to atomic hydrogen line cooling around ~1017 cm-3, and H ionization becomes the main atomic cooling channel. The second factor, f2, ensures that the gas becomes almost completely optically thick to radiative cooling around ~1020 cm-3, so that afterwards the evolution is nearly adiabatic. The fudge factors are calculated as fi = min[1,(1−exp(−τfi)) /τfi] for i = 1,2, using the functional form of the spherical escape probability. The opacities τf1 and τf2 are given by2τf1=max[10-5,(n1017cm-3)5],τf2=max[10-5,(n1020cm-3)8]·\begin{eqnarray} \tau_{\rm f_1} &=& \max{\left[10^{-5}, \left(\frac{n}{10^{17}~{\rm cm^{-3}}}\right)^5\right]}, \\ \tau_{\rm f_2} &=& \max{\left[10^{-5}, \left(\frac{n}{10^{20}~{\rm cm^{-3}}}\right)^8\right]}\cdot \end{eqnarray}

2.2.5. H− cooling

Through radiative association of H and e, H− is formed and a photon is emitted. There will be net cooling if this photon can escape (Omukai 2001; Schleicher et al. 2008). The cooling rate can then be approximated as ΛH−≈kH−nHIneEγ,\begin{equation} \Lambda_{H^-} \approx k_{H^-} n_{\rm HI} n_{\rm e} E_{\gamma}, \end{equation}(6)where Eγ is the approximate energy of the emitted photon. A typical electron undergoing radiative attachment has an energy of the order of kBT, so the average outgoing photon energy can be estimated as Eγ ~ E0 + kBT, where the binding energy E0 of H− is 0.755 eV. Rayleigh scattering (see Sect. 2.2.4), as well as H− bound-free absorption, will suppress this cooling channel, so optical depth approximations for these processes have been taken into account. The cross section for H− bound-free absorption is (John 1988) σH−,bf=10-18λ3(1λ−1λ0)1.5f(λ),\begin{equation} \sigma_{\rm H^-,bf} = 10^{-18} \lambda^3 \left(\frac{1}{\lambda} - \frac{1}{\lambda_0}\right)^{1.5} f\left(\lambda\right), \end{equation}(7)where λ is the wavelength of the scattered radiation in μm, λ0 = 1.6419 μm, and fλ)(\hbox{$f\left(\lambda\right)$} is given by Eq. (5) in John (1988).

2.2.6. Chemical cooling and heating

Various chemical reactions can result in net cooling or heating of the gas (Omukai 2000). In our case, the most important ones are the three-body formation of H2 (Forrey 2013, see Bovino et al. 2014c for a comparison of different rates) and collisional dissociation of H2 (Shapiro & Kang 1987; Martin et al. 1996, 1998). The collisional dissociation process releases 4.48 eV per dissociated H2 molecule (its binding energy), cooling the gas, while the heat deposited by three-body formation is 4.48(1 + ncr/n)-1 eV per H2 molecule. Here, ncr is the critical density, calculated as (Glover & Abel 2008) ncr=(xHIncr,HI+xH2ncr,H2)-1,\begin{equation} n_{\rm cr} = \left(\frac{x_{\rm HI}}{n_{\rm cr,{\rm HI}}} + \frac{x_{{\rm H2}}}{n_{\rm cr,{\rm H2}}}\right)^{-1}, \end{equation}(8)where xHI and xH2 are the number fractions of HI and H2, respectively, and ncr,HI and ncr,H2 are their respective critical densities, given by ncr,HI=dex[3−0.416logT4−0.327(logT4)2],ncr,H2=dex[4.845−1.3logT4+1.62(logT4)2],\begin{eqnarray} n_{\rm cr,{HI}} &=& {\rm dex}{\left[3 - 0.416 \log{T_4} - 0.327 \left(\log{T_4}\right)^2\right]}, \\ n_{\rm cr,H_2} &=& {\rm dex}{\left[4.845 - 1.3 \log{T_4} + 1.62 \left(\log{T_4}\right)^2\right]}, \end{eqnarray}where T4=T104K·\hbox{$T_4 = \frac{T}{10^{4}~{\rm K}}\cdot$}

2.2.7. Radiation background

In our calculations we have used a constant UV background flux with a T5 spectrum below the Lyman limit, which will photo-dissociate H2 and photo-detach H−. The main difference with a T4 spectrum is that lower values of the intensity, J21, are required for the gas to collapse isothermally. We do not expect the choice of spectrum or the specific strength of the UV background to matter, as long as the H2 abundance is kept low so that H2 cooling is unimportant. The difference between the spectra is expressed in the photo-dissociation rate of H2 and photo-detachment rate of H− (see k24 and k25 in Appendix A). We also include H2 self-shielding, using the improved fit described in Wolcott-Green et al. (2011), fsh=0.965(1+xNH2/b5)1.1+0.035(1+xNH2)0.5  ×exp(−8.5×10-4(1+xNH2)0.5),\begin{eqnarray} f_{\rm sh} &=& \frac{0.965}{\left(1 + x_{N_{\rm H_2}}/b_5\right)^{1.1}} + \frac{0.035}{\left(1 + x_{N_{\rm H_2}}\right)^{0.5}} \notag\\ & &\quad\quad \times \exp{\left(-8.5 \times 10^{-4} \left(1+x_{N_{\rm H_2}}\right)^{0.5}\right)}, \end{eqnarray}(11)where xNH2 is given by xNH2=NH25×1014cm-2,\begin{equation} x_{N_{\rm H_2}} = \frac{N_{\rm H_2}}{5 \times 10^{14}~{\rm cm^{-2}}}, \end{equation}(12)with NH2 the column density in cm-2, calculated as NH2 = nH2λJ/ 2. The Doppler broadening parameter for H2, b5, is given by b5=10-5(2kBT2mH)0.5,\begin{equation} b_5 = 10^{-5} \left(\frac{2 k_{\rm B} T}{2 m_{\rm H}}\right)^{0.5}, \end{equation}(13)in units of 105 cm/s. The photo-dissociation rate of H2 is multiplied by the self-shielding factor fsh.

3. Results

3.1. One-zone calculations

thumbnail Fig. 1

Physical quantities as a function of number density in a one-zone calculation, irradiated by a strong T5 background, using our modification of KROME. A: temperature; B: self-shielding factor for H2, fsh (fsh = 1 means no shielding, and the smaller fsh, the stronger the shielding); C: number fractions of H species (and electrons); D: number fractions of He species.

We have performed a one-zone freefall collapse test, already included with the KROME package, to verify our chemical model. The results for a T5 background with J21 = 105 and initial conditions similar to those of the 3D simulations can be seen in Fig. 1, calculated up to a number density of 1021 cm-3. Panel A shows the temperature evolution, panel B shows the self-shielding factor for molecular hydrogen, panel C shows the number fractions of the different H species (and electrons), and panel D shows the number fractions of the different He species. We note here that at densities above 1017 cm-3 an equilibrium approximation could be adopted and might speed up the calculations. Nevertheless, we preferred to follow a complete non-equilibrium evolution.

Initially, the temperature increases adiabatically, due to strong compressional heating. Because the molecular hydrogen is strongly dissociated by the UV background, the gas cannot cool through H2 and instead cools via other processes. During the initial adiabatic phase, H− cooling is the dominant cooling process, but it is not efficient enough to counter the strong heating. When the temperature reaches ~8000 K, around ~105 cm-3, Lyα cooling starts to become dominant and the temperature slope flattens off, now evolving nearly isothermally, though still decreasing slowly. Both chemical cooling and H− cooling are also important during this phase. Around a number density of ~108 cm-3, both of these rates become higher than the atomic cooling. The H− cooling channel becomes strongly suppressed around ~1016 cm-3 as the cloud becomes optically thick to both Rayleigh scattering and H− bound-free absorption. Chemical cooling still maintains the near-isothermal evolution briefly, but then chemical heating cancels out the cooling and the temperature starts rising. Collisional ionization of H starts at ~1017 cm-3, resulting in a slowdown of the temperature rise up to ~1019 cm-3. From this point on, the cloud collapses adiabatically, and after sufficient contraction a protostar is expected to form in the centre.

During the whole collapse, the molecular hydrogen fraction never becomes larger than 10-3, and as a result H2 cooling is not important (except for densities between 104−105 cm-3). Starting from ~1010 cm-3, three-body formation increases the H2 abundance, peaking just before the rise in temperature at ~1017 cm-3, after which strong collisional dissociation drastically decreases the abundance again. At low densities the self-shielding is too weak to prevent H2 from being photo-dissociated (the smaller the factor, the stronger the shielding). At densities above 1010 cm-3 the gas starts to become well-shielded, however, due to the high temperature, collisional dissociation of H2 becomes effective. Additionally, H2 cooling starts to become optically thick at these densities.

3.2. 3D simulation results

thumbnail Fig. 2

Physical quantities, weighted by mass, spherically averaged and radially binned, as a function of radius at the peak density output for the different simulations. A) Number density; B) temperature; C) H2 number fraction; D) electron number fraction; E) turbulent velocity; F) radial velocity, plotted as − vrad; G) enclosed mass, and the Jeans mass for T40R10 (it is very similar for other runs); H) radial mass infall rate, calculated from the density and the radial velocity. The simulation details and abbreviations are listed in Table 1.

thumbnail Fig. 3

Density projections along the x-axis for all simulations, showing the integrated number density for various scales at the peak density output. The simulation details and abbreviations are listed in Table 1.

thumbnail Fig. 4

Density projections along the z-axis for all simulations, showing the integrated number density for various scales at the peak density output.The simulation details and abbreviations are listed in Table 1.

An overview of the different 3D simulations and their abbreviations are listed in Table 1. We have performed one simulation for each set of initial conditions, five in total. After reaching the highest refinement level of 29, the simulations evolved for another ~1.67 × 10-2 years, corresponding to ~5 freefall times, with the freefall time (~3.5 × 10-3 yr) calculated at the moment when the highest refinement level is first reached.

Figure 2 shows several spherically averaged, radially binned profiles of various quantities for all simulations, centred on the peak density location (hereafter referred to as the central clump). The data shown has been obtained at the end of each simulation, when a peak density of 1021 cm-3 was reached. From the density profile (shown in panel A) it can be seen that in general, the density increases with decreasing radius, so that overall the evolution of quantities with decreasing radius corresponds to an evolution with increasing density. Specifically, the density increases approximately as ∝r-2, as is typical of an isothermal collapse. Deviations from this behaviour are caused by local over- or underdensities, resulting from the turbulent nature of the gas. In the very centre of the cloud, inside of the radius corresponding to the minimum Jeans mass, the density profile flattens off, indicating the central clump. This clump can also be seen in the enclosed mass profile (shown in panel G), which steeply decreases inside 10-7 pc, due to enhanced pressure support.

The density profile in simulation T20R10 deviates somewhat from isothermal, with a peak in the density profile around ~4 × 10-5 pc. After close inspection of density projections at different scales (see Figs. 3 and 4, particularly at the 50 pc scale), this appears to be due to the presence of a second concentration of mass containing two additional clumps, which have not collapsed as far as the main clump. However, from comparison runs with the same initial conditions, though with a different random seed for the initialization of the turbulent velocity field, we have found that such additional clumps are only sometimes present for the T20R10 initial conditions. Additionally, a second clump is also sometimes found for the other initial conditions discussed, though always with a lower peak density than the main clump. Hence, this fragmentation is likely not related to the amount of initial turbulence or rotation. It is not yet clear whether these additional clumps will continue to collapse, or instead accrete onto the main clump. Based on a simplified “toy” model of fragmentation in the accretion disk around a protostar, Inayoshi & Haiman (2014) argue that some of the clumps formed in the disk may evolve to zero-age main sequence stars, but that most of these clumps can migrate inward and merge with the central protostar.

The temperature evolution of the gas cloud (shown in panel B of Fig. 2) is very similar in all simulations. In the final stage displayed in the plots, the outer layers of the cloud are at a temperature of ~8000 K. Further inwards, the temperature evolves nearly isothermally, though still gradually drops to about 4000 K, until reaching a radius of about 3 × 10-6 pc. Inside this radius, the evolution proceeds nearly adiabatically and the temperature reaches ~7 × 104 K by the time the peak density is reached. This behaviour is expected based on the one-zone calculations, and for a more detailed description of the involved processes, see Sect. 3.1.

The molecular hydrogen number fraction (shown in panel C) exhibits some variation between simulations at larger radii, but converges for small radii, with only the T40R0 case deviating slightly from the others. At large radii, corresponding to low densities, the fraction slowly decreases as H2 is dissociated by the UV background. Next, a steep increase occurs at the radius corresponding to a density of 1010 cm-3, due to 3-body formation becoming efficient, after which formation and dissociation approximately balance each other for a broad density range. For the highest densities, where the gas is heating up, collisional dissociation starts to dominate and the fraction drops drastically. The H2 number fraction never gets much larger than 10-3, in agreement with the one-zone test, which means that there is never enough H2 for molecular cooling to be important.

The overall evolution of the electron number fraction (shown in panel D) is again quite similar to what is expected from one-zone calculations. The T80R10 case deviates slightly from the others, in that the electron fraction starts to increase already at somewhat larger radii. This is again due to some mass buildup around that radius, reaching slightly higher temperatures than the surrounding matter. The T20R10 case deviates quite strongly around the radius where for the other simulations the minimum occurs, which is due to the aforementioned second mass concentration at that radius.

The RMS turbulent velocity (shown in the panel E) increases slowly with radius from ~1 km s-1 to ~10 km s-1 for all simulations. It is interesting that although initially the amount of turbulence is varied, later in the runs this difference is smoothed out and at least in the turbulent velocities there is no longer a clear difference between the high and low turbulence cases. The radial velocity (shown in panel F) is similar for all simulations as well, and stays around ~11 km s-1 throughout most of the cloud. Only for the smallest radii, inside the minimum Jeans radius, does the radial velocity decrease down to 1−0.1 km s-1, due to the pressure support in the clump.

The radial accretion rate (the rate of mass flow towards the central clump; shown in panel H) is calculated as dM/ dr = 4πr2ρvrad, where ρ is the density and vrad the radial velocity. The rate varies somewhat between different simulations, although there does not seem to be a trend with either turbulence or rotation. The large peak in the accretion rate for the T20R10 run around ~4 × 10-5 pc is due to the close connection of the second mass concentration to the central clump, locally boosting the accretion rate. It can be seen that both the density and radial velocity show a peak at the same location, causing the enhanced accretion rate. Similar features in the accretion rate were found by Regan & Haehnelt (2009a), who also attribute them to clumps of high-density gas. Overall accretion rates of a few solar masses per year are observed in all cases. Given such high accretion rates, a supermassive star of 105 M⊙ is expected to form within 105 years.

From the density projections (Figs. 3 and 4), it can be seen that for the simulations including rotation a disk has formed. Stronger rotation leads to a flatter, more extended disk, and more pronounced spiral structures. The differences in turbulent structures show on the 0.1 pc scale, with an increased amount of structure for higher initial turbulent velocities, also enhanced by stronger rotation. On smaller scales, the differences are no longer clear, as can also be seen from the turbulent velocity profiles in panel E of Fig. 2.

thumbnail Fig. 5

Density and temperature slices along the x-axis for all simulations and for two different scales at the peak density output. The simulation details and abbreviations are listed in Table 1.

thumbnail Fig. 6

Density and temperature slices along the z-axis for all simulations and for two different scales at the peak density output. The simulation details and abbreviations are listed in Table 1.

Figures 5 and 6 show temperature slices for two different scales, next to density slices of the same area. It can be seen that there are hot regions of gas surrounding slightly cooler patches. Such warmer and cooler patches result from local compression and expansion of the gas due to turbulent motions.

thumbnail Fig. 7

Physical quantities, weighted by mass, disk averaged and radially binned, as a function of radius at the peak density output for the different simulations. A) Toomre Q parameter; B) ratio of the rotational to the Keplerian velocity; C) sound speed; D) turbulent velocity; E) rotational frequency; F): surface density. The simulation details and abbreviations are listed in Table 1.

In Fig. 7 we explore the properties of the disk, by displaying several disk averaged, radially binned quantities (using the radius in the x − y plane) for all simulations except T40R0 (as there is no disk present), centred on the peak density location (hereafter referred to as the central clump). The data shown has been obtained at the end of each simulation, when a peak density of 1021 cm-3 was reached.

In panel A, the Toomre Q parameter is shown. This parameter is calculated as Q=σΩπGΣ\hbox{$Q = \frac{\sigma \Omega}{\pi G \Sigma}$}, where σ is the RMS of the sound speed and the turbulent velocity (as both thermal and turbulent motions will play a role in stabilizing the disk; plotted in panels C and D), Ω is the rotation frequency (plotted in panel E), G is the gravitational constant, and Σ is the surface density (plotted in panel F). The disk is stable when Q is larger than a critical value, which is of order unity; when Q approaches this threshold, the disk will become gravitationally unstable. It can be seen that the disk is mildly unstable for all simulations. Only for the smallest radii Q becomes decidedly larger than one, which is due to the close-to-adiabatic core in the central region. On the smallest scales, the adiabatic heating thus stabilizes the protostar against further collapse. We note that the adiabatic equation of state is however only an approximation, while real systems may evolve further via Kelvin-Helmholtz contraction.

The ratio of the rotational velocity to the Keplerian velocity (shown in panel B) follows roughly the same behaviour as the Toomre Q parameter. The ratio is more or less constant over most radii. Some imprint of the initial amount of rotation remains, as the T40R20 run has the highest ratio over nearly all radii. However, for all runs the ratio has increased compared to the initial value, as a spin-up occurs during collapse.

It is interesting to note that we do not find a clear trend with either turbulence or rotation in any of the measured quantities on the smaller scales. It appears that whatever the initial conditions are, at later stages the initial difference in turbulence and rotation is smoothed out on these scales. Of course, on larger scales the presence and size of a disk does vary according to the initial amount of rotation, and on intermediate scales there are more turbulent structures for an increasing amount of initial turbulence, but this does not affect the overall evolution of density, temperature, accretion rate, and other quantities on scales smaller than 1 pc.

Whether one or more clumps are present does not depend on the initial amount of turbulence or rotation either, as we have concluded from comparison runs with the same initial conditions and a different random seed, in which usually one, sometimes two, and in a single case three of these clumps form. However, we never find more than three clumps, none of which have collapsed as far as the main clump, meaning that there is not much fragmentation, irrespective of turbulence or rotation. As mentioned, the simulations evolved for another ~1.67 × 10-2 years after the highest refinement level was reached. Given that no fragmentation occurs in most of our simulations during this time, it can be considered as a lower limit on the fragmentation time scale.

3.2.1. The central object

A quantification of the properties of the central clump in each simulation is listed in Table 2. We find only one of such collapsed clumps in each simulation. The location of the “knee” in the enclosed mass profile is taken as the clump radius. The mass enclosed in this radius corresponds approximately to the minimum Jeans mass (see also panel G in Fig. 2), and thus the clumps are gravitationally bound. This object marks the onset of protostar formation. Due to computational constraints simulations we cannot evolve the simulations further, though we expect the gas in the surroundings to collapse to form a massive protostar.

Table 2

Properties of the central bound clump found in each simulation.

Given the radial accretion rates shown in Fig. 2, a supermassive star of 105 M⊙ is expected to form within 105 years, assuming that the gas reservoir to accrete from is large enough. If the accretion rate remains higher than 10-2 M⊙/ yr, Hosokawa et al. (2012) found that the star will have a bloated envelope and lower surface temperatures, which inhibits the emission of ionizing radiation. In this case, radiative feedback will not be able to interfere with the accretion process. If accretion rates higher than 0.14 M⊙/yr can be maintained until the core has exhausted its hydrogen content through nuclear burning (after ~7 × 106 yr), it is likely that the core of the star will collapse into a black hole, resulting in a quasi-star (Schleicher et al. 2013).

4. Discussion and conclusions

We have performed 3D adaptive mesh refinement simulations using the ENZO code, simulating the formation of a protostar up to unprecedented high central densities of 1021 cm-3, and spatial scales of a few solar radii. To achieve this goal, we have employed the KROME package to improve the modelling of the chemical and thermal processes. Particularly, we have investigated how the fragmentation behaviour of collapsing primordial gas in the presence of a strong Lyman-Werner radiation background is influenced by varying amounts of turbulence and rotation.

We found that in the runs including rotation, a mildly unstable disk forms on scales of ~0.5 pc, with a more extended disk for the stronger rotating case, run T40R20. On somewhat smaller scales, ~0.1 pc, the amount of turbulent structures increases with increasing initial turbulent velocities, as one would expect. However, on even smaller scales, ≲0.01 pc, the differences between the runs disappear, and radial profiles of the density, temperature, accretion rate, and other quantities are all very similar, with no dependence on the initial amount of turbulence or rotation. The thermal evolution of all runs is consistent with the one-zone result from Omukai (2001).

In each simulation we have found a single bound clump collapsed up to a density of 1021 cm-3, with a radius of 2 × 10-2 AU and a mass of ~7 × 10-2 M⊙, corresponding to the minimum Jeans mass. This clump marks the onset of protostar formation. Given the observed accretion rates of a few solar masses per year, the protostar is expected to become a quasi-star with a mass of 105 M⊙ within 105 years, assuming a high accretion rate can be maintained. Ferrara et al. (2014) have derived a detailed prediction for the initial mass function (IMF) of the first massive black holes formed in atomic cooling haloes, combining the physics of SMS evolution and direct collapse black hole formation and growth with cosmological merger-tree simulations. They have found that in the case that minihaloes can form stars and pollute the gas, the IMF is bimodal and spans a broad mass range, M ≈ (0.5−20) × 105 M⊙; while in the case that they cannot form stars, the IMF spans a narrower range, M ≈ (1−2.8) × 106 M⊙. However, they predominantly consider larger scales (several kpc) on a longer-term evolution than the study presented in this paper, as their focus is on modelling the dynamics of halo mergers and the implications for accretion.

In a single run presented in this study (T20R10), the gas fragments into three clumps instead of one. From comparison runs with the same initial conditions and a different random seed for the realization of the turbulent velocity field, we have concluded that this fragmentation does not depend on the initial amount of rotation or turbulence, as usually one, sometimes two, and in a single case three clumps are found, though never more than three. Thus, we do not find much fragmentation, irrespective of turbulence or rotation. It is not clear whether the additional clumps will continue to collapse, or instead accrete onto the main clump, though based on a simplified model of fragmentation in the accretion disk around a protostar, Inayoshi & Haiman (2014) argue that most of these clumps can migrate inward and merge with the central protostar. The simulations have been evolved for another ~1.67 × 10-2 years (~5 freefall times) after the highest refinement level was reached. As no fragmentation occurs in most of our simulations during this time, it can be considered as a lower limit on the fragmentation time scale. To quantify the amount of fragmentation with greater certainty, the simulations should be evolved for longer, though our findings at least hint at the robustness of the direct collapse scenario.

For our simulations, we have used a Lyman-Werner background with a T5 spectrum. However, we expect to find the final result to be similarly independent of turbulence or rotation, as long as the intensity of the UV background is above the critical value, regardless of the stellar spectrum.

Recently, Inayoshi et al. (2014) have done a similar simulation to attempt to resolve protostar formation, starting from equally simplified, though somewhat different initial conditions. They start from a marginally supported isothermal sphere with an initial density of 104 cm-3 and temperature of 8000 K. The mass and radius of their cloud are slightly smaller than ours, 1.17 × 105 M⊙ and 10.8 pc, respectively, though they are of the same order of magnitude. They also resolve the Jeans length by at least 64 grid cells, and their limiting resolution is 0.1 AU. There are a few differences between our chemical models. Concerning three-body rates, they do not take reaction H + H + H2 into account, and for reaction H + H + H they use the rate by Shapiro & Kang (1987), while we have used the updated rate by Forrey (2013). Additionally, their opacity corrections at high density are calculated based on the Rosseland mean opacity, while we have used a different treatment for each cooling process, as described in Sect. 2.2. For atomic cooling, we consider opacity from Rayleigh scattering, and a fudge opacity to reduce cooling at high densities, in accordance with findings from previous one-zone studies. For H− cooling, opacity from both Rayleigh scattering and H− bound-free absorption is taken into account.

They have stopped their simulation when a temperature just in excess of 104 K was reached during the adiabatic phase, and find a hydrostatic core with a mass of 1 M⊙ and a radius of 2 AU, at a peak density of ~5 × 1016 cm-3. In our simulations, the onset of the adiabatic phase occurs approximately one order of magnitude in density later, likely due to differences in the chemical model, which results in our central clump being less massive, ~7 × 10-2 M⊙, and smaller, 2 × 10-2 AU. Findings in agreement with ours include the resulting isothermal profile of the collapsing cloud, the fact that H2 cooling remains inefficient, and the accretion rate. They did not include initial rotation, and do not find a disk, similar to our run without rotation (T40R0).

Presently, there are no simulations that resolve the formation of the protostar starting from cosmological initial conditions. However, our mass infall rates are in agreement with those found in cosmological direct collapse simulations reaching lower peak densities (Latif et al. 2013a,e).

Wise et al. (2008) have performed cosmological simulations following the collapse of two atomic cooling haloes up to densities comparable to our peak density, though using a much simpler chemical model, neglecting, for example, H2 and H− chemistry and cooling. Particularly, they did not consider optical depth effects, overestimating the cooling above column densities of ~1013 cm-2, and thus did not obtain a transition towards an adiabatic equation of state. Therefore, the formation of a quasi static object, like a protostar, was not observed. This can be clearly seen from their Fig. 5: there is no increase in temperature for the inner regions, nor a significant change in slope in the density or enclosed mass profiles. These important differences make it difficult to directly compare their results to ours. However, we can say that their radial velocity and density profiles are, except in the innermost region (R< 10-6 pc), comparable to ours, meaning that the radial accretion rate should be of the same order of magnitude as well.

In this study, we did not take turbulence on subgrid scales into account (for more details on subgrid scale turbulence, see e.g. Schmidt et al. 2006; Schmidt & Federrath 2011). As has been shown by Latif et al. (2013b,a,e), turbulence on unresolved scales affects the morphology of the collapsing gas. Thus, it would be interesting to investigate whether this affects our results. Another caveat is the absence of magnetic fields. Latif et al. (2013d, 2014d) have demonstrated that even a very small seed field can be effectively amplified by the small-scale dynamo mechanism. The resulting strong magnetic field provides additional support against gravity and helps to suppress fragmentation. A further caveat is the simplification of the cooling functions, and more importantly, the opacities. Future work should include a more detailed treatment of optical depth effects, possibly even employing radiative transfer, although this will be computationally more expensive. In the future, it would also be useful to implement equilibrium chemistry, as then it will be possible to follow the evolution for longer, and study the accretion onto the protostar in more detail. Additionally, simulations that resolve protostar formation starting from cosmological initial conditions are needed, to rule out possible effects caused by an idealized setup.

Online material

Appendix A: Reaction rates

Table A.1

Reaction rate coefficients.

1

Publicly available at http://kromepackage.org/

2

The exponents of 5 and 8 do not have a specific physical meaning, but are instead intended to provide a sharp enough cutoff, as in this regime the atomic cooling functions increase steeply with both density and temperature.

Acknowledgments

Computations described in this work were performed using the ENZO code (http://enzo-project.org), which is the product of a collaborative effort of scientists at many universities and national laboratories. The simulation results are analysed using yt, a multi-code analysis toolkit for astrophysical simulation data (Turk et al. 2011). C.V.B., D.R.G.S. and M.A.L. acknowledge funding by the Deutsche Forschungsgemeinschaft (DFG) under grant SFB 963/1 (project A12). D.R.G.S., S.B., and C.V.B. thank the DFG for funding and computing time via the Schwerpunktprogram SPP 1573 “Physics of the Interstellar Medium” (grant SCHL 1964/1 – 1). T.G. acknowledges the Centre for Star and Planet Formation funded by the Danish National Research Foundation.

References

  1. Abel, T., Anninos, P., Zhang, Y., & Norman, M. L. 1997, New Astron., 2, 181 [Google Scholar]
  2. Abel, T., Bryan, G. L., & Norman, M. L. 2002, Science, 295, 93 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  3. Agarwal, B., & Khochfar, S. 2014, MNRAS, submitted [arXiv:1407.4115] [Google Scholar]
  4. Agarwal, B., Dal la Vecchia, C., Johnson, J. L., Khochfar, S., & Paardekooper, J.-P. 2014, MNRAS, 443, 648 [NASA ADS] [CrossRef] [Google Scholar]
  5. Ahn, K., Shapiro, P. R., Iliev, I. T., Mellema, G., & Pen, U.-L. 2009, ApJ, 695, 1430 [NASA ADS] [CrossRef] [Google Scholar]
  6. Aldrovandi, S. M. V., & Pequignot, D. 1973, A&A, 25, 137 [NASA ADS] [Google Scholar]
  7. Alvarez, M. A., Wise, J. H., & Abel, T. 2009, ApJ, 701, L133 [NASA ADS] [CrossRef] [Google Scholar]
  8. Ball, W. H., Tout, C. A., Żytkow, A. N., & Eldridge, J. J. 2011, MNRAS, 414, 2751 [NASA ADS] [CrossRef] [Google Scholar]
  9. Ball, W. H., Tout, C. A., & Żytkow, A. N. 2012, MNRAS, 421, 2713 [NASA ADS] [CrossRef] [Google Scholar]
  10. Begelman, M. C. 2010, MNRAS, 402, 673 [NASA ADS] [CrossRef] [Google Scholar]
  11. Begelman, M. C., & Rees, M. J. 1978, MNRAS, 185, 847 [NASA ADS] [CrossRef] [Google Scholar]
  12. Begelman, M. C., Volonteri, M., & Rees, M. J. 2006, MNRAS, 370, 289 [NASA ADS] [CrossRef] [Google Scholar]
  13. Begelman, M. C., Rossi, E. M., & Armitage, P. J. 2008, MNRAS, 387, 1649 [NASA ADS] [CrossRef] [Google Scholar]
  14. Bovino, S., Grassi, T., Latif, M. A., & Schleicher, D. R. G. 2013, MNRAS, 434, L36 [NASA ADS] [CrossRef] [Google Scholar]
  15. Bovino, S., Grassi, T., Schleicher, D. R. G., & Latif, M. A. 2014a, ApJ, 790, L35 [NASA ADS] [CrossRef] [Google Scholar]
  16. Bovino, S., Latif, M. A., Grassi, T., & Schleicher, D. R. G. 2014b, MNRAS, 441, 2181 [NASA ADS] [CrossRef] [Google Scholar]
  17. Bovino, S., Schleicher, D. R. G., & Grassi, T. 2014c, A&A, 561, A13 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  18. Bromm, V., & Loeb, A. 2003, ApJ, 596, 34 [NASA ADS] [CrossRef] [Google Scholar]
  19. Bromm, V., Kudritzki, R. P., & Loeb, A. 2001, ApJ, 552, 464 [NASA ADS] [CrossRef] [Google Scholar]
  20. Bromm, V., Coppi, P. S., & Larson, R. B. 2002, ApJ, 564, 23 [NASA ADS] [CrossRef] [Google Scholar]
  21. Bryan, G. L., Norman, M. L., O’Shea, B. W., et al. 2014, ApJS, 211, 19 [NASA ADS] [CrossRef] [Google Scholar]
  22. Capitelli, M., Coppola, C. M., Diomede, P., & Longo, S. 2007, A&A, 470, 811 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  23. Cen, R. 1992, ApJS, 78, 341 [NASA ADS] [CrossRef] [Google Scholar]
  24. Clark, P. C., Glover, S. C. O., & Klessen, R. S. 2008, ApJ, 672, 757 [NASA ADS] [CrossRef] [Google Scholar]
  25. Clark, P. C., Glover, S. C. O., Klessen, R. S., & Bromm, V. 2011a, ApJ, 727, 110 [NASA ADS] [CrossRef] [Google Scholar]
  26. Clark, P. C., Glover, S. C. O., Smith, R. J., et al. 2011b, Science, 331, 1040 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  27. Dalgarno, A., & Lepp, S. 1987, in Astrochemistry, eds. M. S. Vardya, & S. P. Tarafdar, IAU Symp., 120, 109 [Google Scholar]
  28. de Jong, T. 1972, A&A, 20, 263 [NASA ADS] [Google Scholar]
  29. Devecchi, B., & Volonteri, M. 2009, ApJ, 694, 302 [NASA ADS] [CrossRef] [Google Scholar]
  30. Devecchi, B., Volonteri, M., Colpi, M., & Haardt, F. 2010, MNRAS, 409, 1057 [NASA ADS] [CrossRef] [Google Scholar]
  31. Devecchi, B., Volonteri, M., Rossi, E. M., Colpi, M., & Portegies Zwart, S. 2012, MNRAS, 421, 1465 [NASA ADS] [CrossRef] [Google Scholar]
  32. Dijkstra, M., Haiman, Z., Mesinger, A., & Wyithe, J. S. B. 2008, MNRAS, 391, 1961 [NASA ADS] [CrossRef] [Google Scholar]
  33. Dijkstra, M., Ferrara, A., & Mesinger, A. 2014, MNRAS, 442, 2036 [NASA ADS] [CrossRef] [Google Scholar]
  34. Donahue, M., & Shull, J. M. 1991, ApJ, 383, 511 [NASA ADS] [CrossRef] [Google Scholar]
  35. Dopcke, G., Glover, S. C. O., Clark, P. C., & Klessen, R. S. 2011, ApJ, 729, L3 [NASA ADS] [CrossRef] [Google Scholar]
  36. Dopcke, G., Glover, S. C. O., Clark, P. C., & Klessen, R. S. 2013, ApJ, 766, 103 [NASA ADS] [CrossRef] [Google Scholar]
  37. Fan, X. 2006, New Astron. Rev, 50, 665 [Google Scholar]
  38. Federrath, C., & Klessen, R. S. 2012, ApJ, 761, 156 [NASA ADS] [CrossRef] [Google Scholar]
  39. Federrath, C., Sur, S., Schleicher, D. R. G., Banerjee, R., & Klessen, R. S. 2011, ApJ, 731, 62 [NASA ADS] [CrossRef] [Google Scholar]
  40. Ferrara, A., Salvadori, S., Yue, B., & Schleicher, D. 2014, MNRAS, 443, 2410 [NASA ADS] [CrossRef] [Google Scholar]
  41. Forrey, R. C. 2013, ApJ, 773, L25 [NASA ADS] [CrossRef] [Google Scholar]
  42. Glover, S. C. O., & Abel, T. 2008, MNRAS, 388, 1627 [NASA ADS] [CrossRef] [Google Scholar]
  43. Grassi, T., Bovino, S., Schleicher, D., & Gianturco, F. A. 2013, MNRAS, 431, 1659 [NASA ADS] [CrossRef] [Google Scholar]
  44. Grassi, T., Bovino, S., Schleicher, D. R. G., et al. 2014, MNRAS, 439, 2386 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
  45. Greif, T. H. 2014, MNRAS, 444, 1566 [NASA ADS] [CrossRef] [Google Scholar]
  46. Greif, T. H., Johnson, J. L., Klessen, R. S., & Bromm, V. 2008, MNRAS, 387, 1021 [NASA ADS] [CrossRef] [Google Scholar]
  47. Greif, T. H., Springel, V., White, S. D. M., et al. 2011, ApJ, 737, 75 [NASA ADS] [CrossRef] [Google Scholar]
  48. Haiman, Z. 2013, in Astrophys. Space Sci. Libr. 396, eds. T. Wiklind, B. Mobasher, & V. Bromm, 293 [Google Scholar]
  49. Hartwig, T., Clark, P. C., Glover, S. C. O., Klessen, R. S., & Sasaki, M. 2014, ApJ, submitted [arXiv:1407.2102] [Google Scholar]
  50. Hirano, S., Hosokawa, T., Yoshida, N., et al. 2014, ApJ, 781, 60 [NASA ADS] [CrossRef] [Google Scholar]
  51. Hocuk, S., & Spaans, M. 2010, A&A, 510, A110 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  52. Holzbauer, L. N., & Furlanetto, S. R. 2012, MNRAS, 419, 718 [NASA ADS] [CrossRef] [Google Scholar]
  53. Hosokawa, T., Omukai, K., Yoshida, N., & Yorke, H. W. 2011, Science, 334, 1250 [Google Scholar]
  54. Hosokawa, T., Omukai, K., & Yorke, H. W. 2012a, ApJ, 756, 93 [Google Scholar]
  55. Hosokawa, T., Yoshida, N., Omukai, K., & Yorke, H. W. 2012b, ApJ, 760, L37 [NASA ADS] [CrossRef] [Google Scholar]
  56. Hosokawa, T., Yorke, H. W., Inayoshi, K., Omukai, K., & Yoshida, N. 2013, ApJ, 778, 178 [NASA ADS] [CrossRef] [Google Scholar]
  57. Inayoshi, K., & Omukai, K. 2012, MNRAS, 422, 2539 [NASA ADS] [CrossRef] [Google Scholar]
  58. Inayoshi, K., & Haiman, Z. 2014, MNRAS, 445, 1549 [NASA ADS] [CrossRef] [Google Scholar]
  59. Inayoshi, K., Omukai, K., & Tasker, E. J. 2014, MNRAS, 445, L109 [NASA ADS] [CrossRef] [Google Scholar]
  60. Janev, R. K., Langer, W. D., & Evans, K. 1987, Elementary processes in Hydrogen-Helium plasmas – Cross sections and reaction rate coefficients (Springer) [Google Scholar]
  61. John, T. L. 1988, A&A, 193, 189 [NASA ADS] [Google Scholar]
  62. Johnson, J. L., Khochfar, S., Greif, T. H., & Durier, F. 2011, MNRAS, 410, 919 [NASA ADS] [CrossRef] [Google Scholar]
  63. Karpas, Z., Anicich, V., & Huntress, W. T. 1979, J. Chem. Phys., 70, 2877 [NASA ADS] [CrossRef] [Google Scholar]
  64. Klessen, R. S., Glover, S. C. O., & Clark, P. C. 2012, MNRAS, 421, 3217 [NASA ADS] [CrossRef] [Google Scholar]
  65. Koushiappas, S. M., Bullock, J. S., & Dekel, A. 2004, MNRAS, 354, 292 [NASA ADS] [CrossRef] [Google Scholar]
  66. Kurucz, R. L. 1970, SAO Spec. Rep., 309 [Google Scholar]
  67. Larson, R. B. 1981, MNRAS, 194, 809 [NASA ADS] [CrossRef] [Google Scholar]
  68. Latif, M. A., Schleicher, D. R. G., Schmidt, W., & Niemeyer, J. 2013a, MNRAS, 433, 1607 [NASA ADS] [CrossRef] [Google Scholar]
  69. Latif, M. A., Schleicher, D. R. G., Schmidt, W., & Niemeyer, J. 2013b, MNRAS, 430, 588 [NASA ADS] [CrossRef] [Google Scholar]
  70. Latif, M. A., Schleicher, D. R. G., Schmidt, W., & Niemeyer, J. 2013c, ApJ, 772, L3 [NASA ADS] [CrossRef] [Google Scholar]
  71. Latif, M. A., Schleicher, D. R. G., Schmidt, W., & Niemeyer, J. 2013d, MNRAS, 432, 668 [NASA ADS] [CrossRef] [Google Scholar]
  72. Latif, M. A., Schleicher, D. R. G., Schmidt, W., & Niemeyer, J. C. 2013e, MNRAS, 436, 2989 [NASA ADS] [CrossRef] [Google Scholar]
  73. Latif, M. A., Bovino, S., Grassi, T., Schleicher, D. R. G., & Spaans, M. 2014a, MNRAS, submitted [arXiv:1408.3061] [Google Scholar]
  74. Latif, M. A., Bovino, S., Van Borm, C., et al. 2014b, MNRAS, 443, 1979 [NASA ADS] [CrossRef] [Google Scholar]
  75. Latif, M. A., Schleicher, D. R. G., Bovino, S., Grassi, T., & Spaans, M. 2014c, ApJ, 792, 78 [NASA ADS] [CrossRef] [Google Scholar]
  76. Latif, M. A., Schleicher, D. R. G., & Schmidt, W. 2014d, MNRAS, 440, 1551 [NASA ADS] [CrossRef] [Google Scholar]
  77. Lodato, G., & Natarajan, P. 2006, MNRAS, 371, 1813 [NASA ADS] [CrossRef] [Google Scholar]
  78. Lupi, A., Colpi, M., Devecchi, B., Galanti, G., & Volonteri, M. 2014, MNRAS, 442, 3616 [NASA ADS] [CrossRef] [Google Scholar]
  79. Mac Low, M.-M., & Klessen, R. S. 2004, Rev. Mod. Phys., 76, 125 [NASA ADS] [CrossRef] [Google Scholar]
  80. Machida, M. N. 2008, ApJ, 682, L1 [NASA ADS] [CrossRef] [Google Scholar]
  81. Madau, P., Haardt, F., & Dotti, M. 2014, ApJ, 784, L38 [NASA ADS] [CrossRef] [Google Scholar]
  82. Martin, P. G., Schwarz, D. H., & Mandy, M. E. 1996, ApJ, 461, 265 [NASA ADS] [CrossRef] [Google Scholar]
  83. Martin, P. G., Keogh, W. J., & Mandy, M. E. 1998, ApJ, 499, 793 [NASA ADS] [CrossRef] [Google Scholar]
  84. McKee, C. F., & Ostriker, E. C. 2007, ARA&A, 45, 565 [NASA ADS] [CrossRef] [Google Scholar]
  85. Milosavljević, M., Bromm, V., Couch, S. M., & Oh, S. P. 2009a, ApJ, 698, 766 [NASA ADS] [CrossRef] [Google Scholar]
  86. Milosavljević, M., Couch, S. M., & Bromm, V. 2009b, ApJ, 696, L146 [NASA ADS] [CrossRef] [Google Scholar]
  87. Mortlock, D. J., Warren, S. J., Venemans, B. P., et al. 2011, Nature, 474, 616 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  88. Omukai, K. 2000, ApJ, 534, 809 [NASA ADS] [CrossRef] [Google Scholar]
  89. Omukai, K. 2001, ApJ, 546, 635 [NASA ADS] [CrossRef] [Google Scholar]
  90. Omukai, K. 2012, PASJ, 64, 114 [NASA ADS] [Google Scholar]
  91. Omukai, K., Tsuribe, T., Schneider, R., & Ferrara, A. 2005, ApJ, 626, 627 [NASA ADS] [CrossRef] [Google Scholar]
  92. Omukai, K., Schneider, R., & Haiman, Z. 2008, ApJ, 686, 801 [NASA ADS] [CrossRef] [Google Scholar]
  93. Palla, F., Salpeter, E. E., & Stahler, S. W. 1983, ApJ, 271, 632 [NASA ADS] [CrossRef] [Google Scholar]
  94. Portegies Zwart, S. F., Baumgardt, H., Hut, P., Makino, J., & McMillan, S. L. W. 2004, Nature, 428, 724 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  95. Poulaert, G., Brouillard, F., Claeys, W., McGowan, J. W., & Van Wassenhove, G. 1978, J. Phys. B At. Mol. Phys., 11, L671 [NASA ADS] [CrossRef] [Google Scholar]
  96. Regan, J. A., & Haehnelt, M. G. 2009a, MNRAS, 396, 343 [NASA ADS] [CrossRef] [Google Scholar]
  97. Regan, J. A., & Haehnelt, M. G. 2009b, MNRAS, 393, 858 [NASA ADS] [CrossRef] [Google Scholar]
  98. Regan, J. A., Johansson, P. H., & Haehnelt, M. G. 2014, MNRAS, 439, 1160 [NASA ADS] [CrossRef] [Google Scholar]
  99. Ripamonti, E., & Abel, T. 2004, MNRAS, 348, 1019 [NASA ADS] [CrossRef] [Google Scholar]
  100. Safranek-Shrader, C., Milosavljević, M., & Bromm, V. 2014, MNRAS, 438, 1669 [NASA ADS] [CrossRef] [Google Scholar]
  101. Schaerer, D. 2002, A&A, 382, 28 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  102. Schleicher, D. R. G., Galli, D., Palla, F., et al. 2008, A&A, 490, 521 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  103. Schleicher, D. R. G., Galli, D., Glover, S. C. O., et al. 2009, ApJ, 703, 1096 [NASA ADS] [CrossRef] [Google Scholar]
  104. Schleicher, D. R. G., Spaans, M., & Glover, S. C. O. 2010, ApJ, 712, L69 [NASA ADS] [CrossRef] [Google Scholar]
  105. Schleicher, D. R. G., Palla, F., Ferrara, A., Galli, D., & Latif, M. 2013, A&A, 558, A59 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  106. Schmidt, W., & Federrath, C. 2011, A&A, 528, A106 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  107. Schmidt, W., Niemeyer, J. C., & Hillebrandt, W. 2006, A&A, 450, 265 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  108. Schneider, R., Ferrara, A., Salvaterra, R., Omukai, K., & Bromm, V. 2003, Nature, 422, 869 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  109. Schneider, R., Omukai, K., Bianchi, S., & Valiante, R. 2012a, MNRAS, 419, 1566 [NASA ADS] [CrossRef] [Google Scholar]
  110. Schneider, R., Omukai, K., Limongi, M., et al. 2012b, MNRAS, 423, L60 [NASA ADS] [CrossRef] [Google Scholar]
  111. Sethi, S., Haiman, Z., & Pandey, K. 2010, ApJ, 721, 615 [NASA ADS] [CrossRef] [Google Scholar]
  112. Shang, C., Bryan, G. L., & Haiman, Z. 2010, MNRAS, 402, 1249 [NASA ADS] [CrossRef] [Google Scholar]
  113. Shapiro, P. R., & Kang, H. 1987, ApJ, 318, 32 [NASA ADS] [CrossRef] [Google Scholar]
  114. Smith, R. J., Glover, S. C. O., Clark, P. C., Greif, T., & Klessen, R. S. 2011, MNRAS, 414, 3633 [NASA ADS] [CrossRef] [Google Scholar]
  115. Smith, R. J., Hosokawa, T., Omukai, K., Glover, S. C. O., & Klessen, R. S. 2012, MNRAS, 424, 457 [NASA ADS] [CrossRef] [Google Scholar]
  116. Spaans, M., & Silk, J. 2006, ApJ, 652, 902 [NASA ADS] [CrossRef] [Google Scholar]
  117. Stacy, A., Greif, T. H., & Bromm, V. 2010, MNRAS, 403, 45 [NASA ADS] [CrossRef] [Google Scholar]
  118. Stacy, A., Greif, T. H., & Bromm, V. 2012, MNRAS, 422, 290 [NASA ADS] [CrossRef] [Google Scholar]
  119. Sugimura, K., Omukai, K., & Inoue, A. K. 2014, MNRAS, 445, 544 [NASA ADS] [CrossRef] [Google Scholar]
  120. Sur, S., Schleicher, D. R. G., Banerjee, R., Federrath, C., & Klessen, R. S. 2010, ApJ, 721, L134 [NASA ADS] [CrossRef] [Google Scholar]
  121. Susa, H. 2013, ApJ, 773, 185 [NASA ADS] [CrossRef] [Google Scholar]
  122. Susa, H., Hasegawa, K., & Tominaga, N. 2014, ApJ, 792, 32 [NASA ADS] [CrossRef] [Google Scholar]
  123. Tegmark, M., Silk, J., Rees, M. J., et al. 1997, ApJ, 474, 1 [NASA ADS] [CrossRef] [Google Scholar]
  124. Tumlinson, J., & Shull, J. M. 2000, ApJ, 528, L65 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  125. Turk, M. J., Abel, T., & O’Shea, B. 2009, Science, 325, 601 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  126. Turk, M. J., Smith, B. D., Oishi, J. S., et al. 2011, ApJS, 192, 9 [NASA ADS] [CrossRef] [Google Scholar]
  127. Turk, M. J., Oishi, J. S., Abel, T., & Bryan, G. L. 2012, ApJ, 745, 154 [NASA ADS] [CrossRef] [Google Scholar]
  128. Van Borm, C., & Spaans, M. 2013, A&A, 553, L9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  129. Venemans, B. P., Findlay, J. R., Sutherland, W. J., et al. 2013, ApJ, 779, 24 [NASA ADS] [CrossRef] [Google Scholar]
  130. Visbal, E., Haiman, Z., & Bryan, G. L. 2014, MNRAS, 445, 1056 [NASA ADS] [CrossRef] [Google Scholar]
  131. Volonteri, M. 2010, A&ARv, 18, 279 [NASA ADS] [CrossRef] [Google Scholar]
  132. Volonteri, M., & Begelman, M. C. 2010, MNRAS, 409, 1022 [NASA ADS] [CrossRef] [Google Scholar]
  133. Whalen, D. J., Johnson, J. L., Smidt, J., et al. 2013, ApJ, 774, 64 [NASA ADS] [CrossRef] [Google Scholar]
  134. Willott, C. J., Delorme, P., Reylé, C., et al. 2010, AJ, 139, 906 [NASA ADS] [CrossRef] [Google Scholar]
  135. Wise, J. H., Turk, M. J., & Abel, T. 2008, ApJ, 682, 745 [NASA ADS] [CrossRef] [Google Scholar]
  136. Wolcott-Green, J., Haiman, Z., & Bryan, G. L. 2011, MNRAS, 418, 838 [NASA ADS] [CrossRef] [Google Scholar]

All Tables

Table 1

Overview of the different simulations and their initial turbulent and rotational velocities.

In the text
Table 2

Properties of the central bound clump found in each simulation.

In the text
Table A.1

Reaction rate coefficients.

In the text

All Figures

thumbnail Fig. 1

Physical quantities as a function of number density in a one-zone calculation, irradiated by a strong T5 background, using our modification of KROME. A: temperature; B: self-shielding factor for H2, fsh (fsh = 1 means no shielding, and the smaller fsh, the stronger the shielding); C: number fractions of H species (and electrons); D: number fractions of He species.
In the text
thumbnail Fig. 2

Physical quantities, weighted by mass, spherically averaged and radially binned, as a function of radius at the peak density output for the different simulations. A) Number density; B) temperature; C) H2 number fraction; D) electron number fraction; E) turbulent velocity; F) radial velocity, plotted as − vrad; G) enclosed mass, and the Jeans mass for T40R10 (it is very similar for other runs); H) radial mass infall rate, calculated from the density and the radial velocity. The simulation details and abbreviations are listed in Table 1.
In the text
thumbnail Fig. 3

Density projections along the x-axis for all simulations, showing the integrated number density for various scales at the peak density output. The simulation details and abbreviations are listed in Table 1.
In the text
thumbnail Fig. 4

Density projections along the z-axis for all simulations, showing the integrated number density for various scales at the peak density output.The simulation details and abbreviations are listed in Table 1.
In the text
thumbnail Fig. 5

Density and temperature slices along the x-axis for all simulations and for two different scales at the peak density output. The simulation details and abbreviations are listed in Table 1.
In the text
thumbnail Fig. 6

Density and temperature slices along the z-axis for all simulations and for two different scales at the peak density output. The simulation details and abbreviations are listed in Table 1.
In the text
thumbnail Fig. 7

Physical quantities, weighted by mass, disk averaged and radially binned, as a function of radius at the peak density output for the different simulations. A) Toomre Q parameter; B) ratio of the rotational to the Keplerian velocity; C) sound speed; D) turbulent velocity; E) rotational frequency; F): surface density. The simulation details and abbreviations are listed in Table 1.
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.