Direct formation of supermassive black holes in metal-enriched gas at the heart of high-redshift galaxy mergers
We present novel 3D multi-scale SPH simulations of gas-rich galaxy mergers between the most massive galaxies at , designed to scrutinize the direct collapse formation scenario for massive black hole seeds proposed in Mayer et al. (2010). The simulations achieve a resolution of 0.1 pc, and include both metallicity-dependent optically-thin cooling and a model for thermal balance at high optical depth. We consider different formulations of the SPH hydrodynamical equations, including thermal and metal diffusion. When the two merging galaxy cores collide, gas infall produces a compact, optically thick nuclear disk with densities exceeding g cm. The disk rapidly accretes higher angular momentum gas from its surroundings reaching pc and a mass of in only a few yr. Outside pc it fragments into massive clumps. Instead, supersonic turbulence prevents fragmentation in the inner parsec region, which remains warm ( K) and develops strong non-axisymmetric modes that cause prominent radial gas inflows ( yr), forming an ultra-dense massive disky core. Angular momentum transport by non-axisymmetric modes should continue below our spatial resolution limit, quickly turning the disky core into a supermassive protostar which can collapse directly into a massive black hole of mass via the relativistic radial instability. Such a “cold direct collapse” explains naturally the early emergence of high-z QSOs. Its telltale signature would be a burst of gravitational waves in the frequency range Hz, possibly detectable by the planned eLISA interferometer.
Subject headings:Black hole physics – Galaxies: nuclei – Hydrodynamics – Methods: numerical
The origin of supermassive black holes (SMBHs) is still a puzzle. At least three different classes of models have been proposed: light seeds from Pop III stars, massive seeds from direct gas collapse and light to intermediate mass seeds formed inside star clusters via runaway mergers of stars and stellar remnants. Light seeds from Pop III stars originate at and are the natural consequence of primordial star formation (Madau & Rees, 2001; Volonteri et al., 2003; Tanaka & Haiman, 2009). However, they may face difficulties in growing at a fast enough rate, at least comparable to the Eddington rate, to produce black holes exceeding in less than half a billion years (Johnson & Bromm, 2007; Pelupessy et al., 2007; Wise et al., 2008; Alvarez et al., 2009; Milosavljević et al., 2009), as required by the existence of high- quasars (QSOs; Fan et al. 2001; Mortlock et al. 2011; Treister et al. 2013). Pop III seeds indeed tend to grow at rates well below Eddington as they are surrounded by low density ionized gas and are hosted in small halos that cannot sustain fast inflows of fresh gas from the cosmic web due to their shallow potential wells. Seeds from dense star clusters form when the cluster is subject to rapid mass segregation and consequent runaway mergers that lead to the formation of a central very massive star (VMS; Gürkan et al. 2004; Freitag et al. 2006). The VMS is then expected to end its life leaving behind a BH seed with (see e.g. Portegies Zwart et al., 2004; Belkus et al., 2007). However, intense stellar wind may strongly suppress the growth of the VMS (e.g. Vink, 2008; Glebbeek et al., 2009), unless it evolves in low metallicity environments such as high- nuclear star cluster (Devecchi & Volonteri, 2009; Devecchi et al., 2012). The resulting BH seed might be more massive if the entire star cluster is pushed towards rapid collapse by mass loading from a large scale gas inflow (Davies et al., 2011; Lupi et al., 2014). Finally, direct collapse models rely on gravitational instabilities, such as spiral instabilities and bars-in-bars, to transport gas all the way from galactic disk scales to sub-pc scales (Begelman et al., 2006; Lodato & Natarajan, 2006; Begelman & Shlosman, 2009). As a result, a Jeans unstable cloud forms and it may later evolve into a supermassive star (SMS), which in turn produces a quasi-star (QS) nurturing a massive BH seed ( ) at its center (Begelman et al., 2008; Volonteri & Begelman, 2010; Dotan et al., 2011). Alternatively, if the cloud or supermassive star is massive enough ( ), and has negligible rotational support, it can collapse via the general relativistic radial instability all the way into a black hole encompassing more than half of its original mass (Hoyle et al., 1963; Baumgarte & Shapiro, 1999; Shibata & Shapiro, 2002; Saijo & Hawke, 2009; Montero et al., 2012). However, conventional direct collapse scenarios posit suppression of radiative cooling in order to avoid fragmentation and form a single supermassive object (Lodato & Natarajan, 2006; Dijkstra et al., 2006). This in turn requires primordial gas composition and efficient dissociation of one of the most efficient primordial coolants, namely molecular hydrogen. These requirements can be fulfilled under special environmental conditions such as strong external Lyman-Werner ionizing flux (see Ferrara & Loeb, 2013; Dijkstra et al., 2014). The conditions arising in direct collapse models may also restart the growth of a dormant Pop III seed that has been previously accreting inefficiently (Begelman, 2012; Madau et al., 2014). Interestingly, cosmological simulations have shown convincingly that cold gas accretion flows originating from the cosmic web can sustain a mass accretion rate high enough to produce the high- QSOs on the required short timescale in rare peaks at (halos with virial mass ), provided that BH seeds in such halos start with a mass and that they accrete constantly near the Eddington rate (Di Matteo et al., 2012).
An alternative scenario for direct collapse has been proposed by Mayer et al. (2010, hereafter MA10), and later tested against observational constraints of QSOs using a semi-analytical model attached to the Millennium Simulation by Bonoli et al. (2014). In such a scenario, massive BH seeds originate from the collapse of a supermassive cloud assembled at the end of gas-rich galaxy mergers by multi-scale gas inflows driven by gravoturbulence and gravitational torques. The model does not require primordial gas composition, rather it is designed to be at work in very massive galaxies that have already been enriched to solar metallicity at . MA10 measure gas inflow rates yr at parsec scales and below. The merger-driven model postulates that the formed seeds grow fast but are rare, since they require gas-rich, major mergers between highly biased host halos corresponding to fluctuation (Bonoli et al., 2014). On the small scale side, the model predicts inflows that well exceed those necessary to build a SMS (see Bonoli et al., 2014, and references therein), while other direct collapse models based on instabilities in protogalaxies at fall short of the required inflow rates (e.g. Regan & Haehnelt, 2009a, b; Latif et al., 2012).
The model, however, has been questioned since it is based on simulations with idealized thermodynamical conditions. Indeed, the simulations of MA10 employed an effective equation of state (EoS) with a variable polytropic index across a wide range of gas densities, valid for metal-enriched gas. It has been pointed out that the adopted EoS may artificially promote the formation of a central supermassive gas cloud by suppressing gas fragmentation and star formation. Instead, the expected high cooling rates should lead to widespread fragmentation, disrupting the gas inflow and leading to a BH of only in the most optimistic scenario (Ferrara et al., 2013). The latter arguments, however, are based on simple one-dimensional disk models evolved in the isochoric limit, which forces the density to remain constant at a given radius, neglecting the effect of shocks. These models also neglect self-gravity. Therefore they lack the self-consistent coupling between thermodynamics and gasdynamics. It is conceivable that self-gravity is a crucial ingredient in the MA10 model since the inflow is generated by multi-scale gravitational torques and it is ultimately mediated by a massive, self-gravitating, marginally unstable nuclear gas disk at scales pc.
Here we present 3D hydrodynamic simulations with initial conditions identical to MA10 but with a much richer inventory of physical processes. Most notably, we include radiative cooling in both the optically-thin and optically-thick regimes in place of a prescribed EoS. We even explore the effect of different formulations of SPH, including a novel implementation of the SPH hydro force along with thermal diffusion which successfully obviate to the notorious difficulties of capturing fluid instabilities and mixing (Keller et al. 2014 – see also Ritchie & Thomas 2001 and Hopkins 2013).
As we show in this paper, the basic picture proposed in MA10 is confirmed and further strengthened. Not only the inflow occurs and forms a supermassive cloud even under these more realistic thermodynamical conditions, but it is enhanced by cooling, producing a much more compact and denser cloud on a shorter timescale. As a result, we are left with the intriguing possibility that the collapse proceeds all the way to direct SMBH formation via the relativistic radial instability. Despite the metal-line cooling, fragmentation is limited due to shock heating and gravoturbulence. We conclude that, contrary to previous claims, our scenario for direct collapse is tenable. It offers an explanation for the rapid emergence of high- QSOs without requiring any special conditions in the interstellar medium (ISM), such as primordial gas composition and dissociation of molecular hydrogen.
|Label||SPH force||cooling type||thermal diffusion||metal diffusion|
|RUN3||standard||metal cooling + thermal balance model||yes||yes|
The paper is organized as follows: in Section 2 we describe the setup of the simulations and the implementation of the physical processes included. We present our findings in Section 3, while we discuss possible evolution scenarios of our final configuration towards a massive BH seed in Section 4. In Section 5 we highlight the limitations that could affect our results and we present our conclusions.
2. Numerical Simulations
2.1. Initial conditions and particle splitting
The main purpose of this work is to study the gas dynamics of a gas-rich major galaxy merger from kiloparsec scales down to pc. To this aim, we start by running a large scale galaxy merger simulation and then we apply particle splitting in the later phases of the evolution to reach sub-pc resolution, following MA10. Since we want to compare with our previous work, we adopt the same initial conditions of the reference large-scale merger simulation used in MA10, namely a prograde coplanar 1:1 merger. The merger simulation was initially run with a mass resolution of for the gas particles and a gravitational softening of 100 pc for all galaxy components. The simulation includes radiative cooling, star formation and supernova feedback (see Kazantzidis et al. 2005; Mayer et al. 2007; MA10 for further details). The galaxies approach on a typical parabolic orbit inferred from cosmological simulations (Khochfar & Burkert, 2006). The significant computational burden demanded by these simulations favors the choice of a coplanar prograde merger, which leads to full coalescence of the two galaxy cores on a timescale somewhat shorter than in the other configurations (Kazantzidis et al., 2005; Capelo et al., 2014). Nonetheless, the mass of gas that concentrates at scales pc, which is a key feature in our scenario, has been shown to be rather insensitive to the merger geometry for a fixed mass ratio between the two galaxies (Kazantzidis et al., 2005).
The two galaxy models represent identical multi-component systems with virial mass . They are constructed following the method by Hernquist (1993) and choosing the structural parameters in agreement with the scaling laws predicted by the -CDM model (Mo et al., 1998; Kazantzidis et al., 2005). The models include a dark matter halo that follows the Navarro-Frank-White (NFW; Navarro et al., 1996) profile with concentration , an exponential stellar and gaseous disk and a stellar bulge. The stellar disk has mass , scale radius kpc and scale height pc. The gaseous disk has the same scale lengths of the stellar disk. The mass is about at the time of the last pericenter passage, when star formation has already consumed part of the gas. The bulge is a Hernquist (1990) model with mass and scale length kpc. We refer to MA10 for more details on the numerical parameters of the large scale merger simulation (see also Kazantzidis et al., 2005; Mayer et al., 2007).
We note that the adopted virial mass is comparable to that of a Milky Way-sized galaxy at , and it corresponds to fluctuations at in the WMAP9 cosmology. This is consistent with the very low abundance and high clustering amplitudes of high- QSOs (Volonteri & Rees, 2006; Bonoli et al., 2014). We have chosen model parameters that are similar to those of a present-day, massive, disk-dominated galaxy, which is supported by recent cosmological simulations showing that typical disk dominated galaxies at have structural properties and morphologies already akin to their counterparts at low , with a disk and a sizable bulge component (Fiacconi et al., 2014). However, gas fractions are often higher at high-, in the range , as suggested by both observations and simulations (Tacconi et al., 2010; Moody et al., 2014; Fiacconi et al., 2014), and the most massive among high- galaxies might have hosted clumpy, turbulent disks (Tacconi et al., 2010; Ceverino et al., 2012; Förster Schreiber et al., 2014; Moody et al., 2014) rather than smooth exponential disks as in our initial conditions. Therefore our initial conditions follow a conservative approach since a more gas-rich and gravoturbulent disk will promote loss of angular momentum and central inflow even before the merger occurs (e.g. Bournaud et al., 2012).
As in MA10, we split individual gas particles within a spherical volume of radius 30 kpc into 8 child particles, when galaxy cores are 6 kpc apart at their last apocenter. We reduce correspondingly their gravitational softening, achieving a resolution of and 0.1 pc. Pre-existing star particles and dark matter particles are not split, and their softening is left unvaried. They essentially provide a smooth background potential to avoid spurious two-body heating against the much lighter gas particles. The momentum conserving splitting procedure, which improves over that used in MA10, is described in Roškar et al. (2014). The choices of when to split the gas particles and the size of splitting volume are identical to that in MA10 and are discussed in the Supplementary Information of the latter paper. They are motivated by inducing minimal numerical fluctuations by introducing a refined region large enough to avoid any contamination of low-res particles for the entire duration of the simulations. After the two cores reached a separation pc, most of the gaseous mass is collected in the inner 100 pc volume and it is traced with as many as 1.5 million gas particles. By performing numerical tests we have verified that, owing to the fact that gas dominates the mass and dynamics of the nuclear region, the large softening adopted for the dark matter particles does not affect significantly the density profile of the inner dark halo that surrounds the nuclear disk. These and other numerical tests showing the robustness of the technique were presented in Mayer et al. (2007).
2.2. Features of the simulation suite
We have run four different versions of the second part of the simulation (i.e. after particle splitting) using GASOLINE2, a new version of the GASOLINE code (Wadsley et al., 2004; Keller et al., 2014). Individual runs differ in the sub-grid model for radiative cooling as well as for the specific implementation of the SPH equations of hydrodynamics. We describe below the different features of the runs and we summarize them in Table 1. In the following, we refer to each run with its own label according to Table 1. In all the runs we adopt the metal-dependent, optically-thin cooling introduced in Shen et al. (2010). It considers tabulated cooling rates in ionization equilibrium, while for H and He we compute directly the rates without assuming equilibrium (Wadsley et al., 2004). As in MA10, we assume solar metallicity gas, consistent with observational constraints on the metallicity of the hosts of high- QSOs (Walter et al., 2004). In the following, we call this scheme “metal cooling”. In one of the runs (RUN3), the code switches to an equilibrium temperature-density relation above H cm. This relation is calibrated on 2D radiative transfer calculations based on an improved version of the Spaans & Silk (2000) model, as described in detail in Roškar et al. (2014). In the following, we call this additional feature “thermal balance model”. The “thermal balance model” includes: (i) metallicity-dependent opacity effects due to absorption and scattering of photons by dust; (ii) IR dust radiation; (iii) photoelectric effect on dust; (iv) atomic and molecular line trapping in an ISM irradiated by stellar light; and finally (v) heating by cosmic rays. The thermal balance model thus accounts for self-shielding effects in the dense ISM.
The “thermal balance model” assumes the presence of a uniform UV photon radiation field produced by a starburst in gas of a specified metallicity. Before the final merger, the large scale simulation show that a starburst with a strength of yr takes place in the inner kiloparsec (see MA10). Therefore, we assume such a star formation rate as input to determine the stellar UV flux required by the thermal balance model. While the intensity of the starburst is actually a free parameter, the assumed value is on the low side. A starburst with ten times higher strength, which is likely at high-, would produce a higher UV flux intensity, which in turn would heat up the dust and enhance photoelectric heating of the gas; this would produce a higher gas temperature floor and moderate net radiative losses, strengthening further the main findings of this paper concerning stability to fragmentation (see Section 3).
A pressure floor ensures that the Jeans mass is resolved locally by resolution element, in order to avoid spurious fragmentation. Following Agertz et al. (2009), the minimum pressure is set to , where is a fudge factor, is the gravitational constant, is the density of the particle, and is the local resolution element, i.e. the maximum between the smoothing length and the softening length of the gas particle. All the simulations (including the the large-scale merger simulation before splitting is applied) adopt star formation and blast-wave supernova feedback following Stinson et al. (2006). In particular, we allow stars to form from gas following a Schmidt law in regions above the density threshold H cm and below the temperature threshold K.
We also compare runs with different implementations of SPH that should capture hydrodynamic instabilities and multi-phase fluid mixing with increasing degree of realism. In particular: (i) RUN1 employs the original SPH implementation of GASOLINE (dubbed in the following “standard” SPH; Wadsley et al. 2004); (ii) RUN2 and RUN3 employ the “standard” SPH with the thermal and metal diffusion terms in the momentum and energy equations based on the sub-grid turbulence prescription described by Shen et al. (2010); and (iii) RUN4 employs thermal and metal diffusion as well as the new implementation of the hydrodynamical force equation described by Keller et al. (2014). The new approach is based on the geometric density average (GDSPH) in the SPH force expression in place of the usual , where and are particle pressures and densities respectively. This modification leads to smoother gradients and removes artificial surface tension (Keller et al., 2014). Detailed tests of the GDSPH implementation combined with diffusion as in RUN4 will be presented in Wadsley et al. (in preparation), showing that it can successfully capture Kelvin-Helmoltz and Rayleigh-Taylor instabilities. Note that other new SPH implementations, while they often track entropy rather than energy as we do, also use a geometric density mean for the forces (e.g. Hopkins 2013; Saitoh & Makino 2013).
3.1. Global evolution of the nuclear region
We begin our analysis when the two galaxy cores are only pc apart and kyr away from final coalescence, as shown by the gas density projection in Figure 1. As the two galaxy cores approach their last encounter, efficient cooling leads to some fragmentation within the colliding gas filaments surrounding them. However, most of the gas flows inward unimpeded as the two cores sink owing to dynamical friction against the background matter. Hereafter we define the reference merging time of the cores, , as the time at which two separate density peaks cannot be identified any longer. Our analysis will focus on .
The gas flow is radial, with velocities exceeding 1000 km s, and highly supersonic, as shown by the Mach number map at in the top row of Figure 2 for RUN4. We define the Mach number as the ratio between the module of the gas 3D velocity and the local thermal sound speed . In the inner pc, , leading to strong shocks and the development of supersonic turbulence, as shown by the 3D velocity dispersion maps in the bottom row of Figure 2. The local velocity dispersion is measured as , where is the local 3D gas velocity and the average is intended as SPH average on the smoothing kernel.
After the final coalescence of the two cores, the residual angular momentum in the infalling gas establishes a centrifugal barrier and leads to the formation of a self-gravitating nuclear disk with radius pc, as shown in the gas density projections of Figure 3. At the center, the gas accumulates in a compact, disky core that shows non-axisymmetric features such as modes. This final configuration appears to be largely independent on the specific SPH implementation and cooling module adopted in the various runs (Figure 3) and is maintained till the end of the simulation. We emphasize that the nuclear disk is extremely well resolved owing to particle splitting, with nearly half a million particles.
The central disk is surrounded by filamentary rings of accreting material which exhibit substantial rotation out to several pc during the whole simulation. This is shown by the radial profile of in Figure 4.
quantifies the balance between rotation and random motions by means of the ratio between the azimuthal velocity and the maximum between the radial velocity and the local velocity dispersion . In all runs the central disk undergoes only moderate fragmentation despite being massive and self-gravitating. Indeed, the rotating disk and surrounding infalling filamentary rings are stable except along the densest spiral arms, where gravitationally bound clumps are produced via fragmentation (the stability of the nuclear region is described with detail in Section 3.2).
Figure 5 shows the time evolution of the gas mass inflow, which we interpret as a gas accretion flow onto the central disky core.
The mass inflow rate is measured inside cylindrical shells with vertical thickness of 2 pc as , where is the radial width of the shell, and are the mass and the radial velocity, respectively, of the -th gas particle inside the annulus. Figure 5 shows that accretion rates of yr are sustained from pc, where the radial motions dominate the dynamics of the gas, down to the scale of the inner compact disc111The outflow occurring at the final stage of the evolution in the inner region are due to the presence of a second massive core; see Section 4 for details. pc, where the gas is rotationally supported (see Figure 4).
The interplay between the strong shock heating associated with the highly supersonic inflow and radiative cooling is such that the nuclear disk maintains a temperature K for . This is shown in Figure 6, where we compare the phase diagrams of all the runs at kyr, highlighting the nuclear disk region with black contours.
The phase diagrams are very similar in the high density region occupied by the nuclear disk, regardless of the particular cooling implementation (i.e. “metal cooling” with or without the “thermal balance model”) and the specific formulation of the SPH hydrodynamic force. In the absence of any cooling, an adiabatic, strong shock at the infall velocity found in the simulations ( km s, see the text above and Figure 2) would induce a post-shock temperature K. However, the presence of cooling regulates the temperature in the whole nuclear disk ( pc) to a few thousand K (note that the cooling rate drops significantly below K even for solar metallicity gas). The similarity between the various runs suggests that shock heating during supersonic infall, which is equally present in all of them, plays a major role in balancing cooling and setting the temperature of the core. Likewise, the diagrams show that the “metal cooling” is more effective at lowering the temperature in the low density, low Mach number regions far from the core. The latter appears less effective when GDPSH is employed (RUN4), likely because the new SPH implementation facilitates mixing between different gas phases when used in conjunction with diffusion (Wadsley et al., in preparation). Finally, cooling is enhanced at both low and intermediate densities when the thermal balance model is used (RUN3) because in such a model the cooling rate for solar metallicity gas is increased at low and intermediate densities, before absorption and re-radiation by dust suppresses cooling at high optical depths, as shown in Roskar et al. (2014). We conclude that the warm temperature maintained by the central disk-like core seems to be a robust product of the supersonic gravitational infall established by the larger scale dynamics of a major galaxy merger such as the one simulated here. We also note that no star formation occurs in the warm core simply because the gas never cools enough to meet the temperature condition of our star formation prescription (see Section 2.2). Instead, sporadic star formation can occur in the densest gas pockets further away from the center, where metal-line cooling is effective.
Since the the structure and thermodynamics of the inner region within pc does not change appreciably among the different runs, in the following we will focus on the results of RUN4, which employs both GDSPH and diffusion. This simulation should capture better mixing as well as the multi-phase structure of the flow (Keller et al., 2014).
3.2. A supermassive compact disk in a gravoturbulent, optically thick inflow
Here we discuss in detail the structural properties and evolution of the central nuclear disk. Figure 7 shows the cumulative mass measured inside cylindrical shells with vertical thickness of 3 pc.
Immediately after formation, the compact disk-like core forms and weighs , with variations of less than a factor of 2 between the different runs, and is barely pc in size (see Figure 3). It has has a mean temperature of K and an aspect ratio . It is born with marked non-axisymmetric structure due to its high self-gravity and mass accretion rates yr, which are almost an order of magnitude higher than those found in MA10 with an effective EoS. Then, the inner compact disk grows inside-out as it accretes infalling gas with increasingly higher angular momentum, reaching a mass of of within pc and forming a nuclear disk with an extension of pc after 30 kyr. Such a disk size is a factor of 10 smaller than the radius reached by the disk in MA10 soon after the merger is completed. The disk is much denser than in MA10 because the infalling gas here is much colder than with an effective EoS. The cold infall also enhances the free-fall velocity since the baryonic mass density increases more and faster relative to MA10, an effect analogous to inside-out protostellar collapse (see e.g. Stahler & Palla, 2005). This in turn increases the mass inflow rate since (at least until the effect of angular momentum becomes important). The free-fall time is also proportionally reduced as , where is the size of the collapsing region. Not surprisingly, the inner compact disk is established with a central density exceeding g cm, comparable to the density that in MA10 was reached only after kyr as a result of a secondary infall inside the disk triggered by spiral instabilities. Therefore, owing to radiative cooling the dynamical evolution of the nuclear region is both stronger and accelerated in these new simulations relative to MA10.
The fast accretion rate on to the nuclear disk/inner core continuously maintains a high level of self-gravity (see e.g. Boley 2009 and Hayfield et al. 2011 for an analogous behavior in protostellar disks). This induces the nuclear disk to develop global spiral modes at a few parsec scale and an oval distortion within the central parsec (see Figure 3). In particular, we expect that a further increase in spatial resolution would turn the oval distortion into a strong bar-like mode since we measure within 0.5 pc, where is the rotational kinetic energy and the gravitational binding energy. This value is indeed close to the condition for bar formation in 3D rotationally flattened, self-gravitating axisymmetric configurations (Durisen & Tohline, 1985; Pickett et al., 1996; Christodoulou et al., 1995), where the exact threshold value depends on the equilibrium structure of the configuration (Christodoulou et al., 1995). The value of that we measure is likely reduced by gravitational softening, which damps the growth of non-axisymmetric instabilities, a well documented fact in galactic disk-scale simulations (Mayer & Wadsley, 2004; Debattista et al., 2006; Kaufmann et al., 2007).
We expect that the unresolved bar mode in the core will sustain an inflow down to pc, as found in circumnuclear disk simulations of Choi et al. (2013). This is also in line with the general picture of bars-in-bars proposed by Begelman et al. (2006) (see also Begelman & Shlosman 2009), except for the much higher inflow rates relative to isolated protogalactic disks owing to the trigger from larger scales provided by the merger dynamics. As the system will lower its rotational energy due to angular momentum transport the bar mode mode will weaken, but other instabilities can intervene at lower , for example one-armed modes (Pickett et al., 1996). Considering the conservative case of one-armed modes, less effective than bars at shedding angular momentum, a decrease of a factor of 2 in the specific angular momentum is still expected after a few core rotations (Pickett et al., 1996).
We now examine the conditions of the self-gravitating gas flow in the central disk more closely. Fragmentation in a (laminar) rotating gas disk requires two conditions: (i) the Toomre parameter222The conventional stability threshold is derived using linear perturbation theory for an infinitesimally thin, axisymmetric disk. However, higher order perturbation theory and numerical simulations in 2D and 3D favor the phenomenological threshold in more general conditions (e.g. Durisen et al., 2007). ; and (ii) the local cooling timescale333The cooling timescale condition is well established in simulations of self-gravitating gas disks (e.g. self-gravitating proto-planetary disks; Gammie, 2001; Rice et al., 2003; Mayer et al., 2005; Meru & Bate, 2011a; Helled et al., 2013). However, the exact value of the ratio between the cooling time and the orbital time is still a matter of debate since proof of convergence as a function of resolution and cooling law implementation has not been achieved yet (see e.g. Meru & Bate, 2011b; Helled et al., 2013). Therefore, we will be intentionally conservative and assume that a cooling time equal to the orbital time is a sufficient second condition for fragmentation. has to be shorter than the local orbital time (e.g. Lodato & Natarajan, 2006; Durisen et al., 2007). Figure 8 shows the time evolution of the profiles of our simulation RUN4.
We compute the Toomre parameter taking into account supersonic turbulence in the gas as , where is the epicyclic frequency, is the gas surface density, and . Here, is the 3D velocity dispersion of the gas (measured as in Section 3.1) and the thermal sound speed. The stability threshold is satisfied in most of the nuclear disk after a few thousand year since its emergence, except for some minima below 1.5 typically arising at along spiral arms in the nuclear disk (see Figure 3), which indeed is where gaseous clumps are observed to form via fragmentation. Supersonic turbulence provides support against gravitational collapse, both directly (kinematically) and via heating through small scale shocks as typically observed in star formation simulations (Padoan et al., 1997). Therefore, turbulence influences also the cooling time criterion increasing the effective cooling time. Although we can not probe directly the non-linear coupling between turbulence heating and radiative cooling, we argue that supersonic turbulence plays an important role in regulating the temperature of the nuclear region, since and the runs appear rather insensitive to the details of the cooling within a few parsecs. The key role of turbulence in regulating disk fragmentation was noted also by Choi et al. (2013) and Latif et al. (2013).
We do not include radiative transfer in our simulations, although the “thermal balance model” can be viewed as an approximate way to introduce radiative transfer effects. Nevertheless, in order to support further our results we can show analytically that the intrinsic optically-thick nature of the nuclear disk would suppress fragmentation in the inner parsec. Indeed, the typical optical depth in the disky core within 1 pc is , where cm is the Thomson scattering cross section and H cm is the mean gas column density within pc (corresponding to a surface density pc). Note that varies by an order of magnitude above an below the quoted value in the very inner region and in the low density gaps between spiral arms and rings, respectively. At the temperature of the core, which is K, the opacity due to electron capture by H might indeed be up to ten 10 times higher than for the densities in the inner parsec ( g cm), making our estimate of the optical depth conservative. Finally, the gas would then cool on the photon diffusion timescale yr, where is the speed of light and pc is the vertical scale height of the disk. This timescale is much longer than the orbital time, which is always yr for pc, meaning that no fragmentation should occur. Note for comparison that Ferrara et al. (2013) determined a similar diffusion timescale444They assumed a lower optical depth but a higher value for ., but they compared it with a free-fall time two orders of magnitude longer then the orbital time we find in our simulation, concluding that fragmentation would be ubiquitous in the post-merger core. This is because they underestimated the central density of the disk due to their isochoric assumption. In other words, they explored the effect of radiative cooling without considering how both cooling and shocks in the infalling gas would modify the density structure of the disk in the first place.
The conditions in the disk outside pc change because the average velocity dispersion decreases, the orbital time increases and fragmentation can take place (see upper panel of Figure 3). The co-existence of an inner region stable to fragmentation and an outer region unstable to fragmentation is reminiscent of massive self-gravitating protostellar and protoplanetary disks (e.g. Boley et al., 2010; Helled et al., 2013). Large scale spiral instabilities are present, and massive clumps in the range form along the arms. Most of them appear to be only marginally bound, being dispersed and reforming on a few orbital timescales. The typical temperature of the clumps is K, hence no conventional star formation can take place inside them. Clumps that are dense enough to survive tides will spiral towards the center as a result of dynamical friction on a few orbital times (Noguchi, 1998; Bournaud et al., 2009; Fiacconi et al., 2013).
The gravoturbulent nature of the flow is reflected also by the density profile and the density probability distribution function (PDF) acquired by the gas, shown respectively in the upper and lower panel of Figure 9 at different times.
The mass density profile achieved in the inner 10 pc as a result of the coincident action of self-gravity, turbulence and thermodynamics is close to . Such a density slope has been shown to favor central collapse relative to fragmentation in the envelope for gravoturbulent molecular clouds, a condition that is often invoked to form a massive star rather than a cluster of lower mass stars (Girichidis et al., 2011). This qualitatively corresponds to the behavior observed in our simulations. The density PDF of the gas in the nuclear region carries the imprint of supersonic turbulence in presence of self-gravity. It exhibits a slope close to , in nice agreement with the simulations of Choi et al. (2013), as opposed to the slope of typically found before self-gravitating collapse ensues in simulations of supersonic turbulence in the ISM (e.g. Scalo et al., 1998; Kritsuk et al., 2011). It has been suggested that a slope may result from dynamically important angular momentum (Kritsuk et al., 2011), and indeed we verified that such region corresponds to the region where the nuclear disk assembles, again in agreement with Choi et al. (2013). However, our PDF cannot be simply described by a log-normal distribution function at low densities as in Choi et al. (2013). We argue that the excess amplitude relative to a log-normal PDF at densities g cm is due to the large amount of infalling gas exterior to the inner compact disky-core. Presumably less gas is present at large radii and lower density in Choi et al. (2013) since these authors model an isolated protogalaxy rather than a galaxy merger. Our interpretation is supported by the fact that the low density tail tends to become slightly lower with time as less infalling continues to reach the inner few parsecs.
4. Possible evolutionary pathways of the inner compact core
We stop our simulations between 30 and 50 kyr after because we reached prohibitively small timesteps and the central inflow is saturated at the resolution scale set by the gravitational softening. Although we cannot ascertain the ultimate fate of the gas below the resolution limit, the simulations provide a wealth of information that allows us to delineate possible scenarios to form a massive BH seed given the final conditions inside the pc-size compact disky core present in all of our runs. Hereafter we will refer to the final conditions of RUN4 described in the previous section.
The mass accretion rate has been shown to be the key factor deciding whether a SMS will form or the system will remain in a protostar-like phase, on the Hayashi track, and produce directly a quasi-star (Begelman, 2010; Schleicher et al., 2013; Choi et al., 2013). Schleicher et al. (2013) found that the critical mass necessary for forming a SMS is , where is the mass accretion rate normalized such that for an accretion rate yr. Such critical mass is the mass above which the accretion timescale becomes longer than the Kelvin-Helmoltz timescale in the stellar interior, implying that the protostellar core will contract, the star will leave the Hayashi track and will become a full fledged SMS. The accretion rates found in all our runs in the inner compact disk are typically yr. They correspond to and to a critical mass to form a SMS . This mass is higher than the entire mass reservoir contained within several hundred parsec, suggesting that a SMS would not form. Instead, a QS could form in yr, after the central region has collapsed into a small BH (Schleicher et al., 2013). We caution that we cannot determine the actual radius of the protostar from of our simulations, which should be smaller than our gravitational softening (Hosokawa et al., 2012, 2013), hence we cannot exclude a reduction of the accretion rate at the scale of the protostar. Moreover, we recall that Schleicher et al. (2013) derived their results under the assumption of a spherical, adiabatic, steady contraction without rotation as opposed to a highly dynamical, angular momentum loaded flow as that in our simulations.
Conversely, another condition may become relevant on much shorter timescales if the high inflow rates measured in our simulations can be sustained at scales below our resolution. This is the condition for global relativistic radial collapse in the strong field regime. Numerical general-relativistic (GR) simulations show that the instability for a rotating fluid configuration (a polytrope with ) is reached for: (i) a compactness threshold , and (ii) a dimensionless spin parameter (Baumgarte & Shapiro, 1999; Shibata & Shapiro, 2002; Saijo & Hawke, 2009; Reisswig et al., 2013), where is the mass and the total angular momentum of the cloud. This has been shown to apply to both uniformly rotating and differentially rotating clouds. The inner disky core in our simulations reaches quickly within pc and it has typically . On the other hand, the compactness threshold for is , which is times smaller than the characteristic radius of pc. Therefore, the disky core would seem to be stable to the general relativistic (GR) radial instability. However, as we outlined in section 3.2 we expect that it will become bar-unstable, which would transport angular momentum outward and cause further contraction, possibly pushing it closer to the verge of the instability.
An alternative, perhaps more effective way to discuss stability is to rely on another result of relativistic simulations. Indeed, starting from 2D and 3D rotating polytropic clouds with masses exceeding , these simulations show that or lower is a sufficient condition to bring the cloud to the radial collapse stage under a small initial perturbation 555Recently Montero et al. (2012) have shown that if nuclear burning via the CNO cycle begins at the very center it can induce a shock wave that deflagrates the protostar but so far this has been obtained only in 2D simulations, and for systems with mass . (Shibata & Shapiro, 2002). This is of course a phenomenological criterion, and may change with a different EoS, but offers a useful guideline. If we adopt the latter condition, the angular momentum in the inner compact disk has to decrease substantially to enter radial collapse since at radii pc. Drawing from the calculations of bar-unstable protostellar clouds, which can apply here since the eventual subsequent contraction will mostly be in the newtonian regime, one expects a decrease of the specific angular momentum by a factor of 2 over a few dynamical times, (e.g. Pickett et al., 1996), i.e. over yr. At the same time, the mass of the system can grow up to a factor of , if accretion rates /yr are sustained by the bar down to scales pc for yr. Therefore, would decrease by almost an order of magnitude at fixed radius on relatively short timescales (since it scales as ), reaching the critical threshold for radial collapse.
Although we cannot reach firm quantitative conclusions, our simple estimates suggest the direct massive BH formation via GR radial instability might occur if gas inflow rates as those measured un our simulations are maintained down to scales close to our resolution limit for a time only slightly longer than we could probe here. Since up to of the progenitor mass can be retained during GR radial collapse (e.g. Saijo & Hawke, 2009; Reisswig et al., 2013), the emerging BH seed would have a mass roughly between and , namely in the SMBH mass range and, most importantly, already close the mass inferred for the SMBHs powering the high- QSOs.
Finally, another intriguing aspect of the simulations is that, at least in RUN4, a binary system of two compact disk-like cores eventually arises. This is shown in the gas density projection in Figure 10 at time kyr.
The secondary clump forms after about 25 kyr from fragmentation along a prominent spiral arm. It has a mass only 5 times smaller than that of the primary, which is . Its appearance is also highlighted by the bump at pc in the cumulative mass profile shown in Figure 7. Simulations of massive black hole binaries which evolve in a similar environment suggest that the separations of the two objects will be reduced to pc by dynamical friction against the background in less than 1 Myr (e.g. Mayer, 2013). If a massive BH seed forms in the meantime at the center of both objects one may speculate that a merger between two massive BHs will occur, leading to a burst of gravitational waves (GWs). The frequency of such signal may vary a lot depending on the mass ratio of the BH seeds and on their orbital eccentricity. Assuming conservatively the BH seed mass found in models of direct collapse that posit a quasi-star stage ( ; Dotan et al. 2011), the signal could be detected by the eLISA probe (Amaro-Seoane et al., 2013). Alternatively, if a secondary massive clump is not formed during the stage studied here, the central core might still fragment later during the relativistic phase of the collapse. Indeed in the latter phase a relativistic bar instability could be triggered, which may lead to fission into two or more massive BHs. This was demonstrated by Reisswig et al. (2013) with 3D relativistic hydrodynamics simulations starting from polytropic, marginally unstable, differentially rotating SMS.
We caution that the arguments outlined in this section do not address the thermodynamical evolution during the final collapse stage, which will not necessarily produce a simple polytropic configuration. Hence, only by evolving the end state of our system with a numerical GR code along with a proper treatment of radiation hydrodynamics we will be able to ascertain whether or not direct formation of a SMBH is possible.
5. Discussion and Caveats
We have presented simulations with greater physical realism than those in our previous work (MA10). Our results lend stronger support to the scenario presented in MA10 and Bonoli et al. (2014), in which the SMBHs powering the bright QSOs at form by direct collapse driven by multi-scale gas inflows in major mergers between the most massive galaxies already present at . A key difference with all other direct collapse models presented in the literature is that here we do not require metal-poor gas or dissociation of H to suppress cooling, rather radiative cooling itself fosters the rapid formation of a compact and dense nuclear disk. Fragmentation can happen, but the inner pc-scale core is stabilized by shock heating, gravoturbulence and the high optical depth of the gas. The gas collapses faster and it reaches higher densities relative to MA10 simulations, which adopted and effective EoS rather than incorporating cooling. While the subsequent transport of angular momentum below pc will have to be explored further, we have argued, based on simple analytical arguments, that rapid removal of angular momentum by non-axisymmetric instabilities might continue below the resolution of our simulations. Most importantly, we have highlighted the possibility that the compact disky-core, which encompasses a mass of nearly , may reach the condition for the onset of the GR radial instability by collapsing just an extra decade in radius.
When such a condition is reached, GR simulations tell us that the inner core would leave behind a BH seed exceeding half of its mass. Being really conservative, we can posit that the BH seed will have a mass of . With such a seed, little accretion will be needed to reach rapidly the BH masses inferred for high- QSOs. Indeed, the BH would take yr to grow from to . If we assume an Eddington ratio all the time, the typical value obtained in radiative simulations of gas accretion onto BH seeds (Milosavljević et al., 2009), we obtain yr, comfortably shorter than the time elapsed to in a -CDM cosmology. Note that the Eddington accretion rate would change from to yr between and . The gas reservoir in the inner 100 pc of the merger remnant is and undergoes infall, hence there is plenty of gas to accrete from, a point already emphasized in MA10.
The thermodynamics of the gas will change drastically after a BH seed forms at the center. If the seed emerges already with the mass of a SMBH, any energetic feedback from subsequent accretion, although it can limit its further growth is of little relevance for the origin of the high-z QSOs. Instead, such an early strong feedback mode might have implications for galaxy formation as it could affect gas cooling in nearby halos (see e.g. Dijkstra et al. 2006, 2014) as well as stifle bulge formation in the host galaxy. Feedback can surely have a more important impact on the growth of lighter seeds in the range . However, we have shown that the photon diffusion timescale is longer than the orbital time below parsec scales and thus longer than the free-fall timescale. This means that radiation generated by the accretion process will be trapped in the very inner region on characteristic timescale of the inflow itself. Therefore, the luminosity produced by accretion on to the BH can be integrated over the time interval to yield an estimate of the energy deposited in the inner region, say roughly below 1 pc. Such integrated energy release is , where parametrizes uncertainties on the coupling between gas and radiation, and is the accretion luminosity on to the BH, which we roughly estimate as the Eddington luminosity for a BH. Note is highly uncertain, but we take it as usually done in effective models of BH feedback (see e.g. Di Matteo et al. 2005; Van Wassenhove et al. 2014). Over a timescale of yr we obtain erg. This should be compared to the binding energy of the gas, which is erg assuming pc and . Therefore feedback should not be strong enough to stifle accretion even if the BH seed accretes up to five orders of magnitude above Eddington.
Based on the scenario outlined in the previous section, the massive BH formation process takes less than yr following the merger if GR radial collapse occurs. This implies that bright QSOs with SMBHs can appear essentially as soon as galaxies massive enough to host sufficiently abundant gas reservoirs and multi-scale scale inflows arise. In MA10 and Bonoli et al. (2014) we studied the dependence of the gas inflow on galaxy mass (but fixed gas fraction), finding that a major merger between galaxies having virial mass is a necessary condition for the inflow to produce a central collapse. We can estimate roughly how the inflow rate should scale with virial mass by assuming that the inflow occurs in a self-similar way. Thus, the inflow rate would scale as , since , where is the virial circular velocity of the halo (Mo et al., 1998). In our simulations yr and we expect yr in galaxies ten times lower in mass. This means that the mass accumulated at the center would be an order of magnitude lower on the same timescale and the inner compact disky-core would not be able to reach the onset of the GR radial instability, unless accretion is sustained at least ten times longer ( yr). However, this longer timescale increases the chances that the inflow rate would decrease significantly before reaching the radial instability conditions and opens the path to the formation of an SMS and, eventually, of a QS. In this case the scenario would be identical to that adopted in Bonoli et al. (2014), and the seed that would form has a mass in the range rather than . We calculated that about 200 major mergers between halos above take place at in the Millennium Simulation, whose volume is . If all these events lead to a SMBH-like seed via GR radial instability with mass , the number density of quasars would be . If we further assume that only halos above produce seeds because of the scaling argument for outlined above, then the number of events drops to only about 10, corresponding to a number density Mpc Interestingly, even in this latter case the predicted number density is still high enough to explain all the QSOs, which have an estimated number density of Mpc (Willott et al., 2010; Treister et al., 2013). We note however that this discussion is based on the scaling arguments derived above, which do not necessarily apply to major mergers. Therefore, the dependence on galaxy mass will have to be revisited exploring a wider range of conditions for galaxy mergers, possibly drawn from cosmological simulations.
Some shortcomings are also present in the simulations themselves. A major limitation is the nature of the adopted initial conditions. They are still based on binary mergers rather than on mergers drawn from cosmological simulations. However, a vast body of literature has recently shown that the most massive high- galaxies have turbulent, clumpy disks (e.g. Tacconi et al. 2010; Ceverino et al. 2012; Förster Schreiber et al. 2014; but see also Fiacconi et al. 2014 for lower mass galaxies). These conditions might favor prominent gas inflows during the mergers since non-axisymmetric torquing of the gas will be aided by gas turbulence. This could lead to higher central gas densities in the galaxy cores prior to merging. Clumps formed at galactic disk scales would rain down to the center as the merger comes to completion due to dynamical friction, ultimately aiding the build up of central self-gravitating mass. However, since angular momentum dissipation is crucial, changing the orbit, geometry and properties of the flow by merging misaligned galaxies or galaxies with a more turbulent ISM may change some of the detailed properties of the inflow. This in turn may affect the properties of the compact disk-like core emerging at the center of the remnant. However, it is reassuring that recent studies of the effect of merger geometry on gas inflows down to pc find very weak variations (Capelo et al., 2014).
Another limitation is that the simulations start from merging galaxies that had no star formation happening at scales below 100 pc because of the initial lack of resolution (i.e. before particle splitting is applied). Pre-existing star formation could have stirred the gas in the galaxy cores via stellar feedback prior to the merger, leading to a more fluffy gas distribution at the center of the remnant. We made an attempt to explore this effect by imposing instantaneous feedback in the central 100 pc region after particle splitting is applied. To achieve this we force all star particle located within this volume, each of them sampling a Kroupa IMF, to have an age of about 6 Myr, so that massive stars readily go off as SN Type II immediately nearly everywhere after the beginning of the simulation. This was done in RUN4 and the the resulting effect on gasdynamics was found to be minor. The gas closest to the sites of the explosions is heated and remains hot for the whole duration of the simulation666The blastwave feedback scheme by construction enforces shut-off of radiative cooling for a timescale of Myr (Stinson et al., 2006)., but most of the gas is not directly affected and keeps cooling and flowing inward supersonically.
A potentially more important shortcoming is the lack of feedback provided by the accretion luminosity of the central clump (see Schleicher et al. 2013), which is the analog of protostellar accretion luminosity in star formation (Krumholz et al., 2006). Drawing from the example of massive stars, radiation pressure associated with accretion luminosity may eventually stop the inflow and generate outflows. The largest accretion luminosity is generated in the inner region of the nuclear disk, where the potential well deepens the most. At pc the core generates erg s (using pc, , yr; all these quantities evolve during the simulations but by less than an order of magnitude). This is smaller than the Eddington luminosity, erg s, for an object with mass , suggesting that radiation pressure should not stop the inflow. Nonetheless, during the contraction required in order to reach the stage of the GR radial instability, might increase and become comparable to . At the same, the gas optical depth will increase further and might lead to super-Eddington, radiatively inefficient accretion as mentioned above (Madau et al., 2014). Therefore, the dynamical effect of the accretion radiation has to be probed further in order to understand its impact.
The environment in the nuclear disk is highly optically thick, hence detection of the initial stages of this SMBH formation process in the electromagnetic domain may be difficult. At a temperature of a few thousand K, the inner disk will shine blackbody photons at visible wavelengths into the dense infalling gas, which will re-radiate in the infrared, resembling a protostellar envelope. The outgoing luminosity might be much less than the accretion luminosity estimated above due to absorption by the larger scale gas-rich, dusty nuclear environment. While the James Webb Space Telescope might be a good candidate instrument to detect the infrared signal, extinction and possibly re-emission to even longer wavelengths by the envelope of the nuclear disk might make it more accessible by the Atacama Large Millimeter Array. The detection of a powerful nuclear FIR source associated with highly supersonic gas infall velocities revealed with high resolution spectroscopy could represent the “smoking gun” signature of a direct collapse event. Yet the spectroscopic accuracy required to detect strong velocity gradients at scales of hundreds of parsecs at is hardly attainable with current or upcoming instruments.
Ultimately GWs will provide the cleanest signature of the flavour of direct collapse that we are proposing, namely a “cold direct collapse” into a SMBH driven by the radial instability. Strong GW bursts with a characteristic wave-form and amplitude do indeed arise during the radial collapse phase of (axisymmetric) supermassive objects in numerical GR simulations, having a frequency Hz that falls within the expected eLISA band (Saijo & Hawke, 2009). Furthermore, if binaries of SMBHs arise commonly “at birth” in this scenario, as our simulations suggest, one would detect two relatively low frequency bursts of GWs as the two SMBHs form via the relativistic collapse, followed by an even lower frequency signal as the two holes spiral-in towards coalescence a few Myr later (Mayer, 2013). A detailed study of the expected overall pattern and features of the composite signal, involving considerations of the timescale that separates birth and coalescence of the two SMBHs, with the aim at assessing detectability by the planned eLISA interferometer, warrants investigation in future work.
- Agertz et al. (2009) Agertz, O., Teyssier, R., & Moore, B. 2009, MNRAS, 397, L64
- Alvarez et al. (2009) Alvarez, M. A., Wise, J. H., & Abel, T. 2009, ApJ, 701, L133
- Amaro-Seoane et al. (2013) Amaro-Seoane, P., Aoudia, S., Babak, S., et al. 2013, GW Notes, Vol. 6, p. 4-110, 6, 4
- Baumgarte & Shapiro (1999) Baumgarte, T. W., & Shapiro, S. L. 1999, ApJ, 526, 941
- Begelman et al. (2006) Begelman, M. C., Volonteri, M., & Rees, M. J. 2006, MNRAS, 370, 289
- Begelman et al. (2008) Begelman, M. C., Rossi, E. M., & Armitage, P. J. 2008, MNRAS, 387, 1649
- Begelman & Shlosman (2009) Begelman, M. C., & Shlosman, I. 2009, ApJ, 702, L5
- Begelman (2010) Begelman, M. C. 2010, MNRAS, 402, 673
- Begelman (2012) Begelman, M. C. 2012, ApJ, 749, L3
- Belkus et al. (2007) Belkus, H., Van Bever, J., & Vanbeveren, D. 2007, ApJ, 659, 1576
- Boley (2009) Boley, A. C. 2009, ApJ, 695, L53
- Boley et al. (2010) Boley, A. C., Hayfield, T., Mayer, L., & Durisen, R. H. 2010, Icarus, 207, 509
- Bonoli et al. (2014) Bonoli, S., Mayer, L., & Callegari, S. 2014, MNRAS, 437, 1576
- Bournaud et al. (2009) Bournaud, F., Elmegreen, B. G., & Martig, M. 2009, ApJ, 707, L1
- Bournaud et al. (2012) Bournaud, F., Juneau, S., Le Floc’h, E., et al. 2012, ApJ, 757, 81
- Capelo et al. (2014) Capelo, P. R., Volonteri, M., Dotti, M., et al. 2014, arXiv:1409.0004
- Ceverino et al. (2012) Ceverino, D., Dekel, A., Mandelker, N., et al. 2012, MNRAS, 420, 3490
- Choi et al. (2013) Choi, J.-H., Shlosman, I., & Begelman, M. C. 2013, ApJ, 774, 149
- Christodoulou et al. (1995) Christodoulou, D.M., Shlosman, I, & Tohline, J.E. 1995, ApJ, 443, 551
- Davies et al. (2011) Davies, M. B., Miller, M. C., & Bellovary, J. M. 2011, ApJ, 740, L42
- Devecchi & Volonteri (2009) Devecchi, B., & Volonteri, M. 2009, ApJ, 694, 302
- Devecchi et al. (2012) Devecchi, B., Volonteri, M., Rossi, E. M., Colpi, M., & Portegies Zwart, S. 2012, MNRAS, 421, 1465
- Debattista et al. (2006) Debattista, V. P., Mayer, L., Carollo, C. M., et al. 2006, ApJ, 645, 209
- Di Matteo et al. (2005) di Matteo, T., Springel, V., & Hernquist, L., 2005, Nature, 433, 604
- Di Matteo et al. (2012) di Matteo, Khandai, N, de Graf, C. , Feng, Y., Croft, R.A.C., Lopez, J., & Springel, V., 2012, ApJ, 745, L29
- Dijkstra et al. (2006) Dijkstra, M., Haiman, Z., & Spaans, M. 2006, ApJ, 649, 14
- Dijkstra et al. (2014) Dijkstra, M., Ferrara, A., & Mesinger, A. 2014, MNRAS, 442, 2036
- Dotan et al. (2011) Dotan, C., Rossi, E. M., & Shaviv, N. J. 2011, MNRAS, 417, 3035
- Durisen & Tohline (1985) Durisen, R. H., & Tohline, J. E. 1985, Protostars and Planets II, 534
- Durisen et al. (2007) Durisen, R. H., Boss, A. P., Mayer, L., et al. 2007, Protostars and Planets V, 607
- Fan et al. (2001) Fan, X., Narayanan, V. K., Lupton, R. H., et al. 2001, AJ, 122, 2833
- Ferrara & Loeb (2013) Ferrara, A., & Loeb, A. 2013, MNRAS, 431, 2826
- Ferrara et al. (2013) Ferrara, A., Haardt, F., & Salvaterra, R. 2013, MNRAS, 434, 2600
- Fiacconi et al. (2013) Fiacconi, D., Mayer, L., Roškar, R., & Colpi, M. 2013, ApJ, 777, L14
- Fiacconi et al. (2014) Fiacconi, D., Feldmann, R., & Mayer, L. 2014, arXiv:1410.6818
- Förster Schreiber et al. (2014) Förster Schreiber, N. M., Genzel, R., Newman, S. F., et al. 2014, ApJ, 787, 38
- Freitag et al. (2006) Freitag, M., Gürkan, M. A., & Rasio, F. A. 2006, MNRAS, 368, 141
- Gammie (2001) Gammie, C. F. 2001, ApJ, 553, 174
- Girichidis et al. (2011) Girichidis, P., Federrath, C., Banerjee, R. & Klessen, R.S., 2011, MNRAS, 413, 2741
- Glebbeek et al. (2009) Glebbeek, E., Gaburov, E., de Mink, S. E., Pols, O. R., & Portegies Zwart, S. F. 2009, A&A, 497, 255
- Gürkan et al. (2004) Gürkan, M. A., Freitag, M., & Rasio, F. A. 2004, ApJ, 604, 632
- Hayfield et al. (2011) Hayfield, T., Mayer, L., Wadsley, J., & Boley, A. C. 2011, MNRAS, 417, 1839
- Helled et al. (2013) Helled, R., Bodenheimer, P., Podolak, M., et al. 2013, arXiv:1311.1142
- Hoyle et al. (1963) Hoyle, F. & Fowler, W.A., 1963. MNRAS, 125, 169
- Hernquist (1990) Hernquist, L. 1990, ApJ, 356, 359
- Hernquist (1993) Hernquist, L. 1993, ApJS, 86, 389
- Hopkins (2013) Hopkins, P. F. 2013, MNRAS, 428, 2840
- Hosokawa et al. (2012) Hosokawa, T., Omukai, K., & Yorke, H. W. 2012, ApJ, 756, 93
- Hosokawa et al. (2013) Hosokawa, T., Yorke, H.W., Inayoshi. K., Omukai, K., & Yoshida, N., 2013, ApJ, 778, 178
- Johnson & Bromm (2007) Johnson, J.L., & Bromm, V., 2007, MNRAS, 374, 1557
- Kaufmann et al. (2007) Kaufmann, T., Mayer, L., Wadsley, J., Stadel, J., & Moore, B. 2007, MNRAS, 375, 53
- Kazantzidis et al. (2005) Kazantzidis, S., Mayer, L., Colpi, M., et al. 2005, ApJ, 623, L67
- Keller et al. (2014) Keller, B. W., Wadsley, J., Benincasa, S. M., & Couchman, H. M. P. 2014, MNRAS, 442, 3013
- Khochfar & Burkert (2006) Khochfar, S., & Burkert, A. 2006, A&A, 445, 403
- Kritsuk et al. (2011) Kritsuk, A. G., Norman, M. L., & Wagner, R. 2011, ApJ, 727, L20
- Krumholz et al. (2006) Krumholz, M. R., Matzner, C. D., & McKee, C. F. 2006, ApJ, 653, 361
- Latif et al. (2012) Latif, M. A., Schleicher, D. R. G., & Spaans, M. 2012, A&A, 540, A101
- Latif et al. (2013) Latif, M. A., Schleicher, D. R. G., Schmidt, W., & Niemeyer, J. 2013, MNRAS, 433, 1607
- Lodato & Natarajan (2006) Lodato, G., & Natarajan, P. 2006, MNRAS, 371, 1813
- Lupi et al. (2014) Lupi, A., Colpi, M., Devecchi, B., Galanti, G., & Volonteri, M. 2014, MNRAS, 442, 3616
- Madau & Rees (2001) Madau, P., & Rees, M. J. 2001, ApJ, 551, L27
- Madau et al. (2014) Madau, P., Haardt, F., & Dotti, M. 2014, ApJ, 784, L38
- Mayer & Wadsley (2004) Mayer, L., & Wadsley, J. 2004, MNRAS, 347, 277
- Mayer et al. (2005) Mayer, L., Wadsley, J., Quinn, T., & Stadel, J. 2005, MNRAS, 363, 641
- Mayer et al. (2007) Mayer, L., Kazantzidis, S., Madau. P., Colpi, M., Quinn, T. & Wadsley, J. 2007, Science, 316, 1874
- Mayer et al. (2010) Mayer, L., Kazantzidis, S., Escala, A., & Callegari, S. 2010, Nature, 466, 1082
- Mayer (2013) Mayer, L. 2013, Classical and Quantum Gravity, 30, 244008
- Meru & Bate (2011a) Meru, F., & Bate, M. R. 2011a, MNRAS, 410, 559
- Meru & Bate (2011b) Meru, F., & Bate, M. R. 2011b, MNRAS, 411, L1
- Milosavljević et al. (2009) Milosavljević, M., Bromm, V., Couch, S. M., & Oh, S. P. 2009, ApJ, 698, 766
- Mo et al. (1998) Mo, H. J., Mao, S., & White, S. D. M. 1998, MNRAS, 295, 319
- Moody et al. (2014) Moody, C. E., Guo, Y., Mandelker, N., et al. 2014, MNRAS, 444, 1389
- Montero et al. (2012) Montero, P., Janka, H.T. & Müller, E. 2012, ApJ, 749, 37
- Mortlock et al. (2011) Mortlock, D. J., et al., 2011, Nature, 464, 616
- Navarro et al. (1996) Navarro, J. F., Frenk, C. S., & White, S. D. M. 1996, ApJ, 462, 563
- Noguchi (1998) Noguchi, M. 1998, Nature, 392, 253
- Padoan et al. (1997) Padoan, P., Nordlund, A., & Jones, B.J.T., 1997, MNRAS, 288, 145
- Pelupessy et al. (2007) Pelupessy, F. I., Di Matteo, T., & Ciardi, B. 2007, ApJ, 665, 107
- Pickett et al. (1996) Pickett, B. K., Durisen, R. H., & Davis, G. A. 1996, ApJ, 458, 714
- Portegies Zwart et al. (2004) Portegies Zwart, S. F., Baumgardt, H., Hut, P., Makino, J., & McMillan, S. L. W. 2004, Nature, 428, 724
- Regan & Haehnelt (2009a) Regan, J. A., & Haehnelt, M. G. 2009a, MNRAS, 393, 858
- Regan & Haehnelt (2009b) Regan, J. A., & Haehnelt, M. G. 2009b, MNRAS, 396, 343
- Reisswig et al. (2013) Reisswig, C., Ott, C. D., Abdikamalov, E., et al. 2013, Physical Review Letters, 111, 151101
- Rice et al. (2003) Rice, W. K. M., Armitage, P. J., Bate, M. R., & Bonnell, I. A. 2003, MNRAS, 339, 1025
- Ritchie & Thomas (2001) Ritchie, B. W., & Thomas, P. A. 2001, MNRAS, 323, 743
- Roškar et al. (2014) Roškar, R., Mayer, L., Fiacconi, D., et al. 2014, arXiv:1406.4505
- Saijo & Hawke (2009) Saijo, M., & Hawke, I. 2009, Phys. Rev. D, 80, 064001
- Saitoh & Makino (2013) Saitoh, T.R. & Makino, J. 2013, ApJ, 768, 44
- Scalo et al. (1998) Scalo, J., Vázquez-Semadeni, E., Chappell, D., & Passot, T. 1998, ApJ, 504, 835
- Schleicher et al. (2013) Schleicher, D. R. G., Palla, F., Ferrara, A., Galli, D., & Latif, M. 2013, A&A, 558, A59
- Shen et al. (2010) Shen, S., Wadsley, J., & Stinson, G. 2010, MNRAS, 407, 1581
- Shibata & Shapiro (2002) Shibata, M., & Shapiro, S. L. 2002, ApJ, 572, L39
- Spaans & Silk (2000) Spaans, M., & Silk, J. 2000, ApJ, 538, 115
- Stahler & Palla (2005) Stahler, S. W., & Palla, F. 2005, The Formation of Stars, by Steven W. Stahler, Francesco Palla, pp. 865. ISBN 3-527-40559-3. Wiley-VCH , January 2005.
- Stinson et al. (2006) Stinson, G., Seth, A., Katz, N., et al. 2006, MNRAS, 373, 1074
- Tacconi et al. (2010) Tacconi, L. J., Genzel, R., Neri, R., Cox, P., Cooper, M. C., Shapiro, K., Bolatto, A., Bouché, N., Bournaud, F., Burkert, A., Combes, F., Comerford, J., Davis, M., Förster Schreiber, N. M., Garcia-Burillo, S., Gracia-Carpio, J., Lutz, D., Naab, T., Omont, A., Shapley, A., Sternberg, A., & Weiner, B. 2010, Nature, 463, 781
- Tanaka & Haiman (2009) Tanaka, T., & Haiman, Z. 2009, ApJ, 696, 1798
- Treister et al. (2013) Treister, E., Schawinski, K., Volonteri, M., & Natarajan, P. 2013, ApJ, 778, 130
- Van Wassenhove et al. (2014) Van Wassenhove, S., Capelo, P. R., Volonteri, M., et al. 2014, MNRAS, 439, 474
- Vink (2008) Vink, J. S. 2008, New Astron. Rev., 52, 419
- Volonteri et al. (2003) Volonteri, M., Haardt, F., & Madau, P. 2003, ApJ, 582, 559
- Volonteri & Rees (2006) Volonteri, M., & Rees, M. J. 2006, ApJ, 650, 669
- Volonteri & Begelman (2010) Volonteri, M., & Begelman, M. C. 2010, MNRAS, 409, 1022
- Walter et al. (2004) Walter, F., et al., 2004, ApJ, 615, L17
- Wadsley et al. (2004) Wadsley, J. W., Stadel, J., & Quinn, T. 2004, New A, 9, 137
- Wise et al. (2008) Wise, J. H., Turk, M. J., & Abel, T. 2008, ApJ, 682, 745
- Willott et al. (2010) Willott, C., et al., 2010, AJ, 139, 906