The Collapse of Atomically-Cooled Primordial Haloes. I. High Lyman-Werner Backgrounds

Pristine, atomically-cooled haloes may be the sites of primordial quasar formation because atomic cooling triggers rapid baryon collapse that can create 10$^4$ - 10$^5$ M$_{\odot}$ black hole seeds. However, no numerical simulation has ever followed the collapse of these haloes for the times required to form supermassive stars and direct-collapse black holes (DCBHs). We have now modeled baryon collapse in atomically-cooled haloes with a wide range of spin parameters and assembly histories for times that are sufficient for DCBH formation. Fragmentation of accretion disks after $\sim$ 500 kyr is nearly ubiquitous in these haloes and in most cases leads to the formation of binary or multiple supermassive stellar systems. They also confirm that rapid baryon collapse proceeds for the times required for these stars to form DCBHs. Our simulations suggest that binary or even multiple DCBH formation was the rule rather than the exception in the primordial Universe.


INTRODUCTION
Hot, atomically-cooled primordial haloes at z ∼ 15 -20 may have been the birthplaces of the earliest quasars in the universe, more than 200 of which have now been discovered at z > 6 (e.g., Fan et al. 2003), including nine at z > 7 (Mortlock et al. 2011;Bañados et al. 2018;Matsuoka et al. 2019;Yang et al. 2020;Wang et al. 2021).In this picture, primordial haloes grow to masses of 10 7 -10 8 M and reach virial temperatures of ∼ 10 4 K without ever having formed a primordial (Pop III) star, either by being immersed in strong Lyman-Werner (LW) UV backgrounds that destroy all their H2 (e.g., Agarwal et al. 2012;Yue et al. 2014;Dijkstra et al. 2014;Johnson et al. 2014;Schauer et al. 2015Schauer et al. , 2017a) ) or in highly supersonic baryon streaming motions that delay the collapse of the halo even if H2 is present (Tseliakhovich & Hirata 2010;Greif et al. 2011;Stacy et al. 2011;Schauer et al. 2017b;Hirano et al. 2017).Virial temperatures of 10 4 K trigger atomic cooling that leads to catastrophic baryon collapse at infall rates of ∼ 0.1 -1 M yr −1 .Stellar evolution models predict that such flows, if they persist, would E-mail: spatrick@ed.ac.uk build up cool, red supermassive stars (SMSs) before dying as 100,000 -300,000 M direct-collapse black holes (DCBHs; e.g., Hosokawa et al. 2013;Umeda et al. 2016;Woods et al. 2017;Haemmerlé et al. 2018a;Woods et al. 2021b;Herrington et al. 2023).However, it has now been shown that the rare, turbulent haloes that have been shown to form quasars by z 6 in large-scale cosmological simulations (Di Matteo et al. 2012, 2017;Lupi et al. 2019;Valentini et al. 2021;Lupi et al. 2021) produced DCBHs without UV backgrounds, supersonic baryon streaming motions, or even atomic cooling (Latif et al. 2022).
DCBHs are currently the leading candidates for the seeds of the first quasars because ordinary Pop III star BHs are only a few tens to hundreds of solar masses at birth (e.g., Hirano et al. 2014Hirano et al. , 2015) ) and form in low densities that preclude rapid initial growth (Whalen et al. 2004;Kitayama et al. 2004;Alvarez et al. 2009;Whalen & Fryer 2012;Smith et al. 2018).In contrast, DCBHs can grow much more rapidly because they are born in dense environments in massive host haloes capable of retaining their fuel supply, even when it is heated by X-rays (Johnson et al. 2013).Runaway stellar collisions in dense, marginally-enriched clusters can create BHs of up to a few thousand solar masses (Devec-chi & Volonteri 2009;Latif et al. 2016;Sakurai et al. 2017;Reinoso et al. 2018;Boekholt et al. 2018) but even these objects may not be massive enough to become quasars by z > 6 (Smidt et al. 2018;Latif & Khochfar 2020;Zhu et al. 2020 -see also Woods et al. 2019;Maio et al. 2019).
Although numerical simulations of rapid baryon collapse in primordial atomically-cooling haloes have steadily improved over the past decade, they remain a trade-off between resolution and evolution time.The first simulations either resolved sub-AU scales that could capture the formation of the protostar but not the accretion disk around it on parsec (pc) scales (Wise et al. 2008) or larger scales that could follow the evolution of the disk for a few dynamical times but not the SMS at its center (Regan & Haehnelt 2009a,b).These studies found large central infall rates at the onset of collapse but could not determine how long they lasted.Later work at high resolution and somewhat longer evolution times found that these infall rates continued down to scales approaching those of the star but did not run for nearly enough times to follow its evolution (Latif et al. 2013a;Regan et al. 2014).
The implementation of pressure floors (Machacek et al. 2001) and sink particles at high resolution confirmed that accretion rates could remain large at the smallest scales and persist for a few tens of thousands of years (Latif et al. 2013c;Shlosman et al. 2016;Chon et al. 2016).Chon et al. (2018) and Regan & Downes (2018b) used sink particles with radiative feedback from the star and a simple prescription for its evolution to follow collapse for 100 kyr and 250 kyr, respectively.They found that radiation from the star was unable to suppress accretion, as did Ardaneh et al. (2018) and Luo et al. (2018), and reported some small-scale fragmentation in the disk at later times (see also Becerra et al. 2015Becerra et al. , 2018)).Chon et al. (2018) found that almost half of the fragments in one of their haloes form in binaries that they suspected could later become SMSs, corroborating idealised simulations by Bromm & Loeb (2003).However, neither study evolved the disks for long enough times to determine if any of the fragments became stars or were simply subsumed into the central object later on.Suazo et al. (2019) examined the collapse of atomically cooled haloes at intermediate resolutions in high LW backgrounds for ∼ 600 kyr, longer than previous studies but still well short of the formation of a DCBH.
The large inflow rates lasting for millions of years required to form SMSs and DCBHs have only recently been confirmed to occur in numerical simulations (Latif et al. 2020;Regan et al. 2020).These models sacrificed the extreme resolution of earlier studies in order to capture the evolution of the gas at the center of the halo over many dynamical times.Latif et al. (2020) found that the disk can fragment into binary and multiple SMS systems, but only in a few haloes with high spin parameters that were special cases.Here, we follow the collapse of an ensemble of atomically-cooling haloes over the full range of spin parameters expected for these objects and for a variety of merger histories.We evolve them for up to 3 Myr, enough time for SMS formation and collapse at their centers.For simplicity, we approximate high LW backgrounds by deactivating H2 chemistry in our runs and defer intermediate LW backgrounds with H2 to later studies.This should be considered to be a strong but useful upper limit to actual LW backgrounds in the primordial Universe because H2 could self-shield in the cores of massive haloes and survive even intense UV fluxes from nearby Pop III stars (Shang et al. 2010;Latif et al. 2014bLatif et al. , 2015)).Even small amounts of H2 could cool gas in the core and result in much less massive stars.In Section 2 we describe our numerical simulations and discuss halo evolution, disk dynamics and central accretion rates in Section 3. We conclude in Section 4.

NUMERICAL METHOD
We first performed a series of low-resolution dark-matter (DM) only runs to identify haloes with a range of assembly histories and spin parameters.These haloes are then resimulated with full baryonic physics at much higher resolution to follow their collapse and the dynamics of the disks at their centers.

Enzo
We use the Enzo adaptive mesh refinement (AMR) cosmology code to model the haloes in our study (Bryan et al. 2014).Enzo has an N −body adaptive particle-mesh scheme for evolving DM (Efstathiou et al. 1985;Couchman 1991;Bryan & Norman 1997) that is self-consistently coupled to hydrodynamics and nonequilibrium primordial gas chemistry (Anninos et al. 1997;Glover & Abel 2008).We employ the piecewise parabolic method for gas dynamics (PPM; Woodward & Colella 1984;Bryan et al. 1995) and apply the HLLC scheme for enhanced stability with strong shocks and rarefaction waves (Toro et al. 1994).Our simulations include six species (H, He, e − , H + , He + , He 2+ ) with no H2 in order to ensure isothermal cooling and collapse.Our chemistry model includes collisional excitational and ionizational cooling by H and He, recombination cooling, bremsstrahlung cooling, and inverse Compton cooling by the cosmic microwave background (CMB).

Simulation Setup
Haloes of interest are identified from a series of unigrid DMonly simulations in 1.5 h −1 Mpc cosmological boxes with a resolution of 256 3 and initial conditions generated at z = 200 with MUSIC (Hahn & Abel 2011).We use the second-year Planck cosmological parameters: ΩM = 0.308, ΩΛ = 0.691, Ω b h 2 = 0.0223, σ8 = 0.816, h = 0.677 and n = 0.968 (Planck Collaboration et al. 2016).The Rockstar halo finder (Behroozi et al. 2013) is used to find haloes with masses of a few times 10 7 M in these boxes.We restart our simulations at z = 200 with a top grid resolution of 256 3 and three additional nested grids centered on the halo (each with resolution of 256 3 ) which yield effective DM and baryon mass resolutions of 28 M and 34 M , respectively.Up to 15 levels of refinement are allowed for a maximum physical resolution of 0.014 pc.A grid is flagged for refinement when baryon or DM overdensities exceed 8.0 times the mean density (δρ/ρ > 8.0).We resolve the Jeans length by 32 cells to capture the formation of turbulent structures and prevent artificial fragmentation, and we activate refinement on Jeans length at z = 30 (Truelove et al. 1997;Latif et al. 2013a).
DM particles are smoothened at the tenth refinement level, which corresponds to a comoving resolution of 5.72 pc, to avoid spurious effects caused by the discrete sampling of the DM potential.Below pc scales, collapse is dominated by baryons, and if there is an insufficient number of DM particles at higher resolutions they can unphysically accelerate and heat gas in their vicinity.Smoothing mitigates such effects.We also turn on a pressure floor at the maximum level of refinement in our simulations that prevents collapse on scales below our maximum resolution by ensuring that the cells on that level are Jeans stable.The pressure floor in our runs smoothens the gravitational potential over 4∆x, where ∆x is the length of a cell at the highest level of refinement (Machacek et al. 2001), so clumps can be smeared out on those scales.Pressure floors never reverse collapse in our models, just prevent it at the highest level of refinement, and they never result in local expansion of the gas that drives shocks, especially at temperatures of 8000 K.These points are corroborated by previous simulations of collapse in atomically-cooled haloes with pressure floors and with sink particles without floors that produced essentially identical results (Latif et al. 2013c(Latif et al. , 2022)).We output the data every 10 kyr during the run and follow the gravitational collapse for about 3 Myr.

Halo Selection Criteria
We initially selected 20 haloes from our DM-only runs that reached masses of at least 10 7 M by z ∼ 14 -20.From this set we chose only eight haloes for resimulation because of their merger histories and spin parameters at the onset of atomic cooling.They primarily grew through accretion, major mergers or some combination of the two, where a major merger is defined to be the collision of two haloes with a mass ratio of 1/5 or more.They were also chosen so that they spanned the likely range of spin parameters, λ, for atomically-cooled haloes, where λ is (Peebles 1969;Bullock et al. 2001) Here, L, R, and M are the angular momentum, virial radius and virial mass of the halo, respectively.Values for λ for all eight haloes are listed in Table 1.The remainder of the 20 haloes were discarded for having merger histories similar to those in our final set.We save outputs at redshift intervals ∆z = 0.5 from z = 25 -15, use Rockstar to locate positions, masses and velocities of the haloes in each redshift bin, and then port them to ytree (Smith & Lang 2019) to build merger trees.Assembly histories for these haloes, designated 1, 2, 8, 10, 12, 16, 19 and 20, are shown in Figure 1.The trees extend from z = 25 at the bottom to z = 15 at the top in increments ∆z = 0.5.Each circle represents an ancestor of the final halo and has a radius proportional to its logarithmically-scaled mass.This proportionality varies between haloes in Figure 1 so the same circle radius in different trees does not imply the same mass.The main progenitor line is marked in red.Halo 12 is an example of a halo that has grown primarily by accretion while halo 16 has three major mergers.The others lie somewhere between these two extremes.Most of the haloes have lines of low-mass haloes in their ancestry, visible as the long vertical black tails with small black dots along them.Halo 1 grows slowly at early times but then   suddenly rises in mass at z = 17.5 because of a surge in accretion and a merger.Halo 2 has the highest accretion rate of the sample but does not begin to cool before reaching a mass of 8.47 × 10 7 M at z = 14.5, likely because of turbulence driven by rapid infall.Halo 8 begins to cool at the earliest times at a mass of 1.15 × 10 7 M at z = 20.4.Halo 10 has a major merger early in its history and has the highest spin parameter at the time it begins to cool at z = 17.3 at a mass of 2.6 × 10 7 M .Halo 12 has the next highest spin when it begins to cool at z = 16.8 at a mass of 1.93 × 10 7 M , which is interesting because it is almost entirely the product of accretion, not mergers that could have spun it up by tidal torquing.Halo 16 undergoes three major mergers at z = 18.5, 17.5, and 16.5 and minor mergers at z=18.5 and 15.It has the lowest accretion rates in our sample.Halo 19 begins to cool at the lowest redshift, z = 13.9 at a mass of 2.12 × 10 7 M , and has the lowest spin at the onset of cooling.Like halo 12, which has the second highest spin parameter, it too is mostly the result of accretion, not mergers.Halo 20, which has the second lowest spin parameter, is also mostly built up by accretion but has a number of minor mergers in its history.We summarise these results in Table 1. In sum, haloes 1, 10 and 16 have  10 -2 10 -1 10 0 10 1 10 2 10 3 radius (pc)  undergone major mergers during their formation, haloes 2, 12 and 19 have grown mainly via accretion, and haloes 8 and 20 are products of both minor mergers and accretion.Redshifts at the onset of atomic cooling are determined in our runs with both gas and DM for accuracy.

RESULTS
Rapid baryon collapse triggered by atomic cooling leads to the formation of large, rotationally-supported disks at the centers of all eight haloes.We now discuss the evolution of these disks, their accretion rates, and their stability over time.

Initial Baryon Collapse
We show spherically-averaged density and temperature profiles for our haloes at the onset of cooling and collapse in Figure 2. The densities all closely approximate the Larson-Penston solution, ρ ∝ r −2.2 , which slightly deviates from the classic r −2 density profile expected for self-gravitating isothermal spheres because of the presence of DM in the halo.The virial shock of the halo is visible as the bump in density and temperature at 200 pc -800 pc in the profiles.Densities range from 10 6 -10 7 cm −3 in the core to 0.1 -1 cm −3 at the virial shock.They level off at small radii because of the absence of torques within the Jeans sphere (not because of pressure floors, which are imposed just after this stage of collapse).The large bump in density at ∼ 200 pc in halo 2 is due to another halo that is about to merge with it.
The temperature profiles show that collapse due to atomic cooling is nearly isothermal, with flat ∼ 8000 K temperatures extending out to ∼ 100 pc.Gas falling onto the halo is first heated to 10,000 K as it passes through the virial shock and then cools to ∼ 8000 K as it settles deeper into the gravitational potential well of the halo.Because of our resolution limits, gas never reaches densities where it becomes optically thick to Lyα photons (e.g., Smith et al. 2017;Smith & Bromm 2019).At higher levels of refinement these photons could become partly trapped by frequent resonant scatterings, reducing cooling efficiencies and increasing temperatures deep in the core.We also note that at number densities n 10 8 cm −3 H − cooling becomes important and would reduce gas temperatures to ∼ 5000 K before the gas becomes optically thick to Lyα, further increasing densities and potentially trapping more cooling photons at later times.

Disk Formation / Evolution
Rapid baryon collapse leads first to the formation of a dense core and then a rotationally-supported, self-gravitating disk with an initial radius of ∼ 0.2 pc by ∼ 200 kyr.A bar instability then appears and creates two spiral arms.The disk thereafter can exhibit a range of morphologies, as we discuss below, but continues to grow in mass and radius and generally becomes more turbulent over time.This turbulence increases the effective alpha parameter of the disk and it grows to over 1 pc in radius in all eight haloes.
Density projections of the disks in haloes 10, 12 and 1 are shown in Figure 3.They represent the range of disk evolution in our runs: relatively quiescent ones, as in halo 10; more turbulent disks that undergo multiple episodes of fragmentation, as in halo 12; and the prompt formation of interacting binary disks at early times, as in halo 1.In the initial stages of disk formation the collapse of low-angular momentum gas creates the small, dense cores like the one at 51 kyr in halo 12.They are soon followed by the appearance of barred spiral arms as gas at higher angular momentum rains down onto the disk, as shown at 226 kyr in halo 10.Since radiation from the star cannot prevent gas from con- tinuing to fall onto the disk at high rates, it grows rapidly in radius, in some cases to ∼ 3 pc by the end of the runs.Neither the bars nor viscous processes can transport angular momentum out from the center of the disk as fast as gas rains down on it so its spiral arms cannot maintain a single angular velocity and they begin to fracture into clumps, generally by 500 -700 kyr but as early as ∼ 250 kyr in halo 1. Disk fragmentation thus proceeds here via the spiral arm instability (SAI) mechanism discussed in Inoue & Yoshida (2020) in which spiral arms form first and then become locally unstable, overcoming the stabilizing effects such as gas pressure and Coriolis force.Most fragments spiral into the center of the disk but some persist for times that are sufficient to form satellite disks, as we discuss below.Tidal torquing between fragments and the disk often cause the center of the disk to eject tidal tails when the clumps crash into it.Turbulence in the disks also breaks up spiral arms to some degree.
The disk in halo 10 is the most stable one and does not begin to fragment until 1.14 Myr, when it has grown to ∼ 2 pc in radius and become more turbulent.It thereafter forms a number of clumps, but they are quickly torn apart by gravitational torques and taken back up into the spiral arms except for one that forms at 1.26 Myr about 1.2 pc from the center.This clump grows in mass over time, collapsing into a second, smaller disk that is in a highly elliptical orbit with the first.This second disk exchanges mass with the first in a series of close encounters by the end of the run at 2.01 Myr.Because the second disk appears well after the first one, either a DCBH -SMS binary system or a binary DCBH could form after the second star collapses.
There is more turbulence in the disk in halo 12 at early times and it produces two episodes of violent fragmentation, shown at 1.17 Myr and 1.85 Myr in Figure 3. Between episodes, most of the short-lived clumps are torn apart by gravitational torques and ejected in the tidal streams visible at 1.43 Myr, but one clump flattens into a secondary disk that persists from 889 kyr -1.39 Myr.Consequently, it is likely that a second SMS will form in this system, as discussed below.The star at the center of the main disk likely collapses to a DCBH before the second episode of clumping, so some of the clumps would rapidly fuel the growth of the BH when they crash into the center of the disk.
The disk in halo 1 exhibits morphologies similar to the others at first but begins to fragment at early times, 289 kyr after formation.This episode is short-lived, however, because of the formation of a second, satellite disk at 366 kyr at an initial distance of ∼ 1 pc from its more massive companion.It orbits its companion until being torn apart and subsumed into the main disk at 947 kyr, which results in the ejection of the multiple tidal tails shown in Figure 3 at 995 kyr.The orbit is highly eccentric and precesses around the main disk, with disk-disk separations varying from ∼ 0.5 -1.2 pc.Mass exchange between the two disks occurs a number of times during close encounters, as shown at 879 kyr.Tidal torques between the two disks suppress subsequent fragmentation in both, with only two short-lived clumps forming after the appearance of the second disk (fragmentation can also be pre-empted as the satellite disk sweeps up gas along its orbit).The survival of the second disk for ∼ 600 kyr leads to the likely formation of a second SMS that is later drawn into a close binary with the star in the main disk, with final separations of less than 0.1 pc at the time they collapse to DCBHs.Large accretion rates in both disks corroborate the likely eventual formation of a binary DCBH, as discussed below.

Binary and Multiple Disk Systems
Although we noted the appearance of binary accretion disks in haloes 1, 10 and 12 in Section 3.2, they form in all the haloes and in some cases multiple disks appear.For example, three additional disks appear in halo 2: from 666 kyr -1.04 Myr, 937 kyr -1.48 Myr and 1.18 Myr -1.32 Myr (right panel of Figure 4).There is almost no overlap in time between the first satellite disk and the other two, which coevolve and even briefly exchange mass with each other at one point.In halo 8 just one binary disk appears in the run from 685 kyr -1.0 Myr.In halo 16, which hosts a relatively stable disk like that in halo 10, a fragment forms at 1.75 Myr, is flung into a highly elliptical orbit by three-body interactions with the main disk and another clump at 2.01 Myr, and then flattens into a satellite accretion disk by 2.7 Myr that is still orbiting the original disk at the end of the simulation at 3.05 Myr.Two additional disks form in halo 19, from 0.675 kyr -907 kyr and 907 kyr -1.41 Myr.Finally, in halo 20, two disk fragments collide at 639 kyr and form a dense clump that flattens into a disk by 774 kyr.It grows and soon rivals the original disk in mass (see the left panel of Figure 4) but is then almost destroyed in close encounters with this disk twice, at 1.36 Myr and at 1.61 Myr.In both cases tidal streams are ejected from the site of the interaction and the satellite disk retains only a small fraction of its mass after 1.61 Myr.
We show center-to-center separations between the main disk and its satellite over time in haloes 8 and 20 in falls somewhere in between these two extremes.Separation distances vary from ∼ 0.3 pc -2 pc, and the satellite disks typically make 5 -10 orbits before either crashing into the center of the main disk or being torn apart by gravitational torques.
The binary disks in our models live for shorter times than those in Latif et al. (2020).Their disks survive for longer times because their simulations probed the upper extremes in spin parameter for the host haloes.Their large angular momenta allowed the satellite disks to orbit the main disks for longer times before crashing into them.In contrast, the haloes in our study have a range of spin parameters that do not reach their upper or lower extremes so satellite disks are taken up into the main disks in less time.

Accretion Rates
We plot accretion rates for all the disks in our haloes in Figure 6.They are calculated by dividing the change in mass in a 0.135 pc sphere centered on the densest cell at the center of the disk over 10 kyr intervals by 10 kyr, the time between data dumps.We use this radius to ensure that the tally sphere is always resolved by at least 10 zones while excluding the spiral arms of the disk, which could produce spurious contributions to the infall rates.We also found that the accretion rates converged at this radius as we decreased the radius.The initial time t = 0.0 is when the densest cell first reaches the maximum refinement level and catastrophic collapse has begun.

Main Disk
In each case there is an initial jump in accretion of 0.3 -1 M yr −1 for 200 -300 kyr that is due to the formation of the dense core of the disk, like the one shown at 51.5 kyr in halo 12 in Figure 3.Rates then fall as conservation of angular momentum of the gas increases its rotational velocity and it flattens into a self-gravitating, rotationally-supported disk, which hinders further collapse.Over the next few hundred kyr, accretion can be relatively smooth as in haloes 10 and 16, highly turbulent as in haloes 1, 8 and 19, or clumpy as in haloes 2, 12 and 20 because of fragmentation of the disk and collision of the fragments with its center.The larger jumps in accretion rate at 300 -600 kyr are from turbulence (200 -400 kyr in halo 8), mass exchange with a smaller, satellite disk (500 -600 kyr in halo 1), or clumps crashing into its center (500 -600 kyr in halo 12).The largest peaks, like those at 1.0 Myr in halo 1 and 1.4 Myr in halo 19, are from the collision of a satellite disk with the main disk.
After the initial jump in accretion due to the formation of the dense core of the disk, the rates in some of the disks fall to lower values and flatten out as gas flows circling the center of the disk become Keplerian.This effect is most prominent in the disks in haloes 10 and 16.Other disks are too turbulent for the flows to ever become Keplerian and tidal interactions with satellite disks or collisions of fragments with the center of the disk disrupt circularised flows at later times, even in disks 10 and 16.
At times, the accretion rate can become negative if turbulence drives gas out of the center of the disk (and thus out of the tally sphere) or during mass exchange with another disk in a close encounter.This process occurs on much larger scales than infall onto the star itself, which would not be expected to lose mass during these episodes.During such times accretion onto the star would at most be temporarily halted.Average accretion rates after the initial surge vary from ∼ 0.1 -0.5 M yr −1 but can peak at up to 2 M yr −1 for up to 100 kyr.Our rates are consistent with those reported at earlier times in previous studies ( 0.1 M yr −1 ; e.g., Latif & Volonteri 2015).Note that these average infall rates would produce SMSs that collapse via the general-relativistic instability rather than depletion of their hydrogen fuel ( 0.04 M yr −1 ; Woods et al. 2017).
We show the mass accumulated at the center of the main disk in each halo in the left panel of Figure 7.The final masses at the end of the runs vary from 2.4 × 10 5 M to over 5 × 10 5 M .The occasional dips in mass correspond to the episodes of negative accretion rate visible in Figure 6.Comparison of these masses to the final SMS masses plotted as a function of accretion rate in Figure 4 of Woods et al. (2017) indicates that all the disks have been evolved for sufficient times for SMSs at their centers to have collapsed to DCBHs prior to the end of the run.

Satellite Disks
Accretion rates for the longest-lived binary disks in our haloes are shown in red in Figure 3.Although they only live for a fraction of the time that the main disk evolves, their infall rates can rival those of the main disk.Peaks in accretion in the satellite disks often coincide with dips in the main disk and peaks in the main disk often coincide with dips in the satellite.These correlations arise from mass exchange between the disks during close encounters, like the ones at 540 kyr and 879 kyr in halo 1 in Figure 3.Other peaks in accretion rate are due to collisions between the satellite disk and other clumps orbiting the main disk.Both mass exchange and collisions can perturb the orbits of binary disks.Accretion in the binary companion disks ends when they collide with the center of the main disk or are destroyed by tidal forces.
The masses accumulated at the centers of the binary disks range from approximately 75,000 M to nearly 200,000 M and are shown in the right panel of Figure 7.At our resolution (and without a detailed stellar evolution model) it is not clear how much of this mass is taken up into a star at the center of the disk, but it is likely that at least some of these disks host SMSs.As noted earlier, in some cases such as halo 19 the second SMS would form early and coevolve with the star at the center of the main disk, possibly producing an SMS binary when the second disk later crashes into the first.In other cases, such as halo 16, the second SMS would likely form after the first collapses to a DCBH and an SMS -DCBH binary would form.

Angular Momentum
The evolution of the angular momentum, L, in the central 0.136 pc of the main disk in all eight haloes is shown in Figure 8.It mostly tracks the enclosed mass (compare with Figure 7), implying that no large external torques are exerted on the disk that could change the angular momentum without mass transfer.The large jumps in L coincide with the collisions of massive clumps or satellite disk with the main disk, which also coincide with the large jumps in accretion rates in the disks.The subsequent drops in L are due to the ejection of tidal streams after the collision that transport L out of the center of the disk.There is no apparent correlation between halo spin at the onset of atomic cooling and large jumps in L because they are present in halo 19, which has the lowest of the spin parameters, and halo 12, which has the second-highest spin.Likewise, there is no obvious connection between previous mergers and sudden increases in L because haloes 1 and 16 both exhibit mergers during assembly but not jumps in L.
Note that the angular momenta shown in Figure 8 are likely not taken up by the SMS at the center of the disk, which we do not resolve in our models.Previous studies going to higher resolution for shorter times report that barswithin-bars instabilities arise on scales that are intermediate to those of the star and center of the disk that rapidly transport angular momentum outward (e.g., Wise et al. 2008).If such mechamisms did not arise SMSs could not form at the center of the disk because they can only accrete at most a few percent of the Keplerian angular momentum without being destroyed by the ΩΓ radiative instability (Haemmerlé et al. 2018b).

Disk Stability
The gravitational stability of disks is often parametrised by the radially-dependent Toomre Q parameter (Toomre 1964), where cs is the sound speed, v turb is the turbulent speed, σ is the surface density, and κ is the epicyclic frequency (Binney & Tremaine 1987;Oh & Haiman 2002), If the rotation of the disk is Keplerian, κ reduces to the rotational frequency.The Q parameter represents the ability of restorative forces in the disk (rotation and pressure gradients) to counteract gravitational perturbations.If Q falls below unity, perturbations in the disk grow rather than dampen, and it is prone to fragmentation at that radius.Radial profiles of the Toomre parameter for the disks in haloes 10, 12, and 1 are shown in Figure 9 for the times in Figure 3.At 0.226 kyr in halo 10, Q is greater than 1 at all radii and, as shown in Figure 3, the disk is stable.At 1.029 Myr Q dips below 1 at ∼ 0.25 pc and 0.75 pc and the disk becomes prone to gravitiational instabilities that distort the spiral arms near their base and create turbulence but do not break them up or produce long-lived clumps at this time.However, the instabilities persist and later cause the disk to fragment at 1.14 Myr and produce the clumps visible at 1.261 Myr.Q falls below 1 over a greater range of radii by this time because of the growth of the disk in mass and surface area but it again exceeds 1 everywhere by 2.006 Myr, when the fragments have coalesced into a satellite disk whose orbit around the main disk has stabilised it against fragmentation.
In contrast to halo 10, the disk in halo 12 exhibits multiple episodes of fragmentation.Q falls below 1 at several radii at 1.174 Myr and again at 1.851 Myr in halo 12 and clumps are visible at those radii in Figure 3. Outward transport of angular momentum due to the destruction of clumps by tidal torques from the disk stabilise the disk somewhat at 1.435 Myr, when tidal tails are visible instead of fragments, but Q has again fallen below 1 at 0.4 pc and beyond 0.9 pc at this time, portending the next round of fragmentation at 1.851 Myr.As in halo 10, as the disk in halo 12 grows in mass more of it becomes subject to gravitational instabilities.
The Toomre parameter falls below 1 by 0.289 Myr at radii of 0.5 pc -1.1 pc in halo 1, when the disk has indeed produced the two massive clumps in Figure 3.The disk continues to be unstable at these radii at 0.540 Myr but then becomes stable at 1.261 Myr after the clumps have merged to produce a satellite disk who mass rivals that of the main disk.Their mutual orbit suppresses fragmentation in both disks.The satellite then merges with the main disk by 0.995 Myr, which ejects the tidal tails in Figure 3, and it again becomes prone to fragmentation as Q again falls below 1 at ∼ 0.4 pc and beyond 1.1 pc.

Disk Rotation / Turbulence
Rotational velocities for the main disks in all eight haloes at the end of the simulations are shown in the left panel of Figure 10.After multiple dynamical times the disks have all settled into rotationally supported structures with mostly Keplerian velocities.Rotational velocities are greatest at 0.1 -0.2 pc and reach 75 -110 km s −1 .The highest rotation speeds occur in the disks with the largest accretion rates, such as halo 8, while disks exhibiting the the lowest rotation rates have the lowest accretion rates, such as haloes 10 and 16.These velocities are consistent with those in Regan & Haehnelt (2009b), who report rotation rates of up to 60 km s −1 , with higher velocities corresponding to more massive disks.Their velocities are somewhat lower than ours because their disks are evolved for much shorter times, when they have only grown to 0.3 -0.6 pc in radius.
The peaks and dips in velocity visible in some of the haloes at the end of the run are due to satellite disks that are orbiting the main disk at faster or slower rates.The absence of these features in the central 0.3 -0.5 pc, along with our Toomre profiles, indicate that fragmentation does not occur near the centers of the disks at later times.The bump at ∼ 200 pc in halo 2 is due to an interloper that later merges with the halo.The widths of the peaks in velocity mark the outer boundaries of the disks, as gas at larger radii falling onto the disks has much lower rotation rates.
We calculate turbulent velocities in the disks by subtracting the infall and rotational velocities from the absolute velocities in quadrature: r .These profiles are shown for all eight disks at the end of the simulation in the right panel of Figure 10.They peak at slightly larger radii than the rotational velocities, 0.2 -0.3 pc, as gas from the spiral arms crashes into the center of the disk and the kinetic energy of bulk inflow builds up turbulent cascades and the gas becomes chaotic.Turbulent motions transport angular momentum outward through shocks that dissipate kinetic energy and they help support the disks against further collapse.

Disk Evolution vs. Halo Assembly History
Although we find a variety of disk evolution ranging from relatively stable disks to violently fragmenting ones, for the most part there is no apparent correlation between their evolution and the assembly histories or spin parameters of our eight haloes.For example, we find stable disks in haloes with high and low spin parameters and with as few as one and as many as three major mergers in their assembly histories (haloes 10 and 16, respectively).Likewise, disks in haloes 1 and 20 both fragment at early times but halo 1 has twice the spin parameter of halo 20.Furthermore, there is no obvious connection between longevity of binary disks and halo spin.
Halo 1 has a much longer-lived satellite disk than halo 2 but a similar spin parameter.Baryons become gravitationally decoupled from DM in atomically-cooled haloes on scales of 5 -10 pc so the evolution of the disk is governed by local gas flows and turbulence rather than DM dynamics on larger scales.DM dynamics does appear to affect the evolution of the disks in haloes 10 and 16.From Figure 7 it can be seen that they are only half as massive as those in the other haloes and have correspondingly lower central accretion rates.These haloes begin to atomically cool at z = 17.2 and 16.5, just after major mergers at z = 18 and 17.5 that give them the highest spin parameters of all the haloes.The gas inherits some of this angular momentum as it collapses to form the disks, which have more rotational support due to the larger centrifugal forces on the gas.The disks consequently grow more slowly and accumulate less mass at their centers.Suazo et al. (2019) studied fragmentation in three primordial haloes that grew primarily by accretion, major mergers, or some mix of the two in a variety of LW backgrounds.They found no fragmentation in the highest backgrounds, which approximated atomic cooling, but only evolved these haloes for a few hundred kyr.Our results are consistent with their models because fragmentation generally happened after these times in our disks.Binary and multiple SMSs have previously been shown to form in atomically-cooled disks at high redshift, but only in the special case of haloes with high spin parameters (Latif et al. 2020).Our new simulations prove that multiple SMS formation was the rule rather than the exception in atomically-cooled haloes because it occurred over their full range of spin parameters and for a variety of assembly histories.Fragmentation in the vicinity of the star on AU scales has also recently been used to invoke the formation of multiple massive stars instead of a single SMS (e.g., Regan & Downes 2018a).However, they could only follow the evolution of the clumps for at most a few hundred kyr, so it is not clear if they grow to large enough masses to become stars or are just taken up into the central object, as are most of the fragments in our models on larger scales.Regardless of the final fates of fragments on small scales, our simulations show that binaries or small multiples will eventually form in the halo on larger scales.

DISCUSSION AND CONCLUSION
Because DCBHs in our models typically end up in close proximity to each other after formation (0.1 -0.2 pc) and in high ambient densities ( 10 7 cm −3 ), radiative drag forces would rapidly decay their orbits and produce mergers whose GW emission could be detected by future observatories such as LISA (Hartwig et al. 2018) in addition to NIR or radio from the BHs themselves (Pacucci et al. 2015(Pacucci et al. , 2017;;Natarajan et al. 2017;Barrow et al. 2018;Whalen et al. 2020aWhalen et al. ,b, 2021;;Vikaeus et al. 2022) or their progenitor stars (Surace et al. 2018(Surace et al. , 2019)).If a SMS assumes a highly elliptical or-bit around a DCBH companion it could be destroyed in a tidal disruption event (TDE).Luminosities for TDEs of 10 -40 M Pop III stars with DCBHs have been calculated (Kashiyama & Inayoshi 2016) but not those for an SMS, which could produce an extended afterglow because of its large mass, so it is not known if such events would produce enough NIR flux to be detected today.Both types of TDEs may be possible in a given halo because X-rays from the DCBH are known to trigger lower-mass Pop III star formation in its vicinity at later times (Machacek et al. 2003;Aykutalp et al. 2014Aykutalp et al. , 2020)).
Using the accretion rates in Figure 6 in the Kepler stellar evolution code (Weaver et al. 1978;Woosley et al. 2002), we have determined that the SMSs in our disks would all evolve as cool, red supergiants with little ionizing UV that could affect their growth rates (Woods et al. 2021a,b).Consequently, radiation transport was never required in our Enzo models to obtain their true final masses.We do not include magnetic fields in our simulations, which are thought to suppress fragmentation on small scales and result in fewer, more massive objects (Turk et al. 2012;Latif et al. 2013bLatif et al. , 2014a;;Sharda et al. 2020).However, magnetic dynamos typically affect gas flows on much smaller spatial scales than those of clump formation in our models so their absence in our simulations probably did not affect fragmentation.Woods et al. (2021a) have now shown that SMS -SMS, SMS -DCBH, and DCBH -DCBH binaries could form in our haloes.Supermassive X-ray binaries with unique spectra could result from SMS -DCBH systems if they are in close enough proximity.Because we do not resolve flows close to the surface of the star, our accretion rates should be taken to be upper limits.Nevertheless, as noted in the Introduction, numerical simulations that follow the collapse of atomicallycooled flows down to protostellar scales for short times find that accretion proceeds at rates similar to ours even down to these radii because of efficient angular momentum transport by bar instabilities.
We only considered isothermal collapse, in which no H2 is present.In reality, it is difficult to destroy all H2 in the cores of massive haloes because of self shielding, even with very high LW backgrounds, and low to intermediate LW backgrounds were far more common in the early Universe.Studies have shown that even small amounts of H2 can dramatically reduce central infall rates that only produce 3000 -30,000 M stars, not 100,000 M stars (see Latif et al. 2021).There may thus have been two or three tiers of DCBH mass that led to populations of less massive quasars at high redshifts, not just the 10 9 M SMBHs at z > 6. H2 cooling also promotes fragmentation of the gas that could create even less massive stars.Given that primordial haloes fragment even in the absence of H2 cooling, it was likely a common occurence in all primordial haloes.We are investigating longterm disk evolution in a grid of LW backgrounds in all eight of our haloes now.

Figure 2 .
Figure 2. Spherically-averaged density (left) and temperature (right) profiles of the eight haloes at the onset of atomic cooling and rapid baryon collapse.

Figure 3 .
Figure 3. Projections of gas density in 3 of the 8 disks in our simulations.Top row: the relatively stable disk in halo 10.Left to right: 226 kyr, 1.02 Myr, 1.26 Myr, and 2.01 Myr.Center row: the more turbulent, rapidly fragmenting disk in halo 12. Left to right: 51 kyr, 1.17 Myr, 1.43 Myr and 1.85 Myr.Bottom row: early binary disk formation in halo 1. Left to right: 289 kyr, 540 kyr, 879 kyr and 995 kyr.

Figure 4 .
Figure 4. Binary and multiple disk systems in our runs.Left: the main disk and its satellite in halo 20 at 1.19 Myr.Right: Two smaller disks orbiting the main disk at 1.29 Myr in halo 2.

Figure 5 .
Figure 5. Center-to-center separations between the main disk and satellite disk in haloes 8 and 20.

Figure 6 .Figure 7 .
Figure6.Accretion rates at the centers of the main and binary disks in all eight haloes.

Figure 8 .
Figure 8.The evolution of angular momentum in the central 0.136 pc of the disks over time.

Table 1 .
Halo properties at the onset of atomic cooling.From left to right: halo, collapse redshift, mass at the onset of atomic cooling, spin parameter and number of major mergers.