11institutetext: Centre for Astrophysics and Space Science Maynooth,
Department of Theoretical Physics, Maynooth University, Maynooth, Ireland.
Universität Heidelberg, Zentrum für Astronomie, Institut für Theoretische
Astrophysik, Albert-Ueberle-Straße 2, D-69120 Heidelberg, Germany.
Universität Heidelberg, Interdisziplinäres Zentrum für Wissenschaftliches
Rechnen, Im Neuenheimer Feld 205, D-69120 Heidelberg, Germany.
Cardiff University School of Physics and Astronomy, United Kingdom.

# Heavy Black Hole Seed Formation in High-z Atomic Cooling Halos

Lewis R. Prole, John A. Regan E-mail<EMAIL_ADDRESS>Simon C. O. Glover
22 Ralf S. Klessen 2233 Felix D. Priestley 44 Paul C. Clark 44 11223344

(Received ; accepted )

###### Abstract

Context. Halos with masses in excess of the atomic limit are believed to be
ideal environments in which to form heavy black hole seeds with masses above
$10^{3}$ M⊙. In cases where the $\rm{H_{2}}$ fraction is suppressed this is
expected to lead to reduced fragmentation of the gas and the generation of a
top heavy initial mass function. In extreme cases this can result in the
formation of massive black hole seeds. Resolving the initial fragmentation
scale and the resulting protostellar masses has, until now, not been robustly
tested.

Aims. We run zoom-in simulations of atomically cooled halos in which the
formation of $\rm{H_{2}}$ is suppressed to assess whether they can truly
resist fragmentation at high densities and tilt the initial mass function
towards a more top heavy form and the formation of massive black hole seeds.

Methods. Cosmological simulations were performed with the moving mesh code
Arepo using a primordial chemistry network until $z\sim 11$. Three haloes with
masses in excess of the atomic cooling mass were then selected for detailed
examination via zoom-ins. A series of zoom-in simulations, with varying levels
of maximum spatial resolution, capture the resulting fragmentation and
formation of metal free stars using the sink particle technique. The highest
resolution simulations resolve densities up to $10^{-6}$ g cm-3 (1018 cm-3)
and capture a further 100 yr of fragmentation behaviour at the center of the
halo. Lower resolution simulations were then used to model the future
accretion behaviour of the sinks over longer timescales.

Results. Our simulations show intense fragmentation in the central region of
the halos, leading to a large number of near-solar mass protostars. Even in
the presence of a super-critical Lyman-Werner radiation field
($J_{LW}>10^{5}J_{21}$) $\rm{H_{2}}$ continues to form within the inner
$\sim$2000 au of the halo. Despite the increased fragmentation the halos
produce a protostellar mass spectrum that peaks at higher masses relative to
standard Population III star forming halos. The most massive protostars have
accretion rates of 10-3-10-1 M⊙ yr-1 after the first 100 years of evolution,
while the total mass of the central region grows at 1 M⊙ yr-1. Lower
resolution zoom-ins show that the total mass of the system continues to
accrete at $\sim$1 M⊙ yr-1 for at least 104 yr, although how this mass is
distributed amongst the rapidly growing number of protostars is unclear.
However, assuming that a fraction of stars can continue to accrete rapidly the
formation of a sub-population of stars with masses in excess of $10^{3}$
$\rm{M_{\odot}}$ is likely in these halos. In the most optimistic case we
predict the formation of heavy black hole seeds with masses in excess of
$10^{4}$ $\rm{M_{\odot}}$ assuming an accretion behaviour in line with
expectations from super-competitive accretion and/or frequent mergers with
secondary protostars.

###### Key Words.:

Stars: Population III – (Galaxies:) quasars: supermassive black holes – Stars:
black holes – (Cosmology:) dark ages, reionization, first stars –
Hydrodynamics

##  1 Introduction

Quasi-stellar radio sources (quasars) are thought to be powered by
supermassive black holes (SMBHs) which populate the centers of most, if not
all, massive galaxies. Quasars have been detected out to redshifts of $z>7$
with black hole (BH) masses believed to be in excess of $10^{9}$ M⊙ (e.g.
Mortlock et al., 2011; Matsuoka et al., 2019), implying that the seeds for
these SMBHs formed in the early Universe. The existence of SMBHs with masses
in excess of $10^{9}$ M⊙ within the first billion years of the Universe poses
a challenge to our understanding of both BH formation and BH accretion. How
could such massive objects appear so early in cosmic history? Two mainstream
pathways have emerged over the last four decades; SMBHs may originate from so-
called light seeds with masses less than 103 $\rm{M_{\odot}}$, or from heavy
BH seeds with masses significantly in excess of 103 $\rm{M_{\odot}}$.

Light seeds are typically thought to form from the remnants of the first stars
- Population III (Pop III) stars. Initial ab-initio modelling of Pop III star
formation predicted an extremely top heavy IMF for the first stars (e.g. Bromm
et al., 1999, 2001; Abel et al., 2002) with characteristic masses of order
1000 $\rm{M_{\odot}}$. However, more recent studies throughout the last decade
have found that Pop III stellar masses are lower than initially suggested
(e.g. Stacy et al. 2010; Clark et al. 2011a; Hirano et al. 2014; Susa 2019;
Prole et al. 2023; for a recent review see Klessen & Glover 2023) with
characteristic masses of a few tens of M⊙. Depending on the exact onset of BH
formation these light seeds would need to maintain accretion at the Eddington
rate for several hundred million years in order to bridge the approximately
seven orders of magnitude in mass required to reach the upper limits of the
SMBH threshold. Alternatively periods of super-Eddington accretion may offer a
solution. In this case the light seed BH can grow at supra-exponential rates
(Alexander & Natarajan, 2014) which only need to last for brief periods of
time (Lupi et al., 2016; Inayoshi et al., 2016). However, in both Eddington
limited and super-Eddington cases, radiative feedback from the accretion of
matter onto BHs heats the surrounding gas and lowers the accretion rate
(Johnson & Bromm, 2007; Milosavljević et al., 2009; Alvarez et al., 2009;
Smith et al., 2018) \- making maintaining either Eddington accretion and/or
super-Eddington accretion unlikely over sustained periods (Regan et al., 2019;
Su et al., 2023; Massonneau et al., 2023). Finally, both Eddington and super-
Eddington limited growth requires that the light seed sit at the centre of a
powerful gas inflow and that the embryonic BH can readily accrete the
surrounding dense gas. However, dynamical studies of BHs have shown that this
is also a challenge, with light seeds tending to walk random trajectories
around the host halo centres (Beckmann et al., 2019; Pfister et al., 2019). As
a result of these obstacles to light seed growth, the possibility of much
heavier BH seeds has also been studied, beginning with Rees (1978).

Forming heavy seeds has hinged on two separate but not necessarily distinct
pathways. On the one hand dynamical processes have been invoked to explain
heavy BH seed formation; runaway collisions in dense young star clusters could
produce massive black holes (MBHs) of $\sim 10^{3}$ M⊙ (e.g. Portegies Zwart
et al., 2004; Glebbeek et al., 2009; Katz et al., 2015; Reinoso et al., 2023)
that are candidates for the ultraluminous X-ray sources observed in young star
forming regions (Ptak & Colbert, 2004). These MBHs may grow into SMBHs through
binary mergers and/or gas accretion (Micic et al., 2007). Additionally, BH
mergers within a dense BH cluster may achieve the same outcome (Stone et al.,
2017; Schleicher et al., 2022) through either collisions and/or gas accretion,
although the growth prospects are far from certain (Arca Sedda et al., 2023).
On the other hand conditions may exist in the early Universe conducive to the
formation of truly massive stars with final masses well in excess of $10^{3}$
$\rm{M_{\odot}}$ (Regan et al., 2020; Latif et al., 2022; Regan et al., 2023)
and possibly up as high as a few times $10^{5}$ $\rm{M_{\odot}}$(Woods et al.,
2017). The formation of these massive primordial Pop III stars requires high
inflow rates onto the stellar surface (Haemmerlé et al., 2018; Woods et al.,
2020) but repeated numerical experiments have shown the stars to be stable for
at least 2 Myr until their inevitable direct collapse into a MBH (e.g.
Hosokawa & Omukai, 2009; Hosokawa et al., 2012a). Finally, it may also be that
the processes that lead to a dense stellar cluster or a massive Pop III star
form part of a continuum with a massive star forming at the very centre of a
dense stellar cluster where collisions drive the formation of a very massive
star (Boekholt et al., 2018; Chon & Omukai, 2020; Schleicher et al., 2023).
The goal of this paper will be to test whether so-called atomic cooling halos
can lead to the formation of a heavy seed Pop III star using state-of-the-art,
high resolution, hydrodynamic simulations.

Pristine atomic cooling halos have long been suggested as the ideal
environment in which to seed MBHs (Loeb & Rasio, 1994; Spaans & Silk, 2006;
Prieto et al., 2013). If the $\rm{H_{2}}$ abundance inside the massive halo
can be suppressed then the gas must cool predominatly through atomic hydrogen
line emission and H- free-bound emission. The thermal pathway then taken by
the gas inside an atomic cooling halo in the absence of effective $\rm{H_{2}}$
cooling deviates significantly from the standard Pop III star formation
scenario (e.g. Omukai 2000; Klessen & Glover 2023).
Pop III star formation in minihalos is facilitated by molecular hydrogen. The
DM halo potential well pulls in the gas and shock heats it up to $\sim 1000$
K, At these temperatures, the H2 abundance increases to $\sim 10^{-4}$ (e.g.
Tegmark et al. 1997; Greif et al. 2008) where rotational transitions can occur
via electrical quadrupole radiation, which allows the gas to cool down to
minimum temperatures of $\sim$200 K (e.g. Abel et al. 1997; Bromm et al. 1999;
Glover & Abel 2008) and collapse, decoupling from the DM halo.
One potential way in which H2 abundances can be reduced is via a nearby source
of Lyman-Werner (LW) radiation. These far-ultraviolet (FUV) photons in the
Lyman and Werner bands of H2 (11.2 – 13.6 eV) can dissociate H2 via the two-
step Solomon process (Field et al., 1966; Stecher & Williams, 1967).
Additionally, photons of above 0.76eV can photodissociate H-, disrupting the
primary H2 formation channel (e.g. Chuzhoy et al. 2007). Before the Strömgren
spheres of Pop III stars overlap, the UV background below the ionization
threshold is able to penetrate large clouds and suppress their H2 abundance
(Haiman et al., 1997). This photodissociation of H2 suppresses further star
formation inside small halos and delays reionization until larger halos form
(Haiman et al., 2000). While other physical mechanism can have a similar
effect, we are focused here in studying the gas collapse inside atomically
cooling halos, which are both pristine and have had their $\rm{H_{2}}$ cooling
efficiency suppressed.

In halos with virial temperatures below $T_{\rm vir}\sim 8000$ K, star
formation is suppressed entirely if LW radiation reduces the H2 abundance
below the level at which gas can cool within a Hubble time. For example, an
intense burst of LW radiation from a neighbouring star-bursting protogalaxy
just before the gas cloud undergoes gravitational collapse is proposed to
prevent the cloud from collapsing or forming stars (Regan et al., 2017).
However, the halo will continue to grow through hierarchical mergers (e.g.
Chon et al. 2016; Dong et al. 2022). Once the virial temperature reaches $\sim
8000$ K, it becomes possible for the gas to cool via Lyman-$\alpha$ emission.
The required virial temperature is related to the virial mass through the
relation (Fernandez et al., 2014)

$\rm T_{vir}=0.75\times
1800\Big{(}\frac{M}{10^{6}M_{\odot}}\Big{)}^{2/3}\Big{(}\frac{1+z}{21}\Big{)},$
(1)

giving a virial mass of $\sim 3\times 10^{7}$ M⊙ at $z=12$ for a virial
temperature of 8000 K. If a halo can grow to this mass it will become hot
enough to cool via atomic line emission and begin to collapse (Oh & Haiman,
2002; Bromm & Loeb, 2003; Bromm & Yoshida, 2011). In this scenario, the
collapse occurs almost isothermally, and fragmentation of the gas is thought
to be suppressed throughout.

The maximum density that can be reached in simulations of this process is
related to the resolution of the simulation. To avoid artificial
fragmentation, it is necessary to resolve the Jeans length, which
progressively shrinks as the gas collapses to higher densities (Truelove et
al., 1997). Therefore, the better the resolution of the simulation, the higher
the density that it can reach. We know from simulations of the standard Pop
III star formation scenario that the formation of the primordial protostar
occurs at densities of 10-6 \- 10-4 g cm-3 where the gas becomes adiabatic
(Omukai, 2000; Machida & Nakamura, 2015). The Jeans length at this stage is
0.01-0.1 au, and so the required resolution is roughly a factor of ten smaller
than this. Despite this, the resolution of most atomically cooled halo
simulations is relatively poor in comparison, with most studies not resolving
their gas past densities of 10-17 g cm-3 (e.g. Shang et al. 2010; Sugimura et
al. 2014; Regan et al. 2014; Hartwig et al. 2015; Agarwal & Khochfar 2015;
Glover 2015; Agarwal et al. 2016; Regan et al. 2017; Dunn et al. 2018). This
resolution is sufficient to determine whether or not H2 cooling is important
during the initial collapse of the gas, but does not allow one to draw
conclusions about the later stages of the collapse. Some studies find evidence
for small-scale fragmentation, even in the absence of effective H2 cooling
(e.g. Becerra et al. 2015, 2018a; Chon et al. 2018; Latif et al. 2020; Patrick
et al. 2023). The number of fragments formed is generally much smaller than
the number found in recent simulations of the standard Pop III star formation
scenario. However, this may be a consequence of the limited peak density: in
most of these studies, the gas density never exceeds $\rho\sim 10^{-10}\>{\rm
g\>cm^{-3}}$, four orders of magnitude smaller than the point at which we
expect the collapse to become adiabatic (Becerra et al., 2018b).

This study aims to provide the most accurate picture of atomically cooled halo
collapse at high densities to date, answering whether atomic cooling halos do
experience reduced fragmentation and higher stellar/BH seed masses compared to
the Pop III minihalo scenario, or whether the fragmentation at high densities
produces protostellar masses similar to what we have seen in simulations of
Pop III star forming minihalos. To that end, we simulate the collapse of
atomically cooled halos from cosmological initial conditions with zoom-in
simulations running up to the lower limit of protostellar formation at 10-6 g
cm-3 and capturing a further $\sim$100 yr of disc fragmentation. We directly
compare our results to a recent study, Prole et al. (2023) (hereafter LP23)
which examined H2 cooling minihalos with the same simulation code, chemical
set-up and maximum resolution as the simulations presented in this work.

The format of this paper is as follows: in §2 we outline the numerical
technique used including the simulation code and chemical network. In §3 we
discuss the collapse of the gas up to point immediately prior to sink
formation, while in §4 we discuss the fragmentation of the gas after the
insertion of sink particles. In §5 we analyse the growth of the stellar system
and compare, via a convergence study, the evolution of the star particles into
main sequence stars and discuss their eventual evolution into (massive) BHs.
In §6 we discuss some caveats before concluding in §7.

##  2 Numerical method

###  2.1 Arepo

The simulations presented here were performed with the moving mesh code Arepo
(Springel, 2010) with a primordial chemistry set-up described in §2.2. Arepo
combines the advantages of adaptive mesh refinement (AMR: Berger & Colella
1989) and smoothed particle hydrodynamics (SPH: Monaghan 1992) with a mesh
made up of a moving, unstructured, Voronoi tessellation of discrete points.
Arepo solves hyperbolic conservation laws of ideal hydrodynamics with a finite
volume approach, based on a second-order unsplit Godunov scheme with an exact
Riemann solver. Automatic and continuous refinement overcome the challenge of
structure growth associated with AMR (e.g. Heitmann et al. 2008).

###  2.2 Chemistry

Collapse of primordial gas is closely linked to the chemistry involved (e.g.
Glover et al. 2006; Yoshida et al. 2007; Glover & Abel 2008; Turk et al.
2011). We therefore use a fully time-dependent chemical network to model the
gas. We use the treatment of primordial chemistry and cooling originally
described in Clark et al. (2011b), but with updated values for some of the
rate coefficients, as summarised in Schauer et al. (2019). The network has 45
chemical reactions to model primordial gas made up of 12 species: H, H+, H-,
H${}^{+}_{2}$ , H2, He, He+, He++, D, D+, HD and free electrons. Optically
thin H2 cooling is modelled as described in Glover & Abel (2008): we first
calculate the rates in the low density ($n\rightarrow 0$) and LTE limits, and
the smoothly interpolate between them as a function of $n/n_{\rm cr}$, where
$n_{\rm cr}$ is the H2 critical number density above which collisions are so
frequent that they keep the populations close to their LTE values. To compute
the H2 cooling rate in the low density limit, we account for the collisions
with H, H2, He, H+ and electrons. To calculate the H2 cooling rate in the
optically thick limit, we use an approach based on the Sobolev approximation
(Yoshida et al., 2006; Clark et al., 2011b). Prior to the simulation, we
compute a grid of optically thick H2 cooling rates as a function of the gas
temperature and H2 column density. During the simulation, if the gas is dense
enough for the H2 cooling to potentially be in the optically thick regime
($\rho>2\times 10^{-16}\>{\rm g\>cm^{-3}}$), we interpolate the H2 cooling
rate from this table, using the local gas temperature and an estimate of the
effective H2 column density computed using the Sobolev approximation. In
addition to H2 cooling, we also account for several other heating and cooling
processes: cooling from atomic hydrogen and helium, collisionally-induced H2
emission, HD cooling, ionisation and recombination, heating and cooling from
changes in the chemical make-up of the gas and from shocks, compression and
expansion of the gas, three-body H2 formation and heating from accretion
luminosity. For reasons of computational efficiency, the network switches off
tracking of deuterium chemistry111Note that HD cooling continues to be
included in the model. at densities above 10-16 g cm-3, instead assuming that
the ratio of HD to H2 at these densities is given by the cosmological D to H
ratio of 2.6$\times$10-5. Note that although our treatment of H2 cooling
accounts for the opacity of the gas at high densities, our treatment of the
effects of other cooling processes, such as H- free-bound emission, does not
currently account for the continuum opacity of the gas. At densities below
$\sim 10^{-8}\>{\rm g\,cm^{-3}}$, this makes little difference to the thermal
evolution of the gas, but it means that we will tend to overestimate the
cooling rate at densities above this value. The adiabatic index of the gas is
computed as a function of chemical composition and temperature with the Arepo
HLLD Riemann solver.

###  2.3 Simulation Setup

As discussed in the Introduction the goal of this study is to investigate the
formation of primordial Pop III stars at the centre of a metal-free atomically
cooling halo. Such halos have previously been investigated as ideal sites for
heavy seed formation (Haiman et al., 1996; Regan & Haehnelt, 2009; Latif et
al., 2013b; Latif & Volonteri, 2015; Latif et al., 2015; Latif et al., 2020;
Wise et al., 2019). To generate the appropriate initial conditions we
generated cosmological initial conditions using MUSIC (Hahn & Abel, 2011). An
initial cosmological simulation was performed within a co-moving box of side
length $1\>h^{-1}\>{\rm Mpc}$ using a $\Lambda$CDM cosmology with parameters
$h=0.6774$, $\Omega_{0}=0.3089$, $\Omega_{\rm b}=0.04864$,
$\Omega_{\Lambda}=0.6911$, $n=0.96$ and $\sigma_{8}=0.8159$ (Planck-
Collaboration et al., 2020). The simulations were initialized at $z=127$ with
an initial dark matter distribution using the transfer functions of Eisenstein
& Hu (1998). The gas distribution was set within Arepo to initially follow the
dark matter (i.e. GENERATE_GAS_IN_ICS = 1). We model the dark matter with
$512^{3}$ particles and the gas was modelled with $512^{3}$ grid cells (prior
to refinement). During the simulation, an additional Jeans refinement
criterion was applied such that the Jeans length of the gas is always resolved
with at least 4 grid cells. Hence, the gas is able to dynamically refine
during the simulation allowing maximum resolution where required.
During this initial phase, we disabled the molecular chemistry functionality
and hence utilised a simpler 6 species model such that H2 and HD abundances
remained at their initial value. This eliminated the need for a LW background
radiation field at this initial stage of the calculation since the goal of
this study is to investigate the idealised case of Pop III formation inside of
a pristine, atomically cooling halo at a virial temperature of approximately
8000 K. We note that applying a super-critical LW field has the same result on
the H2 abundance. As the dominant HD formation pathway relies on the H2
abundance (Nakamura & Umemura, 2002), a super-critical LW field also
suppresses the HD abundance.

At $z\sim 13.3$ the first atomic halo begins to collapse within our box. Two
more, physically distinct halos collapse at z = 11.5. These first three halos
to begin atomically cooling and collapsing were then extracted. The selection
criteria was simple; while H2 rich halos begin cooling down from 1000 K from
$\sim$10-23 g cm-3, reaching 200 K by $\sim$10-22 g cm-3, we have disabled H2
formation, hence the only way for the gas to collapse to these densities is
via atomic cooling. A density threshold for selecting a collapsing halo was
therefore chosen as slightly higher than this density at 10-21 g cm-3. Only
halos containing gas cells with a density exceeding this threshold were
selected for extraction and resimulation. The central coordinate of the halos
was found using a friends of friends (FoF) algorithm and a new box length of 2
kpc (physical units) was cut around it. Using the new box cut from the parent
box, the simulations were restarted as ”zoom-in” simulations. The units were
converted into physical units for the remainder of the zoom-in calculation and
the simulation essentially run as an isolated galaxy simulation with periodic
boundary conditions. Table 1 shows the virial mass, radius and temperature of
each halo as calculated by the FoF algorithm as well as the redshift when it
was extracted, while Figure 1 shows the temperature-density profiles at that
stage and Figure 2 shows the radial profiles of the enclosed mass for both the
gas and the DM.

Figure 1: 2D histograms of the temperature - density profiles for the halos at
the point at which they are extracted from the initial cosmological
simulation, weighted by total gas mass within each 2D bin. The halos have
begun to gravitationally collapse via atomic cooling as seen from the close to
isothermal temperature profile. They differ from the regular H2 minihalo case
by the absence of a sharp drop in temperature to $\sim$200 K beginning at
$\sim 10^{-23}$ g cm-3.

Figure 2: Radial profiles of enclosed gas and DM mass for the three atomic
halos at the point where they are extracted from the initial cosmological
simulation. The baryonic component becomes comparable to the DM within the
central $\sim 100$ pc where it begins to decouple from the DM. In Halo 3, the
baryonic component is dominant over DM on these scales.

For the zoom-in simulations, the full chemistry model (12 species) is again
enabled to get an accurate picture of $\rm{H_{2}}$ formation at very high
densities inside the core of the collapsing halo. We invoke a super-critical
LW radiation field of $J_{LW}=10^{5}J_{21}$ (where $J_{21}$ is in units of
10-21 erg s-1 cm-2 Hz-1 sr-1) to suppress H2 formation. We model the effects
of H2 self-shielding using the TreeCol algorithm (Clark et al., 2012).
We turn-on the 12 species model to examine the impact of $\rm{H_{2}}$
production which is inevitable at the highest densities. In our highest
resolution simulations, we evolve the collapse up to a density of 10-6 g cm-3
before inserting sink particles (see §2.5) and capturing a further $\sim$100
yr of disc fragmentation and accretion behaviour, while in the lowest
resolution simulation we capture $\sim$104 yr of accretion. See Table 2 for a
list of simulation realisations and resolutions employed. The structure of
Halo 3 at the end of the high resolution simulation is visualised in Figure 3,
which shows the density field along with H2 fraction and temperatures at
various zoom-in scales.

Figure 3: Projection images of density, H2 fraction and temperature for Halo 3
at the end of the simulation. From left to right, the zoom-in plots show at
radius of 1 kpc, 10 pc and 0.1 pc, respectively.

Figure 4: 2D histograms of the characteristics of the halos at the point just
before the formation of the first sink particle i.e. when the simulations
approach their maximum density, weighted by total gas mass within each 2D bin.
Top Panel: Temperature-density relation. The collapse is close to isothermal
at approximately 8000 K, with the temperature decreasing by only a factor of
two over more than ten orders of magnitude in density. Bottom Panel: Radial
profiles of density and H2 abundance. The density-radius relationship follows
the $\rho\propto R^{-2}$ relationship expected for an isothermal collapse. The
$\rm{H_{2}}$ abundances are initially negligible at the halo virial radius (R
$\sim$ 100 pc) but increases in the centre of the halo once the density
becomes high enough for self-shielding and three-body H2 formation to become
effective).

Figure 5: Top - Comparison of the build up of H2 in high density regions when a $J_{21}=10^{5}$ and 1010 LW radiation field are used. Bottom - H2 shielding factor calculated using equation 12 from Wolcott-Green et al. (2011). Table 1: Halo virial masses, radii and temperatures along with redshifts at the point where they are extracted from the initial cosmological simulation. Halo | M200 [M⊙] | R200 [kpc] | Tvir [K] | $z$
---|---|---|---|---
1 | 6.89 $\times$ 107 | 1.02 | 13536.60 | 11.5
2 | 5.86 $\times$ 107 | 0.85 | 13890.77 | 13.3
3 | 5.61 $\times$ 107 | 0.95 | 11794.74 | 11.5

###  2.4 Low resolution re-simulations

As our high resolution simulations are only able to evolve for approximately
one hundred years after the formation of the first protostar due to their
extreme computational cost, it is not possible to directly determine the zero
age main sequence mass of the stars formed in the system. We therefore re-run
the simulation of Halo 3 with lower resolution to make an estimate of
protostellar growth on longer timescales. The resolution study presented in
Prole et al. (2022a) found that the total mass accreted across all sink
particles is well estimated by low resolution simulations, with less
fragmentation yielding higher mass protostars. Our lower resolution
simulations can therefore estimate the total mass in sink particles long after
the high resolution simulations end. Since the three atomically cooled halos
examined here have a similar total sink particle mass evolution (see §4), we
only re-ran lower resolution simulations of Halo 3, which should nonetheless
be a good proxy for all three halos.
We chose to re-run the simulation with threshold sink particle creation
densities of 10-13 and 10-10 g cm-3, changing the minimum cell volume and
gravitational softening lengths appropriately (see §2.5 for details). The
three different density threshold simulations will be referred to as high,
medium and low resolution from here on, although we emphasise that even the
low resolution simulations have an extremely high spatial resolution of
$\sim$300 au.

Table 2: Simulation parameters. Left to right - sink particle creation density, accretion radius, gas minimum cell volume and minimum gravitational softening length, DM mass resolution and gravitational softening length. resolution |  $\rho_{\rm sink}$ [g cm-3] | Rsink [au] | Vmin [au3] | Lsoft [au] | MDM [M⊙] | LDM [pc]
---|---|---|---|---|---|---
low | 10-13 | 264.11 | 2.88$\times$ 106 | 132.02 | 793.25 | 1.18
medium | 10-10 | 9.20 | 12.18 | 4.60 | 793.25 | 1.18
high | 10-6 | 0.11 | 2.33$\times$ 10-5 | 0.057 | 793.25 | 1.18

###  2.5 Sink particles

The simulation mesh must be refined during a gravitational collapse to ensure
the local Jeans length is resolved. Once the mesh refines down to its minimum
cell volume, sink particles must be introduced to represent the dense gas,
preventing artificial instability in cells larger than their Jeans length. Our
sink particle implementation was introduced in Wollenberg et al. (2020) and
Tress et al. (2020). A cell is converted into a sink particle if it satisfies
three criteria:

  1. 1.

The cell reaches a threshold density.

  2. 2.

It is sufficiently far away from pre-existing sink particles so that their
accretion radii do not overlap.

  3. 3.

The gas occupying the region inside the sink is gravitationally bound and
collapsing.

Likewise, for the sink particle to accrete mass from surrounding cells, the
cell must meet two criteria:

  1. 1.

The cell lies within the accretion radius.

  2. 2.

It is gravitationally bound to the sink particle.

A sink particle can accrete up to 90$\%$ of a cell’s mass, above which the
cell is removed and the total cell mass is transferred to the sink.

Increasing the threshold density for sink particle creation drastically
increases the degree of fragmentation, reducing the masses of subsequent
secondary protostars (Prole et al., 2022a). However, increasing the sink
particle threshold density also increases the computational challenge beyond
which it is currently intractable. Ideally, sink particles would be introduced
when the gas becomes adiabatic at $\sim 10^{-4}$ g cm-3 (Omukai, 2000).
However, the zero metallicity protostellar model of Machida & Nakamura (2015)
suggests that stellar feedback kicks in to halt collapse at $\sim 10^{-6}$ g
cm-3 (1018 cm-3), so we choose this as our sink particle creation density for
our highest density simulations.

The accretion radius of a sink particle $R_{\text{sink}}$ is chosen to be the
Jeans length $\lambda_{\text{J}}$ corresponding to the sink particle creation
density and corresponding temperature. Taking our high resolution simulation
as an example, at a density of $10^{-6}$ g cm-3, we take the temperature value
from Prole et al. (2022a) of 4460 K to give a Jeans length of $1.67\times
10^{12}$ cm (0.11 au). We set the minimum cell length to fit 8 cells across
the sink particle accretion radius in compliance with the Truelove condition,
by setting a minimum cell volume $V_{\text{min}}=(R_{\text{sink}}/4)^{3}$. The
minimum gravitational softening length for cells and sink particles
$L_{\text{soft}}$ is set to $R_{\text{sink}}/2$. The simulation parameters for
the low, medium and high resolution simulations are summarised in Table 2.

The sink particle treatment also includes the accretion luminosity feedback
from Smith et al. (2011), as implemented in Arepo by Wollenberg et al. (2020).
Stellar internal luminosity is not included in this work, which is not a
problem in our high resolution simulations because the Kelvin-Helmholtz times
of the protostars formed are much longer than the period simulated, meaning
that no stars will have yet begun nuclear burning. This however would affect
our lower resolution simulations, which run for significantly longer periods,
although the sink particles in those simulations (potentially) represent
groups of protostars rather than individual stars. A comprehensive treatment
of protostellar feedback is nonetheless outside the scope of this study.

Lastly, as protostellar mergers are often reported in primordial star forming
simulations (e.g. Greif et al. 2012; Hirano & Bromm 2017; Susa 2019), we
include the treatment of sink particle mergers originally implemented in Prole
et al. (2022a). Sink particles are merged if:

  1. 1.

They lie within each other’s accretion radius.

  2. 2.

They are moving towards each other.

  3. 3.

Their relative accelerations are ¡0.

  4. 4.

They are gravitationally bound to each other.

Since sink particles carry no thermal data, the last criteria simply require
that their gravitational potential exceeds the kinetic energy of the system.
When these criteria are met, the larger of the sinks gains the mass and linear
momentum of a smaller sink, and its position is shifted to the centre of mass
of the system. We allow multiple mergers per time-step. For example, if sink A
is flagged to merge into sink B, and sink B is flagged to merge into sink C,
then both A and B will be merged into sink C simultaneously.

##  3 Initial collapse

Figure 4 summarises the state of the collapse at a point just before the
formation of the first sink particle. The temperature-density relationship in
the top panel shows that the initial contraction up to densities of $10^{-15}$
g cm-3 follows a near isothermal collapse as expected (e.g. Omukai 2000;
Klessen & Glover 2023). Above this density, the temperature drop steepens, but
remains close to isothermal. Here the temperature drops by around a factor of
2 over a density range of $10^{4}$, giving an effective gamma of $\sim 0.95$.
This steepening of the temperature profile was first reported on in the one-
zone calculations of Omukai (2001); at this density, the H2 abundance is still
too low to provide significant cooling, instead the temperature drop
corresponds to the point where cooling becomes dominated by H- free-bound
cooling. The near-isothermal collapse of the gas is further evidenced in the
bottom panel of Figure 4 where we clearly see a $\rho\propto R^{-2}$ scaling
of the density versus radius over more than six orders of magnitude in radius.

From the very bottom panel of Figure 4, the H2 abundance begins to build up,
albeit from very low abundance levels, within the inner 10-2 pc ($\sim$ 2000
au) of the halo. While our LW field strength of $J_{21}=10^{5}$ is already
quite extreme, we run a second realisation of Halo 3 with an extremely high
value of $10^{10}$ (as shown in the top panel of Figure 5) to demonstrate that
above roughly 10-15 g cm-3 LW photo-dissociation can not prevent an increase
in the H2 abundance (nor probably can any other physical process). The three-
body H2 formation rate per unit volume is proportional to $n^{3}$ (where $n$
is number density), whereas the corresponding scalings for the
photodissociation rate and the H2 collisional dissociation rate are $n$ and
$n^{2}$, respectively. Therefore H2 formation will inevitably overcome its
destruction at high densities. In practice then the formation of $\rm{H_{2}}$
at high densities is inevitable.

Furthermore, a build in the H2 fraction at high densities owes to the
exponential nature of the self-shielding. To demonstrate this, we calculate
the H2 shielding factor given in Wolcott-Green et al. (2011) as

$f_{\rm H_{2}}=\frac{1}{(1+\frac{N_{\rm H_{2}}}{2.34\times
10^{19}})^{2.38\times 10^{-1}}}{\rm exp}\left(\frac{-5.2\times 10^{-3}N_{\rm
H_{2}}}{2.34\times 10^{19}}\right),$ (2)

where $N_{\rm H_{2}}$ is the H2 column density. The shielding factor acts as a
transmission factor for the LW radiation, with lower values representing more
effective shielding. We calculate values of $N_{\rm H_{2}}$ by extrapolating
from the column density-number density power law relation shown in Figure 3 of
Wolcott-Green et al. (2011). The shielding factor is given as a function of
density in the bottom panel of Figure 5, which quickly falls to 0 above
densities of 10-15 g cm-3, independently of the LW intensity. This shows that
the core will always be shielded from LW radiation above these densities
regardless of the value of $J_{21}$. A build up of H2 fraction within the
dense core has been seen to various degrees in the radial profiles of previous
studies. For example, the one-zone calculations of Omukai (2001) show a
$\rm{H_{2}}$ fraction of $\sim 0.01$ at a density of $\sim 10^{-8}$ g cm-3,
while simulations by Becerra et al. (2015) reach a $\rm{H_{2}}$ fraction of
$\sim 0.1$ by $\sim 10^{-6}$ g cm-3 and Regan et al. (2017) reaches a
$\rm{H_{2}}$ fraction of 10-4 by their maximum density of $10^{-15}$ g cm-3.
What the results here show is that the very central regions of atomic cooling
halos are highly likely to have small pockets of fully molecular gas.

Figure 6: Radial profile of the mass flux through consecutive shells i.e. the
accretion rate onto the center of the halo. For comparison, results from the
15 Pop III star forming minihalos from LP23 are shown in green.

Figure 7: Time evolution of the ratio of mass of the largest sink to total
mass accreted accross all sinks, mass of the largest sink particle, total mass
accreted onto sink particles and total number of sink particles. For
comparison, results from the 15 Pop III star forming minihalos from LP23 are
shown in green.

Figure 8: Number of enclosed Jeans masses as a function of density, shown at a
time 100 yr after the formation of the first sink particle. The Jeans mass at
each density bin is calculated using the mass weighted average temperature and
density within the bin. The calculated Jeans mass is then compared to the
total mass of gas at or above the density of the bin.

##  4 Fragmentation and the Initial Mass Function

We now move onto to discuss the build of the initial mass function and the
resulting protostellar masses that are found in our simulations at the very
highest densities when we introduce sink particles. In Figure 6 we show the
mass flux into concentric shells just before the formation of the first sink
particle, comparing the standard Pop III star forming minihalos ($\sim 10^{6}$
M⊙) of LP23 with the more massive, atomically cooled halos simulated here.
Within the inner $\sim$100 pc, the accretion rate into the atomic halos
exceeds the minihalo case by a factor of between 10-1000 due to the combined
effects of a stronger DM gravitational well and a larger available reservoir
of baryonic gas. In principle this should lead to the formation of more
massive sinks given the higher accretion rates possible.

Figure 7 shows how the total mass in sink particles, the mass of the most
massive sink particle and the total number of sinks evolves in time. The total
mass in sinks of the system exceeds the minihalo case by almost 2 order of
magnitude due to the higher mass flux. While the initial accretion onto the
most massive sink particle far exceeds the same quantity in the minihalos, the
most massive sink is ejected from the system in all three atomic halos within
the first $\sim$50 yr after its formation. This results in the most massive
sink particle aligning with the upper limit seen in the minihalos. The most
surprising result is that the fragmentation of the gas in the center of the
atomic halos far exceeds the minihalo cases.

As discussed earlier, H- formation becomes the dominant cooling process above
10-15 g cm-3, leading to a slightly steeper decline in temperature with
density. However, the temperature is still higher in the atomically cooled
halos when compared to the minihalo case at the same density (see e.g. Omukai
2000; Yoshida et al. 2006; Prole et al. 2022a), corresponding to a higher
Jeans mass. The increase in fragmentation is therefore attributable to the
higher mass infall rate. For example, assuming fragmentation is perfectly
efficient, the maximum number of fragments that can form is the number of
Jeans masses present within the disc. We can therefore get an estimate of the
fragmentation from Figure 8, which shows the number of enclosed Jeans masses
of gas at or above a given density. At high densities the number of Jeans
masses present in the atomically cooled halos exceeds the minihalo case by an
order of magnitude and roughly matches the number of sink particles formed in
each halo (see Figure 7), which explains the increase in fragmentation when
compared to the minihalos.

The combined sink particle mass functions from all three halos are shown in
Figure 9 at $\sim 100$ yr after the formation of the first sink particle. We
also show the mass functions from the minihalos at the end of the LP23
simulations ($\sim$300 yr). Despite increased fragmentation and a factor of 3
less in time to accrete, the atomic halo mass function peaks at a higher
protostellar mass of 3 M⊙ compared to the 0.3 M⊙ peak in the minihalos. We
note that these protostars will continue to accrete for roughly 104 yr before
the end of the pre-main sequence. While the zero-age main sequence (ZAMS)
masses of these stars are unknown (see §5), we have confirmed that atomically
cooled halos do produce a population of higher mass protostars when compared
to regular Pop III star forming minihalos, at least after the first 100 yr
after the formation of the first protostar. These protostars can then accrete
the available gas in competition with further star formation.

Figure 9: Comparison of the sink particle mass function from the 3 atomic
cooling halos at $\sim$100 yr versus the 15 H2 cooling minihalos from LP23 at
$\sim$300 yr. Power laws of M0.85 and M-2 are superimposed to give the reader
an idea of the slopes involved.

Figure 10: Comparison between the high resolution and low resolution
simulations. Top - time evolution of the total mass of the system. The times
at which new sink particles formed are indicated with star-shaped markers in
the medium and low resolution simulations. We do not show the formation of
sink particles for the high resolution simulations as fragmentation occurs
almost instantly. Bottom - accretion rates onto the most massive surviving
(non-ejected) sink particle. Also shown are the masses of the most massive
surviving sink at the end of the simulations.

##  5 Final stellar and black hole masses

In order to estimate the ZAMS masses of the sinks we run lower resolution
simulations of Halo 3 which allows us to evolve the system over longer
timescales. The top panel of Figure 10 shows how the total mass of the system
of sink particles grows as a function of time across the different resolutions
tested. We show the formation of new sink particles in the 10-13 and 10-10 g
cm-3 simulations as star shaped markers. We do not show the formation of sinks
for the 10-6 g cm-3 cases as fragmentation occurs almost immediately. Clearly
the higher the resolution used, the earlier fragmentation occurs. If we take
$t=0$ to be the time at which the first sink forms in each case then the
second sink forms almost immediately in the highest resolution case, after
$\sim 10$ yr in the $10^{-10}\>{\rm g\,cm^{-3}}$ case and only after $\sim
1000$ yr in the $10^{-13}\>{\rm g\,cm^{-3}}$ case. As the resolution used does
not significantly affect the growth of the total system, the lowest resolution
run shows that the system will continue to grow linearly through the pre-main
sequence phase to reach a mass of $\sim 10^{4}$ M⊙ by 104 yr i.e. there will
be $\sim 10^{4}$ M⊙ available for star formation within the central 264 au.
The increased fragmentation in the higher resolution simulations complicates
how this mass will be distributed amongst the protostars. The bottom panel of
Figure 10 shows the accretion rate onto the most massive surviving sink
particle i.e. the most massive non-ejected protostar. At the end of the high
resolution simulations, the largest survivors in the three halos have
accretion rates in the range 10-3-10-1 M⊙ yr-1 and have masses in the range of
7.5-12 M⊙. Among the two lower resolution simulations the final masses of the
protostars are approximately 1500 M⊙ (after 1000 years) and 25000 M⊙ (after
20,000 years) respectively. If the most massive sink particle from the highest
resolution simulations can accrete or grow at or near the rates found in the
lower resolution simulations then super-massive star formation may be
realised.

The zero-age main sequence (ZAMS) mass of a star depends on the pathway
between the formation of the protostar and the eventual beginning of core
hydrogen burning, which in turn depends strongly on the evolution of the
protostellar accretion rate. When a protostar’s Kelvin–Helmholtz (KH)
timescale (time to radiate away its own gravitational energy) is shorter than
its accretion timescale (time to double its mass by accretion), it will
radiate away its energy and begin to contract, which typically occurs at
masses $\sim$10 M⊙ (Palla & Stahler, 1991; Omukai & Palla, 2003) though this
depends sensitively on the assumed accretion rate (e.g. Nandal et al., 2023).
The contraction causes an increased extreme ultraviolet (EUV) luminosity and
surface temperature, ionizing infalling gas and accelerating it outwards. This
triggers a runaway collapse as the decreased accretion rate causes further
contraction until hydrogen burning begins and the protostar reaches the ZAMS.
If the accretion rate onto the largest sink particle in our high resolution
simulations falls and remains below $\sim 4\times 10^{-3}$ M⊙ yr-1, the KH
contraction will begin at $\sim$10 M⊙ and the growth of the protostar will be
self-regulated by these radiative feedback effects, limiting the final mass to
a few tens of M⊙ (Hosokawa et al., 2011, 2012b). However, numerous studies
have examined the impact that accretion has on the contraction of a Pop III
star (e.g. Hirano et al., 2014). A key quantity here is the critical accretion
rate - this is the accretion rate onto the stellar surface that prevents
contraction of the star. Recently Nandal et al. (2023) investigated in detail
this quantity in terms of episodic accretion rates using the Geneva Stellar
Evolution code (Genec: Eggenberger et al. 2008) . They found that during the
pre-main sequence, which is of most relevance here, the critical accretion
rate is $\sim 2.5\times 10^{-2}$ M⊙ yr-1. In this case the effective surface
temperature remains below $10^{4}$ K, hence feedback can not form a HII region
and accretion is not prevented, leading to the formation of a super-giant
protostar, which can grow up to several thousand solar masses depending on the
details of the accretion (Umeda et al., 2016; Woods et al., 2017; Haemmerlé et
al., 2018). The end of accretion onto a super-giant protostar is caused by a
fast contraction when the accretion rate falls below $\sim 7\times 10^{-3}$ M⊙
yr-1.

At the end of the high resolution simulations, the accretion rates onto the
largest surviving sink particles vary significantly in time. As the accretion
rates onto $\sim$10 M⊙ protostars appear to be a good indicator of future
contraction/accretion evolution (Hosokawa et al., 2011, 2012b; Hirano & Bromm,
2017) and the largest survivor in Halos 2 and 3 have accretion rates between
10-2-10-1 M⊙ yr-1, they could feasibly go on to form (super-massive) Pop III
stars with masses anywhere in the range of 10 M⊙ up to $10^{4}$ M⊙ within
their main sequence (MS) lifetime. If they grow in excess of 25 M⊙ the type II
supernova explosion would be too weak to eject much of the star and the
subsequent fallback of material causes the resulting neutron star to collapse
into a BH (MacFadyen et al., 2001), while above 40 M⊙ the neutron star is
unable to launch a supernova shock and collapses directly to form a BH with no
mass loss (Fryer, 1999). If these stars avoid the disruptive pair instability
supernova mass range of 140-260 M⊙, they will result in a massive BH equal in
mass to the stellar progenitor. However, whether the protostars will maintain
their accretion rates is uncertain. In the optimistic case that they maintain
accretion rates of 10-1 M⊙ yr-1, these results are in line with many previous,
lower resolution simulations of heavy seed black hole formation (e.g. Johnson
et al. 2011; Latif et al. 2013a, b; Regan et al. 2014; Latif & Volonteri 2015;
Choi et al. 2015; Regan et al. 2017; Smidt et al. 2017; Ardaneh et al. 2018)
despite the increased gas fragmentation. The largest survivor in Halo 1 ends
with an accretion rate of $\sim$10-3 M⊙ yr-1, which if maintained will likely
lead to an early KH contraction and limit the final mass to 10-30 M⊙. It is
unclear if the accretion rate onto Halo 1 will remain below 10-3 or if it will
experience an increase similar to that of Halo 3. Certainly the formation of
massive BH seeds of 104 M⊙ needed to explain $z\sim 7$ quasar observations
would rely on frequent mergers with other protostars.

The sink particle masses in the lower resolution simulations represent whole
groups of protostars in the high resolution simulations. Frequent mergers of
secondary protostars into the main protostar would push the stellar and later
BH masses up to close to what is achieved in the low resolution simulations,
with an upper limit of 104 M⊙ within a 10 kyr period. Recent studies have
pointed out the importance of super competitive accretion, in which a central
few stars grow supermassive while a large number of other stars are competing
for the gas reservoir (Chon & Omukai, 2020), with central objects growing
significantly through mergers (Sassano et al., 2021; Vergara et al., 2021;
Trinca et al., 2022; Schleicher et al., 2022; Zwick et al., 2023). Here the
mass growth by collisions can be comparable to its growth via accretion
(Schleicher et al., 2023). The super-competitive accretion scenario was
initially invoked to explore the ’low metallicity’ regime ($\rm{Z}\sim
10^{-6}-10^{-3}$ Z⊙) where fragmentation is expected to be more active
compared to the metal-free case. What we find here is that fragmentation is
already vigorous within the central 2000 au of the halo. It therefore appears
that for all metallicities below $\sim 10^{-3}$ Z⊙ we can expect a scenario
where vigorous fragmentation competes with a rapidly growing protostar.
Differentiating the environments which produce a dense stellar cluster versus
those that produce a central massive Pop III star may be the next frontier.

To investigate the merger behaviour further, we plot the total number of sink
mergers, ratio of total merged mass to accreted sink mass and the ratio of the
number of mergers to the number of surviving sink particles as a function of
time in Figure 11. From this we see that the mass gained through mergers
constitutes between 1-10$\%$ of the total sink particle mass. The frequency of
mergers is high, with the number of mergers totalling between 10-50$\%$ of the
total number of surviving sink particles. If this holds throughout the
following 104 yr, the total mass gained through mergers alone could be as high
as 103 M⊙, although it is unclear if this mass would be shared between the
growing number of protostars randomly or be preferentially received by the
largest growing BH seed.

Figure 11: Time evolution of the merger data. Top - cumulative number of sink
particle mergers. Middle - ratio of cumulative merged mass to the total
accreted sink particle mass (M${}_{\rm tot}-$M0) where M0 is the initial mass
of a sink particle. Bottom - ratio of cumulative number of sink mergers
against the current number of surviving sink particles.

##  6 Caveats

We have not included magnetic fields in these simulations. While studies of
primordial magnetic fields suggest that they can increase the mass of
protostars (e.g. Saad et al. 2022; Hirano & Machida 2022; Stacy et al. 2022),
the fields have no effect when they are properly resolved, distributing the
magnetic energy from the small-scale turbulent dynamo to smaller spatial
scales (Prole et al., 2022b).

We have also not included radiative feedback from our protostars, which has
the effect of heating the surrounding gas and lowering accretion rates. While
our high resolution simulations end at a point when the protostars would not
produce significant levels of feedback. The low resolution simulations however
run for much longer and would be subject to feedback effects.

As our high resolution simulations have only captured the first $\sim 100$ yr
of accretion after the formation of the first protostar and the lower
resolution simulations lack the capacity to resolve individual star foratmion,
to what extent competitive accretion and mergers affect the subsequent final
stellar/BH masses is unclear.

Ideally we would resolve up to the protostellar formation density of 10-4 g
cm-3, which is currently computationally unfeasible. Failure to resolve up to
that density means we have not achieved numerical convergence. Despite this,
our implemented threshold density of 10-6 g cm-3 and post-sink particle run
time of 100 yr represents the current state of the art in this regime.

As mentioned in §2.2, although we account self-consistently for the effects of
opacity when computing the H2 cooling rate (due to line cooling and collision-
induced emission), we do not yet account for the continuum opacity of the gas
when computing the cooling rate due to other radiative processes, e.g. H-
free-bound cooling. At most densities in our models, the optical depth of the
gas in the continuum is small and this simplification has little impact – for
example, it should be valid at all of the densities traced in our medium and
low resolution models. However, above a density of $\sim 10^{-8}\>{\rm
g\>cm^{-3}}$, we expect the optical depth of the gas in the continuum to
become significant, and hence our current simulations over-estimate the
cooling rate of the gas above this density. This will likely give us slightly
more effective fragmentation at the highest densities than we would find in
reality. However, we do not expect this to significantly impact the main
conclusions of our study. As 8 demonstrates, the number of Jeans masses
present in the gas at $\rho=10^{-8}\>{\rm g\>cm^{-3}}$ in the atomic cooling
halos already substantially exceeds the corresponding value in H2-cooled
minihalos, demonstrating that this result is not a consequence of our
simplified treatment of very high density cooling. Our finding that the gas
fragments extensively in the centre of the atomic-cooled halos should
therefore be a robust result, although some of the details (e.g. the precise
number of fragments and their formation masses) will depend on the treatment
of the cooling at very high densities. Our finding that the fragments quickly
grow to become more massive than their Pop III counterparts should also be a
robust result: properly accounting for the continuum opacity will likely give
us fragments with slightly larger initial masses, which if anything will
increase the rate at which they later gain mass, exacerbating the difference
between the atomic-cooled and H2-cooled results.

We have assumed that the background radiation field from Pop III stars takes
the form of a black body spectrum with an effective temperature of 105 K,
peaking in the UV. However, this relies on the stars having masses of $\sim$
300 M⊙ (Schaerer, 2002). The masses of Pop III stars are currently uncertain,
so our choice of effective temperature affects the photodissociation of the
molecules present in the gas. For example, the effective temperature for a 1
M⊙ Pop III star is $\sim$ 7000 K (Larkin et al., 2022), which peaks at lower
frequencies, resulting in significantly less H2 photodissociating LW radiation
but more infrared radiation, which can still photodissociate H- to disrupt H2
formation (Chuzhoy et al., 2007).

The methodology here was idealised in that we used a six species chemical
model to construct a pristine atomic cooling halo before switching to a full
12 species chemical model and employing an intense LW background. In future
work we will apply the same methodology to rapidly assembling halos without
the reliance on super-critical LW values.

##  7 Conclusions

We have performed high resolution simulations of three pristine atomically
cooled halos of mass $\sim 10^{8}$ M⊙ which begin to collapse at $z\sim
11.5-13.3$. The goal of our simulation suite was to examine the fragmentation
and protostellar formation within the core of such systems as they reach
protostellar densities. For the highest resolution we could tractably simulate
($\rho\sim 10^{-6}$ g cm-3) we followed the protostar formation for
approximately 100 years by introducing sink particles representing individual
protostars. The main results are as follows:

  1. 1.

Resolving protstellar densities reveals that the center of (pristine)
atomically cooled halos is subject to intense fragmentation, even more-so than
in the canonical Pop III minihalo case. This is driven by H- free-bound
cooling and an increased number of Jeans masses within the core, owing to
enhanced accretion rates.

  2. 2.

Despite increased fragmentation, the atomically cooled halos formed a
population of higher mass protostars compared to Pop III star forming
minihalos.

  3. 3.

The initial accretion rates onto the most massive surviving protostars
indicates that their zero-age main sequence masses could range from
$10^{2}-10^{4}$ M⊙ depending on subsequent accretion and/or mergers.

  4. 4.

Our coarser resolution simulations show that the total mass of the system
continues to grow steadily for at least 104 years after the formation of the
first protostar (achieving central sink particle masses in excess of $10^{4}$
$\rm{M_{\odot}}$), although how this is distributed amongst the rapidly
growing number of protostars is unknown and requires further study.

  5. 5.

The formation of a massive Pop III star (M${}_{*}\gtrsim 1000$
$\rm{M_{\odot}}$) is therefore realistic and achievable in a high-z setting
but relies on (super-competitive) accretion and frequent mergers with
secondary protostars

  6. 6.

The H2 fraction within the inner $\sim$2000 au of all three halos was able to
build up independently of the strength of the external LW radiation field.
This strongly indicates that all collapsing halos (even those subject to
strong LW fields and likely other physical processes which could suppress star
formation up to the atomic limit e.g. streaming velocities between DM and
baryonic gas) contain small pockets of cold, H2 rich gas at their centre.

The findings here show that the formation of massive Pop III stars (i.e. heavy
seed MBHs) is entirely plausible in sufficiently massive halos but will depend
sensitively on the future evolution of the protostars and on the balance
between stellar accretion and mergers in the dense cluster that forms within
the center of the atomic cooling halo. While we model the impact of a super-
critical LW radiation field here, any physical process which results in the
initial suppression of star formation up to the atomic cooling limit should
have a similar impact and allow for the formation of more massive Pop III
stars.

###### Acknowledgements.

LP and JR acknowledge support from the Irish Research Council Laureate
programme under grant number IRCLA/2022/1165. JR also acknowledges support
from the Royal Society and Science Foundation Ireland under grant number
URF\R1\191132. This work used the DiRAC@Durham facility managed by the
Institute for Computational Cosmology on behalf of the STFC DiRAC HPC Facility
(www.dirac.ac.uk). The equipment was funded by BEIS capital funding via STFC
capital grants ST/P002293/1, ST/R002371/1 and ST/S002502/1, Durham University
and STFC operations grant ST/R000832/1. DiRAC is part of the National
e-Infrastructure. We also acknowledge the support of the Supercomputing Wales
project, which is part-funded by the European Regional Development Fund (ERDF)
via Welsh Government. SCOG and RSK acknowledge financial support from the
European Research Council via the ERC Synergy Grant “ECOGAL” (project ID
855130), from the Heidelberg Cluster of Excellence (EXC 2181 – 390900948)
“STRUCTURES”, funded by the German Excellence Strategy, and by the German
Ministry for Economic Affairs and Climate Action in project “MAINN” (funding
ID 50OO2206). They also thank for computing resources provided by the Ministry
of Science, Research and the Arts (MWK) of the State of Baden-Württemberg and
Deutsche Forschungsgemeinschaft (DFG) through grant INST 35/1134-1 FUGG and
for data storage at SDS@hd through grant INST 35/1314-1 FUGG. Finally, we
acknowledge Advanced Research Computing at Cardiff (ARCCA) for providing
resources for the project.

## References

  * Abel et al. (1997) Abel T., Anninos P., Zhang Y., Norman M. L., 1997, New Astronomy, 2, 181
  * Abel et al. (2002) Abel T., Bryan G. L., Norman M. L., 2002, Science (New York, N.Y.), 295, 93
  * Agarwal & Khochfar (2015) Agarwal B., Khochfar S., 2015, Monthly Notices of the Royal Astronomical Society, 446, 160
  * Agarwal et al. (2016) Agarwal B., Smith B., Glover S., Natarajan P., Khochfar S., 2016, Monthly Notices of the Royal Astronomical Society, 459, 4209
  * Alexander & Natarajan (2014) Alexander T., Natarajan P., 2014, Science (New York, N.Y.), 345, 1330
  * Alvarez et al. (2009) Alvarez M. A., Wise J. H., Abel T., 2009, The Astrophysical Journal, 701, L133
  * Arca Sedda et al. (2023) Arca Sedda M., Mapelli M., Benacquista M., Spera M., 2023, Monthly Notices of the Royal Astronomical Society, 520, 5259
  * Ardaneh et al. (2018) Ardaneh K., Luo Y., Shlosman I., Nagamine K., Wise J. H., Begelman M. C., 2018, Monthly Notices of the Royal Astronomical Society, 479, 2277
  * Becerra et al. (2015) Becerra F., Greif T. H., Springel V., Hernquist L. E., 2015, Monthly Notices of the Royal Astronomical Society, 446, 2380
  * Becerra et al. (2018a) Becerra F., Marinacci F., Bromm V., Hernquist L. E., 2018a, Monthly Notices of the Royal Astronomical Society, 480, 5029
  * Becerra et al. (2018b) Becerra F., Marinacci F., Inayoshi K., Bromm V., Hernquist L. E., 2018b, The Astrophysical Journal, 857, 138
  * Beckmann et al. (2019) Beckmann R. S., et al., 2019, Astronomy and Astrophysics, 631, A60
  * Berger & Colella (1989) Berger M. J., Colella P., 1989, Journal of Computational Physics, 82, 64
  * Boekholt et al. (2018) Boekholt T. C. N., Schleicher D. R. G., Fellhauer M., Klessen R. S., Reinoso B., Stutz A. M., Haemmerlé L., 2018, Monthly Notices of the Royal Astronomical Society, 476, 366
  * Bromm & Loeb (2003) Bromm V., Loeb A., 2003, The Astrophysical Journal, 596, 34
  * Bromm & Yoshida (2011) Bromm V., Yoshida N., 2011, Annual Review of Astronomy and Astrophysics, 49, 373
  * Bromm et al. (1999) Bromm V., Coppi P. S., Larson R. B., 1999, The Astrophysical Journal, 527, L5
  * Bromm et al. (2001) Bromm V., Kudritzki R. P., Loeb A., 2001, The Astrophysical Journal, 552, 464
  * Choi et al. (2015) Choi J.-H., Shlosman I., Begelman M. C., 2015, Monthly Notices of the Royal Astronomical Society, 450, 4411
  * Chon & Omukai (2020) Chon S., Omukai K., 2020, Monthly Notices of the Royal Astronomical Society, 494, 2851
  * Chon et al. (2016) Chon S., Hirano S., Hosokawa T., Yoshida N., 2016, The Astrophysical Journal, 832, 134
  * Chon et al. (2018) Chon S., Hosokawa T., Yoshida N., 2018, Monthly Notices of the Royal Astronomical Society, 475, 4104
  * Chuzhoy et al. (2007) Chuzhoy L., Kuhlen M., Shapiro P. R., 2007, The Astrophysical Journal, 665, L85
  * Clark et al. (2011a) Clark P. C., Glover S. C. O., Smith R. J., Greif T. H., Klessen R. S., Bromm V., 2011a, Science (New York, N.Y.), 331, 1040
  * Clark et al. (2011b) Clark P. C., Glover S. C. O., Smith R. J., Greif T. H., Klessen R. S., Bromm V., 2011b, Science (New York, N.Y.), 331, 1040
  * Clark et al. (2012) Clark P. C., Glover S. C. O., Klessen R. S., 2012, Monthly Notices of the Royal Astronomical Society, 420, 745
  * Dong et al. (2022) Dong F., Zhao D., Han J., Li Z., Jing Y., Yang X., 2022, The Astrophysical Journal, 929, 120
  * Dunn et al. (2018) Dunn G., Bellovary J., Holley-Bockelmann K., Christensen C., Quinn T., 2018, The Astrophysical Journal, 861, 39
  * Eggenberger et al. (2008) Eggenberger P., Meynet G., Maeder A., Hirschi R., Charbonnel C., Talon S., Ekström S., 2008, Astrophysics and Space Science, 316, 43
  * Eisenstein & Hu (1998) Eisenstein D. J., Hu W., 1998, The Astrophysical Journal, 496, 605
  * Fernandez et al. (2014) Fernandez R., Bryan G. L., Haiman Z., Li M., 2014, Monthly Notices of the Royal Astronomical Society, 439, 3798
  * Field et al. (1966) Field G. B., Somerville W. B., Dressler K., 1966, Annual Review of Astronomy and Astrophysics, 4, 207
  * Fryer (1999) Fryer C. L., 1999, The Astrophysical Journal, 522, 413
  * Glebbeek et al. (2009) Glebbeek E., Gaburov E., de Mink S. E., Pols O. R., Zwart S. F. P., 2009, Astronomy & Astrophysics, 497, 255
  * Glover (2015) Glover S. C. O., 2015, Monthly Notices of the Royal Astronomical Society, 451, 2082
  * Glover & Abel (2008) Glover S. C. O., Abel T., 2008, Monthly Notices of the Royal Astronomical Society, 388, 1627
  * Glover et al. (2006) Glover S. C., Savin D. W., Jappsen A.-K., 2006, The Astrophysical Journal, 640, 553
  * Greif et al. (2008) Greif T. H., Johnson J. L., Klessen R. S., Bromm V., 2008, Monthly Notices of the Royal Astronomical Society, 387, 1021
  * Greif et al. (2012) Greif T. H., Bromm V., Clark P. C., Glover S. C. O., Smith R. J., Klessen R. S., Yoshida N., Springel V., 2012, Monthly Notices of the Royal Astronomical Society, 424, 399
  * Haemmerlé et al. (2018) Haemmerlé L., Woods T. E., Klessen R. S., Heger A., Whalen D. J., 2018, Monthly Notices of the Royal Astronomical Society, 474, 2757
  * Hahn & Abel (2011) Hahn O., Abel T., 2011, Monthly Notices of the Royal Astronomical Society, 415, 2101
  * Haiman et al. (1996) Haiman Z., Thoul A. A., Loeb A., 1996, The Astrophysical Journal, 464, 523
  * Haiman et al. (1997) Haiman Z., Rees M. J., Loeb A., 1997, The Astrophysical Journal, 476, 458
  * Haiman et al. (2000) Haiman Z., Abel T., Rees M. J., 2000, The Astrophysical Journal, 534, 11
  * Hartwig et al. (2015) Hartwig T., Glover S. C. O., Klessen R. S., Latif M. A., Volonteri M., 2015, Monthly Notices of the Royal Astronomical Society, 452, 1233
  * Heitmann et al. (2008) Heitmann K., et al., 2008, Computational Science and Discovery, 1, 015003
  * Hirano & Bromm (2017) Hirano S., Bromm V., 2017, Monthly Notices of the Royal Astronomical Society, 470, 898
  * Hirano & Machida (2022) Hirano S., Machida M. N., 2022, The Astrophysical Journal, 935, L16
  * Hirano et al. (2014) Hirano S., Hosokawa T., Yoshida N., Umeda H., Omukai K., Chiaki G., Yorke H. W., 2014, The Astrophysical Journal, 781, 60
  * Hosokawa & Omukai (2009) Hosokawa T., Omukai K., 2009, The Astrophysical Journal, 691, 823
  * Hosokawa et al. (2011) Hosokawa T., Omukai K., Yoshida N., Yorke H. W., 2011, Science (New York, N.Y.), 334, 1250
  * Hosokawa et al. (2012a) Hosokawa T., Omukai K., Yorke H. W., 2012a, The Astrophysical Journal, 756, 93
  * Hosokawa et al. (2012b) Hosokawa T., Omukai K., Yorke H. W., 2012b, The Astrophysical Journal, 756, 93
  * Inayoshi et al. (2016) Inayoshi K., Haiman Z., Ostriker J. P., 2016, Monthly Notices of the Royal Astronomical Society, 459, 3738
  * Johnson & Bromm (2007) Johnson J. L., Bromm V., 2007, Monthly Notices of the Royal Astronomical Society, 374, 1557
  * Johnson et al. (2011) Johnson J. L., Khochfar S., Greif T. H., Durier F., 2011, Monthly Notices of the Royal Astronomical Society, 410, 919
  * Katz et al. (2015) Katz H., Sijacki D., Haehnelt M. G., 2015, Monthly Notices of the Royal Astronomical Society, 451, 2352
  * Klessen & Glover (2023) Klessen R. S., Glover S. C. O., 2023, Annual Review of Astronomy and Astrophysics, 61, 65
  * Larkin et al. (2022) Larkin M. M., Gerasimov R., Burgasser A. J., 2022, The Astronomical Journal, 165, 2
  * Latif & Volonteri (2015) Latif M. A., Volonteri M., 2015, Monthly Notices of the Royal Astronomical Society, 452, 1026
  * Latif et al. (2013a) Latif M. A., Schleicher D. R. G., Schmidt W., Niemeyer J., 2013a, Monthly Notices of the Royal Astronomical Society, 433, 1607
  * Latif et al. (2013b) Latif M. A., Schleicher D. R. G., Schmidt W., Niemeyer J. C., 2013b, Monthly Notices of the Royal Astronomical Society, 436, 2989
  * Latif et al. (2015) Latif M. A., Bovino S., Grassi T., Schleicher D. R. G., Spaans M., 2015, Monthly Notices of the Royal Astronomical Society, 446, 3163
  * Latif et al. (2020) Latif M. A., Khochfar S., Whalen D., 2020, The Astrophysical Journal, 892, L4
  * Latif et al. (2022) Latif M. A., Whalen D. J., Khochfar S., Herrington N. P., Woods T. E., 2022, Nature, 607, 48
  * Loeb & Rasio (1994) Loeb A., Rasio F. A., 1994, The Astrophysical Journal, 432, 52
  * Lupi et al. (2016) Lupi A., Haardt F., Dotti M., Fiacconi D., Mayer L., Madau P., 2016, Monthly Notices of the Royal Astronomical Society, 456, 2993
  * MacFadyen et al. (2001) MacFadyen A. I., Woosley S. E., Heger A., 2001, The Astrophysical Journal, 550, 410
  * Machida & Nakamura (2015) Machida M. N., Nakamura T., 2015, Monthly Notices of the Royal Astronomical Society, 448, 1405
  * Massonneau et al. (2023) Massonneau W., Volonteri M., Dubois Y., Beckmann R. S., 2023, Astronomy and Astrophysics, 670, A180
  * Matsuoka et al. (2019) Matsuoka Y., et al., 2019, The Astrophysical Journal, 872, L2
  * Micic et al. (2007) Micic M., Holley-Bockelmann K., Sigurdsson S., Abel T., 2007, Monthly Notices of the Royal Astronomical Society, 380, 1533
  * Milosavljević et al. (2009) Milosavljević M., Couch S. M., Bromm V., 2009, The Astrophysical Journal, 696, L146
  * Monaghan (1992) Monaghan J. J., 1992, Annual Review of Astronomy and Astrophysics, 30, 543
  * Mortlock et al. (2011) Mortlock D. J., et al., 2011, Nature, 474, 616
  * Nakamura & Umemura (2002) Nakamura F., Umemura M., 2002, The Astrophysical Journal, 569, 549
  * Nandal et al. (2023) Nandal D., Regan J., Woods T., Farrell E., Ekström S., Meynet G., 2023, Astronomy & Astrophysics
  * Oh & Haiman (2002) Oh S. P., Haiman Z., 2002, The Astrophysical Journal, 569, 558
  * Omukai (2000) Omukai K., 2000, The Astrophysical Journal, 534, 809
  * Omukai (2001) Omukai K., 2001, The Astrophysical Journal, 546, 635
  * Omukai & Palla (2003) Omukai K., Palla F., 2003, The Astrophysical Journal, 589, 677
  * Palla & Stahler (1991) Palla F., Stahler W S., 1991, Astrophysical Journal, 375, 288
  * Patrick et al. (2023) Patrick S. J., Whalen D. J., Elford J. S., Latif M. A., 2023, Monthly Notices of the Royal Astronomical Society, 522, 3795
  * Pfister et al. (2019) Pfister H., Volonteri M., Dubois Y., Dotti M., Colpi M., 2019, Monthly Notices of the Royal Astronomical Society, 486, 101
  * Planck-Collaboration et al. (2020) Planck-Collaboration et al., 2020, Astronomy and Astrophysics, 641, A6
  * Portegies Zwart et al. (2004) Portegies Zwart S. F., Baumgardt H., Hut P., Makino J., McMillan S. L. W., 2004, Nature, 428, 724
  * Prieto et al. (2013) Prieto J., Jimenez R., Haiman Z., 2013, Monthly Notices of the Royal Astronomical Society, 436, 2301
  * Prole et al. (2022a) Prole L. R., Clark P. C., Klessen R. S., Glover S. C. O., 2022a, Monthly Notices of the Royal Astronomical Society, 510, 4019
  * Prole et al. (2022b) Prole L. R., Clark P. C., Klessen R. S., Glover S. C. O., Pakmor R., 2022b, Monthly Notices of the Royal Astronomical Society, 516, 2223
  * Prole et al. (2023) Prole L. R., Schauer A. T. P., Clark P. C., Glover S. C. O., Priestley F. D., Klessen R. S., 2023, Monthly Notices of the Royal Astronomical Society, 520, 2081
  * Ptak & Colbert (2004) Ptak A., Colbert E., 2004, The Astrophysical Journal, 606, 291
  * Rees (1978) Rees M. J., 1978, Physica Scripta, 17, 193
  * Regan & Haehnelt (2009) Regan J. A., Haehnelt M. G., 2009, Monthly Notices of the Royal Astronomical Society, 396, 343
  * Regan et al. (2014) Regan J. A., Johansson P. H., Wise J. H., 2014, The Astrophysical Journal, 795, 137
  * Regan et al. (2017) Regan J. A., Visbal E., Wise J. H., Haiman Z., Johansson P. H., Bryan G. L., 2017, Nature Astronomy, 1, 1
  * Regan et al. (2019) Regan J. A., Downes T. P., Volonteri M., Beckmann R., Lupi A., Trebitsch M., Dubois Y., 2019, Monthly Notices of the Royal Astronomical Society, 486, 3892
  * Regan et al. (2020) Regan J. A., Wise J. H., Woods T. E., Downes T. P., O’Shea B. W., Norman M. L., 2020, The Open Journal of Astrophysics, 3, 10.21105/astro.2008.08090
  * Regan et al. (2023) Regan J. A., Pacucci F., Bustamante-Rosell M. J., 2023, Monthly Notices of the Royal Astronomical Society, 518, 5997
  * Reinoso et al. (2023) Reinoso B., Klessen R. S., Schleicher D., Glover S. C. O., Solar P., 2023, Monthly Notices of the Royal Astronomical Society, 521, 3553
  * Saad et al. (2022) Saad C. R., Bromm V., El Eid M., 2022, Monthly Notices of the Royal Astronomical Society, 516, 3130
  * Sassano et al. (2021) Sassano F., Schneider R., Valiante R., Inayoshi K., Chon S., Omukai K., Mayer L., Capelo P. R., 2021, Monthly Notices of the Royal Astronomical Society, 506, 613
  * Schaerer (2002) Schaerer D., 2002, Astronomy & Astrophysics, 382, 28
  * Schauer et al. (2019) Schauer A. T. P., Glover S. C. O., Klessen R. S., Ceverino D., 2019, Monthly Notices of the Royal Astronomical Society, 484, 3510
  * Schleicher et al. (2022) Schleicher D. R. G., et al., 2022, Monthly Notices of the Royal Astronomical Society, 512, 6192
  * Schleicher et al. (2023) Schleicher D. R. G., Reinoso B., Klessen R. S., 2023, Monthly Notices of the Royal Astronomical Society, 521, 3972
  * Shang et al. (2010) Shang C., Bryan G. L., Haiman Z., 2010, Monthly Notices of the Royal Astronomical Society, 402, 1249
  * Smidt et al. (2017) Smidt J., Whalen D. J., Johnson J. L., Li H., 2017, preprint, 1703
  * Smith et al. (2011) Smith R. J., Glover S. C. O., Clark P. C., Greif T., Klessen R. S., 2011, Monthly Notices of the Royal Astronomical Society, 414, 3633
  * Smith et al. (2018) Smith B. D., Regan J. A., Downes T. P., Norman M. L., O’Shea B. W., Wise J. H., 2018, Monthly Notices of the Royal Astronomical Society, 480, 3762
  * Spaans & Silk (2006) Spaans M., Silk J., 2006, The Astrophysical Journal, 652, 902
  * Springel (2010) Springel V., 2010, Monthly Notices of the Royal Astronomical Society, 401, 791
  * Stacy et al. (2010) Stacy A., Greif T. H., Bromm V., 2010, Monthly Notices of the Royal Astronomical Society, 403, 45
  * Stacy et al. (2022) Stacy A., McKee C. F., Lee A. T., Klein R. I., Li P. S., 2022, Monthly Notices of the Royal Astronomical Society, 511, 5042
  * Stecher & Williams (1967) Stecher T. P., Williams D. A., 1967, The Astrophysical Journal, 149, L29
  * Stone et al. (2017) Stone N. C., Küpper A. H. W., Ostriker J. P., 2017, Monthly Notices of the Royal Astronomical Society, 467, 4180
  * Su et al. (2023) Su K.-Y., Bryan G. L., Haiman Z., Somerville R. S., Hayward C. C., Faucher-Giguère C.-A., 2023, Monthly Notices of the Royal Astronomical Society, 520, 4258
  * Sugimura et al. (2014) Sugimura K., Omukai K., Inoue A. K., 2014, Monthly Notices of the Royal Astronomical Society, 445, 544
  * Susa (2019) Susa H., 2019, The Astrophysical Journal, 877, 99
  * Tegmark et al. (1997) Tegmark M., Silk J., Rees M. J., Blanchard A., Abel T., Palla F., 1997, The Astrophysical Journal, 474, 1
  * Tress et al. (2020) Tress R. G., Smith R. J., Sormani M. C., Glover S. C. O., Klessen R. S., Mac Low M.-M., Clark P. C., 2020, Monthly Notices of the Royal Astronomical Society, 492, 2973
  * Trinca et al. (2022) Trinca A., Schneider R., Valiante R., Graziani L., Zappacosta L., Shankar F., 2022, Monthly Notices of the Royal Astronomical Society, 511, 616
  * Truelove et al. (1997) Truelove J. K., Klein R. I., McKee C. F., Holliman II J. H., Howell L. H., Greenough J. A., 1997, The Astrophysical Journal Letters, 489, L179
  * Turk et al. (2011) Turk M. J., Clark P., Glover S. C. O., Greif T. H., Abel T., Klessen R., Bromm V., 2011, Astrophys.J.726:55,2011, 726
  * Umeda et al. (2016) Umeda H., Hosokawa T., Omukai K., Yoshida N., 2016, The Astrophysical Journal Letters, 830, L34
  * Vergara et al. (2021) Vergara M. Z. C., Schleicher D. R. G., Boekholt T. C. N., Reinoso B., Fellhauer M., Klessen R. S., Leigh N. W. C., 2021, Astronomy and Astrophysics, 649, A160
  * Wise et al. (2019) Wise J. H., Regan J. A., O’Shea B. W., Norman M. L., Downes T. P., Xu H., 2019, Nature, 566, 85
  * Wolcott-Green et al. (2011) Wolcott-Green J., Haiman Z., Bryan G. L., 2011, Monthly Notices of the Royal Astronomical Society, 418, 838
  * Wollenberg et al. (2020) Wollenberg K. M. J., Glover S. C. O., Clark P. C., Klessen R. S., 2020, Monthly Notices of the Royal Astronomical Society, 494, 1871
  * Woods et al. (2017) Woods T. E., Heger A., Whalen D. J., Haemmerlé L., Klessen R. S., 2017, The Astrophysical Journal Letters, 842, L6
  * Woods et al. (2020) Woods T. E., Heger A., Haemmerlé L., 2020, Monthly Notices of the Royal Astronomical Society, 494, 2236
  * Yoshida et al. (2006) Yoshida N., Omukai K., Hernquist L., Abel T., 2006, The Astrophysical Journal, 652, 6
  * Yoshida et al. (2007) Yoshida N., Omukai K., Hernquist L., 2007, The Astrophysical Journal, 667, L117
  * Zwick et al. (2023) Zwick L., Mayer L., Haemmerlé L., Klessen R. S., 2023, Monthly Notices of the Royal Astronomical Society, 518, 2076