Supermassive Black Hole Formation at High Redshifts via Direct Collapse in a Cosmological Context

We study the early stage of the formation of seed supermassive black holes via direct collapse in dark matter (DM) halos, in the cosmological context. We perform high-resolution zoom-in simulations of such collapse at high-$z$. Using the adaptive mesh refinement code ENZO, we resolve the formation and growth of a DM halo, until its virial temperature reaches $\sim 10^4$K, atomic cooling turns on, and collapse ensues. We demonstrate that direct collapse proceeds in two stages, although they are not well separated. The first stage is triggered by the onset of atomic cooling, and leads to rapidly increasing accretion rate with radius, from $\dot M\sim 0.1\,M_\odot {\rm yr^{-1}}$ at the halo virial radius to few $M_\odot \,{\rm yr^{-1}}$, around the scale radius $R_{\rm s}\sim 30$pc of the NFW DM density profile. The second stage of the collapse commences when the gas density takes precedence over the DM density. This is associated with the gas decoupling from the DM gravitational potential. The ensuing collapse approximates that of an isothermal sphere with $\dot M ( r )\sim $const. We confirm that the gas loses its angular momentum through non-axisymmetric perturbations and gravitational torques, to overcome the centrifugal barrier. During the course of the collapse, this angular momentum transfer process happens on nearly all spatial scales, and the angular momentum vector of the gas varies with position and time. Collapsing gas also exhibits supersonic turbulent motions which suppress gas fragmentation, and are characterized by density PDF consisting of a lognormal part and a high-density power law tail.


INTRODUCTION
galaxy co-evolution. Therefore, the formation of SMBHs should be considered alongside overall structure formation in the universe.
In addition to primordial SMBH seeds (e.g., Carr et al. 2010), several scenarios for early SMBH formation exist, such as growth from Population III remnants (e.g., Haiman & Loeb 2001;Abel et al. 2002;Bromm & Larson 2004;Yoo & Miralda-Escudé 2004;Li et al. 2007;Pelupessy et al. 2007;Tanaka & Haiman 2009), and collapse of stellar clusters (e.g., Devecchi & Volonteri 2009;Lupi et al. 2014). Population III progenitors, originally estimated to be as massive as ∼ 1000 M⊙, have been recently downsized by a factor of 10, due to radiation feedback-limited mass and fragmentation (e.g., Turk  To explain the detected high-z quasars, Pop III SMBH seeds are required to grow from ∼ 10 M⊙ to > ∼ 10 9 M⊙ in less than ∼ 7 × 10 8 yrs, which is uncomfortably close to the age of the Universe at this redshift. Is it plausible that super-Eddington accretion rates persist for ∼ 1 Gyr? Of course, Pop III remnants can lead to less massive SMBHs. On the other hand, models involving relativistic instabilities in stellar clusters must explain the origin of these clusters in the first place, and a substantial metal enrichment at these high redshifts.
Models, where early SMBHs seeds formed via direct collapse of gas into dark matter (DM) halos at z ∼ 10 − 20, provide an attractive alternative (e.g., Oh & Haiman 2002;Bromm & Loeb 2003;Haehnelt & Rees 1993;Volonteri & Rees 2005;Begelman et al. 2006;Wise et al. 2008;Regan & Haehnelt 2009;Begelman & Shlosman 2009;Milosavljević et al. 2009;Mayer et al. 2010;Schleicher et al. 2010;Hosokawa et al. 2011;Johnson et al. 2011;Prieto et al. 2013;Choi et al. 2013;Latif et al. 2013a,b). Gas collapse is triggered by atomic gas cooling, and occurs when the halo virial temperature surpasses Tvir ∼ 10 4 K. Several recent studies have dealt with the halo population hosting massive SMBH seeds (Prieto et al. 2013;Agarwal et al. 2013Agarwal et al. , 2014. Within the direct collapse framework, two alternative pathways have been proposed. First, a massive central object -a supermassive star (SMS) -forms at the center of the DM halo and is powered by a combination of core nuclear burning and Kelvin-Helmholtz contraction (e.g., Begelman et al. 2006Begelman et al. , 2008Begelman 2010). Following collapse of the stellar core and formation of the SMBH seed of ∼ 10 − few × 100 M⊙, depending on the angular momentum distribution (Begelman 2010), its convective envelope is powered by super-Eddington accretion onto this seedsuch a configuration has been termed a 'quasistar.' The seed SMBH can grow to ∼ 10 5−6 M⊙ in less than few Myrs. An important ingredient in the formation of the SMS is the very high gas accretion rate, > ∼ 0.1 M⊙ yr −1 (Begelman et al. 2006;Hosokawa et al. 2013;Schleicher et al. 2013). A basic ingredient of this model is trapping of the escaping energy and momentum within the SMS and the quasistar.
Recent work by Becerra et al. (2015) has used cosmological simulations with the AREPO moving mesh code (Springel 2010). Probably the most important assumption made was that the cooling rate of the collapsing gas has been artificially and exponentially suppressed at densities above 10 16 cm −3 . This led to a sharp increase in the gas temperature above these densities and the truncation of the gravitational collapse. As a result, the formation of the SMS has been supplemented by subsequent accretion from the surrounding disk, and by the disk fragmentation into low-mass stars.
According to the second pathway, gravitational collapse can retain a disky character even at the innermost scales, and, in tandem with the existence of a preferred channel for momentum and energy release (e.g., jets or winds), can bypass the SMS and quasistar stages, and the associated thermonuclear reactions (e.g., Begelman & Shlosman 2009;Choi et al. 2013).
Two caveats to the direct collapse scenario have been singled out because of their importance (e.g., Begelman & Shlosman 2009). First, the angular momentum barrier, in principle, can terminate the collapse well before it reaches ∼ 1 − 10 AU scale. Second, the gas could fragment, depleting the accretion stream by forming clumps and ultimately stars, and disturbing the accretion pattern. Choi et al. (2013) have performed a baseline study of direct collapse under idealized conditions, when the DM halo is isolated and the gas is in rotational equilibrium with a cosmological spin λ ∼ 0.05. They confirmed that gas within slowly tumbling DM halos of Mvir ∼ 10 8 M⊙ and Rvir ∼ 1 kpc can collapse to ∼ 10 AU scale. The collapsing flow overcomes the centrifugal barrier, losing its angular momentum by breaking axial symmetry and forming nested gaseous bars. Ultimately, its angular momentum is removed by gravitational torques from the gas and the DM. In the collapse phase, virial supersonic turbulence develops and fragmentation is damped. The gas accretion rate exceeds ∼ 1M⊙ yr −1 at various spatial scales, and allows for the formation of SMBH seeds at high redshifts, at least in principle.
In this work, we explore direct collapse in a full cosmological setup, and study the physical processes associated with its early stages within DM halos with virial temperatures of ∼ 10 4 K. We apply self-consistent zoom-in cosmological simulations. Our goal is to compare the gravitational collapse model in cosmological and isolated halos and to target the angular momentum and fragmentation problems. In Section 2, we explain the numerical details and the initial conditions. Sections 3.1 and 3.2 provide the results. The general DM halo properties are discussed in Section 3.1, and the dynamical aspects of the gas collapse are analyzed in Section 3.2. This is followed by discussion and conclusions.

Numerical Resolution
In this study, we use the Eulerian adaptive mesh refinement (AMR) code ENZO-2.3, which has been tested extensively and is publicly available (Bryan & Norman 1997;Norman & Bryan 1999;Bryan et al. 2014). ENZO uses a particle-mesh N -body method to calculate the gravitational dynamics, including collisionless DM particles, and a secondorder 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 lengthscale, if either the gas or DM density become greater than ρ 0,gas,dm N l , where ρ 0,gas,dm is the cosmic mean density for the gas or DM respectively, N = 2 is the refinement factor, and l = 35 is the maximal AMR refinement level. This refinement corresponds to a spatial resolution of ∼ 0.005 AU.

Zoom-in simulations
ENZO follows the non-equilibrium evolution of six species: H, H + , He, He + , He ++ , and e − Anninos et al. 1997) in a gas with a 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 discussed in Choi et al. (2013, and references therein), several physical processes have been suggested to prevent H2 formation, such as a very stong Lyman-Werner background radiation (e.g. Dijkstra et al. 2008;Ahn et al. 2009), Ly-α photon trapping (e.g. Spaans & Silk 2006;Choi et al. 2013), and the collisional dissociation in the shocked gas (e.g., Inayoshi & Omukai 2012). In this work, we neglect the H2 formation and destruction processes altogether, as well as exclude the chemistry and cooling related to H2. Understanding of H2 suppression mechanism can be important in estimating the population of the high-z SMBH seeds (e.g. Agarwal et al. 2014).
In this study, we are interested in the detailed dynamical evolution of the collapsing gas within a DM halo in the fully cosmological environment and subject to atomic cooling. To satisfy the resolution requirement, we use the MUSIC code (Hahn & Abel 2011) to generate the cosmological zoom-in initial conditions (ICs). MUSIC uses a realspace convolution approach in conjunction with an adaptive multi-grid Poisson solver to generate highly accurate nested density, particle displacement, and velocity fields suitable for multi-scale zoom-in simulations of structure formation in the universe. Generating a set of "zoom-in" ICs is a two-step process. First, we generate 1h −1 Mpc comoving 128 3 DM-only ICs for the pathfinder simulation and run it without AMR until z = 10. Using the HOP group finder (Eisenstein & Hut 1998), we select an appropriate DM halo, whose mass is ∼ 10 8 h −1 M⊙ at z = 10. Second, we generate 1 Mpc/h ICs with 512 3 resolution in DM and gas. Since we use the same random seeds of these ICs as the ICs at the first step, the phases of both ICs are identical. The zoomin region is centered on the selected halo position and is set to be large enough to cover the initial positions of all selected halo particles (see Figures 1 and 2). We set the DM particle smoothing length at 0.24h −1 pc in comoving coordinates. The ICs are generated using WMAP5 cosmology: ΩΛ = 0.721, Ωm = 0.279, Ω b = 0.0445, h = 0.701, σ8 = 0.807, and ns = 0.961. In the following we use R for spherical coordinates and r for cylindrical ones.

RESULTS
The computational box contains a number of DM halos whose virial temperatures exceed 10 4 K at z = 10. We follow a representative DM halo of our choice down to z ∼ 12 and observe its growth from mergers and accretion. Around t ∼ 350 Myr, the DM halo has reached the virial mass and radius of M h ∼ 2 × 10 7 h −1 M⊙ and Rvir ∼ 10h −1 kpc in comoving coordinates, and has acquired cosmological spin λ ∼ 0.03. The DM density profile is well approximated by the NFW profile (Navarro et al. 1997) with the characteristic radius of Rs ∼ 30 pc (in physical coordinates), beyond which the DM density profile steepens gradually to ∼ −3.
The DM halo concentration parameter is c ∼ 25. This halo and its environment are shown on various (comoving) spa-tial scales, from 250h −1 kpc down to 10h −1 kpc, in Figures 1 and 2, at the end of the simulation, t ∼ 360.13 Myr.
The filamentary structure of DM is evident on all scales in Figures 1 and 2. The gas distribution follows that of the DM. The targeted halo becomes dominant on the smallest scales shown (bottom frames of Fig. 2). It is connected to the DM web via three filaments whose width is nearly comparable to the halo virial diameter.   Figure 3. Evolution of spherically-averaged gas temperature (left) and density (right) profiles. The time in the legend shows the age of the universe. The initial profile corresponds to the time when the halo virial mass is sufficient to trigger the atomic cooling in the gas. The last profile corresponds to time when the gas collapse has reached down to ∼ 10 AU scale. The x-axis provides the distance from the densest cell, and is measured in physical units. The density profiles confirm that the halo gas experiences nearly isothermal central runaway collapse.
The most important feature of the target DM halo is its triaxiality, i.e., its three major axes all differ from each other, in all three projection planes, xy, yz, and xz. This is not suprising, as DM halos are universally triaxial as they form in numerical simulations (e.g., Allgood et al. 2006) before dissipative processes axisymmetrize them in the subsequent evolution (e.g., Berentzen et al. 2006;Shlosman 2007). Moreover, the rotation of the halo figure, i.e., its rate of tumbling, is extremely slow, and seems to be the general property of DM halos, as pointed out by Romano-Díaz et al. (2009).

Properties of the collapsing gas in DM halo
Since we are interested in details of the dynamical process within this halo, we switch to physical units. At the early stage of halo growth, the gas mostly follows the DM assembly. This means that it is accumulating within the growing halo, and, during the quiescent time periods of no major mergers, the gas is largely hydrostatic with a small degree of the rotational support. Of course 'hydrostatic' has a very approximate meaning here, as the gas joins the halo partly via penetrating filaments which result in large-scale motions within the halo, i.e., the streamers. The halo gas stops increasing its temperature and starts to cool via atomic cooling, when the halo becomes massive enough and its virial temperature has reached ∼ 10 4 K. Note that we ignore the H2 cooling and the associated chemistry, and implement only atomic cooling (see Section 2.2).
Cooling allows the gas to be driven into the gravitational potential minimum. This is demonstrated in Figure 3 which shows the evolution of gas temperature and density profiles within the halo. When Tvir reaches 10 4 K, atomic cooling becomes important and collapse is triggered. Low-T gas continues to be accreted from outside Rvir along the filaments as well as smooth accretion from arbitrary directions.
The smooth accretion experiences a shock at R ∼ 800 pc, which virializes it. The gas being accreted along the filaments virializes about a decade deeper in R. The collapsing gas roughly maintains isothermality at the cooling floor, with a small radial decline in T , as the heating appears to be inefficient. The gas collapse leads to the establishment of an isothermal density profile ρ ∝ R −2 inside Rs. Finally, the collapse reaches ∼ 10 AU scale.
The collapse clearly proceeds from outside in (Figure 3, right frame). Figure 4, which shows the gas-to-DM density ratio profiles at various times, reveals the critical detail of this collapse. At R > ∼ 100 pc, the gas closely follows the DM density profile, and the baryon-to-DM ratio is at about its cosmic average. Once the gas temperature is reduced below the virial temperature, this ratio increases in the inner halo, and eventually exceeds unity. The dotted-gashed line at ∼ 355 Myr corresponds to the time when the inner gas density, inside ∼ 5 − 10 pc, exceeds that of the DM. After this stage, the inner gas density rapidly increases and establishes the isothermal density profile -as the second stage of the collapse develops, with a very high inflow rate. The gas essentially decouples from the background DM potential. Similar behavior has been found in the evolution of isolated models (Choi et al. 2013). Figure 5 shows the development of the gas inflow rate. From Rvir and down to a ∼ 10 pc scale (just inside Rs), it increases fromṀ ∼ 0.03 M⊙ yr −1 toṀ ∼ few M⊙ yr −1 , i.e., by nearly two orders of magnitude. Inside Rs,Ṁ becomes approximately constant with radius. The dramatic increase inṀ is a reflection of the steep DM density distribution outside Rs and associated gas distribution, with a logarithmic slope varying from ∼ −3 down to ∼ −2 at Rs. As long as this DM profile persists, the mass flux into the inner halo ( < ∼ Rs) will not change: the DM distribution creates a kind of bottleneck for the gas supply rate. Only growth in the DM virial mass will allow a higherṀ outside Rs, but the growth   . Spherically-averaged gas-to-DM density ratio profiles. The time in the legend gives the age of the universe. Initially, the DM and the gas density profiles exhibit very similar shapes reflecting the cosmic mean where ρgas/ρ DM ∼ 0.17. The cooling allows the gas to collapse, increasing the ratio, and to reach the DM density at Rs ∼ 30 pc. This radius movies only slightly inward with time. The outer region, r ∼ 200 − 500 pc becomes gas deficient (compared to the cosmological mean) as the gas inflow across r vir cannot replenish the collapsing gas on a short timescale.
of the DM halo proceeds on a timescale much longer then the gravitational collapse at its center. An interesting corollary of this effect is that after the halo gas has collapsed to the center, the DM halo will be largely depopulated of gas. The radial (infall) velocity and andṀ actually decrease from Rvir (800 pc) to ∼ 100 pc, where they have either global or local minima (see the upper frame of Figure 5 and Figure 6). As the gas moves toward R ∼ 100 pc, it experiences an increase in density, R −2.5−3 , which is offset by a sharper decrease in vR. To fully understand this behavior ofṀ (R) and vR(R), one should also note a simultaneous increase in the tangential velocity at the same radii -this means that the gas angular momentum becomes more significant. Inside the minimum at R ∼ 100 pc, both the density and the radial velocity increase inward, which results in the sharp increase ofṀ .
Inside Rs,Ṁ reaches ∼ few M⊙ yr −1 , and stays approximately constant down to R ∼ 10 −4 pc. Rs is roughly the radius at which the gas density approaches the DM density, and the gas decouples from the background DM potential. The value ofṀ can be estimated from the isothermal spherical gas collapse model, ∼ (Mgas/Mtot)v 3 ff /G, where the free-fall velocity, v ff , is a measure of the local gravitational potential, and Mgas/Mtot is the gas-to-total mass ratio within R (Choi et al. 2013). The virial temperature inside Rs is Tgas ∼ 8000 K, which gives the observed value ofṀ ∼ few M⊙ yr −1 . Note that we continue this simulation until the collapse has been established all the way to R ∼ 10 −4 pc, the radius where we estimate that the optical depth for the radiation produced internally by the gas will reach unity (e.g., Choi et al. 2013  The R-axis is the distance from the densest cell measured in physical units. The accretion rate,Ṁ, increases from the virial radius to the inner ∼ 10 pc scale, which reflects the increase in the free-fall velocity within the NFW DM potential.Ṁ saturates at around ∼ few M ⊙ yr −1 in the inner halo, which extends from r ∼ 10 pc (i.e., just inside ∼ Rs) all the way inward. The location of the Rs is indicated with the black arrow.
this radius and Rs. When the collapse is allowed to proceed further, we expect that the virial temperature of the decoupled gas will exceed the halo virial temperature by a large factor (even if the gas density profile flattens), as the gas will determine the depth of the potential well. The actual temperature of the optically-thick collapsing gas, of course, depends on radiative transfer effects (and other cooling mechanisms), first for the bound-bound transitions and then for the continuum. This rapid gas accretion is one of the key differences between the physical conditions during star formation and SMBH seed formation (Begelman et al. 2006;Begelman & Shlosman 2009). The cosmological simulation in this paper demonstrates that the gas collapse, which is facilitated through atomic cooling, reproduces the central runaway with highṀ and without significant fragmentation. Our results show that the gas collapse proceeds roughly in two stages, reflecting the shape of the gravitational potential dominated by the DM initially and by the gas thereafter.

Dynamics of the collapsing gas
Evolution of radial and tangential velocities in the collapsing gas is shown in Figure 6. Note that the highest velocities are achieved inside the central pc, where the gas has decoupled from the background DM potential and increases its virial velocity above that of the DM. The tangential velocity approaches ∼ 20 km s −1 and so does its radial counterpart -a clear indication that the gas has substantial rotational support, albeit sub-Keplerian, as it continues to collapse. We shall return to this issue again. Figure 7 displays gas density slice maps at the end of the simulation, on four characteristic (physical) scales, from 1 kpc down to 10 −4 pc, on the xy, xz and yz planes. The  filamentary gas distribution, which can be observed on the largest scales, has its origin in the DM distribution. The top row exhibits the overall environment of the growing halo, on scales of ∼ 1 kpc. Multiple filaments which fuel the accretion can be observed on this scale. The second row (∼ 1 pc) shows the inner halo where the gas density becomes higher than the DM density. The third and last rows, ∼ 0.01 − 10 −4 pc, display regions fully dominated by the gas and very asymmetric in distribution and filamentary. the last row, at ∼ 200 AU, represents the environment of the expected optically thinto-thick transition in the collapsing gas. The continuity of the filamentary structure to small radii hints that its origin does not lie in shocks but rather that this is inherently a signature of a cosmological accretion flow before it virializes. Note that we stop the simulation very early, but already at this time we see the non-axisymmetric gas response to a non-axisymmetric background potential completely dominated by the DM in the top row frames. Already in the second row of frames, the gas potential dominates, but the DM on larger scales can still provide the finite amplitude non-axisymmetric perturbation on smaller scales 1 , assisted by the non-axisymmetric distribution of the gas on all radii. The gas density slices in the bottom row frames clearly show the presence of a nonaxisymmetric feature centered on the highest gas density cell. The apparent presence of low-level Fourier harmonics is verified by the mode analysis of the gas density described below. These non-axisymmetric features play a central role in the transfer of angular momentum outward by gravitational torques on small scales -a process that allows for the continuous gas collapse.
Despite the fact that the halo has only a small fraction of rotational support (λ ∼ 0.03), as does the gas, the centrifugal support of the gas would increase quickly if angular 1 In a non-axisymmetric density distribution, the material at larger radii exerts gravitational torques on the interior material. momentum were conserved during the collapse. Without efficient transfer of angular momentum, the collapse would be halted when the angular momentum reached its Keplerian value. In the presence of the background DM, this corresponds to collapse by a factor of 10 in radius from the largest scales (e.g., Shlosman 2013). As gas does not accumulate at any radius in our simulation (although it does slow down its radial motion at various radii), it is clear that it does not reach the centrifugal barrier, and, therefore, its angular momentum is not conserved. That non-axisymmetric perturbations can facilitate angular momentum transfer in steady state systems is well known -these can be spiral arms (e.g., Lynden-Bell & Kalnajs 1972;Tremaine & Weinberg 1984), large-scale bars or a hierarchical bars-in-bars structure (e.g., Shlosman et al. 1989;Shlosman 2005).
While spontaneous bar instability typically requires a few rotation periods to develop, bars and bars-in-bars can also be triggered by a finite amplitude perturbation, e.g., such as provided by triaxial DM halos (Shlosman 2011). Strong gravitational torques which accompany finite amplitude perturbations can transfer angular momentum on the short dynamical timescale encountered in direct collapse, and are more efficient as their rise time is negligible. Such torques can follow from the low-m Fourier non-axisymmetric modes, like m = 1 and m = 2, which can be associated with the displacement of the center of mass of the gas with respect to the DM, and/or the development of nested gaseous bars (Shlosman et al. 1989(Shlosman et al. , 1990Englmaier & Shlosman 2004;Begelman & Shlosman 2009). We, therefore, analyze the prevailing non-axisymmetric modes arising in our simulations (e.g., Figure 7). Although Figure 7 clearly indicates the existence of nonaxisymmetric density features in the inner halo, it is not clear what mode dominates the structure. We, therefore, perform a Fourier analysis for two density modes, m = 1 and 2 -the fastest growing non-axisymmetric modes (see e.g., Long et al. (2014) for technical details). The m = 1 mode requires displacement of the center of mass of the system. Its presence is evident in the fact that the densest cell in the simulation separates from the center of mass of the collapsing gas, measured within a sphere with a radius of 0.08 pc. The m = 2 mode is a barlike mode and usually plays the dominant role in angular momentum transfer. Figure 8 shows the evolution of the Fourier amplitudes of m = 1 and 2 modes normalized by the amplitude of the m = 0 mode. The mode evolution shown starts after the time when atomic cooling has been triggered in the halo, ∼ 310 Myr. We display the mode evolution only for one plane, but it is representative of the overall behavior. The DM halo shape is that of triaxial ellipsoid and exerts gravitational torques on the gas. The latter responds to this finite amplitude perturbation in a nonlinear fashion and develops m = 1, 2 and higher modes.
The top panel in Figure 8 shows the mode evolution for the cylindrical annulus defined by the minimal and maximal radii of rmin = 10 pc and rmax = 50 pc in the xy-plane, with the vertical thickness of the slice ∆z = 5 pc. This region is dominated by the DM, and the behavior of both modes is similar to that of the isolated halo studied by Choi et al. (2013), although the amplitude of both modes is higher. Before the gas decouples from the DM, the m = 1 mode dominates over the m = 2 mode. However, around the start of the central runaway, ∼ 355 Myr, the m = 2 mode starts to dominate over m = 1. The bottom panel shows mode evolution for a smaller annulus defined by rmin = 0.01 pc, rmax = 0.1 pc and ∆z = 0.01 pc, respectively. This scale is dominated by the gas. It shows only the last stage of the central runaway, where m = 2 dominates over the m = 1, until the very last moment where both of them increase dramatically, reflecting the formation of the elongated structure -the nonlinear response of the very central gas. Comparing the large and small scale evolution of these modes, we note that last rise of the m = 1 mode is important only at the very center. It has a much smaller amplitude on the larger scale, indicating that the cause for its rise lies in the off-center motion on scales < ∼ 0.1 pc.

Collapse and angular momentum
The angular momentum axis of the collapsing gas varies with R -an effect discussed by Romano-Díaz et al. (2009, see their Fig. 19) for the growth of DM halos in the cosmological context. It reflects the variability in the angular momentum axis of the continuous gas inflow. Since the observer's frame is fixed in its position, the gas moving to smaller radii is replaced by gas inflow from larger radii. This 'replacement' gas has typically a different orientation of the angular momentum axis, and can in principle drive variability of the Fourier density mode amplitude. Figure 9 provides two snapshots of the specific angular momentum profile of the gas, jgas(R). These profiles are superposed on the circular specific angular momentum jc(R). At each radius, jc is the maximum possible angular momentum allowed for the bound gas at a fixed energy. The total mass included within a sphere of radius R is about linear with R in the halo, the circular velocity is about constant, and jc ∼ R. jgas(R) is typically smaller than jc, by a factor that varies with R and with time. This means that the collapse is typically not prohibited by the angular momentum barrier, for a considerable range in radii. At smallest radii, it lingers not far from the centrifugal barrier, at each time. This can happen only if there is a constant flow of angular momentum away from the gas as it moves inward.
As the DM distribution is triaxial, it exerts gravitational torques on the gas at all radii, as we noted above, but its contribution is gradually washed out in the region where the gas dominates. The gas is losing its angular momentum at a rate of djgas/dt ∼ τ , where τ is the torque per unit mass. It depends on the offset in the position angle between the gas and DM, i.e., on the asymmetric parts of the DM and gas density distributions, and the asymmetric part of the background gravitational potential (e.g., Berentzen et al. 2007). The angular momentum, therefore, can flow from the inner gas to the outer gas, and from the gas to the DM.
The efficiency of angular momentum extraction by gravitational torques can be measured by η ≡ 1 − jgas/jc, i.e., by the 'separation' between the centrifugal barrier and the actual angular momentum in the gas at each R. The right frame of Figure 9 displays the typical situation in the later stages of the gravitational collapse. The efficiency η is about constant for a wide range in radii. Note, that if jgas moved away from jc as the gas moves in, η would increase and the residual jgas would be less important dynamically. On the other hand, if jgas decreased more slowly than jc, the angular momentum would become more important. The fact that η ≈ const., means that there exists a tight balance between the characteristic timescale of j loss by the gas and the inflow timescale.
The gas distribution is substantially asymmetric on all scales (Figure 9) and is subject to strong gravitational torques. The bottom panel of Figure 8 confirms that the m = 1 and 2 modes are the key contributors to the transfer of the angular momentum. Figure 8 and scrupulous inspection of the full evolutionary feature of the angular momentum profile at various times confirm that the angular momentum transfer happens multiple times on multiple spatial scales. The figure rotation of the DM halo can be neglected, as it tumbles extremely slowly (Romano-Díaz et al. 2009). The outer torques are dominated by the DM distribution, while at smaller radii, it is the nonaxisymmetric shape of the gas distribution that provides the torques.
Unlike in the case of collapse from idealized conditions within an isolated halo, the cosmological collapse, expect- edly, does not show the formation of a disk, because rotational support never becomes high enough and because the vector of the specific angular momentum exhibits temporal and spatial variability.

Collapse and turbulence
Figure 10 displays slices of the collapsing gas with the gas velocities and the degree of turbulence given by the magnitude of the vorticity, defined by w = ∇ × v, where v is the velocity field. The evolutional time, spatial scales, and viewing angles are the same as in Figure 7. The arrows represent the direction of gas flow and their sizes give the relative flow speeds. The top panels exhibit the gas motion on the scale of the halo, 1 kpc. The middle panels show the gas flow on scales where it decouples from the DM -rotational motions are clearly visible and so some degree of rotational support is present. Similar phenomena can be observed on the smaller scales (bottom panels). Figure 10 demonstrates the degree of turbulence. The collapsing halo gas develops the supersonic turbulent motions (e.g., Wise et al. 2008;Begelman & Shlosman 2009;Regan & Haehnelt 2009), which suppress fragmentation (Begelman & Shlosman 2009;Choi et al. 2013). The supersonic turbulence works to both damp and trigger fragmentation. The damping is provided via turbulent pressure, above the level of the thermal pressure, thus acting against the self-gravity. On the other hand, shocks associated with supersonic turbulent motions induce fragmentation (e.g., Krumholz & McKee 2005). The fragments, however, must collapse before the passage of the next shock front, which otherwise will destroy the fragments. We also note that fragmentation in the DM-dominated phase is suppressed by the DM background, because its action dilutes the gas selfgravity, and, therefore, increases its Jeans mass -the gas cannot collapse until its density surpasses that of the background DM (Choi et al. 2013).
In Figure 10, the smaller scales exhibit larger vorticity magnitude than larger scales (note the shifting colors in the color palettes). This implies that the turbulent motions increase and will continue to increase with the gas collapse, as the potential well generated by the gas deepens. This is consistent with Figure 6: vR and vt increase at smaller R. Figures 7 and 10 show that no major gas fragmentation occurs. A central density maximum is well-identified at all times. Although the density slices do not display major fragmentation, multiple shocks are present as the analysis shows. These shocks result from the supersonic turbulent motion of the collapsing halo gas. The existence of the shocks without major fragmentation suggests that incipient fragments in the collapsing halo gas are destroyed by the next incoming shock, before the fragmentation proceeds into a strongly nonlinear regime. If fragmentation occured, it would deplete the available gas supply to the center and would disturb the developing low-m Fourier density modes which facilitate the angular momentum transfer. More quantitative analysis of this is given elsewhere.
Our simulation has been stopped when the collapse reaches the required scale of ∼ 10 −4 pc, which happens when the refinement level has reached 35. In order to compare directly with the simulations of Choi et al. (2013), we have restarted the simulation when it reaches the comparable re-  Figure 11. Evolution of the volume-averaged, gas density PDF as a function of log 10 ρ measured at the end of the simulation, at t ∼ 360.13 Myr, and sampled with > ∼ 10 6 AMR cells at the fixed refinement level (see text). The sampling shows the PDFs of the central sphere of radius 200 AU as blue histogram. Shown also are the lognormal fit (thick blue solide line) presented in Equation 1 with σ ∼ 1.17M. The average density for the sampled sphere is ρ = 1.76 × 10 −12 g cm −3 (for the blue histogram). It also shows the continuous evolution of the high density side, ρ > 10 −10 g cm −3 , of the blue histogram at consecutive times separated by ∼ 10 yrs as green dotted, red dot-dashed, and cyan dashed lines. The evolution shows that the lines have reached the slope of ∼ −1.18 in the power law tail. The collapsing gas has been sampled at the resolution of ∼0.7 AU and the density fluctuations extend over 7 decades. The refinement levels for these figures is kept at 28, in order to compare with identical conditions in Figure 15 of Choi et al. (2013).
finement level of 28. We sample the gas density cells within a distance of 200 AU from the center of the collapse with > ∼ 10 6 cells, and fit the PDF (by least squares) to a pure lognormal distribution with σ ∼ 1.17M (blue histogram) from Equation 1. We have also examined several different sampling scales and locations and obtained similar results. Figure 11 confirms that the obtained blue histogram PDF is nicely fit with the lognormal distribution, without any visible power law tail. We, therefore, continue the simulation for additional time at a fixed refinement level. Figure 11 shows the evolution of the density PDF and one can clearly observe the formation of a progressively shallow power law tail. The terminal slope (cyan dashed line) is about −1.18. As the power law tail appears at high densities, its origin must be related to the onset of self-gravity as the gas starts to pile up in the center. However, the power law slope attained in this simulation is still evolving, and its terminal value most probably would steepen as the collapse proceeds.

DISCUSSION
In this paper, we have investigated some aspects of direct collapse that can lead to the formation of SMBH seeds at high redshifts. We have used the cosmological zoom-in initial conditions and an Eulerian AMR code to resolve the gas and DM dynamics inside and in the vicinity of the targeted DM halo. We have shown that the gas atomic cooling in the DM halo will trigger central runaway collapse without significant fragmentation. The central runaway can be divided into two stages: outer collapse in the DM-dominated potential, where the density flattens gradually from a log slope of −3 to −2 at a few pc from the center, i.e., the NFW profile, and inner collapse at smaller R. The characteristic radius where the second stage is triggered corresponds roughly to the NFW scale radius Rs, where the background DM density gradually becomes shallower and reaches the (log) slope of −2.
The associated mass accretion rate reaches a few M⊙ yr −1 at this radius. The second stage of the collapse represents the gas-dominated region, where the gas decouples from the DM background potential and continues its collapse. As the physical conditions in the central regions require the introduction of radiative transfer on-the-fly, due to the buildup of optical depth in the hydrogen lines, we have terminated the collapse in the early stage when it reaches the required resolution scale of 10 −4 pc. We confirm the Choi et al. (2013) finding that the lowm Fourier non-axisymmetric modes are responsible for the angular momentum flow outward which results in overcoming the angular momentum barrier at various radii. The specifics of direct collapse in the cosmological context involves the variability of the angular momentum vector as a function of radius and time. Our analysis of the mode and angular momentum evolution demonstrates that gravitational torques clearly dominate over the hydrodynamical torques over a wide dynamic range during the collapse. Efficient loss of angular momentum allows for continuous collapse over 7 decades in radius, compared to one decade if the angular momentum were conserved in the gas. We also demonstrate that the collapsing flow develops supersonic turbulent motion. The degree of the turbulence increases as the collapse proceeds to small radii, and the supersonic turbulent motion suppresses gas fragmentation. Latif et al. (2013a) have argued that turbulent eddies will be under-resolved at small R, and have added the turbulent driving on all scales. However, we achieve the highest resolution at the smallest scales and have tested both the rms velocities there, which appear to be supersonic, and the density PDF, which extends over 7 decades in density. We find no indication that the turbulence is unresolved there. One can compare the final state of the simulation achieved with and without turbulent driving. Indeed, the driving leads to a more regular rotation on small scales, but one questions whether this is realistic. This issues must be studied further.
Simulations of direct collapse in the cosmological context have been also performed by Prieto et al. (2013), but their spatial resolution has been limited to ∼ 1% of the halo virial radius. This is insufficient to resolve the second stage of the gravitational collapse, where the gas decouples from the DM background.
Additional processes are known to help to suppress fragmentation in the collapsing gas. The most important process is the dissociation of H2 by UV background continuum from external sources. It produces the temperature floor and maintains isothermality with T ∼ 10 4 K. We assume that the UV background radiation prevents H2 formation.
There are several recent studies arguing that this is a plausible case, when the Lyman-Werner background continuum is present at sufficient levels (e.g., Omukai 2001;Shang et al. 2010;Regan et al. 2014). This radiation background can in principle constrain the population of massive SMBH seeds at high-z.
Another issue is related to the shrinking numerical timestep which allows us to follow the evolution of the central regions but basically freezes the outer regions. A possible solution to circumvent the small timestep is to impose a sink particle mechanism for a given resolution level (e.g. Latif et al. 2013b). In this case, the spatial resolution is limited by the scale imposed for the sink particles. The simulation can continue further and one can follow the longterm evolution of the direct collapse. This approximation can provide some details of the evolution of the central objects (Shlosman, Choi, Begelman & Nagamine 2015 in preparation).
Given that we have confirmed that the central runaway collapse can proceed without significant fragmentation, we conclude that high-z SMBH seeds can form, at least in principle, through a process of direct collapse within a DM halo. SMBH seeds formed in this way, while relatively massive, are not expected to follow the low-redshift M-σ relation. This suggests that the co-evolution of massive SMBH seeds and their host protogalaxies may follow a very different evolutionary path, particularly at high-z, than their later counterparts.