SUPERMASSIVE BLACK HOLE FORMATION AT HIGH REDSHIFTS VIA DIRECT COLLAPSE: PHYSICAL PROCESSES IN THE EARLY STAGE

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 $\lambda = J/\sqrt{2}M_{\rm vir}v_{\rm c}R_{\rm vir}$, where J is the halo angular momentum, and M_vir and v_c are its virial mass and the circular velocity at the virial radius R_vir (e.g., Peebles 1969; Barnes & Efstathiou 1987; Bullock et al. 2001). The mean spin of DM halos is λ ∼ 0.05. 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 ∼100. When the DM potential dominates, the collapse will stop after a decade in radius (e.g., Mo et al. 1998; Shlosman 2013, and references 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, m = 1–3, prevail in both the gas and the DM, because their rise times are shortest. The m = 1 mode requires perturbation of the center of mass of the gaseous disk—this is not always possible. The next fastest growing mode is m = 2—the bar-like mode. Typically it will have the highest amplitude, in the absence of m = 1. 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 1992). 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, T/|W|, 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 α ≡ 0.5fT/|W| ⩾ 0.34 (Christodoulou et al. 1995), where f = 1 for disks and f = (2/3) 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. 1997; 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 sphere⁴ρ∝R⁻². 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

$$\begin{eqnarray} &&\dot{M}(R) = \frac{M_{\rm gas}}{t_{\rm ff}}\sim \frac{M_{\rm gas}}{M_{\rm tot}} \frac{v_{\rm ff}^3}{G}\sim 0.4\,v_{30}^3\,M_\odot \,{\rm yr^{-1}}, \end{eqnarray} \tag{ 1 }$$

where v₃₀ ≡ v_ff/30 km s⁻¹ for a 2 × 10⁸ M_☉ 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 the order of ${\sim} \alpha _{\rm v} c_{\rm s}^{3}/G \propto T^{3/2}/G$, where c_s is the gas sound speed and α_v 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, α_v < 1, but for nonlocal viscosities, e.g., magnetic or gravitational torques, α_v ≳ 1, especially in the self-gravitating case where the effective α_v ≫ 1 (Shlosman et al. 1990). In this latter case, the inflow rate is still given by Equation (1), where v_ff is the virial velocity of the gravitational potential which is now dominated by the gas, i.e., M_gas/M_tot = 1 in Equation (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 ∼1 M_☉ yr⁻¹ is necessary in order to form a 10⁶ M_☉ SMS, which has a nuclear burning timescale of only a few million years. The initial SMBH seed of 100–1000 M_☉, 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 ∼10⁵ M_☉ or more in less than a few megayears.

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):

$$\begin{eqnarray} &&p(x) = \frac{1}{\left(2 \pi \sigma ^{2}_{\rm p} \right)^{0.5}} \frac{1}{x} \exp \left[ -\frac{(\ln x - \overline{\ln x})^2}{2\sigma ^{2}_{\rm p}} \right], \end{eqnarray} \tag{ 2 }$$

where the distribution mean $\overline{\ln x} = -0.5 \sigma ^{2}_{\rm p}$, x ≡ ρ/ρ₀ (ρ₀ being mean density), and its dispersion is $\sigma ^{2}_{\rm p} \sim [\ln (1 + 3\mathcal {M}^{2}/4)]$ ($\mathcal {M}$ is the flow Mach number). This distribution is the result of x being a random variable, which is itself a product of independent random variables. Note that the coefficient b = 3/4 in the definition of $\sigma ^2_{\rm p}$ is not a constant and depends on the driving mode of the turbulence. Federrath et al. (2008) found that b 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 x_crit ≡ (λ_J0/λ_s)², where λ_J0 is the Jeans length at the mean density and λ_s 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 x_crit are unstable to fragmentation. Estimating the mass fraction in Jeans-unstable fragments that collapse on a timescale shorter than the free-fall timescale as $f = \int ^{\infty }_{x_{\rm crit}} x ({dp(x)}/{dx}) dx$, we find that f ⩽ 0.02 for $\mathcal {M}\ge 3$. 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 ∼7 × 10⁻³X pc, where X is the ionization fraction. For temperatures encountered in this simulation, X ≪ 1 and the electron scattering optical depth is negligible, as is the free–free absorption. For primordial chemical composition and the range of temperatures ∼3000–10, 000 K, the opacity is expected to be dominated by bound–bound transitions in hydrogen, and especially by the opacity for Lyα 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 >8000 K (e.g., Omukai 2001; Shang et al. 2010). The prime agent of H₂ dissociation can be an externally produced Lyman–Werner (∼11 eV–13 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 T ∼ 10⁴ K 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 ∼8000 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 τ_α ≳ 10¹⁴ (e.g., Harrington 1973; Neufeld 1990), which can be attained at radii ≲ 10⁻³ 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 10⁻⁸ 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 $\mathcal {M}$.

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 B-fields are expected to be weak at high z, and because of the high temperature of the flow, ≳ 1000 K. At these high temperatures, the value of the critical B-field strength needed to limit the compression in supersonic isothermal shocks is also higher (Padoan et al. 2007).

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, 1998; Norman & Bryan 1999). ENZO uses a particle-mesh N-body method to calculate the gravitational dynamics including collisionless DM particles, and a second-order piecewise parabolic method (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 two in length scale if the gas density is greater than ρ₀N^l, where ρ₀ is the minimum density above which refinement occurs, N = 2 is the refinement factor and l = 25 is the maximal AMR refinement level. This refinement corresponds to a spatial resolution of ∼10 AU.

The Truelove et al. (1997) requirement for resolution of the Jeans length, i.e., at least four 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, 10–10³ pc, the Jeans length is resolved with ≳ 100 cells. Between 1–10 pc, it is resolved by ∼32 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 difference in the evolution has been found, except a slight, ∼10%, increase in the onset time of the second 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: H, H⁺, He, He⁺, He⁺⁺, and e⁻ (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 H₂ formation. In order to include these effects without implementing a full radiative transfer calculation, we exclude the chemistry and cooling related to H₂ in this paper.

For a gas with primordial composition and no H₂, the earliest collapse can occur in DM halos whose virial temperatures exceed ∼10⁴ K. For this reason we set up a spherical DM halo with M_vir ∼ 2 × 10⁸ h⁻¹ M_☉ at z = 15. According to the top-hat model, such a halo will have a virial radius of R_vir ∼ 945 h⁻¹ pc and a virial temperature T_vir ∼ 32, 000 K, according to the WMAP5 cosmology (Ω_m = 0.279, Ω_b = 0.0445, and h = 0.701; Komatsu et al. 2009). We simulate this DM halo within a (6 kpc)³ computational domain.

For the DM halo, we assume an isothermal sphere. This halo model is similar to the universal DM halo model (Navarro–Frenk–White (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 10⁶ particles. We introduce a DM core radius, R_dm, 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, R_dm plays the role of the King radius for a nonsingular isothermal sphere. We vary R_dm 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 an R_core = 100 pc 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, Ω_b/Ω_m ∼ 0.16. The angular momentum of the halo gas is computed from $J = \lambda \sqrt{2} M_{\rm vir} v_{\rm c} R_{\rm vir}$, defined in Section 2. In large-scale N-body cosmological simulations, Bullock et al. (2001) has found that the specific angular momentum in DM halos roughly follows j(R) ∼ R^1.1, 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 z-axis. For R ≲ R_core, we assume solid body rotation for the gas, i.e., j(R) ∼ R² (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 ∼R⁻³ beyond R_vir.

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 j of the inner shells and compensated the outer shells by a slight increase in j. 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.

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 r and larger gravitational accelerations there. This leads to the establishment of ρ ∼ R⁻² 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 R_dm, and, therefore, in the depth of the DM potential well. Figure 1 shows face-on density-weighted projections of the density⁵ of the gas disk at t = 4.6 Myr in models B and D, with R_dm ∼ 0.4 pc and 1.5 pc, respectively. Unless mentioned, other models behave similarly.

Zoom In Zoom Out Reset image size

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 r ∼ 1 pc. The disk grows quickly in radius and mass, and by t ∼ 4 Myr, the snapshots display well-developed gas disks with radii ∼10 pc. The disk boundaries, both in r and in z, are delineated by standing shocks, discussed later. The surface densities of all models are well-approximated by a power-law Σ ∼ r⁻ⁿ with n ∼ 1.3. Figure 2 displays the evolution of the surface density for model B, as well as comparisons of the surface densities at t = 4.7 Myr among models B, D, and E. The disks are truncated at r ∼ 10 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 ∼1–2 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.

Zoom In Zoom Out Reset image size

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 ∼10 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 R_dm of less cuspy halos.

Zoom In Zoom Out Reset image size

Over the next time period of a few × 10⁵ 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, j increases, moving the barrier outward. Inside the centrifugal barrier, the disk develops a global instability, dominated by the Fourier component m = 2—a bar-like mode, which transfers j to the outer gas and to the DM. For example, model B exhibits a well defined gaseous bar in the central ∼1–2 pc of the disk at ∼4.7 Myr. The appearance of this bar-like mode is characteristic also of smaller spatial scales of this model (Figure 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 m = 2 axial symmetry breaking. At a later stage, these models exhibit off-center fragmentation.

Zoom In Zoom Out Reset image size

Figure 1(a) 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, m = 1, 2, and occasionally m = 3. Only m = 1 exhibits a mild displacement of the center of mass of the disk, while m = 2 has the appearance of a strong gaseous bar which drives spiral structure, and m = 3 shows up as tri-armed spirals. We follow the runaway collapse to a spatial scale of ∼10⁻⁴ 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 ∼20 AU at the end of the simulation, while model D gas only reaches a resolution of ∼0.1 pc. Resolving the ∼20 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 ∼1 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 t = 4.0 Myr to t = 4.7 Myr, a truly runaway collapse with a dynamical time of ≲ 10⁶ yr.

For the idealized case of spherical accretion within an external gravitational potential given by ρ∝R⁻², the gas mass flux is given by $\dot{M}\sim R^2 \rho v_{\rm R}$ with v_R weakly dependent on R. Figure 3 confirms that the density profile of the collapsing gas in model B tends to the isothermal density profile, ρ∝R⁻², for R ≲ 10 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 5(a) shows the radial profile of the inflow velocity measured in the disk plane. We can divide the r-range into the collapsing gas in the region outside 10 pc, the disk plateau at r ∼ 1–10 pc, and the inner region of runaway collapse—all at the time of the central runaway, t ∼ 4.7 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 v_R toward the center, as the result of the runaway collapse of a finite gas mass.

Zoom In Zoom Out Reset image size

Thus the outermost region is dominated by free-fall kinematics because the angular momentum is low. The density profile is maintained as ρ ∼ R⁻² if one neglects the weak variation of v_R. The intermediate disky region has v_R ∼ constant and the same −2 slope of log ρ. The innermost collapse exhibits a rapidly increasing and steepening infall velocity (Figures 5(a) and (b)). The accretion flow is mostly rotationally supported in the disky region of 1–10 pc, as well as in the central region, where the self-gravitating collapse proceeds with a non-negligible angular momentum. The time evolution of v_R 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, $\dot{M}(R)$, appears to be approximately flat, ∼0.1–0.2 M_☉ yr⁻¹, exterior to the radial shock, in agreement with Equation (1) (Figure 6). In the disk region, $\dot{M}(R)$ decreases abruptly to ∼0.01 M_☉ yr⁻¹ at t ∼ 3 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, t = 4.7 Myr, the central runaway is accompanied by peak values of $\dot{M}(R)\sim 2\,M_\odot \,{\rm yr^{-1}}$. The decrease in the accretion rate in the range of r ∼ 1–5 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 $\dot{M}(R)$ 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.

Zoom In Zoom Out Reset image size

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 Equation (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 M_gas/M_tot ∼ 1 in Equation (1), and the free-fall velocity, v_ff, becomes a rapidly growing function of time. As a result, the characteristic r-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 (z-axis) shocks, as is seen in Figures 5(a), 7, and 13. But there is no associated shock in the azimuthal motion. The face-on and edge-on velocity fields on scales of ∼10 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.

Zoom In Zoom Out Reset image size

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 t ∼ 4.7 Myr, along two directions in the equatorial plane, ±x and ±y, and along the ±z-axis (Figure 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 t ∼ 4.7 Myr, the radial shock in the equatorial plane is positioned at r ∼ 10 pc (see also Figure 5(a)).

Inside the radial shock marking the centrifugal barrier, i.e., between r ∼ 1 pc and 10 pc, the equatorial density increases substantially while the density along the z-axis is almost constant. At ∼1 pc, the density along the z-axis is ∼10⁻³ times the density measured along the x- or y-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 z) in the disk. While the centrifugal barrier stops the cold inflow toward the rotation axis at r ∼ 10 pc, the position of the disk surface shock, i.e., the disk thickness, is determined by the postshock temperature which never exceeds ∼10⁴ K, as seen in Figure 9.

A strong surface shock is maintained by the infall along the z-axis at ∼0.5 pc (Figure 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 x = y = 0. The "peanut" shape of the surface shock at this time can be inferred from Figure 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 r are associated with different dynamical timescales. Hence the fact that ρ(z)/ρ(r) ∼ constant for different r 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 (Figure 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 ∼1–2 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 0.01 pc. The resulting values of the pattern speeds, Ω_bar, confirm that the m = 2 mode tumbles with substantially different speeds. On ∼1 pc scale, Ω_bar ∼ 4 × 10⁻⁵ yr⁻¹, while on ∼0.01 pc scale, Ω_bar ∼ 5 × 10⁻³ yr⁻¹, about two 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, Ω_bar(r) ∼ r^−α, with α ∼ 1. The parameter α is determined by the mass distribution M(r) inside the collapsing gas at the onset of the central runaway. The gas volume density scales as ρ ∼ r⁻² in the central region (Figure 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 Ω_bar(r) ∼ [M(r)/r³]^1/2, where M(r) ∼ r, which explains the instantaneous value of α above.

To verify that the bar-like (m = 2) mode dominates over other modes, e.g., m = 1 and 3, we have Fourier analyzed the gas response to the asymmetric potential. The dominant mode at early times is the m = 1 mode—the disk center of mass is initially perturbed by the asymmetry developed in the overall mass distribution (Figure 8). This mode, however, decays quickly to a negligible amplitude. An additional odd mode, m = 3, is also present, although at lower amplitude—it decays similarly. On the other hand, the even mode m = 2 grows to a substantial amplitude with time. At about t ∼ 4 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 α < 0.34 (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.

Zoom In Zoom Out Reset image size

The temperature variations with density for model B are shown in Figure 9. At the start of disk formation, a radial shock forms at r ∼ 2 pc and then gradually propagates out to 10 pc. At this time the central runaway is triggered. As shown by Figure 9, the posthock gas rapidly cools down to low T in the postshock region because of the high density. Within the growing disk, r ≲ 10 pc, we observe an increasing range in T ∼ 1.5 × 10³–10⁴ K and ρ ∼ 10⁻¹⁷–10⁻²¹ g cm⁻³. The upper limit to the temperature is dropping with density. The central runaway at ≲ 1 pc narrows the T-range to mostly T ∼ 3 × 10³ K but increases the range in ρ ∼ 10⁻¹⁵–10⁻¹⁰ g cm⁻³. The maximal postshock T in the disk is slightly increasing with time, which reflects the increasing shock-impact velocity of the gas which comes from progressively larger distances.

Zoom In Zoom Out Reset image size

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 t = 4.7 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 ∼2–3 pc at t = 4.7 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 ∼8 pc at this time. When the gas distributions in models B and D are compared before the central runaway in model B, at t ∼ 4.0 yr, the core density in B is higher by about two orders of magnitude than in D (Figure 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 m = 1 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 m = 2 mode, which is responsible for the central runaway collapse, can only grow after m = 1 has decayed significantly (Figure 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 m = 1, 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 j, and steadily growing free-fall velocity, arriving from larger r. After ∼6–7 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 (Figure 10). Figure 11 displays the formation and evolution of this torus at r ∼ 20–100 pc. The density profile at t = 13.2 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., Figure 10(a)). Finally, at t ∼ 13.4 Myr, one of the fragments becomes substantially compact and begins gravitational collapse and the simulation is stopped (Figure 12). Model E is similar to D, with off-center fragmentation in the torus occurring at the same time, t ∼ 13.4 Myr.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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 10⁴ M_☉, 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 m = 1 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 Figure 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.

ENZO has been extensively tested to handle supersonic turbulent motions in a uniform density background (e.g., Kritsuk et al. 2007, 2011a; Kitsionas et al. 2009; Padoan et al. 2009) and in stratified densities (Kritsuk et al. 2011b).

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 $\mathcal {M}\sim 3\hbox{--}5$. 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 z-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.

Zoom In Zoom Out Reset image size

While turbulence is notoriously difficult to define, we follow the definition which relies on the vorticity ${\rm \bf w} = {\rm \bf \nabla \times v}$, and its cross product with the velocity field, i.e., the inertial vortex force ${\rm \bf v \times \nabla \times v}$. To quantify the turbulence within the disk, we have followed the vorticity field within the computational box (lower frames of Figure 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.

Table 2. Sampling the Turbulent Velocities in Model B

Center at rc	Sphere Radius	$\mathcal {M}$
(pc)	(pc)
7.2 × 10−3	7.0 × 10−3	1.85
7.5 × 10−3	7.2 × 10−3	1.62
1.0 × 10−1	5.0 × 10−2	0.72
1.0	0.3	0.48
2.5	1.0	0.63
5.0	2.0	0.94

Notes. r_c is the distance from the rotation axis of the center of the sampling sphere (see text for additional explanations). Last column: the estimated Mach number within the sampling sphere.

Download table as:  ASCIITypeset image

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 r_c 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.

Zoom In Zoom Out Reset image size

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 ∼1500 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 t = 4.7 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 log ρ > −12, the best fit is a power-law tail with the slope of ∼ − 1 for about 3.5 decades in log ρ. No such tail is detected for other locations of sampling spheres.

Zoom In Zoom Out Reset image size

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 ∼ − 2 before the self-gravitating clumps have formed, and a slope of −1 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. 2011b).

In a more recent work, Kritsuk et al. (2011b) 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: −1.7 for Kritsuk et al. (2011b) and −1 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., −1.7 versus −1? Kritsuk et al. (2011b) comment that a shallower, −1 slope does appear at high densities and then 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 −1 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-Alfvénic 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., T ∼ 10⁴ K, within a 2 × 10⁸ M_☉ DM halo shows that the ratio of the gas mass within some spherical radius R within this halo to the Jeans mass at this temperature and density is of the order of ∼3–4. Because the Jeans mass depends on the rms velocity to the third power, this ratio will decrease to ∼1 in a flow with transonic turbulent velocities.

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 j slightly decreased and their orientation randomized. To keep the total angular momentum unchanged, the outer shells have to compensate for this and their j 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 r, unlike in Figures 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 z-axis, and that the density peak in the disk experiences some weak, low-amplitude m = 1 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.

We now turn to estimating the seed SMBH masses. Table 3 shows the onset time of the central runaway, t_coll, 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, $\dot{M}(R)$, for all models (e.g., Figure 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 $\dot{M}(R)$ 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 M_coll ∼ 2 × 10⁴ M_☉ and 6 × 10⁵ M_☉.

We now attempt to assess the validity of these estimates. Plotting M_coll as function of the onset of the central runaway collapse time, t_coll, gives a nearly perfect log-linear dependence, log M_coll ∼ t_coll, namely,

$$\begin{eqnarray} &&\log \,(M_{\rm coll}/M_\odot) = a (t_{\rm coll}/1\,{\rm Myr}) + b, \end{eqnarray} \tag{ 3 }$$

where a ∼ 0.18 and b ∼ 3.95. On the other hand, t_coll depends linearly on the size of the DM density core in models A–C, R_dm, given in Table 1. Assuming that M_coll has a direct relationship to the mass of the SMBH seed, M_•, models with larger t_coll, which lead to larger M_coll, should also result in larger M_•. However, an upper limit on t_coll 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 ∼13.4 Myr. Hence t_coll ≲ 13.4 Myr, which is rather a conservative estimate as the fragments will need some time to affect the central runaway dynamics.

The upper limit of 13.4 Myr intersects the t_coll(R_s) line at R_s ∼ 1.27 pc. Models with larger R_dm should exhibit off-center fragmentation. We test this on models D and E—both lie to the right of 1.27 pc (Table 1). Therefore, our simplistic argument has passed the first test successfully. What have we learned from this reasoning?

The central runaway drains baryons within the radius ∼R_coll, when the gas accumulation inside this radius roughly exceeds that of the DM. The collapse time can be estimated roughly as ∼2–3 × t_ff, where t_ff is the local, i.e., within ∼R_coll, free-fall timescale. Thus the collapse timescale is ${\sim } 3\times 10^6\,R_{\rm coll,10}^{3/2}M_{\rm coll,6}^{-1/2}$ yr, where R_{coll, 10} ≡ R_coll/10 pc and M_{coll, 6} ≡ M_coll/10⁶ M_☉. One should consider that baryons inside R_coll 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 log M_coll ∼ t_coll to obtain the maximal M_coll ∼ 2 × 10⁶ M_☉, which is about 10% of the amount of baryons in DM halos of interest, M_vir ∼ 1–2 × 10⁸ M_☉. So the mass range for SMBH seeds appears to be 2 × 10⁴ M_☉ ≲ M_• ≲ 2 × 10⁶ M_☉. 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 solutions 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 correlates 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 M_• ∼ 10⁶ M_☉ for a DM halo of M_vir ∼ 2 × 10⁸ M_☉ and R_vir ∼ 1 kpc, which will have a DM density core of R_dm ∼ 1 pc. This is close to the upper limit on M_• 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, R_dm does not correlate with R_vir, the above estimate of an SMBH mass range 2 × 10⁴ M_☉ ≲ M_• ≲ 2 × 10⁶ M_☉ 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 $\dot{M}\sim 1\hbox{--}2\,M_\odot$. Moreover, the proposed range of M_• 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 second 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 T_gas ∼ 10⁴ K. 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 10⁴ K. This happens because the flow becomes much more supersonic and the fraction of fragmenting gas decreases with increasing $\mathcal {M}$. 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 second 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 second 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 second stage of the collapse. Their Figure 4 clearly exhibits the dominant m = 2 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 disk-like 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.