The formation of compact massive self-gravitating discs in metal-free haloes with virial temperatures of 13000-30000 K
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
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 &
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 Enzo
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.
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,
while the circular velocity, , can be written as
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
The Hubble constant and density parameter are related to their present-day values by
where E(z) is given by
We can now define the virial temperature of the halo as
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)
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 ,
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
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.
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
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
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
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
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,
where is the angular velocity. Using the square root of the largest eigenvalue of the inertia tensor, , we then estimate the rotation velocity as
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.
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
with scale height
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
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
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.
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.
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.
- Abel T., Anninos P., Norman M. L., Zhang Y., 1998, ApJ, 508, 518
- Abel T., Bryan G. L., Norman M. L., 2000, ApJ, 540, 39
- Abel T., Bryan G. L., Norman M. L., 2002, Science, 295, 93
- Barnes J., Efstathiou G., 1987, ApJ, 319, 575
- Begelman M. C., Volonteri M., Rees M. J., 2006, MNRAS, 370, 289
- Berger M. J., Colella P., 1989, Journal of Computational Physics, 82, 64
- Berger M. J., Oliger J., 1984, Journal of Computational Physics, 53, 484
- Bett P., Eke V., Frenk C. S., Jenkins A., Helly J., Navarro J., 2007, MNRAS, 376, 215
- Bromm V., Coppi P. S., Larson R. B., 1999, ApJ, 527, L5
- Bromm V., Coppi P. S., Larson R. B., 2002, ApJ, 564, 23
- Bromm V., Larson R. B., 2004, ARA&A, 42, 79
- Bromm V., Loeb A., 2003, ApJ, 596, 34
- Bryan G. L., Norman M. L., 1995, Bulletin of the American Astronomical Society, 27, 1421
- 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–+
- Bullock J. S., Dekel A., Kolatt T. S., Kravtsov A. V., Klypin A. A., Porciani C., Primack J. R., 2001, ApJ, 555, 240
- Carlberg R. G., 1990, ApJ, 350, 505
- Cavaliere A., Vittorini V., 2002, ApJ, 570, 114
- Colbert E. J. M., Heckman T. M., Ptak A. F., Strickland D. K., Weaver K. A., 2004, ApJ, 602, 231
- Croton D. J., 2006, MNRAS, 369, 1808
- Dekel A., Silk J., 1986, ApJ, 303, 39
- Di Matteo T., Colberg J., Springel V., Hernquist L., Sijacki D., 2008, ApJ, 676, 33
- Efstathiou G., Davis M., White S. D. M., Frenk C. S., 1985, ApJS, 57, 241
- Efstathiou G., Rees M. J., 1988, MNRAS, 230, 5P
- 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
- Fabbiano G., King A. R., Zezas A., Ponman T. J., Rots A., Schweizer F., 2003, ApJ, 591, 843
- Fan X., Carilli C. L., Keating B., 2006, ARA&A, 44, 415
- Fan X. e. a., 2001, AJ, 122, 2833
- Ferrarese L., Merritt D., 2000, ApJ, 539, L9
- Haehnelt M. G., Rees M. J., 1993, MNRAS, 263, 168
- Haiman Z., 2006, New Astronomy Review, 50, 672
- Hockney R. W., Eastwood J. W., 1988, Computer simulation using particles. Bristol: Hilger, 1988
- Hopkins P. F., Hernquist L., Cox T. J., Kereš D., 2008, ApJS, 175, 356
- 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
- Kauffmann G., Haehnelt M., 2000, MNRAS, 311, 576
- Koushiappas S. M., Bullock J. S., Dekel A., 2004, MNRAS, 354, 292
- Kuhlen M., Madau P., 2005, MNRAS, 363, 1069
- Larson R. B., 1969, MNRAS, 145, 271
- Lodato G., Natarajan P., 2006, MNRAS, 371, 1813
- Lukic Z., Heitmann K., Habib S., Bashinsky S., Ricker P. M., 2007, ArXiv Astrophysics e-prints
- Lynden-Bell D., 1969, Nature, 223, 690
- 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
- Merloni A., Heinz S., 2008, ArXiv e-prints, 805
- Merloni A., Rudnick G., Di Matteo T., 2004, MNRAS, 354, L37
- Mo H. J., White S. D. M., 2002, MNRAS, 336, 112
- Monaco P., Salucci P., Danese L., 2000, MNRAS, 311, 279
- 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–+
- Oh S. P., Haiman Z., 2002, ApJ, 569, 558
- O’Shea B. W., Bryan G., Bordner J., Norman M. L., Abel T., Harkness R., Kritsuk A., 2004, ArXiv Astrophysics e-prints
- O’Shea B. W., Norman M. L., 2007, ApJ, 654, 66
- O’Shea B. W., Norman M. L., 2008, ApJ, 673, 14
- Rees M. J., 1978, Phys. Scr, 17, 193
- Rees M. J., 1984, ARA&A, 22, 471
- 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
- Salpeter E. E., 1964, ApJ, 140, 796
- Shankar F., Weinberg D. H., Miralda-Escude’ J., 2007, ArXiv e-prints, 710
- Sijacki D., Springel V., di Matteo T., Hernquist L., 2007, MNRAS, 380, 877
- 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
- Spitzer L. J., 1942, ApJ, 95, 329
- Truelove J. K., Klein R. I., McKee C. F., Holliman II J. H., Howell L. H., Greenough J. A., 1997, ApJ, 489, L179+
- Volonteri M., 2007a, ArXiv e-prints, 709
- Volonteri M., 2007b, ApJ, 663, L5
- Volonteri M., Haardt F., Madau P., 2003, ApJ, 582, 559
- Volonteri M., Lodato G., Natarajan P., 2008, MNRAS, 383, 1079
- Volonteri M., Perna R., 2005, MNRAS, 358, 913
- Volonteri M., Rees M. J., 2005, ApJ, 633, 624
- Wise J. H., Abel T., 2007a, ArXiv e-prints, 710
- Wise J. H., Abel T., 2007b, ApJ, 671, 1559
- Wise J. H., Turk M. J., Abel T., 2007, ArXiv e-prints, 710
- Zel’Dovich Y. B., Novikov I. D., 1964, Soviet Physics Doklady, 9, 246