TY - JOUR AU - Begelman, Mitchell, C AB - ABSTRACT We have modelled direct collapse of a primordial gas within dark matter haloes in the presence of radiative transfer, in high-resolution zoom-in simulations in a cosmological framework, down to the formation of the photosphere and the central object. Radiative transfer has been implemented in the flux-limited diffusion (FLD) approximation. Adiabatic models were run for comparison. We find that (a) the FLD flow forms an irregular central structure and does not exhibit fragmentation, contrary to adiabatic flow which forms a thick disc, driving a pair of spiral shocks, subject to Kelvin–Helmholtz shear instability forming fragments; (b) the growing central core in the FLD flow quickly reaches |${\sim } 10\, \text{M}_\odot$| and a highly variable luminosity of |$10^{38}\text{--}10^{39}\, {\rm erg\, s^{-1}}$|⁠, comparable to the Eddington luminosity. It experiences massive recurrent outflows driven by radiation force and thermal pressure gradients, which mix with the accretion flow and transfer the angular momentum outwards; and (c) the interplay between these processes and a massive accretion, results in photosphere at ∼10 au. We conclude that in the FLD model (1) the central object exhibits dynamically insignificant rotation and slower than adiabatic temperature rise with density; (2) does not experience fragmentation leading to star formation, thus promoting the fast track formation of a supermassive black hole (SMBH) seed; (3) inclusion of radiation force leads to outflows, resulting in the mass accumulation within the central 10−3 pc, which is ∼100 times larger than characteristic scale of star formation. The inclusion of radiative transfer reveals complex early stages of formation and growth of the central structure in the direct collapse scenario of SMBH seed formation. methods: numerical, galaxies: formation, galaxies: high-redshift, cosmology: theory, dark ages, reionization, first stars, quasars:supermassive black holes 1 INTRODUCTION Supermassive black holes (SMBHs) of |$\mathrel {\raise.3ex\hbox{}\gt {\rm }\lower0.6ex\hbox{}{\sim }{\rm }}10^9\, \text{M}_\odot$| are increasingly found at high redshifts, |$z\mathrel {\raise.3ex\hbox{}\gt {\rm }\lower0.6ex\hbox{}{\sim }{\rm }}6$| (e.g. Fan et al. 2003; Willott et al. 2010a; Mortlock et al. 2011; Wu et al. 2015), up to |$z$| ∼ 7.54 at present (e.g. Venemans et al. 2017; Banados et al. 2018). They reside in very luminous active galactic nuclei (AGN) and appear to form the extreme of the overall population of such objects. While two or three specific alternatives have been explored for explaining the development of such massive SMBHs at these redshifts, the broader issue of SMBH formation is equally important for our understanding of structure evolution in the universe and galaxy evolution (e.g. Shlosman 2013, for review). The main difficulty in forming the SMBHs in galaxy centres is their long growth time, if the initial seed mass is small compared to the final product. At high redshift, the only realistic options include direct collapse leading to massive seeds (e.g. Haehnelt & Rees 1993; Bromm & Loeb 2003; Begelman, Volonteri & Rees 2006; Wise, Turk & Abel 2008; Begelman & Shlosman 2009; Milosavljević et al. 2009; Regan & Haehnelt 2009; Schleicher, Spaans & Glover 2010; Hosokawa et al. 2011; Choi, Shlosman & Begelman 2013, 2015; Latif et al. 2013a,b; Shlosman et al. 2016), or remnants of Population III stars, which, if fed vigorously and merged, can form less massive SMBHs of Seyfert galaxies, inprinciple. Conditions for direct collapse with a gas of pristine composition are relatively simple. Namely, if the gas can form molecular hydrogen and cool below ∼1000 K, it will collapse into dark matter (DM) minihaloes of |${\sim } 10^6\, \text{M}_\odot$|⁠, while if the gas remains atomic, it requires DM haloes with virial temperatures reaching ∼104 K. In the former case, the process is expected to lead to the formation of a Pop III star or a few stars (e.g. Abel, Bryan & Norman 2002; Bromm, Coppi & Larson 2002; O’Shea & Norman 2007; Turk, Abel & O’Shea 2009). In the presence of radiative feedback, a rather normal initial mass function (IMF) has been found (e.g. Hosokawa et al. 2011; Hirano et al. 2015; Hosokawa et al. 2016). In the latter case, the result can be the formation of a supermassive star (SMS; Begelman et al. 2006; Begelman 2010), which evolves to an SMBH seed after its core collapse – a so-called quasi-star that accretes in the hyper-Eddington regime. Supersonic streaming velocities – remnants from the recombination epoch – can suppress formation of Pop III stars, resulting in more massive central objects (Tanaka, Li & Haiman 2013; Hirano et al. 2017). The quasi-star mass has been argued to lie in the range of |${\sim } 10^{1}\text{--}10^{6}\, \text{M}_\odot$|⁠, or even higher (Begelman, Rossi & Armitage 2008), if fragmentation can be suppressed and the angular momentum can be efficiently transferred outwards (Begelman & Shlosman 2009; Begelman 2010; Choi et al. 2013, 2015). Gravitational torques assisted by shocks have been verified to transfer the angular momentum of the collapsing gas to the DM and the outer gas (Choi et al. 2013, 2015). Supersonic turbulent motions, induced by the torques and resulting shear, damp fragmentation in the atomic gas. This evolution differs from that described by the self-similar analysis, which was necessarily limited to a linear stage (Hanawa & Matsumoto 2000); the growing bar-like m = 2 mode in its nonlinear stage does not lead to fragmentation, but induces gas inflow. Alternatively, it is possible that the stellar evolution stage can be by-passed completely, e.g. if the temperature does not rise sufficiently high to trigger thermonuclear reactions (e.g. Begelman & Shlosman 2009; Choi et al. 2013; Shlosman et al. 2016). Direct collapse within DM haloes has been investigated in the optically thin regime, allowing the gas to cool down efficiently. Difficulties in performing on-the-fly radiative transfer in the collapsing gas have motivated models that switch to an adiabatic equation of state, by cutting off cooling within a region having a density above some critical value (e.g. Becerra et al. 2015, 2018). These models lead to a rotationally dominated flow in the inner region – a disc and fragmentation due to Jeans instability. However, this approximation neglects a long list of processes that operate in the central region, such as radiation pressure and its force, and artificially suppresses cooling in regions where it should be able to operate. In a recent paper (Luo et al. 2018), we have dealt with the optically thick stage of direct collapse within isolated DM haloes, including radiative transfer in the flux-limited diffusion (FLD) approximation. Furthermore, we have compared the FLD models with models assuming an adiabatic equation of state, and found that their evolutions diverge. The main goal of this work is to model the optically thick stage of direct collapse with radiation transfer in the FLD regime, in a fully self-consistent cosmological framework. The formation of the central object under direct collapse has not been simulated so far, except under simplified adiabatic conditions, with the exception of Luo et al. (2018) that dealt with isolated DM haloes. In this paper, we take an additional step towards realistic modelling of the outcome. If the collapse leads to an SMS, we wish to determine its properties, namely, its mass, size, and spin. What opacities dominate in its interior and near its photosphere? If it is rotation-dominated, does it rotate differentially, or more like a solid body? How does the massive accretion flow affect its internal structure, including convection, and its thermal and dynamical relaxation? What is its central temperature and is it sufficiently high to trigger thermonuclear reactions? Are outflows associated with this stage and what drives them? Inclusion of radiation transfer should in principle answer these questions. This paper is structured as follows. The next section describes the numerical code used in our simulation, as well as the details of radiative transfer implemented here, and the initial cosmological conditions used. Section 3 presents our results for adiabatic flows, and Section 4 for FLD flows. This is followed by a discussion section and conclusions of this work. 2 NUMERICAL TECHNIQUES For simulations of direct collapse within DM haloes, we invoke the modified Eulerian adaptive mesh refinement (AMR) code enzo-2.4 (Bryan & Norman 1997; Norman & Bryan 1999). enzo implements a particle-mesh N-body method to calculate the gravitational dynamics, to follow collisionless DM particles, and a second-order piecewise parabolic method to solve hydrodynamics (Colella & Woodward 1984; Bryan et al. 1995). Supplementary inner meshes are allowed by the structured AMR, in order to enhance the resolution in the pre-specified regions. The number of rectangular grids that cover some region at a given refinement level, and the number of refinement levels are not subject to limitation (Berger & Colella 1989). When densities in the DM or gas exceed ρ0Nl, the simulation grid is refined by a factor of 2 in length-scale, where N = 2 is the refinement factor, l is the maximal AMR level of refinement, and ρ0 is the threshold density for refinement. The force resolution corresponds to twice the minimal cell size in adaptive PM codes (e.g. Kravtsov, Klypin & Khokhlov 1997) The Truelove et al. (1997) condition to resolve the Jeans length, i.e. to have at least four cells per Jeans length, is exceeded, in order to avoid spurious fragmentation. We have imposed the condition of at least eight cells per Jeans length in our simulations, in line with recent numerical experiments requiring even higher numerical resolution in order to properly resolve the turbulent motions (e.g. Sur et al. 2010; Federrath et al. 2011; Turk et al. 2012; Latif et al. 2013a). 2.1 Radiation hydrodynamics and radiative transfer formalism We use the FLD approximation in order to model radiation transport. Local thermodynamic equilibrium (LTE) is imposed in optically thick regions of the flow, where the Planck intensity is used for the gas emissivity, and the Saha equation determines the gas ionization level (e.g. Rybicki & Lightman 1979). In enzo-2.4, in order to solve the radiation transport equation, a fully implicit inexact Newton method has been applied and coupled to the AMR cosmological hydro solver with an explicit, operator-split algorithm, but only at the end of the top level timestep (Reynolds et al. 2009). We have modified this by updating each level of refinement at the end of respective timestep. This makes the FLD fully consistent with the hydrodynamics (Luo et al. 2018). In Luo et al. (2018) and here, we have introduced the radiation force and |$v$|/c order terms, with c and |$v$| being the speed of light and the gas velocity, respectively. In this work, we have included the cosmological terms, with a being the cosmological expansion parameter. The modified Euler equation is \begin{eqnarray*} \frac{\partial \rho {\boldsymbol v}}{\partial t} + \frac{1}{a} \nabla \cdot (\rho {\boldsymbol v} {\boldsymbol v} + \mathbb {I} p) = - \frac{\dot{a}}{a} \rho {\boldsymbol v} - \frac{1}{a}\rho \nabla \phi + \frac{1}{a}\frac{\kappa _{\rm R}}{c}{{\boldsymbol F}}. \end{eqnarray*} (1) Here, p, ρ, and |${{\boldsymbol F}}$| are the thermal pressure, comoving baryon density, and the radiation energy flux, respectively. |$\mathbb {I}$| is the identity matrix, and κR is the Rosseland mean opacity (Section 2.2). The gravitational potential ϕ is calculated from the DM density ρDM and the baryon density ρ. The radiative flux vector in the FLD approximation can be presented by the gradient of radiation energy density (Levermore 1984), i.e. using the Fick’s diffusion law, \begin{eqnarray*} {{\boldsymbol F}} = - \frac{c\lambda _{\rm F}}{\kappa _{\rm R}}\nabla E, \end{eqnarray*} (2) where λF is the flux limiter, |$\lambda _{\rm F} = \lambda _{\rm F}(E, \nabla E \kappa _{\rm R}) = (9 + R_\lambda ^2)^{-1/2}$|⁠. The auxiliary function Rλ is defined as Rλ = |∇E|/(κRE). The comoving radiation energy density in a cosmological medium, omitting the frequency-dependence by integration over the radiation energy spectrum, is given by \begin{eqnarray*} &&\frac{\partial E}{\partial t} + \frac{1}{a} \nabla \cdot (E {{\boldsymbol v}})\nonumber\\ &&\quad= - \frac{1}{a^2}\nabla \cdot {{\boldsymbol F}} - \frac{\dot{a}}{a} E - {\boldsymbol P}:\nabla {{\boldsymbol v}}- c\kappa _{\rm P} E + \eta - \frac{1}{a}\frac{\kappa _{\rm R}}{c}{{\boldsymbol F}}\cdot {\boldsymbol v} \end{eqnarray*} (3) (Reynolds et al. 2009; Bryan et al. 2014). Here, parameter η is the blackbody emissivity given by |$\eta = 4\kappa _{\rm P}\, \sigma _{\rm SB}\, T^4$|⁠, where T is the gas temperature, σSB is the Stefan–Boltzmann constant, and κP is the Planck mean opacity. The term has been added to the above equation, and the radiation pressure tensor, |${\boldsymbol P}$|⁠, has been written using auxiliary functions, \begin{eqnarray*} {\boldsymbol P} &=& {\boldsymbol D}E\nonumber\\ {\boldsymbol D}&=&\frac{1-\chi }{2}\mathbb {I}+\frac{3\chi -1}{2}{\boldsymbol n}\otimes {\boldsymbol n}\nonumber\\ \chi &=&\lambda _{\rm F}+\lambda _{\rm F}^2R^2\nonumber\\ {\boldsymbol n}&=&\frac{\nabla E}{|\nabla E|}. \end{eqnarray*} (4) The equation of the gas energy density e has been modified by introducing the O(⁠|$v$|/c) term, \begin{eqnarray*} &&\frac{\partial e}{\partial t} + \frac{1}{a} \nabla \cdot [(e+p){\boldsymbol v}]\nonumber\\ &&\quad= - \frac{2\dot{a}}{a} e - \frac{1}{a}\rho {\boldsymbol v}\cdot \nabla \phi + c \kappa _{\rm P} E - \eta + \frac{1}{a}\frac{\kappa _{\rm R}}{c}{{\boldsymbol F}}\cdot {\boldsymbol v}. \end{eqnarray*} (5) 2.2 Opacities We use tabulated opacities (Mayer & Duschl 2005) in the form of Planck and Rosseland mean opacities for matter with a primordial composition (Fig. 1). Three elements have been included, namely, H, He, and Li. The following H species have been accounted for: H, H−, H+, H2, H|$_2^+$|⁠, H|$_3^+$|⁠, as well as D, He, and Li. Figure 1. Open in new tabDownload slide (a) Planck mean and (b) Rosseland mean opacities as a function of temperature for different densities, adopted from Mayer & Duschl (2005). These opacities are defined as area per unit volume, and hence have units per unit length. The opacity table in Mayer & Duschl (2005) covers a density range from |$-16 \lt \log \rho \ {\rm [g\, cm^{-3}]} \lt -2$| and a temperature range of |$1.8 \lt \, \log \, T\, {\rm [K]}\, \lt 4.5$|⁠. The table has been extrapolated for temperatures above 2 × 104 K. Figure 1. Open in new tabDownload slide (a) Planck mean and (b) Rosseland mean opacities as a function of temperature for different densities, adopted from Mayer & Duschl (2005). These opacities are defined as area per unit volume, and hence have units per unit length. The opacity table in Mayer & Duschl (2005) covers a density range from |$-16 \lt \log \rho \ {\rm [g\, cm^{-3}]} \lt -2$| and a temperature range of |$1.8 \lt \, \log \, T\, {\rm [K]}\, \lt 4.5$|⁠. The table has been extrapolated for temperatures above 2 × 104 K. The tables cover the ranges |$-16 \lt {\rm log}\, \rho \, {\rm [g\, cm^{-3}]} \lt -2$| in density, and |$1.8 \lt {\rm log}\, T\, {\rm [K]} \lt 4.5$| in temperature (Fig. 1). The temperature-dependence of the opacity has been extrapolated by using the analytic expressions for the free-free, bound-free, and electron-scattering opacities. 2.3 Cooling and heating rates Cooling rates from Luo, Nagamine & Shlosman (2016) have been adopted for the optically thin part of the collapse, solving for internal energy and radiative cooling. A non-equilibrium primordial chemistry network has been invoked to calculate pressure, temperature, ratio of specific heats, and cooling time, for atomic H and He (Abel et al. 1997; Anninos et al. 1997). For this purpose, we used grackle 1.1 package (Bryan et al. 2014; Kim et al. 2014; Smith et al. 2017, https://grackle.readthedocs.org/). We have assumed dust-free gas, and calculated radiative cooling and heating rates accounting for collisional excitation and free-free transitions, recombination, and atomic line excitation. For comparison models with adiabatic equation of state, we have used the local Jeans length, λJ, for each cell, i.e. τ = κPλJ, to calculate the optical depth. An exponential cut off in the cooling rate has been imposed for the optically thin-to-thick transition region, \begin{eqnarray*} \Lambda = \Lambda _{\rm thin}{\rm e}^{-\tau }, \end{eqnarray*} (6) where Λthin is the optically thin cooling rate. LTE has been assumed for the optically thick part of the adiabatic and FLD flows. A number of alternative options have been used to determine the position of the photosphere in the FLD flow (Section 4.1). 2.4 Cosmological initial conditions We use fully cosmological initial conditions (ICs) for our current runs and invoke zoom-in simulations (e.g. Choi et al. 2015; Luo et al. 2016; Shlosman et al. 2016). The gas density exceeds the DM density on spatial scales smaller than ∼0.3–3 pc, and hence the gas decouples from the DM. For the gas, the gravitational softening is adaptive and varies with refinement level. While the cosmological evolution is of course tied to the time since the big bang, we find it helpful to renormalize the time in our simulations to t = 0 when the optical depth in the collapsing gas just exceeds unity. Time before this benchmark is considered negative. We use the WMAP5 cosmology (Komatsu et al. 2009). ΩΛ = 0.721, Ωm = 0.279, Ωb = 0.0445, h = 0.701, σ8 = 0.807, and ns = 0.961. The cosmological ICs are initialized at |$z$| = 199 with the music code (Hahn & Abel 2011), as described in Luo et al. (2016). Generating a set of zoom-in ICs is a two-step process. First, we produce |$1\, h^{-1}$| Mpc comoving 1283 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 |$\mathrel {\raise.3ex\hbox{}\gt {\rm }\lower0.6ex\hbox{}{\sim }{\rm }}10^8\, h^{-1}\, {\rm M_\odot }$| at |$z$| = 10. Secondly, we generate |$0.18\, h^{-1}$| Mpc ICs with 5123 effective resolution in DM and gas embedded in the zoom-in region. Since we use the same random seeds for these ICs as the first step, the phases of both ICs are identical. The zoom-in region is centered on the selected halo position and is set to be large enough to cover the Lagrangian volume of the selected halo and the immediate neighborhood. We have measured the spin parameters of DM haloes within the relevant mass and cosmological spin ranges, λ ∼ 0.01–0.07, in the zoom region at |$z$| ∼ 10. Out of this range, we chose a DM halo with λ ∼ 0.05, which is close to the average spin (e.g. Bullock et al. 2001; Collier, Shlosman & Heller 2018). The gravitational collapse happens at |$z$| ∼ 15.8 in both the adiabatic and FLD runs. At this time, the DM halo has a virial mass of |$1.4\times 10^7\, h^{-1}\,$|M⊙ and virial radius of |${\sim } 0.5\, h^{-1}$| kpc. 3 RESULTS: ADIABATIC FLOW In contrast to the isolated halo model (Luo et al. 2018), cosmological collapse is filamentary and we observe the filaments penetrating inside the growing DM haloes. The adiabatic accretion flow is traced deep inside the parent halo, where the character of the flow changes within ∼10−3 pc, as shown in Fig. 2. We describe evolution of this flow on scales |$\mathrel {\raise.3ex\hbox{}\lt {\rm }\lower0.6ex\hbox{}{\sim }{\rm }}10^{-3}$| pc below using Figs 2 and 3. On these scales, the angular momentum becomes more important, the flow acquires a disc-like character and is fed by a couple of filaments extending from larger scales. The rotational axis of this flow is slightly tilted with respect to the zx-plane. Figure 2. Open in new tabDownload slide Adiabatic collapse: density-weighted map of the gas density projected onto the yz-plane (left-hand column), zx-plane (middle column), and xy-plane (right-hand column), on four spatial scales: 5 kpc (top row), 10−2 pc (second row), 2 × 10−3 pc (third row), and 5 × 10−4 pc (bottom row). All panels shown at the end of the simulation at t = 220.8 yr. Figure 2. Open in new tabDownload slide Adiabatic collapse: density-weighted map of the gas density projected onto the yz-plane (left-hand column), zx-plane (middle column), and xy-plane (right-hand column), on four spatial scales: 5 kpc (top row), 10−2 pc (second row), 2 × 10−3 pc (third row), and 5 × 10−4 pc (bottom row). All panels shown at the end of the simulation at t = 220.8 yr. Figure 3. Open in new tabDownload slide Adiabatic collapse: evolution of the gas density-weighted projection on the zx-plane (face-on view along the y-axis) for the central 2 × 10−3 pc. The black continuous contour represents the position of the ‘photosphere’, at τ = 1 (see Section  2.3). The time t = 0 is defined by the appearance of the photosphere. Figure 3. Open in new tabDownload slide Adiabatic collapse: evolution of the gas density-weighted projection on the zx-plane (face-on view along the y-axis) for the central 2 × 10−3 pc. The black continuous contour represents the position of the ‘photosphere’, at τ = 1 (see Section  2.3). The time t = 0 is defined by the appearance of the photosphere. We display snapshots of disc evolution in Fig. 3, approximately along its rotation axis. The cooling has been artificially and exponentially suppressed at optical depths |$\tau \mathrel {\raise.3ex\hbox{}\gt {\rm }\lower0.6ex\hbox{}{\sim }{\rm }}1$|⁠. At about 250.264 Myr after the big bang, the ‘photosphere’ begins to form, and we re-normalize this time to t = 0. Hence, the photosphere appears within an already discy flow. By t ∼ 20 yr, the disc becomes stratified, both in radial and vertical directions. As the disc becomes visible, it experiences azimuthal distortions, and forms an oval that evolves into a gaseous bar. The amplitude of this bar varies in time, occasionally becoming very strong, and the bar drives a pair of open spiral arms – a sign that angular momentum redistribution continues. The disc size reaches ∼(1–2) × 10−4 pc and ∼(0.5–1) × 10−4 pc, in radius and in height, respectively. So this disc is not geometrically thin. After t ∼ 40 yr, we observe clump formation in the spiral arms, while the central disc and its bar do not experience fragmentation. Luo et al. (2018) argued that Kelvin–Helmholtz (K-H) shear instability (e.g. Chandrasekhar 1961), rather than Toomre/Jeans instability is responsible for similar fragmentation seen in simulations of adiabatic inflow into an isolated DM halo. The main argument is that the clumps form in the spiral arms or shocks, outside the disc. Such a configuration induces shear in the flow, which may therefore be subject to K-H instability, when the Richardson number, Ri < 0.25. This instability will cause the shock front to ‘wiggle,’ and clumps will form and grow (e.g. Balbus 1988; Kim & Ostriker 2002). The gas self-gravity will act as a stabilizer. The gravitational force across the shock front and post-shock layer can be estimated from the value of the non-axisymmetric term in the gravitational potential, perpendicularly to the shock front. To estimate Ri, we assume that the non-axisymmetric force induced by a spiral arm is a fraction β ∼ 0.05 of the radial potential measured by the centrifugal acceleration, |$v_{\rm t}^2/R$|⁠, where |$v$|t is the tangential velocity. This value of β is typical for spiral arms. To project this acceleration perpendicularly to the shock front, we use the pitch angle i of the spiral. Adopting values from the current run that are similar to those in Luo et al. (2018), we obtain Ri ∼ 0.1, confirming that clumps can form as a result of the K-H shear instability. Initially, a single clump appears, but additional clumps continue to form in the spiral arms. Most of the clumps spiral in and merge with the central disc, while two to three outer clumps acquire angular momentum from the bar and move out. These outer clumps grow faster by accretion, especially after t ∼ 60 yr. The gas density profile is shown for a number of representative times (Fig. 4a), and is centered on the densest cell of the most massive object, i.e. the disc. Note the sharp increase in the central density, shortly after formation of the photosphere, to |${\sim } 10^{-6}\, {\rm g\, cm^{-3}}$|⁠. It decreases thereafter, and increases again. By t ∼ 90 yr, the disc becomes heavily distorted by the outer clump, its geometry is complicated but remains planar. Two clumps drive spirals and shocks of their own, and by t ∼ 150 yr, there are basically two cores that grow by accretion. There is in excess of |$300\, \text{M}_\odot$| within the central 10−3 pc at the end of the run, and much of it is found in the diffuse state, the most massive core is |${\sim } 60\, \text{M}_\odot$| (Fig. 4b). Figure 4. Open in new tabDownload slide Adiabatic collapse: evolution of basic parameters of the flow at few representative times: (a) mass density profiles; (b) mass within a fixed spherical radius; (c) temperature profiles; (d) mass accretion rate profiles. Figure 4. Open in new tabDownload slide Adiabatic collapse: evolution of basic parameters of the flow at few representative times: (a) mass density profiles; (b) mass within a fixed spherical radius; (c) temperature profiles; (d) mass accretion rate profiles. The large-scale filaments that feed the central discy flow are relatively cold – their temperature is close to the temperature floor of the atomic gas at ∼3000 K. The disc, outside the photosphere is cold as well. At the same time, the clumps are warmer, T ∼ 104 K, with the most massive clump being the hottest, and so are the spiral arms driven by the asymmetric disc. This higher temperature is associated with larger optical depth for the escaping radiation, as well with increased compression. As the gas shocks and enters the photosphere, its temperature rises sharply by almost a decade. By the end of the run, the central temperature has reached ∼3 × 104 K, when averaged over spherical shells (Fig. 4c). More careful analysis shows that the central temperature fluctuates in the range of T ∼ few × 104–105 K. After the first appearance of the photosphere at few × 10−6 pc, it expands out to ∼2 × 10−4 pc by the end of the run. The mass accretion rate drops inside the photosphere from |${\sim } 1\, \text{M}_\odot \, {\rm yr^{-1}}$| by two to three orders of magnitude, meaning that rotation and gas pressure gradients terminate the inflow below this radius (Fig. 4d), the radiation pressure being much smaller at this time. We calculate the accretion rate, |$\dot{M}$|⁠, by measuring the mass difference within spherical radii per timestep. The importance of angular momentum in the central region of ∼10−3 pc is obvious also from the local velocity field shown in Fig. 5 in three projections. The zx-projection plane shows that the dominant rotation has flattened the object along a plane that is tilted with respect to the zx-plane. It also confirms that the clumps have an intrinsic temperature higher than the accreting gas, including the disc itself. The clumps show a tendency of spiralling in and are expected to merge in a couple of orbital periods. Figure 5. Open in new tabDownload slide Adiabatic collapse: density-weighted map of the temperature projected onto the three principal planes at the end of the simulation at t = 220.8 yr, positioned on the centre of mass. Overplotted arrows represent the velocity field and their length is proportional to the velocity magnitude (see Fig. 6). The colours represent the gas temperature. Figure 5. Open in new tabDownload slide Adiabatic collapse: density-weighted map of the temperature projected onto the three principal planes at the end of the simulation at t = 220.8 yr, positioned on the centre of mass. Overplotted arrows represent the velocity field and their length is proportional to the velocity magnitude (see Fig. 6). The colours represent the gas temperature. The kinematics of the innermost flow has been quantified in Fig. 6, where radial (⁠|$v$|r), tangential (⁠|$v$|t), and circular (⁠|$v$|c) velocity profiles, as well as the sound speed (cs), are shown at four representative times. The maximum value of |$v$|c(r) is increasing with time and moving out until t ∼ 30 yr, then stagnates in position but continues its growth in value – a sign of accumulating mass and the binary nature of the flow before the two clumps merge. The sharp maximum in |$v$|r of the last frame in this figure, at t ∼ 221 yr, reflects the instantaneous radial velocity of the neighbouring clump on an elliptical trajectory. This is also true for the tangential velocity, which is the result of a complex flow driven by two massive clumps, and is nearly cancelled out in the region between them. Note that these clumps have parallel spins. Figure 6. Open in new tabDownload slide Adiabatic collapse: radial profiles of the radial (⁠|$v$|r), tangential (⁠|$v$|t), and circular (⁠|$v$|c) velocities, and sound speed (cs) are shown at four representative times, (a) t = −3.8 yr, (b) t = 1.7 yr, (c) t = 32.6 yr, and (d) t = 220.8 yr. Negative values of radial velocity represent the inflow. These profiles are centered on the densest cell of the most massive object, i.e. the disc. Figure 6. Open in new tabDownload slide Adiabatic collapse: radial profiles of the radial (⁠|$v$|r), tangential (⁠|$v$|t), and circular (⁠|$v$|c) velocities, and sound speed (cs) are shown at four representative times, (a) t = −3.8 yr, (b) t = 1.7 yr, (c) t = 32.6 yr, and (d) t = 220.8 yr. Negative values of radial velocity represent the inflow. These profiles are centered on the densest cell of the most massive object, i.e. the disc. The distribution of the clump masses at the end of the run is given in Fig. 7. There are two massive clumps in the |$20\text{--}90\, \text{M}_\odot$| bin. The distribution peaks between ∼0.1 and |$0.3\, \text{M}_\odot$|⁠. Figure 7. Open in new tabDownload slide Adiabatic collapse: distribution of the clump masses at the end of the simulation, t = 220.8 yr. Figure 7. Open in new tabDownload slide Adiabatic collapse: distribution of the clump masses at the end of the simulation, t = 220.8 yr. To summarize the adiabatic run, the final configuration consists of a rotationally supported, geometrically thick disc. Being asymmetric, it drives a pair of spiral shocks at larger radii. The associated fragmentation is limited to the shocked material and appears to be the result of K-H shear instability and not Toomre gravitational instability. Importantly, the discy character of the inner flow is only perturbed by the clumps but not destroyed. The forming clumps constitute a transient phenomenon. The central mass concentration is of the order of |${\sim } 300\, \text{M}_\odot$| within the central10−3 pc. An important question is whether the adiabatic approximation and the resulting discy flow provide an adequate description for the innermost flow by the FLD approximation, as well. This issue, and the comparison between isolated and cosmological runs, are addressed in Sections 4 and 5. 4 RESULTS: FLD FLOW The FLD run does not differ from the adiabatic run at early time and on large spatial scales. We, therefore, limit our discussion to the central region, |$R\mathrel {\raise.3ex\hbox{}\lt {\rm }\lower0.6ex\hbox{}{\sim }{\rm }}10^{-2}$| pc. The filamentary inflow extends to ∼10−4 pc. The inflow is channelled along the main filament, and the outside material joins the filament after experiencing an oblique shock on its surface. The photosphere forms at ∼250.264 Myr after the big bang, which is very similar to that in the adiabatic run. We define this moment as t = 0, which is used in the subsequent evolution. 4.1 Determining the position of the photosphere For the FLD runs, the position of the photosphere is calculated using a number of alternative options, following Luo et al. (2018). First, we trace rays away from the densest cell in the growing core to a distance of 1 pc, then integrate inwards along the ray to the point where τ = 1, using the Planck mean opacity coefficient, κP. We use 9000 rays equally spaced on a spherical surface. The position of the photosphere is determined for each ray separately. It has no particular symmetry and its shape evolves each timestep. Secondly, the position of the photosphere is calculated using the flux limiter (Section 2.1), as a trace of the optically thick region and shown in Fig. 3. We assume that the photosphere lies at λF = 0.33, where λF decreases sharply with radius. LTE has been assumed for the optically thick part of the flow. The surface τ = 1 is slightly offset outwards from the surface of λF = 0.33. The main problem in calculating the optical depth is that the ray can intersect the fixed τ surface a number of times, the resulting photospheric radius is often overestimated and we avoid using this method here. We plot the surface of constant ionization fraction at the level of 50 per cent in Fig. 8. This surface follows the λF = 0.33 surface sufficiently closely, and lies just outside it. Hence both of these latter surfaces can be used for the purpose of determining the photospheric shape and radius, which lies between them. Figure 8. Open in new tabDownload slide FLD flow: the flux limiter, λF = 0.33, contour (black line) and ionization fraction of 50 per cent contour (red line) superposed on the density slice map (top frame), and ionization slice map (bottom frame) at t = 37.4 yr. The colour palettes are given for each frame. Note the anisotropic ionization map and the generally neutral gas outside the photosphere that lies immediately outside the contours, meaning that the Stromgren sphere size is small and the accreting gas is neutral. Each frame is 2 × 10−4 pc on the side. Figure 8. Open in new tabDownload slide FLD flow: the flux limiter, λF = 0.33, contour (black line) and ionization fraction of 50 per cent contour (red line) superposed on the density slice map (top frame), and ionization slice map (bottom frame) at t = 37.4 yr. The colour palettes are given for each frame. Note the anisotropic ionization map and the generally neutral gas outside the photosphere that lies immediately outside the contours, meaning that the Stromgren sphere size is small and the accreting gas is neutral. Each frame is 2 × 10−4 pc on the side. The evolution of the photosphere has been traced in Fig. 9. Its position has been calculated using the flux limiter. Its shape differs from averaging in spherical shells, and we avoid using the latter, except in Section 5.2. Figure 9. Open in new tabDownload slide FLD collapse: evolution of the central core and its photosphere, superposed on the projected density map. The photosphere is defined using the flux limiter value (Section 2.1) and is given by the thin solid lines. The colour palette provides the mass density of the gas. The size of each frame is 2 × 10−4 pc. Note that the field of view in this figure is a factor of 10 smaller than shown in Fig. 3 for the adiabatic run. Figure 9. Open in new tabDownload slide FLD collapse: evolution of the central core and its photosphere, superposed on the projected density map. The photosphere is defined using the flux limiter value (Section 2.1) and is given by the thin solid lines. The colour palette provides the mass density of the gas. The size of each frame is 2 × 10−4 pc. Note that the field of view in this figure is a factor of 10 smaller than shown in Fig. 3 for the adiabatic run. 4.2 Central core evolution The overall evolution proceeds as follows. The main filament is sheared by rotation on scales of ∼10−5 pc. The mass starts to accumulate where the inflow from the opposite sides of the filament encounters itself. As the mass accumulates in the centre, the filament itself thickens and continues to be sheared in the xz-plane, i.e. the rotation axis of the flow on this spatial scale is nearly aligned with the y-axis, in agreement with the adiabatic run. The central core of ∼10−5 pc becomes visible immediately. The temperatures of the filament and of the core appear to be ∼(1–2) × 104 K, above that of the surrounding gas, which remains near the floor temperature of the atomic gas. The shear, which originates in the rotation of the flow, wraps the filament around the core. The core stays at its original temperature, while the gas around it becomes slightly hotter. A few solar masses of gas are found within the photosphere by t ∼ 5 yr (Fig. 9). In a short time of ∼4 yr after the formation of the photosphere, the first outflow develops and is driven by the gas pressure gradients, with some contribution from the radiation force that becomes dominant occasionally (Fig. 10). The outflow is not symmetric with respect to the core, nor is the temperature of the gas surrounding it. The solid angle of the outflow is about 1 steradian. The resulting hot bubble expands rapidly, driving a dense shell, reaching ∼2 × 10−4 pc in a few years, where it stagnates. The associated velocity of expansion is |${\sim } 50\, {\rm km\, s^{-1}}$|⁠, and hence is supersonic with respect to the ambient gas. Indeed, we observe an expanding shock wave. Figure 10. Open in new tabDownload slide FLD collapse: evolution of the velocity vector field of an inner collapsing flow superposed onto the density slice map at various representative times. The colour palette provides the mass density of the gas. Overplotted arrows represent the velocity field and their size is proportional to the velocity magnitude. Note the alternate phases of inflow and outflow, and their coexistence on various snapshots. The size of each frame is 10−4 pc, twice smaller than in Fig. 9. Figure 10. Open in new tabDownload slide FLD collapse: evolution of the velocity vector field of an inner collapsing flow superposed onto the density slice map at various representative times. The colour palette provides the mass density of the gas. Overplotted arrows represent the velocity field and their size is proportional to the velocity magnitude. Note the alternate phases of inflow and outflow, and their coexistence on various snapshots. The size of each frame is 10−4 pc, twice smaller than in Fig. 9. Even with the outflow present, the inflow continues unabated along the filaments. With time, the outflow region starts to envelop the core, still asymmetry remains strong. This affects the filamentary inflow, which becomes progressively cut off from the core, at least on the side of the outflow. We observe its effect on the mass accretion across the photosphere, (see Section 5.2). At t ∼ 7–9 yr, the mass accretion rate dives by about an order of magnitude, and the photospheric radius shrinks visibly. The mass of the core decreases as well. The temperature of the hot expanding bubble is T ∼ 105 K. As the outflow envelops the central core and evacuates material in its vicinity, we observe that the actual shape of the core resembles that of a triaxial ellipsoid, i.e. it is bar-like, and tumbles. During t ∼ 8–11 yr, the core appears to be cut off completely from the feeding filaments. As the outflow ceases by t ∼ 10 yr, the accretion resumes, as observed also in velocity maps. By this time, much of the core mass is ‘eaten away,’ but it resumes its growth. This pattern of evolution is followed by another round of outflow, which ‘eats up’ the core visibly. By t ∼ 12 yr, the central structure has lost its discy appearance completely, while the core grows. Interestingly, the shell driven by the outflow increases its surface density and forms a mass accumulation by t ∼ 20 yr, which behaves like a fragment. Fig. 9 shows this fragment at t ∼ 22.2 yr, already within the photosphere, falling inwards and dissolving within the central core shortly thereafter. The time periods of recurrent outflows appear to come and go, e.g. at t ∼ 17, ∼22, and especially at ∼24–25 yr (not shown here), then after ∼29 yr (seen only in other projections), and after t ∼ 36 yr. These outflows have a profound effect on the growth of the mass in the central region. The final snapshot of the central core in the FLD evolution is shown in Fig. 9 at t ∼ 37.4 yr after the formation of the photosphere. The gas density radial profile is shown in Fig. 11a. The central density has increased substantially, and reached |${\sim } 10^{-6}\, {\rm g\, cm^{-3}}$| at the end of the run. The density peak at R ∼ 3 × 10−4 pc, visible at t ∼ 37.4 yr, is the mass accumulation at the position of the standing shock of the stagnating bubble. The density profile of the flow departs from the r−2 law within the photosphere bylevelling off. Figure 11. Open in new tabDownload slide FLD collapse: radial profiles of the (a) gas density, (b) enclosed gas mass, (c) gas temperature, and (d) accretion rate at different times. The vertical arrow shows the final position of the photospheric radius. The photospheric temperature fluctuates but, on average, increases monotonically. Figure 11. Open in new tabDownload slide FLD collapse: radial profiles of the (a) gas density, (b) enclosed gas mass, (c) gas temperature, and (d) accretion rate at different times. The vertical arrow shows the final position of the photospheric radius. The photospheric temperature fluctuates but, on average, increases monotonically. Fig. 9 shows that the growth of the central core is not monotonic, and the shape of the photosphere is heavily affected by both inflow and outflow. Within the photosphere, about |$10\, \text{M}_\odot$| have accumulated by the end of the run, and in excess of |$100\, \text{M}_\odot$| are found inside 10−3 pc (Fig. 11b). The photospheric temperature is about 3 × 104 K by the end, up sharply from before the region became optically thick (Fig. 11c). The mass accretion rate fluctuates, as we shall discuss later, and at the end of the run is just below |$1\, \text{M}_\odot \, {\rm yr^{-1}}$| (Fig. 11d). It drops sharply inside the photospheric radius, but also experiences a local minimum around R ∼ 10−2–10−1 pc, due to the increased importance of the angular momentum there. This outer minimum in |${\dot{M}}$| is especially pronounced in isolated models, which have more axisymmetric DM haloes, compared to cosmological models (Choi et al. 2013, 2015). The characteristic velocity profiles of the collapsing flow are given in Fig. 12 at four representative times. Circular velocity provides a measure of the radial mass distribution. We observe that its maximum value does not increase monotonically – another signature of alternating inflows and outflows that affect the central mass accumulation substantially. Radial velocity profiles reflect the same tendency, e.g. at t ∼ 25.7 and 37.4 yr, when the inflow at large radii is reversed at R ∼ 10−4 pc and ∼10−5 pc, respectively, and experiences an outflow. Averaging in spherical shells neglects that both inflow and outflow proceed along particular directions and often coexist at the same radius. The tangential velocity reaches its maximum around the photosphere, but it must be remembered that it is averaged in cylindrical shells and that counterstreams arepresent. Figure 12. Open in new tabDownload slide FLD collapse: radial profiles of the radial, tangential, and circular velocities and the sound speed shown at four representative times. Negative radial velocities correspond to inflow. The arrows in each panel denotes the position of the photospheric radius calculated based on the averaging over spherical shells, except that the tangential velocity averaged on cylindrical shells. The variability arises from the appearance of outflows that disrupt the accretion flow temporarily. Figure 12. Open in new tabDownload slide FLD collapse: radial profiles of the radial, tangential, and circular velocities and the sound speed shown at four representative times. Negative radial velocities correspond to inflow. The arrows in each panel denotes the position of the photospheric radius calculated based on the averaging over spherical shells, except that the tangential velocity averaged on cylindrical shells. The variability arises from the appearance of outflows that disrupt the accretion flow temporarily. Next, we quantify the effect of various forces on the gas, namely, gravity, hydrodynamical force, and radiation force (Fig. 13). During the formation of the photosphere, at t = 0, the radiation force is less important by a factor of a few than the other two forces, up to about a decade in radius outside the photosphere. At larger radii, the radiation force is completely negligible. Gravity dominates outside the photosphere, but inside gravity is balanced by radiation and hydrodynamical forces and some angular momentum. At subsequent times, the radiation force exceeds the hydro force from time to time, especially around the photosphere and at larger radii. The radial profiles of radiation force oscillate widely and correlate with the outflow periods. Towards the end of the run, the radiation force dominates nearly everywhere inside thephotosphere. Figure 13. Open in new tabDownload slide FLD collapse: radial profiles of the different types of acceleration: acceleration due to the gradient of gas thermal pressure, ath, radiative acceleration, ar, and gravitational acceleration, ag, presented at four representative times. The arrow in each panel denotes the photospheric radius. t = 0 represents the appearance of the photosphere in the simulation. Figure 13. Open in new tabDownload slide FLD collapse: radial profiles of the different types of acceleration: acceleration due to the gradient of gas thermal pressure, ath, radiative acceleration, ar, and gravitational acceleration, ag, presented at four representative times. The arrow in each panel denotes the photospheric radius. t = 0 represents the appearance of the photosphere in the simulation. The thermodynamic state of the collapsing gas can be described by averaging properties in spherical shells, but this is not always representative of the actual evolution. Therefore, we show the correlation between the gas temperature and its density using mass averages and displaying all grid cells at the end of the run (Fig. 14). The gas appears nearly isothermal near the cooling floor of the atomic gas, for |$\rho \mathrel {\raise.3ex\hbox{}\lt {\rm }\lower0.6ex\hbox{}{\sim }{\rm }}10^{-12}\, {\rm g\, cm^{-3}}$|⁠. Above this density, the temperature distribution for a given density broadens due to the increasing opacity, and exhibits a upward-pointing ‘hot leg’ distribution, roughly following ρ0.2 (see the inset in Fig. 14), due to the gas becoming optically thick. The important point here is that the temperature is rising slower than in the adiabatic case, where one would expect T ∝ ρ2/3. Figure 14. Open in new tabDownload slide FLD collapse: gas temperature as a function of the gas density, at the end of the simulation, t = 37.4 yr. The colour corresponds to the total mass of all cells for the specific temperature and density. The upward-directed branch reflects the temperature rise due to an increased opacity in the flow. Inset: displays the mass-weighted average profile. The red line shows a T∝ ρ0.2 profile for comparison. Note that adiabatic flow should have T ∝ ρ2/3. Figure 14. Open in new tabDownload slide FLD collapse: gas temperature as a function of the gas density, at the end of the simulation, t = 37.4 yr. The colour corresponds to the total mass of all cells for the specific temperature and density. The upward-directed branch reflects the temperature rise due to an increased opacity in the flow. Inset: displays the mass-weighted average profile. The red line shows a T∝ ρ0.2 profile for comparison. Note that adiabatic flow should have T ∝ ρ2/3. To summarize the results of modelling the FLD flow, the character and evolution of this flow differ fundamentally from the adiabatic flow presented in the previous section and elsewhere in the literature. The central object that forms in the FLD flow is not dominated by rotation. After the initial stage, the radiation force determines its dynamics in tandem with gas pressure gradients. Consequently, no fragmentation occurs. Instead, strong outflows develop but are contained in the central ∼10−3 pc. These outflows slow down the growth of the central mass accumulation, which reaches about |$100\, \text{M}_\odot$| within this radius. Moreover, these outflows mix with the massive accretion flow and transfer angular momentum outwards, lowering the spin of the central object, within its photosphere. This characteristic scale of 10−3 pc is determined by the ability of the accretion flow to contain the outflow, and is about a factor of 100 larger than characterstic scales for star formation. Comparison of the FLD flow with the adiabatic flow is discussed in the next section. 5 DISCUSSION We have followed direct baryonic collapse within DM haloes using high-resolution zoom-in cosmological hydrodynamic simulations. The inclusion of radiative transfer has allowed us to reach the spatial scales of ∼0.01 au, or ∼10−7 pc, for the first time taking into account the associated physical processes involving radiative fluxes and forces in optically thick and partially thick regions. The radiative transfer has been performed in the FLD approximation, and LTE has been assumed for the optically thick collapsing region. The adiabatic model has been evolved for comparison. We find that the collapse proceeds in a filamentary way, and is nearly isothermal in the outer part, down to ∼10−4–10−3 pc from the centre. The gas is channelled along the filaments, with shocks formed by the material joining the filaments. Inside the optically thick region, a central object forms in response to the converging flow and is delineated by its photosphere, initially ∼10−6 pc and expanding thereafter. Reassuringly, it has a similar size as in the isolated collapse (Luo et al. 2018). Moreover, the adiabatic collapse forms a ∼10−6 pc ‘photosphere’ in the cosmological run (Section 3), and ∼2 × 10−6 pc in the isolated run. This core is well resolved throughout the simulations. Its mass |${\sim } 10\, \text{M}_\odot$| within its photosphere is well above the local cell mass of |${\sim } 10^{-6}\, \text{M}_\odot$|⁠, and its central density is about |$10^{-6}\, {\rm g\, cm^{-3}}$|⁠, similar to the adiabatic run. About |$100\, \text{M}_\odot$| have assembled within the central 10−3 pc, and |${\sim } 3\, 000\, \text{M}_\odot$| within 0.1 pc – again, similar to the isolated halo runs. The adiabatic run has about three times larger mass accumulation in a similar time. While the central object is clearly identified in the FLD simulation, it is not expected and indeed is not found to be in hydrostatic or rotational equilibrium. As the flow enters the τ > 1 region, the radial velocity drops, but its internal structure is not relaxed and exhibits streams and turbulent motions. The tangential velocity increases with radius towards the photosphere. The angular momentum profile shows only partial, ∼10 per cent rotational support in this region at the end of the run. The remaining support is provided by the buildup of thermal and radiation pressure gradients. Important new ingredients in the FLD model are the ability of the gas to lose and gain its radiation energy along the radiation energy gradients and addition of the radiation force. This results in slower than adiabatic rise of temperature with density. In the isolated and cosmological models, we find that both the kinematics and the dynamics of the FLD flow differs from the adiabaticcase. A number of important issues must be resolved in order to understand the differences between the adiabatic and FLD runs in particular, and the advanced stage of direct collapse in general. An incomplete list of issues includes the following. Why are the cores obtained in adiabatic and FLD runs so inflated compared to the protostellar stage of massive stars discussed in the literature? Why does angular momentum dominate the central region kinematics in the adiabatic collapse and not in the FLD collapse? Why does fragmentation fail to occur in the central region of the FLD model, contrary to that observed in the adiabatic model? We start by discussing the first question. The forming cores in the adiabatic and FLD runs are not relaxed, neither thermally nor dynamically. Density and temperature profiles as well as all major acceleration profiles, i.e. due to gravity, radiation pressure, and gas pressure are variable. The reason for this is the large accretion rate, which exceeds that encountered in ‘normal’ star formation by orders of magnitude, including formation of the Pop III stars. Moreover, the accretion rate is strongly variable. An important signature of being out of equilibrium is that the maximal central temperature of the gas is |$T\mathrel {\raise.3ex\hbox{}\lt {\rm }\lower0.6ex\hbox{}{\sim }{\rm }}7\times 10^4$| K in the adiabatic run, and |$T\mathrel {\raise.3ex\hbox{}\lt {\rm }\lower0.6ex\hbox{}{\sim }{\rm }}4\times 10^4$| K in the FLD run, when the core reaches a mass of |${\sim } 10\, \text{M}_\odot$|⁠. While this temperature is insufficient to provide for a hydrostatic support due to the gas pressure, this is enough to provide radiation pressure support, because the opacity exceeds that of electron scattering opacity by more than an order of magnitude. In the adiabatic model, this support is provided mostly by the angular momentum, with some contribution from the pressure gradients. 5.1 Outflows and the angular momentum problem Next, we deal with the question of angular momentum redistribution in the adiabatic and FLD runs. The main difference between these runs is the appearance and even dominance of energetic outflows driven by a combination of radiation and gas pressure gradients. These outflows are supersonic and drive shocks into the accretion flow. Luo et al. (2018) observed the complete core dissolution in the FLD model, and this phenomenon has reappeared a number of times in the cosmological runs, as can be seen in Fig. 9 at t ∼ 10.1, 27.9, and 30–33 yr. In these cases, the core did not disappear completely, but lost a substantial mass. Why are the outflows so powerful in the FLD runs? They are so powerful, indeed, that they prevent the core from growing at the full rate provided by the accretion flow. And what is the fate of the expanding gas? These outflows break out in specific directions. Typically, as they evolve, they tend to envelop the photosphere after some time, becoming quasi-isotropic. The outflows are stopped around ∼10−4 pc by the accretion flow, and mix, presumably via a Rayleigh–Taylor instability. Hence, accreting gas accumulates in a shell at ∼10−4–10−3 pc, where it is stirred by the outflows. This phenomenon is important in that it has no known counterpart in star formation, where stellar winds from massive stars disperse the surrounding gas as well as the star-forming cloud itself. In direct collapse considered here, the accretion rate is so high that is capable of containing the outflow from the central core. The ultimate driving force behind these outflows is of course the potential energy of the accretion flow that is converted into kinetic energy and deposited below the photosphere. This process is ∼103–105 times more energetic than during formation of massive stars, e.g. OB type and Pop III stars, because the mass accretion flow is smaller by this factor in the latter cases compared to a direct collapse within DM haloes. Namely, the accretion rate is |${\sim } 10^{-5}\text{--}10^{-3}\, \text{M}_\odot \, {\rm yr^{-1}}$| for massive stars and |${\sim } 1\, \text{M}_\odot \, {\rm yr^{-1}}$| is encountered here. We defer analysis and discussion about the nature of these winds to a later publication. Outflows play an important role in redistributing the angular momentum in the central region, and we elaborate on this below. An important question emerging from comparison of adiabatic and FLD runs is the dominant role of the angular momentum in adiabatic models and its secondary role in the FLD models. Clearly this difference appears only for the innermost flow, roughly within the central ∼10−4 pc. In this region, the two flow realizations differ from each other profoundly. While we do observe some rotation in the FLD flow prior to the formation of the photosphere, it is marginalized shortly thereafter. Here, we attempt to address this important issue. On larger scales, the angular momentum is redistributed by gravitational torques and induced shocks. On scales ∼10−4–10−3 pc of the adiabatic flow, the angular momentum is transferred outwards by the recurrent bar-like perturbation that drives strong spiral shocks (Fig. 3). No such configuration forms in the FLD case, despite the initial sheared flow that is present. As expected, the molecular or ion viscosity has an enormously long timescale and can be neglected. But the forming core drives strong anisotropic outflows repeatedly (e.g. Figs 10 and 15), which extend to ∼10−4 pc at t ∼ 17 yr (see associated snapshot in Fig. 10), and stir the gas within this region, mixing with the accretion flow, and concurrently mixing its angular momentum. Note that the specific angular momentum of the accretion flow is substantially higher than the angular momentum of the outflow. Such a turbulent region drives shocks and is capable of transferring angular momentum outwards. Figure 15. Open in new tabDownload slide FLD collapse: projection snapshots along the three major axes at the end of the simulation, t = 37.4 yr. The velocity field is overplotted on the projected density (top frames), and the temperature map is given in the bottom frames. The arrow size is proportional to the velocity value (see Fig. 12), and the colour palettes are shown on the right margin. The frame size is 10−4 pc. Note the anisotropy of the outflows. Figure 15. Open in new tabDownload slide FLD collapse: projection snapshots along the three major axes at the end of the simulation, t = 37.4 yr. The velocity field is overplotted on the projected density (top frames), and the temperature map is given in the bottom frames. The arrow size is proportional to the velocity value (see Fig. 12), and the colour palettes are shown on the right margin. The frame size is 10−4 pc. Note the anisotropy of the outflows. We, therefore, estimate the timescale invoking a turbulent viscosity in order to extract the angular momentum from the core and vicinity. The timescale for viscosity to have an effect is |$t_{\rm turb}\sim R_{\rm turb}^2/\nu$|⁠, where |$\nu \sim l\, v_{\rm turb}$| is kinematic turbulent viscosity, l is the mean free path for a turbulent cell, Rturb is the size of the turbulent region, and |$v$|turb is the typical turbulent velocity. We take Rturb ∼ l ∼ 10−4 pc, |$v_{\rm turb}\sim c_{\rm s}\sim 2\times 10^6\, {\rm cm\, s^{-1}}$| – all from the FLD runs. The resulting turbulent viscosity timescale is \begin{eqnarray*} t_{\rm turb}\sim 4.8\, \bigg (\frac{l}{10^{-4}\, {\rm pc}}\bigg)\, \bigg (\frac{v_{\rm turb}}{20\, {\rm km\, s^{-1}}}\bigg)^{-1}\, {\rm yr}, \end{eqnarray*} (7) meaning that a turbulent flow can redistribute the angular momentum on a timescale that we observe in the FLD run. This short timescale explains one of the main differences between the adiabatic and FLD flows. 5.2 Luminosity of the central pre-SMBH object One of the differences between the FLD and adiabatic runs is the ability of the former cells to radiate energy in addition to the expansion and contraction available to the cells in the adiabatic approximation. We calculate the radiative luminosity of the pre-SMBH object when it forms a photosphere. This is done by using Fick’s law (equation 2). The emerging photospheric radiation luminosity after t ∼ 5 yr is |$L_{\rm rad}\sim {\rm few}\times 10^{38}-{\rm few}\times 10^{39}\, {\rm erg\, s^{-1}}$| – about the Eddington luminosity for this mass, and exhibits strong variability above and below this range, as shown in Fig. 16c. As the central core grows in mass and in size, its radiation luminosity grows in tandem. Figure 16. Open in new tabDownload slide FLD collapse: time evolution of (a) the mass within the photospheric radius, Rph, enclosed by the mass within 10−5 and 10−4 pc, respectively, for clarity; (b) the photospheric radius, Rph; and (c) the mass accretion luminosity (line with full circles) and radiative luminosity (line with open circles), Lacc and Lrad, respectively. All quantities have been calculated at the photospheric radius that has been obtained averaging in the spherical shells. Note that this method is more approximate than the other methods used in this work, and, therefore, this radius differs somewhat from those used elsewhere in the text. Figure 16. Open in new tabDownload slide FLD collapse: time evolution of (a) the mass within the photospheric radius, Rph, enclosed by the mass within 10−5 and 10−4 pc, respectively, for clarity; (b) the photospheric radius, Rph; and (c) the mass accretion luminosity (line with full circles) and radiative luminosity (line with open circles), Lacc and Lrad, respectively. All quantities have been calculated at the photospheric radius that has been obtained averaging in the spherical shells. Note that this method is more approximate than the other methods used in this work, and, therefore, this radius differs somewhat from those used elsewhere in the text. The energy source of this radiation comes from the potential energy of the accretion flow, Lacc, that is converted into kinetic energy and deposited inside the photosphere, where it is randomized and thermalized. The kinetic luminosity of accretion can be estimated as |${\sim }0.5\dot{\ }M v_{\rm R}^2\sim 10^{38}\, (\dot{M}/1\, \text{M}_\odot \, {\rm yr^{-1}})\, (v_{\rm R}/20\, {\rm km\, s^{-1}})^2\, {\rm erg\, s^{-1}}$|⁠, and varies with time as well. But one should not overlook the associated accretion of the gas thermal energy, |${\sim } 2\pi \rho c_{\rm s}^2\, R_{\rm ph}^2\, v_{\rm R}\sim 3\times 10^{38}\, (R_{\rm ph}/4\times 10^{-5}\, {\rm pc})^2\, (v_{\rm R}/20\, {\rm km\, s^{-1}})\, (c_{\rm s}/15\, {\rm km\, s^{-1}})^2\, {\rm erg\, s^{-1}}$|⁠. Here, all the values are taken at the end of the FLD run, and Rph is the photospheric radius. Hence, unlike accretion on compact objects, the contribution of the second term can play an important role in direct collapse. Obviously, this estimate indicates that the forming core in direct collapse should be relatively loosely bound. In other words, the core material rather ‘levitates’ within the photosphere due to a dominant radiation force – this is supported by Fig. 13d. This situation is expected to be maintained at least until the core grows substantially above the characteristic stellar mass. The multiple periods of outflows, which we have observed during the FLD run, confirm this expectation. Fig. 16 displays the evolution of the photospheric mass, radius and radiation, and accretion luminosities. Strong variability characterizes the evolution of all of these quantities. This variability is about a factor of 100 in amplitude for Lacc and about 1000 for Lrad, superposed on steady growth with an average accretion rate of |${\sim } 0.3\, \text{M}_\odot \, {\rm yr^{-1}}$|⁠. Hence, the rate of supplied accretion energy varies substantially less than the emerging radiationluminosity. The photospheric radius, which shows a slow growth after t ∼ 5 yr (Fig. 16b), exhibits a deep minimum around t ∼ 30–33 yr – a consequence of a strong outflow, which appears to be a response to the peak |$L_{\rm rad}\sim 10^{40}\, {\rm erg\, s^{-1}}$| prior to this time. This alternating growth and decrease in the mass of the core leads to a complicated behavior of Lrad and Lacc in Fig. 16c. There is only a remote correspondence between their oscillations. The reason for this lies in the ability of the mass accretion flow to deposit and store energy deep within the photosphere. The typical diffusion time of the photons from the depth of |$l\sim 0.1\, R_{\rm ph}$| is |$t_{\rm diff}\sim l^2/2D\sim 0.2\, (R/4\times 10^{-6}\, {\rm pc})^2$| yr, where l is the characteristic diffusion radius. This timescale has been Fourier analysed in Luo et al. (2018) for a |$1\, \text{M}_\odot$| core to be at ∼0.12 yr. The same analysis repeated here for Fig. 16c results in an identical timescale of ∼0.10 yr timescale. Accretion luminosity exhibits a characteristic timescale of ∼5 yr, compared to 10 yr in the isolated case. This timescale is related to that of the accretion flow itself and hence can be affected by the outflow feedback. 5.3 No fragmentation in the FLD flow The kinematics of the FLD flow exhibits no dominant role for angular momentum in the core and its vicinity. The initial disc-like flow within the central region does not acquire rotational support because it is capable of transferring its angular momentum via turbulent mixing between the radiation-driven outflow and accretion inflow. This turbulent flow outside the photosphere and very optically thick flow inside the photosphere are not prone tofragmentation. On the other hand, the accretion flow on similar spatial scales in the adiabatic flow forms a disc, partially supported by internal pressure, which also thickens it in the vertical direction. This disc is subject to non-axisymmetric perturbations, mainly the m = 2 mode, and drives spiral shocks (Fig. 3). We have argued here and in Luo et al. (2018) that shocks in the sheared accretion flow will drive the K-H instability and form clumps. Indeed, these clumps form in the spiral shocks and not in the disc itself, which supports our argument. Therefore, we do not agree with the view that the forming clumps are result of gravitational fragmentation in the disc (e.g. Becerra et al. 2015, 2018). In fact, the disc is geometrically thick, which damps the fragmentation exponentially (Toomre1964). 6 CONCLUSIONS We have modelled gravitational collapse of a primordial gas within DM haloes, including radiative transfer following the establishment of a photosphere. This corresponds to the most advanced stage of direct collapse to form seeds of SMBHs at high redshift in a cosmological framework performed so far. Using high-resolution zoom-in cosmological simulations, we have compared runs with an adiabatic equation of state with those in the FLD approximation. We have observed the formation of central cores surrounded by an irregularly shaped photosphere, nearly simultaneously, 250 Myr, after the big bang, i.e. at |$z$| ∼ 15.8, in both approaches. Yet the properties of the cores appear to be quite different. A rotationally dominated core, in the form of a geometrically thick disc, forms in the adiabatic run, supplemented by smaller fragments forming as a result of K-H shear instability in the spiral arms driven by an asymmetric disc. These fragments are transient and eventually merged with the disc, which has a mass of |${\sim } 100\, \text{M}_\odot$|⁠. The central mass concentration achieved at the end of the adiabatic run is about |$300\, \text{M}_\odot$| assembled within the central 10−3 pc, and |${\sim } 3000\, \text{M}_\odot$| within the central 0.1 pc. The central region of the FLD flow evolves differently. The initially discy flow within the central 10−3 pc quickly loses its angular momentum and an amorphous core develops and grows to |${\sim } 10\, \text{M}_\odot$| within a photosphere of close to 10−4 pc. No fragmentation is observed to occur because the central region has lost its angular momentum rapidly and the K-H instability does not operate there. The mass concentration on larger, |$\mathrel {\raise.3ex\hbox{}\gt {\rm }\lower0.6ex\hbox{}{\sim }{\rm }}10^{-3}$| pc scales is similar to that in the adiabatic flow, but its dynamics is fundamentally different. The absence of dominant rotational support of the central object in the FLD run, is due to the development of massive outflows, triggered by the presence of radiation force and gas pressure gradients. These recurrent supersonic outflows are found to drive dense shells of gas by about a factor of 10 in radius, in essence cutting off the core from the accretion flow. They also are responsible for redistributing the angular momentum away from the core. The core radiation luminosity in the FLD run is of the order of the Eddington luminosity, and highly variable, i.e. |$L_{\rm rad}\sim 10^{38}\text{--}10^{39}\, {\rm erg\, s^{-1}}$|⁠. The cores in both runs are substantially inflated in comparison to expected protostellar sizes of comparable masses. The reason for this is the dominant radiation pressure within the photosphere, which results in the gas essentially levitating at the Eddington limit. We confirm our previous result, obtained for direct collapse in isolated haloes (Luo et al. 2018), that radiation transfer allows the gas in the central structure to cool due to anisotropic density and thermal gas and radiation gradients, in the presence of an irregularly shaped photosphere. The resulting rise of temperature with density is substantially shallower than the adiabatic rise. This result is in a strong contrast with the adiabatic approximation for the equation of state used currently in the literature. ACKNOWLEDGEMENTS We thank the enzo and yt support team for help. All analysis has been conducted using yt (Turk et al. 2011), http://yt-project.org/ and VisIt (Childs et al. 2012). We thank Daniel Reynolds for help with the FLD solver, and Kazuyuki Omukai for providing the updated opacities for a comparison. Discussions with Pengfei Chen, Michael Norman, Kazuyuki Omukai, Daniel Reynolds, and Kengo Tomida are gratefully acknowledged. This work has been partially supported by the Hubble Theory grant HST-AR-14584, and by JSPS KAKENHI grant 16H02163 (to IS). IS and KN are grateful for a generous support from the International Joint Research Promotion Program at Osaka University. JHW acknowledges support from NSF grant AST-1614333, Hubble Theory grants HST-AR-13895 and HST-AR-14326, and NASA grant NNX-17AG23G. MB acknowledges NASA ATP grants NNX14AB37G and NNX17AK55G and NSF grant AST-1411879. The STScI is operated by the AURA, Inc., under NASA contract NAS5-26555. Numerical simulations have been performed on XC30 at the Center for Computational Astrophysics, National Astronomical Observatory of Japan, on the KDK computer system at Research Institute for Sustainable Humanosphere, Kyoto University, on VCC at the Cybermedia Center at Osaka University, as well as on the DLX Cluster of the University of Kentucky. REFERENCES Abel T. , Anninos P. , Zhang Y. , Norman M. L. , 1997 , New Astron. , 2 , 181 https://doi.org/10.1016/S1384-1076(97)00010-9 Crossref Search ADS Abel T. , Bryan G. L. , Norman M. L. , 2002 , Science , 295 , 93 https://doi.org/10.1126/science.1063991 Crossref Search ADS PubMed Anninos P. , Zhang Y. , Abel T. , Norman M. L. , 1997 , New Astron. , 2 , 209 https://doi.org/10.1016/S1384-1076(97)00009-2 Crossref Search ADS Balbus S. A , 1988 , ApJ , 324 , 60 https://doi.org/10.1086/165880 Crossref Search ADS Banados E. et al. , 2018 , Nature , 553 , 473 https://doi.org/10.1038/nature25180 Crossref Search ADS PubMed Becerra F. , Greif T. H. , Springel V. , Hernquist L. , 2015 , MNRAS , 446 , 2380 https://doi.org/10.1093/mnras/stu2284 Crossref Search ADS Becerra F. , Marinacci F. , Inayoshi K. , Bromm V. , Hernquist L. , 2018 , ApJ , 857 , 138 https://doi.org/10.3847/1538-4357/aab8f4 Crossref Search ADS Begelman M. C. , Volonteri M. , Rees M. J. , 2006 , MNRAS , 370 , 289 https://doi.org/10.1111/j.1365-2966.2006.10467.x Crossref Search ADS Begelman M. C. , Rossi E. , Armitage P. J. , 2008 , MNRAS , 387 , 1649 https://doi.org/10.1111/j.1365-2966.2008.13344.x Crossref Search ADS Begelman M. C. , Shlosman I. , 2009 , ApJ , 702 , L5 https://doi.org/10.1088/0004-637X/702/1/L5 Crossref Search ADS Begelman M. C. , 2010 , MNRAS , 402 , 673 https://doi.org/10.1111/j.1365-2966.2009.15916.x Crossref Search ADS Berger M. J. , Colella P. , 1989 , J. Comput. Phys. , 82 , 64 https://doi.org/10.1016/0021-9991(89)90035-1 Crossref Search ADS Bromm V. , Loeb A. , 2003 , ApJ , 596 , 34 https://doi.org/10.1086/377529 Crossref Search ADS Bromm V. , Coppi P. S. Larson R. B. , 2002 , ApJ , 564 , 23 https://doi.org/10.1086/323947 Crossref Search ADS Bryan G. L. , Norman M. L. , 1997 , in Clarke D. A. , West M. J. , eds, ASP Conf. Ser. Vol. 123, Computational Astrophysics; 12th Kingston Meeting on Theoretical Astrophysics . Astron. Soc. Pac , San Francisco , p. 363 Google Preview WorldCat COPAC Bryan G. L. , Norman M. L. , Stone J. M. , Cen R. , Ostriker J. P. , 1995 , Comput. Phys. Commun. , 89 , 149 https://doi.org/10.1016/0010-4655(94)00191-4 Crossref Search ADS Bryan G. L. et al. , 2014 , ApJS , 211 , 19 https://doi.org/10.1088/0067-0049/211/2/19 Crossref Search ADS Bullock J. S. , Dekel A. , Kolatt T. S. , Kravtsov A. V. , Klypin A. A. , Porciani C. , Primack J. R. , 2001 , ApJ , 555 , 240 https://doi.org/10.1086/321477 Crossref Search ADS Chandrasekhar S. , 1961 , Hydrodynamic & Hydromagnetic Stability . Dover Press , New York Google Preview WorldCat COPAC Childs H. et al. , 2012 , VisIt: An end-user Tool for Visualizing and Analyzing Very Large Data, No. LBNL-6320E . Ernest Orlando Lawrence Berkeley National Laboratory , Berkeley, CA Google Preview WorldCat COPAC Choi J.-H. , Shlosman I. , Begelman M. C. , 2013 , ApJ , 774 , 149 https://doi.org/10.1088/0004-637X/774/2/149 Crossref Search ADS Choi J.-H. , Shlosman I. , Begelman M. C. , 2015 , MNRAS , 450 , 4411 https://doi.org/10.1093/mnras/stv694 Crossref Search ADS Colella P. , Woodward P. R. , 1984 , J. Comput. Phys. , 54 , 174 https://doi.org/10.1016/0021-9991(84)90143-8 Crossref Search ADS Collier A. , Shlosman I. , Heller C. H. , 2018 , MNRAS , 476 , 1331 https://doi.org/10.1093/mnras/sty270 Crossref Search ADS Eisenstein D. J. , Hut P. , 1998 , ApJ , 498 , 137 https://doi.org/10.1086/305535 Crossref Search ADS Fan X. et al. , 2003 , AJ , 125 , 1649 https://doi.org/10.1086/368246 Crossref Search ADS Federrath C. , Sur S. , Schleicher D. R. G. , Banerjee R. , Klessen R. S. , 2011 , ApJ , 731 , 62 https://doi.org/10.1088/0004-637X/731/1/62 Crossref Search ADS Haehnelt M. G. , Rees M. J. , 1993 , MNRAS , 263 , 168 https://doi.org/10.1093/mnras/263.1.168 Crossref Search ADS Hahn O. , Abel T. , 2011 , MNRAS , 415 , 2101 https://doi.org/10.1111/j.1365-2966.2011.18820.x Crossref Search ADS Hanawa T. , Matsumoto T. , 2000 , PASJ , 52 , 241 https://doi.org/10.1093/pasj/52.2.241 Crossref Search ADS Hirano S. , Hosokawa T. , Yoshida N. , Omukai K. , Yorke H. W. , 2015 , MNRAS , 448 , 568 https://doi.org/10.1093/mnras/stv044 Crossref Search ADS Hirano S. , Hosokawa T. , Yoshida N. , Kuiper R. , 2017 , Science , 357 , 1375 https://doi.org/10.1126/science.aai9119 Crossref Search ADS PubMed Hosokawa T. , Omukai K. , Yoshida N. , Yorke H. W. , 2011 , Science , 334 , 1250 https://doi.org/10.1126/science.1207433 Crossref Search ADS PubMed Hosokawa T. , Hirano S. , Kuiper R. , Yorke H. W. , Omukai K. , Yoshida N. , 2016 , ApJ , 824 , 119 https://doi.org/10.3847/0004-637X/824/2/119 Crossref Search ADS Kim W. , Ostriker E. C. , 2002 , ApJ , 570 , 132 https://doi.org/10.1086/339352 Crossref Search ADS Kim J.-H. et al. , 2014 , ApJS , 210 , 14 https://doi.org/10.1088/0067-0049/210/1/14 Crossref Search ADS Komatsu E. et al. , 2009 , ApJS , 180 , 330 https://doi.org/10.1088/0067-0049/180/2/330 Crossref Search ADS Kravtsov A. V. , Klypin A. A. , Khokhlov A. M. , 1997 , ApJS , 111 , 73 https://doi.org/10.1086/313015 Crossref Search ADS Latif M. A. , Schleicher D. R. G. , Schmidt W. , Niemeyer J. , 2013a , MNRAS , 433 , 1607 https://doi.org/10.1093/mnras/stt834 Crossref Search ADS Latif M. A. , Schleicher D. R. G. , Schmidt W. , Niemeyer J. C. , 2013b , MNRAS , 436 , 2989 https://doi.org/10.1093/mnras/stt1786 Crossref Search ADS Levermore C. D. , 1984 , J. Quant. Spectrosc. Radiat. Tranfer , 31 , 149 https://doi.org/10.1016/0022-4073(84)90112-2 Crossref Search ADS Luo Y. , Nagamine K. , Shlosman I. , 2016 , MNRAS , 459 , 3217 https://doi.org/10.1093/mnras/stw698 Crossref Search ADS Luo Y. , Ardaneh K. , Shlosman I. , Nagamine K. , Wise J. , Begelman M. C. , 2018 , MNRAS , 476 , 3523 https://doi.org/10.1093/mnras/sty362 Crossref Search ADS Mayer M. , Duschl W. J. , 2005 , MNRAS , 358 , 614 https://doi.org/10.1111/j.1365-2966.2005.08826.x Crossref Search ADS Milosavljević M. , Bromm V. , Couch S. M. , Oh S. P. , 2009 , ApJ , 698 , 766 https://doi.org/10.1088/0004-637X/698/1/766 Crossref Search ADS Mortlock D. J. et al. , 2011 , Nature , 474 , 616 https://doi.org/10.1038/nature10159 Crossref Search ADS PubMed Norman M. L. , Bryan G. L. , 1999 , in Miyama S. M. , Tomisaka K. , Hanawa T. , eds, Astrophysics & Space Science Library, Vol. 240, Numerical Astrophysics . Kluwer , Dordrecht , p. 19 Google Preview WorldCat COPAC O’Shea B. W. Norman M. L. , 2007 , ApJ , 654 , 66 https://doi.org/10.1086/509250 Crossref Search ADS Regan J. A. , Haehnelt M. G. , 2009 , MNRAS , 396 , 343 https://doi.org/10.1111/j.1365-2966.2009.14579.x Crossref Search ADS Reynolds D. R. , Hayes J. C. , Paschos P. , Norman M. L. , 2009 , J. Comput. Phys. , 228 , 6833 https://doi.org/10.1016/j.jcp.2009.06.006 Crossref Search ADS Rybicki G. B. , Lightman A. P. , 1979 , Radiative Processes in Astrophysics . Wiley , New York Google Preview WorldCat COPAC Schleicher D. R. G. , Spaans M. , Glover S. C. O. , 2010 , ApJ , 712 , L69 https://doi.org/10.1088/2041-8205/712/1/L69 Crossref Search ADS Shlosman I. , 2013 , Falcon-Barroso J. , Knapen J. H. , eds, Secular Evolution of Galaxies . Cambridge Univ. Press , Cambridge , p. 555 Google Scholar Crossref Search ADS Google Preview WorldCat COPAC Shlosman I. , Choi J.-H. , Begelman M. C. , Nagamine K. , 2016 , MNRAS , 456 , 500 https://doi.org/10.1093/mnras/stv2700 Crossref Search ADS Smith B. D. et al. , 2017 , MNRAS , 466 , 2217 https://doi.org/10.1093/mnras/stw3291 Crossref Search ADS Sur S. , Schleicher D. R. G. , Banerjee R. , Federrath C. , Klessen R. S. , 2010 , ApJ , 721 , L134 https://doi.org/10.1088/2041-8205/721/2/L134 Crossref Search ADS Tanaka T. L. , Li M. , Haiman Z. , 2013 , MNRAS , 435 , 3559 https://doi.org/10.1093/mnras/stt1553 Crossref Search ADS Toomre A. , 1964 , ApJ , 139 , 1217 https://doi.org/10.1086/147861 Crossref Search ADS Truelove J. K. , Klein R. I. , McKee C. F. , Holliman J. H. , II , Howell L. H. , Greenough J. A. , 1997 , ApJ , 489 , L179 https://doi.org/10.1086/310975 Crossref Search ADS Turk M. J. , Abel T. , O’Shea B. , 2009 , Science , 325 , 601 https://doi.org/10.1126/science.1173540 Crossref Search ADS PubMed Turk M. J. , Smith B. D. , Oishi J. S. , Skory S. , Skillman S. W. , Abel T. , Norman M. L. , 2011 , ApJS , 192 , 9 https://doi.org/10.1088/0067-0049/192/1/9 Crossref Search ADS Turk M. J. , Oishi J. S. , Abel T. , Bryan G. L. , 2012 , ApJ , 745 , 154 https://doi.org/10.1088/0004-637X/745/2/154 Crossref Search ADS Venemans B. P. et al. , 2017 , ApJ , 851 , L8 https://doi.org/10.3847/2041-8213/aa943a Crossref Search ADS Willott C. et al. , 2010a , AJ , 139 , 906 https://doi.org/10.1088/0004-6256/139/3/906 Crossref Search ADS Wise J. , Turk A. , Abel T. , 2008 , ApJ , 682 , 745 https://doi.org/10.1086/588209 Crossref Search ADS Wu X.-B. et al. , 2015 , Nature , 518 , 512 https://doi.org/10.1038/nature14241 Crossref Search ADS PubMed © 2018 The Author(s) Published by Oxford University Press on behalf of the Royal Astronomical Society This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/about_us/legal/notices) TI - Direct collapse to supermassive black hole seeds with radiation transfer: cosmological haloes JF - Monthly Notices of the Royal Astronomical Society DO - 10.1093/mnras/sty1657 DA - 2018-09-11 UR - https://www.deepdyve.com/lp/oxford-university-press/direct-collapse-to-supermassive-black-hole-seeds-with-radiation-SDUwYwhZJQ SP - 2277 VL - 479 IS - 2 DP - DeepDyve ER -