The formation of compact massive self-gravitating discs in metal-free haloes with virial temperatures of \sim 13000-30000 K

The formation of compact massive self-gravitating discs in metal-free haloes with virial temperatures of 13000-30000 K

Abstract

We have used the hydrodynamical AMR code ENZO to investigate the dynamical evolution of the gas at the centre of dark matter haloes with virial velocities of and virial temperatures of at in a cosmological context. The virial temperature of the dark matter haloes is above the threshold where atomic cooling by hydrogen allows the gas to cool and collapse. We neglect cooling by molecular hydrogen and metals, as may be plausible if cooling is suppressed by a meta-galactic Lyman-Werner background or an internal source of Lyman-Werner photons, and metal enrichment has not progressed very far. The gas in the haloes becomes gravitationally unstable and develops turbulent velocities comparable to the virial velocities of the dark matter haloes. Within a few dynamical times it settles into a nearly isothermal density profile over many decades in radius losing most of its angular momentum in the process. About 0.1 - 1 % of the baryons, at the centre of the dark matter haloes, collapse into a self-gravitating, fat, ellipsoidal, centrifugally supported exponential disc with scale-length of pc and rotation velocities of . We are able to follow the settling of the gas into centrifugal support and the dynamical evolution of the compact disc in each dark matter halo for a few dynamical times. The dynamical evolution of the gas at the centre of the haloes is complex. In one of the haloes the gas at the centre fragments into a triple system leading to strong tidal perturbations and eventually to the in-fall of a secondary smaller clump into the most massive primary clump. The formation of centrifugally supported self-gravitating massive discs is likely to be an important intermediary stage en route to the formation of a massive black hole seed.

keywords:
Cosmology: theory – large-scale structure – black holes physics – methods: numerical
Sim Boxsize DM mass R
[Comoving Mpc] [pc]
A 5.0 200.0 15.0
B 2.0 250.0 14.0
C 2.0 250.0 15.0
[kpc] [km ] [K] [K]
1.28 29.9
0.64 19.0
0.59 19.3
Table 1: Basic properties of the three simulations (A, B and C): boxsize(comoving kpc), starting redshift, collapse redshift, DM particle mass (), spatial resolution (pc), total mass of the halo (), virial radius ( kpc), circular velocity (km ), virial temperature (K), maximum baryon density in the halo (), angular momentum parameter , temperature at the core of the halo (K). All units are physical units, unless explicitly stated otherwise.

1 Introduction

Supermassive black holes (SMBH)s were invoked as the central engine powering quasi-stellar objects (QSO)s soon after the first QSOs were discovered Zel’Dovich & Novikov (1964); Salpeter (1964). Early predictions that supermassive black holes are ubiquitous and maybe found in many if not most galaxies as the remnants of dead QSOs Lynden-Bell (1969) have stood the test of time. SMBHs are now believed to reside in most if not all galactic bulges and their masses correlate tightly with the stellar velocity dispersion and the mass of the bulges of their host galaxies (Gebhard et al. 2000; Ferrarese & Merritt, 2000). The luminosity functions of active galactic nuclei across the electromagnetic spectrum and its evolution with redshift have been used to constrain the growth history of supermassive black holes (e.g., Merloni et al., 2004; Shankar et al., 2007; Merloni & Heinz, 2008). Little is known however, about how the first (super-)massive black holes came into being. The discovery of very luminous, bright QSOs at appears to suggest that at least some of the most massive black holes were already in place when the Universe was less than 1 billion years old Fan (2001); Fan et al. (2006); Haiman (2006).
A wide range of generic pathways to a supermassive black hole are possible ranging from direct collapse of gas into rather massive seed black holes, to Eddington-limited accretion onto stellar mass black holes and the dynamical evolution of dense star clusters Rees (1978, 1984). The tight relation between galactic bulges and their central mass suggests that they have formed and grown in a connected way. In the well-established CDM paradigm of structure formation, dark matter haloes grow by hierarchical merging. The last decade has seen extensive and detailed modeling of the build-up of galaxies and supermassive black holes in this context using semi-analytic descriptions as well as numerical simulations Efstathiou & Rees (1988); Carlberg (1990); Haehnelt & Rees (1993); Kauffmann & Haehnelt (2000); Monaco et al. (2000); Cavaliere & Vittorini (2002); Volonteri et al. (2003); Volonteri & Perna (2005); Volonteri & Rees (2005); Volonteri (2007a); Croton (2006); Sijacki et al. (2007); Di Matteo et al. (2008); Hopkins et al. (2008). Due to the lack of observational constraints it is, however, still an open question at what mass scale this hierarchical build-up starts.
The modeling of the growth of supermassive black holes by hierarchical merging of dark matter haloes/galaxies from the remnant black holes of the first generation of stars has identified a number of problems. The shallow potentials hosting the first generation of stars are rather sensitive to the effect of stellar feedback Dekel & Silk (1986); Wise & Abel (2007a); O’Shea & Norman (2008) and may thus not provide sufficient amounts of fuel for a rapid growth of stellar mass black holes. The predicted frequent merging may lead to the expulsion due to gravitational recoil and/or gravitational slingshot unless stellar seed black holes are sufficiently rare (e.g. Volonteri, 2007b). The available time for the formation of the most massive black holes requires almost continuous accretion at the Eddington rate or super-Eddington accretion rates. There is also circumstantial evidence that a certain class of X-ray sources, Ultra-luminous X-ray sources or ULXs for short, are powered by intermediate mass black holes of masses up to Makishima et al. (2000); Fabbiano et al. (2003); Colbert et al. (2004); Fabbiano et al. (2006). There is thus plenty of motivation to consider seriously the possibility that the growth of supermassive black holes has started from massive seed black holes with masses substantially larger than those forming as stellar remnants (see Begelman et al. (2006) for a recent discussion).
The centres of high-redshift dark matter haloes with virial temperatures K, not yet significantly enriched with metals, have been identified as a promising environment to form such massive seed black holes. If cooling is suppressed in these haloes either due to an external UV background or more likely especially in the later stages of the collapse due to internal sources of UV radiation (cf. Wise & Abel, 2007b) early fragmentation should not occur favouring the formation of a compact massive rotationally supported disc Oh & Haiman (2002); Bromm & Loeb (2003); Volonteri et al. (2003); Koushiappas et al. (2004); Begelman et al. (2006); Lodato & Natarajan (2006); Rees & Volonteri (2007); Volonteri et al. (2008). This is rather different from the situation in the lower mass haloes studied extensively as possible sites for the formation of the “first” stars where cooling is the dominant cooling mechanism (e.g. Bromm et al., 1999, 2002; Bromm & Larson, 2004; Abel et al., 1998, 2000, 2002; O’Shea & Norman, 2007, 2008).
We use the publicly available code Enzo1 to perform adaptive mesh simulations of the dynamical evolution of the gas in three dark matter haloes with virial velocities between and and virial temperatures between and at in a cosmological context. The work is similar in spirit to that of Wise et al. (2007), hereafter W07, who simulated two haloes somewhat less massive than our haloes. W07 find that the gas at the centre of their haloes does not reach rotational support and does not fragment. Our simulations have a similar set-up but unlike W07 we do not push for resolution in the central region where a very small fraction of the gas can collapse to very high density but rather follow the dynamical evolution of a significant fraction of the gas settling into a rotationally supported disc for several dynamical times. The paper is structured as follows. In §2 we describe the details of the numerical simulations. In §3 we summarise some basic formulae that we use to characterise the properties of the haloes. In §4 we describe the results of our numerical simulations and in §5 we summarise our conclusions. Throughout this paper we assume a standard CDM cosmology with the following parameters (Spergel et al., 2003, based on the WMAP 1st year data), = 0.74, = 0.26, = 0.0463, = 0.9 and = 0.72. We further assume a spectral index of primordial density fluctuations of .

2 The Setup of the numerical simulations

2.1 The adaptive-mesh refinement code Enzo 

We have used the publicly available adaptive mesh refinement (AMR) code Enzo. Enzo was originally developed by Greg Bryan and Mike Norman at the University of Illinois (Bryan & Norman 1995b, Bryan & Norman 1997, Norman & Bryan 1999, O’Shea et al. 2004). The gravity solver in Enzo uses an N-Body particle mesh technique (Efstathiou et al. 1985, Hockney & Eastwood 1988). Additional finer meshes can then be added in regions of high density to calculate the dynamics of the dark matter (DM) particles more accurately.
The hydrodynamics solver employs the Piecewise Parabolic method combined with a non-linear Riemann solver for shock capturing. The Eulerian AMR scheme was first developed by Berger & Oliger (1984) and later refined by Berger & Colella (1989) to solve the hydrodynamical equations for an ideal gas. Bryan & Norman (1997) adopted such a scheme for cosmological simulations. In addition to this there are also modules available which compute the radiative cooling of the gas together with a multispecies chemical reaction network. There are two versions of the chemistry solver available, one with 6 species () and one with 9 species (same as before plus ). As stated previously the simulations conducted here do not include cooling and chemistry. Our simulations make extensive use of Enzo’s capability to employ nested grids which we describe in the next section.

2.2 Nested grids & initial conditions

Initial conditions were generated with the initial conditions generator supplied with the Enzo code. The nested grids are introduced at the initial conditions stage. We have first run exploratory DM only simulations with coarse resolution, setting the maximum refinement level to 4. These DM only simulations have a root grid size of and no nested grids. In these exploratory simulations we have identified the most massive halo at a redshift of 10 and then rerun the simulations, including the hydrodynamics module. We also introduce nested grids at this point. The nested grids are placed around the region of interest, as identified from the coarse DM simulation. We have used four levels of nested grids in our simulations with a maximum effective resolution of . The introduction of nested grids is accompanied by a corresponding increase in the DM resolution by increasing the number of particles in the region of interest. We only use the highest resolution nested grid to refine further thus economising our computational output .This corresponds to a region of (comoving) kpc for simulation A and (comoving) kpc for simulations B & C. The total number of particles in our simulation is 4935680, with of these in our highest resolution region. The nested grids are distributed as follows in root grid sizes, . Table 1 gives the details of the simulations discussed here.

Figure 1: A sequence of visualisations of the density distribution in simulation A. The gas is collapsing in a halo of DM mass of , virial velocity of and virial temperature of . At the end of the simulation the gas at the centre of the halo has settled into a rotationally supported disc. Each plot shows a thin slice of the density distribution centered on the grid cell with the highest density in our halo of interest. Between panels either the density range or the time of the simulation output change. The colour range represents approximately an order of magnitude range in density in each plot. The scales are proper distances.
Figure 2: Same as figure 1 for simulation B. The gas is collapsing in a halo of DM mass of , virial velocity of and virial temperature of . Note how the gas fragments into three clumps at Myrs which tidally distort each other and undergo a violent dynamical interaction.
Figure 3: Same as figure 1 for simulation C but only for the later stages of the evolution. The gas is collapsing in a halo of DM mass of , virial velocity of and virial temperature of . Panels 9 to 15 illustrate the dynamical evolution of the rotationally supported gas at the centre of the halo for a few dynamical times.
Figure 4: Left-hand Panels: The temperature profile for the halo in simulations A and B as a function of radius. Middle Panels: The time evolution of the radial velocity (black), thermal velocity (blue) and the turbulent velocity (red) as a function of total enclosed gas mass. Right-hand Panels: Ratio of the gravitational acceleration to the acceleration due to thermal pressure as a function of total enclosed mass. Note that initially the central regions are supported by thermal pressure.
Figure 5: Left-hand Panels: The evolution of the enclosed (gas) mass within spherical shells around the centre of the halo as a function of radius. The enclosed DM mass profile at the end of the simulation is shown by the filled diamonds. Right-hand Panels: The evolution of the (gas) density profile during the collapse. The gas settles into an close to isothermal () density profile (thick dot-dashed line) over many decades in radius. The DM density at the end of the simulation is shown by the filled diamonds.

Figure 6: : Left-hand Panel: The evolution of the total specific angular momentum () of the gas as a function of enclosed mass for eight of the simulation outputs in figure 1. Right-hand Panel: The same for simulation B for eight simulation outputs shown in figure 2.
Figure 7: Left-hand Panels: The circular velocity (blue), and our two estimates of the rotation velocity (black) and (red) are plotted as a function of radius. As the collapse proceeds the inner part of the halo settles into centrifugal support. Right-hand Panels: The ratio of to , for the same six output times plotted against radius. The central region attains rotational support, after Myrs and Myrs, in simulations A and B, respectively.

2.3 Following the collapse of the gas and dark matter with Enzo 

As mentioned previously Enzo uses adaptive grids to provide increased resolution where it is required. For the simulations discussed in this paper we have used four refinement criteria implemented in enzo: DM over-density, baryon over-density, Jeans length and cooling time. The first two criteria introduce additional meshes when the over-density of a grid cell with respect to the mean density exceeds 3.0 for baryons and/or DM. The third criterion is the Truelove criterion Truelove et al. (1997) which in its original form states that at least 4 grid cells should be used to resolve the Jeans length to ensure that no artificial fragmentation takes place. The cooling time is often shorter than the dynamical time and the cooling time refinement criterion helps to ensure that the Jeans mass is properly resolved. As other authors (e.g. O’Shea & Norman, 2008, W07) we are conservative here and set this criterion to grid cells in our simulations. The fourth criterion ensures that the cooling time in a given grid cell is always longer than the sound crossing time of that cell. We also set the MinimumMassForRefinementExponent parameter to making the simulation super-Lagrangian and therefore reducing the criteria for refinement as higher densities are reached O’Shea & Norman (2008). We furthermore set the MinimumPressureSupportParameter equal to ten as we have restricted the maximum refinement level in our simulations (e.g. Kuhlen & Madau, 2005). With the MinimumPressureSupport option the code assumes in each computational cell the minimum temperature necessary to make the cell Jeans stable at the highest refinement level.
We first advance our simulation forward in time using a maximum refinement level of 4. If we would have run the simulation with a high maximum refinement level from the outset the simulation would have stalled once the virial temperature of the halo has reached 10000K in the attempt to follow the dynamical evolution of a small dense region before a sufficiently massive DM halo had formed. We are effectively delaying the collapse by limiting the resolution. We use a friends of friends (FoF) algorithm to keep track of how our halo is progressing through time. Once a sufficiently massive halo (corresponding to a DM particle number of about 30000 particles in simulation A and 120000 particles in simulations B and C) has formed we halt the simulation. We then restart the simulation with the maximum level of refinement set to 18 (16 in simulation B and C). We will refer to this point as the collapse redshift.
The simulation then continues at a much higher resolution albeit at a significantly slower speed. We allow the simulation to evolve until the code no longer tracks the hydrodynamic evolution of the gas accurately due to a lack of resolution at the highest densities caused by restricting the highest allowed refinement level. The Enzo code issues warnings alerting the user when the code begins to produce spurious results due to a lack of resolution. At this point we terminate the simulation.
We have also followed the recommendation in Lukic et al. (2007) for the starting redshift of our simulations. As a result the initial redshift is quite a bit higher than that of other simulations described in the literature. Experimenting with the value of the initial redshift suggests, however, that this increased initial redshift has little effect on our results.

3 Characteristic Formulae

As mentioned before we use a FoF algorithm to identify DM haloes in our simulation outputs. We adopt a linking length of 0.2 for all of our analysis (see Jenkins et al., 2001; Lukic et al., 2007). We compute physical properties of the virialised haloes using standard formulae (see e.g. Mo & White, 2002). For convenience we summarise some of the formulae below. We calculate characteristic properties for a sphere for which the mean enclosed density is 200 times the mean cosmic value . The characteristic “virial” radius is then related to the mass as,

(1)

while the circular velocity, , can be written as

(2)

where is Hubble’s constant at redshift and is the density parameter of non-relativistic matter, is the mass of our halo as determined by FoF group finder. Rewriting the equation for the virial velocity in terms of the mass, , only gives

(3)

The Hubble constant and density parameter are related to their present-day values by

(4)

and

(5)

where E(z) is given by

(6)

We can now define the virial temperature of the halo as

(7)

where , with being the proton mass, is the mean molecular weight.
The level of rotational support is often characterized by the dimensionless angular momentum parameter which can be written as (e.g. Bullock et al., 2001)

(8)

where is the angular momentum inside a sphere of radius containing mass and is the virial velocity of the halo. The mean value of is for typical haloes in cosmological simulations Barnes & Efstathiou (1987); Bett et al. (2007). We have included the above properties for the three virialised haloes in table 1 for reference.
Also of interest is the dynamical time which we express in terms of the enclosed mass at radius ,

(9)

4 Results of the numerical simulation

4.1 The global dynamics of the collapse of the gas and dark matter

We first discuss the global dynamical evolution of the gas and dark matter in the three haloes we have simulated. Note that the dynamical evolution in simulation C is similar to that in simulation A and we will thus show only a subset of the plots for simulation C in the following. We have simulated boxes with Mpc and Mpc on a side, respectively. Figure 1 and 2 show the location of the haloes in simulation A and B within the typical characteristic web of filaments and sheets which is already well developed at (simulation C looks very similar in this respect). The haloes are located at the intersection of several filaments. In the top row of figure 1 and 2 we zoom in several steps onto the haloes which are approximately spherical and have a virial radius of about 1.5 and 0.5 kpc. In Figure 3 we focus on the later stages of the evolution of the gas at the centre of simulation C. From one panel to the next we either change the output time or the density/length scale. We have selected time 0 as the time at which we re-start the simulation with increased resolution, with a maximum refinement level of 18 for simulation A and a maximum refinement level of 16 for simulations B & C. Recall that up to now we had run the simulation limiting the refinement to 4 levels in order to allow the halo to build up. We will discuss possible limitations of this approach later. Once the resolution is increased the inner part of the haloes starts to collapse.
As discussed by W07 the collapse is highly turbulent and dynamically complex. In each halo we are able to follow the collapse of the gravitationally unstable gas at the centre of the halo into a centrifugally supported disc (or disc-like object) and follow the dynamical evolution of the disc for several dynamical times. In simulation A this disc has a mass of . In the second of our haloes the gas at the centre of the halo fragments into three gravitationally bound clumps with masses of a few times and a mass ratio of approximately 3:1:1. The clumps tidally distort each other and subsequently undergo a violent dynamical interaction. One of the smaller clumps eventually merges with the most massive clump and forms a disc of . The second small clump is tidally stripped of most of its mass and has a mass of a few times at the end of the simulation. The (temporary) formation of multiple systems appears, however, not to be generic. Neither simulation A with its five times more massive halo nor simulation C with its halo of similar mass show the formation of a similar multiple system. The disc forming in simulation C has a mass of .

4.2 The onset of gravitational instability, mass inflow and the evolution of the density profile

The gas at the centre of our haloes attains an almost isothermal temperature profile with a characteristic temperature of 6000-7000K the temperature below which atomic cooling due to hydrogen cooling becomes inefficient. The temperature rises gently to 9000-10000 K at larger radii. The left-hand panels in figure 4 show the temperature profile of the gas in the haloes and its evolution for simulation A (top panel) & B (bottom panel), respectively.
In order to understand the dynamical evolution of the gas it is illustrative to consider the thermal and turbulent velocities of the gas. In the middle panel of figure 4 we show the evolution of the radial velocity (black), the thermal velocity (blue) and the turbulent velocity (red) of the gas. As the collapse gets under way the gas develops turbulent velocities comparable to the virial velocities of the halo starting in the outer parts of the halo and progressively moving inwards. The turbulent velocities are calculated by computing the root mean square velocity of the gas after subtracting the centre of mass velocity of the halo and the velocity due to the radial inflow of the gas. We also subtract an estimate of the rotation velocity (based on the specific angular momentum and the inertia tensor as described in section 4.3) in quadrature. The turbulent velocities in the later stages of the collapse become significantly larger than the thermal velocities of the gas which we calculate as , where is the mean molecular weight, chosen to be and is the mass of the hydrogen atom. As in the simulations of W07 the turbulent velocities are supersonic. There is a net inflow of gas with velocities about a factor three smaller than the turbulent velocity component.
In the right hand panel of figure 4 we compare the ratio of the inward acceleration due to gravity to the outward acceleration due to the thermal pressure (gradient) for both simulation A (upper panel) and simulation B (bottom panel). Once the DM halo has built up, the mass of the halo is above the Jeans mass and the inward gravitational acceleration comfortably exceeds the outward acceleration due to the thermal pressure. The peak which builds up at an enclosed mass of a few times for simulation A and for simulation B is due to the formation of a centrifugally supported disc at the centre. The double peak at certain output times for simulation B is due to the presence of multiple gas clumps within the system.
In figure 5 we show the evolution of the enclosed gas and DM mass as a function of radius and the density profiles of the gas and DM distribution. The DM density (in units of ) and the DM fraction are shown for the final output times only as filled diamonds. The density profiles are calculated by averaging over spherical shells centred on the densest point in the halo. Note that in simulation B the densest point is always found in the most massive of the three clumps. The haloes collapse within 12 and 7 Myrs, respectively, after we have restarted the simulation with the increased refinement level. This corresponds to about a few free-fall times for the gas in the inner few hundred parsecs which initially has pretty much constant density. As discussed already, we initially held the refinement level at a maximum of four until the halo had a mass corresponding to ( ) DM particles. The gas completely decouples from the dark matter during the collapse and becomes self-gravitating. The density profile of the gas in the inner part steepens dramatically and settles into a close to density profile over many decades in radius as expected for an isothermal collapse (e.g. Larson, 1969, W07). The inner few times collapse further and settle eventually into rotational support. In the bottom right hand panel of figure 5 (simulation B) the secondary peak in the density profile is due to the secondary clump which eventually merges with the primary clump. The peak is absent in the final two output times as one of the smaller clumps has merged with the primary clump and the other smaller clump has had most of its mass stripped away. Note that the continuous accretion onto the virialised haloes from the filaments occurs on much longer time scales and thus has no visible effect.

Figure 8: The ratio of the square root of the smallest eigenvalue, , of the inertia tensor to the distance to the centre of the halo. Note the dip at a few times where the gas settles into a rotationally supported fat disc in simulation A. For simulation B the dip at the end of the simulation occurs at a much higher mass of .

4.3 Angular momentum loss and settling into rotational support

In Figure 6 we show the evolution of the specific angular momentum, , as a function of enclosed mass for both simulation A and simulation B. The angular momentum vector is calculated in the usual way via the cross product where is the momentum vector. We then obtain the specific angular momentum vector as a function of enclosed mass centred on the centre of mass by dividing by the enclosed mass. In both haloes there is very significant angular momentum loss. In simulation A the angular momentum drops by a factor of 20 or more between the initial output and the final output in the inner few times . The evolution in simulation B is more complicated due to the complex dynamical interaction of the triple system. Initially the angular momentum of the gas drops by a factor 20-100 in the primary clump and then increases again by a factor three to five when one of the smaller clumps merges with the primary clump.
We now want to discuss in more detail to what extent the gas settles into rotational support. During the turbulent collapse of the gas, it is not obvious how best to define rotation velocities. We follow W07 who use the ratio , where is the specific angular momentum vector and is the position vector from the centre of mass as a simple but rough measure of the rotation velocity of the gas. We also calculate a second estimate of the rotation velocity based on the angular momentum and the inertia tensor, .

The nine components of the inertia tensor are given by

(10)
(11)

and the corresponding cyclic permutations, where the sum is over computational cells. The off-diagonal components are called the products of inertia while the diagonal components are referred to as the moments of inertia. The matrix is symmetric which guarantees that the eigenvalues are real.

The angular momentum and the inertia tensor are related as,

(12)

where is the angular velocity. Using the square root of the largest eigenvalue of the inertia tensor, , we then estimate the rotation velocity as

(13)

In figure 7 we compare our two estimates of the rotation velocity (red) and (black) to the circular velocity (blue), , for simulation A (top panel) and simulation B (bottom panel). Early on in the collapse the gravitational potential is dominated by the dark matter halo and the gas is only slowly rotating with a ratio of rotation to circular velocity of about . As the gas in the inner part of the haloes collapses and becomes self-gravitating both the circular velocities and the rotation velocities rise. As shown in the right hand panel, the inner few times reach (approximate) rotational support with slightly larger than after about 12 Myrs in simulation A. In simulation B the increase in specific angular momentum at around 6.5 Myrs due to the merger of one of the small clumps with the primary clump leads to a sharp increase in the rotation velocity. As we will see in more detail later the gas nevertheless settles into a regular rotationally supported disc despite this rather violent dynamical evolution. Our estimate of the rotation velocity based on the inertia tensor exceeds the circular velocity by a factor of up to two. We will come back to this point in section 4.5.

Figure 9: A close up of the final panels of Figures 1 and 2 showing the 2D density projection of the rotationally supported gas at the centre of haloes at the end of simulation A and B.
Figure 10: The density distribution of the rotationally supported fat discs from all three simulations shown ’face-on’ (left panels) and edge on (right panels). In all cases (approximate) iso-density surfaces are plotted. Note the prominent tidal tail(s) in simulation A and the secondary clump in simulation B. The mass of the disc in simulation A is , while the disc in simulation B has a mass of , with the smaller clump having a mass of . The mass of the disc in simulation C is .
Figure 11: Left-hand Panels: The surface mass density profile of the rotationally supported gas at the centre of each halo shown in figure 10. Middle Panel: The black solid lines are the contours of the azimuthally averaged density in the plane of the simulated discs. The red contours and the blue dashed curve are the density and vertical scale height an isothermal exponential disc with the same central density, scale length and temperature would have. Right Panel: The four velocities; , , and as a function of radius.

Note that our two estimates of the rotation velocity trace each other, but that is systematically lower than by a factor 1.5-3. As the gas settles into rotational support a disc-like flattened structure should gradually form. For a turbulent collapse like the one considered here we found this to be most easily studied by looking at the evolution of the smallest eigenvalue of the inertia tensor which for a disc should be a good proxy for the thickness of the disc. In figure 8 we show the ratio of the square root of the smallest eigenvalue of the inertia tensor, , to the radius in which the inertia tensor is calculated as a function of the enclosed mass. The ratio reaches a minimum value at a few times and a few times for simulation A and B, respectively, as expected for the formation of a rotationally supported disc. The disc appears to fatten towards the centre and is surrounded by a more spherical in-fall. In simulation B the minimum moves to higher mass values as the merger of the initially three clumps progresses. At the end of the simulation the minimum occurs at .

4.4 Formation of a self-gravitating massive fat disc

About 0.1-1% of the gas in the inner part of both haloes has settled into a rotationally-supported self-gravitating object with rotation velocities of . We now progress to inspect these structures in somewhat more detail.
In figure 9 we show the final panels of figures 1 and 2 at a larger scale. The mass of the central object in simulation A in the left hand panel is . The radius of the central object is pc with a compact inner core with a radius of pc. In the right hand panel we show the central object in simulation B at the final output time. The mass of the central object (large clump) is and the radius is pc. The second clump on the right hand side has a mass of and a radius of pc. The disc in simulation C (shown in figure 3) has a final mass of and a radius of pc.
In figure 10 we show the rotationally supported objects face-on and edge-on. We have plotted iso-density surfaces using overdensities within a relatively narrow range (approximately an order of magnitude). These visualisations give a clear picture of the shape of the central objects that form in all three simulations. Each object is quite clearly a disc. The feature in the edge on visualisation of simulation A (top panel) which points north-west is the tidal tail seen at the bottom of the face-on view. There is another, less prominent, tidal tail at the top of the disc in simulation A which appears at the bottom left in the edge on view. The disc in simulation B (middle panel) is less obstructed, the smaller clump is also visible in both visualisations. The disc in simulation C (bottom panel) is the cleanest example of a disc. Like the disc in simulation B it has a somewhat ellipsoidal shape.

4.5 Properties of the centrifugally supported disc(-like object)

We now have a closer look at the surface mass density profile and the level of rotational support of the disc(-like) object for each simulation. As can be seen in the left panel of figure 11 in all three cases the surface mass density profile is exponential over several scale lengths. The scale lengths are 0.075, 0.20 and 0.27 pc in simulations A-C, respectively. At the outer “edge” of the disc at about 0.3 - 1 pc the surface mass density profile reverts to that expected for an isothermal spherical density distribution in which the discs are embedded.
In the middle panel of figure 11 we compare the density structure of the discs to that of an exponential disc which is given by

(14)

with scale height

(15)

where is the sound speed, is the mean molecular weight, is the central density, is the distance from the centre of mass and is the scale-radius of the disc Spitzer (1942); Oh & Haiman (2002). We show contours of the azimuthally averaged density in the plane of the discs. The discs (especially that in simulation A) are somewhat fatter than expected if they were supported by pressure only. This is not surprising considering the rather unrelaxed dynamical state of the discs. The general agreement of the density structure with that of an exponential isothermal disc is nevertheless remarkable. In simulation B the structure in the top right corner is due to the secondary less massive clump.

In the right panel of figure 11 we compare the actual rotation velocities of the gas in the discs (shown by the black solid curve) to our two estimates for the rotational velocity (black dot-dashed curve) and (red curve) and the circular velocity (blue curve). We compute the actual rotational velocities by rotating our coordinate system into the coordinate system of the disc using the matrix of eigenvectors obtained from the inertia tensor (which we then checked visually). The rotation velocity in the plane of the disc is then easily calculated using trigonometic transforms. The peak of the actual rotation velocities occurs at radii corresponding to, two to three times the scale radius of the exponential discs and range from 25 to 60 . The actual rotation velocities agree with our estimate from the angular momentum vector and the largest eigenvalue of the inertia tensor, to within 10-20%, despite the fact that the latter was calculated for the enclosed mass in spheres centred on the densest cell. As already mentioned the other estimate for the rotation velocities is systematically lower by a factor 1.5-3.
The actual rotation velocities exceed the circular velocities by up to a factor of 2. This is probably due to a combination of reasons. The peak rotation velocities of a thin exponential disc is about 15 % larger than that of a spherical mass configuration with the same enclosed mass. More importantly the discs in simulation B and C and probably also that in simulation A have an ellipsoidal shape which should significantly raise the rotation velocities compared to the circular velocity of a spherical mass distribution. Finally the discs have probably not yet reached centrifugal equilibrium. The gas is still settling into rotational support and has probably fallen to somewhat smaller radius initially than expected for centrifugal equilibrium and will take some time to reach centrifugal equilibrium.

4.6 Numerical Limitations

Probing the collapse of gas at the centre of dark matter haloes in a cosmological context is an extremely challenging exercise when we wish to follow the collapse to very high densities. In this work we have reached radii as small as 0.01% of the virial radius while at the same time following the dynamical evolution of a substantial fraction of the gas in the halo. This presents a considerable challenge even for an AMR code like Enzo. The main problem is that the high refinement levels necessary to achieve such a large dynamic range mean that the code will normally grind to a halt following the dynamical evolution of whatever is the first high density region to form. As discussed, in order to ameliorate this problem we have changed the refinement level during the simulation in order to allow our three haloes of choice to build up to their full mass. To what level may this have affected our results?
In order to address this we have run simulations with twice the initial refinement level and sixty four times the initial refinement level and found no systematic difference in the initial value of the angular momentum of the halo. While the detailed dynamical evolution is obviously different especially in the later stages of the collapse when the initial distribution of matter in the halo is different the qualitative behavior of the dynamical evolution did not change. We would also like to point out that our results are in most aspects similar to those of W07 who did not change the refinement level during the simulation. This meant, however, that they were not able to follow the gas at the centre of the halo settling into rotational support. We feel that this is the most interesting aspect of our work and well worth exploring the limits of the code with some manual intervention. Obviously further work is needed to address these questions in more detail.

5 Conclusions

We have investigated the dynamical evolution of the gas in haloes with virial temperatures of between K and K assuming that cooling by atomic hydrogen and helium are the dominant cooling processes. The dynamical evolution of the gas in the haloes is complex and highly turbulent with turbulent velocities approaching the virial velocity of the haloes. The gas in the inner part of our haloes collapses isothermally with a temperature of 6000-7000K on a somewhat slower time scale than the free-fall time and settles into a close to isothermal () density profile. We find no signs of efficient fragmentation confirming suggestions that the isothermal collapse of gas in a dark matter halo at temperatures close to the virial temperature of the halo leads to only modest fragmentation. The inner 0.1-1 % of the gas in the virialised haloes loses as much as 95% of its initial angular momentum during the collapse. The gas thereby collapses by a factor of 300-1000 in radius and eventually settles into a very compact rotationally supported self-gravitating disc with peak rotation velocities of 25 - 60 and “radii” of 0.3 pc - 0.6 pc (0.05% - 0.1% of the virial radius of the host halo).
The discs have an exponential surface mass density profile with scale length in the range pc which extends over several scale lengths. The vertical structure of the disc is somewhat more extended than expected for a purely pressure supported isothermal axisymmetric exponential disc.
Massive compact self-gravitating discs such as those found in our simulations have been suggested to evolve into massive seed black holes which later in the hierarchical build-up of galaxies will grow into the supermassive black holes found at the centre of the bulges of present-day galaxies. Unfortunately we were not yet able to follow the further dynamical evolution of the discs or the gas further out in the haloes in our simulations. However, independent of whether the gas in these discs will continue to efficiently lose angular momentum and contract further or will fragment and form an ultra-compact star cluster we most probably have identified an important intermediary stage en route to the formation of a massive seed black hole.

Acknowledgments

The simulations were run on the Cosmos (SGI Altix 3700) supercomputer at DAMTP in Cambridge and on the Cambridge High Performance Computing Cluster Darwin in Cambridge. Cosmos is a UK-CCC facility which is supported by HEFCE and PPARC. Darwin is the primary supercomputer of the University of Cambridge High Performance Computing Service (http://www.hpc.cam.ac.uk/), provided by Dell Inc. using Strategic Research Infrastructure Funding from the Higher Education Funding Council for England. We used the Visit (https://wci.llnl.gov/codes/visit/home.html) visualisation software to generate the visualisation plots shown here. We are grateful to Brian O’Shea, John Wise, Volker Springel, Debora Sijacki, Tirthankar Roy Choudhury, Tom Abel, Mitch Begelman, Giuseppe Lodato, Priya Natarajan, Jim Pringle and Martin Rees for useful discussions. We are also grateful to Greg Bryan, Mike Norman, Brian O’Shea and the rest of the Enzo team for making the Enzo code publicly available.

Footnotes

  1. http://cosmos.ucsd.edu/enzo/

References

  1. Abel T., Anninos P., Norman M. L., Zhang Y., 1998, ApJ, 508, 518
  2. Abel T., Bryan G. L., Norman M. L., 2000, ApJ, 540, 39
  3. Abel T., Bryan G. L., Norman M. L., 2002, Science, 295, 93
  4. Barnes J., Efstathiou G., 1987, ApJ, 319, 575
  5. Begelman M. C., Volonteri M., Rees M. J., 2006, MNRAS, 370, 289
  6. Berger M. J., Colella P., 1989, Journal of Computational Physics, 82, 64
  7. Berger M. J., Oliger J., 1984, Journal of Computational Physics, 53, 484
  8. Bett P., Eke V., Frenk C. S., Jenkins A., Helly J., Navarro J., 2007, MNRAS, 376, 215
  9. Bromm V., Coppi P. S., Larson R. B., 1999, ApJ, 527, L5
  10. Bromm V., Coppi P. S., Larson R. B., 2002, ApJ, 564, 23
  11. Bromm V., Larson R. B., 2004, ARA&A, 42, 79
  12. Bromm V., Loeb A., 2003, ApJ, 596, 34
  13. Bryan G. L., Norman M. L., 1995, Bulletin of the American Astronomical Society, 27, 1421
  14. Bryan G. L., Norman M. L., 1997, in ASP Conf. Ser. 123: Computational Astrophysics; 12th Kingston Meeting on Theoretical Astrophysics Simulating X-Ray Clusters with Adaptive Mesh Refinement. pp 363–+
  15. Bullock J. S., Dekel A., Kolatt T. S., Kravtsov A. V., Klypin A. A., Porciani C., Primack J. R., 2001, ApJ, 555, 240
  16. Carlberg R. G., 1990, ApJ, 350, 505
  17. Cavaliere A., Vittorini V., 2002, ApJ, 570, 114
  18. Colbert E. J. M., Heckman T. M., Ptak A. F., Strickland D. K., Weaver K. A., 2004, ApJ, 602, 231
  19. Croton D. J., 2006, MNRAS, 369, 1808
  20. Dekel A., Silk J., 1986, ApJ, 303, 39
  21. Di Matteo T., Colberg J., Springel V., Hernquist L., Sijacki D., 2008, ApJ, 676, 33
  22. Efstathiou G., Davis M., White S. D. M., Frenk C. S., 1985, ApJS, 57, 241
  23. Efstathiou G., Rees M. J., 1988, MNRAS, 230, 5P
  24. Fabbiano G., Kim D.-W., Fragos T., Kalogera V., King A. R., Angelini L., Davies R. L., Gallagher J. S., Pellegrini S., Trinchieri G., Zepf S. E., Zezas A., 2006, ApJ, 650, 879
  25. Fabbiano G., King A. R., Zezas A., Ponman T. J., Rots A., Schweizer F., 2003, ApJ, 591, 843
  26. Fan X., Carilli C. L., Keating B., 2006, ARA&A, 44, 415
  27. Fan X. e. a., 2001, AJ, 122, 2833
  28. Ferrarese L., Merritt D., 2000, ApJ, 539, L9
  29. Haehnelt M. G., Rees M. J., 1993, MNRAS, 263, 168
  30. Haiman Z., 2006, New Astronomy Review, 50, 672
  31. Hockney R. W., Eastwood J. W., 1988, Computer simulation using particles. Bristol: Hilger, 1988
  32. Hopkins P. F., Hernquist L., Cox T. J., Kereš D., 2008, ApJS, 175, 356
  33. Jenkins A., Frenk C. S., White S. D. M., Colberg J. M., Cole S., Evrard A. E., Couchman H. M. P., Yoshida N., 2001, MNRAS, 321, 372
  34. Kauffmann G., Haehnelt M., 2000, MNRAS, 311, 576
  35. Koushiappas S. M., Bullock J. S., Dekel A., 2004, MNRAS, 354, 292
  36. Kuhlen M., Madau P., 2005, MNRAS, 363, 1069
  37. Larson R. B., 1969, MNRAS, 145, 271
  38. Lodato G., Natarajan P., 2006, MNRAS, 371, 1813
  39. Lukic Z., Heitmann K., Habib S., Bashinsky S., Ricker P. M., 2007, ArXiv Astrophysics e-prints
  40. Lynden-Bell D., 1969, Nature, 223, 690
  41. Makishima K., Kubota A., Mizuno T., Ohnishi T., Tashiro M., Aruga Y., Asai K., Dotani T., Mitsuda K., Ueda Y., Uno S., Yamaoka K., Ebisawa K., Kohmura Y., Okada K., 2000, ApJ, 535, 632
  42. Merloni A., Heinz S., 2008, ArXiv e-prints, 805
  43. Merloni A., Rudnick G., Di Matteo T., 2004, MNRAS, 354, L37
  44. Mo H. J., White S. D. M., 2002, MNRAS, 336, 112
  45. Monaco P., Salucci P., Danese L., 2000, MNRAS, 311, 279
  46. Norman M. L., Bryan G. L., 1999, in Miyama S. M., Tomisaka K., Hanawa T., eds, ASSL Vol. 240: Numerical Astrophysics Cosmological Adaptive Mesh Refinement. pp 19–+
  47. Oh S. P., Haiman Z., 2002, ApJ, 569, 558
  48. O’Shea B. W., Bryan G., Bordner J., Norman M. L., Abel T., Harkness R., Kritsuk A., 2004, ArXiv Astrophysics e-prints
  49. O’Shea B. W., Norman M. L., 2007, ApJ, 654, 66
  50. O’Shea B. W., Norman M. L., 2008, ApJ, 673, 14
  51. Rees M. J., 1978, Phys. Scr, 17, 193
  52. Rees M. J., 1984, ARA&A, 22, 471
  53. Rees M. J., Volonteri M., 2007, in Karas V., Matt G., eds, IAU Symposium Vol. 238 of IAU Symposium, Massive black holes: formation and evolution. pp 51–58
  54. Salpeter E. E., 1964, ApJ, 140, 796
  55. Shankar F., Weinberg D. H., Miralda-Escude’ J., 2007, ArXiv e-prints, 710
  56. Sijacki D., Springel V., di Matteo T., Hernquist L., 2007, MNRAS, 380, 877
  57. Spergel D. N., Verde L., Peiris H. V., Komatsu E., Nolta M. R., Bennett C. L., Halpern M., Hinshaw G., Jarosik N., Kogut A., Limon M., Meyer S. S., Page L., Tucker G. S., Weiland J. L., Wollack E., Wright E. L., 2003, ApJS, 148, 175
  58. Spitzer L. J., 1942, ApJ, 95, 329
  59. Truelove J. K., Klein R. I., McKee C. F., Holliman II J. H., Howell L. H., Greenough J. A., 1997, ApJ, 489, L179+
  60. Volonteri M., 2007a, ArXiv e-prints, 709
  61. Volonteri M., 2007b, ApJ, 663, L5
  62. Volonteri M., Haardt F., Madau P., 2003, ApJ, 582, 559
  63. Volonteri M., Lodato G., Natarajan P., 2008, MNRAS, 383, 1079
  64. Volonteri M., Perna R., 2005, MNRAS, 358, 913
  65. Volonteri M., Rees M. J., 2005, ApJ, 633, 624
  66. Wise J. H., Abel T., 2007a, ArXiv e-prints, 710
  67. Wise J. H., Abel T., 2007b, ApJ, 671, 1559
  68. Wise J. H., Turk M. J., Abel T., 2007, ArXiv e-prints, 710
  69. Zel’Dovich Y. B., Novikov I. D., 1964, Soviet Physics Doklady, 9, 246
103750
This is a comment super asjknd jkasnjk adsnkj
Upvote
Downvote
Edit
-  
Unpublish
""
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters
Submit
Cancel
Comments 0
Request answer
""
The feedback must be of minumum 40 characters
Add comment
Cancel
Loading ...