Massive protostars in atomic cooling haloes

Formation of massive protostars in atomic cooling haloes

Abstract

We present the highest-resolution three-dimensional simulation to date of the collapse of an atomic cooling halo in the early Universe. We use the moving-mesh code arepo with the primordial chemistry module introduced in Greif (2014), which evolves the chemical and thermal rate equations for over more than 20 orders of magnitude in density. Molecular hydrogen cooling is suppressed by a strong Lyman-Werner background, which facilitates the near-isothermal collapse of the gas at a temperature of about K. Once the central gas cloud becomes optically thick to continuum emission, it settles into a Keplerian disc around the primary protostar. The initial mass of the protostar is about , which is an order of magnitude higher than in minihaloes that cool via molecular hydrogen. The high accretion rate and efficient cooling of the gas catalyse the fragmentation of the disc into a small protostellar system with members. After about yr, strong gravitational interactions disrupt the disc and temporarily eject the primary protostar from the centre of the cloud. By the end of the simulation, a secondary clump has collapsed at a distance of au from the primary clump. If this clump undergoes a similar evolution as the first, the central gas cloud may evolve into a wide binary system. High accretion rates of both the primary and secondary clumps suggest that fragmentation is not a significant barrier for forming at least one massive black hole seed.

keywords:
hydrodynamics – stars: formation – galaxies: formation – galaxies: high-redshift– cosmology: theory – early Universe.

1 Introduction

Black holes (BHs) are a key ingredient in the formation and evolution of galaxies. In the local Universe, the stellar velocity dispersion in galaxy bulges is correlated with the mass of the BH at their centre (Ferrarese & Merritt, 2000; Gebhardt et al., 2000). BHs also power luminous quasars by accreting gas from their host galaxies. Recent observations suggest that quasars powered by BHs with masses were already present when the Universe was less than one billion year old (Fan et al., 2003, 2006). These supermassive black holes most likely grew from smaller seed BHs that formed earlier, but the origin of these seeds remains unclear (Haiman, 2006, 2009; Greene, 2012; Volonteri, 2012; Volonteri & Bellovary, 2012). One possible candidate are the remnants of massive Population III stars (Madau & Rees, 2001; Li et al., 2007; Johnson et al., 2012), or the direct collapse of primordial gas in haloes with virial temperatures K, so-called atomic cooling haloes (Bromm & Loeb, 2003; Bromm & Yoshida, 2011). In the former case, the seeds have initial masses of the order of , and grow at or above the Eddington limit for the remaining Myr between seed formation and . However, numerical simulations have shown that accretion on to early BHs is inefficient, due to the low density of the gas surrounding the BH remnant, which is caused by photoionization heating from the progenitor star (Johnson & Bromm, 2007; Alvarez, Wise & Abel, 2009). Accretion rates are thus not high enough to allow efficient growth of the seed, which poses a serious complication for the Population III stellar remnant scenario.

In the direct collapse scenario, haloes with virial temperatures K may host seed BHs that are substantially more massive. A prerequisite is that the accretion rate on to the central object is high enough that radiative feedback does not severely impede the accretion flow (Johnson et al., 2011, 2012; Hosokawa, Omukai & Yorke, 2012; Hosokawa et al., 2013). In this case, a supermassive star or ‘quasi-star’ forms, which may collapse into a BH of mass (Heger et al., 2003; Begelman, Volonteri & Rees, 2006; Begelman, Rossi & Armitage, 2008; Begelman, 2010; Volonteri & Begelman, 2010; Montero, Janka & Müller, 2012; Volonteri, 2012; Inayoshi, Hosokawa & Omukai, 2013; Schleicher et al., 2013; Chen et al., 2014). Since the accretion rate in a Jeans-unstable cloud scales as , molecular hydrogen cooling must be suppressed until the virial temperature of the halo is high enough that Ly cooling becomes important. This may be achieved by a Lyman-Werner (LW) radiation background (Omukai, 2001; Bromm & Loeb, 2003; Volonteri & Rees, 2005; Spaans & Silk, 2006; Schleicher, Spaans & Glover, 2010; Johnson et al., 2013). Simple one-zone models have found that the critical flux is of the order of in units of for a blackbody spectrum with K (Omukai, 2001). For Population I/II stars, recent studies have found that the critical flux may be somewhat lower (Shang, Bryan & Haiman, 2010; Wolcott-Green & Haiman, 2012; Van Borm & Spaans, 2013; Agarwal & Khochfar, 2014; Latif et al., 2014b, a; Regan, Johansson & Wise, 2014; Sugimura, Omukai & Inoue, 2014). Even though the LW flux on cosmological scales is well below this value, local star formation may raise the flux to supercritical levels (Dijkstra et al., 2008; Dijkstra, Ferrara & Mesinger, 2014; Agarwal et al., 2012, 2014; Visbal, Haiman & Bryan, 2014).

If the LW flux is high enough, the halo gas collapses nearly isothermally at K up to a density of , where the gas becomes optically thick to Ly emission (Omukai, 2001). At this point, continuum cooling via free-bound emission of H takes over, and allows the gas to again contract nearly isothermally up to a density of . Once the continuum emission becomes trapped, the gas evolves nearly adiabatically and a protostar forms at the centre of the halo. During the initial collapse, the angular momentum is constantly redistributed by turbulence and bar-like instabilities, such that the cloud contracts nearly unhindered (Oh & Haiman, 2002; Koushiappas, Bullock & Dekel, 2004; Begelman, Volonteri & Rees, 2006; Lodato & Natarajan, 2006; Wise, Turk & Abel, 2008; Begelman & Shlosman, 2009; Choi, Shlosman & Begelman, 2013; Latif et al., 2013a; Prieto, Jimenez & Haiman, 2013).

The subsequent accretion phase was investigated by Regan & Haehnelt (2009) and Latif et al. (2013b, a). They found that a Keplerian disc forms around the primary protostar, which becomes gravitationally unstable and fragments into a small system of protostars. The secondary protostars merge on a short time-scale and do not prevent the growth of the primary protostar. These studies employed a pressure floor beyond a certain refinement level, such that the maximum density was limited to . The simulations of Regan, Johansson & Haehnelt (2014) also displayed the formation of a disc-like object at the centre of the halo, which in some cases fragmented on a scale of au. However, these simulations also suffered from limited resolution, and did not include the relevant cooling and chemistry. Recently, Inayoshi, Omukai & Tasker (2014) used the most detailed chemical and thermal model to date, but stopped the simulation once the primary protostar had formed. In addition, they did not use cosmological initial conditions. We here attempt to improve upon these studies by carrying out a simulation that starts from cosmological initial conditions and is not resolution-limited. We use a slightly less sophisticated chemical model as Inayoshi, Omukai & Tasker (2014), but evolve the simulation well beyond the formation of the first protostar at the centre of the halo.

Our paper is organized as follows. In Section 2, we describe the simulation setup and the chemistry and cooling network. In Section 3, we analyse the simulation and discuss the collapse of the central gas cloud, the formation and fragmentation of the disc, the development of the protostellar system, and the collapse of a secondary clump towards the end of the simulation. Finally, in Section 4 we summarize and draw conclusions. All distances are quoted in proper units, unless noted otherwise.

2 Simulations

We perform three-dimensional, cosmological hydrodynamical simulations to investigate the collapse of gas in atomic cooling haloes in which the formation of has been suppressed by a LW background. For this purpose we employ the moving-mesh code arepo (Springel, 2010). We also include the recently developed primordial chemistry and cooling network of Greif (2014). In the following, we briefly describe the initialization of the simulations, the extraction procedure and refinement criteria used to achieve densities , and the chemistry and cooling network.

2.1 Dark matter simulations

We first initialize a dark matter (DM)-only simulation at a redshift of in a standard cold dark matter (CDM) cosmology. We adopt cosmological parameters based on the Wilkinson Microwave Anisotropy Probe results (Komatsu et al., 2009). We use a matter density , baryon density , Hubble parameter (where is the present Hubble expansion rate), spectral index , and normalization . The simulation is initialized in a box of side length Mpc (comoving) with a total of DM particles of mass . The gravitational softening length is set to pc (comoving), which corresponds to 5% of the initial mean inter-particle separation. We stop the simulation when the first halo with virial mass exceeding collapses. This occurs at , when the first halo reaches . At this point the halo has a virial radius of kpc and a spin parameter .

2.2 Resimulations

The second step is to locate the target halo and flag it for further refinement. We select the particles belonging to that halo and a sufficiently large boundary region around it, and trace them back to their initial conditions. Once the particle locations have been determined, we reinitialize the simulation centred on the target halo. In order to acquire higher resolution we replace each DM particle by 64 less-massive DM particles and 64 mesh-generating points. The resolution is gradually decreased as the distance from the high-resolution region increases, replacing cells and DM particles by higher-mass particles outside the target region. The resimulation has lower resolution towards the edges of the box than the original DM-only simulation, but the accuracy of the gravitational tidal field around the target halo is preserved. The refined DM particle mass is given by , and the gravitational softening length is set to pc (comoving). The refined mass of each cell is given by .

We stop the resimulation once the first cell has exceeded a density of . We then proceed to extract the particles in the central pc and reinitialize the simulation with reflective boundary conditions. Hence, the central region of the final output in the first resimulation becomes the initial condition for a second resimulation with a box size of pc. Furthermore, at those densities the gas component is already well decoupled from the DM component (Greif et al., 2011), so we discard the DM and keep only the gas particles. We evolve the second resimulation until it exceeds a density of , after which we conduct a second extraction similar in nature to the first, but cut out the central pc of the second resimulation. We then use this as the side length for the third resimulation. This approach has the risk that perturbations from the edges of the box might influence the central regions. However, we explicitly avoid this issue by assuring that the sound crossing time through the box is much longer than the free-fall time of the central high-density cloud.

2.3 Refinement

An essential refinement criterion that grid codes have to fulfill to resolve gravitational instability and avoid artificial fragmentation is the so-called Truelove criterion (Truelove et al., 1997). This criterion states that the local Jeans length needs to be resolved by at least four cells, where the cell size is approximately given by , and is the volume of the cell. In order to adequately resolve turbulence, recent studies using grid codes with a fixed mesh have found that the Jeans length must be resolved by at least 32 cells (Federrath et al., 2011; Turk et al., 2012; Latif et al., 2013a). A disadvantage of using refinement based on the Jeans length is that shock-heated regions may be much less resolved than adjacent cold regions. In order to avoid this problem, we follow the refinement criterion proposed by Turk, Norman & Abel (2010), who suggest using the minimum temperature of the gas to evaluate the Jeans length. We slightly modify this criterion by using K for cells with , but the correct temperature for cells with . This ensures that the initial collapse phase is adequately resolved, while at high densities the resolution does not become excessively high and slows down the calculation. Below , we employ 64 cells per Jeans length, which is degraded to 8 cells above . We use a linear extrapolation between these densities. The maximum spatial resolution achieved with this refinement strategy is au. Next to the Jeans refinement, we refine a cell if its mass increases to more than twice its initial mass.

2.4 Chemistry and cooling

A detailed description of the chemical and thermal model used here can be found in Greif (2014). Here, we only briefly describe the most important reactions and cooling processes. The chemical network employs a non-equilibrium solver at low densities and an equilibrium solver at high densities for the species H, , , , and . The transition from non-equilibrium to equilibrium H chemistry occurs at , since three-body reactions depend on the cube of the density and would otherwise prohibitively decrease the time-step of the non-equilibrium solver. For densities above , the electron and abundances are also considered to be in equilibrium. The main reactions include the formation of via associative detachment as well as three-body reactions, the destruction of via collisions and photodissociation, and the formation and destruction of by collisional ionizations and recombinations.

Figure 1: Zoom-in on the gas cloud that forms at the centre of the atomic cooling halo. The number density of hydrogen nuclei is weighted with the square of the density along the line of sight, which is perpendicular to the plane of the disc. Clockwise from the top left, the width of the individual cubes are pc, pc, pc, au, au, and 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.

The relevant cooling processes are line cooling, collision-induced emission, Ly cooling, and inverse Compton cooling. cooling plays a substantial role up to , where the gas becomes optically thick to the line emission, while collision-induced emission becomes important at and provides the last radiative cooling channel (Omukai & Nishi, 1998; Ripamonti & Abel, 2004). Although we include molecular hydrogen cooling, its effect does not become important during the evolution of the simulation due to the presence of a strong LW background that dissociates via the Solomon process (Abel et al., 1997). Previous studies found that a strong LW flux with is required to dissociate molecular hydrogen in the progenitors of an atomic cooling halo (Omukai, 2001; Johnson & Bromm, 2007; Dijkstra et al., 2008; Latif et al., 2013b; Wolcott-Green, Haiman & Bryan, 2011). Here, we assume a constant LW flux of for a blackbody spectrum with , which is commonly used to estimate the spectra of Population III stars. In this case, the photodissociation rate is much smaller than the photodissociation rate (Sugimura, Omukai & Inoue, 2014). We approximate the combined effects of Ly cooling and continuum cooling by assuming that Ly cooling remains optically thin up to densities . The cooling rate is exponentially suppressed at densities to approximately reproduce the density-temperature relation found in Omukai (2001). Due to this simplification, we may somewhat underestimate the true cooling rate.

3 Results

3.1 Collapse of central gas cloud

A number of studies have discussed the properties of the collapse of primordial gas clouds in atomic cooling haloes (e.g., Bromm & Loeb, 2003; Regan & Haehnelt, 2009; Choi, Shlosman & Begelman, 2013; Latif et al., 2013a; Inayoshi, Omukai & Tasker, 2014; Regan, Johansson & Haehnelt, 2014). Here, we investigate the collapse of the gas over an unprecedented range in scale, as shown in Fig. 1. The six panels show a zoom-in on the central gas cloud, ranging from pc down to scales of au. The panel on the bottom-left side of the figure shows the primary protostar surrounded by an accretion disc. The cloud shows an irregular morphology and changes shape as it collapses. Its filamentary structure is indicative of turbulence, which is especially pronounced during the later stages of the collapse. On the largest scales, the cloud shows less substructure and is more spherically symmetric.

Figure 2: Radial profiles for the mass-weighted number density of hydrogen nuclei, enclosed gas mass, temperature, and , , and Hii abundances between about au and kpc. The various line styles and colours denote different epochs of the evolution of the halo, labelled by their lookback time as measured from the end of the run according to the legend. The blue dashed line corresponds to redshift , the green dash-dotted line to , the red dotted line to when the number density first exceeds , the cyan dashed line to yr after that, the purple dash-dotted line to when the number density first exceeds , the yellow dotted line to yr after the formation of first protostar, and the black solid line to the end of the simulation after approximately yr. The halo follows several evolutionary stages from large to small scales: shock-heating to the virial temperature, onset of cooling, Jeans instability, isothermal contraction, formation of the primary protostar and disc, and fragmentation of the disc (see Section 3 for details).
Figure 3: Clockwise from the top left panel: distribution of the gas in temperature, , Hii , and abundances versus number density of hydrogen nuclei at the end of the simulation. The mass per bin over the total mass in the computational domain is colour-coded from blue (lowest) to red (highest). The solid black lines show the mass-weighted average values. After shock-heating to the virial temperature of K, the gas collapses nearly isothermally to densities of . The gas then becomes optically thick to continuum emission and evolves nearly adiabatically. At this point, the Hii abundance dramatically increases from to unity. The H abundance stays below due to the LW background, but then increase to as three-body reactions set in. However, the H abundance never becomes high enough for H cooling to become important. The ‘fingers’ visible in the various distributions show the evolutionary paths of individual protostars.

Fig. 2 shows various physical quantities as a function of distance from the densest cell in the halo. The radial profiles are constructed from data of the three resimulations. We proceed by extracting the inner au from the last resimulation, while the range between and au is taken from the second resimulation. To complete the profiles, the outer region corresponds to data from the first resimulation up to au. Due to the self-similarity of the collapse, moving from large to small radii is equivalent to moving from early to late times. Properties plotted in the figure are the number density of hydrogen nuclei, enclosed gas mass, temperature, abundance, abundance, and Hii abundance. These profiles have been calculated using mass-weighted averages of the cells contributing to the radial bins. Colours and line styles represent different evolutionary stages of the gas cloud as described in the legend and the caption of the figure.

Figure 4: Mass-weighted average radial velocity, rotational velocity, Keplerian velocity, turbulent velocity, and sound speed versus radius (left-hand panel), as well as the corresponding Mach numbers (right-hand panel). The turbulent velocity is supersonic with a Mach number of throughout the collapse, while the radial velocity briefly becomes subsonic. Once the disc forms, the rotational velocity becomes nearly equal to the Keplerian velocity, and comparable to the turbulent velocity.
Figure 5: Left: from top left to bottom right, the panels show the mass-weighted average surface density, sound speed, orbital frequency, and Toomre parameter versus radius just before the disc fragments. Above au, the power-law profiles of the surface density and rotation speed yield a Toomre parameter that is close to unity, which indicates that perturbations in the disc can grow. Right: effective equation of state, root-mean-squared density contrast, cooling time over free-fall time, and free-fall time over sound-crossing time. The isothermal collapse of the gas on scales au results in , while the increasing optical depth of the gas to continuum emission on smaller scales results in an exponent that is closer to that of an adiabatic gas. The cooling time over the free fall time has a local minimum on a scale of au: this is approximately the radius at which the first fragment forms. The density contrast created by the supersonic turbulence is between and . The free-fall time exceeds the sound crossing time on a scale of au, which shows the size of the central, Jeans-unstable clump.

As the gas collapses into the DM halo, it is shock-heated to the virial temperature. In the central parts of the halo, Ly cooling becomes important and keeps the gas nearly isothermal at K (Wise & Abel, 2007). During this period, the abundance builds up from to , with small spikes due to the existence of shocks in the outer regions of the halo. The Hii abundance increases by about one order of magnitude. The collapse then approximately follows the Larson-Penston solution for an isothermal, self-gravitating gas cloud (Larson, 1969; Penston, 1969), which is described by a density profile following . The strong LW background suppresses the formation of and maintains an abundance of down to scales of au. This prevents cooling, which in turn leads to a roughly isothermal collapse between and au. Over these scales, the and Hii abundances drop by many orders of magnitude due to recombination. On a scale of au, the fraction increases due to three-body reactions. Up to this point, the radial profiles agree well with those of previous studies (Wise, Turk & Abel, 2008; Latif et al., 2013a; Regan, Johansson & Haehnelt, 2014; Inayoshi, Omukai & Tasker, 2014). In the final stage of the collapse, when the primary protostar forms, the gas becomes optically thick to continuum cooling at au, which results in a rise in temperature of more than two orders of magnitude to K. This is accompanied by a drop in the abundance and an increase in both the and Hii abundances. At the end of the simulation, two pronounced spikes in the density profile are clearly visible in the central au. These correspond to secondary protostars that have formed due to fragmentation in the disc. This will be discussed in detail in Section 3.2.

In Fig. 3, we show the temperature-density distribution of the gas at the end of the simulation, next to those of the , , and Hii abundances. At low densities, the temperature distribution spans almost six orders of magnitude, reaching as high as . A similarly high scatter is present in the and abundances, while the Hii abundance varies only by two orders of magnitude. Up to , the temperature distribution becomes much narrower, showing the near-isothermal collapse of the gas. Once three-body reactions become important, the distribution of the fraction widens for densities in the range , with particles reaching abundances as high as . The resulting temperature dispersion leads to an increasing dispersion in the and Hii abundances, while their average values continue to decrease due to recombinations. The values of the abundance are somewhat smaller than those found in Inayoshi, Omukai & Tasker (2014), but agree with Latif et al. (2013a). We therefore do not distinguish between two thermal phases of the gas as in Inayoshi, Omukai & Tasker (2014). For densities , the formation of the primary and secondary protostars can be recognized as ‘fingers’ of gas in the individual panels, which evolve nearly adiabatically. The high temperatures in the interior of the protostars results in a decrease of the and abundances, and an increase of the Hii abundance to unity.

The radial profiles of the magnitude of the radial velocity, rotational velocity, Keplerian velocity, turbulent velocity, and sound speed at the end of the simulation are shown in the left-hand panel of Fig. 4. In addition, the right-hand panel shows the Mach numbers of each velocity component. The turbulent Mach number is given by

(1)

where denotes the sound speed of the radial bin, the total mass, the index of a cell contributing to the bin, its mass, the velocity, the radial velocity vector, and the rotational velocity vector. During the initial free-fall phase, the turbulent component is supersonic with . In contrast, the Mach number of the rotational velocity remains below unity, indicating the poor rotational support of the cloud at that stage. The trend for each component is roughly maintained once the halo has entered the isothermal collapse phase, with the exception of the Mach number of the radial velocity, which briefly drops to below unity. Down to au, the rotational velocity oscillates between 0.2 and 0.5 of the Keplerian velocity, indicating a substantial degree of rotational support. It reaches its peak at the edge of the disc on scales au, where . On smaller scales, the primary protostar is characterized by an increase in temperature and thus sound speed, such that all velocity components become subsonic. In addition, drops precipitously, which shows that the infall rate decreases rapidly within the primary protostar. Similar values for the velocity have been found in previous studies (e.g. Regan, Johansson & Haehnelt, 2014).

Overall, we find good agreement between our results and previous work. However, some differences exist. The morphology of the halo between and is similar to that of Inayoshi, Omukai & Tasker (2014), but we do not find clumps on larger scales as pointed out by Regan & Haehnelt (2009), Latif et al. (2013a), and Regan, Johansson & Haehnelt (2014). However, this is not surprising, since in our case the gas has not yet had time to settle into a disc on these scales. The radial profiles resemble those of Latif et al. (2013a) quite well, while Inayoshi, Omukai & Tasker (2014) found a slightly higher abundance, which is also reflected in lower temperatures during the isothermal collapse phase. These differences may be caused by the different chemical networks used in our study (see Section 3.6).

Figure 6: Enclosed gas mass over the mass-weighted average BE mass as a function of enclosed gas mass. Colours and line styles are the same as in Fig. 2. As the halo grows in mass, the BE mass increases due to the rise in the virial temperature, which reduces . Once the atomic cooling halo is assembled, this ratio exceeds unity on a scale of . Following the onset of runaway cooling due to Ly emission, the central become Jeans-unstable (red dotted line). The minimum Jeans mass of the cloud is indicated by the purple dash-dotted line as , which coincides with the initial mass of the primary protostar.

3.2 Disc formation and fragmentation

Figure 7: Evolution of the protostellar system that forms at the centre of the atomic cooling halo. The number density of hydrogen nuclei is weighted with the square of the density along the line of sight, which is perpendicular to the plane of the disc. The top three rows show cubes with a side length of au, centred on the position of the first protostar. The bottom row shows the later evolution of the protostellar system on a somewhat larger scale of au, where the centre has been fixed on the position of the primary protostar after yr. The time is measured from the instant when the density first exceeds . The formation of a Keplerian disc around the primary protostar is clearly visible. Shortly thereafter, the disc becomes Toomre-unstable and spiral arms form that transport mass inwards and angular momentum outwards. After yr, the disc becomes gravitationally unstable and fragments due to the high mass accretion rate from the surrounding cloud on to the disc, and the efficient cooling of the disc by continuum emission. Over the next yr, an additional five protostars form before three-body interactions lead to the temporary ejection of the primary protostar from the cloud, which disrupts the disc.

After the formation of the first protostar, the gas becomes fully rotationally supported in a Keplerian disc. We study its stability by computing Toomre’s parameter (Toomre, 1964):

(2)

where is the sound speed of the gas, the epicyclic frequency of the disc, the gravitational constant, and the surface density of the gas. For the case of a Keplerian disc, the epicyclic frequency may be replaced by the orbital frequency . The parameter was originally proposed to determine whether perturbations can grow in an infinitely thin, isothermal disc. Later studies have extended this criterion for thick discs, finding that it only deviates by a factor of order unity from the above equation (Wang et al., 2010). For values greater than , the system is stable due to gas pressure and shear by the differential rotation of the disc, while for lower values the system is unstable and hence susceptible to the growth of perturbations. These lead to the formation of spiral arms that transport mass inwards and angular momentum outwards.

Figure 8: Relative distance (left-hand panel) and velocity (right-hand panel, blue solid) between the primary and a secondary protostar. The latter initially orbits around the primary protostar at a distance of au, but strong gravitational forces due to three-body interactions temporarily eject both protostars from the centre of the cloud after yr. The relative velocity reaches a peak value of about , which declines to towards the end of the simulation. For comparison, we also show the escape velocities of both protostars. The green dashed line corresponds to the primary protostar, and the red dotted line to the secondary protostar. Since the relative velocity decreases to well below the escape velocity, both protostars will likely return to the centre of the cloud.

Radial profiles for the gas surface density, sound speed, orbital frequency, and Toomre parameter are shown in the left-hand panel of Fig. 5. We compute the profiles using mass-weighted spherical shells centred on the densest cell in the halo of the final resimulation. The surface density increases as in the interior of the protostar, where the density is almost constant. On larger scales, the radial dependence changes to , as deduced from the relation for isothermal collapse, while the orbital frequency roughly follows . On scales between 0.1 and au, the radial dependence of and thus cancel each other, such that remains roughly constant around unity. In the interior of the protostar, increases due to the increase in the sound speed and the different radial scaling between and . Since the value of is roughly equal to the critical value, the disc is prone to perturbation growth.

Further properties at the time when the primary protostar has just formed and is surrounded by a disc that has not yet fragmented are shown in the right-hand panel of Fig. 5. From the top left to the bottom right, the panels show the effective equation of state, root-mean-squared density contrast, cooling time over free-fall time, and free-fall time over sound-crossing time. In the outer region of the disc, on scales au, the equation of state is characterized by , as expected for isothermal collapse. The density contrast is roughly constant around unity, while the interior of the primary protostar is characterized by values close to 0.1. Here, increases to as a reflection of the temperature increase by almost two orders of magnitude in the central au. The cooling time remains well above the free-fall time in the inner au of the halo and down to au, the scale at which the gas becomes optically thick to continuum cooling. The free-fall time remains below the sound-crossing time down to au, showing the gravitational instability of the cloud down to this scale.

3.3 Minimum fragment mass

Further evidence for the gravitational instability of the gas is presented in Fig. 6, where we plot the enclosed gas mass over the locally estimated Bonnor-Ebert (BE; Ebert, 1955; Bonnor, 1956) mass as a function of enclosed gas mass. Colours and line styles are the same as in Fig. 2. The profiles have been computed using spherical shells centred on the densest cell, where the BE mass is calculated as the mass-weighted average among cells within a given radius according to:

(3)

During the initial collapse, the ratio of enclosed gas mass to BE mass decreases as a consequence of the rise in temperature as the gas is shock-heated. The enclosed gas mass surpasses at , which is in agreement with the mass of the halo. As the halo keeps accreting, another region where the ratio exceeds unity emerges at about . This marks the initial Jeans instability of the cloud. From , the point where surpasses moves down to when the densest cell first reaches . This is the minimum fragment mass and coincides with the initial mass of the protostar formed at the centre of the halo. From then on the temperature of the central object increases, which is translated into an increase of the BE mass, and hence a decrease of the ratio. As a result, the point at which this ratio equals unity briefly moves up to and always stays above .

3.4 Protostellar system

The fragmentation of the disc into a small protostellar system is shown in Fig. 7. The top three rows show cubes of side length au centred on the position of the primary protostar, while the cubes of the last row are au wide with the centre fixed on the position of the primary protostar after yr. In total, we present 16 different output times, which are measured with respect to the point in time at which the densest cell first exceeds . During the first yr, perturbations grow between and au in the form of spiral arms. After yr, they become gravitationally unstable and the first secondary protostar forms. In the next yr, the efficient cooling of the gas results in the formation of additional protostars, and after a small protostellar system with six members has emerged. Shortly thereafter, three-body interactions and strong tidal forces during a close passage of a secondary protostar and the primary results in the disruption of the disc. Both protostars are ejected from the centre of the cloud. The sequence in the bottom row of Fig. 7 shows the evolution of this interaction and how both protostars move away from each other.

Figure 9: Stellar mass, radius, and accretion rate of all protostars formed in the simulation. Each line corresponds to an individual protostar, and the black dashed lines shows the total mass and accretion rate, respectively. Initially, the mass budget is entirely dominated by the primary protostar (blue line), which grows from to at a rate of , while its radius swells to well over au. Once the primary protostar is expelled from the centre, its accretion rate and size drop significantly. The protostar formed in the secondary clump (red line) grows to about at a rate of . The other protostars stay below , and thus do not contribute significantly to the total protostellar mass.

To quantify the interaction between both protostars, Fig. 8 shows the relative distance and velocity between the protostars over time. For comparison, we also include the escape velocity of both protostars, using the enclosed mass in a spherical region around their respective centres, with radii equal to their separation. Once the secondary protostar forms, it orbits at a roughly constant distance of au from the primary protostar, but they soon move together and their separation decreases to au. Shortly thereafter, three-body interactions eject both protostars from the centre of the halo. This is reflected by a high relative velocity with a peak value of , which is followed by a gradual drop in the relative velocity. The parabolic shape of the relative distance suggests that it may reach a point of turnaround and the protostars will begin to re-collapse towards the centre. This trend is supported by the declining profile of the relative velocity and the fact that it has fallen well below the escape velocity by the end of the simulation.

In Fig. 9, we show the mass, radius, and accretion rate of all protostars over time. The solid lines correspond to individual protostars, while the black dashed lines denote the total mass and accretion rate, respectively. The radius of a protostar is calculated as the distance at which the Rosseland mean opacity reaches its maximum value (Stacy et al., 2013). The protostellar mass is given by the mass enclosed within that radius, and the accretion rate by the time derivative of the enclosed mass. A total of eight secondary protostars form during the evolution and fragmentation of the disc. Out of these, four survive until the end of the simulation. The rest merge with other protostars or are tidally disrupted. During the first yr after the formation of the primary protostar, its mass builds up from to at a rate of roughly , in agreement with previous work (Latif et al., 2013a; Inayoshi, Omukai & Tasker, 2014). Its radius increases from to . The second protostar forms after yr with an initial mass of and a radius of . Most of the gas is accreted by the primary protostar, while the second protostar only accretes at a rate of before it is tidally disrupted.

Figure 10: Simultaneous collapse of a secondary gas clump at a distance of au from the primary clump in the atomic cooling halo. The number density of hydrogen nuclei is weighted with the square of the density along the line of sight in a cube with a side length of au. The panels on the right show zoom-ins on the primary and secondary clumps with a side length of 60 and au, respectively. While strong interactions occur in the central protostellar system, a second clump has collapsed and is in the early stages of its evolution. Ultimately, the clumps may evolve into a wide binary system.

Shortly thereafter, the disc fragments vigorously and gives rise to a protostellar system characterized by a massive primary protostar with and a radius of , while the secondary protostars only have masses between and , and radii in the range . The accretion on to the secondary protostars results in a slight decrease of the accretion rate on to the primary protostar of , while the total accretion rate remains roughly constant at . After , the primary protostar is expelled from the centre of the halo and the disc is disrupted (see bottom row of Fig. 7). As the primary protostar is deprived of gas, its accretion rate drops to , and its radius decreases from to . Its final mass is .

The formation of a protostellar system has recently been reported in studies of minihaloes that cool via lines (Clark, Glover & Klessen, 2008; Clark et al., 2011; Greif et al., 2011). Initial masses of protostars in atomic cooling haloes are an order of magnitude higher, while the accretion rates exceed those in minihaloes by about three orders of magnitude. Other studies have found similar values for the initial protostellar masses and accretion rates (Regan & Haehnelt, 2009; Latif et al., 2013a; Regan, Johansson & Haehnelt, 2014; Inayoshi, Omukai & Tasker, 2014).

3.5 Secondary clump

About yr into the evolution of the protostellar system, a second clump collapses at a distance of about au from the primary clump. Fig. 10 shows both clumps in a cube with a side length of au at the end of the simulation, with smaller cubes showing zoom-ins on to the individual clumps. The zoom-in of the secondary clump shows a protostar with a disc and spiral arms similar to the early evolutionary stages of the primary clump. The protostar in the secondary clump is denoted by the red line in Fig. 9. Its mass quickly grows from to , and its radius increases from to . Despite its later formation, it accretes more rapidly than the first protostar. Ultimately, both clumps may evolve into a wide binary system.

3.6 Caveats

Previous studies that investigated the collapse and fragmentation of gas in atomic cooling haloes did not have sufficiently high resolution to self-consistently follow the formation of protostars at the centre of the cloud. We have attempted to address this shortcoming by performing a simulation that is not resolution-limited. Nevertheless, we have neglected some physical processes that might affect the fragmentation of the cloud. In particular, we have assumed that the optically thin regime for atomic hydrogen cooling extends up to densities . In reality, the gas becomes optically thick to Ly radiation at densities of , and then free-bound continuum emission of becomes the main cooling agent. Previous studies have found that this kind of cooling may lower the temperature by up to a factor of in the range compared to our study (Omukai, 2001; Inayoshi, Omukai & Tasker, 2014). A lower temperature should translate into a lower Toomre parameter, which would enhance the fragmentation seen in our simulation. In addition, at we have introduced an artificial cut-off for continuum cooling in order to approximately reproduce the density-temperature relation found in Omukai (2001). This simplification may also affect the thermal and gravitational stability of the gas.

Another factor that might influence the temperature of the disc is the heating from the accretion luminosity of the primary protostar. The accretion luminosity is given by

(4)

where is the Planck mean opacity, the distance from the source, and the accretion luminosity. The effects of the accretion luminosity have been discussed in similar studies that focused on minihaloes (Greif et al., 2011; Smith et al., 2011). They found that the additional heating of the gas may slightly delay fragmentation, but does not prevent it. The photospheric temperature of the protostar of K during the early stages of the collapse is too low to produce significant amounts of ionizing radiation. Latif et al. (2013a) investigated the influence of accretion luminosity in atomic cooling haloes. Assuming a power-law relation between the mass and the radius of the star, and an accretion rate of , they computed an accretion luminosity of for a clump with a size of au and a temperature of . This value is comparable to the energy emitted by Ly cooling, and may exceed it once the mass of the clump reaches . However, Latif et al. (2013a) found that this difference only translates into an increase of the temperature by . Since we investigate the evolution of the protostellar system at even earlier times, when the mass of the protostar is much lower, the effects of the accretion luminosity are expected to be even smaller.

Next to the aforementioned cooling and heating processes, we do not include the effects of magnetic fields. These are expected to become dynamically important in minihaloes as well as atomic cooling haloes (e.g. Xu et al., 2008; Schleicher, Spaans & Glover, 2010; Sur et al., 2010; Peters et al., 2012, 2014; Schober et al., 2012; Turk et al., 2012; Latif et al., 2013c). Indeed, Latif, Schleicher & Schmidt (2014) found that the magnetic pressure provides additional support against gravity and delays or suppresses fragmentation. Future simulations should therefore include magnetic fields as well as a more detailed chemical and thermal model.

4 Summary and Conclusions

We have performed the highest-resolution cosmological simulation to date of the formation and evolution of a protostellar system in an atomic cooling halo. We follow the collapse of the gas from a few Mpc down to au, spanning almost 13 orders of magnitude in scale, and reaching densities as high as . The simulation includes an equilibrium/non-equilibrium primordial chemistry solver that evolves five species (H, , , , and ), and includes line emission, collision-induced emission, Ly cooling, and inverse Compton cooling. Additionally, we have included a uniform LW background of strength to prevent star formation in progenitor haloes.

During the initial collapse, the gas is shock-heated to the virial temperature of about . The molecular hydrogen abundance briefly increases due to the presence of supersonic shocks, but the external radiation background photodissociates H to a level of within the halo. As a result, runaway collapse due to Ly cooling ensues once the virial mass has risen to . The central gas cloud becomes Jeans-unstable with a mass of and collapses nearly isothermally over many orders of magnitude in density, characterized by a profile of the form . At densities , the gas becomes optically thick to Ly emission and effectively cools via free-bound continuum emission of up to a density of , where the continuum emission is trapped. The average abundance increases to at due to three-body reactions, but never becomes high enough for H line emission to become important. The abundance declines to due to recombinations before increasing to unity for densities , where the gas evolves nearly adiabatically and a protostar with an initial mass of is formed.

Following the formation of the primary protostar, the gas settles into a Keplerian disc. The Toomre parameter within the disc is close to unity, such that perturbations can grow. The emerging spiral arms feed gas on to the primary protostar at a rate of . However, this is not sufficient to process the mass that accretes from the surrounding cloud on to the disc. In combinations with the efficient cooling of the gas via continuum emission, the disc becomes gravitationally unstable and a secondary protostar forms after only yr. The disc continues to fragment, such that after yr a total of eight secondary protostars have formed. By the end of the simulation, four of these have survived, while the rest have merged away or are tidally disrupted. The primary protostar has grown to a mass of , while all other secondary protostars have masses . Three-body interactions lead to the temporary ejection of the primary protostar from the disc after yr, which is disrupted in the process. However, an analysis of the relative velocity of the protostars shows that it is well below the escape velocity. It will therefore likely return to the centre of the cloud. After , a second clumps collapses at a distance of au from the primary clump. It has not yet fragmented and contains a single protostar that rapidly grows to . If this clump show a similar pattern of rapid migration and merging, the cloud may evolve into a wide binary system.

Despite the temporary ejection of the primary protostar from the centre of the cloud, subfragmentation likely does not substantially impede its growth. Once it returns to the centre of the cloud, its accretion rate will likely again increase to . In addition, the secondary protostars formed in the disc quickly migrate to the centre of the cloud, where they merge with the primary protostar. They are also typically times less massive than the primary protostar, which has accreted by the end of the simulation, while the most massive secondary protostar has only grown to . Most of the accreted material thus does not stem from other protostars, but from the bar-like instabilities in the disc. The secondary clump may be a much more potent candidate for accreting mass that may have otherwise been accreted by the primary clump, but even in this case the growth of the most massive protostar would be reduced by at most a factor of . It thus appears that fragmentation is not a significant barrier for forming at least one massive BH seed per atomic cooling halo, assuming that the LW background is high enough to prevent H cooling. Recent simulations have shown that this may indeed be the much more limiting factor (Latif et al., 2014b, a; Regan, Johansson & Wise, 2014).

One of the main caveats of this study is the simplified chemistry and cooling network. Future work should include a more detailed chemical model, such as that used in Inayoshi, Omukai & Tasker (2014). It may also become possible to treat the radiative transfer of the various line and continuum processes (e.g. Greif, 2014). Finally, the influence of magnetic fields may be investigated with modules that have already been implemented in arepo (Pakmor, Bauer & Springel, 2011). The influence of the radiation may not be that strong, since it is difficult to heat the gas above K, while magnetic fields may have a substantial effect on the thermal and gravitational stability of the cloud (e.g. Latif et al., 2013c; Latif, Schleicher & Schmidt, 2014). The additional support provided by magnetic fields may reduce the ability of the gas to fragment, and further increase the accretion rate of the primary protostar.

Acknowledgements

FB would like to thank Simon Glover, Paul Clark, John Regan, Muhammad Latif, and Kohei Inayoshi for stimulating discussions and feedback during the conference ‘The physics of first stars and galaxy formation’ at The Higgs Centre for Theoretical Physics of the University of Edinburgh. VS acknowledges support by the European Research Council under ERC-StG EXAGAL-308037. The simulations were carried out at the Texas Advanced Computing Center (TACC) under XSEDE allocation AST130020.

References

  1. Abel T., Anninos P., Zhang Y., Norman M. L., 1997, New Astron., 2, 181
  2. Agarwal B., Dalla Vecchia C., Johnson J. L., Khochfar S., Paardekooper J.-P., 2014, MNRAS, 443, 648
  3. Agarwal B., Khochfar S., 2014, preprint (arXiv:1407.4115)
  4. Agarwal B., Khochfar S., Johnson J. L., Neistein E., Dalla Vecchia C., Livio M., 2012, MNRAS, 425, 2854
  5. Alvarez M. A., Wise J. H., Abel T., 2009, ApJL, 701, L133
  6. Begelman M. C., 2010, MNRAS, 402, 673
  7. Begelman M. C., Rossi E. M., Armitage P. J., 2008, MNRAS, 387, 1649
  8. Begelman M. C., Shlosman I., 2009, ApJL, 702, L5
  9. Begelman M. C., Volonteri M., Rees M. J., 2006, MNRAS, 370, 289
  10. Bonnor W. B., 1956, MNRAS, 116, 351
  11. Bromm V., Loeb A., 2003, ApJ, 596, 34
  12. Bromm V., Yoshida N., 2011, ARA&A, 49, 373
  13. Chen K.-J., Heger A., Woosley S., Almgren A., Whalen D. J., Johnson J. L., 2014, ApJ, 790, 162
  14. Choi J.-H., Shlosman I., Begelman M. C., 2013, ApJ, 774, 149
  15. Clark P. C., Glover S. C. O., Klessen R. S., 2008, ApJ, 672, 757
  16. Clark P. C., Glover S. C. O., Smith R. J., Greif T. H., Klessen R. S., Bromm V., 2011, Science, 331, 1040
  17. Dijkstra M., Ferrara A., Mesinger A., 2014, MNRAS, 442, 2036
  18. Dijkstra M., Haiman Z., Mesinger A., Wyithe J. S. B., 2008, MNRAS, 391, 1961
  19. Ebert R., 1955, Z. Astrophys., 37, 217
  20. Fan X. et al., 2006, AJ, 132, 117
  21. Fan X. et al., 2003, AJ, 125, 1649
  22. Federrath C., Sur S., Schleicher D. R. G., Banerjee R., Klessen R. S., 2011, ApJ, 731, 62
  23. Ferrarese L., Merritt D., 2000, ApJL, 539, L9
  24. Gebhardt K. et al., 2000, ApJL, 539, L13
  25. Greene J. E., 2012, Nature Communications, 3, 1304
  26. Greif T. H., 2014, MNRAS, 444, 1566
  27. Greif T. H., Springel V., White S. D. M., Glover S. C. O., Clark P. C., Smith R. J., Klessen R. S., Bromm V., 2011, ApJ, 737, 75
  28. Haiman Z., 2006, New Astron. Rev., 50, 672
  29. Haiman Z., 2009, Observing the First Stars and Black Holes, Thronson H. A., Stiavelli M., Tielens A., eds., p. 385
  30. Heger A., Fryer C. L., Woosley S. E., Langer N., Hartmann D. H., 2003, ApJ, 591, 288
  31. Hosokawa T., Omukai K., Yorke H. W., 2012, ApJ, 756, 93
  32. Hosokawa T., Yorke H. W., Inayoshi K., Omukai K., Yoshida N., 2013, ApJ, 778, 178
  33. Inayoshi K., Hosokawa T., Omukai K., 2013, MNRAS, 431, 3036
  34. Inayoshi K., Omukai K., Tasker E. J., 2014, MNRAS, 445, L109
  35. Johnson J. L., Bromm V., 2007, MNRAS, 374, 1557
  36. Johnson J. L., Khochfar S., Greif T. H., Durier F., 2011, MNRAS, 410, 919
  37. Johnson J. L., Whalen D. J., Fryer C. L., Li H., 2012, ApJ, 750, 66
  38. Johnson J. L., Whalen D. J., Li H., Holz D. E., 2013, ApJ, 771, 116
  39. Komatsu E. et al., 2009, ApJS, 180, 330
  40. Koushiappas S. M., Bullock J. S., Dekel A., 2004, MNRAS, 354, 292
  41. Larson R. B., 1969, MNRAS, 145, 271
  42. Latif M. A., Bovino S., Grassi T., Schleicher D. R. G., Spaans M., 2014a, preprint (arXiv:1408.3061)
  43. Latif M. A., Schleicher D. R. G., Bovino S., Grassi T., Spaans M., 2014b, ApJ, 792, 78
  44. Latif M. A., Schleicher D. R. G., Schmidt W., 2014, MNRAS, 440, 1551
  45. Latif M. A., Schleicher D. R. G., Schmidt W., Niemeyer J., 2013a, MNRAS, 433, 1607
  46. Latif M. A., Schleicher D. R. G., Schmidt W., Niemeyer J., 2013b, MNRAS, 430, 588
  47. Latif M. A., Schleicher D. R. G., Schmidt W., Niemeyer J., 2013c, MNRAS, 432, 668
  48. Li Y. et al., 2007, ApJ, 665, 187
  49. Lodato G., Natarajan P., 2006, MNRAS, 371, 1813
  50. Madau P., Rees M. J., 2001, ApJL, 551, L27
  51. Montero P. J., Janka H.-T., Müller E., 2012, ApJ, 749, 37
  52. Oh S. P., Haiman Z., 2002, ApJ, 569, 558
  53. Omukai K., 2001, ApJ, 546, 635
  54. Omukai K., Nishi R., 1998, ApJ, 508, 141
  55. Pakmor R., Bauer A., Springel V., 2011, MNRAS, 418, 1392
  56. Penston M. V., 1969, MNRAS, 144, 425
  57. Peters T., Schleicher D. R. G., Klessen R. S., Banerjee R., Federrath C., Smith R. J., Sur S., 2012, ApJL, 760, L28
  58. Peters T., Schleicher D. R. G., Smith R. J., Schmidt W., Klessen R. S., 2014, MNRAS, 442, 3112
  59. Prieto J., Jimenez R., Haiman Z., 2013, MNRAS, 436, 2301
  60. Regan J. A., Haehnelt M. G., 2009, MNRAS, 396, 343
  61. Regan J. A., Johansson P. H., Haehnelt M. G., 2014, MNRAS, 439, 1160
  62. Regan J. A., Johansson P. H., Wise J. H., 2014, ApJ, 795, 137
  63. Ripamonti E., Abel T., 2004, MNRAS, 348, 1019
  64. Schleicher D. R. G., Palla F., Ferrara A., Galli D., Latif M., 2013, A&A, 558, A59
  65. Schleicher D. R. G., Spaans M., Glover S. C. O., 2010, ApJL, 712, L69
  66. Schober J., Schleicher D., Federrath C., Glover S., Klessen R. S., Banerjee R., 2012, ApJ, 754, 99
  67. Shang C., Bryan G. L., Haiman Z., 2010, MNRAS, 402, 1249
  68. Smith R. J., Glover S. C. O., Clark P. C., Greif T., Klessen R. S., 2011, MNRAS, 414, 3633
  69. Spaans M., Silk J., 2006, ApJ, 652, 902
  70. Springel V., 2010, MNRAS, 401, 791
  71. Stacy A., Greif T. H., Klessen R. S., Bromm V., Loeb A., 2013, MNRAS, 431, 1470
  72. Sugimura K., Omukai K., Inoue A. K., 2014, MNRAS, 445, 544
  73. Sur S., Schleicher D. R. G., Banerjee R., Federrath C., Klessen R. S., 2010, ApJL, 721, L134
  74. Toomre A., 1964, ApJ, 139, 1217
  75. Truelove J. K., Klein R. I., McKee C. F., Holliman, II J. H., Howell L. H., Greenough J. A., 1997, ApJL, 489, L179
  76. Turk M. J., Norman M. L., Abel T., 2010, ApJL, 725, L140
  77. Turk M. J., Oishi J. S., Abel T., Bryan G. L., 2012, ApJ, 745, 154
  78. Van Borm C., Spaans M., 2013, A&A, 553, L9
  79. Visbal E., Haiman Z., Bryan G. L., 2014, MNRAS, 442, L100
  80. Volonteri M., 2012, Science, 337, 544
  81. Volonteri M., Begelman M. C., 2010, MNRAS, 409, 1022
  82. Volonteri M., Bellovary J., 2012, Reports on Progress in Physics, 75, 124901
  83. Volonteri M., Rees M. J., 2005, ApJ, 633, 624
  84. Wang H.-H., Klessen R. S., Dullemond C. P., van den Bosch F. C., Fuchs B., 2010, MNRAS, 407, 705
  85. Wise J. H., Abel T., 2007, ApJ, 665, 899
  86. Wise J. H., Turk M. J., Abel T., 2008, ApJ, 682, 745
  87. Wolcott-Green J., Haiman Z., 2012, MNRAS, 425, L51
  88. Wolcott-Green J., Haiman Z., Bryan G. L., 2011, MNRAS, 418, 838
  89. Xu H., O’Shea B. W., Collins D. C., Norman M. L., Li H., Li S., 2008, ApJL, 688, L57
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
""
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
   
Add comment
Cancel
Loading ...
105170
This is a comment super asjknd jkasnjk adsnkj
Upvote
Downvote
""
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters
Submit
Cancel

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test
Test description