# Turbulent cold flows gave birth to the first quasars

Latif  M. A.1 Whalen  D. J.∗2 Khochfar  S.3 Herrington  N. P.4 & Woods  T. E.5

###### Abstract

How quasars powered by supermassive black holes (SMBHs) formed less than a
billion years after the Big Bang is still one of the outstanding problems in
astrophysics 20 years after their discovery[1, 2, 3, 4]. Cosmological
simulations suggest that rare cold flows converging on primordial haloes in
low-shear environments could have created these quasars if they were 104 \-
105 solar masses at birth but could not resolve their formation[5, 6, 7, 8].
Semianalytical studies of the progenitor halo of a primordial quasar found
that it favours the formation of such seeds but could not verify if one
actually appeared[9]. Here we show that a halo at the rare convergence of
strong, cold accretion flows creates massive BH seeds without the need for UV
backgrounds, supersonic streaming motions, or even atomic cooling. Cold flows
drive violent, supersonic turbulence in the halo that prevents star formation
until it reaches a mass that triggers sudden, catastrophic baryon collapse
that forms 31,000 and 40,000 solar-mass stars. This simple, robust process
ensures that haloes capable of forming quasars by z $>$ 6 produce massive
seeds. The first quasars were thus a natural consequence of structure
formation in cold dark matter cosmologies, not exotic, finely-tuned
environments as previously thought[10, 11, 12, 13, 14].

Physics Department, College of Science, United Arab Emirates University, PO
Box 15551, Al-Ain, UAE<EMAIL_ADDRESS>
Institute of Cosmology and Gravitation, Portsmouth University, Dennis Sciama
Building, Portsmouth PO1 3FX<EMAIL_ADDRESS>
Institute for Astronomy, University of Edinburgh, Royal Observatory, Blackford
Hill, Edinburgh EH9 3HJ, UK

School of Physics and Astronomy, University of Exeter, Stocker Road, Exeter,
EX4 4QL, UK

National Research Council of Canada, Herzberg Astronomy & Astrophysics
Research Centre, 5071 West Saanich Road, Victoria, BC V9E 2E7, Canada

We investigated SMBH seed formation during the earliest stages of collapse of
the progenitor halo of a z $>$ 6 quasar with the Enzo cosmology code. This
halo grows to 1.4 $\times$ 1012 $\mathrm{M}_{\odot}$ ($\mathrm{M}_{\odot}$,
solar mass) by z $=$ 6 mostly by cold flow accretion[15, 16, 17] rather than
by major mergers with other gas-rich haloes[18], and it is not in close
proximity to any external sources of Lyman-Werner (LW) UV flux or in any LW
backgrounds when it first collapses. Our simulation resolves the halo when it
reaches a mass of 3.9 $\times$ 105 $\mathrm{M}_{\odot}$ at z $=$ 35, the
approximate minimum mass at which it could form Population III (Pop III)
stars. However, highly supersonic turbulence driven by the cold flows from
large scales prevent gas in the halo from collapsing into stars at this mass.
As shown at z $=$ 29 in Figure 1, turbulent velocities are 30 - 40 km s-1, a
factor of a few higher than the infall velocities and more than an order of
magnitude greater than the sound speed of the gas.

Over the next 45 Myr the halo grows by two orders of magnitude in mass, to 4
$\times$ 107 $\mathrm{M}_{\odot}$ by z $=$ 25.7 at an average rate of $\sim$
0.9 $\mathrm{M}_{\odot}$ yr-1. Although this is roughly four times the mass at
which atomic cooling normally begins[19] ($\rm 2.2\times
10^{7}~{}{\mathrm{M}_{\odot}}~{}[(1+z)/16]^{-3/2}$), the gas does not collapse
by Ly$\alpha$ emission because temperatures in the halo never exceed 2,000 K.
Collapse is instead triggered when infall velocities finally exceed turbulent
velocities because the halo has grown so much in mass, as shown in Figure 1,
and gravitational compression overcomes turbulent pressure support. We show
density and temperature projections of the halo at the onset of collapse in
Figure 2. The four cold accretion streams seen in the 100 pc image carry gas
at densities of 103 \- 104 cm-3 deep into the halo and create the highly
turbulent core visible in the 10 pc panel. Energy transfer from scales of tens
of kpc down to pc scales drives this turbulence. A dense clump appears in the
core 5 kyr after the onset of collapse, as shown in the 0.1 pc image. It forms
a small self-gravitating disk with two prominent spiral arms and a mass of
about 1,000 $\mathrm{M}_{\odot}$, as shown in the 0.02 pc panel. Analysis
indicates that it is stable to fragmentation, and it rotates at nearly the
Keplerian velocity with influx rates of about 1 $\mathrm{M}_{\odot}$ yr-1.
However, before this clump can form a star it becomes part of a much larger
clump that forms a supermassive star later on (C2, as discussed below).

Unlike collapse in normal Pop III star-forming haloes, which creates quasi-
static settling flows[20], or in atomically-cooled haloes in which there is at
most transonic turbulence, collapse in this halo proceeds on two scales via
two, distinct physical processes. As shown in Figure 3, it begins on scales of
20 - 200 pc when cold flows finally overcome the turbulent pressure they
originally created and settle onto the core of the halo at influx rates of up
to 40 $\mathrm{M}_{\odot}$ yr-1. The pileup of gas at the outer layers of the
core is visible in the peaks in the density and infall velocity profiles at 20
pc. Collapse on these scales is driven mostly by dynamics. Infall at larger
scales compresses gas to densities of 10-16 g cm-3 at radii below 0.01 pc,
activating rapid H2 formation via three-body reactions that almost fully
molecularise the core by 11 kyr as shown in Extended Data Figure 1. Rapid
cooling by these large H2 mass fractions causes a second stage of collapse at
radii below 0.1 pc, with inflow rates of up to 5 $\mathrm{M}_{\odot}$ yr-1
even down to 20 AU, the smallest scales resolved by our simulation at this
stage. Densities at the center of the halo rise by four orders of magnitude to
$\sim$ 10-11 g cm-3 by 11 kyr. Gravitational compression briefly raises
temperatures at the center of the halo from 750 K to 1,200 K before rapid H2
cooling drops them back down to 700 K. The dramatic rise in density at the
center of the halo causes gas to pile up at radii of 0.01 pc, where a shock
forms and heats the gas to 1,700 K at 11 kyr (this shock is also visible in
the dip in infall rates at the same radius).

While resolutions of 20 AU capture the initial stages of collapse on all
relevant scales they are not practical for following massive seed formation
for 1 - 2 Myr because of short Courant times. We halted the simulation,
restarted it from a slightly higher redshift (z $=$ 25.73), when the halo was
only resolved by 18 levels of refinement (0.01 pc), and then evolved the halo
at this new limit of refinement for 1.6 Myr to determine if it produces a
supermassive star (SMS). As shown in Figure 4, turbulence prevents the
formation of any coherent structures in the halo like the large accretion
disks that form during the collapse of atomically-cooled haloes[21, 22]. The
core of the halo instead develops a wispy, filamentary structure that is
interspersed with dense clumps. The two most massive ones, C1 and C2, appear
in the first 100 kyr. C1 initially grows mainly through accretion but merges
with several other clumps 200 kyr later. As shown in Extended Data Figure 2,
these collisions raise its mass to 104 $\mathrm{M}_{\odot}$ over 400 kyr and
it then doubles in mass over the next 600 kyr. Its growth later falls off
because it migrates into regions with lower gas densities. C1 reaches a mass
of 2.6 $\times$ 104 $\mathrm{M}_{\odot}$ just before merging with another
cluster of clumps at the end of the run.

C2 grows to 1,000 $\mathrm{M}_{\odot}$ in 100 kyr before multiple encounters
with dense filaments increase its mass by a factor of four over the next 100
kyr. It then grows to 8,000 $\mathrm{M}_{\odot}$ in about 900 kyr primarily by
accretion. The four cold accretion streams impart angular momentum to the core
and create the transient annular structure shown in Fig. 4 that persists for
about a Myr. An elongated dense filament forms in the annulus and fragments
into multiple clumps. C2 merges with a few of these clumps and grows by a
factor of five in mass in about 200 kyr. The other clumps orbit C2, which
briefly forms spiral arms due to tidal interactions with them, but most of
them eventually merge with C2 except for one that is ejected from the cluster
and two that remain in its orbit. C2’s mass at the end of the simulation is
4.1 $\times$ 104 $\mathrm{M}_{\odot}$.

C1 and C2 collapse before any stars form in the halo, at infall rates of 0.01
- 0.55 $\mathrm{M}_{\odot}$ yr-1 (left panel of Extended Data Figure 3). These
rates are on par with those in atomically-cooled haloes that form SMSs in 1 -
2 Myr[23]. We modeled the evolution of the stars that form in these clumps
with the MESA code and show their Hertzsprung-Russell diagrams in the right
panel of Extended Data Figure 3. Both initially evolve along the Hayashi track
as cool, red hypergiants because of H- opacity in their outer layers, which
causes them to expand and remain at approximately 5,500 K as the stars grow in
mass. However, during periods of quiescence when there is little or no
accretion at later times, the stars can contract onto hotter and bluer tracks
and both eventually reach temperatures of 100,000 - 120,000 K. SMSs 1 and 2
collapse at masses of 31,316 $\mathrm{M}_{\odot}$ and 39,765
$\mathrm{M}_{\odot}$, respectively. Although both stars are millions of times
more luminous than the Sun, neither will be visible to the James Webb Space
Telescope[24] but their BHs could be detected at z $\gtrsim$ 20 after growing
to 105 $\mathrm{M}_{\odot}$[25]. Since both stars become blue and hot later in
their lives, their ionizing UV flux could reduce their accretion rates[26] so
these masses should be taken to be upper limits. Nevertheless, it is clear
that dynamically-driven collapse creates two DCBHs that are viable seeds for a
quasar forming by z $\sim$ 6 - 7.

Our simulations show that the number of quasars at z $>$ 6 does not depend on
LW backgrounds so it is not subject to the large uncertainties in the critical
flux, Jcrit, required to form DCBHs in atomically cooled haloes, which can
vary several orders of magnitude[27, 28]. The number of haloes that form
quasars by z $\sim$ 6 in large-scale numerical simulations[29] is 1 - 2
$\times$ 10-8 cMpc-3, within an order of magnitude of the number of quasars
observed at z $>$ 6, about 10-9 cMpc-3 or 1 cGpc-3. Delayed collapse due to
turbulence driven by cold flows therefore only has to occur in about 10% of
these haloes to account for the number of quasars at z $>$ 6 observed to date.
Because the topology of cold flows in the primordial Universe governed the
numbers of these unusual reservoirs, and they could form their own massive
seeds regardless of proximity to LW sources, the demographics of the first
quasars may ultimately have been determined by the structure of the cosmic
web.

## References

  * [1] Fan, X. _et al._ A Survey of $z>$ 5.7 Quasars in the Sloan Digital Sky Survey. II. Discovery of Three Additional Quasars at $z>$ 6\.  _AJ_ 125, 1649–1659 (2003).  arXiv:astro-ph/0301135.
  * [2] Mortlock, D. J. _et al._ A luminous quasar at a redshift of z = 7.085.  _Nature_ 474, 616–619 (2011).  1106.6088.
  * [3] Volonteri, M.  The Formation and Evolution of Massive Black Holes.  _Science_ 337, 544– (2012).  1208.1106.
  * [4] Wang, F. _et al._ A Luminous Quasar at Redshift 7.642.  _ApJ_ 907, L1 (2021).  2101.03179.
  * [5] Tenneti, A., Di Matteo, T., Croft, R., Garcia, T. & Feng, Y.  The descendants of the first quasars in the BlueTides simulation.  _MNRAS_ 474, 597–603 (2018).  1708.03373.
  * [6] Smidt, J., Whalen, D. J., Johnson, J. L., Surace, M. & Li, H.  Radiation Hydrodynamical Simulations of the First Quasars.  _ApJ_ 865, 126 (2018).  1703.00449.
  * [7] Huang, K.-W., Di Matteo, T., Bhowmick, A. K., Feng, Y. & Ma, C.-P.  BLUETIDES simulation: establishing black hole-galaxy relations at high redshift.  _MNRAS_ 478, 5063–5073 (2018).  1801.04951.
  * [8] Zhu, Q. _et al._ The Formation of the First Quasars. I. The Black Hole Seeds, Accretion and Feedback Models.  _arXiv e-prints_ arXiv:2012.01458 (2020).  2012.01458.
  * [9] Lupi, A., Haiman, Z. & Volonteri, M.  Forming massive seed black holes in high-redshift quasar host progenitors.  _MNRAS_ 503, 5046–5060 (2021).  2102.05051.
  * [10] Alexander, T. & Natarajan, P.  Rapid growth of seed black holes in the early universe by supra-exponential accretion.  _Science_ 345, 1330–1333 (2014).  1408.1718.
  * [11] Latif, M. A., Bovino, S., Grassi, T., Schleicher, D. R. G. & Spaans, M.  How realistic UV spectra and X-rays suppress the abundance of direct collapse black holes.  _MNRAS_ 446, 3163–3177 (2015).  1408.3061.
  * [12] Hirano, S., Hosokawa, T., Yoshida, N. & Kuiper, R.  Supersonic gas streams enhance the formation of massive black holes in the early universe.  _Science_ 357, 1375–1378 (2017).  1709.09863.
  * [13] Woods, T. E., Heger, A., Whalen, D. J., Haemmerlé, L. & Klessen, R. S.  On the Maximum Mass of Accreting Primordial Supermassive Stars.  _ApJ_ 842, L6 (2017).  1703.07480.
  * [14] Woods, T. E., Patrick, S., Elford, J. S., Whalen, D. J. & Heger, A.  On the Evolution of Supermassive Primordial Stars in Cosmological Flows.  _ApJ_ 915, 110 (2021).  2102.08963.
  * [15] Feng, Y., Di Matteo, T., Croft, R. & Khandai, N.  High-redshift supermassive black holes: accretion through cold flows.  _MNRAS_ 440, 1865–1879 (2014).  1312.1391.
  * [16] Lupi, A. _et al._ High-redshift quasars and their host galaxies - I. Kinematical and dynamical properties and their tracers.  _MNRAS_ 488, 4004–4022 (2019).  1901.02464.
  * [17] Valentini, M., Gallerani, S. & Ferrara, A.  Host galaxies of high-redshift quasars: SMBH growth and feedback.  _MNRAS_ 507, 1–26 (2021).  2107.05638.
  * [18] Li, Y. _et al._ Formation of z~6 Quasars from Hierarchical Galaxy Mergers.  _ApJ_ 665, 187–208 (2007).  astro-ph/0608190.
  * [19] Wise, J. H. _et al._ Formation of massive black holes in rapidly growing pre-galactic gas clouds.  _Nature_ 566, 85–88 (2019).  1901.07563.
  * [20] Bromm, V. & Larson, R. B.  The First Stars.  _ARA &A_ 42, 79–118 (2004).  arXiv:astro-ph/0311019.
  * [21] Lodato, G. & Natarajan, P.  Supermassive black hole formation during the assembly of pre-galactic discs.  _MNRAS_ 371, 1813–1823 (2006).  arXiv:astro-ph/0606159.
  * [22] Regan, J. A. & Haehnelt, M. G.  Pathways to massive black holes and compact star clusters in pre-galactic dark matter haloes with virial temperatures $>$ 10000K.  _MNRAS_ 396, 343–353 (2009).  0810.2802.
  * [23] Patrick, S. J., Whalen, D. J., Elford, J. S. & Latif, M. A.  The Collapse of Atomically-Cooled Primordial Haloes. I. High Lyman-Werner Backgrounds.  _arXiv e-prints_ arXiv:2012.11612 (2020).  2012.11612.
  * [24] Surace, M. _et al._ On the Detection of Supermassive Primordial Stars.  _ApJ_ 869, L39 (2018).  1811.08911.
  * [25] Whalen, D. J. _et al._ Finding Direct-collapse Black Holes at Birth.  _ApJ_ 897, L16 (2020).  2005.03018.
  * [26] Latif, M. A., Khochfar, S., Schleicher, D. & Whalen, D. J.  Radiation hydrodynamical simulations of the birth of intermediate-mass black holes in the first galaxies.  _MNRAS_ 508, 1756–1767 (2021).  2012.09177.
  * [27] Agarwal, B., Smith, B., Glover, S., Natarajan, P. & Khochfar, S.  New constraints on direct collapse black hole formation in the early Universe.  _MNRAS_ 459, 4209–4217 (2016).  1504.04042.
  * [28] Valiante, R., Agarwal, B., Habouzit, M. & Pezzulli, E.  On the Formation of the First Quasars.  _Publications of the Astronomical Society of Australia_ 34, e031 (2017).  1703.03808.
  * [29] Di Matteo, T., Croft, R. A. C., Feng, Y., Waters, D. & Wilkins, S.  The origin of the most massive black holes at high-z: BlueTides and the next quasar frontier.  _MNRAS_ 467, 4243–4251 (2017).  1606.08871.

Figure 1: Spherically-averaged velocity profiles of the halo. The infall
velocities (solid), turbulent velocities (dashed) and sound speeds (dotted)
are at just before collapse (z $=$ 29, green), the onset of collapse (z $=$
25.74, blue), and well into collapse (z $=$ 25.73, red). Enclosed gas masses
of 104 $\mathrm{M}_{\odot}$ and 106 $\mathrm{M}_{\odot}$ correspond to radii
of about 1 pc and 100 pc, respectively. Figure 2: Catastrophic baryon collapse
in the halo at z $=$ 25.7. Densities are in a - d and temperatures are in e -
h where a and e are at the onset of collapse, b and f are at 4.5 kyr, c and g
are at 5.5 kyr, and d and h are at 11 kyr. Figure 3: Spherically-averaged gas
profiles of the halo. The enclosed masses (a), temperatures (b), densities
(c), and inflow rates (d) are for the onset of collapse (green), 5 kyr (blue),
and 11 kyr (red). Figure 4: Evolution of the turbulent core. C1 and C2 are
marked by the white circles and the images are 25 pc on a side. The times in a
- d are 0.6 Myr - 1.7 Myr after collapse.

Methods

Enzo[2] is an adaptive mesh refinement (AMR) cosmology code with an $N-$body
particle-mesh scheme for evolving dark matter (DM)[3, 4] that is self-
consistently coupled to hydrodynamics, nonequilibrium primordial gas
chemistry[5], and a multigrid Poisson solver for self-gravity. Our simulations
utilise the piecewise parabolic method for hydrodynamics[6, 7] and the HLLC
scheme for enhanced stability with strong shocks and rarefaction waves[8]. The
reaction network follows the thermal and chemical evolution of the gas by
evolving nine primordial species (H, He, e-, H+, He+, He2+, H-, H2 and
H${}_{2}^{+}$) with updates to the gas energy that include collisional
excitational and ionizational cooling by H and He, recombination cooling, H2
cooling, bremsstrahlung cooling, and inverse Compton cooling by the cosmic
microwave background (CMB). Our H2 cooling rates[9] are valid at low densities
and high densities, at which local thermodynamic equilibrium (LTE) effects
become important. They also include the transition from optically-thin to
optically-thick cooling at densities when H2 becomes opaque to its own cooling
lines[10]. Our chemistry model omits collisional ionizations of H by H and
contributions to the H2 collisional dissociation rate by dissociative
tunnelling, which are important when strong LW fields are present[11, 12].
However, since there are no LW fluxes in our simulation and the temperatures
at which these two reactions become significant are 8000 K and 2000 - 4500 K,
respectively (which are greater than those in our halo), their omission has
little impact on our results. The choice of chemistry network, reaction
coefficients and update scheme can lead to some variations in the thermal
evolution of the gas.

##  0.1 Halo selection and simulation setup.

We initialise our simulations at z $=$ 150 with cosmological initial
conditions generated with MUSIC[13] using the second-year Planck[14] best fit
lowP + lensing + BAO + JLA + H0 cosmological parameters:
$\Omega_{\mathrm{M}}=$ 0.3089, $\Omega_{\Lambda}=$ 0.691,
$\Omega_{\mathrm{b}}=$ 0.0486, h $=$ 0.677, $\sigma_{8}=$ 0.816, and n $=$
0.967. The simulation was performed in a 37.5 Mpc box. Given that there are
only a few dozen haloes per Gpc3 at the nexus of converging cold flows that
sustain rapid early quasar growth in large-scale cosmological simulations, no
single 37.5 Mpc box at random would be expected to enclose one[15]. We
therefore evolved dozens of such boxes[16] with different random seeds at a
low top-grid resolution of $\rm 256^{3}$ and DM mass resolution of 108
$\mathrm{M}_{\odot}$ until we found a halo that reached 1.4 $\times$ 1012
$\mathrm{M}_{\odot}$ by z $=$ 6 primarily by accretion rather than major
mergers (i.e., collisions with other haloes with masses greater than 10% of
the mass of the host halo).

To capture the assembly history of this halo in greater detail the simulation
was rerun with a DM particle mass of 2.3 $\times$ 105 $\mathrm{M}_{\odot}$ and
outputs at 10 Myr intervals. A halo catalogue for each output was then
computed with Rockstar[17] and the merger history of the halo was constructed
from the catalogues with the consistent-tree algorithm. The mass of the halo
over time is shown in Extended Data Figure 4. Analysis shows that the halo
grows mainly by accretion and minor mergers with mass ratios below 1:10,
having experienced only one major merger with a mass ratio of 1:3 at z $=$
11.3 and two mergers with mass ratios of 1:10 at z $=$ 12.

We then identified the DM particles in the Lagrangian volume of the halo at z
$=$ 6 with Rockstar, traced their masses and positions back to the
initialization redshift, added one level of refinement in the region that
encloses these particles (twice the Lagrange radius of the halo), and reran
the simulation down to z $=$ 6\. We again identified the DM particles of the
now better-resolved halo, traced their masses and positions back to z $=$ 150,
applied another level of refinement to the volume enclosing this new particle
set, and evolved the simulation to the final redshift. This procedure[18] was
repeated until it produced five nested grids with a maximum DM resolution of
3,636 $\mathrm{M}_{\odot}$, an effective initial spatial resolution of 8,1923,
and 975 million DM particles. The minimum DM particle mass allows us to
resolve the halo at the lowest mass at which it could have formed normal Pop
III stars, a few 105 $\mathrm{M}_{\odot}$.

Our simulation is partitioned into two stages. In the first, we evolved the
halo from z = 150 with DM and gas down to the onset of collapse at z = 25.7
while allowing up to 25 additional levels of refinement in the innermost
nested grid for a maximum resolution of 20 AU. Our refinement criteria are
based on baryonic overdensity, Jeans refinement and the DM mass resolution. A
cell is flagged for refinement if it exceeds 2.5 times the cosmic mean density
or a DM particle density of 0.0391 $\rho_{\rm{DM}}r^{\ell\alpha}$, where
$\rho_{\rm{DM}}$ is the dark matter density, r $=$ 2 is the refinement factor,
$\ell$ is the refinement level, and $\alpha=-0.2$, which enforces super-
Lagrangian refinement. The Jeans length in the simulations was always covered
by at least 64 cells in order to resolve turbulent eddies in the flows[19,
20]. We also used the must-refine particle flag, which guarantees that DM
particles are refined to at least their generation level.

When the simulation reached 25 levels of refinement, we halted it and imposed
a pressure floor on the highest level to prevent collapse to even smaller
scales. This was done by artificially setting temperatures in the cells to
ensure they are Jeans stable and cannot fragment. We then evolved the halo
with this floor for 20 kyr to follow its collapse at a resolution of 20 AU. In
the second stage of the run we restarted the simulation from a slightly
earlier time (z $=$ 25.73), when it was only at a maximum of 18 levels of
refinement (a resolution of $\sim$ 0.01 pc), and limited the AMR to 18 levels
thereafter to be able to evolve the halo for the times required for the two
SMSs to form and collapse, about 1.7 Myr. We set a new pressure floor on the
18th level so gas would not collapse below our resolution limit.

##  0.2 H2 Cooling.

H2 cooling can dissipate virial shock heating and cool a halo of this mass to
less than 2,000 K because its H2 mass fractions remain high in the absence of
nearby sources of LW photons, and at 1,000 - 2,000 K the H2 cooling rates can
be 10 - 100 times higher than those in normal Pop III star-forming haloes.
Densities at the virial radius of this halo are also several times those of
other haloes of comparable mass and also enhance its cooling. A few cells
reach 8,000 - 10,000 K and cool via Ly$\alpha$ emission but their total mass
is small and they have no effect on the collapse. As shown in Extended Data
Figure 1, H2 mass fractions vary from 0.1 - 0.8 at radii less than 0.01 pc to
a few 10-3 out to the boundary of the turbulent core at 10 - 20 pc and then
down to 10-4 at 100 pc. H2 cooling generally keeps temperatures in the halo
below 1,000 K, the 2,000 K peak at 100 pc is due to virial heating by the
accretion shock. The densities exhibit significant departures from classic
$r^{-2}$ isothermal collapse profiles because temperatures vary by an order of
magnitude in the halo and collapse proceeds in two stages at two different
scales. The gas enclosed by the halo is 4 $\times$ 106 $\mathrm{M}_{\odot}$,
two orders of magnitude higher than in normal Pop III star-forming haloes and
a factor of a few higher than conventional direct-collapse black hole (DCBH)
host haloes.

##  0.3 Pop III Star Formation.

Although we did not have star or sink particle formation turned on, no normal
Pop III stars formed in the halo during the first stage of the simulation
because (extremely) large turbulent velocities prevented the collapse of the
gas, even though it was cooled by H2. In molecular clouds today, supersonic
turbulent shocks can create density fluctuations that can collapse to stars
via efficient cooling by dust, metals and molecules[21, 22]. Such shocks
produce similar structures in our simulations but they do not form stars
because H2 cooling is much less efficient. Gas in the halo did not reach high
enough densities to trigger Pop III star formation before global collapse in
the halo took place.

No Pop III stars formed during the collapse of the halo in the second stage of
the run before the two supermassive stars appeared because supersonic
turbulence driven by the cold flows modified the Jeans mass of the gas:

$M_{\rm Jeans}\,=\,\frac{\pi}{6}\frac{v_{\rm eff}^{3}}{G^{3/2}\rho^{1/2}},$
(1)

where

$v_{\rm eff}\,=\,\sqrt{c_{\rm s}^{2}+v_{\rm t}^{2}}.$ (2)

Here, vt is the turbulent velocity[23] and cs is the sound speed. After the
onset of collapse, densities in the core rise to over 1 $\times$ 10-16 g cm-3
and vt and cs are $\sim$ 20 km s-1 and 2 km s-1, respectively, yielding an
effective Jeans mass of about 104 $\mathrm{M}_{\odot}$. Consequently, in spite
of H2 cooling, the gas could only fracture on large mass scales and no normal
Pop III stars formed, just the two massive clumps that created the 31,000
$\mathrm{M}_{\odot}$ and 39,000 $\mathrm{M}_{\odot}$ stars. No star or sink
formation was turned on during the second stage of the simulation either, so
to identify potential sites of SMS formation we used the clump finder in the
yt analysis package[24] to search for objects that are gravitationally bound,
are resolved by a minimum of 20 cells, and have threshold densities $\geq$
10-17 g cm-3. C1 and C2 were the first self-gravitating clumps to form SMSs
after the halo began to collapse, no other clumps produced stars before the
SMSs became DCBHs.

The pressure floors we imposed on the highest level of refinement at the end
of stage 1 and during stage 2 suppress fragmentation on scales of at most 10
times their resolution limits (200 AU and 0.1 pc, respectively[25]).
Fragmentation can occur in accretion disks on scales below our resolution[26,
27, 28], but these low-mass clumps are expected to migrate inwards and merge
with the central object on timescales that are much shorter than its Kelvin-
Helmholz contraction times onto the main sequence. This has been found to be
true in 3D simulations[29] and analytical studies of Pop III disk
fragmentation[30, 31].

We confirmed that no Pop III stars formed in the halo by turning on sink
particles instead of pressure floors in stage 1 and for the first Myr of stage
2. Sinks were created if five criteria were met:

  1. 1.

The grid cell is at the maximum refinement level.

  2. 2.

The gas density in the cell is higher than the Jeans density.

  3. 3.

The flow is convergent ($\nabla\cdot\mathbf{v}<$ 0).

  4. 4.

The cooling time is shorter than the free-fall time.

  5. 5.

The cell is at a local minimum of the gravitational potential.

Besides accreting gas, these sinks, once created, also act as star particles
so they can be sources of ionizing UV radiation. This algorithm has been
tested and used in a number of cosmological simulations[32, 33, 26, 34] and is
similar to those used in studies of star formation in molecular clouds
today[35, 36]. No sink particles formed during either stage of the run,
confirming that no Pop III stars formed before the two SMSs. This outcome is
consistent with previous work that found essentially no differences when using
sink particles or pressure floors in simulations of DCBH formation[37]. We
show images of the halo 0.6 Myr and 1.0 Myr after the onset of collapse in our
sink particle run in Extended Data Figure 5. They exhibit morphologies that
are quite similar to those in Figure 4. As shown in Extended Data Figure 2, C1
and C2 have similar masses in the two runs and were the first to collapse and
form SMSs in them. They also exhibit central infall rates that are similar to
those in the original run, showing that pressure floors had little effect on
the evolution of the core of the halo.

Our results are also consistent with a recent simulation in which supersonic
gas streaming motions relative to the DM potential of a primordial halo[38]
modify the Jeans mass of the gas and prevent star formation in the halo until
it grows to over 107 $\mathrm{M}_{\odot}$ and creates a SMS[12]. Turbulence
due to 3$\sigma$ fluctuation streaming velocities suppressed collapse until
the halo grew to 2 $\times$ 107 $\mathrm{M}_{\odot}$ and formed the star via
H2 and atomic cooling. As in our simulation, no Pop III stars formed before
the halo reached this mass. Although the origin of the supersonic turbulence
in our study is quite different, its suppression of star formation prior to
collapse is similar.

Pop III stars are expected to eventually appear in the halo because X-rays
from the DCBH can catalyze H2 formation and trigger their formation[39, 40,
41], and they play an important role in regulating the growth of the BH at
later times[6]. But even if normal Pop III stars had formed before the two
SMSs they probably would not have been able to interfere with their growth.
Recent radiation hydrodynamical simulations indicate that such stars would
have had to be at least 1,000 $\mathrm{M}_{\odot}$ before their ionizing UV
flux could disrupt massive inflows in their vicinity in similar
environments[26].

H2 cooling has previously been shown to play a role in SMS formation when
large reservoirs of gas and high halo temperatures are present. Only the most
extreme LW backgrounds destroy all the H2 in a halo because it can self-shield
in the high densities at the center. Not even energy injection by violent
mergers destroy it all, even in moderate LW backgrounds, and such haloes have
still been found to form SMSs[19, 33]. A 104 $\mathrm{M}_{\odot}$ star even
formed in the presence of metals in one simulation[42], in which Z $<$ 0.001
$\mathrm{Z}_{\odot}$ ($\mathrm{Z}_{\odot}$, solar metallicity). These studies
show that any rapid dynamical collapse can produce SMSs, regardless of
origin[43, 32].

##  0.4 Jeans masses and H2 mass fractions at the onset of collapse.

In Extended Data Figure 1 we plot the ratio of the enclosed gas mass to the
Jeans mass, $\rm M_{en}/M_{J}$, at the onset of dynamically-driven collapse,
at 5 kyr, and at 11 kyr. Wherever this ratio exceeds one the gas is prone to
collapse. The three plots show that the gas in the halo becomes increasingly
unstable to collapse over time, with the Jeans instability operating at the
smallest radii resolved in the simulation by just 11 kyr. The peaks in the
ratio indicate that gravitational instabilities are present at multiple scales
and that almost all the gas in the halo is destined to collapse. This trend is
also visible in the infall rates, $\rm 4\pi R^{2}\rho v$, plotted as a
function of enclosed mass at the same times in Extended Data Figure 1. At 5
kyr collapse rates at the center of the halo are modest, a few 0.01
$\mathrm{M}_{\odot}$ yr-1, but rise to nearly 1 $\mathrm{M}_{\odot}$ yr-1 by
11 kyr.

##  0.5 Dynamical heating.

Recent studies have examined the formation of SMSs in atomically-cooled haloes
in marginal LW backgrounds via dynamical heating due to rapid growth driven by
major mergers[19, 33]. In the first study (W19), a halo reaches 2 $\times$ 107
$\mathrm{M}_{\odot}$ by $z\sim$ 16.5, having grown by a factor of 30 in mass
in just 30 Myr from $z=$ 21 - 19 because of a major merger. It was exposed to
a LW flux of 3 $J_{21}$ and collapsed at a mass of $\sim$ 2 $\times$ 106
$\mathrm{M}_{\odot}$. Dynamical heating of the gas in the halo due to rapid
inflows suppressed normal Pop III star formation even though the LW background
was small.

Although our halo also hosts rapid infall, dynamical heating plays little or
no role in the prevention of normal Pop III SF in it prior to catastrophic
collapse. While its average growth rate is larger than in W19, below 3
$\times$ 106 $\mathrm{M}_{\odot}$ its instantaneous growth rate is a few 105
$\mathrm{M}_{\odot}$ per unit redshift, an order of magnitude smaller than in
W19. The maximum growth rate of our halo is still a factor of 2 - 3 times
below that in W19 and peaks when the halo has already reached 107
$\mathrm{M}_{\odot}$. Typical H2 mass fractions in W19 are 10-5 \- 10-6, more
than an order of magnitude lower than those in our halo, $\sim$ 10-4. At these
H2 fractions the growth rate of our halo would have to exceed those in W19 by
an order of magnitude for dynamical heating to be effective (see their
Extended Data Figure 1) so it has little or no impact on collapse. We also
note that turbulent and radial infall velocities in W19 are about a factor of
ten lower than in our halo so turbulence has little effect on the collapse of
their halos. The conditions leading to the formation of SMSs in our halo are
thus quite distinct from those in W19.

##  0.6 SMS evolution and DCBH birth.

We use the Modules for Experiments in Stellar Astrophysics (MESA) code[44] to
follow the growth of the SMSs in the halo. MESA is a one-dimensional
lagrangian stellar evolution and hydrodynamics code with nuclear burning that
is implicitly coupled to updates of the equations of stellar structure[45,
46]. We initialise the stars as 10 $\mathrm{M}_{\odot}$, fully convective (n
$=$ 1.5) polytropes with hydrogen and helium mass fractions of 0.76 and 0.24,
respectively, and central densities $\rho_{\mathrm{c}}=$ 7.6 $\times$10-6 g
cm-3 and temperatures T${}_{\mathrm{c}}=$ 9.0 $\times$ 105 K.

MESA adaptively rezones the stars over their evolution. The mesh is refined
when gradients in pressure, temperature and 4He abundances between cells
exceed preset limits: $\delta\rm\log P/P>$ 1/30, $\delta\rm\log T/T>$ 1/80,
and $\delta\log(\chi+0.01)>$ 1/20 $\log\chi$, where $\chi$ is the 4He mass
fraction. This approach typically partitions the stars into 750 - 1,200 mass
zones, with larger numbers of zones being allocated to the center of the star
to resolve nuclear burning and convective processes. We use the 21-isotope
APPROX21 nuclear reaction network, which includes the p-p chain, CNO cycle,
helium burning, and the CO and alpha chains.

We apply the Ledoux criterion for convection and use Henyey mixing-length
theory[47]. The equation of state (EOS) in our models is a composite of
several datasets, including the OPAL/SCVH tables[48, 49], which are used at
lower densities and temperatures like the outer regions of the star and its
atmosphere, and the HELM and PC EOSs,[50, 51] which are used in the core of
the star. Above 104 $\mathrm{M}_{\odot}$ stars can encounter the
Chandrasekhar-Feynman general relativistic (GR) instability, in which thermal
photons in their cores reach densities at which they can enhance gravity via
the stress-energy tensor and cause the core to collapse.[52] This effect can
be approximated by the post-Newtonian Tolman-Oppenheimer-Volkoff (TOV)
correction to the equation of hydrostatic equilibrium,

$\frac{dP}{dr}=-\frac{Gm\rho}{r^{2}},\vspace{0.1in}$ (3)

by replacing the gravitational constant G with Grel,

$G_{\textrm{rel}}=G\left(1+\frac{P}{\rho c^{2}}+\frac{4\pi
Pr^{3}}{m_{r}c^{2}}\right)\left(1-\frac{2Gm_{r}}{rc^{2}}\right)^{-1}.$ (4)

The stars are evolved until they reach advanced stages of collapse, with
infall speeds above 1,000 km s-1.

Masses for C1 and C2 over time are shown in Extended Data Figure 2. Accretion
rates in the two clumps were estimated by dividing the change in mass of each
clump (as found by our halo finder) over 20 kyr intervals by this interval. At
times the rates can become negative if tidal interactions with other objects
strip mass from the clump, but these outflows occur on much larger scales than
accretion onto the star, which is not expected to lose mass during these
episodes. At most, accretion onto the star would temporarily halt so we set
the rates to zero when the two clumps lose mass. The rates are shown in the
left panel of Extended Data Figure 3, where t $=$ 0 marks the onset of the
collapse of the clumps. Most of the peaks are from turbulent flows but the
large dip at 1.45 Myr and subsequent spike at 1.55 Myr are due to tidal
interactions that cause C2 to exchange mass with another smaller clump. The
entropy of the accreted material is matched to that of the outer layer of the
star, which neglects the luminosity due to the accretion shock. However, it is
at most 104 $\mathrm{L}_{\odot}$ ($\mathrm{L}_{\odot}$, solar luminosity) even
at 1 $\mathrm{M}_{\odot}$ yr-1 so it is small in comparison to that of the
star. We ignore stellar rotation because the spiral arms of the disks that are
known to form in the clumps efficiently transport most of the angular momentum
from their centers[53, 54, 55]. There is no mass loss due to stellar winds
because of low internal opacities at zero metallicity[56, 57], and we neglect
pulsational mass losses because they are expected to saturate at accretion
rates that are an order of magnitude below those considered here[58].

Central p-p chain burning begins immediately in both stars but at very low
rates that cannot halt stellar contraction, as is the case in all massive
stars. They contract until they ignite low-level triple-alpha burning at
$\sim$ 750 yr and 1,400 yr in SMSs 1 and 2, respectively. Although the triple-
alpha process only provides a small contribution to the net energy generation
rate, it quickly produces the carbon, nitrogen, and oxygen required to
catalyze hydrogen burning via the CNO chain, which becomes the dominant energy
source in SMS 1 by 9,000 yr and 10,000 yr in SMS 2. CNO burning halts further
contraction in both stars and signals the onset of main sequence burning.

The radii of the two SMSs vary from 30 – 6,000 $\mathrm{R}_{\odot}$
($\mathrm{R}_{\odot}$, solar radii), or 0.14 – 27.9 AU, over their lifetimes
as seen in Extended Data Figure 3. Their internal structures are shown in the
Kippenhahn diagrams in Extended Data Figure 6. Most of the mass of the stars
lies within their convective interiors, relatively little mass is in their
higher-entropy outer envelopes. Modest mean accretion rates in the first Myr
allow SMS 2 to become largely thermally relaxed, remaining hot and blue even
with much stronger accretion after 1 Myr and becoming mostly convective after
1.2 Myr. Both stars die when they encounter the GR instability while still on
the hydrogen-burning main sequence. SMS 1 collapses with a central hydrogen
fraction of just below 20% while SMS 2 collapses earlier its evolution (after
the onset of very rapid accretion at 1 Myr) with a central hydrogen fraction
of $\sim$50%. Internal velocities for both stars are shown at the final stages
of collapse in Extended Data Figure 7.

Data Availability The Enzo parameter files and initial conditions files
generated by MUSIC that are required to perform the simulations are available
at https://doi.org/10.5281/zenodo.5853118 (DOI: 10.5281/zenodo.5853118). The
MUSIC input files required to generate the initial conditions are available at
https://sites.google.com/site/latifmaastro/ics. The MESA inlist files required
to evolve the two stars can be found at
https://cococubed.com/mesa_market/inlists.html.

Code Availability The code used to produce our cosmological simulations, Enzo
2.6, can be found at https://bitbucket.org/enzo/enzo-dev/tree/enzo-2.6.1. The
MESA code used to evolve the two stars, version 12778, can be found at
https://zenodo.org/record/3698354#.YgkaN-TfUlQ.

## References

  * [1]
  * [2] Bryan, G. L. _et al._ ENZO: An Adaptive Mesh Refinement Code for Astrophysics.  _ApJS_ 211, 19 (2014).
  * [3] Efstathiou, G., Davis, M., White, S. D. M. & Frenk, C. S.  Numerical techniques for large cosmological N-body simulations.  _ApJS_ 57, 241–260 (1985).
  * [4] Couchman, H. M. P.  Mesh-refined P3M - A fast adaptive N-body algorithm.  _ApJ_ 368, L23–L26 (1991).
  * [5] Anninos, P., Zhang, Y., Abel, T. & Norman, M. L.  Cosmological hydrodynamics with multi-species chemistry and nonequilibrium ionization and cooling.  _New Astronomy_ 2, 209–224 (1997).
  * [6] Woodward, P. & Colella, P.  The numerical simulation of two-dimensional fluid flow with strong shocks.  _Journal of Computational Physics_ 54, 115–173 (1984).
  * [7] Bryan, G. L., Norman, M. L., Stone, J. M., Cen, R. & Ostriker, J. P.  A piecewise parabolic method for cosmological hydrodynamics.  _Computer Physics Communications_ 89, 149–168 (1995).
  * [8] Toro, E. F., Spruce, M. & Speares, W.  Restoration of the contact surface in the HLL-Riemann solver.  _Shock Waves_ 4, 25–34 (1994).
  * [9] Glover, S. C. O. & Abel, T.  Uncertainties in H2 and HD chemistry and cooling and their role in early structure formation.  _MNRAS_ 388, 1627–1651 (2008).
  * [10] Ripamonti, E. & Abel, T.  Fragmentation and the formation of primordial protostars: the possible role of collision-induced emission.  _MNRAS_ 348, 1019–1034 (2004).
  * [11] Glover, S. C. O.  Simulating the formation of massive seed black holes in the early Universe - I. An improved chemical model.  _MNRAS_ 451, 2082–2096 (2015).
  * [12] Glover, S. C. O.  Simulating the formation of massive seed black holes in the early Universe - II. Impact of rate coefficient uncertainties.  _MNRAS_ 453, 2901–2918 (2015).
  * [13] Hahn, O. & Abel, T.  Multi-scale initial conditions for cosmological simulations.  _MNRAS_ 415, 2101–2121 (2011).
  * [14] Planck Collaboration _et al._ Planck 2015 results. XIII. Cosmological parameters.  _A &A_ 594, A13 (2016).
  * [15] Khandai, N., Feng, Y., DeGraf, C., Di Matteo, T. & Croft, R. A. C.  The formation of galaxies hosting z $\sim$ 6 quasars.  _MNRAS_ 423, 2397–2406 (2012).
  * [16] Trenti, M., Santos, M. R. & Stiavelli, M.  Where Can We Really Find the First Stars’ Remnants Today?  _ApJ_ 687, 1–6 (2008).
  * [17] Behroozi, P. S., Wechsler, R. H. & Wu, H.-Y.  The ROCKSTAR Phase-space Temporal Halo Finder and the Velocity Offsets of Cluster Cores.  _ApJ_ 762, 109 (2013).
  * [18] Wise, J. H.  Enzo-MRP-Music, https://github.com/jwise77/enzo-mrp-music.
  * [19] Federrath, C., Sur, S., Schleicher, D. R. G., Banerjee, R. & Klessen, R. S.  A New Jeans Resolution Criterion for (M)HD Simulations of Self-gravitating Gas: Application to Magnetic Field Amplification by Gravity-driven Turbulence.  _ApJ_ 731, 62 (2011).
  * [20] Latif, M. A., Schleicher, D. R. G., Schmidt, W. & Niemeyer, J.  Black hole formation in the early Universe.  _MNRAS_ 433, 1607–1618 (2013).
  * [21] Hennebelle, P. & Chabrier, G.  Analytical Theory for the Initial Mass Function: CO Clumps and Prestellar Cores.  _ApJ_ 684, 395–410 (2008).
  * [22] Federrath, C.  The turbulent formation of stars.  _Physics Today_ 71, 38–42 (2018).
  * [23] Stacy, A., Bromm, V. & Loeb, A.  Rotation speed of the first stars.  _MNRAS_ 413, 543–553 (2011).
  * [24] Turk, M. J. _et al._ yt: A Multi-code Analysis Toolkit for Astrophysical Simulation Data.  _ApJS_ 192, 9 (2011).
  * [25] Machacek, M. E., Bryan, G. L. & Abel, T.  Simulations of Pregalactic Structure Formation with Radiative Feedback.  _ApJ_ 548, 509–521 (2003).
  * [26] Greif, T. H. _et al._ Simulations on a Moving Mesh: The Clustered Formation of Population III Protostars.  _ApJ_ 737, 75 (2011).
  * [27] Becerra, F., Greif, T. H., Springel, V. & Hernquist, L. E.  Formation of massive protostars in atomic cooling haloes.  _MNRAS_ 446, 2380–2393 (2015).
  * [28] Becerra, F., Marinacci, F., Bromm, V. & Hernquist, L. E.  Assembly of supermassive black hole seeds.  _MNRAS_ 480, 5029–5045 (2018).
  * [29] Hosokawa, T. _et al._ Formation of Massive Primordial Stars: Intermittent UV Feedback with Episodic Mass Accretion.  _ApJ_ 824, 119 (2016).
  * [30] Latif, M. A. & Schleicher, D. R. G.  Magnetic fields in primordial accretion disks.  _A &A_ 585, A151 (2016).
  * [31] Inayoshi, K. & Haiman, Z.  Does disc fragmentation prevent the formation of supermassive stars in protogalaxies?  _MNRAS_ 445, 1549–1557 (2014).
  * [32] Regan, J. A. & Downes, T. P.  Rise of the first supermassive stars.  _MNRAS_ 478, 5037–5049 (2018).
  * [33] Regan, J. A. _et al._ The Formation of Very Massive Stars in Early Galaxies and Implications for Intermediate Mass Black Holes.  _The Open Journal of Astrophysics_ 3, 15 (2020).
  * [34] Latif, M. A., Whalen, D. & Khochfar, S.  The Birth Mass Function of Pop III Stars.  arXiv:2109.10655 (2021).
  * [35] Krumholz, M. R., McKee, C. F. & Klein, R. I.  Embedding Lagrangian Sink Particles in Eulerian Grids.  _ApJ_ 611, 399–412 (2004).
  * [36] Federrath, C., Banerjee, R., Clark, P. C. & Klessen, R. S.  Modeling Collapse and Accretion in Turbulent Gas Clouds: Implementation and Comparison of Sink Particles in AMR and SPH.  _ApJ_ 713, 269–290 (2010).
  * [37] Latif, M. A., Schleicher, D. R. G., Schmidt, W. & Niemeyer, J. C.  The characteristic black hole mass resulting from direct collapse in the early Universe.  _MNRAS_ 436, 2989–2996 (2013).
  * [38] Tseliakhovich, D. & Hirata, C.  Relative velocity of dark matter and baryonic fluids and the formation of the first structures.  _Phys. Rev. D_ 82, 083520 (2010).
  * [39] Machacek, M. E., Bryan, G. L. & Abel, T.  Effects of a soft X-ray background on structure formation at high redshift.  _MNRAS_ 338, 273–286 (2003).
  * [40] Aykutalp, A., Wise, J. H., Spaans, M. & Meijerink, R.  Songlines from Direct Collapse Seed Black Holes: Effects of X-Rays on Black Hole Growth and Stellar Populations.  _ApJ_ 797, 139 (2014).
  * [41] Aykutalp, A., Barrow, K. S. S., Wise, J. H. & Johnson, J. L.  Induced Metal-free Star Formation around a Massive Black Hole Seed.  _ApJ_ 898, L53 (2020).
  * [42] Chon, S. & Omukai, K.  Supermassive star formation via super competitive accretion in slightly metal-enriched clouds.  _MNRAS_ 494, 2851–2860 (2020).
  * [43] Latif, M. A. & Volonteri, M.  Assessing inflow rates in atomic cooling haloes: implications for direct collapse black holes.  _MNRAS_ 452, 1026–1044 (2015).
  * [44] Paxton, B. _et al._ Modules for Experiments in Stellar Astrophysics (MESA).  _ApJS_ 192, 3 (2011).
  * [45] Paxton, B. _et al._ Modules for Experiments in Stellar Astrophysics (MESA): Planets, Oscillations, Rotation, and Massive Stars.  _ApJS_ 208, 4 (2013).
  * [46] Paxton, B. _et al._ Modules for Experiments in Stellar Astrophysics (MESA): Convective Boundaries, Element Diffusion, and Massive Star Explosions.  _ApJS_ 234, 34 (2018).
  * [47] Henyey, L., Vardya, M. S. & Bodenheimer, P.  Studies in Stellar Evolution. III. The Calculation of Model Envelopes.  _ApJ_ 142, 841 (1965).
  * [48] Rogers, F. J. & Nayfonov, A.  Updated and Expanded OPAL Equation-of-State Tables: Implications for Helioseismology.  _ApJ_ 576, 1064–1074 (2002).
  * [49] Saumon, D., Chabrier, G. & van Horn, H. M.  An Equation of State for Low-Mass Stars and Giant Planets.  _ApJS_ 99, 713 (1995).
  * [50] Timmes, F. X. & Swesty, F. D.  The Accuracy, Consistency, and Speed of an Electron-Positron Equation of State Based on Table Interpolation of the Helmholtz Free Energy.  _ApJS_ 126, 501–516 (2000).
  * [51] Potekhin, A. Y. & Chabrier, G.  Thermodynamic Functions of Dense Plasmas: Analytic Approximations for Astrophysical Applications.  _Contributions to Plasma Physics_ 50, 82–87 (2010).
  * [52] Chandrasekhar, S.  The Dynamical Instability of Gaseous Masses Approaching the Schwarzschild Limit in General Relativity.  _ApJ_ 140, 417 (1964).
  * [53] Haemmerlé, L., Woods, T. E., Klessen, R. S., Heger, A. & Whalen, D. J.  On the Rotation of Supermassive Stars.  _ApJ_ 853, L3 (2018).
  * [54] Haemmerlé, L., Woods, T. E., Klessen, R. S., Heger, A. & Whalen, D. J.  The evolution of supermassive Population III stars.  _MNRAS_ 474, 2757–2773 (2018).
  * [55] Haemmerlé, L. & Meynet, G.  Magnetic braking of supermassive stars through winds.  _A &A_ 623, L7 (2019).
  * [56] Vink, J. S., de Koter, A. & Lamers, H. J. G. L. M.  Mass-loss predictions for O and B stars as a function of metallicity.  _A &A_ 369, 574–588 (2001).
  * [57] Baraffe, I., Heger, A. & Woosley, S. E.  On the Stability of Very Massive Primordial Stars.  _ApJ_ 550, 890–896 (2001).
  * [58] Hosokawa, T., Yorke, H. W., Inayoshi, K., Omukai, K. & Yoshida, N.  Formation of Primordial Supermassive Stars by Rapid Mass Accretion.  _ApJ_ 778, 178 (2013).

MAL was supported by UAEU UPAR grant No. 31S390. SK acknowledges support from
the UK STFC via grant ST/V000594/1. NPH acknowledges funding from the European
Research Council for the Horizon 2020 ERC Consolidator Grant project ICYBOB,
grant number 818940. TEW acknowledges support from the NRC-Canada Plaskett
Fellowship. The Enzo simulations were performed on HPC resources at UAEU and
the Institute of Cosmology and Gravitation at the University of Portsmouth
(Sciama).

DJW proposed this study, assisted in its development, and wrote the paper. MAL
developed this study, performed the Enzo simulations, and analyzed its data.
SK contributed to the development of this study and the interpretation of its
results. NPH performed the MESA calculations and analyzed its results with
TEW.

The authors declare no competing interests.

are available in the online version of this paper.

is available at www.nature.com/reprints.

for materials should be addressed to DJW or MAL.

Extended Data Figure 1: Additional radial gas profiles at collapse. The ratio
of the enclosed gas mass to the Jeans mass (a) and gas infall rates (b) are at
the onset of collapse (green), 5 kyr (blue), and at 11 kyr (red). Regions
where the ratio exceeds 1 are where gas is prone to collapse. H2 mass
fractions (c) are at the same times. At 11 kyr they approach unity at the
smallest scales, indicating the nearly complete molecularization of the core.
Extended Data Figure 2: C1 and C2 masses vs. time in the pressure floor run
(solid) and the sink particle run (dashed). The occasional slight dips in C1
mass are due to tidal stripping with other structures in the core. Extended
Data Figure 3: Evolution of SMS 1 and SMS 2. Accretion rates are shown in (a)
and Hertzsprung-Russell tracks are shown in (b). Extended Data Figure 4: Halo
mass vs. redshift. The masses in (a) are from z $=$ 35 (80 Myr after the Big
Bang), when it crosses the threshold for normal Pop III star formation, to z
$=$ 25.7 (126 Myr after the Big Bang), when it begins to dynamically collapse.
The subsequent growth of the halo from z $=$ 15 to z $=$ 6 is shown in (b).
Extended Data Figure 5: Structure of the turbulent core in the sink particle
run. Images are 25 pc on a side and the times in a and b are 0.6 Myr and 1.0
Myr after collapse, respectively. Extended Data Figure 6: Kippenhahn diagrams
for the two stars. Internal structures for SMS 1 and 2 are shown in (a) and
(b), respectively. The x- and y-axes are time and mass coordinate,
respectively, so the interior structure of the star and energy generation at a
given time can be read from the vertical line through the diagram at that
time. Colours show energy generation rates, ${\varepsilon}_{\mathrm{nuc}}$, at
each mass coordinate and red contours mark regions where convection dominates
energy transport (all other regions are radiative). Extended Data Figure 7:
The collapse of SMS 1 and SMS 2 to DCBHs. Infall velocities are shown vs. mass
coordinate in (a) and vs. radius in (b).