Supermassive Black Hole Formation at High Redshifts via Direct Collapse:
Physical Processes in the Early Stage
We use numerical simulations to explore whether direct collapse can lead to the formation of supermassive black hole (SMBH) seeds at high redshifts. Using the adaptive mesh refinement code ENZO, we follow the evolution of gas within slowly tumbling dark matter (DM) halos of and kpc. For our idealized simulations, we adopt cosmologically motivated DM and baryon density profiles and angular momentum distributions. Our principal goal is to understand how the collapsing flow overcomes the centrifugal barrier and whether it is subject to fragmentation which can potentially lead to star formation, decreasing the seed SMBH mass. We find that the collapse proceeds from inside out and leads either to a central runaway or to off-center fragmentation. A disk-like configuration is formed inside the centrifugal barrier, growing via accretion. For models with a more cuspy DM distribution, the gas collapses more and experiences a bar-like perturbation and a central runaway on scales of pc. We have followed this inflow down to pc ( AU), where it is estimated to become optically thick. The flow remains isothermal and the specific angular momentum, , is efficiently transferred by gravitational torques in a cascade of nested bars. This cascade is triggered by finite perturbations from the large-scale mass distribution and by gas self-gravity, and supports a self-similar, disk-like collapse where the axial ratios remain constant. The mass accretion rate shows a global minimum on scales of pc at the time of the central runaway. In the collapsing phase, virial supersonic turbulence develops and fragmentation is damped. Models with progressively larger initial DM cores evolve similarly, but the timescales become longer. In models with more organized initial rotation — when the rotation of spherical shells is constrained to be coplanar — a torus forms on scales pc outside the disk, and appears to be supported by turbulent motions driven by accretion from the outside. The overall evolution of the models depends on the competition between two timescales, corresponding to the onset of the central runaway and of off-center fragmentation. In models with less organized rotation — when the rotation of spherical shells is randomized (but the total angular momentum remains unchanged) — the torus is greatly weakened, the central accretion timescale is shortened, and off-center fragmentation is suppressed — triggering the central runaway even in previously ‘stable’ models. The resulting seed SMBH masses is found in the range , substantially higher than the mass range of Population III remnants. We argue that the above upper limit on appears to be more realistic, and lies close to the cutoff mass of detected SMBHs. Corollaries of this model include a possible correlation between SMBH and DM halo masses, and similarity between the SMBH and halo mass functions, at time of formation.
Subject headings:methods: numerical — galaxies: formation — galaxies: high-redshift — cosmology: theory — cosmology: dark ages, reionization, first stars
Observational evidence points to most large galaxies hosting supermassive black holes (SMBHs) in their centers (e.g., Ferrarese, 2005) and to the possible coevolution of galaxies and SMBHs, whether causal (e.g., Gebhardt et al., 2000; Ferrarese & Merritt, 2000; Tremaine et al., 2002) or not (Jahnke & Maccio, 2011). A growing number of QSOs have been found at , including a bright QSO at (e.g., Fan et al., 2003; Mortlock et al., 2011). Hence, some SMBHs with few had already formed before the universe was Myr old, raising the issue of how such large black holes (BHs) could have grown so quickly.
The early formation and growth of SMBHs must be understood in a cosmological context. Early SMBHs could have formed from remnants of the first generation of stars, Population III, which subsequently grew by gas accretion and mergers (e.g., Haiman & Loeb, 2001; Yoo & Miralda-Escudé, 2004; Volonteri & Rees, 2006; Li et al., 2007; Pelupessy et al., 2007; Tanaka & Haiman, 2009). Early studies of the initial mass function (IMF) of these objects argued in favor of very massive objects (e.g., Abel et al., 2002; Bromm et al., 2002; O’Shea & Norman, 2007), but more recent work indicates a rather normal IMF (e.g., Turk et al., 2009; Clark et al., 2011; Hosokawa et al., 2011; Stacy et al., 2012; Wise et al., 2012). Even if there were a sufficient supply of remnant BHs which grew via gas accretion at the Eddington rate with 10% radiation efficiency, it would have taken them at least yr to reach (Salpeter, 1964). This timescale barely matches the available time required to make bright QSOs at .
Such efficient growth by accretion is unlikely, for two reasons. First, Population III stars in the vicinity of the seed, and the BH progenitor itself, ionize the surrounding matter, causing expansion of the H II region and likely suppression of subsequent accretion onto the SMBH (e.g., Safranek-Shrader et al., 2012). Second, the possible supernova explosion of a Pop III star can evacuate the ambient gas in the vicinity of the seed. This suggests that mergers must play an essential role in growing SMBHs from Population III seeds, but this is difficult as well. According to the cold dark matter (CDM) cosmology, dark matter (DM) halos frequently merge in hierarchical structure formation. However, too frequent halo mergers would result in the formation of small -body systems of BH seeds, which would be prone to SMBH slingshot ejection(s) from the DM halo. Given the relatively short time available for growing SMBH seeds to quasar masses, the long list of delaying processes can seriously hurt the Pop III seed scenario.
An appealing alternative is that early SMBHs formed via direct collapse of the halo gas at (e.g., Oh & Haiman, 2002; Bromm & Loeb, 2003; Volonteri & Rees, 2005; Begelman et al., 2006; Wise et al., 2008; Levine et al., 2008; Regan & Haehnelt, 2009; Begelman & Shlosman, 2009; Mayer et al., 2010; Schleicher et al., 2010; Hosokawa et al., 2011; Johnson et al., 2011; Prieto et al., 2013). The direct collapse paradigm assumes that an SMBH seed, , forms from the cold halo gas. This process favors a high density, high inflow rate environment in which star formation is suppressed. A massive central object — a supermassive star (SMS), forms at the center of the collapsing region. This object is powered by a combination of core nuclear burning and Kelvin-Helmholtz contraction (e.g., Begelman et al., 2006, 2008; Begelman, 2010). After the stellar core collapses and forms the SMBH seed, its convective envelope is powered by the accretion onto this seed — this configuration is termed a quasistar. Consequently, the seed BH, at the beginning, can grow to in less than a few Myrs.
In this context, constraints on the accretion rate onto the SMBH become weaker, and additional problems anticipated in the Pop III seed scenario do not play a major role. In order to make this scenario work, the gas needs to inflow rapidly from kpc to 100 AU scales. The most probable site for such runaway gas collapse is expected in DM halos with virial temperatures K. If the halo gas can cool only via atomic cooling, the necessary condition for the collapse to be triggered is for the gas temperature to be lower than the virial temperature of the host DM halo.
However, two caveats to the direct collapse scenario must be addressed (e.g., Begelman & Shlosman, 2009). First, the angular momentum barrier, in principle, can terminate the collapse well before it reaches AU scales. Owing to angular momentum conservation, collapse is halted when the rotational velocity reaches the circular velocity. In order to overcome this barrier, there should be a physical process to redistribute the gas angular momentum outward. Second, gas could fragment during the collapse, depleting the accretion stream by forming gas clumps and ultimately stars. The gas clumps could also disturb the accretion pattern.
In this paper, we study the physical processes that accompany the initial and intermediate stages of gas collapse in a DM halo. We focus on the dynamical processes: spontaneous and induced symmetry breaking, nested gaseous bar formation, and the role of supersonic turbulence. Section 2 discusses the relevant processes. In Section 3 we describe the numerical setups for a set of simulations intended to quantify these processes, and in Section 4 we present our results. We summarize and discuss our results in Section 5.
2.1. Angular Momentum Transfer
It is generally understood that local viscous transport mechanisms are inefficient on galactic scales. Even on scales of a few parsecs, their characteristic timescales are prohibitively long and angular momentum transport therefore requires nonlocal mechanisms, e.g., large-scale magnetic fields (e.g., Blandford & Payne, 1982), or by gravitational torques. Turbulent motions, e.g., driven by the gravitational collapse (e.g., Hoyle, 1953; Klessen & Hennebelle, 2010; Sur et al., 2010; Vazquez-Semadeni, 2010) and other mechanisms, can amplify the seed magnetic fields (e.g., Balbus & Hawley, 1998; Subramanian, 2010, for a review). On larger spatial scales, the angular momentum redistribution most likely depends on long-range gravitational torques (e.g., Lynden-Bell & Kalnajs, 1972; Shlosman et al., 1989, 1990). When angular momentum is conserved, gas collapse is quickly stopped by the centrifugal barrier. Prior to collapse within DM halos, the gas can exhibit the same angular momentum distribution as the DM. The typical spin acquired by a halo during maximal expansion is given by a dimensionless parameter , where is the halo angular momentum, and and are its virial mass and the circular velocity at the virial radius (e.g., Peebles, 1969; Barnes & Efstathiou, 1987; Bullock et al., 2001). The mean spin of DM halos is . If the gravitational potential of the gas dominates, this will bring the gas to a centrifugal barrier when it has collapsed by a factor of . When the DM potential dominates, the collapse will stop after a decade in radius (e.g., Mo et al., 1998; Shlosman, 2013, and refs. therein).
Inflow beyond the centrifugal barrier requires substantial and continuous angular momentum loss by the gas. What are the options? If the gas disk becomes nonaxisymmetric, gravitational torques between the inner and outer gas, as well as between the gas and DM, can drain angular momentum away from the collapsing gas. The lowest nonaxisymmetric modes, , prevail in both the gas and the DM, because their rise times are shortest. The mode requires perturbation of the center of mass of the gaseous disk — this is not always possible. The next fastest growing mode is — the bar-like mode. Typically it will have the highest amplitude, in the absence of . For the inflow to extend to smaller spatial scales, the gas disk must be self-gravitating. Under these conditions, a cascade of bars can be maintained over a substantial dynamic range in radius. Contraction of the gas via such an avalanche can be related to self-organized criticality (Lichtenberg & Lieberman, 1995). An analogy can be made between the chaos driven when the bar exceeds a certain strength and a sandpile whose slope is increased until the sand slides off (Shlosman, 2005).
The bar cascade can be triggered in a number of ways. First, if the disk becomes self-gravitating and cold, its axial symmetry is known to be broken spontaneously when the ratio of bulk kinetic energy to absolute potential energy, , is larger than a certain critical value. The gas bar exerts torque on the gas disk and transfers angular momentum outward (e.g., Shlosman et al., 1989, 1990; Englmaier & Shlosman, 2004). This process allows the collapsing gas to cross the centrifugal barrier and to continue the collapse to smaller scales, if angular momentum is continuously extracted. The threshold value for gaseous bar formation depends also on the shape of the gas distribution and is (Christodoulou et al., 1995), where for disks and for spheres. This instability is spontaneous and does not require any trigger — it will grow exponentially from any infinitesimal perturbation.
Several additional effects can drive a cascade through finite amplitude perturbations. First, the shape and dynamics of a self-gravitating gas disk will respond to the potential of the DM halo in which it is embedded. DM halos are known to be triaxial (e.g., Allgood et al., 2006; Berentzen et al., 2006; Berentzen & Shlosman, 2006; Heller et al., 2007), and therefore exert gravitational torques on the disk, which will respond with low Fourier modes. Second, finite bar perturbations can also be driven by tidal effects of massive DM and baryon clumps, i.e., halo substructure (e.g., Romano-Diaz et al., 2008) and from mergers. In addition, asymmetry in the background gravitational potential can be generated by the spontaneous break of axial symmetry in the collapsing flow, DM, gas or both. It is well-known that nonradial perturbations can grow in supersonic accretion flows onto compact objects (e.g., Hunter, 1962; Garlick, 1979; Moncrief, 1980; Goldreich et al., 1996; Lai & Goldreich, 2000), without any assistance from self-gravity. Gas collapse inside an axisymmetric DM halo is similar — in Section 4.1 of the present work, we show the appearance of nonradial perturbations in the equatorial plane of the disk, although we start from idealized conditions of an isolated halo and embedded baryons which possess axial symmetry. In a comprehensive cosmological simulation, all these factors would provide finite amplitude perturbations that do not require the -parameter to cross the threshold.
At the spatial scale on which the gas disk forms, the potential contribution from the gas is already comparable with that of the DM, and the gas cooling time is shorter than the dynamical time. These conditions imply that the collapsing gas becomes self-gravitating and nearly isothermal, and the density profile evolves toward a singular isothermal sphere111Throughout this paper, we use for spherical radii, and for cylindrical radii. . Before the gas reaches the centrifugal barrier, we can assume it flows in with the free-fall speed, whose dependence on radius is weaker than logarithmic. Away from the disk and the centrifugal barrier, the gravitational potential is still dominated by the DM, so the accretion rate can be estimated at
where for a DM halo, and we used the universal baryon fraction for the gas-to-total mass ratio. At these radii, only the gas participates in the collapse, as the DM is in equilibrium due to the initial conditions.
On the other hand, when the gas dominates the potential, the mass accretion rate is of order , where is the gas sound speed and is the viscosity parameter in the disk (e.g. Shu, 1977; Shlosman & Begelman, 1989). As we shall see in Section 4.1, this situation in our simulations occurs inside the centrifugal barrier. For local viscosity in the disk, , but for nonlocal viscosities, e.g., magnetic or gravitational torques, , especially in the self-gravitating case where the effective (Shlosman et al., 1990). In this latter case, the inflow rate is still given by Equation 1, where is the virial velocity of the gravitational potential which is now dominated by the gas, i.e., in Eq. 1. The inflow velocity will depend on the compactness of the gas distribution and can substantially exceed the virial speed of the DM halo.
A high accretion rate is crucial for the formation of a SMS or disk that can develop into an SMBH seed. For example, Begelman (2010) pointed out that an infall rate exceeding yr is necessary in order to form a SMS, which has a nuclear burning timescale of only a few million years. The initial SMBH seed of , formed by core collapse of such a star, then grows at a highly super-Eddington rate inside the remaining convective envelope (the ‘quasistar’ phase: Begelman et al. (2008)), reaching or more in less than a few Myr.
2.2. Fragmentation of Accretion Flows
Gas fragmentation can terminate the collapse by consuming the gas supply. In addition, the clumps formed during fragmentation can excite odd-mode perturbations in the gas disk, damping the bar instability which plays the key role in angular momentum redistribution. In order to continue collapse to very small spatial scales, fragmentation should be suppressed.
There are several ways to suppress gas fragmentation. Supersonic turbulence can suppress fragmentation via shocks that sweep the density fluctuations in a crossing time (e.g., Padoan, 1995; Padoan & Nordlund, 2002; Krumholz & McKee, 2005). Studies of supersonic turbulence find that the probability distribution function (PDF) of the density fluctuations is lognormal (i.e., Gaussian in the log) under a wide range of conditions (e.g., Vazquez-Semadeni, 1994; Padoan, 1995; Scalo et al., 1998; Ostriker et al., 1999; Padoan & Nordlund, 2002; Krumholz & McKee, 2005):
where the distribution mean , ( being mean density), and its dispersion is ( is the flow Mach number). This distribution is the result of being a random variable, which is itself a product of independent random variables. Note that the coefficient in the definition of is not a constant and depends on the driving mode of the turbulence. Federrath et al. (2008) found that varies between 1/3 to unity, for solenoidal and compressive driving, correspondingly.
The lognormal density PDF has been observed in simulations of turbulent motions on the scales of giant molecular clouds, where the background (unperturbed) density is uniform. Here, we shall test the development of turbulence in collapsing flows with steep preexisting density gradients.
The critical overdensity required for fragmentation is , where is the Jeans length at the mean density and is the turbulent length scale at which the velocity dispersion reaches the sound speed (Padoan & Nordlund, 2002; Krumholz & McKee, 2005; Begelman & Shlosman, 2009). In a supersonic turbulent medium, overdensities higher than are unstable to fragmentation. Estimating the mass fraction in Jeans-unstable fragments that collapse on a timescale shorter than the free-fall timescale as , we find that for . If the collapsing gas maintains supersonic turbulence, the amount of fragmenting gas is negligible.
In our models, we have assumed optically thin flows. The spherically-symmetric, isothermal density distribution of baryons modeled in this paper becomes opaque at pc, where is the ionization fraction. For temperatures encountered in this simulation, and the electron scattering optical depth is negligible, as is the free-free absorption. For primordial chemical composition and the range of temperatures K, the opacity is expected to be dominated by bound-bound transitions in hydrogen, and especially by the opacity for Lyman photons.
Radiation fields can suppress gas fragmentation by dissociating molecular gas or preventing its formation in the first place. If H is absent, the gas temperature will be maintained at K (e.g., Omukai, 2001; Shang et al., 2010). The prime agent of H dissociation can be an externally-produced Lyman-Werner ( eV eV) continuum, e.g., from neighboring Pop III stars. (e.g., Dijkstra et al., 2008). H formation could be inhibited by the high temperatures likely to obtain in the collapsing gas. The cooling efficiency of Ly produced in the accretion flow at is limited by the high optical depths. For a sufficiently high column density of neutral hydrogen, collisional de-excitation can decrease the escape fraction of Ly, preventing the gas from cooling below K and suppressing H formation thermally (e.g., Spaans & Silk, 2006; Schleicher et al., 2010; Latif et al., 2011).
One might speculate whether Ly photons generated within the accretion flow could be upscattered into the Lyman-Werner continuum before escaping, thus producing a dissociating flux self-consistently. For an isothermal density profile, we estimate the required Ly optical depth (measured in the Doppler core) of (e.g., Harrington, 1973; Neufeld, 1990), which can be attained at radii pc for the parameters of our models. Other processes, however, are likely to allow the energy in these upscattered Ly photons to leak away before they can escape. For one thing, these photons can be absorbed by even small amounts of dust. Extrapolating Figure 18 of Neufeld (1990) to our column densities (which can be inferred from the next section), we find that a metallicity as low as solar will do the job. An additional process that can destroy the Ly photons is the resonant pumping of H, followed by fluorescent decay to vibrational levels.
We note that fragmentation might be avoidable in the presence of supersonic turbulence, even if H is present (e.g., Begelman & Shlosman, 2009). This happens because the fraction of fragmenting gas decreases with .
In our simulations, we assume that H is destroyed and do not analyze its contribution to the gas cooling. We also neglect magnetic fields and their effects on the turbulent flow, and model gravitational collapse within the DM halo hydrodynamically. This is justified because -fields are expected to be weak at high , and because of the high temperature of the flow, K. At these high temperatures, the value of the critical -field strength needed to limit the compression in supersonic isothermal shocks is also higher (Padoan et al., 2007).
3. Numerical technique
3.1. The numerical code ENZO
In order to test the theoretical arguments made in Section 2, we use an Eulerian adaptive mesh refinement (AMR) code ENZO-2.1, which has been tested extensively and is publicly available (Bryan & Norman, 1997; Norman & Bryan, 1999). ENZO uses a particle-mesh -body method to calculate the gravitational dynamics including collisionless DM particles, and a second-order piecewise parabolic method (PPM, Bryan et al., 1995) to solve hydrodynamics. The structured AMR used in ENZO places no fundamental restrictions on the number of rectangular grids used to cover some region of space at a given level of refinement, or on the number of levels of refinement (Berger & Colella, 1989). A region of the simulation grid is refined by a factor of 2 in length scale if the gas density is greater than , where is the minimum density above which refinement occurs, is the refinement factor and is the maximal AMR refinement level. This refinement corresponds to a spatial resolution of AU.
The Truelove et al. (1997) requirement for resolution of the Jeans length, i.e., at least 4 cells, has been verified. However, the actual resolution of the Jeans length in our simulations exceeds this criterion, depending on the distance from the center, because of the baryon density dependence. Specifically, throughout most of the spherical collapse region, pc, the Jeans length is resolved with cells. Between 1 – 10 pc, it is resolved by cells, for 0.001 – 1 pc by 10-30 cells, and inside 0.001 pc by 4 cells. This estimate is based on the baryon properties only, and ignores the effect of DM. If the latter is taken into account, the Jeans length increases by a factor of a few. Therefore, our resolution of the Jeans length increases correspondingly.
We have also run test models resolving the Jeans length with 64 cells, as required by e.g., Sur et al. (2010), Federrath et al. (2011), Turk et al. (2012) and Latif et al. (2013) in the MHD simulations. No qualitative diference in the evolution has been found, except a slight, , increase in the onset time of the 2nd stage of the gravitational collapse, i.e., of the central runaway (see Section 4 for the definition).
ENZO follows the non-equilibrium evolution of six species: , and (Abel et al., 1997; Anninos et al., 1997) in a gas with primordial composition. It calculates radiative heating and cooling following atomic line excitation, recombination, collisional excitation and free-free transitions. Radiative losses from atomic cooling are computed in the optically-thin limit. As we discussed in Section 2, there are several radiation transfer processes that have been suggested to prevent formation. In order to include these effects without implementing a full radiative transfer calculation, we exclude the chemistry and cooling related to in this paper.
3.2. Simulation Setups
For a gas with primordial composition and no , the earliest collapse can occur in DM halos whose virial temperatures exceed . For this reason we set up a spherical DM halo with at . According to the top-hat model, such a halo will have a virial radius of and a virial temperature , according to the WMAP5 cosmology (, , and ) (Komatsu et al., 2009). We simulate this DM halo within a computational domain.
For the DM halo, we assume an isothermal sphere. This halo model is similar to the universal DM halo model (NFW, Navarro et al., 1996, 1997) at the spatial scales of interest. The isothermal model also provides a simple analytical form for the velocity distribution function which allows us to generate a live isotropic DM particle distribution that maintains a stable equilibrium with particles. We introduce a DM core radius, , within which the DM density is approximately constant (e.g., El-Zant et al., 2001; Tonini et al., 2006; Romano-Diaz et al., 2008; Primack, 2009). Therefore, plays the role of the King radius for a nonsingular isothermal sphere. We vary to see the effect of DM halo structure on the development of gas collapse. In particular, we implement a number of simulations, models A – E, with different DM density core sizes, as given in Table 1.
The DM halo and the gas have been laid down at virial and pressure (for the gas) equilibrium, so no transient adjustment occurs. The gas halo has a constant density core at the start of the simulations. The gas, however, instantly cools down from the virial temperature to about 10,000 K, and so its pressure support is negligible. The total gas mass in the halo is estimated from the cosmological baryon fraction, . The angular momentum of the halo gas is computed from , defined in Section 2. In large-scale -body cosmological simulations, Bullock et al. (2001) has found that the specific angular momentum in DM halos roughly follows , close to a flat rotation curve. We therefore impose a constant rotation velocity on the gas outside the core, assuming that the rotational axis runs parallel to the -axis. For , we assume solid body rotation for the gas, i.e., (e.g., Samland & Gerhard, 2003; Heller et al., 2007). The outer boundary of the model we smooth the sharp density cutoff at the virial radius by adopting the gas density profile beyond .
We also ran a number of models with modified initial conditions. One such representative model is discussed as model Dmod, i.e., modified model D, in Section 4.4. In this model we randomized the orientation of the initial angular momentum in spherical shells, keeping the mean-square angular momentum unchanged. For this to happen, we have slightly decreased of the inner shells and compensated the outer shells by a slight increase in . This was done in order to mimic cosmological initial conditions, based on the simulations by Romano-Diaz et al. (2009), in particular on their Figure 19.
To verify that the DM halo is indeed formed in equilibrium, we tested DM-only models and confirmed their stability, especially in the central region. Additional tests have been run for model B with an isothermal equation of state for the gas; these show very similar evolution. We also tested a model with an adiabatic equation of state — the gravitational collapse did not proceed far in this model, as expected.
We would like to emphasize several important points about the gravitational collapse modeled here, before we show the results of numerical simulations. First, the evolution described here truly proceeds from inside out. Development of the central mass accumulation and, therefore, of the prospective seed SMBH, happens on timescales much shorter, Myr, than the free-fall timescale for the DM halo, Myr. Hence, the dynamics of the outer parts on scales of kpc, although interesting in itself, appears to be irrelevant for the central regions of pc, which dominate the formation of the SMBH seed. Next, the gravitational collapse proceeds in two stages. The first stage involves infall, braking and the formation of a gaseous disk-like configuration inside the centrifugal barrier. All the models exhibit this phase. Models A – C show the second stage of gravitational collapse in which the inner part of the disk, pc, develops a runaway which proceeds to the point where the numerical evolution has been terminated, at pc AU. Our task will be to explain this evolution. We shall emphasize models B and D as representative of 2-stage and 1-stage collapse, respectively. Model Dmod is discussed separately.
4.1. Loss of axial symmetry and central runaway
As the gas has little rotational support, it goes into nearly free-fall collapse, which develops from inside out because of the shorter dynamical timescales at small and larger gravitational accelerations there. This leads to the establishment of density profile and to a largely homologous collapse, except when and where the angular momentum becomes important. There is little difference among models at this stage, as they differ only in the value of the DM core radius , and, therefore, in the depth of the DM potential well. Figure 1 shows face-on density-weighted projections of the density222, where the weight function is of the gas disk at Myr in models B and D, with pc and 1.5 pc, respectively. Unless mentioned, other models behave similarly.
During the first stage of collapse, the baryon angular momentum flows inward. The inner gas is the first to reach the centrifugal barrier, and develops a disk-like configuration at pc. The disk grows quickly in radius and mass, and by Myr, the snapshots display well-developed gas disks with radii pc. The disk boundaries, both in and in , are delineated by standing shocks, discussed later. The surface densities of all models are well-approximated by a power-law with . Figure 2 displays the evolution of the surface density for model B, as well as comparisons of the surface densities at Myr among models B, D and E. The disks are truncated at pc, as clearly seen in this figure. The outer surface density profiles of all disks, in models A to E, are very similar, but in the central pc the density profiles differ: the density in model B within the central 0.1 pc is more than an order of magnitude higher than the density in model D.
The DM density profiles show little evolution even in the central regions. In models B, D and E (Figure 2), and also in all other models, we observe the trend that more cuspy DM halos develop higher surface and volume density gas disks. The disk in D has a surface density about a factor of times higher than E. The disk volume density also appears higher in more cuspy DM halos (Figure 3). This difference is clearly observed within the DM core radius of less cuspy halos.
Over the next time period of a few yr, the disk surface density grows as a result of accretion. The second stage of the collapse is reached only in models A – C, and is characterized by both baryonic angular momentum outflow and continuing mass influx. On larger scales, just outside the centrifugal barrier, increases, moving the barrier outward. Inside the centrifugal barrier, the disk develops a global instability, dominated by the Fourier component — a bar-like mode, which transfers to the outer gas and to the DM. For example, model B exhibits a well defined gaseous bar in the central pc of the disk at Myr. The appearance of this bar-like mode is characteristic also of smaller spatial scales of this model (Fig. 4), as we shall analyze below. The timescale for the onset of this stage of the central runaway is given in Table 1, and is increasing along the sequence A C, i.e., toward less cuspy halos. On the contrary, models D and E display a weak central depression which has the appearance of a ring-like structure around the center and no axial symmetry-breaking. At a later stage, these models exhibit off-center fragmentation.
Figure 1a confirms the crucial role of gravitational torques in shaping the inner mass distribution via angular momentum transfer, which requires a substantial departure from axial symmetry. This symmetry breaking gives rise to the lowest modes, , 2, and occasionally . Only exhibits a mild displacement of the center of mass of the disk, while has the appearance of a strong gaseous bar which drives spiral structure, and shows up as tri-armed spirals. We follow the runaway collapse to a spatial scale of pc (20 AU), and stop the calculation there because the flow is expected to become opaque and enter a dynamical regime where radiation pressure must be taken into account (see Section 2.2).
Owing to the AMR nature of ENZO, the small scale structure in Figures 2 and 3 is only resolved when a given region exceeds the density threshold for a new refinement. The density profile in model B is resolved to AU at the end of the simulation, while model D gas only reaches a resolution of 0.1 pc. Resolving the AU scale means that the gas inflow has reached this scale by overcoming the angular momentum barrier. It demonstrates that the gravitational torques resulting from axial symmetry-breaking at pc in model B trigger continuous gas inflow on progressively smaller scales. Finally, Figure 3 demonstrates that the timescale of this inflow is very short, from Myr to Myr, a truly runaway collapse with a dynamical time of yrs.
For the idealized case of spherical accretion within an external gravitational potential given by , the gas mass flux is given by with weakly dependent on . Figure 3 confirms that the density profile of the collapsing gas in model B tends to the isothermal density profile, , for pc. It decreases sharply outside this radius due to the shock at the centrifugal barrier, then continues with the same slope at larger radii. (Note that all models start with a flat gas density core of 100 pc.) Figure 5a shows the radial profile of the inflow velocity measured in the disk plane. We can divide the -range into the collapsing gas in the region outside 10 pc, the disk plateau at pc, and the inner region of runaway collapse — all at the time of the central runaway, Myr. The stronger than expected increase of Mach number with decreasing radius has its origin in a combination of temperature decrease (due to an increased gas density) and the steepening increase of toward the center, as the result of the runaway collapse of a finite gas mass.
Thus the outermost region is dominated by free-fall kinematics because the angular momentum is low. The density profile is maintained as if one neglects the weak variation of . The intermediate disky region has and the same slope of log . The innermost collapse exhibits a rapidly increasing and steepening infall velocity (Fig. 5a,b). The accretion flow is mostly rotationally supported in the disky region of pc, as well as in the central region, where the self-gravitating collapse proceeds with a non-negligible angular momentum. The time evolution of in the central region reflects the self-gravitating collapse developing there, and a dynamical decoupling of the gas from the DM background potential.
The measured mass accretion rate, , appears to be approximately flat, , exterior to the radial shock, in agreement with Eq. 1 (Fig. 6). In the disk region, decreases abruptly to at Myr, and drops by an additional order of magnitude during the second stage of the collapse, when the material is drained by the central runaway. At this latter time, Myr, the central runaway is accompanied by peak values of . The decrease in the accretion rate in the range of pc is a reflection of mass conservation, as the bar-like mode channels gas inward at a higher rate than it can be resupplied by disk accretion, and thus drains this region. This minimum in is found in the region where the turbulence induced by the radial shock decays. In Section 5.1, we use this radius in model B to estimate the baryon mass that participates in the second stage of the collapse, and therefore, the mass of the seed SMBH.
The high mass accretion rate is an important ingredient of early SMBH formation in the direct collapse model (e.g., Begelman et al., 2006, 2008; Begelman & Shlosman, 2009; Begelman, 2010). The central runaway region exceeds the mass accretion rate given in Eq. 1 by about an order of magnitude. The reason for this lies in the fact that during the second stage of the collapse, the inflow velocities are determined by the compactness of the gas distribution (Section 2.1) and not by the DM halo virial velocity. Under these conditions, the ratio in Eq. 1, and the free-fall velocity, , becomes a rapidly growing function of time. As a result, the characteristic -profile of radial infall velocity during this stage is characteristic of the ‘avalanche’ behavior of a dynamically decoupled gas. This radial velocity is expected to grow sharply with time. The timescale of this process depends only on the mass involved in the collapse, and on the ability of the gas to cool ahead of the free-fall time. If the latter condition is fulfilled, the gas will collapse to an infinite density in a finite time.
We note that the gas joining the disk experiences strong radial and vertical (-axis) shocks, as is seen in Figure 5a (and Figs. 7 and 13). But there is no associated shock in the azimuthal motion. The face-on and edge-on velocity fields on scales of pc are shown in Figure 13. While the face-on view is dominated by rotation, the edge-on view exhibits turbulent motions dominated by vortices, as we discuss in the next section.
We next examine the geometry of the collapsing gas during the runaway stage: is it best approximated as one-dimensional collapse with a spherical symmetry, where the angular momentum is completely unimportant, or does it exhibits a cylindrical or disk-like geometry? To answer this, we compiled the density profiles of the collapsing gas in model B at Myr, along two directions in the equatorial plane, and , and along the -axis (Fig. 7). We verify that the collapse outside the centrifugal barrier proceeds in a spherical fashion, which is understandable because the angular momentum in this region is not important. At Myr, the radial shock in the equatorial plane is positioned at pc (see also Figure 5a).
Inside the radial shock marking the centrifugal barrier, i.e., between pc and 10 pc, the equatorial density increases substantially while the density along the axis is almost constant. At pc, the density along the -axis is times the density measured along the - or -axis at the same radial distance. The reason for this abrupt change in density ratio lies in the relative positions of the radial shock and the surface shock (i.e., the shock in the vertical velocity, at roughly constant ) in the disk. While the centrifugal barrier stops the cold inflow toward the rotation axis at pc, the position of the disk surface shock, i.e., the disk thickness, is determined by the postshock temperature which never exceeds K, as seen in Fig. 9.
A strong surface shock is maintained by the infall along the -axis at pc (Fig. 7). Both radial and surface shocks are slowly moving outward with time. Most interestingly, at the time of the central runaway, the surface shock collapses toward the equatorial plane around . The ‘peanut’ shape of the surface shock at this time can be inferred from Fig. 13. Figure 7 reveals that within the central region the second stage of the collapse proceeds in the self-similar fashion, preserving the disk-like configuration. We have tested this by plotting these distributions at various times. But even as a single snapshot, Figure 7 nevertheless reflects the evolutionary trend because different are associated with different dynamical timescales. Hence the fact that for different means that this ratio also does not evolve in time — a clear sign of self-similarity during the central runaway. Hence the gas configuration remains dynamically cold, sustaining the bar cascade that efficiently transfers angular momentum outward (Fig. 4).
We now return to Figure 4, which displays the non-axisymmetric bar-like perturbations on a number of spatial scales. We first observe this instability developing on scales pc, and then propagating inward. The dynamics of gaseous bars has been analyzed by Englmaier & Shlosman (2004), who focused on the decoupling of a gaseous bar from large-scale stellar bars. This decoupling is associated with the rapidly increasing pattern speed of the gaseous bar and its radial contraction.
We measured the pattern speed of the gaseous bar in model B on scales of 1 pc and and 0.01 pc. The resulting values of the pattern speeds, , confirm that the mode tumbles with substantially different speeds. On pc scale, , while on pc scale, , about 2 orders of magnitude faster.
While these pattern speeds are well-defined, there is no clear ‘material’ separation between the spatial scales, i.e., no separation between bars. This continuity of the flow properties is associated with continuity of pattern speeds — smaller scales correspond to larger pattern speeds, , with . The parameter is determined by the mass distribution inside the collapsing gas at the onset of the central runaway. The gas volume density scales as in the central region (Fig. 7) and dominates over the DM density distribution there, causing the gas to decouple dynamically from the DM background in the region of the central runaway. The pattern speed is given by , where , which explains the instantaneous value of above.
To verify that the bar-like () mode dominates over other modes, e.g., and 3, we have Fourier analyzed the gas response to the asymmetric potential. The dominant mode at early times is the mode — the disk center of mass is initially perturbed by the asymmetry developed in the overall mass distribution (Fig. 8). This mode, however, decays quickly to a negligible amplitude. An additional odd mode, , is also present, although at lower amplitude — it decays similarly. On the other hand, the even mode grows to a substantial amplitude with time. At about Myr, it enters exponential growth — this explains the appearance of the gaseous bar at this stage. Clearly, this bar-like mode has formed early in the evolution, when the gas component does not dominate the potential at any radius, and when the global stability parameter (see Section 2.1). This mode is stimulated by the overall mass distribution. The exponential growth at later time happens exactly when and where the gravitational potential of the gas becomes the dominant one. Thus the bar-like mode appears to be driven initially by the shape of the overall potential, but subsequently runs away due to the self-gravity of the gas.
The temperature variations with density for model B are shown in Fig. 9. At the start of disk formation, a radial shock forms at pc and then gradually propagates out to 10 pc. At this time the central runaway is triggered. As shown by Fig. 9, the posthock gas rapidly cools down to low in the postshock region because of the high density. Within the growing disk, pc, we observe an increasing range in K and . The upper limit to the temperature is dropping with density. The central runaway at pc narrows the -range to mostly K but increases the range in . The maximal postshock in the disk is slightly increasing with time, which reflects the increasing shock-impact velocity of the gas which comes from progressively larger distances.
4.2. Interplay between off-center fragmentation and central runaway
As we noted in Section 4.1, models A – C exhibit a central runaway collapse, while models D and E do not show it — not at Myr nor at anytime later on. We now describe the evolution of the representative model D and analyze the reasons for its differences from model B. We shall discuss additional models as well.
The first stage of collapse in model D proceeds similarly to model B, but leads to a disk within the centrifugal barrier that has a much lower surface density at the center. Figure 2 shows that this disk possesses a density core that grows to pc at Myr — the time of the onset of the central runaway in model B, in sharp contrast with this model. In model E, the density core extends to pc at this time. When the gas distributions in models B and D are compared before the central runaway in model B, at yr, the core density in B is higher by about 2 orders of magnitude than in D (Fig. 3), and by even more compared to model E.
Clearly, a cuspier DM distribution leads to a more centrally concentrated gas distribution. This has a clear dynamical consequence: the dynamical timescale prior to the central runaway is shorter in more cuspy potentials. We observe that, at early times of the simulations, the mode is the dominant mode for all the models. This mode is damped faster in models with cuspier DM distributions, via dynamical friction, presumably of the gas against the DM distributions. The mode, which is responsible for the central runaway collapse, can only grow after has decayed significantly (Fig. 8). Therefore, our models A – E form a sequence along which the dynamical timescale becomes longer. Any dynamical instability will have a tendency to increase its amplitude more slowly along this sequence. This includes the growth of , 2 and 3 modes that we observe, in different combinations, in these models.
Further evolution of model D shows the growth of the disk behind the radial shock. The central density increases by about an order of magnitude, but still falls short of the beginning of the central runaway in model B (Figure 3). This disk growth is related to gas with an increasing angular momentum , and steadily growing free-fall velocity, arriving from larger . After Myr, the turbulent pressure in the shocked gas substantially exceeds the thermal pressure behind the shock. This leads to the formation of a geometrically-thick torus with a growing surface (and volume) density (Fig. 10). Figure 11 displays the formation and evolution of this torus at pc. The density profile at Myr shows that a significant fraction of the collapsing gas has stagnated in the torus. At this time, the gas in this region is mainly in circular motion — the radial motion essentially ceases and the turbulent motions decay. The torus becomes azimuthally inhomogeneous, and increasing density leads initially to mild off-center fragmentation, but most of the fragments are immediately sheared and destroyed (e.g., Fig. 10a). Finally, at Myr, one of the fragments becomes substantially compact and begins gravitational collapse and the simulation is stopped (Fig. 12). Model E is similar to D, with off-center fragmentation in the torus occurring at the same time, Myr.
We have followed the evolution of these fragments in models D and E. It differs substantially from the evolution of the central runaway (i.e., collapse of the central fragment) in models A – C. The typical mass of the fragments is about , but it is not clear if this mass is characteristic of the off-center runaway. Each of the fragments in D and E collapses and exhibits a fission into two fragments. We did not follow their evolution further.
The evolution of model D thus diverges from that of model B. The appearance of fragments perturbs the central region of the gas disk. We observe that these perturbations lead to a loss of symmetry in the disk as well as displacing its center of mass, depending on the number of fragments and their azimuthal distribution. The formation of an even-mode bar is suppressed. Therefore, the effect of fragmentation in the torus is that the bar instability is damped.
This result reveals competition between two timescales — that of the central runaway and that of off-center fragmentation in the models. The development of the bar cascade and of off-center fragmentation are facilitated by different instabilities. Gas inflow, of course, drives both instabilities, but the DM plays an important role by regulating the timescale of the central runaway, as well as the decay timescale for the mode. Gaseous bar formation is a global instability in the gas disk and appears to be triggered by both the gas gravity and the global distribution of the collapsing matter. In more centrally-concentrated disks, this instability happens earlier. On the other hand, the fragmentation in the torus is determined by the local surface density and the local pressure — this is a local instability. If the gas turbulent velocities decay first in the torus (i.e., in the outer disk), the torus fragments first. The competition between these two instabilities, local and global, has been investigated in Begelman & Shlosman (2009).
To summarize, we find that the DM density profile, specifically its cuspiness, plays the role of the discriminator between models that experience the second stage of gravitational collapse past the centrifugal barrier and models that do not exhibit this dynamical stage. This dichotomy is related directly to the surface density profile of the growing disk within the centrifugal barrier. But, as we have shown in Section 4.1, it also depends on the asymmetry of the outer mass distribution that develops in the collapsing matter, triggering the bar cascade. Essentially, this corresponds to the dynamical decoupling of the gas from the underlying DM potential. During the initial stage of gas collapse, the DM potential dominates at all radii and the baryon density is lower than the DM density everywhere. However, since the DM potential of model B is more cuspy than that of model D (Table 1 and Fig. 2), its gas contracts more, and the surface density of the forming disk is substantially higher. The disk continues to grow as a result of ongoing accretion and the settling down of the turbulent gas. We, therefore, turn to the details of this turbulent dynamics.
4.3. The role of turbulence
ENZO has been extensively tested to handle supersonic turbulent motions in a uniform density background (e.g., Kritsuk et al., 2007; Kitsionas et al., 2009; Padoan et al., 2009; Kritsuk et al., 2011b) and in stratified densities (Kritsuk et al., 2011a).
Figure 5 shows the radial and tangential velocity profiles of the flow in Mach numbers, calculated using the velocity and temperature maps for model B. The right-hand maximum corresponds to the initial stage of the gravitational collapse and exceeds . The left-hand maximum reflects the accelerated runaway in the central region and exhibits the same range of Mach numbers. Both the radial and tangential flows are clearly supersonic.
The upper frames of Figure 13 show the velocity field for model B within the central 10 pc at the time of the central runaway. The face-on disk displays spiral shocks in the gas, while the velocity field is dominated by rotation. At the same time, the edge-on disk shows a turbulent velocity field and a number of eddies. The geometrically-thin disk appears substantially puffed-up behind the radial shock. Because the temperature maps reveal that the shocks are nearly isothermal, the concurrent increase in the vertical thickness of the disk can come only from the turbulent flow behind the standing radial and vertical shocks. Indeed, our estimates of thermal pressure gradients in the -direction reveal that thermal pressure gradients are insufficient to puff up the disk to the extent seen in the upper right frame of Figure 13.
While turbulence is notoriously difficult to define, we follow the definition which relies on the vorticity , and its cross product with the velocity field, i.e., the inertial vortex force . To quantify the turbulence within the disk, we have followed the vorticity field within the computational box (lower frames of Fig. 13). As expected, on larger spatial scales the velocity field is irrotational due to the relaxed initial conditions used in our simulations. As the inflow velocity grows with time and the velocity field becomes less regular, the vorticity increases and exhibits a discontinuity at the standing shock which envelops the disk. The postshock flow shows a sharp increase in the vorticity, which decays toward the equatorial plane. In the postshock region, the turbulence is transonic (Table 2). It provides support for the vertical disk structure. Its decay, when moving radially-inward from the shock, results in the sharp decrease of the disk vertical thickness at smaller radii.
The vorticity in the disk appears to be driven by the spiral shocks and by the bar-like perturbation at the center (at later times) — the spiral arms around this perturbation are turbulent as well, as can be seen in the face-on disk region. The central region, which experiences the runaway collapse, exhibits the highest vorticity.
We have also sampled the turbulent velocity field (in terms of the Mach number) on smaller spatial scales, at the time of the central runaway, using spherical sampling volumes whose positions in the disk plane at radii are given in Table 2. The evolution of the turbulent velocity within the central sphere (i.e., first line in Table 2) is also displayed in Figure 14.
An alternative method for measuring the properties of supersonic turbulence is the PDF of the gas density. As discussed in Section 2, the density PDF for homogeneous supersonic turbulence is expected to follow a lognormal distribution (e.g., Vazquez-Semadeni, 1994; Padoan, 1995; Padoan & Nordlund, 2002; Krumholz & McKee, 2005; Federrath et al., 2011). To estimate the PDF in different regions of model B, we measure the grid densities within sampling spheres centered at various points of interest. The sampling of the central region is centered at (0, 0, 0). For each sampling center, we choose the radius of a sphere to enclose 1,500 AMR cells. These spheres sample different dynamical states of the collapsing gas: the central runaway, the outskirts of the disk, the transition region between the disk and the halo, and the gas at large. We carefully choose the size of these spheres to sample a large enough number of grid cells to obtain reasonable statistics, without mixing different dynamical regions in a single sphere.
Figure 15 shows the volume-averaged density PDF in the central region of model B during the central runaway, at Myr. We are unable to fit a single analytic lognormal PDF, for a given mean density (Equation 2), to the entire density range. Instead, we use a combination of a lognormal and a power-law distribution, with the power-law tail dominating at higher densities. At lower densities, comparison between the measured density PDF and the analytical lognormal PDF shows excellent agreement, with the lognormal PDF fit extending over 4 decades in . For , the best fit is a power law tail with the slope of for about 3.5 decades in . No such tail is detected for other locations of sampling spheres.
Such a power-law tail has been observed previously in two-dimensional simulations of homogeneous, supersonic hydrodynamical turbulence in its early stages, before the formation of self-gravitating clumps, i.e., before fragmentation (Scalo et al., 1998). This stage is similar to the present models, albeit with an important difference — our central region is also in a state of a supersonic gravitational collapse and therefore has developed a substantial density gradient. In comparison, the Scalo et al. (1998) power-law tail has a measured slope of before the self-gravitating clumps have formed, and a slope of when these clumps are present. In this respect, the central runaway in our models agrees well with the similar dynamic state in the gas of Scalo et al. — it corresponds to the gravitational collapse of a “fragment.” Such a power law is predicted for the solutions of Burgers equation (which is pressureless) at high densities (e.g., Gotoh & Kraichnan, 1993). Power-law tails have also been observed in three-dimensional simulations (e.g., Federrath et al., 2008; Kritsuk et al., 2011a).
In a more recent work, Kritsuk et al. (2011a) has targeted isothermal supersonic turbulent flows in the presence of gas self-gravity, for the purpose of determining the mass density PDF. A random force has been used to drive the turbulence on large spatial scales, in contrast with our models where turbulence is driven by gravitational collapse only. Despite these differences, we are in agreement that the power-law tail in the PDF appears at the time of the central runaway, when and where the local gravity in the gas becomes important, albeit the slopes are different at the end of the simulations: for Kritsuk et al. (2011a) and for our models A – C. Sampling away from the central runaway site in our models does not show the power-law tail, indicating that there are no other self-gravitating fragments within the computational box.
Why is the power-law slope shallower in our simulations, i.e., vs ? Kritsuk et al. (2011a) comment that a shallower, slope does appear at high densities and associate this with the mass pile-up resulting from dynamically important angular momentum in the region. Recall that the central runaway in our simulations is angular momentum-dominated, as we show in Section 4.1.
We have tested the origin of this power-law PDF by sampling the region with concentric spheres having progressively smaller radii. We find that the range fit by the lognormal PDF shrinks along this sequence, primarily because of cutting off the lower densities, while the high-density end remains untouched (compare the green and blue histograms in Figure 15). The high-density end of the power-law PDF, therefore, corresponds to the very high central density in our simulations, and the density stratification in the central region is responsible for the slope power-law PDF. Clearly, understanding of the power tail in the PDF requires additional theoretical analysis.
In our models discussed here, the initial conditions are simplified and do not include turbulent motions. This delays the onset of turbulence, which takes time to fully develop in the outer part of the flow. In more realistic models, the infalling gas is expected to be already partly turbulent. However, one can also argue in favor of initially laminar and mildly sub-Alfvenic flow in minihalos (e.g., Abel et al., 2002; Bromm & Larson, 2004; Yoshida et al., 2006; Greif et al., 2009; Schleicher et al., 2009). As the infall develops, turbulence sets in spontaneously and greatly increases after the gas has crossed the standing shock enveloping the disk within the central few parsecs. Turbulence which has developed during gravitational collapse has been noticed by other authors (Levine et al., 2008; Wise et al., 2008; Regan & Haehnelt, 2009), and its impact on the gravitational stability of the flow has been analyzed by Begelman & Shlosman (2009). As discussed in Section 2, the effect of turbulence developing during gravitational collapse is to suppress gas fragmentation.
The absence of fragmentation in our simulations of collapsing flow is evident in Figure 1, which provides a snapshot of the face-on disk in model B. A simple estimate based on the floor temperature of the gas with a primordial abundance, e.g., K, within a DM halo shows that the ratio of the gas mass within some spherical radius within this halo to the Jeans mass at this temperature and density is of order . Because the Jeans mass depends on the rms velocity to the third power, this ratio will decrease to in a flow with transonic turbulent velocities.
4.4. Randomizing gas motions: model Dmod
We have rerun model D with less organized baryonic initial momenta (Section 4.2 and Table 1). The purpose of this run is to introduce more realistic initial conditions expected in the cosmological context. In model Dmod, the inner spherical shells have their slightly decreased and their orientation randomized. To keep the total angular momentum unchanged, the outer shells have to compensate for this and their has been slightly increased. The large-scale evolution is somewhat different than that of model D. Namely, the centrifugal barrier has moved inward somewhat. The disk forms with higher surface density and is geometrically thicker. The disk thickness also appears independent of , unlike in Figs. 10 and 13. We also observe that the position of the disk’s equator is less stable relative to the equatorial plane of the DM halo, moving periodically along the -axis, and that the density peak in the disk experiences some weak, low-amplitude perturbations.
Probably the most important difference is the dramatic weakening of the outer torus. This evolution has a dual effect — the timescale for the onset of the central runaway is shortened to 4.5 Myr, and the off-center fragmentation does not materialize because the surface density of the torus is very low. As a result, unlike in model D, model Dmod experiences a ‘classical’ central runaway, similar to models A – C.
5. Discussion and Summary
We have studied some of the physical processes that can operate during the first stages of SMBH formation at high within the direct collapse paradigm. The initial conditions used in this work have been substantially simplified and consist of an isolated, responsive, spherical DM halo of and a virial radius of kpc, having a spin parameter ; embedded baryons with an average cosmological fraction; a universal angular momentum distribution; and nonsingular isothermal density profiles for DM and gas. We limit ourselves to the primordial composition and the absence of molecular cooling.
In a set of high-resolution numerical models we have only varied the size of the DM density core, i.e., the region where the DM density is constant. The collapsing flow is followed down to spatial scales pc (20 AU), over a dynamic range of decades. In these simulations, we have assumed that the inflow is optically thin. This is tested a posteriori — the flow remains optically thin for electron scattering and free-free opacities. The Ly photons produced in the inflow have a large optical depth, but our estimates show that they do not have an effect on the dynamics of the gravitational collapse, in agreement with Rees & Ostriker (1977). Our modeling has been stopped at the approximate radius where the physical conditions of radiation transfer are expected to modify the flow substantially. The inner boundary of the flow should also be investigated in connection with additional physical processes, e.g., magnetic torques. We expect the weak magnetic field advected across the virial radius to be amplified by the dynamo, and the compression to become significant here. MHD processes may also trigger an outflow — all this remains outside the scope of this work.
Two processes can dramatically affect the outcome of direct collapse and prevent a seed SMBH with substantial mass from forming: efficient fragmentation and an angular momentum barrier. The former process can lead to star formation, which will sap the available mass supply and can also expel baryons from the DM halo altogether, e.g., by supernovae feedback and massive stellar winds. The latter process can stop the collapse at the centrifugal barrier, which has been estimated to lie at (e.g., Mo et al., 1998). Our results show that efficient fragmentation has been damped by the development of supersonic turbulence, as suggested by Begelman & Shlosman (2009) — this is especially true during the 2nd stage of the collapse. We also find that the angular momentum is not conserved within the centrifugal barrier — both the outer baryonic collapse and increasing self-gravity in the interior flow trigger the growth of the bar-like mode, which channels the gas inward. With more realistic initial conditions, the typical triaxiality of the DM halo will amplify the non-conservation of angular momentum.
To confirm that fragmentation is indeed damped in the disk forming within the central pc, we have estimated the fragmentation timescale using the analytical approximation provided by Hopkins & Christiansen (2013) for fragmentation in the proto-planetary disks (section 3 there). The resulting characteristic timescale for the disk fragmentation in model B can be expressed as , where is the disk thickness-to-radius ratio, is the Toomre’s parameter, and erfc() is the complementary error function of . The typical values in the model B are and . This results in typical and erfc() . Hence Hubble time. So indeed such disks are not expected to fragment during their lifetime of a few yrs, confirming our numerical results. We note, that the same estimates for the tori bring down to yrs, again as observed in the simulations.
Overall, we find that direct collapse within DM halos depends on the competition between two timescales — that of the central runaway within the centrifugal barrier, and that of off-center fragmentation. If we take a broader view, the centrifugal barrier appears to be a typical rather than exceptional feature of such a collapse. However, in all models it is situated initially much deeper than anticipated, at pc compared to the expected pc for such DM halos. In models with larger DM cores, the central runaway time is delayed, and the radial shock has time to propagate much farther out to pc.
Why does the centrifugal barrier lie so much deeper than anticipated? The reason for this is the inside-out development of the collapse, where only the inner gas has time to reach the barrier. For example, in model B, the centrifugal barrier and the radial shock, which delineate it, form at about 1 pc and advance to about 10 pc by the end of the simulation. We note an important detail — our simulations extend over a timescale which is much shorter than the global free-fall time, Myr, in these DM halos. The central runaway is triggered within the first of this time, and lasts for Myr. This explains why the supersonic turbulence did not decay in our models and why the fragmentation process is so inefficient.
The most intriguing characteristics of the central runaway are that it is self-similar and disky. This means that the angular momentum is dynamically important down to the optically-thick boundary of about 10 AU. The collapsing region is partially supported by rotation in the radial direction and pressure supported in the vertical direction. Published models of SMS (e.g., Begelman, 2010), do assume a degree of rotation in order to ensure stability, but it is well below the dynamically important rotation encountered near the inner boundary of our models (e.g., Montero et al., 2012). It is, therefore, natural to assume that the transition from a rotationally-supported entity to a pressure-supported SMS happens close to this boundary. However, the details of this transition are completely unclear. Moreover, a possibility exists that the runaway collapse dominated by will bypass the state of a hydrostatic equilibrium, that thermonuclear fusion will not play a major role, and that the collapse will proceed directly to forming the SMBH horizon.
To properly follow up the gravitational collapse it is critical to resolve the centrifugal barrier and the associated radial shock, i.e., to have a spatial resolution of a fraction of a parsec. At lower resolution the evolution can diverge from the one we observe.
Usage of a randomized orientation leads to the following corollaries. The centrifugal barrier moves slightly inward, the disk becomes somewhat smaller — these changes do not appear significant. But the disk surface density increases substantially. As a result, the central runaway in model Dmod happens after Myr. The corresponding model D exhibits off-center fragmentation instead. The largest difference between these models is the absence of the massive torus that dominates the outer disk in D. Only a trace of this configuration remains in Dmod, and it is stable against fragmentation. So more random initial conditions move model D into the ‘mainstream’ of models A – C. Clearly, cosmological initial conditions will show the most realistic solution.
Another important requirement is to resolve the supersonic turbulence which develops in various parts of the collapsing flow, and especially behind the radial shock and during the central runaway. The relative absence of turbulence between the shock and the virial radius comes from the quiescent flow in this region — a direct consequence of our initial conditions of an isolated DM halo. In more realistic cosmological initial conditions the laminar flow around may be already turbulent.
We have analyzed the density PDF in the central runaway and found that it is not the lognormal PDF typically encountered in simulations of non-self-gravitating isothermal supersonic flows. The PDFs in models A – C consist of the usual lognormal part as well as a high-density power-law part. The slope of the power-law is found to be at the time and position of the central runaway. Limiting the sampling of the density fluctuations to a smaller region close to the very center, we find that the lognormal PDF fades away but the power law part remains intact (Fig. 15). The lognormal density PDF extends over 4 decades in density and the power-law extends over 3.5 additional decades at higher density. Comparison with two-dimensional hydrodynamical simulations with gas self-gravity (Scalo et al., 1998, see also Kritsuk et al. 2011 for 3-D simulations) confirms that the formation of self-gravitating clumps in the presence of supersonic turbulence depends on the position of the velocity sampling and shows the same structure of a lognormal power-law PDF.
5.1. Estimating the seed SMBH mass range
We now turn to estimating the seed SMBH masses. Table 3 shows the onset time of the central runaway, , in models A – C, increasing from 2.1 Myr (A), to 4.7 Myr (B), to 8.7 Myr (C). Models D and E do not show central collapse before we observe off-center fragmentation. We have calculated the radial dependence of the mass accretion rate, , for all models (e.g., Fig. 6 for model B). In all cases, the central runaway extends radially over a fraction of the central disk. We use the radial mass accretion rate profiles to estimate the baryon mass that participates in the central runaway. In Figure 6 for model B, as well as for models A and C, we observe that these baryons are the ones located within about the half-radius of the disk at the runaway time, as given in Table 3. Note that the disk radius is given by the position of the radial shock at the centrifugal barrier. This also corresponds to the global minimum of at that time — baryons within this radius effectively decouple from the DM background and are dumped onto the center. Baryon masses which are part of this breakaway are also given in Table 3 and range between and .
We now attempt to assess the validity of these estimates. Plotting as function of the onset of the central runaway collapse time, , gives a nearly perfect log-linear dependence, , namely,
where and . On the other hand, depends linearly on the size of the DM density core in models A – C, , given in Table 1. Assuming that has a direct relationship to the mass of the SMBH seed, , models with larger , which lead to larger , should also result in larger . However, an upper limit on appears to come from the condition for off-center fragmentation in the torus — this limit comes from models D and E, which both show off-center fragmentation at Myr. Hence Myr, which is rather a conservative estimate as the fragments will need some time to affect the central runaway dynamics.
The upper limit of Myr intersects the line at pc. Models with larger should exhibit off-center fragmentation. We test this on models D and E — both lie to the right of 1.27 pc (Table 1). So our simplistic argument has passed the first test successfully. What have we learned from this reasoning?
The central runaway drains baryons within the radius , when the gas accumulation inside this radius roughly exceeds that of the DM. The collapse time can be estimated roughly as , where is the local, i.e., within , free-fall timescale. Thus the collapse timescale is yrs, where and . One should consider that baryons inside can be replenished, in principle, as the material flows in across the radial and surface shocks. This would determine a characteristic timescale which may be an order of magnitude above the estimated collapse time.
Probably the most intriguing consequence of this argument is the emerging mass range for the SMBH seeds. If a large fraction of the overall inflow goes into formation of the SMBH seed, we can extrapolate to obtain the maximal , which is about 10% of the amount of baryons in DM halos of interest, . So the mass range for SMBH seeds appears to be . If, in addition, the size of the flat DM density core correlates with the halo virial radius, the mass of the SMBH seed is expected to correlate with the DM halo mass, at the time of formation. This also hints at the possible correlation between the DM halo mass function and the SMBH seed mass function.
Our models relate the properties of SMBHs formed through direct collapse to the sizes of the flat density cores of DM halos. Pure DM simulations (e.g., NFW) have claimed universal density profiles with a cusp, while observations hint rather at the existence of flat density cores (e.g., Flores & Primack, 1994; de Blok, 2005; Primack, 2009, for review). A possible explanation for the flattening of NFW density cusps, appealing to the action of clumpy baryons (El-Zant et al., 2001, 2004), has been verified in numerical simulations (Romano-Diaz et al., 2008). Other solution within the CDM paradigm rely on baryon energy feedback (e.g., Mashchenko et al., 2006).
The size of the DM density cusp in the NFW profiles, and, therefore, the size of the DM density core replacing the cusp, strictly correlate with the halo virial radius. This assumption is probably overly optimistic, and relies heavily on the fragility of the cusp due to its thermodynamic improbability (El-Zant et al., 2001). Nevertheless, we point out that such a correlation will lead to an SMBH mass which initially depends linearly on the DM halo mass. In this case the SMBH seed mass can be inferred from Table 3 to lie at for a DM halo of and kpc, which will have a DM density core of pc. This is close to the upper limit on we have estimated above. It is tantalizing that this upper limit lies so close to the characteristic lowest detected mass of the SMBHs in galaxies at present (e.g., Kormendy & Ho, 2013), and can explain this cutoff. If, however, does not correlate with , the above estimate of an SMBH mass range appears to be more realistic.
The above conclusions might be modified if a substantial amount of the collapsing baryons is expelled via some feedback from the SMS (e.g., Hosokawa et al., 2011; Johnson et al., 2013) or wind mechanism, effectively decreasing the peak accretion below its nominal value of . Moreover, the proposed range of is attributed only to a single cycle in the accretion process, by which we mean one central runaway resulting in the formation of the SMBH seed. The conditions leading to the second cycle will differ because of the existence of the SMBH. However, it is not clear if the mechanical and radiative feedback from this seed will have an effect on the next runaway, i.e., on the 2nd cycle. One can envision that the feedback is directed along the rotation axis, while the next central runaway proceeds in the equatorial plane.
Finally, we note that while our initial calculations of SMBH formation in the direct collapse scenario have emphasized some interesting outcomes, a long list of issues to be resolved remains. One such issue is whether molecular hydrogen can affect the outcome of this process by inducing fragmentation in the collapsing gas. Nearby stars can contribute to the UV background which have an adverse effect on H formation (e.g., Dijkstra et al., 2008). In this work we assume that the UV background will damp H formation and, therefore, will maintain the gas temperature not much below . Several physical mechanisms have been proposed that can support this state of the gas (Omukai, 2001; Oh & Haiman, 2002; Spaans & Silk, 2006; Shang et al., 2010; Schleicher et al., 2010; Latif et al., 2011). Implementation of radiative transfer calculations to study the details of this process is a next logical step. We note that Begelman & Shlosman (2009) have argued that fragmentation will be suppressed even if the cooling floor moves substantially below K. This happens because the flow becomes much more supersonic and the fraction of fragmenting gas decreases with increasing . Alternatively, if the collapse happens inside more massive halos, the ratio of the virial velocity to the sound speed will increase and lead to the same result.
Our results can be compared to some extent with the concurrent works available in the literature, using ENZO (Wise et al., 2008; Latif et al., 2013) and RAMSES (Prieto et al., 2013). All these works used cosmological initial conditions which provide a more realistic setting for the gravitational collapse, but provide less leverage when studying its detail, and are also more time-consuming, limiting the number of models run. The outcome of these models agrees generally with our results of rotationally-dominated disks, and turbulence-suppressed or delayed fragmentation. Prieto et al. (2013), obtains rotationally-supported “cores” in only 3 out of 19 cases. This can be simply explained by the maximal resolution of their models limited to 8 pc (but typically larger, e.g., 14, 15 and 22 pc). At this resolution the disk-like structure, and even the radial shock positioned at the centrifugal barrier, obtained in our simulations would remain unresolved, and the 2nd stage of the collapse will be missed. Based on our model Dmod (Section 4.4), we expect that cosmological initial conditions will lead to the formation of a rotationally-supported disk at somewhat smaller radii than in our models. This would explain why Prieto et al. (2013) missed this runaway stage altogether.
Latif et al. (2013) has imposed a specific subgrid turbulence model of Schmidt et al. (2006). This model has been calibrated against the subsonic turbulence regime. We find that the supersonic turbulence regime operates in various places of the DM halo, especially during the 2nd stage of the collapse. Furthermore, Latif et al. (2013) does not consider the role of gravitational torques in the angular momentum transfer during the collapse, while we find it to be of a prime importance, especially when triggering the 2nd stage of the collapse. Their figure 4 clearly exhibits the dominant barlike or spiral mode, and the gravitational torques are expected to play at least some role in the gas inflow. Despite these differences in the interpretation, our results broadly agree, especially regarding the product of the collapse — the central disklike configuration. The same applies to Wise et al. (2008), who also have concluded that gravitational torques appear as a main mechanism for the angular momentum redistribution in the system.
Although our initial conditions have been motivated by the current cosmology framework, they are significantly simplified and idealized. Owing to this simplification, the simulation results can qualitatively demonstrate the physical processes that work but cannot quantitatively predict the physical timescales discussed here. Varying the DM halo profile, gas density profiles, and angular momentum distribution can affect the bar formation timescale and the torus fragmentation timescale. For example, Koushiappas et al. (2004) and Lodato & Natarajan (2006) suggested that early SMBHs tend to be formed in halos with a low angular momentum. It will be interesting to predict the environments of early SMBHs formed through direct collapse, using full cosmological simulations with many different halo conditions.
- Abel et al. (1997) Abel, T., Anninos, P., Zhang, Y., & Norman, M. L. 1997, New Astron., 2, 181
- Abel et al. (2002) Abel, T., Bryan, G. L., & Norman, M. L. 2002, Science, 295, 93
- Allgood et al. (2006) Allgood, B., Flores, R.A., Primack, J.R., Kravtsov, A.V., Wechsler, R.H., Faltenbacher, A., & Bullock, J.S. 2006, MNRAS, 367, 1781
- Anninos et al. (1997) Anninos, P., Zhang, Y., Abel, T., & Norman, M. L. 1997, New Astron., 2, 209
- Balbus & Hawley (1998) Balbus, S. A., & Hawley, J. F. 1998, Reviews of Modern Physics, 70, 1
- Barnes & Efstathiou (1987) Barnes, J., & Efstathiou, G. 1987, ApJ, 319, 575
- Begelman (2010) Begelman, M. C. 2010, MNRAS, 402, 673
- Begelman et al. (2008) Begelman, M. C., Rossi, E. M., & Armitage, P. J. 2008, MNRAS, 387, 1649
- Begelman & Shlosman (2009) Begelman, M. C., & Shlosman, I. 2009, ApJ, 702, L5
- Begelman et al. (2006) Begelman, M. C., Volonteri, M., & Rees, M. J. 2006, MNRAS, 370, 289
- Blandford & Payne (1982) Blandford, R. D., & Payne, D. G. 1982, MNRAS, 199, 883
- Berentzen et al. (2006) Berentzen, I., Shlosman, I. & Jogee, S. 2006, ApJ, 637, 582
- Berentzen & Shlosman (2006) Berentzen, I., & Shlosman, I. 2006, ApJ, 648, 807
- Berger & Colella (1989) Berger, M. J., & Colella, P. 1989, Journal of Computational Physics, 82, 64
- Bromm et al. (2002) Bromm, V., Coppi, P. S., & Larson, R. B. 2002, ApJ, 564, 23
- Bromm & Loeb (2003) Bromm, V., & Loeb, A. 2003, ApJ, 596, 34
- Bromm & Larson (2004) Bromm, V., & Larson, R.B. 2004, ARA&A, 42, 79
- Bryan & Norman (1998) Bryan, G. L., & Norman, M. L. 1998, ApJ, 495, 80
- Bryan & Norman (1997) Bryan, G. L., & Norman, M. L. 1997, in Astronomical Society of the Pacific Conference Series, Vol. 123, Computational Astrophysics; 12th Kingston Meeting on Theoretical Astrophysics, ed. D. A. Clarke & M. J. West, 363
- Bryan et al. (1995) Bryan, G. L., Norman, M. L., Stone, J. M., Cen, R., & Ostriker, J. P. 1995, Computer Physics Communications, 89, 149
- Bullock et al. (2001) Bullock, J. S., Dekel, A., Kolatt, T. S., Kravtsov, A. V., Klypin, A. A., Porciani, C., & Primack, J. R. 2001, ApJ, 555, 240
- Christodoulou et al. (1995) Christodoulou, D. M., Shlosman, I., & Tohline, J. E. 1995, ApJ, 443, 551
- Clark et al. (2011) Clark, P. C., Glover, S. C. O., Smith, R. J., et al. 2011, Science, 331, 1040
- de Blok (2005) de Blok, W.J.G. 2005, ApJ, 634, 227
- Dijkstra et al. (2008) Dijkstra, M., Haiman, Z., Mesinger, A., & Wyithe, J.S.E. 2008, MNRAS, 391, 1961
- El-Zant et al. (2001) El-Zant, A., Shlosman, I., & Hoffman, Y. 2001, ApJ, 560, 636
- El-Zant et al. (2004) El-Zant, A., Hoffman, Y., Primack, J., Combes, F., & Shlosman, I. 2004, ApJ, 607, L75
- Englmaier & Shlosman (2004) Englmaier, P., & Shlosman, I. 2004, ApJ, 617, L115
- Fan et al. (2003) Fan, X., et al. 2003, AJ, 125, 1649
- Federrath et al. (2008) Federrath, C., Glover, S.C.O., Klessen, R.S., & Schmidt, W. 2008, Phys. Scr. T, 132, 014025
- Federrath et al. (2011) Federrath, C., Sur, S., Schleicher, D. R. G., Banerjee, R., & Klessen, R. S. 2011, ApJ, 731, 62
- Ferrarese & Merritt (2000) Ferrarese, L., & Merritt, D. 2000, ApJ, 539, L9
- Ferrarese (2005) Ferrarese, L. 2005, Space Sci. Rev., 116, 523
- Flores & Primack (1994) Flores, R.A., & Primack, J.R. 1994, ApJ, 427, L1
- Garlick (1979) Garlick, A.R. 1979, A&A, 73, 171
- Gebhardt et al. (2000) Gebhardt, K., et al. 2000, ApJ, 539, L13
- Goldreich et al. (1996) Goldreich, P., Lai, D., & Sahrling, M. 1996, Unsolved Problems in Astrophysics, eds. J.N.Bahcall & J.P. Ostriker (Princeton: Princeton Univ. Press)
- Gotoh & Kraichnan (1993) Gotoh, T., & Kraichnan, R.H. 1993, Phys.Fluids A, 5, 445
- Greif et al. (2009) Greif, T. H., Johnson, J. L., Klessen, R. S., & Bromm, V. 2009, MNRAS, 399, 639
- Haiman & Loeb (2001) Haiman, Z., & Loeb, A. 2001, ApJ, 552, 459
- Harrington (1973) Harrington, P.J., 1973, MNRAS, 162, 43
- Heller et al. (2007) Heller, C.H., Shlosman, I., & Athanassoula, E. 2007, ApJ, 671, 226
- Hopkins & Christiansen (2013) Hopkins, P. F., & Christiansen, J. L. 2013, arXiv:1301.2600
- Hosokawa et al. (2011) Hosokawa, T., Omukai, K., Yoshida, N., & Yorke, H. W. 2011, Science, 334, 1250
- Hoyle (1953) Hoyle, F. 1953, ApJ, 118, 513
- Hunter (1962) Hunter, C. 1962, ApJ, 136, 594
- Jahnke & Maccio (2011) Jahnke, K., & Maccio, A.V. 2011, ApJ, 734, 92
- Johnson et al. (2011) Johnson, J. L., Khochfar, S., Greif, T. H., & Durier, F. 2011, MNRAS, 410, 919
- Johnson et al. (2013) Johnson, J. L., Whalen, D. J., Even, W., et al. 2013, arXiv:1304.4601
- Kitsionas et al. (2009) Kitsionas, S., Federrath, C., Klessen, R. S., et al. 2009, A&A, 508, 541
- Klessen & Hennebelle (2010) Klessen, R. S., & Hennebelle, P. 2010, A&A, 520, A17
- Komatsu et al. (2009) Komatsu, E., et al. 2009, ApJS, 180, 330
- Kormendy & Ho (2013) Kormendy, J., & Ho, L. C. 2013, ARAA, in press, arXiv:1304.7762
- Koushiappas et al. (2004) Koushiappas, S. M., Bullock, J. S., & Dekel, A. 2004, MNRAS, 354, 292
- Kritsuk et al. (2007) Kritsuk, A.G., Norman, M.L., Padoan, P., & Wagner, R. 2007, ApJ, 665, 416
- Kritsuk et al. (2011a) Kritsuk, A.G., Norman, M.L., & Wagner, R. 2011, ApJ, 727, L20
- Kritsuk et al. (2011b) Kritsuk, A. G., Nordlund, Å., Collins, D., et al. 2011, ApJ, 737, 13
- Krumholz & McKee (2005) Krumholz, M. R., & McKee, C. F. 2005, ApJ, 630, 250
- Lai & Goldreich (2000) Lai, D., & Goldreich, P. 2000, ApJ, 535, 402
- Latif et al. (2011) Latif, M. A., Zaroubi, S., & Spaans, M. 2011, MNRAS, 411, 1659
- Latif et al. (2013) Latif, M. A., Schleicher, D. R. G., Schmidt, W., & Niemeyer, J. 2013, MNRAS, 1155
- Levine et al. (2008) Levine, R., Gnedin, N. Y., Hamilton, A. J. S., & Kravtsov, A. V. 2008, ApJ, 678, 154
- Li et al. (2007) Li, Y., et al. 2007, ApJ, 665, 187
- Lichtenberg & Lieberman (1995) Lichtenberg, A.J., & Lieberman, M.A 1995, Regular & Chaotic Dynamics, Springer
- Lodato & Natarajan (2006) Lodato, G., & Natarajan, P. 2006, MNRAS, 371, 1813
- Lynden-Bell & Kalnajs (1972) Lynden-Bell, D., & Kalnajs, A. J. 1972, MNRAS, 157, 1
- Mashchenko et al. (2006) Mashchenko, S., Couchman, H.M.P., & Wadsley, J. 2006, Nature, 442, 539
- Mayer et al. (2010) Mayer, L., Kazantzidis, S., Escala, A., & Callegari, S. 2010, Nature, 466, 1082
- Mo et al. (1998) Mo, H,J., Mao, S., & White, S.D.M. 1998, MNRAS, 295, 319
- Moncrief (1980) Moncrief, V. 1980, ApJ, 235, 1038
- Montero et al. (2012) Montero, P. J., Janka, H.-T., & Müller, E. 2012, ApJ, 749, 37
- Mortlock et al. (2011) Mortlock, D. J., et al. 2011, Nature, 474, 616
- Navarro et al. (1996) Navarro, J. F., Frenk, C. S., & White, S. D. M. 1996, ApJ, 462, 563
- Navarro et al. (1997) —. 1997, ApJ, 490, 493
- Neufeld (1990) Neufeld, D.A., 1990, ApJ, 350, 216
- Norman & Bryan (1999) Norman, M. L., & Bryan, G. L. 1999, in Astrophysics and Space Science Library, Vol. 240, Numerical Astrophysics, ed. S. M. Miyama, K. Tomisaka, & T. Hanawa, 19
- Oh & Haiman (2002) Oh, S. P., & Haiman, Z. 2002, ApJ, 569, 558
- Omukai (2001) Omukai, K. 2001, ApJ, 546, 635
- O’Shea & Norman (2007) O’Shea, B. W., & Norman, M. L. 2007, ApJ, 654, 66
- Ostriker et al. (1999) Ostriker, E. C., Gammie, C. F., & Stone, J. M. 1999, ApJ, 513, 259
- Padoan (1995) Padoan, P. 1995, MNRAS, 277, 377
- Padoan & Nordlund (2002) Padoan, P., & Nordlund, A. 2002, ApJ, 576, 870
- Padoan et al. (2009) Padoan, P., Juvela, M., Kritsuk, A.G., & Norman, M.L. 2009, ApJ, 707, 153
- Padoan et al. (2007) Padoan, P., Nordlund, A., Kritsuk, A.G., Norman, M.L., & Li, P.S. 2007, ApJ, 661, 972
- Peebles (1969) Peebles, P. J. E. 1969, ApJ, 155, 393
- Pelupessy et al. (2007) Pelupessy, F. I., Di Matteo, T., & Ciardi, B. 2007, ApJ, 665, 107
- Prieto et al. (2013) Prieto, J., Jimenez, R., & Haiman, Z. 2013, MNRAS, submitted, arXiv:1301.5567
- Primack (2009) Primack, J.R. 2009, New J. of Phys., 11, 105029
- Rees & Ostriker (1977) Rees, M.J., & Ostriker, J.P. 1977, MNRAS, 179, 541
- Regan & Haehnelt (2009) Regan, J. A., & Haehnelt, M. G. 2009, MNRAS, 396, 343
- Romano-Diaz et al. (2008) Romano-Diaz, E., Shlosman, I., Hoffman, Y., & Heller, C. 2008, ApJ, 685, L105
- Romano-Diaz et al. (2009) Romano-Diaz, E., Shlosman, I., Heller, C., & Hoffman, Y. 2009, ApJ, 702, 1250
- Safranek-Shrader et al. (2012) Safranek-Shrader, C., Agarwal, M., Federrath, C., Dubey, A., Milosavljevic, M., & Bromm, V. 2012, MNRAS, 426, 1159
- Salpeter (1964) Salpeter, E. E. 1964, ApJ, 140, 796
- Samland & Gerhard (2003) Samland, M., & Gerhard, O.E. 2003, A&A, 399, 961
- Scalo et al. (1998) Scalo, J., Vazquez-Semadeni, E., Chappell, D., & Passot, T. 1998, ApJ, 504, 835
- Schleicher et al. (2009) Schleicher, D. R. G., Galli, D., Glover, S. C. O., et al. 2009, ApJ, 703, 1096
- Schleicher et al. (2010) Schleicher, D. R. G., Spaans, M., & Glover, S. C. O. 2010, ApJ, 712, L69
- Schmidt et al. (2006) Schmidt, W., Niemeyer, J. C., & Hillebrandt, W. 2006, A&A, 450, 265
- Shang et al. (2010) Shang, C., Bryan, G. L., & Haiman, Z. 2010, MNRAS, 402, 1249
- Shlosman & Begelman (1989) Shlosman, I., & Begelman, M. C. 1989, ApJ, 341, 685
- Shlosman et al. (1989) Shlosman, I., Frank, J., & Begelman, M. C. 1989, Nature, 338, 45
- Shlosman et al. (1990) Shlosman, I., Begelman, M. C., & Frank, J. 1990, Nature, 345, 679
- Shlosman (2005) Shlosman, I. 2005, The Evolution of Starbursts, (eds.) S. Hüttermeister et al., AIP 783, 223
- Shlosman (2013) Shlosman, I. 2013, Secular Evolution of Galaxies, (eds.) J.Falcon-Barroso & J.H.Knapen, CUP, in press, arXiv:1212.1463
- Shu (1977) Shu, F. H. 1977, ApJ, 214, 488
- Spaans & Silk (2006) Spaans, M., & Silk, J. 2006, ApJ, 652, 902
- Stacy et al. (2012) Stacy, A., Greif, T. H., & Bromm, V. 2012, MNRAS, 422, 290
- Subramanian (2010) Subramanian, K. 2010, Astronomische Nachrichten, 331, 110
- Sur et al. (2010) Sur, S., Schleicher, D. R. G., Banerjee, R., Federrath, C., & Klessen, R. S. 2010, ApJ, 721, L134
- Tanaka & Haiman (2009) Tanaka, T., & Haiman, Z. 2009, ApJ, 696, 1798
- Tonini et al. (2006) Tonini, C., Lapi, A., & Salucci, P. 2006, ApJ, 649, 591
- Tremaine et al. (2002) Tremaine, S., et al. 2002, ApJ, 574, 740
- Truelove et al. (1997) Truelove, J.K., Klein, R.I., McKee, C.F., Holliman, J.H., Howell, L.H., & Greenough, J.A. 1997, ApJ, 489, L179
- Turk et al. (2009) Turk, M. J., Abel, T., & O’Shea, B. 2009, Science, 325, 601
- Turk et al. (2011) Turk, M. J., Smith, B. D., Oishi, J. S., et al. 2011, ApJS, 192, 9
- Turk et al. (2012) Turk, M. J., Oishi, J. S., Abel, T., & Bryan, G. L. 2012, ApJ, 745, 154
- Vazquez-Semadeni (1994) Vazquez-Semadeni, E. 1994, ApJ, 423, 681
- Vazquez-Semadeni (2010) Vazquez-Semadeni, E. 2010, in IAU Symp. 270, Computational Star Formation, ed. J. Alves et al. (Cambridge: Cambridge Univ. Press
- Volonteri & Rees (2005) Volonteri, M., & Rees, M. J. 2005, ApJ, 633, 624
- Volonteri & Rees (2006) Volonteri, M., & Rees, M. J. 2006, ApJ, 650, 669
- Wise et al. (2008) Wise, J. H., Turk, M. J., & Abel, T. 2008, ApJ, 682, 745
- Wise et al. (2012) Wise, J.H., Turk, M.J., Norman, M.L., & Abel, T. 2012, ApJ, 745, 50
- Yoo & Miralda-Escudé (2004) Yoo, J., & Miralda-Escudé, J. 2004, ApJ, 614, L25
- Yoshida et al. (2006) Yoshida, N., Omukai, K., Hernquist, L., & Abel, T. 2006, ApJ, 652, 6