Numerical resolution effects on simulations of massive black hole seeds

We have performed high-resolution numerical simulations with the hydrodynamical AMR code Enzo to investigate the formation of massive seed black holes in a sample of six dark matter haloes above the atomic cooling threshold. The aim of this study is to illustrate the effects of varying the maximum refinement level on the final object formed. The virial temperatures of the simulated haloes range from $\rm{T} \sim 10000\ \rm{K} - 16000\ \rm{K}$ and they have virial masses in the range $\rm{M} \sim 2 \times 10^7 \rm{M_{\odot}}$ to $\rm{M} \sim 7 \times 10^7 \rm{M_{\odot}}$ at $z \sim 15$. The outcome of our six fiducial simulations is both generic and robust. A rotationally supported, marginally gravitationally stable, disk forms with an exponential profile. The mass and scale length of this disk depends strongly on the maximum refinement level used. Varying the maximum refinement level by factors between 1 / 64 to 256 times the fiducial level illustrates the care that must be taken in interpreting the results. The lower resolution simulations show tentative evidence that the gas may become rotationally supported out to 20 pc while the highest resolution simulations show only weak evidence of rotational support due to the shorter dynamical times for which the simulation runs. The higher resolution simulations do, however, point to fragmentation at small scales of the order of $\sim 100$ AU. In the highest resolution simulations a central object of a few times $10^2\ \rm{M_{\odot}}$ forms with multiple strongly bound, Jeans unstable, clumps of $\sim 10\ \rm{M_{\odot}}$ and radii of 10 - 20 AU suggesting the formation of dense star clusters in these haloes.


INTRODUCTION
It is now widely accepted that supermassive black holes (SMBHs) populate the centres of most if not all galaxies. SMBHs were invoked early on in the literature to explain the powering of extremely luminous, extra-galactic sources (Lynden-Bell 1969) initially dubbed quasi-stellar objects (QSOs) (Zel'Dovich & Novikov 1964;Salpeter 1964). Since their discovery nearly five decades ago much theoretical and observational work has gone into explaining both their presence and their properties. The mass of central SMBHs appears to correlate surprisingly strongly with the luminosity and the stellar velocity dispersion of the galactic bulges hosting them (Magorrian et al. 1998;Gebhardt et al. 2000;Ferrarese & Merritt 2000, see Kormendy & Ho 2013 for a recent review). The physical origin of the claimed tight correlation of stellar velocity dispersion and black hole mass, the "M BH − σ " relation, is still controversial (see Fabian 2012 for a recent review), ⋆ E-mail:john.regan@helsinki.fi but is generally attributed to the joint formation history of (proto-)galaxies and their central black hole . Numerical simulations that include models for the feedback from the SMBHs on the surrounding gas support this picture (e.g. Di Matteo et al. 2005;Sijacki & Springel 2006;Sijacki et al. 2007Sijacki et al. , 2009Johansson et al. 2009a,b;Debuhr et al. 2011;Choi et al. 2013).
Surprisingly, billion solar mass black holes already exist at z > ∼ 6 when the Universe was less than a Gyr old (Fan 2004;Fan et al. 2006;Mortlock et al. 2011;Venemans et al. 2013). This poses problems for models that assume Eddington limited growth of stellar mass seed black holes in the limited time available (e.g. Costa et al. 2013 andsee Volonteri 2010;Haiman 2013 for recent reviews). Stellar mass seed black holes form in shallow potential wells and their growth is hampered by the (negative) feedback due to photo-ionisation heating and the energy and momentum injection due to supernovae (Haiman et  Eddington limited growth have been further compounded by recent simulations reporting more efficient fragmentation than earlier work lowering the expected masses of population III (POP III) stars to ≈ 10M ⊙ − 100M ⊙ Clark et al. 2011;Stacy et al. 2012;Greif et al. 2012). This increased fragmentation is due to a combination of turbulence, H 2 collisional dissociation cooling and collision-induced emission which could not be followed in earlier simulations due to resolution constraints.
Growth from massive seed black holes with masses above 10 4 M ⊙ that form in a rapid phase where the Eddington limit is exceeded by a large factor is therefore often invoked as an alternative (Begelman & Rees 1978;Begelman et al. 1984;Loeb & Rasio 1994;Haiman & Loeb 2001;Begelman 2001;Oh & Haiman 2002;Bromm & Loeb 2003). The most promising route for the formation of such massive seed black holes is probably "direct collapse" in dark matter haloes with virial temperatures above the atomic cooling threshold in which cooling by molecular hydrogen and metals is not important and fragmentation is therefore suppressed (Shlosman et al. 1989;Loeb & Rasio 1994;Begelman et al. 2008;Begelman 2008;Begelman et al. 2006;Volonteri & Begelman 2010). A number of studies has shown that it may indeed be plausible that a sufficient number of such haloes has not been enriched by metals and are subject to sufficiently strong (local) UV radiation that molecular hydrogen is dissociated (Bromm & Loeb 2003;Haiman 2006;Mesinger et al. 2006;Dijkstra et al. 2008;Cen & Riquelme 2008;Ahn et al. 2009;Shang et al. 2010;Wolcott-Green et al. 2011;Tanaka & Li 2013). Note, however, recent work by Aykutalp et al. (2013) on (negative) photo-ionisation feedback driven by X-rays has also shown that the composition of the gas surrounding the seed black hole may play a role in determining the final outcome.
Simulations of the isothermal collapse of the gas in such haloes are numerically challenging, but significant progress has been made mainly, but not exclusively, with grid based adaptive mesh refinement (AMR) codes (e.g. Johnson & Bromm 2007;Wise et al. 2008;Regan & Haehnelt 2009b;Greif et al. 2008;Johnson et al. 2011;Latif et al. 2012Latif et al. , 2013cPrieto et al. 2013;Latif et al. 2013a). In the simulations the gas cools efficiently to temperatures of about 7000-8000K and becomes strongly turbulent. Turbulent angular momentum transport thereby leads to substantial and sustained mass inflow and allows the gas to settle into an isothermal (ρ ∝ r −2 ) density profile. A major difficulty with these simulations thereby stems from the fact that the dynamical timescales get shorter as the collapse progresses and the simulations are generally able to follow only an increasingly smaller fraction of the gas as it collapses to the highest densities. This makes it very difficult to follow gas that has settled into angular momentum support further out in the halo while the inner part continues to collapse. This has led to discussions as to the extent of angular momentum support and the subsequent fragmentation scale (Wise et al. 2008;Regan & Haehnelt 2009b;Latif et al. 2013b). Making progress with these questions is important for judging the the importance of the "direct collapse" model as a route to the formation of supermassive black holes. We therefore follow on here from our previous studies of this problem with a larger suite of simulations performed with increased resolution with a newer version of the ENZO code. In this way we are able to follow the collapse to higher density as well studying to an unprecedented extent the longer term evolution of (marginally) angular momentum supported gas during the collapse. The larger sample of simulations also allowed us to look for systematic trends in the environmental properties of the haloes. The aim being to search for clues within the large scale environment which may favour the formation of massive seed black holes.
The paper is structured as follows. In §2 we describe the details of the numerical simulations. In §3, appendix A & §4 we describe the results of our numerical simulations and in §5 we deliver our conclusions. Throughout this paper we assume a standard ΛCDM cosmology with the following parameters (Planck Collaboration et al. 2013, based on the latest Planck data), Ω Λ,0 = 0.6817, Ω m,0 = 0.3183, Ω b = 0.0463, σ 8 = 0.8347 and h = 0.6704. We further assume a spectral index of primordial density fluctuations of n = 1.

The Adaptive-Mesh Refinement Code ENZO
We have used the publicly available adaptive mesh refinement (AMR) code ENZO 1 . The code has matured significantly over the last few years and as of July 2013 is available as version . Throughout this study we use ENZO version 2.2. ENZO was originally developed by Greg Bryan and Mike Norman at the University of Illinois (Bryan & Norman 1995Norman & Bryan 1999;O'Shea et al. 2004;The Enzo Collaboration et al. 2013). The gravity solver in ENZO uses an N-Body particle mesh technique (Efstathiou et al. 1985;Hockney & Eastwood 1988) while the hydro calculation are performed using the piecewise parabolic method combined with a non-linear Riemann solver for shock capturing. One of ENZO'S greatest strengths lies in the fact that additional finer meshes can be laid down as the simulation runs to enhance the resolution in a given, user defined, region. 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 (1995) 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. Numerous chemistry solvers are available as part of the ENZO package. For our purposes we use only the six species model which includes: H, H + , He, He + , He ++ , e − . In this case we are neglecting the effects of H 2 cooling -see §2.5 for more details and limitations of this approach. We also allow the gas to cool radiatively during the course of the simulation. Our simulations make extensive use of ENZO'S capability to employ nested grids.

Nested Grids & Initial Conditions
For our fiducial simulations the maximum refinement level is set to 18. Our fiducial box size is 2 h −1 Mpc comoving giving a maximum comoving resolution of ∼ 6 × 10 −2 h −1 pc. At the typical redshifts we are interested in (z ∼ 15) this corresponds to a physical resolution of ∼ 4×10 −3 h −1 pc. Furthermore, during the course of the resolution study conducted here ( §4) we increase the resolution up to a maximum refinement level of 26. This leads to a maximum comoving resolution of ∼ 2 × 10 −4 h −1 pc, corresponding to ∼ 1 h −1 AU physical at z ∼ 15. Initial conditions were generated with the "inits" 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 256 3 and no nested grids. In these exploratory simulations we have identified the most massive halo at a redshift of 15 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 1024 3 . 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. Within the highest resolution region we further restrict the refinement region to a comoving region of size 128 h −1 kpc around the region of interest so as to minimise the computational overhead of our simulations. We do this for all of our simulations. The total number of particles in our simulation is 4935680, with 128 3 of these in our highest resolution region.  Table 1 gives further details on the simulations discussed here.

Comparison to Previous Work
In comparison to our previous work in this area (Regan & Haehnelt 2009b) these new simulations differ in several ways. Firstly, the code used is a newer version of ENZO (version 2.2), secondly, the methodology of the simulations is somewhat different. While in Regan & Haehnelt (2009b) we first allowed a large halo to build up before increasing the refinement, we do not do that here and instead allow maximum refinement from the start of the simulations and follow the gravitational collapse consistently. Thirdly, we increase the resolution of our highest resolution simulations significantly -by a factor of 256. Finally, we modified the refinement criteria and some of the runtime options used by ENZO compared to the Regan & Haehnelt (2009b) simulations.

Refinement Criteria
ENZO uses adaptive grids to provide increased resolution where it is required. For the simulations discussed in this paper we have used three refinement criteria implemented in ENZO : DM overdensity, baryon over-density, and Jeans length. The first two criteria introduce additional meshes when the over-density ( ∆ρ ρ mean ) of a grid cell with respect to the mean density exceeds 3.0 for baryons 2 The virial mass is defined as 200 times the mean density of the Universe in this case -see §3.3 for further details. and/or DM. The third criterion has received quite a lot of attention in recent years with several groups arguing for a higher threshold, than the canonical Truelove criterion of 4 cells (Truelove et al. 1997). Federrath et al. (2011), Turk et al. (2012 and Latif et al. (2013c) have all shown that simulations require a minimum Jeans length resolution of 32 cells in order to obtain converged turbulent energy results. We follow Latif et al. (2013c) and Meece et al. (2013) in using a Jeans resolution criteria of 64 cells in order to fully resolve the turbulent energy requirement in the highest density regimes. We set the MinimumMassForRefinementExponent parameter to −0.1 making the simulation super-Lagrangian and therefore reducing the threshold 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). When this option is selected the code defines a minimum temperature for each grid at the highest refinement level. This minimum temperature is that required to make each grid Jeans stable multiplied by the above parameter. This parameter was introduced into ENZO to alleviate artificial fragmentation and angular momentum non-conservation -see Machacek et al. (2003) for further details. We further make use of the RefineRegionAutoAdjust parameter setting it to one. This parameter modifies the refinement region during the course of a simulation so as to allow only the highest resolved dark matter particles into the refinement region. Larger mass dark matter particles are therefore excluded. We then allow our simulations to evolve restricted only by the maximum refinement threshold.

Zero Metallicity and H 2 Dissociation
Assuming zero metallicity and efficient dissociation of H 2 (Wise et al. 2008;Regan & Haehnelt 2009b;Johnson et al. 2011;Prieto et al. 2013;Latif et al. 2013b,c,d,e) strongly reduces fragmentation and therefore optimises conditions for direct collapse of gas into a massive (seed) black hole. As discussed extensively in the literature these conditions may be fulfilled in a small fraction of haloes above the atomic cooling threshold which may nevertheless be sufficiently abundant to act as seeds for the most massive supermassive black holes at high redshift (e.g. Loeb & Rasio 1994;Bromm & Loeb 2003;Spaans & Silk 2006;Dijkstra et al. 2008). Accurate modelling of the metal enrichment history of these haloes is very challenging. Cen & Riquelme (2008) have argued that the mixing of metals produced by early episodes of population III (POP III) star formation is very inefficient, but as pointed out by Omukai et al. (2008) even small amounts of metal contamination may be enough to induce fragmentation. Efficient dissociation is Density slices through the centre of the central object in each halo. The "up" vector is chosen to be the angular momentum vector and so we are looking down onto the central object. Each image is scaled in density and length in the same way. The colour scale for the density runs from 10 −19 g cm −3 to 10 −15 g cm −3 . It is clearly visible that in each case a disk is formed, with well defined spiral arms. The size of the plotted region in each panel is 1pc.
required to allow the build up of a halo to T vir > 10000 K without collapse. In order to efficiently dissociate H 2 a strong UV background is required. Recent work by several groups (Shang et al. 2010;Wolcott-Green et al. 2011;Johnson et al. 2013) has determined that a Lyman-Werner flux of ∼ 1000 J 21 for a UV background with a POP III 3 spectrum or ∼ 30 -300 J 21 for a POP II spectrum is required to allow for the build up of an atomic cooling halo, where J 21 is the canonical background intensity in units of 10 −21 erg cm −2 s −1 Hz −1 sr −1 . Here we also make the assumption that in our simulations both metal and H 2 cooling are inefficient. Such haloes will be rare, but may exist when a local strong UV source exists 10 kpc from the collapsing halo which augments the global UV background in the early Universe (Dijkstra et al. 2008;Agarwal et al. 2012).

The Suite of Simulations
We have performed 6 simulations of haloes with virial temperatures in the range of 9000 K to 16000 K and virial velocities in the range of 15 to 22 km/s collapsing at redshifts 15 < z < 22 all within a 2 Mpc h −1 box starting at z = 100. The details of each simulation are shown in Table 1. Each simulation was run to a maximum refinement level of 18 with an initial root grid of dimension 128 grid cells. This results in a maximum comoving resolution of 5.96 × 10 −2 h −1 pc. The physical resolution achieved in each run is given in Table 1. Once the simulations reach the maximum refinement level the collapse is able to continue adiabatically. Given the short dynamical times at the maximum refinement level the evolution of all, but the central regions has effectively stopped. As gas continues to fall onto the central regions, thereby increasing the densities, the cooling times of the gas continue to decrease. As the maximum refinement level has already been reached the code is no longer able to calculate accurately the gas properties and errors in the hydro solver begin to appear. At this point we terminate the simulation, discarding any spurious outputs. In addition to the above simulations with a maximum refinement level set to 18 we further run simulations A and C with 4 different maximum refinement settings (see Table 2 for more details). We are thus able to see the effect of numerical resolution on the central object obtained. We discuss this further in §4. Note that unlike in our previous simulations (Regan & Haehnelt 2009b), where we only started refinement when a simulated halo was well above the atomic cooling threshold, we have followed here the growth of the virial temperature of the six individual haloes across the atomic cooling threshold with the code always allowed to refine where necessary.

Simulation Morphologies
In Figure 1 we show visualisations of the central pc of our simulations containing about 10 4 M ⊙ at a time when the gas is settling into (marginal) rotational support. In all cases unstable discs have formed which develop prominent spiral features. As in our previous simulations and that of Latif et al. (2013c,b,e) 10 4 M ⊙ appears to be the characteristic mass scale for which the collapse becomes marginally rotationally supported and angular momentum support first sets in. Note, however, that as we will discuss later the innermost part of the halo nevertheless continues to collapse. The panels are density slices taken through a plane perpendicular to the angular momentum vector, calculated at the core of the halo. We are therefore looking down onto the central 10 -3 10 -2 10 -1 10 0 10 1 10 2 10 3 Radius (pc) 10 -3 10 -2 10 -1 10 0 10 1 10 2 10 3 Radius ( show a more complex morphology. As we will see in more detail later, formation of a (marginally) dynamically stable and rotationally supported object is the generic outcome of our simulations.

Properties of the Gas in the Dark Matter Haloes when Collapse Occurs
The collapse of the gas in our simulations starts when the virial temperature of the halo has reached the atomic cooling threshold. The virial temperatures of our haloes are therefore all in the narrow range of 9000 K -16000 K with an average virial temperature at which collapse begins of T vir ∼ 11000 K. The virial radius of the collapsed object is approximately 0.5 kpc in all simulations. The virial quantities are defined such that the density at the virial radius is 200 times the mean density of the Universe at that redshift. For the haloes found here this corresponds to a virial mass of a few times 10 7 M ⊙ (Mo & White 2002). That all of our haloes have similar masses is due to the maximum refinement criteria we have imposed.
Once collapse begins, additional refinement is engaged by the code. This extra refinement allows us to track the collapse down to sub parsec scales, but effectively freezes any further growth of the halo. Hence, in a simulation such as this, with a relatively high maximum refinement level, following the growth of haloes to much higher masses is currently extremely challenging. We will return to this issue in §4 where we conduct further simulations with different maximum resolutions.

Profiling the Central Object
As cooling is facilitated only via atomic hydrogen cooling in the simulations performed here, the gas cannot cool below approximately 7000 K. In Figure 2 (left panel) we show the temperature of the gas over several decades in radius. We plot the temperature out to approximately the virial radius which is well outside the realm of the collapse. The temperature of the gas is found by averaging the temperature of cells in spherical shells outwards from the densest point in the halo. As expected, the temperature remains approximately constant as the density grows towards the centre of the halo. Initially the gas is shock-heated to T ∼ 10 4 K, close to the virial radius, from where it cools via atomic hydrogen transitions to T ∼ 7000 K. Figure 2 (right panel) shows the gas density profile of the halo over the same range. The density profile is initially quite flat at, or outside, the virial radius, but quickly steepens to attain a slope of n ∝ r −2 as expected for an isothermal collapse. The profile is not completely smooth due to the presence of small dense clumps within the halo (as shown in Figure 1). All simulations show similar profiles with the maximum number density obtained in each simulation found to be ≈ 1 × 10 11 cm −3 . The result here are in very good agreement with those of our own previous simulations as well as those of other authors (Wise et al. 2008;Latif et al. 2013c,b,e). The maximum density reached is a function of the maximum resolution of the simulation. As we will see in §4 higher resolution simulation are able to follow the collapse to higher densities at the expense of failing to track the further evolution of the gas at lower resolution (lower density).
In the left hand panel of Figure 3 we have plotted the thermal velocity (blue lines), turbulent velocity (red lines) and the radial velocities (black) lines against enclosed mass. Again we plot out to the virial mass. The thermal velocity is computed as V T H = 3k B T /M, where k B is the Boltzmann constant, T is the temperature of the gas in a given shell and M is the gas mass in that shell. 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. Finally, the radial velocity of the gas is determined by computing the radial component of the Cartesian velocities in each shell.
Near the virial radius the turbulent velocities are stable at approximately 20 km s −1 and similar to the virial velocity of the DM The velocity profile for three distinct velocities. In blue we have plotted the thermal velocity profile. In red the turbulent velocity profile and finally in black the radial velocity profile of the gas. Right Panel: The central object loses angular momentum allowing the gas to fall to the centre. In Regan & Haehnelt (2009b) we showed the evolution of a central object and found that it can lose ≈ 90% of its initial angular momentum.
10 -3 10 -2 10 -1 10 0 10 1 10 2 10 3 Radius (pc)  halo, the thermal velocity is very low and radial in-fall is constant at approximately 15 km s −1 . As the gas falls towards the centre the turbulent velocities grow and exceed the thermal and radial components. All of the simulations show significant turbulent velocities consistent with their complex morphology as seen in Figure 1. The radial velocities have a mean value of close to 10 km s −1 , for the majority of outputs. Simulations C, D and F have radial velocities which display a significant positive contribution at an enclosed mass of about 10 4 M ⊙ . At this point the central object becomes rotationally supported and the radial velocity switches sign. Significant amounts of gas are flung out from the central core along the spiral arms resulting in average positive values of the radial velocity. The same is seen -albeit to a lesser extent -in each of the other simulations.
In the right hand panel of Figure 3 we plot the angular mo-mentum of the central object against the enclosed gas mass in spherical shells. The angular momentum is calculated in spheres centred on the densest point in the simulation. As discussed in detail in Regan & Haehnelt (2009b) the innermost gas shells "lose" up to 90% of their initial angular momentum during the turbulent collapse probably due to a combination of angular momentum redistribution and angular momentum cancellation. The low angular momentum gas that ends up flowing to the centre then settles into a rotationally supported object. We determine the rotational properties of the collapsing gas in a similar way as in Regan & Haehnelt (2009b). We calculate the inertia tensorĨ and calculate its eigenvectors to describe the principal axes of the rotation of the central object. 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, a 1 , we then estimate the rotation velocity as As in Regan & Haehnelt (2009b) we then use the square root of the smallest eigenvalue of the inertia tensor as a proxy for the thickness of the flattened object formed.
In the left panel of Figure 4 we show the enclosed gas mass against radius for the six simulations. In each simulation the enclosed mass first decreases linearly inwards as expected for an isothermal density profile. The gas at these radii is virtually in freefall, collapsing turbulently at a significant fraction (≈ 25 − 30%) of the free-fall speed. Between 1 and 0.1 pc the curve of enclosed mass flattens indicating that the collapse becomes marginally rotationally supported. The onset of rotational support is also reflected in the right panel of Figure 4 where we plot the ratio of the square root of the smallest eigenvalue of the inertia tensor, a3, and the radius at that point. Note that the plot focuses on the range of mass shells where the rotational support sets in between 5 × 10 3 M ⊙ and 5 × 10 5 M ⊙ . All of the simulations show a characteristic dip in "thickness" between ∼ 1 × 10 4 M ⊙ and ∼ 1 × 10 5 M ⊙ . The dip is due the formation of a flattened disk-like structure. Interestingly simulations A, C & E show the most prominent dip at a few times 10 4 M ⊙ . This corresponds to a radius of ∼ 1 pc. Figure 1 clearly shows a flattened object with spiral arms for simulation A, but simulations C & E have a more complex morphology. Nonetheless the calculations indicate the presence of a fat disk (since the ratio of a3/R is relatively large). Simulation B shows no prominent dip compared to the surrounding gas, rather the ratio is constant at ∼ 0.3. Simulation D shows a double dip at enclosed masses of ∼ 3 × 10 4 M ⊙ and ∼ 1 × 10 5 M ⊙ consistent with the visual impression from Figure 1. Simulation F shows a wide dip at ∼ 4 × 10 4 M ⊙ which again is well matched by the visual impression from Figure 1.

Rotational Support and Gravitational Stability
We now investigate the actual level of rotational support in the collapsing gas in our simulations by comparing the rotational velocity to the circular velocity, In the left panel of Figure 5 we show the circular velocity for our six haloes and in the right panel we show the ratio of rotational and circular velocity plotted against radius. The rotational velocity is computed using Eq. 2. Between ∼ 1 pc and ∼ 1 × 10 −2 pc the ratio of rotational to circular velocity reaches a value greater than 1.0 in all simulations -right panel of Figure 5. To investigate the stability of the disks formed we also calculate their surface mass density and the corresponding Toomre parameter as a function of radius. Similar to Regan & Haehnelt (2009b) we find exponential surface mass density profiles (left panel of Figure 6) with scale lengths in the range R d ∼ 0.24 − 0.53 pc. The Toomre stability parameter shown in the right panel of Figure 6 is calculated as (Toomre 1964), where c s is the sound speed, κ = √ 2 (V rot /r) (1+d ln V rot /d ln r) 1/2 (Oh & Haiman 2002;Binney & Tremaine 2008) is the epicyclic frequency, V rot is the rotational velocity, Σ is the surface mass density and r is the radius. For values of Q < 1 the disc is expected to be gravitationally unstable. Simulations C, E & F all have values Q > 1 within a scale length or more. Simulations A, B & D are approaching gravitational stability, but do not have values > 1.0 at the end of the simulation. These results should not be surprising given the rather unrelaxed dynamical state of the discs as they undergo further collapse and evolution. Overall the discs are found to be marginally stable with some fragmentation and formation of clumps during the course of the collapse. Small scale turbulence effects within the disk, where Q 1, may furthermore They found that while the general properties of a collapse are unchanged with the subgrid model, fragmentation is increased slightly and the formation of discrete clumps is more pronounced. We will discuss the formation of bound clumps further in §4.6. The morphology of the marginally rotationally supported gas varies rapidly between different simulations, but the formation of a central rotationally supported (marginally) gravitationally stable disc appears to be a robust outcome of our simulations.

Effects of Environment on Halo Angular Momentum Dissipation
As noted in the introduction we conducted a study using our six fiducial simulations on the environmental dependence of haloes on early angular momentum dissipation. We investigated the dependence of the early dissipation of gas angular momentum on halo "rareness" which we quantified using the peak height relative to the RMS fluctuation of the density field of each halo. As discussed in more detail in Appendix A we could not identify any significant trend with the rareness of the halo.

Varying the Resolution
The maximum spatial (physical) resolution of our fiducial simulations discussed so far is ≈ 10 −3 pc and the minimum dark matter particle mass is M DM = 8.301 × 10 2 M ⊙ . We investigate now the effect that increasing and decreasing the maximum refinement level has on the outcome of our simulations. For this we rerun simulations A & C with different maximum refinement levels. We have also decreased the minimum dark matter particle mass for the highest refinement runs. The parameters of the simulations used for this resolution study are shown in Table 2. In Table 2 simulation A12 refers to simulation A run with a maximum refinement level of 12, A18 is simulation A run with a maximum refinement level of 18 (which was the fiducial value used in this study), etc. The same notation applies to simulation C. Note that the simulations run at a refinement level of 12 are able to run for considerably longer as they do not resolve the collapsing structures as early as the higher resolution simulations. All simulations were run until the maximum refinement level was reached and the simulations were no longer able to properly compute the hydrodynamics due to the lack of resolution.

Splitting of Dark Matter Particles
When we run the higher resolution simulations in Table 2 with higher refinement levels we have also decreased the minimum mass of the dark matter particles in these simulations. This was to ensure that the dark matter particles do not introduce any artificial fragmentation in the higher resolution simulations. In the higher resolution runs, those with a maximum refinement of > 18, the gas mass resolution becomes substantially higher. We initially ran simulations without altering the dark matter particle mass, but found that the dark matter particles were aiding fragmentation at small scales in some cases. We therefore utilised a modified form of the dark matter splitting algorithm in ENZO -which follows the prescription given in Kitsionas & Whitworth (2002). The initial dark matter particle mass was M DM = 8.301 × 10 2 M ⊙ . The splitting algorithm splits the dark matter particle into 13 new particles each of mass M DM = 6.385 × 10 1 M ⊙ . We did this for simulations A22, A26, C22 and C26. Simulations A22 and A26 were split and restarted at z = 22 while simulations C22 and C26 were split and restarted at z = 18. In order to fully utilise this feature of ENZO we had to modify the ENZO code that contains the dark matter splitting algorithm to suit our specific needs. The patch was then upstreamed to the ENZO mainline.

Results of the Resolution Study
In Figures 7 and 8 we show a density projection of the output at the end of each of the highest resolution runs for simulations C and A ( Table 2). The projection is centred on the point of highest density   −1 ), the virial temperature f (K), the DM mass g within the virial radius (M ⊙ ), the maximum number density h in the halo (cm −3 ), the dark matter particle mass i (M ⊙ ) and the spatial resolution j (pc) of the simulation. in each case. The projection is created using the YT analysis suite (Turk et al. 2011). To create the projections the radiative transfer equation along the line of sight is integrated by converting field values to emission and absorption values producing a final image. Fragments are easily identified in this projection. As the maximum refinement level and therefore the maximum resolution of the simulation increases, more and more of the computational resources are directed at the densest clump(s) collapsing first. The outer collapsing mass shells "freeze out" more quickly with increasing maximum refinement level. This strongly affects the resulting morphology. Figure 7 shows density projections of the results from our high resolution run of simulation C with a maximum refinement level of 26 at two different zoom levels. The left panel has a size of 1 pc while the right panel shows the same time output, but zoomed in by a factor of ≈ 200. The right hand panel displays clear evidence for the onset of fragmentation within a marginally stable disc with a mass of a few times 10 2 M ⊙ , which we will explore in more detail in subsection 4.6. Similarly, Figure 8 shows the density projection for simulation A26. Here a single dense clump of gas (marked by the white circle) fragments at the edge of the main structure (Mass ∼ 10 4 M ⊙ ) and collapses to high density early before the rest of the gas can evolve further. ENZO follows this small fragment at the expense of the evolution of the outer mass shells. The width of the left panel is 1 pc. The zoomed in visualisation on the right, with a width of 0.01 pc, is centred on this clump which has a mass of ≈ 10 M ⊙ and is an order of magnitude less massive than the clump which collapses at the centre of simulation C26.
In Figure 9 we show how the radial profiles of enclosed mass, density, disc thickness and ratio of rotational to circular velocity are affected by the choice of maximum refinement level. As expected with increased refinement level the collapse can be followed to higher density. The maximum density reached in simulation C26 is > 1 × 10 16 cm −3 , similar to that in the highest resolution runs presented in Latif et al. (2013b). Note that for the highest densities the gas should have become optically thick to Thompson scattering and our assumption that the gas is optically thin to cooling radiation breaks down even for continuum radiation (see section §4.5 for more discussion). The top left panel of Figure 9 shows the enclosed mass against radius. C12, C14, C18 and to a somewhat smaller extent C22 all show plateaus in the radial enclosed mass profile which indicate that the collapse becomes marginally rotationally supported. The enclosed mass where this happens decreases, however, with increasing resolution as the code follows an increasingly smaller fraction of the gas at increasing resolution. Note also the absence of such a plateau in simulation C26. In C26 the time steps have become so short that the settling into rotational support cannot be tracked for any of the outer mass shells. In C26 a disk does not form at all, because it does not have adequate time to do so. The dips in the top right panel of Figure 9 for C12, C14, C18 and C22 showing the ratio of the eigenvalue a 3 of the inertia tensor to the radius as a proxy for disk thickness suggests that indeed a (fat) disk has formed in the lower resolution simulations. C26 shows no obvious dip. In the bottom right panel of Figure 9 C12 and C14 show that marginal rotational support is achieved at scales of ∼ 10 − 20 pc, C18 and C22 show strong rotational support at scales of ∼ a few times 10 −1 and ∼ a few times 10 −2 pc respectively. C26 shows no sign of rotational support -the ratio goes only above 1.0 near the resolution limit and is therefore not indicative of a settled disk. Figure 10 shows the radial profiles for simulation A. Again with increasing resolution higher densities are reached, but as we can clearly see here it is not necessarily the central regions of the halo which are collapsing first. The decoupling of the small offcentre fragment from the rest of the gas has a noticeable effect on the profile in simulations A22 and A26. However, in A22 the fragment rejoins the collapsing outer structure as the simulation progresses. In the top left panel of Figure 10 we show the enclosed mass as a function of radius. A26 shows no sign of disk formation and the simulation follows only the collapse of the off-centre clump to ever higher densities at late times. It is also worth noting that in the right hand panel of Figure 6 where we plot the Toomre stability parameter, the disk in simulation A18 is observed to be highly gravitationally unstable while the disk in simulation C18 appears to be (marginally) gravitationally stable within a scale length or more. This suggests that at a resolution of ∼ 1 × 10 −3 pc (maximum refinement = 18) the disk in simulation A is already highly gravitationally unstable, but the resolution of the simulations is not sufficient to detect the onset of fragmentation.

Choosing the Right Resolution
The above results demonstrate that choosing the correct maximum resolution for collapse simulations is of crucial importance. Both choosing a resolution that is too low and equally choosing a resolution that is too high can lead to a misleading interpretation of results. Our lower resolution simulations, C14 for example, provide tentative evidence that rotational support may extend out to ∼ 10 pc or more, but they lack the required resolution to probe the dynamics at sub-parsec scales. In C26 we showed that at very high resolution a rotationally supported disk does not form. This is because the simulation was unable to follow the outer mass shells for the required dynamical time in order for the disks to form. The dynamical time required to form the rotationally supported disk in C18 is greater than 10 Myrs, however, C26 only evolves gas above a similar density for the order of 1 to 2 Myrs. The high resolution runs did however, indicate that fragmentation is common at sub-parsec scales in agreement with the results found by Latif et al. (2013b). In addition, running simulations that closely resembles our highest resolution simulations Latif et al. (2013b) also found that a disk does not form. In both cases, the reason is most likely that the simulation does not run for the required dynamical time of the outer shells rather than the fact that a disk does not form at all. This is an inherent limitation of high resolution AMR simulations that can only be overcome by running simulations at varying resolutions.
Increasing the resolution in simulations without a feedback mechanism which will set a characteristic resolution scale, means that an AMR code will follow the densest fluctuation at the expense of lower density gas. Feedback mechanisms such as turbulence (Latif et al. 2013b), magnetic fields (Latif et al. 2013a), or ionising radiation from a previous source naturally provide a mechanism to prevent or reduce fragmentation. In the case of such a feedback mechanism high resolution simulations are extremely useful, because the simulation is able to track the dynamics of the in-falling gas at both very small and intermediate scales. The feedback sets a characteristic resolution scale by acting against density fluctuations at very small scales. Without such feedback mechanisms converged numerical results are impossible to reach in this case. This is due to the nature of the AMR simulations conducted. Multiple simulations at varying resolutions are therefore required to ascertain the complete picture.

Limitations of the Six Species Cooling Model
As already discussed our simulations have included only the six species chemical model ( H, H + , He, He + , He ++ , e − ). We have neglected the presence of H 2 in our simulations and have thus implicitly assumed a dissociating Lyman-Werner background capable of destroying all the H 2 within our halo. As discussed in §2.5 a background capable of keeping a halo H 2 free is certainly plausible at this redshift for selected haloes in the Universe, at least in the initial stages of the collapse. However, at the densities reached in our high resolution simulations it is likely that H 2 will be formed, through both the usual H − and H + channels and via three body reactions involving H (Palla et al. 1983), at the very centre of our collaps-ing object, possibly at a rate faster than the dissociation timescale. Shang et al. (2010) ran several simulations with varying H 2 dissociating backgrounds. Their strongest fields, 10 3 J 21 , effectively reduce the H 2 fraction to 10 −8 , which is several orders of magnitude below the level required to influence the gas thermodynamics. However, the densities reached in their simulations were significantly below those reached in our high resolution runs. Radiative transfer methods will be required to accurately model this environment as we exit the regime in which the optically thin approximation holds.
In simulation C26, for example, the highest densities reached are n ∼ 10 16 cm −3 . As already noted in §4.3 we have now moved outside the realm in which the optically thin approximation is valid and hence we are no longer accurately modelling the cooling processes which are activated at these high densities (e.g. line trapping of cooling radiation, Thomson scattering). This will have the knock-on effect of impacting the subsequent cooling and possibly the evolution of the central object. In a future study we will investigate in detail the effects of a Lyman-Werner dissociating background, fully taking into account the effect a dissociating flux has on the H 2 cooling and therefore on the gas thermodynamics.

The Onset of Fragmentation and the Formation of Bound Clumps
In the absence of metal and molecular hydrogen cooling fragmentation is strongly suppressed in collapsing dark matter haloes at the atomic cooling threshold (Oh & Haiman 2002;Bromm & Loeb 2003;Lodato & Natarajan 2006;Spaans & Silk 2006;Begelman et al. 2006;Schleicher et al. 2010). This is because the cooling effectiveness of atomic hydrogen lines drops below ∼ 10000 K. The temperature remains at close to 7000K for gas outside of the core of the halo, but can drop to between 3000 K and 5000 K at the very centre of the highest density regions. As the collapse becomes marginally rotationally supported and the gas settles into a disc like structure, our simulations begin, however, to show the onset of significant fragmentation and the formation of clumps. In Figure 11 we show a time series of outputs at fixed density and fixed size of our simulation C26. Several clumps are clearly noticeable. Note that the maximum circular velocity of the marginally stable disc at this stage is 40 km/s. The temperature of the gas is therefore well below the "virial temperature" of the disk which may explain the onset of fragmentation at this stage. Fragmentation may also be aided by the interaction of the gas that is flung out with the gas flowing inwards from larger radii (see Bonnell 1994 for a similar discussion in a different context). The clumps have typical masses of 5 -20 M ⊙ and could therefore be potential formation sites of stars. Let us now have a closer look at clump A -identified in Figure 11. The time at which the clump is first visually identified is selected as T = 0 and shown in the top left panel. The clump merges with another clump in the third picture of the sequence and the two clumps coalesce before circling back towards the dynamical centre. The physical properties of clump A are shown in Figure 12. Using the clump finding algorithm (Smith et al. 2009) built into YT we track the clump over several dynamical times. In the bottom left panel of Figure 12 we plot the density against radius for the clump. As the clump moves away from the central density and collapses the profile steepens significantly between T = 0 yrs and T = 212 yrs. The radius of the clump is approximately 20 AU. In the top left panel we show the enclosed mass profile. The enclosed mass grows approximately linearly with the radius out to about 20 AU. The bottom right hand panel shows the circular velocity of the clump against radius. The circular velocity increases strongly between 20 and 40 AU, reaching a maximum value of ≈ 25 − 30 km s −1 . The clump is thus strongly bound. In the top right hand panel we plot the ratio of the enclosed mass divided by the Jeans mass. Initially the clump is Jeans stable, but as the clump evolves, and through its merger with another clump, the clump grows in mass, becomes Jeans unstable, and collapses.

Predicting the Further Evolution
It is very difficult to follow the long-term evolution of the bulk of the gas in the haloes at high resolution due to the prohibitively short dynamical time scales. As we have discussed earlier, there are, however, clear signs for the onset of fragmentation in the simulations, despite the absence of metal and molecular hydrogen cooling. The long term fate of these fragments is uncertain at this point, some of them appear to dissolve quickly, but some of them appear to become strongly bound and are likely to persist and form stars.
Predicting the final outcome of these simulations even for the simplified case of no radiative and supernova feedback and no optical depth effects is thus not yet possible and will probably have to involve simulations making use of sink particles. Judging from our current simulations, it appears that it is rather unlikely that most of the gas at the centre of these haloes will rapidly accrete onto a single star/black hole and that a dense star cluster may form instead which then evolves further (Regan & Haehnelt 2009a). However, the inclusion of radiative transfer and feedback effects may well again change the outcome not only quantitatively, but also qualitatively.

CONCLUSIONS
We have performed a suite of ENZO AMR simulations of the collapse of gas in DM haloes with virial temperatures just above the atomic cooling threshold. As in previous simulations by us and other authors the highly turbulent gas attains an isothermal density profile during the collapse. The gas looses angular momentum efficiently due to angular momentum separation and cancellation and settles towards the centre at a significant fraction of the free-fall velocity. Furthermore, within our suite of simulations, we ran an extended set of simulations at varying maximum refinement level. We found that the gas becomes marginally rotationally supported at characteristic radii, settles and forms a thick often gravitationally unstable disc. The mass at which the onset of rotational support can be studied by such AMR simulations is set by the refinement level and decreases with increasing maximum refinement level. At the highest refinement levels for which we ran our simulations the dynamical times become too short to study the onset of rotational support at any radius. The interpretation of these AMR simulations requires therefore great care. Similar to what was found in the simulations by Latif et al. (2013c,b,e) we found that in our highest resolution simulations (maximum refinement level 26) the gas reaches number densities as high as 10 16 cm −3 and is prone to fragmentation. The fragments have masses of a few tens of solar masses and radii of 10-20 AU. These dense often strongly bound clumps with circular velocities > 10 km/s forming in the highest resolution simulations are, however, not necessarily located at the very centre of the halo. They are probably best interpreted as the onset of fragmentation within the marginally stable, highly self-gravitating, discs which appear to be a generic outcome of these simulations. Predicting the further evolution of the collapse is not yet possible with the current setup, but the formation of a dense stellar cluster, perhaps as an intermediate stage to the formation of a massive seed black hole is a possible outcome. Future progress will require the inclusion of feedback effects due to the radiation of massive stars and supernovae and probably also the employment of appropriately chosen sink particles in order to follow the evolution of the unstable discs over longer timescales.

A THE PROBABILITY DISTRIBUTION OF THE INITIAL ANGULAR MOMENTUM AND ITS DEPENDENCE ON HALO RARENESS AND ENVIRONMENT
As discussed in the introduction the space density of massive seed black holes required to explain the presence of the observed billion solar mass black holes at very high redshift is very much smaller than the space density of DM haloes at the atomic cooling threshold. Low metallicity and a high amplitude of the Lyman-Werner radiation may therefore not necessarily be the only criteria required for the (efficient) formation of a massive seed black hole. As we were able to perform here a larger sample of simulations of these haloes compared to our previous work (Regan & Haehnelt 2009b) we have looked also briefly into the possibility of systematic trends with the "rareness" and the environment of the simulated haloes (see also Prieto et al. 2013). We thereby characterise the "rareness" of the halo (or peak height) in terms of the RMS fluctuation amplitude of the linearly extrapolated density field where δ c is the threshold over-density from Press & Schechter (1974), D(z) is the growth factor and σ (M) is the mass fluctuation inside a halo of mass M: where the integral is over the wavenumber k, P(k) is the power spectrum and W (kR) is the top hat window function. As the haloes undergo collapse, it is the lowest angular momentum gas that can fall radially to the centre (Dubois et al. 2012;Bellovary et al. 2013). As the collapse goes on angular momentum redistribution and cancellation occurs (Wise et al. 2008;Regan & Haehnelt 2009b). In order to increase the range of ν beyond that of our fiducial simulations we ran an extra 150 extra dark matter simulations to find a rarer peak than in our six fiducial simulations. We picked the highest ν peak from these 150 runs and reran the simulation at the same maximum refinement level as the fiducial simulations. This high-ν halo is included in the following results. Figure 13 shows the PDF of the projected angular momentum for our seven simulated haloes with virial temperatures of ∼ 5000 K. The projected angular momentum is calculated by first rotating into the coordinate system defined by the inertia tensor. The left panel of Figure 13 shows the PDF for each halo when it has reached a virial temperature of T vir ∼ 5000 K. At T vir ∼ 5000 K the halo is still growing by mass inflow along the surrounding filamentary structures and is just about to collapse. The angular momentum distribution at this stage should be a good representation of the initial angular momentum distribution with which the gas enters into the collapse. Figure 13 shows no clear trend of lower angular momentum gas with rareness of the halo. The angular momentum distributions are pretty symmetrically distributed at zero with similar amounts of co-and counter-rotating gas which explains why the gas "looses" so efficiently angular momentum during the collapse due to angular momentum cancellation. There is also no evidence of a trend towards more symmetric distribution in rarer haloes as suggested by Dubois et al. (2012), but the range of ν of our haloes is still rather moderate due to the small size of the simulation box. Unfortunately, with our current set-up of the simulations it is not straightforward to extend the range further. The right hand panel of Figure 13 shows the mass fraction of absolute projected angular momentum above a given value. Similarly no trend with the rareness of the halos is seen. Figure 14 shows the PDF and mass fraction plots for the outputs at the end of our simulation runs (i.e. when T vir 10000 K). Again no clear trend with ν is detected. Other properties of the seven haloes we have looked at show likewise no significant trend.
In Figure 15 we illustrate the location of the central object amid the cosmic web for two of the simulations. The left panel of Figure 15 shows the central object in simulation A while the right hand panel shows the location of the central object in simulation B. The object in simulation A is located at the knot of a cosmic web of gas while the central object in simulation B is located midway along a filament. Both objects therefore have formed in relatively different environments with regard to gas inflow and in particular filamentary inflow. The visualisations also show the presence of other dense clumps of gas forming within 5 kpc of the central object in each simulation. We are not able to confirm the trends in different collapse properties of haloes at the atomic cooling threshold found by Prieto et al. for our (admittedly smaller) sample of simulations. The difference in resolution between the studies is, however, significant. The study undertaken by Prieto et al. (2013) used much larger box sizes (up to 8 Mpc 3 ), but with reduced resolution (8 pc at maximum). Therefore, they were able to study the dynamical evolution of the central object over much longer periods, but were not able to resolve the sub parsec evolution of these objects.