# The Formation of Very Massive Stars in Early Galaxies and Implications for
Intermediate Mass Black Holes

John A. Regan1,2,∗ John H. Wise3 Tyrone E. Woods4 Turlough P. Downes2 Brian W.
O’Shea5,6,7,8 Michael L. Norman9 1Department of Theoretical Physics, Maynooth
University, Maynooth, Ireland  2Centre for Astrophysics & Relativity, School
of Mathematical Sciences, Dublin City University, Glasnevin, D09 W6Y4, Ireland
3Center for Relativistic Astrophysics, Georgia Institute of Technology, 837
State Street, Atlanta, GA 30332, USA  4National Research Council of Canada,
Herzberg Astronomy & Astrophysics Research Centre, 5071 West Saanich Road,
Victoria, BC V9E 2E7, Canada  5Department of Computational Mathematics,
Science, and Engineering, Michigan State University, MI, 48823, USA
6Department of Physics and Astronomy, Michigan State University,MI, 48823, USA
7Joint Institute for Nuclear Astrophysics - Center for the Evolution of the
Elements, USA  8National Superconducting Cyclotron Laboratory, Michigan State,
University, MI, 48823, USA  9Center for Astrophysics and Space Sciences,
University of California, San Diego, 9500 Gilman Dr, La Jolla, CA 92093

###### Abstract

We investigate the ab-initio formation of super-massive stars in a pristine
atomic cooling halo. The halo is extracted from a larger self-consistent
parent simulation. The halo remains metal-free and star formation is
suppressed due to a combination of dynamical heating from mergers and a mild
($J_{\rm LW}\sim 2-10\ J_{21}$(z)) Lyman-Werner (LW) background. We find that
more than 20 very massive stars form with stellar masses greater than 1000
$\rm{M_{\odot}}$. The most massive star has a stellar mass of over 6000
$\rm{M_{\odot}}$. However, accretion onto all stars declines significantly
after the first $\sim$ 100 kyr of evolution as the surrounding material is
accreted and the turbulent nature of the gas causes the stars to move to lower
density regions. We post-process the impact of ionising radiation from the
stars and find that ionising radiation is not a limiting factor when
considering SMS formation and growth. Rather the birth environments are highly
turbulent and a steady accretion flow is not maintained within the timescale
(2 Myr) of our simulations. As the massive stars end their lives as direct
collapse black holes this will seed these embryonic haloes with a population
of black holes with masses between approximately 300 $\rm{M_{\odot}}~{}$and
10,000 $\rm{M_{\odot}}$. Afterwards they may sink to the centre of the haloes,
eventually coalescing to form larger intermediate mass black holes whose in-
situ mergers will be detectable by LISA.

###### keywords:

Early Universe, Supermassive Stars, Star Formation, First Galaxies, Numerical
Methods

††thanks<EMAIL_ADDRESS>Royal Society - SFI University Research
Fellow

##  1 Introduction

Supermassive stars (SMSs) with masses between $10^{4}$ and $10^{5}$
$\rm{M_{\odot}}~{}$have, over the past few decades, been invoked (Rees, 1978;
Begelman & Rees, 1978; Begelman, Volonteri & Rees, 2006; Begelman, Rossi &
Armitage, 2008; Latif, Schleicher & Hartwig, 2016; Woods et al., 2019) as an
intermediate phase to explain the existence of supermassive black holes
(SMBHs) at the centres of massive galaxies (Fan et al., 2006; Kormendy & Ho,
2013). The pathways to forming a SMBH are thus far unknown with a number of
theoretical models proposing to explain their existence.
Perhaps the simplest explanation is to invoke the black holes left over from
the first generation of stars as seeds for SMBHs. The first generation of
(metal-free) stars are referred to as Population III (PopIII) stars and
according to current theoretical models (e.g. Turk, Abel & O’Shea, 2009;
Clark, Glover & Klessen, 2008; Hirano et al., 2014; Stacy, Bromm & Lee, 2016)
the initial mass function should be top heavy with a characteristic mass of
tens of solar masses. However, PopIII remnant black holes are expected to form
in low density environments (Whalen, Abel & Norman, 2004; O’Shea et al., 2005;
Milosavljević, Couch & Bromm, 2009) and are not expected to accrete
substantially, at least not initially (Alvarez, Wise & Abel, 2009; Smith et
al., 2018). PopIII stars are therefore not seen as good candidates to explain
the existence of SMBHs without invoking super-Eddington accretion scenarios
which can boost their initial seed masses by an order of magnitude or more
over a short period (Lupi et al., 2014; Pacucci, Volonteri & Ferrara, 2015;
Sakurai, Inayoshi & Haiman, 2016; Inayoshi, Haiman & Ostriker, 2016a; Pacucci
et al., 2017; Inayoshi, Haiman & Ostriker, 2016b).
SMSs provide an alternative path to forming a SMBH by giving the seed black
hole a head-start compared to a black hole formed from a PopIII remnant. Under
a SMS formation scenario the accretion rate onto the protostar must exceed a
critical threshold thought to be around 0.001
$\rm{M_{\odot}~{}yr^{-1}}~{}$(Haemmerlé et al., 2018). When this threshold
accretion rate is reached and maintained the stellar radius inflates reducing
its surface temperature to approximately 5000 K and making the star resemble a
red giant star (Omukai & Palla, 2003; Hosokawa, Omukai & Yorke, 2013; Hosokawa
et al., 2013; Woods et al., 2017). However, the SMS must continue to accrete
above this threshold rate. If the accretion rate falls below the critical
rate, for a time exceeding the Kelvin-Helmholtz time (Sakurai et al., 2016),
the star contracts to the main sequence and becomes a hyper-luminous PopIII
star with a mass set approximately by the mass at which the accretion rate
dropped. When the accretion rate is maintained the star grows rapidly but
emits only weak radiative feedback with the spectrum of the emitted radiation
peaking below the hydrogen ionisation limit (Woods et al., 2019).
As discussed, the key requirement for forming a SMS is that the mass accretion
rate onto the star exceeds approximately 0.001 $\rm{M_{\odot}~{}yr^{-1}}$,
however, a sufficient baryon reservoir is also required and furthermore the
metallicity of the gas being accreted should be below $10^{-3}$
$\rm{Z_{\odot}}$. Detailed high resolution simulations have found that gas
that has been enriched above this threshold fragments into lower mass stars
which do not converge to a single object and in this case the formation of a
SMS is suppressed (e.g. Chon & Omukai, 2020). For these reasons metal-poor
(i.e., Z $\lesssim 10^{-3}$ $\rm{Z_{\odot}}$) atomic cooling haloes are seen
as the most promising candidates in which to form SMSs. Haloes which have
higher levels of metal-enrichment (Z $>10^{-3}$ $\rm{Z_{\odot}}$) may also be
viable candidates for SMS formation in the early universe if metal mixing is
inhomogeneous (Regan et al., 2020a). Atomic cooling haloes which provide the
above requirements for SMS formation were recently investigated by Wise et al.
(2019) and Regan et al. (2020b). In particular Wise et al. (2019) found that
the combination of a mild Lyman-Werner background combined with the impact of
dynamical heating effects due to minor and major mergers can suppress star
formation until a halo crosses the atomic cooling threshold. These haloes are
therefore predominantly metal-poor (with any metal enrichment coming
externally), have large baryon reservoirs and suppressed $\rm{H_{2}}$ content
due to the LW radiation fields. In this study we build on the previous works
cited above.
The goal of this study is to model the formation and evolution of
(super-)massive star formation in haloes that are exposed to moderate LW
backgrounds, which when combined with the effects of dynamical heating can
suppresses star formation below the atomic cooling limit. To pursue this
research we re-simulate two haloes from the original Renaissance simulations
using the zoom technique. We designate these haloes as HaloA and HaloB. Both
haloes were chosen as they exhibited near isothermal collapse of their inner
core in the original Renaissance datasets as shown by Regan et al. (2020b).
They were therefore identified as among the most promising candidates for SMS
formation. Both haloes were exposed to moderate levels of LW radiation from
nearby radiation sources as well as constant mergers which dynamically heated
the gas within the haloes. To further understand the impact of the LW field we
re-simulate HaloB with no LW radiation in this work. This is done to determine
if a halo can remain star-free due to only dynamical heating effects or if the
LW field remains a critical component. HaloA, on the other-hand, is re-
simulated with a LW background composed of both local source contributions and
background contributions.
In the zoom simulations, we find that HaloA forms stars with masses greater
than 6000 $\rm{M_{\odot}}~{}$but that the accretion rate onto individual
proto-stars always declines as the star’s immediate gas supply is depleted. In
the re-simulation of HaloB, without a LW field, we find that the halo
undergoes premature collapse (compared to the original case where a LW field
of JLW$\sim 2$ J21111J21 is shorthand for $1\times 10^{-21}\ \rm{erg\ cm^{-2}\
s^{-1}\ Hz^{-1}\ sr^{-1}}$ existed). In HaloB the most massive star in the
halo has a mass of approximately 173 $\rm{M_{\odot}}$. The impact of ionising
radiation is not considered in these simulations but post-processing of the
stellar feedback using Cloudy (Ferland et al., 2017) is instead used to gauge
the likely impact of ionising sources, particularly for HaloA which forms a
number of hyper-luminous PopIII stars.
The paper is laid out as follows: In §2 we very briefly review the original
Renaissance simulations as well as discussing the zoom-in simulations. In §3
we analyse the results of the zoom-in simulations. In §4 we discuss the
implications of the results and the connection with upcoming gravitational
wave observatories. In §5 we summarize our results and outline our
conclusions.

##  2 Methods

###  2.1 Renaissance Simulation Suite

Enzo has been extensively used to study the formation of structure in the
early universe (Abel, Bryan & Norman, 2002; O’Shea et al., 2005; Turk et al.,
2012; Wise et al., 2012, 2014; Regan, Johansson & Wise, 2015; Regan et al.,
2017). Enzo includes a ray tracing scheme to follow the propagation of
radiation from star formation and black hole formation (Wise & Abel, 2011) as
well as a detailed multi-species chemistry model that tracks the formation and
evolution of nine species (Anninos et al., 1997; Abel et al., 1997). In
particular the photo-dissociation of $\rm{H_{2}}$ is followed, which is a
critical ingredient for determining the formation of the first metal-free
stars (Abel, Bryan & Norman, 2000).
The original Renaissance simulations Xu, Wise & Norman (2013); Xu et al.
(2014); O’Shea et al. (2015) were carried out on the Blue Waters supercomputer
using the adaptive mesh refinement code Enzo (Bryan et al., 2014; Brummel-
Smith et al., 2019)222https://enzo-project.org/. The datasets that formed the
basis for this study were originally derived from a simulation of the universe
in a 40 Mpc on the side box using the WMAP7 best fit cosmology (Komatsu et
al., 2011). For more details on the Renaissance simulation suite see Chen et
al. (2014). Here we outline only the details relevant to this study for
brevity. The simulation suite was broken down into three separate regions,
namely the Rarepeak, Normal and Void regions. Each region was simulated with
an effective initial resolution of $4096^{3}$ grid cells and particles giving
a maximum dark matter particle mass resolution of $2.9\times 10^{4}$
$\rm{M_{\odot}}$. Further refinement was allowed throughout each region up to
a maximum refinement level of 12, which corresponded to 19 pc comoving spatial
resolution. Given that the regions focus on different overdensities each
region was evolved forward in time to different epochs. The Rarepeak region,
being the most overdense and hence the most computationally demanding at
earlier times, was run until $z=15$. The Normal region ran until $z=11.6$, and
the Void region ran until $z=8$. In all of the regions the halo mass function
was very well resolved down to M${}_{\rm halo}\sim 2\times 10^{6}$
$\rm{M_{\odot}}$.
As noted already in §1, in Wise et al. (2019) we examined two metal-free and
star-free haloes from the final output of the Rarepeak simulation and re-
simulated those two haloes at significantly higher resolution (maximum spatial
resolution, $\Delta x\sim 60$ au) until the point of collapse. This re-
simulation allowed us to investigate the evolution of the inner halo and the
mass distribution of the clumps formed. However, no star formation
prescription was employed during this re-simulation. In Regan et al. (2020b)
we subsequently investigated the occurrence of metal-free and star-free atomic
cooling haloes across all of the Renaissance datasets. We found a total of 79
such haloes in the Rarepeak outputs and three such haloes in the Normal
outputs. None were found in the Void outputs. Of the 79 haloes which were
metal-free and star free above the atomic cooling limit four haloes showed
almost ideal isothermal collapse (see Figure 6. in Regan et al. (2020b))
consistent with what has previously been identified as ideal conditions for
forming SMSs (Inayoshi, Omukai & Tasker, 2014; Becerra et al., 2015; Latif,
Schleicher & Hartwig, 2016; Regan et al., 2017; Chon, Hosokawa & Yoshida,
2018; Regan & Downes, 2018b). Of those haloes which collapsed isothermally we
then selected two haloes for re-simulation in this study.

Figure 1: Left Panel: The LW background field imposed on Halo A for the re-
simulation. HaloB is not shown as no LW background is imposed on HaloB. The LW
field values are composed of local sources, plus a LW background from Wise et
al. (2012), extracted from the original Renaissance simulations. The local
sources dominate over the background by at least an order of magnitude at all
times. In the original Renaissance simulations HaloA (and HaloB) remained
metal-free and star-free until $z=15.6$ ($z=15.0$). The redshift of original
collapse is marked with a green circle. A black star marks the redshift at
which star formation occurs in HaloA in the zoom-in simulations. Right Panel:
The merger history of both HaloA and HaloB as determined from the original
Renaissance simulations. Again we mark the redshift of first star formation in
HaloA and HaloB as found in the re-simulations with black stars. Note that
star formation now occurs significantly earlier in HaloB because no LW
background is imposed. In HaloA, on the other hand, star formation occurs 15
Myr after the halo was detected in the original (lower resolution)
simulations. Matom, orange line, is the atomic cooling threshold (Fernandez et
al., 2014) and Mmin,LW is the threshold mass for halo collapse in the presence
of an expected LW background at these redshifts (Machacek, Bryan & Abel, 2001;
O’Shea & Norman, 2008).

###  2.2 Re-simulation of pristine atomic cooling haloes

In order to make the simulation time tractable, we restricted the mesh
refinement to be focused around the target haloes. In order to do this we
first identified the Lagrangian volume (setting this to three times the virial
radius) of the target halo at the redshift at which it was first identified as
a star-free, atomic cooling halo with zero metallicity ($z=15.6$,
M${}_{halo}\sim 7\times 10^{7}$$\rm{M_{\odot}}~{}$for HaloA and $z=15.0$,
M${}_{halo}\sim 5\times 10^{7}$ $\rm{M_{\odot}}~{}$for HaloB) and tagged each
dark matter particle within this volume. This ensured that we captured the
dynamics of the gas and dark matter both within and surrounding the halo. Upon
restarting the simulation subsequently at $z=20$ we then set each of the
tagged dark matter particles as “must-refine-particles.” This meant that any
cell, in the simulation, containing one of these particles was allowed to
refine. Any cells not containing one of these particles could not refine. This
optimisation focuses the refinement solely onto the target halo (and gas
surrounding the halo). In addition to the “must-refine-particle” refinement
criteria, refinement is updated to be also based on the Jeans length of the
gas with the additional criteria that the Jeans length is always refined by at
least 64 cells.
Having now optimised the simulations to focus only on the target halo we next
looked to (dark matter) particle splitting. By splitting the dark matter
particles we increased the mass resolution of the simulation so as to match
the increased spatial resolution. Dark matter particles are split in Enzo
using the prescription of Kitsionas & Whitworth (2002) and was previously
described in Regan, Johansson & Wise (2015). By employing particle splitting
we reduced the dark matter particle mass to $M_{\rm DM}\sim 170$
$\rm{M_{\odot}}~{}$inside the target halo and its surrounding Lagrangian
volume.

###  2.3 The External Radiation Field

With refinement targeted only on a single halo (and its progenitors) star
formation in the surrounding galaxies is therefore neglected. In order to
account for the radiation that would otherwise be emitted by these galaxies we
extracted the LW emission from the original simulations and created a table of
LW values that this target halo is exposed to as a function of redshift. By
doing this the spatial information of the LW field is preserved and so the
target haloes feel the same LW radiation as they did in the original
simulations where star formation and feedback was followed in every galaxy in
the volume. Using the Grackle-2.1 (Smith et al., 2017) software library we
then used these LW tables as the “background” that this halo is exposed to.
Self-shielding to LW radiation is also invoked in the re-simulations based on
the Wolcott-Green, Haiman & Bryan (2011) prescription.
For the re-simulation of HaloA the LW radiation background follows exactly the
radiation field as determined from the original simulations and is shown in
Figure 1. For HaloB we purposely set the background radiation field to zero.
This allows us to study the evolution of HaloB with star formation suppression
due to dynamical heating only. HaloB is therefore the control
simulation333Ideally HaloA without a LW field would have been the control for
the simulation campaign. However, as noted in footnote 4 this was not possible
at this time. and allows us to examine the impact of what happens when no LW
field is present.

###  2.4 Subgrid Star Formation Prescription

In order to resolve star formation in the collapsing target haloes we set the
maximum refinement level of the simulation to 20. This is an increase of a
factor of $2^{8}(256)$ compared to the original Renaissance simulations and
allows us to reach a maximum spatial resolution of $\Delta x\sim 1000$ au.
While this (maximum) resolution is less than what was achieved in Wise et al.
(2019) it was necessary as the goal of this re-simulation was not only to
follow the collapse of the target halo but to also follow the formation of
stars within the collapsing halo for up to 2 Myr following the formation of
the first star. At the resolution used in Wise et al. (2019) this proved
intractable and so we reduced the resolution by a factor of $2^{4}(16)$,
compared to Wise et al. (2019), as a compromise. Reducing the refinement
factor compared to Wise et al. (2019) reduced the computational load while
still allowing us to resolve star formation at an acceptable resolution.
In order to model star formation within the collapsing gas cloud we employed a
star formation criteria using the methodology first described in Krumholz,
McKee & Klein (2004). The implementation in Enzo is described in detail in
Regan & Downes (2018a) and Regan & Downes (2018b) and we give a brief overview
here for completeness. Stars are formed when all of the following conditions
are met:

  1. 1.

The cell is at the highest refinement level

  2. 2.

The cell exceeds the Jeans density

  3. 3.

The flow around the cell is converging

  4. 4.

The cooling time of the cell is less than the freefall time

  5. 5.

The cell is at a local minimum of the gravitational potential

Once the star is formed accretion onto the star is determined by evaluating
the mass flux across a sphere with a radius of 4 cells centered on the star.
Initially all stars are assumed to be stars with low surface temperatures that
are appropriate for main sequence SMSs and less massive proto-stars on the
Hayashi track. The accretion onto the surface of the embryonic star is found
by applying Gauss’s divergence theorem to the volume integral of the accretion
zone (e.g Bleuler & Teyssier, 2014) (i.e. the volume integral of flux inside
the accretion zone)

$\dot{M}=4\pi\int_{\Omega}{\rho v_{r}^{-}r^{2}dr}$ (1)

where $\dot{M}$ is the mass accretion rate, $\Omega$ is the accretion zone
over which we integrate, $\rho$ is the density of the cells intersecting the
surface, $v_{r}^{-}$ is the velocity of cells intersecting the surface which
have negative radial velocities and $r$ is the radius of our surface. As noted
above we set the accretion radius to be 4 cells. The accretion onto the star
is calculated at each timestep, however this is likely to be a very noisy
metric. To alleviate this to some degree we average the accretion rate over
intervals of 1 kyr and use that averaged accretion rate in data outputs. The
accretion rate is added as an attribute to each star and hence a full
accretion history of every star is outputted as part of every snapshot.
Mergers with other stars are also included in the accretion onto the stars.
Simulations of massive star formation (e.g. Meyer et al., 2020) do show that
stellar mergers do increase the “bloatedness” of a star and hence including
the mergers in the accretion rates is likely to be valid. In this case the
more massive star retains its information (e.g., age, type, etc.) after the
merger event - information on the less massive star is lost. The mass of the
less massive star is added to the accretion rate of the more massive star for
that timestep. Stars are merged when they come within 3 times the accretion
radius of each other (i.e. 12 cell lengths). We also note here that the dark
matter particle mass is approximately $M_{\rm DM}\sim 170$ $\rm{M_{\odot}}$.
This mass is comparable to and greater than the mass of some of the stars that
are formed. While the star forming regions are strongly baryon dominated
individual dark matter particles may produce some shot noise. The shot noise
effects are accounted for however by smoothing the dark matter particles over
a scale of approximately 1 pc. This alleviates the impact of individual dark
matter particles on the central gas and star particles (e.g. Regan, Johansson
& Wise, 2015).

Figure 2: Left Panel: The percentage of the inner volume ionised in each halo
due to all stars as a function of time. In each case the volume is centred on
the most massive star at that time. Two radii are considered, a radius of 1 pc
from the most massive star is shown as a solid line, a radius of 5 pc is shown
as a dashed line. HaloA is shown in blue, HaloB in red. For both haloes the 1
pc volume surrounding the most massive star becomes ionised but as we look to
the larger region (5 pc) the medium remains neutral. Note that for HaloA the
most massive star changes over time and this is reflected in the ionisation
spikes. The first spike, where the volume becomes ionised quickly (at
approximately 700 kyr) is at least partially due to the star moving to a lower
density environment and hence ionisation becomes easier. For HaloB, after
approximately 1800 kyr the gas within 1 pc of the most massive star has become
fully ionised. However, at a radius of 5 pc the gas shows virtually no
ionisation. So while the stars are able to ionise their immediate surroundings
the radiative feedback does not penetrate much further than a pc. Right Panel:
The ratio of the ionising photon luminosity and the total integrated
recombination rate of the gas. For HaloA we see a similar trend to the left
panel. The most massive star, at a time less than $\sim 750$ kyr is quickly
able to ionise its surroundings. However, at larger radii (and also at later
times for HaloA when the most massive star is embedded in the central region)
the gas remains neutral. For HaloB the results are again similar to the volume
averaged calculation. At 1 pc the region becomes fully ionised after
approximately 1000 kyr but again at a radius 5 pc the region remains dominated
by recombinations.

Each star also has the ability to provide both radiative and mechanical
feedback, which is most appropriate in the case where the star has
transitioned into a (accreting) black hole. As the accretion rate onto the
star varies the star can transition its type from a SMS, with an inflated
surface, to a PopIII star. This transitioning only occurs if the accretion
rate onto the star either never exceeds the critical rate (set in our
simulations be be 0.04 $\rm{M_{\odot}~{}yr^{-1}}$) or if the accretion rate
onto the star falls below the critical accretion rate. While the star remains
bloated the radiative feedback from the star is primarily below the hydrogen
ionisation limit and is mostly in the form of infrared radiation. However, if
the accretion rate drops and the star contracts to the main sequence then its
surface temperature dramatically increases up to $10^{5}$ K, causing its
spectrum to harden and peak in the UV.
All stars in this simulation emit radiative feedback below the ionisation
threshold of hydrogen. The radiation is followed explicitly using the ray
tracing technique (Wise & Abel, 2011). Pop III stars are modelled assuming a
blackbody spectrum with a characteristic mass of 40 $\rm{M_{\odot}}$(Table 4,
Schaerer (2002)). From that we assign a LuminosityPerSolarMass to the Pop III
star and the star consequently becomes more luminous and the feedback more
intense as the mass of the star increases. SMSs are modelled by assuming a
blackbody spectrum with an effective temperature of Teff = 5500 K (Hosokawa et
al., 2013). The radiation spectrum for a SMS therefore peaks in the infrared
as opposed to the UV for Pop III stars. For the specific luminosity of the SMS
we take a characteristic mass of 500 $\rm{M_{\odot}}~{}$and apply the
contribution from the non-ionising photons only (Schaerer, 2002). As with the
’normal’ Pop III stars the SMS luminosity changes as mass is accreted and the
total luminosity then scales up as the mass increases.
In both cases the radiation from the stars is propagated outwards from the
star using the MORAY radiative transfer package (Wise & Abel, 2011) that is
part of Enzo. MORAY is able to model the ionisation of H, He and He+. It can
also account for the photo-dissociation of $\rm{H_{2}}$ for photons with
energies within the Lyman-Werner band and the photo-detachment of $\rm{H^{-}}$
and $\rm{H_{2}^{+}}$ for photons in the infrared band. For each type of star
we use five energy bins. The first two energy bins (E $<13.6$ eV) are weighted
by the cross section peaks for $\rm{H^{-}}$, $\rm{H_{2}^{+}}$ and $\rm{H_{2}}$
photo detachment/dissociation respectively. The next three energy bins are
determined using the sedop code developed by Mirocha et al. (2012) which
determines the optimum number of energy bins needed to accurately model
radiation with energy above the ionisation threshold of hydrogen. For the
self-shielding of $\rm{H_{2}}$ against LW radiation we use the prescription of
Wolcott-Green, Haiman & Bryan (2011).
As stars in this simulation contract to the main sequence, when the accretion
rate onto the star drops, radiation above the hydrogen ionisation threshold is
not initiated although in principle it should be and our radiative transfer
scheme does support ionising radiation. The problem however is that the shocks
generated near the star particle are too strong and unresolved for the PPM
reconstruction scheme and HLLC Riemann solver to handle, the steep gradients
cause the solution to overshoot to negative densities and energies (see Bryan
et al., 2014, for details). We could have avoided this problem by utilizing a
more diffusive solver at the expense of accuracy, but instead we neglect the
impact of ionising feedback and post-process the impact of ionising feedback
with Cloudy and use the results to estimate the extent to which the stars
(particularly the hyper luminous stars in HaloA) ionise their surroundings and
potentially shutdown further star formation444Initially our simulation plan
had been to run four haloes with and without a background LW field at a
refinement level of 24 (16 times higher than we report here). However, the
simulations proved intractable and computationally expensive and hence we
reduced the scope of the campaign. As a result, we did not have the capacity
to run HaloA without a LW radiation field or HaloB with a radiation field. We
are now in the process of undertaking a new more expansive campaign to address
these shortcomings..

Figure 3: Both panels show the projected number density6 of the region around
which the first star forms. HaloA is in the left hand panel, HaloB in the
right hand panel. The legend in each panel gives the mass of each star at the
first output time following star formation as well as the age of each star at
that time. The extent of each panel is 2 pc (physical). The orange star in
each case represents the most massive star. Stars coloured in red are stars in
which the accretion rate exceeds that required for supermassive star formation
(0.04 $\rm{M_{\odot}~{}yr^{-1}}$in our simulations). Blue stars are those
stars (normal PopIII stars) for which the rate is below the critical rate.
HaloA contains a single initial star. HaloB contains two stars (closely
separated) at the first output following star formation. For illustrative
purposes the size of each star in the projection is scaled as
$\rm{R_{star}}\propto\rm{M_{star}^{0.6}}$.

###  2.5 Post-Processing with Cloudy

Cloudy (Ferland et al., 2017) is a spectral synthesis code which models
radiative transfer through a gas, and its resulting thermal and chemical
equilibrium, under a wide range of conditions encompassing those expected for
interstellar matter. The code relies on a number of databases for computing
the behaviour of atoms and molecules, including tabulated recombination
coefficients obtained from Badnell et al. (2003) and Badnell (2006), with Case
A and B recombination predictions for single-electron systems from Storey &
Hummer (1995) and He I recombination rates from Porter et al. (2012), as well
as ionic emission data from the CHIANTI database (Dere et al., 1997; Landi et
al., 2012). Its chief limitation is that it is largely constrained to
modelling environments in 1-D and/or assuming spherical symmetry (though see,
e.g., Morisset, 2013; Fitzgerald et al., 2020, for recent efforts to extend
its implementation to pseudo-3D problems). For this reason, and because we do
not follow the stellar evolution of each formed star in detail, we must make a
number of simplifying approximations in using cloudy to estimate the impact of
ionising feedback.

For each halo, we consider the ionising feedback for a range of times between
1 kyr and 2000 kyr after the initial star formation has commenced. For each
snapshot, we divide the stars formed between those which are accreting
$>5\times 10^{-3}$ $\rm{M_{\odot}~{}yr^{-1}}$, which will evolve on the
Hayashi track due to $\rm{H}^{-}$ opacity in their inflated envelopes (e.g.,
Hosokawa et al., 2013), and those accreting below this threshold, which will
evolve blue-ward as they contract to become hot ionising sources on the ZAMS
(Haemmerlé et al., 2018). We assume that all rapidly-accreting (“red”) stars
remain negligible ionising sources, and that all slowly- or non-accreting
(“blue”) stars have thermally-relaxed to a main sequence temperature of
$\approx 10^{5}$K (Schaerer, 2002; Woods, Heger & Haemmerlé, 2020). We further
take their luminosities to be approximately Eddington (L $\approx 1.3\times
10^{38}\times(M/M_{\odot})$ erg s-1), and their spectra to be well-
approximated as blackbodies.

For each blue star, we then model the ionisation state of the surrounding gas
assuming spherical symmetry, with the density profile found from the average
gas density in successive shells 0.5pc in width, centred on each star in the
halo, and primordial abundance ratios taken to be 1.0:0.08232:1.6e-10:1e-16
for H:He:Li:Be (consistent with the results of the Planck Collaboration et
al., 2014, table 2, see cloudy documentation for further discussion). We
assume an inner boundary of $10^{15}$ cm ($\sim 3\times 10^{-4}$ pc) and
terminate our calculations at the outer boundary of each nebula once either
the gas temperature falls to 8000 K or the electron fraction falls below 5%.
Modelling the ionised nebula associated with each star in this way does not
account for the overlapping of Strömgren spheres associated with distinct
stars; we address this point and further limitations in $\S$3.2.

Figure 4: Both panels show the projected number density6 in each halo at the
end of each simulation. The left panel shows HaloA while the right panel shows
HaloB. The legend in each figure gives the mass of the five most massive stars
at the final output time, as well as the age of the star at that time. The
extent of each panel is 20 pc (physical). The orange star in each case
represents the most massive star. Blue stars are those stars (normal PopIII
stars) for which the rate is below the critical rate (0.04
$\rm{M_{\odot}~{}yr^{-1}}$) for SMS formation. At the times shown no SMSs
exist in either simulation because the accretion rate onto each star is less
than the critical rate. For illustrative purposes the size of each star in the
projection is scaled as $\rm{R_{star}}\propto\rm{M_{star}^{0.6}}$.

##  3 Results

###  3.1 Conditions for (Super-)Massive Star Formation

As noted in §1 the goal of this study is to model the formation of SMS
formation in haloes which are experiencing both LW feedback from nearby
galaxies and dynamical heating from mergers. Both of these effects suppress
the formation of regular PopIII stars in haloes below the atomic cooling limit
and hence may lead to the formation of a SMS in a larger halo (as was found in
Wise et al., 2019). In this study we include sub-grid star formation
prescriptions to follow the star formation process in two of these candidate
haloes.
In the left hand panel of Figure 1 we show the LW history that HaloA is
exposed to throughout its re-simulation. As discussed in §2 the LW field at
the location of the halo is determined from the original simulation, which
included star formation in all of the surrounding haloes. In the re-simulation
the adaptive mesh refinement grids are focused only around the target halo and
so the LW field must be imposed as a global background, albeit focused on a
single halo. As can be seen in the left hand panel of Figure 1 the LW field is
very flat and sits below JLW $\sim 2$ J21 until a redshift of $z\sim 16.5$ at
which point it increases strongly due to a nearby source which undergoes a
starburst and emits copious amounts of LW radiation (and ionising
radiation555The shorter mean free path of the ionising radiation means that
HaloA does not get photoionised.) in the direction of HaloA. The LW field
reaches its zenith at a redshift of $z\sim 15.8$ and thereafter starts to
decline. The green circle at $z=15.6$ represents the redshift at which the
halo was detected as a metal-free and star-free atomic cooling halo in the
original simulations (see Regan et al., 2020b, for details). The black star
indicates the redshift at which star formation begins in this re-simulation
($z=15.05$).
In the right hand panel we show the merger history of both HaloA and HaloB as
found in the original simulations. The mergers, both major and minor, drive
the dynamical heating which heats the gas and delays star formation (Wise et
al., 2019). HaloA experiences a major merger starting at $z\sim 16.2$ and
lasts until $z\sim 15.7$. This merger, in combination with the LW field,
suppresses star formation until z = 15.05 (again the green circles denotes the
time of collapse in the original, lower resolution, simulations). The merger
history of HaloB is shown as the red line. Recall that HaloB is run without
any LW background as a control case. As a result star formation is triggered
at $z=17.21$. Using the results from the original simulation which allows us
to see the future evolution of that halo we see that immediately after star
formation the halo undergoes a major merger. It is therefore highly likely
that, in the absence of a suppressing LW background, the merger triggers star
formation, which in the original simulation was suppressed until after $z=15$.
Therefore, we can see here that the absence of the background LW radiation
field was sufficient to allow PopIII star formation to take place. Dynamical
heating by itself was not sufficient for this halo to avoid star formation.
Star formation occurs in both HaloA and HaloB when they are at very different
phases in their evolution. In HaloA the mass of the halo is significantly
above the atomic cooling limit (M${}_{\rm{HaloA}}=9.3\times 10^{7}$
$\rm{M_{\odot}}$) and the halo has remained metal-free. HaloB, on the other
hand, collapses early with a mass of M${}_{\rm{HaloB}}=3.7\times 10^{6}$
$\rm{M_{\odot}}~{}$due to the lack of a LW background and thus can be
classified as a mini-halo. Before analysing the star formation in detail we
use results from Cloudy to determine if the negative radiative feedback from
stars can quench further star formation through photoionisation in each halo.

###  3.2 Determining the time of star formation quenching

To assess the impact of photoionising stellar feedback, we must investigate
the ionisation state of the gas in the innermost star-forming regions of each
halo. Here we adopt two different radii at which to calculate the ionisation
state of the gas. We use a radius of 1 pc to estimate the ionisation state of
the gas close to the most massive star and we also use a radius of 5 pc to
estimate the more “large scale” ionisation of the entire halo. We use the same
radii, at which to estimate the ionisation state of the gas, for both haloes
for consistency. We then evaluate the evolution of both the total ionised
volume and the total ionisation budget through each snapshot, in order to
provide subtly distinct but complementary measures of the strength of stellar
feedback.

Figure 5:  The mass function for the stars that formed in HaloA & HaloB over
approximately 2 Myr. Blue bars represent stellar masses from HaloA while red
bars represent stellar masses from HaloB. The dichotomy in masses is clearly
evident with a strong bias towards more massive stars in HaloA. This is due to
the higher temperatures in HaloA compared to HaloB which in turn leads to high
infall rates to the centre and hence more mass available for proto-stars to
accrete. The median mass for a single star in HaloA is 683
$\rm{M_{\odot}}~{}$and the median mass for HaloB is 44 $\rm{M_{\odot}}$.

Figure 6: We examine the radial profiles for both HaloA (blue line) and HaloB
(red line) immediately prior to star formation. The four panels are from
bottom left going clockwise: temperature, $\rm{H_{2}}$ fraction, enclosed gas
mass and infall rate. As HaloB collapsed early due to the lack of a LW
background (compared to the original Renaissance simulation) the temperature
of this halo is significantly lower in the centre. The central temperature
equilibrium value is close to 500 K which is characteristic of PopIII
simulations for the chemical network used in this work. This lower temperature
is due to the high $\rm{H_{2}}$ fraction, as shown in the upper left panel.
This can be contrasted with the radial profiles for HaloA which show
systematically higher temperatures in the core of the halo and lower
$\rm{H_{2}}$ fractions. The lower $\rm{H_{2}}$ fractions are due to a
combination of dynamical heating from mergers and from local radiation
sources. The enclosed mass fractions identify the difference in masses between
the haloes. In the case of HaloB the enclosed masses are systematically lower
compared to HaloA as this is a effectively a minihalo. Infall rates in both
haloes are very different in the centre of each halo. Average infall rates in
the very centre (R $\lesssim 0.1$ pc) of HaloA are above 0.1
$\rm{M_{\odot}~{}yr^{-1}}~{}$compared to below 0.001
$\rm{M_{\odot}~{}yr^{-1}}~{}$in HaloB. These infall rates are consistent with
the accretion rates that are seen onto the stars that subsequently form in
each halo. The flat patches in the infall rate are due to outflow at those
radii (which show up as negative infall rates).

First, we compute the ionised volume associated with each thermally-relaxed
star as described in §2.5, and integrate the intersection of the union of
these Strömgren regions within the central star-forming regions, as defined in
the preceding paragraph for each halo, using a Monte Carlo approach. Note that
for a constant combined ionising luminosity, any reduction in ionised volume
in accounting for overlaps in our spherically-symmetric Strömgren spheres
would in principle be compensated by the ionisation of a greater (though
presumably not spherically symmetric) volume by the remaining available
ionising photons. The extent to which this would raise the total ionised
volume beyond our estimate depends sensitively on the local density
distribution beyond each Strömgren region.

Indeed, an alternative approach is to simply compare the ionising photon
luminosity of all stars within the inner star-forming region with the total
integrated recombination rate $\dot{\mathcal{N}}_{\rm R}$ of the gas if all
hydrogen atoms were ionised; the latter sets the total emission rate of
photons with $E>13.6$ eV needed to maintain ionisation of the star-forming
region in equilibrium. This can be found from integrating over the hydrogen
number density:

$\dot{\mathcal{N}}_{\rm
R}=\int_{\rm{V}}\alpha_{\rm{B}}(\rm{H}^{0},T)n_{\rm{H}}^{2}(r)d\rm{V}$ (2)

where $n_{\rm{H}}(r,t)$ is the spherically-averaged density profile centered
on the most massive star in the Halo at time $t$ (see above), and
$\alpha_{\rm{B}}(\rm{H}^{0},T)$ is the case B recombination coefficient for
hydrogen, which we take as $\approx 1.4\times 10^{-13}\rm{cm}^{3}\rm{s}^{-1}$
for primordial gas at a temperature typical of plasmas photoionised by
Population III stars, $\rm{T}\sim 2\times 10^{4}$K (Osterbrock & Ferland,
2006; Johnson et al., 2012).

These two measures of the strength of ionising stellar feedback are compared
in Figure 2. The left panel shows the volume ionised fraction of the innermost
region as a function of time, while the right panel shows the fractional
ionising budget as a function of time. The ionisation state of the gas is
shown in blue for HaloA and red for HaloB. Calculations with a 1 pc sphere are
shown as solid line, calculations with a 5 pc radius are shown as dashed
lines. Note that the ionised volume (left panel) and the fractional ionising
photon budget (right panel) do not necessarily move perfectly in tandem. This
is partly due to our simplified treatment as described above, including
overlapping Strömgren regions and spherically-averaged density profiles, but
also reflects the variable density distribution, with the lowest density
regions preferentially ionised over high-density knots.
We begin by focusing on HaloA. In the left hand panel, we see a marked rise in
the ionised volume between $\sim 500$ and $\sim 1200$ kyr when measured in a
sphere of 1 pc, this is matched in both cases by a similar rise in the
fractional ionising photon budget (right panel). In HaloA we see a sharp rise
matched by a sharp drop (solid blue line). This is because the most massive
star in the simulation switches at approximately 1000 kyr as another star’s
accretion rate gives it a larger mass. Prior to this time the most massive
star easily ionises its immediate surrounding and this is at least partly
because the star also moves into a lower density region (something similar
occurs in HaloB).
After 1000 kyr the most massive star in HaloA becomes embedded in dense
filaments and its ability to ionise even its immediate surroundings is blunted
(and hence the ionisation fraction hovers between approximately 0.1 and 0.4).
The right hand panel shows a similar result for HaloA with the most massive
star, prior to 1000 kyr, easily able to ionise its immediate surroundings.
However, after this point, with the most massive star more centrally located,
recombinations dominate. When considering the larger 5 pc radius (dashed blue
line) ionisation becomes much more difficult. In the left hand panel we see
that the 5 pc region is approximately 40% ionised by 2000 kyr but if we look
at the right hand panel we see that recombinations are still overwhelmingly
dominant. This is because HaloA contains a plethora of dense knots and
filaments in which recombinations are dominant and into which ionising photons
cannot penetrate.
For HaloB both measures, at 1 pc, show a fully ionised medium after
approximately 1000 kyr. The most massive star has a mass of 173
$\rm{M_{\odot}}~{}$and crucially it moves to a lower density region away from
the central over-density. However, like HaloA, using a radius of 5 pc shows a
medium whose ionisation state remains low reflecting the fact that while the
stars are very massive and easily able to ionise their immediate surroundings
they are not yet numerous enough to ionise the entire halo.
In summary, ionising stellar emission is unlikely to shut down star formation
in either halo prior to the onset of chemical and radiative emission from the
supernovae explosions expected to occur for some of the stars (e.g. see Heger
et al., 2003) approximately 2 Myr after their formation. In the case of HaloA
this is particularly important. HaloA contains a number of hyper-luminous
PopIII stars, however, the much denser gas in the central region of HaloA is
robust against stellar feedback allowing star formation to proceed in these
dense pockets, albeit with the overall volume ionisation growing in importance
as the simulation progresses (see left panel of Figure 2). The results from
our Cloudy models have very recently been confirmed by the 1-D analysis of
Sakurai, Haiman & Inayoshi (2020) who investigated the impact of ionising
feedback from a growing protostar in a very similar environment. They found,
similar to us, that no significant HII region forms during the evolution of
the star.

Figure 7:  Left Panel: The accretion rate onto all stars in HaloA. The mass
accretion onto the most massive star is shown by the blue line. The dashed
blue line gives the accretion rate onto the median mass star in HaloA.
Accretion onto all others stars are shown as grey lines. In all cases initial
accretion rates decline over time, dropping below the critical rate for SMS
formation after, on average, 100 kyr. The black dashed line gives the critical
rate for SMS star formation (Haemmerlé et al., 2018). Right Panel: The mass
evolution of every star HaloA. The solid blue line is for the most massive
star in HaloA. It has a mass of 6127 $\rm{M_{\odot}}~{}$after 143 kyr. We also
look at the mass evolution of an average star in HaloA for comparison - given
by the dashed blue line. In HaloA a median4 star has a mass of 683
$\rm{M_{\odot}}$.

Figure 8:  The same as Figure 7 except this time for the evolution of stars in
HaloB. The most massive star in HaloB has a mass of 173 $\rm{M_{\odot}}$. Note
that the sharp spike observed at approximately 200 kyr in the left hand panel
is caused by a stellar collision.

###  3.3 The onset and end of star formation

Star formation is triggered when all of the criteria set out in §2.4 are
fulfilled. In Figure 3 we show a projection of the total gas number
density666Number density refers here to the total number density from each
separate species weighted by each species atomic mass unit. in the regions
surrounding the formation of the first star in each halo. The projections are
made from the first snapshot following star formation. Each panel is 2 pc
across. In the left hand panel we show the projection from HaloA, which shows
the formation of a single star (coloured in orange at the centre of the green
coloured gas cloud). The legend on the top left gives the mass and age of the
star at this time. In the right-hand panel we show the projection for HaloB.
In this case two stars are formed by the time of the first output following
star formation. The most massive star has a mass of 37 $\rm{M_{\odot}}~{}$and
is 28 kyr old, the second star has a mass of 6 $\rm{M_{\odot}}~{}$and is 4 kyr
old. Stars are coloured in red if the accretion rate exceeds the SMS critical
rate of 0.04 $\rm{M_{\odot}~{}yr^{-1}}$, otherwise stars are coloured in blue
(denoting PopIII stars). The most massive star is coloured in orange.
In Figure 4 we show the extent of star formation by the end of each
simulation. In the left panel we show the results for HaloA and for HaloB in
the right hand panel. Both simulations were terminated approximately 2 Myr
after the formation of the first star. In both cased we terminate the
simulations before the imminent impact of chemical and mechanical feedback
from the first supernovae. At the end of the HaloA simulation there are 99
stars in the simulation with masses ranging from approximately 40
$\rm{M_{\odot}}~{}$up to a maximum mass of over 6000 $\rm{M_{\odot}}$. Star
formation in HaloA is widespread throughout the inner 20 pc of the halo with a
number of different gas clouds giving rise to star formation. Furthermore, as
the simulation develops the interactions between clouds triggers star
formation as individual clouds merge and interact. HaloB, on the other hand,
contains essentially only a single site of star formation due to the
significantly smaller halo mass. At the end of the simulation of HaloB there
are 21 stars with masses ranging from approximately 20 $\rm{M_{\odot}}~{}$up
to approximately 170 $\rm{M_{\odot}}$.
In Figure 5 we show the mass distribution of stars at the final output time.
Stars from HaloA are binned in blue, stars from HaloB are binned in red. The
median mass of stars formed in HaloB is 44 $\rm{M_{\odot}}~{}$while the median
mass of stars formed in HaloA is 683 $\rm{M_{\odot}}$. The most massive star
in HaloA at the final output time is 6127 $\rm{M_{\odot}}~{}$with another star
having a mass of 4477 $\rm{M_{\odot}}$. HaloA has 22 stars with masses
exceeding 1000 $\rm{M_{\odot}}$. The lowest mass star in HaloA has a mass of
40 $\rm{M_{\odot}}$. The mass distribution of stars in HaloB is significantly
smaller, running from 22 $\rm{M_{\odot}}~{}$up to 173 $\rm{M_{\odot}}$.
However, it should be noted that due to our finite resolution ($\Delta x\sim
1000$ au) we cannot accurately probe the lower end of the initial mass
function of either halo and we are likely missing some lower mass stars.

Figure 9: The enclosed gas mass (blue line), enclosed stellar mass (orange
line) and the number of stars (dashed black line) as a function of radius. The
radial profile is centred on the centre of mass of stars in the simulation.
The profile is created from the final dataset of HaloA. The stars are
concentrated in the central 20 pc of the simulation with a total stellar mass
of approximately 90,000 $\rm{M_{\odot}}$. The baryonic mass at 1 pc is
approximately $3\times 10^{4}$ $\rm{M_{\odot}}~{}$rising to approximately
$6\times 10^{5}$ $\rm{M_{\odot}}~{}$at 20 pc.

Figure 10: The evolution of the most massive star in HaloA and four other
stars which cross regions of high density during their evolution resulting in
enhanced accretion rates. The accretion rate for each star (most massive in
green) is shown in the top left panel. The final mass of each star ranges from
approximately 1135 $\rm{M_{\odot}}~{}$up to 6127 $\rm{M_{\odot}}$. Each of the
stars (with the exception of the most massive star) shows a second spike of
$\dot{M}\gtrsim 10^{-}3$ $\rm{M_{\odot}~{}yr^{-1}}~{}$during its evolution.
These times of high accretion are highlighted using circular markers with
embedded black dots. For most stars this happens a few hundred kyr after
formation. This spike in accretion is caused by the star crossing a region of
high density. The gas number density surrounding each star is shown in the top
right panel. The relative velocity of the stars (which also impacts accretion)
is shown in the bottom right panel. Note that the relative velocity values can
be as high as 10 - 15 km $\rm{s^{-1}}$during a high accretion event as the
star can be crossing through a dense clump. Finally, in the bottom left panel
we plot the distance from each star to the nearest high density clump. Times
of high accretion always corresponds to a distance of 0 pc from a clump.
However, as stars cross the inner halo their distances from a high density
clump varies in general. Note that for the most massive star, shown in green,
the stellar age is quite short as this star formed near the end of the
simulation.

###  3.4 The case for super-massive star formation

In Figure 4 we saw that the mass of the most massive star in HaloA was 6127
$\rm{M_{\odot}}$. While this mass is well beyond the mass of ordinary PopIII
stars (Turk, Abel & O’Shea, 2009; Greif et al., 2011; Wise et al., 2012;
Crosby et al., 2013; Susa, Hasegawa & Tominaga, 2014; Hirano et al., 2014;
Stacy, Bromm & Lee, 2016) and also more than 30 times larger than the most
massive star in HaloB it is still well short of the mass often associated with
truly super-massive stars (e.g. Woods et al., 2019) which are expected to have
end stage masses of $\sim 10^{5}$ $\rm{M_{\odot}}$. Furthermore, HaloA was
chosen here to present near ideal initial conditions in which to form a
supermassive star. To illustrate this point we show radial profiles of the
halo properties of both HaloA and HaloB in Figure 6. Figure 6 shows the state
of the gas in each halo immediately prior to star formation. HaloA is denoted
by the blue line, HaloB by the red line. The temperature, $\rm{H_{2}}$
fraction, enclosed gas mass and mass inflow rate are shown in clockwise order
starting from the bottom left panel.
HaloA shows near isothermal collapse, although the gas does show some degree
of cooling at approximately 10 pc. The $\rm{H_{2}}$ fraction is very different
between the two haloes, with HaloA showing a steep decline in $\rm{H_{2}}$
towards the centre due to the combination of the LW background and the rapid
assembly of mass. HaloB on the other-hand shows a more typical $\rm{H_{2}}$
evolution consistent with mini-halo formation. The enclosed mass plot in the
top right panel illustrates the different mass associated with each halo as a
function of radius. HaloA being well inside the atomic cooling mass range. In
the bottom right hand panel we see the mass inflow rate for each halo. HaloA
shows mass infall rates averaging greater than 0.1
$\rm{M_{\odot}~{}yr^{-1}}~{}$all the way into the centre of the halo. HaloB’s
mass infall rates on the other hand fall from approximately 0.1
$\rm{M_{\odot}~{}yr^{-1}}~{}$at 1 pc down to approximately $10^{-3}$
$\rm{M_{\odot}~{}yr^{-1}}~{}$in the very centre. As we will see this result is
reflected in the initial accretion rates observed onto the protostars.
In order to understand why, given apparent high mass infall rates that are
greater than the critical rate, a SMS does not form we need to examine the
accretion rate that is measured onto the stars within each simulation itself.
In Figure 7 we show the mass evolution and the accretion rate onto every star
in HaloA, the most massive star in HaloA is coloured in blue. In the left
panel we show the mass accretion rate while in the right hand panel we show
the mass evolution. In each case we plot as a function of the stellar age. The
dashed line in the panel shows the mass evolution of a star with the
median777We define a median star here as the star at the end of the simulation
which has a final mass closest to the median stellar mass at the end stellar
mass in the simulation. The most massive star in HaloA is formed when two gas
clouds within the central region merge and trigger star formation. The star
quickly grows in mass up to its final mass of 6127 $\rm{M_{\odot}}~{}$within
approximately 20 kyr, after which its accretion rate declines. The shaded grey
lines are the evolutionary tracks for each of the other 100 stars in HaloA.
The trend is similar for every other star.
Initially each star follows approximately the mass accretion rate as found in
the radial profiles (see Figure 6). However, in every case after a few tens of
kyr the accretion rate onto each star drops dramatically and mass growth is
halted. This drop in accretion is not due to feedback since no ionising
feedback is present in these simulations. Instead mass accretion is terminated
as the star accretes all the gas in its surroundings and moves into a less
dense region of the inner halo. All of the stars formed have a small initial
velocity of a few km $\rm{s^{-1}}$, which it takes from the gas that initially
formed the star. While this initial velocity is much smaller than the circular
velocity of the host halo it does cause the stars to move around and decouple
from the high density gas cloud in which they initially formed.
In Figure 8 we show the mass evolution and accretion rate onto each star in
HaloB. Again in the left panel we show the mass accretion rate while in the
right hand panel we show the mass evolution. The mass evolution of the most
massive star (coloured in blue) in HaloB is a little different. In this case
accretion is somewhat lessened at approximately $10^{-3}$
$\rm{M_{\odot}~{}yr^{-1}}~{}$(left panel). The most massive star in HaloB
results from a stellar merger which then drives the star to increase its mass
suddenly. This occurs when the main progenitor star is approximately 200 kyr
old and results in a star of 173 $\rm{M_{\odot}}$. However, the general trends
are similar with initial accretion rates always declining as stars move into
less dense environments.

###  3.5 Stellar Multiplicity

In Figure 9 we plot the enclosed gas (blue line) and stellar mass (orange
line) of HaloA as a function of radius. The radial profiles are centred on the
centre of mass of the stars in the halo. We also show the number of stars
(black dashed line) in each radial bin up to the number of stars in the halo
(101). All of the stars are concentrated in a radius of 20 pc of the halo
centre. By the end of the simulation the total stellar mass in HaloA is 90,000
$\rm{M_{\odot}}~{}$(cf. 1300 $\rm{M_{\odot}}~{}$in HaloB). The star formation
efficiency in the inner 20 pc (5 pc) of HaloA (HaloB) is $\sim$ 23% ($\sim
2$%). The stellar to dark matter ratios are higher in each case since this
inner regime is baryon dominated. 20 pc is the approximate Jeans length of the
gas in HaloA ($\sim 5$ pc for HaloB). In order for the formation of a truly
supermassive star in HaloA all of the stellar mass (i.e. $\sim$90,000
$\rm{M_{\odot}}$) would have needed to have been accreted onto the stellar
surface. Instead what happened was that the outer gas cloud in HaloA
fragmented into a small number of sub-clouds which each generated hotspots of
star formation. These sub-clouds tidally disrupt each other during the
subsequent evolution both triggering further star formation but also tidally
disrupting accretion. From a total baryonic reservoir of approximately 600,000
$\rm{M_{\odot}}~{}$inside 20 pc, approximately 1% of the baryons went into the
most massive star. To form a SMS with a mass of close to 100,000
$\rm{M_{\odot}}~{}$we would need over a fifth of the mass to flow into a
single object. While this is possible it is clearly going to be very
challenging given the turbulent nature of the environment in which these stars
are forming.
The stars in the central 20 pc region of HaloA wander the inner parts of the
halo in a random walk fashion. As stars cross regions (filaments and clumps)
of high density they can experience enhanced accretion rates. In Figure 10 we
show the environmental changes which occur for five stars within HaloA. Four
of the stars were chosen because at some point in their evolution they undergo
a second period of high accretion ($\dot{\rm{M}}>1\times 10^{-3}$
$\rm{M_{\odot}~{}yr^{-1}}$) at least 100 kyr after their formation. We also
include the most massive star, coloured in green. Periods of high accretion
are identified by a circular marker with an embedded black dot. Note that each
of the stars (excluding the most massive) undergoes a “bump” in its accretion
rate independent of the initial high accretion rate value (top left panel of
Figure 10). The high accretion rate is directly correlated with the star
passing through a region of high density which is shown in the top right
panel. The relative velocities of each of the stars are shown in the bottom
right panel while the distance from each star to the nearest high density
clump is shown in the bottom left panel. Average relative velocities are below
10 km $\rm{s^{-1}}$ with some spikes showing velocities between 10 and 20 km
$\rm{s^{-1}}$. A high density clump is defined as any cell that has a density
more than 10% of the Jeans Density. The distance of a star from a high density
clump is highly variable with stars encountering a high density clump only
occasionally. Encounters with filaments and clumps do occur however with an
associated spike in the accretion rate (and density) around the star. What is
also clear is that periods of high accretion are fleeting, lasting only a few
kyr at most. Apart from these short periods of rapid growth the initial mass
of the stars does not increase appreciably and a consistently high accretion
rate is not maintained. This is similar in spirit to what Smith et al. (2018)
found for accreting stellar mass black holes in similar high-z environments.
In summary none of the stars in our simulation are able to sustain critical
accretion rates for SMS formation. Nonetheless, the most massive star(s) in
the halo will form of population of intermediate mass black holes (IMBHs) with
masses ranging from 300 $\rm{M_{\odot}}~{}$to several thousand
$\rm{M_{\odot}}$. The subsequent growth of these IMBHs within the halo then
through mergers and accretion will then slowly build the black hole mass. The
key issue preventing the formation of a monolithic SMS appears to be that
these haloes are highly turbulent and that sustaining accretion onto a single
object in such an environment is hugely challenging (see Chon, Hosokawa &
Yoshida (2018) and §5 for a more detailed discussion). This conclusion is
unlikely to be due to insufficient resolution: it is well-known that low
resolution damps turbulent motion very significantly (e.g. Federrath et al.,
2010; Downes, 2012). Thus insufficient resolution will lead to an over-
estimate of the likelihood of forming an SMS.

Figure 11: The collisional time and relaxation times for the stars in HaloA.
The collisional times is an estimate of the time taken for a collision between
two stellar populations. For the stars in HaloA the time taken for a single
collision to occur exceeds the Hubble time at all radii. The relaxation time
is a measure of the time taken for the stellar system to contract to higher
density, perhaps driving a higher collision rate, due to two body relaxation.
For HaloA the relaxation time is approximately 10 Myr at 1 pc.

##  4 Intermediate Mass Black Hole Merging and Implications for Gravitational
Wave detection

We begin this section by analysing the stellar content of HaloA and in
particular the likelihood that HaloA will undergo two-body collisional
processes akin to the runaway collapse of a dense stellar cluster (Portegies
Zwart et al., 2004; Gürkan, Freitag & Rasio, 2004; Freitag, Gürkan & Rasio,
2006; Regan & Haehnelt, 2009; Katz, Sijacki & Haehnelt, 2015). In Figure 11 we
plot both the collisional time for the stars in HaloA and the relaxation time
of the stellar system in HaloA as a function radius. The collisional time,
tcoll, for a stellar population consisting of two populations can be written
as (Freitag, 2008):

$t_{\rm{coll}}=5\ \rm{Gyr}\Big{(}{10^{6}\over
n_{*}}\Big{)}\Big{(}{\sigma_{vel}\over 10\ \rm{km\
s^{-1}}}\Big{)}\Big{(}{2R_{\odot}\over
R_{1}+R_{2}}\Big{)}\Big{(}{2M_{\odot}\over M_{1}+M_{2}}\Big{)}$ (3)

where $n_{*}$ is the stellar density, $\sigma_{vel}$ is the one dimensional
velocity dispersion, R1 is the stellar radius of the first population of stars
and R2 is the stellar radius of the second population. We choose the median
stellar mass, 683 $\rm{M_{\odot}}$, for M1 and we choose our maximum stellar
mass, 6127 $\rm{M_{\odot}}$, for M2. The stellar radii are calculated using
the standard formulae from Stahler, Palla & Salpeter (1986). The collisional
time, plotted in blue, is a measure of time required for one stellar collision
between these two populations. In orange we plot the relaxation time of the
cluster. The relaxation time is the time required for the system to contract
due to two body interactions in a populous cluster (N $\gg$ 10). The
relaxation time is given by:

$t_{\rm{relax}}=2\ \rm{Myr}\Big{(}{10\over\lambda}\Big{)}\Big{(}{10^{6}\over
n_{*}}\Big{)}\Big{(}{\sigma_{vel}\over 10\ \rm{km\
s^{-1}}}\Big{)}^{3}\Big{(}{1M_{\odot}\over\langle m\rangle}\Big{)}^{2}$ (4)

where $\lambda=\ln(0.02N)$ and N is the number of cluster stars and $\langle
m\rangle$ is the median stellar mass. The collisional time, tcoll, is greater
than the Hubble time at all radii indicating that collisions are very unlikely
for the stars in HaloA. The relaxation time, trelax, is approximately 10 Myr
at 1 pc. However, this time exceeds the lifetime of the stars in our cluster
and contraction to very high densities appears unlikely and hence the runaway
collapse of the cluster is not predicted (or observed here).
At the end of the simulation HaloA has 22 stars with masses greater than 1000
$\rm{M_{\odot}}~{}$with the most massive star having a mass of 6127
$\rm{M_{\odot}}$. Within the central 1 pc of the halo centre there is a total
(stellar + gas) mass of $\sim 3\times 10^{4}$ $\rm{M_{\odot}}$. This puts
HaloA at the lower end of what defines a Nuclear Star Cluster (NSC) (e.g.
Georgiev et al., 2016). However, subsequent, further rapid growth of this halo
beyond this initial phase of star formation will drive the enclosed mass
further upwards. This would have the effect of pushing the environment to a
point where the central region of a galaxy becomes a proto-NSC. At this stage
if the density of the core reaches a point where the enclosed total baryonic
mass within the central parsec reaches approximately $10^{6}$
$\rm{M_{\odot}}~{}$then core collapse can occur (Davies, Miller & Bellovary,
2011) potentially reaching a central black hole mass of M${}_{BH}\sim 10^{5}$
$\rm{M_{\odot}}$. Depending on the densities reached core collapse could occur
on timescales of a few hundred Myr (Davies, Miller & Bellovary, 2011) although
equally the NSC could act as an incubator for the rapid growth of a seed black
hole through gas accretion achieving similar masses over a potentially longer
timescale (Natarajan, 2020).
Our specific case here shows that we will have a proto-galaxy that gets
populated with a plethora of weakly accreting black holes early in the proto-
galactic evolution. It is therefore reasonable to ask what are the prospects
for the merger of this first generation of black holes into a population of
more massive black holes?.
Pfister et al. (2019) investigated the merging of black holes in high
resolution simulations ($\Delta x\sim 10$ pc) with black hole masses in the
range $10^{4}$ $\rm{M_{\odot}}~{}$to $10^{5}$ $\rm{M_{\odot}}~{}$in high-z
galaxies. They found that dynamical friction from stars has a more stabilizing
effect on the black holes than the gas component but that the stellar
distribution in high-z galaxies is highly irregular and hence is not able to
stabilise the black hole orbits. The irregular stellar distributions found by
Pfister et al. (2019) correlates with what we see in our simulations here.
Pfister et al. (2019) found that black holes with masses less than $10^{5}$
$\rm{M_{\odot}}~{}$do not sink to the centre of a galaxy, but instead exhibit
random walk characteristics. However, their simulation neglected the impact of
any radiative feedback from the black hole (which could impede the sinking
timescale) and also lacked sufficient resolution to detect high density gas
which could promote dynamical friction.
Toyouchi et al. (2020), building on previous investigations of Park, Wise &
Bogdanović (2017) and Park & Bogdanović (2019), investigated a more idealised
scenario where they modelled the dynamics of a single $10^{4}$
$\rm{M_{\odot}}~{}$mass black hole incorporating the effects of dynamical
friction from gas drag and radiative feedback. They found that if the gas
density is low that the black hole does not sink to the centre - similar to
what Pfister et al. (2019) found. However, in the case where the black hole
encounters dense gas ($n_{gas}\gtrsim 10^{6}\ \rm{cm}^{-3}$) the ram pressure
of the head wind (due to radiative feedback) causes the black hole to lose
orbital energy and so the black hole can sink to the centre on a timescale of
t${}_{dyn}\sim$ 0.01 Myr, consistent with the dynamical friction timescale in
the Bondi-Hoyle-Littleton case (i.e. no radiative feedback), and much shorter
than the dynamical timescale in galactic gas disks.
Taking the two results together this appears to imply that the specific
environment that the black hole finds itself in will play a central role in
determining the timescale for the black hole to sink towards the centre. If
the black hole can interact with sufficiently dense gas it may sink towards
the centre (even for black holes with masses of $\sim 10^{4}$
$\rm{M_{\odot}}$) while in the case where the black hole does not encounter
sufficiently dense gas then the black hole may wander the central region of
the galaxy for many dynamical times. The galactic centres in which our massive
PopIII stars form (especially HaloA) contain a web of dense gaseous filaments
(with gas clearly approaching that critical level of $n_{gas}\gtrsim 10^{6}\
\rm{cm}^{-3}$ (see Figure 4) in the centre of the halo) which may help to
supply the necessary ram pressure and extract orbital energy. A further build
up of dense gas would potentially help to promote efficient frictional forces
even allowing for radiative feedback.

#### What then does this mean for potential gravitational wave detection?

Within approximately 4 Myr after the formation of the first massive stars here
there will be a population of black hole seeds with masses in the range 300
$\rm{M_{\odot}}~{}$(upper limit of pair instability supernovae) up to 6000
$\rm{M_{\odot}}~{}$(and likely higher). In Figure 12 we plot the
characteristic strain expected from the merger of two black holes with masses
of 5000 $\rm{M_{\odot}}~{}$each at $z=15$ we note that here that gravitational
recoil effects can eject black holes from galactic centres (e.g. Gültekin,
Miller & Hamilton, 2006), however we consider the idealised case of a non-
spinning binary black hole merger with a mass ratio, $q=1$ and hence in this
idealised case ejections can be avoided (Holley-Bockelmann et al., 2008;
Morawski et al., 2018; Dunn, Holley-Bockelmann & Bellovary, 2020). These
masses are consistent with the masses of the stars found in this simulation
and black holes with similar masses are expected from the direct collapse of
these massive PopIII stars (Heger et al., 2003). The merger of two black holes
with these masses will produce a gravitational wave signal detectable by LISA
(eLISA Consortium et al., 2013; Sesana, 2016; Cornish & Shuman, 2020) emitting
gravitational waves between approximately $10^{-4}$ and $10^{-1}$ Hz (detector
frame frequencies). The Signal-to-Noise ratio (SNR) from the merger of two
black holes with masses of 5000 $\rm{M_{\odot}}~{}$at $z=15$ reaches a SNR
$\sim 32$ over its time inside the LISA band and hence is well within LISA’s
detection parameters888The SNR was calculated using the LISA sensitivity
calculator available at
github.com/eXtremeGravityInstitute/LISA_Sensitivity.git (Robson, Cornish &
Liu, 2019). The merger of a 10,000 $\rm{M_{\odot}}~{}$binary black hole with a
mass ratio of 1 at a redshift of 15 will enter the LISA band approximately 2
months before merger and complete thousands of orbits before the final plunge.
The SNR value is calculated by integrating over the entire frequency range
over which the binary overlaps with the LISA band. The SNR exceeds 10
approximately 12 hours before merger. With this number of cycles LISA will be
able to detect the redshifted mass with a precision of close to 1% (Sesana,
2013), with the strongest signal and hence the highest SNR being achieved in
the final few hours and hundreds of orbits (Robson, Cornish & Liu, 2019). The
largest uncertainty will however come from the determination of the redshift
(or luminosity distance) which at this redshift is expected to have an
uncertainty of approximately 30% (Sesana, 2013). Additional uncertainties due
to weak lensing will also contribute at the percent level (Shapiro et al.,
2010; Petiteau, Babak & Sesana, 2011) but the dominant error for these high-z
sources will be LISA intrinsic errors on the measurement of the wave
amplitude.
This redshift uncertainty will translate to a similar uncertainty in the chirp
mass of the binary. Klein et al. (2016) modelled the merger of similar black
hole masses and redshifts and using a SNR = 7 finding that the expected number
of detections would be 358 over a 5 year mission period (their N2A5M5L6 model
in Table II). Their modelling only accounted for the merging of heavy seed
black holes following a galaxy merger and neglected the in-situ models as
described here which could increase the numbers of events potentially
discoverable by LISA significantly. This point was previously investigated by
Hartwig, Agarwal & Regan (2018) as a method to discriminate between seed
formation scenarios. Further, investigation of the number density of in-situ
mergers and how they could influence the number of detections by LISA is now
underway (O’Brennan et al. in prep). The merger of less massive black holes
($\rm{M_{BH}\sim 1000}$ $\rm{M_{\odot}}$) that are produced in these
environments may also be detectable by 3rd generation, ground based, GW
observatories. An exploration of their potential is outside the scope of this
discussion but see Valiante et al. (2020) for more details.

Figure 12: The LISA sensitivity curve with the strain due to two massive black
holes each with a mass of 5000 $\rm{M_{\odot}}$. The merger of two 5000
$\rm{M_{\odot}}~{}$black holes at $z=15$ gives a signal to noise ratio of 32
that is visible by LISA. The inspiral will enter the LISA band at
approximately $5\times 10^{-4}$ Hz, corresponding to a time to merger of
approximately 2 months.

##  5 Discussion and Conclusions

The aim of this study was to examine the prospect of (super-)massive star
formation in haloes that are exposed to moderate LW radiation and also
experience significant dynamical heating effects due to mergers. We identified
79 such haloes in the Renaissance datasets (Regan et al., 2020b) and then
subsequently selected two haloes, which showed no star formation or metal
enrichment prior to reaching the atomic cooling limit, for re-simulation at
enhanced resolution. These halo characteristics make these environments ideal
for forming massive stars.
The two haloes, HaloA and HaloB, were then targeted for re-simulation by
increasing the resolution and by adding additional physics capabilities to
provide a more self-consistent calculation of (super-)massive star formation
in these haloes. Firstly, the spatial resolution was increased by a factor of
256 ($\Delta x\sim 1000$ au) and the mass resolution by a factor of 169
(M${}_{\rm DM}\sim 170$ $\rm{M_{\odot}}$). Additionally, self-shielding of
$\rm{H_{2}}$ from LW was included and the Jeans length was refined by at least
64 cells at all times. To model star formation within the collapsing gas
clouds a star formation prescription was employed, which additionally tracks
accretion onto the stars and allows for the distinction between compact PopIII
stars and massively inflated SMSs.
HaloA was simulated with a LW radiation field impinging on it that followed
exactly what was found in the original simulations and so accounted for star
formation from nearby galaxies. The LW radiation field imposed is shown in
Figure 1. HaloB on the other hand was used as the control halo and no LW field
was imposed to test the impact that turning off the LW field would have.
In the high-resolution re-simulation HaloA formed 22 stars with masses greater
than 1000 $\rm{M_{\odot}}$. The most massive star in HaloA has a mass of 6127
$\rm{M_{\odot}}$. HaloB, re-simulated without a LW radiation field, collapsed
earlier in its evolution before reaching the atomic cooling limit. The lack of
any LW field allows PopIII star formation to occur prior to the onset of
atomic line cooling. In HaloB a total of 21 stars formed by the end of the
simulation with the most massive star in HaloB having a mass of 173
$\rm{M_{\odot}}$. While our results cannot be used to infer a PopIII initial
mass function our results are consistent with the results of other groups who
have found that PopIII mini-haloes tend to host a small number of stars with
peak masses in the range of a few tens of solar masses (Hirano et al., 2014;
Susa, Hasegawa & Tominaga, 2014; Stacy, Bromm & Lee, 2016; Skinner & Wise,
2020).
While the most massive star in HaloA has a mass of more than 6000
$\rm{M_{\odot}}~{}$its accretion rate had fallen to less than $10^{-6}$
$\rm{M_{\odot}~{}yr^{-1}}~{}$by the time the simulation terminated. No
ionising radiative feedback was employed in these simulations and hence
radiative feedback is not the reason for the falloff in accretion. Instead,
while high accretion rates of greater than 0.1
$\rm{M_{\odot}~{}yr^{-1}}~{}$initially fall onto the protostar those rates are
not sustained. The turbulent nature of the collapsing parent cloud causes a
multitude of sub-clouds to form and dissipate. The sub-clouds with masses of
thousands of solar masses form stars but also tidally disrupt other sub-clouds
disrupting accretion onto other (proto-)stars. This phenomenon is observed in
all of the stars formed in the simulation. On average stars accrete above 0.01
$\rm{M_{\odot}~{}yr^{-1}}~{}$for less than 100 kyr before accretion is halted
due to external tidal disruption.
What does this mean for the formation of SMSs and massive black hole seeds?
This is only a single halo simulated for 2 Myr. However, it is likely that the
phenomenon of a turbulent environment is common in the early stages of galaxy
formation and PopIII studies have also converged on the formation of multiple
stars in mini-haloes (e.g. Turk et al., 2012). It is therefore not surprising
that this picture appears in more massive galaxies as well. Chon et al. (2016)
and Chon, Hosokawa & Yoshida (2018) investigated a similar system where, like
here, they selected massive haloes from a larger scale parent simulation.
Using particle splitting techniques they preformed two very high resolution
zoom-in simulations for approximately 100 kyr each. In one case they found
that the central disk fragmented into approximately 15 stars with masses of a
few$\times 10^{3}$ $\rm{M_{\odot}}~{}$\- similar to what we find here and
indeed what was found in the more idealised simulation of Regan & Downes
(2018b) where an external LW field was used to suppress PopIII star formation.
In the other case from Chon, Hosokawa & Yoshida (2018) they found much less
fragmentation leading to an almost spherical collapse and the formation of a
handful of stars with masses in the range of a few $\times 10^{4}$
$\rm{M_{\odot}}$. They concluded that the tidal field surrounding the
collapsing haloes in both cases played a major role. Given that HaloA in this
study experiences a sharp rise in its impinging LW field from a nearby star-
bursting galaxies shortly before its own collapse similar tidal effects are
also expected here. We plan to evaluate the impact of tidal fields on the
formation of haloes hosting massive star formation in an upcoming study (Regan
et al. in prep.).
While the tidal effects can change the dynamics of the large scale disk, below
the parsec scale the tidal effects only act as a seed to further turbulence
which grows through isothermal collapse (Chon, Hosokawa & Yoshida, 2018).
Turbulence is a common feature of present day giant molecular clouds in which
star formation is clearly evident (e.g. Girichidis et al., 2020; Lee et al.,
2020; Krause et al., 2020) and while in both cases the gas is moving super-
sonically the coolants available to the gas in both cases make direct
comparisons difficult. After approximately 2 Myr the stars in HaloA will start
to directly collapse in black holes retaining close to 100% of the stellar
mass (Heger et al., 2003). While this event is unlikely to have any immediate
impact on the mass accretion rate onto the black hole it will leave this halo
with a small population of massive black holes only a few Myr after the onset
of galaxy formation. Allowing the galaxy to evolve further (through a few tens
of dynamical times) may reveal than these intermediate mass black holes can
sink towards the centre of the halo and subsequently merge to form more
massive black holes with masses close to $10^{5}$ $\rm{M_{\odot}}~{}$(similar
to what most canonical galaxy formation simulation use as their seed mass).
Such black hole mergers may even be detectable by LISA as we demonstrate in
Figure 12. However, it is likely that some of this population of IMBHs will
not merge or accrete any significant amount of gas, will display very low duty
cycles and will wander their parent galaxy (Tremmel et al., 2018; Reines et
al., 2020; Barausse et al., 2020).
To develop this model further what will actually be required is further high
resolution simulations which can model the formation of massive stars in high-
infall rate haloes, perhaps independent of metallicity (Chon & Omukai, 2020;
Regan et al., 2020a), which experience a large transfer of baryonic mass
towards their centres. These simulations will equally need to be able to
differentiate between PopIII star formation and SMS star formation, as done
here, but for a longer time period and beyond the formation of the subsequent
black holes. Tracking the evolution of these truly seed black holes is the
next frontier.

## Acknowledgments

JR acknowledges support from the Royal Society and Science Foundation Ireland
under grant number URF$\backslash$R1$\backslash$191132\. JHW is supported by
National Science Foundation grants AST-1614333 and OAC-1835213, and NASA
grants NNX17AG23G and 80NSSC20K0520. BWO acknowledges support from NSF grants
PHY-1430152, AST-1514700, AST-1517908, and OAC-1835213, by NASA grants
NNX12AC98G and NNX15AP39G, and by HST-AR-13261 and HST-AR-14315. JR wishes to
acknowledge the DJEI/DES/SFI/HEA Irish Centre for High-End Computing (ICHEC)
for the provision of computational facilities and support on which the zoom
simulations were run. The authors also thank the Extreme Gravity Institute at
Montana State University for making the LISA sensitivity curve scripts
available on github. The original Renaissance simulations were performed on
Blue Waters, which is operated by the National Center for Supercomputing
Applications (NCSA) with PRAC allocation support by the NSF (awards
ACI-1238993, ACI-1514580, and OAC-1810584). This research is part of the Blue
Waters sustained-petascale computing project, which is supported by the NSF
(awards OCI-0725070, ACI-1238993) and the state of Illinois. Blue Waters is a
joint effort of the University of Illinois at Urbana-Champaign and its NCSA.
The freely available plotting library matplotlib (Hunter, 2007) was used to
construct numerous plots within this paper. Computations and analysis
described in this work were performed using the publicly-available Enzo (Bryan
et al., 2014; Brummel-Smith et al., 2019) and yt  (Turk et al., 2011) codes,
which are the product of a collaborative effort of many independent scientists
from numerous institutions around the world. Their commitment to open science
has helped make this work possible. Finally, the authors thank the anonymous
referee for a detailed and very helpful report.

## References

  * Abel et al. (1997) Abel T., Anninos P., Zhang Y., Norman M. L., 1997, New Astronomy, 2, 181
  * Abel, Bryan & Norman (2000) Abel T., Bryan G. L., Norman M. L., 2000, ApJ, 540, 39
  * Abel, Bryan & Norman (2002) Abel T., Bryan G. L., Norman M. L., 2002, Science, 295, 93
  * Alvarez, Wise & Abel (2009) Alvarez M. A., Wise J. H., Abel T., 2009, ApJ, 701, L133
  * Anninos et al. (1997) Anninos P., Zhang Y., Abel T., Norman M. L., 1997, New Astronomy, 2, 209
  * Badnell (2006) Badnell N. R., 2006, ApJS, 167, 334
  * Badnell et al. (2003) Badnell N. R. et al., 2003, A&A, 406, 1151
  * Barausse et al. (2020) Barausse E., Dvorkin I., Tremmel M., Volonteri M., Bonetti M., 2020, ApJ, 904, 16
  * Becerra et al. (2015) Becerra F., Greif T. H., Springel V., Hernquist L. E., 2015, MNRAS, 446, 2380
  * Begelman & Rees (1978) Begelman M. C., Rees M. J., 1978, MNRAS, 185, 847
  * Begelman, Rossi & Armitage (2008) Begelman M. C., Rossi E. M., Armitage P. J., 2008, MNRAS, 387, 1649
  * Begelman, Volonteri & Rees (2006) Begelman M. C., Volonteri M., Rees M. J., 2006, MNRAS, 370, 289
  * Bleuler & Teyssier (2014) Bleuler A., Teyssier R., 2014, MNRAS, 445, 4015
  * Brummel-Smith et al. (2019) Brummel-Smith C. et al., 2019, The Journal of Open Source Software, 4, 1636
  * Bryan et al. (2014) Bryan G. L., Norman M. L., O’Shea B. W., Abel T., Wise J. H., Turk M. J., The Enzo Collaboration, 2014, ApJS, 211, 19
  * Chen et al. (2014) Chen P., Wise J. H., Norman M. L., Xu H., O’Shea B. W., 2014, ApJ, 795, 144
  * Chon et al. (2016) Chon S., Hirano S., Hosokawa T., Yoshida N., 2016, ApJ, 832, 134
  * Chon, Hosokawa & Yoshida (2018) Chon S., Hosokawa T., Yoshida N., 2018, MNRAS, 475, 4104
  * Chon & Omukai (2020) Chon S., Omukai K., 2020, MNRAS
  * Clark, Glover & Klessen (2008) Clark P. C., Glover S. C. O., Klessen R. S., 2008, ApJ, 672, 757
  * Cornish & Shuman (2020) Cornish N. J., Shuman K., 2020, Phys. Rev. D, 101, 124008
  * Crosby et al. (2013) Crosby B. D., O’Shea B. W., Smith B. D., Turk M. J., Hahn O., 2013, ApJ, 773, 108
  * Davies, Miller & Bellovary (2011) Davies M. B., Miller M. C., Bellovary J. M., 2011, ApJ, 740, L42
  * Dere et al. (1997) Dere K. P., Landi E., Mason H. E., Monsignori Fossi B. C., Young P. R., 1997, A&AS, 125, 149
  * Downes (2012) Downes T. P., 2012, MNRAS, 425, 2277
  * Dunn, Holley-Bockelmann & Bellovary (2020) Dunn G., Holley-Bockelmann K., Bellovary J., 2020, ApJ, 896, 72
  * eLISA Consortium et al. (2013) eLISA Consortium et al., 2013, ArXiv e-prints 1305.5720
  * Fan et al. (2006) Fan X. et al., 2006, AJ, 131, 1203
  * Federrath et al. (2010) Federrath C., Roman-Duval J., Klessen R. S., Schmidt W., Mac Low M. M., 2010, A&A, 512, A81
  * Ferland et al. (2017) Ferland G. J. et al., 2017, Rev. Mexicana Astron. Astrofis., 53, 385
  * Fernandez et al. (2014) Fernandez R., Bryan G. L., Haiman Z., Li M., 2014, MNRAS, 439, 3798
  * Fitzgerald et al. (2020) Fitzgerald K., Harvey E. J., Keaveney N., Redman M. P., 2020, Astronomy and Computing, 32, 100382
  * Freitag (2008) Freitag M., 2008, in Astronomical Society of the Pacific Conference Series, Vol. 387, Massive Star Formation: Observations Confront Theory, Beuther H., Linz H., Henning T., eds., pp. 247–+
  * Freitag, Gürkan & Rasio (2006) Freitag M., Gürkan M. A., Rasio F. A., 2006, MNRAS, 368, 141
  * Georgiev et al. (2016) Georgiev I. Y., Böker T., Leigh N., Lützgendorf N., Neumayer N., 2016, MNRAS, 457, 2122
  * Girichidis et al. (2020) Girichidis P. et al., 2020, Space Sci. Rev., 216, 68
  * Greif et al. (2011) Greif T. H., Springel V., White S. D. M., Glover S. C. O., Clark P. C., Smith R. J., Klessen R. S., Bromm V., 2011, ApJ, 737, 75
  * Gültekin, Miller & Hamilton (2006) Gültekin K., Miller M. C., Hamilton D. P., 2006, ApJ, 640, 156
  * Gürkan, Freitag & Rasio (2004) Gürkan M. A., Freitag M., Rasio F. A., 2004, ApJ, 604, 632
  * Haemmerlé et al. (2018) Haemmerlé L., Woods T. E., Klessen R. S., Heger A., Whalen D. J., 2018, MNRAS, 474, 2757
  * Hartwig, Agarwal & Regan (2018) Hartwig T., Agarwal B., Regan J. A., 2018, MNRAS, 479, L23
  * Heger et al. (2003) Heger A., Fryer C. L., Woosley S. E., Langer N., Hartmann D. H., 2003, ApJ, 591, 288
  * Hirano et al. (2014) Hirano S., Hosokawa T., Yoshida N., Umeda H., Omukai K., Chiaki G., Yorke H. W., 2014, ApJ, 781, 60
  * Holley-Bockelmann et al. (2008) Holley-Bockelmann K., Gültekin K., Shoemaker D., Yunes N., 2008, ApJ, 686, 829
  * Hosokawa, Omukai & Yorke (2013) Hosokawa T., Omukai K., Yorke H. W., 2013, ApJ, 778, 178
  * Hosokawa et al. (2013) Hosokawa T., Yorke H. W., Inayoshi K., Omukai K., Yoshida N., 2013, ApJ, 778, 178
  * Hunter (2007) Hunter J. D., 2007, Computing in Science & Engineering, 9, 90
  * Inayoshi, Haiman & Ostriker (2016a) Inayoshi K., Haiman Z., Ostriker J. P., 2016a, MNRAS, 459, 3738
  * Inayoshi, Haiman & Ostriker (2016b) Inayoshi K., Haiman Z., Ostriker J. P., 2016b, MNRAS, 459, 3738
  * Inayoshi, Omukai & Tasker (2014) Inayoshi K., Omukai K., Tasker E., 2014, MNRAS, 445, L109
  * Johnson et al. (2012) Johnson J. L., Whalen D. J., Fryer C. L., Li H., 2012, ApJ, 750, 66
  * Katz, Sijacki & Haehnelt (2015) Katz H., Sijacki D., Haehnelt M. G., 2015, MNRAS, 451, 2352
  * Kitsionas & Whitworth (2002) Kitsionas S., Whitworth A. P., 2002, MNRAS, 330, 129
  * Klein et al. (2016) Klein A. et al., 2016, Phys. Rev. D, 93, 024003
  * Komatsu et al. (2011) Komatsu E. et al., 2011, ApJS, 192, 18
  * Kormendy & Ho (2013) Kormendy J., Ho L. C., 2013, ARA&A, 51, 511
  * Krause et al. (2020) Krause M. G. H. et al., 2020, Space Sci. Rev., 216, 64
  * Krumholz, McKee & Klein (2004) Krumholz M. R., McKee C. F., Klein R. I., 2004, ApJ, 611, 399
  * Landi et al. (2012) Landi E., Del Zanna G., Young P. R., Dere K. P., Mason H. E., 2012, ApJ, 744, 99
  * Latif, Schleicher & Hartwig (2016) Latif M. A., Schleicher D. R. G., Hartwig T., 2016, MNRAS, 458, 233
  * Lee et al. (2020) Lee Y.-N., Offner S. S. R., Hennebelle P., André P., Zinnecker H., Ballesteros-Paredes J., Inutsuka S.-i., Kruijssen J. M. D., 2020, Space Sci. Rev., 216, 70
  * Lupi et al. (2014) Lupi A., Colpi M., Devecchi B., Galanti G., Volonteri M., 2014, MNRAS, 442, 3616
  * Machacek, Bryan & Abel (2001) Machacek M. E., Bryan G. L., Abel T., 2001, ApJ, 548, 509
  * Meyer et al. (2020) Meyer D. M. A., Vorobyov E. I., Elbakyan V. G., Eislöffel J., Sobolev A. M., Stöhr M., 2020, MNRAS
  * Milosavljević, Couch & Bromm (2009) Milosavljević M., Couch S. M., Bromm V., 2009, ApJ, 696, L146
  * Mirocha et al. (2012) Mirocha J., Skory S., Burns J. O., Wise J. H., 2012, ApJ, 756, 94
  * Morawski et al. (2018) Morawski J., Giersz M., Askar A., Belczynski K., 2018, MNRAS, 481, 2168
  * Morisset (2013) Morisset C., 2013, pyCloudy: Tools to manage astronomical Cloudy photoionization code
  * Natarajan (2020) Natarajan P., 2020, MNRAS
  * Omukai & Palla (2003) Omukai K., Palla F., 2003, ApJ, 589, 677
  * O’Shea et al. (2005) O’Shea B. W., Abel T., Whalen D., Norman M. L., 2005, ApJ, 628, L5
  * O’Shea & Norman (2008) O’Shea B. W., Norman M. L., 2008, ApJ, 673, 14
  * O’Shea et al. (2015) O’Shea B. W., Wise J. H., Xu H., Norman M. L., 2015, ApJ, 807, L12
  * Osterbrock & Ferland (2006) Osterbrock D. E., Ferland G. J., 2006, Astrophysics of gaseous nebulae and active galactic nuclei
  * Pacucci et al. (2017) Pacucci F., Natarajan P., Volonteri M., Cappelluti N., Urry C. M., 2017, ApJ, 850, L42
  * Pacucci, Volonteri & Ferrara (2015) Pacucci F., Volonteri M., Ferrara A., 2015, MNRAS, 452, 1922
  * Park & Bogdanović (2019) Park K., Bogdanović T., 2019, ApJ, 883, 209
  * Park, Wise & Bogdanović (2017) Park K., Wise J. H., Bogdanović T., 2017, ApJ, 847, 70
  * Petiteau, Babak & Sesana (2011) Petiteau A., Babak S., Sesana A., 2011, ApJ, 732, 82
  * Pfister et al. (2019) Pfister H., Volonteri M., Dubois Y., Dotti M., Colpi M., 2019, MNRAS, 486, 101
  * Planck Collaboration et al. (2014) Planck Collaboration et al., 2014, A&A, 571, A16
  * Portegies Zwart et al. (2004) Portegies Zwart S. F., Baumgardt H., Hut P., Makino J., McMillan S. L. W., 2004, Nature, 428, 724
  * Porter et al. (2012) Porter R. L., Ferland G. J., Storey P. J., Detisch M. J., 2012, MNRAS, 425, L28
  * Rees (1978) Rees M. J., 1978, Phys. Scr, 17, 193
  * Regan & Downes (2018a) Regan J. A., Downes T. P., 2018a, MNRAS, 475, 4636
  * Regan & Downes (2018b) Regan J. A., Downes T. P., 2018b, MNRAS, 478, 5037
  * Regan & Haehnelt (2009) Regan J. A., Haehnelt M. G., 2009, MNRAS, 396, 343
  * Regan et al. (2020a) Regan J. A., Haiman Z., Wise J. H., O’Shea B. W., Norman M. L., 2020a, The Open Journal of Astrophysics, 3, E9
  * Regan, Johansson & Wise (2015) Regan J. A., Johansson P. H., Wise J. H., 2015, MNRAS, 449, 3766
  * Regan et al. (2017) Regan J. A., Visbal E., Wise J. H., , Haiman Z., Johansson P. H., Bryan G. L., 2017, Nature Astronomy, 1, 0075
  * Regan et al. (2020b) Regan J. A., Wise J. H., O’Shea B. W., Norman M. L., 2020b, MNRAS, 492, 3021
  * Reines et al. (2020) Reines A. E., Condon J. J., Darling J., Greene J. E., 2020, ApJ, 888, 36
  * Robson, Cornish & Liu (2019) Robson T., Cornish N. J., Liu C., 2019, Classical and Quantum Gravity, 36, 105011
  * Sakurai, Haiman & Inayoshi (2020) Sakurai Y., Haiman Z., Inayoshi K., 2020, MNRAS, 499, 5960
  * Sakurai, Inayoshi & Haiman (2016) Sakurai Y., Inayoshi K., Haiman Z., 2016, MNRAS, 461, 4496
  * Sakurai et al. (2016) Sakurai Y., Vorobyov E. I., Hosokawa T., Yoshida N., Omukai K., Yorke H. W., 2016, MNRAS, 459, 1137
  * Schaerer (2002) Schaerer D., 2002, A&A, 382, 28
  * Sesana (2013) Sesana A., 2013, in Astronomical Society of the Pacific Conference Series, Vol. 467, 9th LISA Symposium, Auger G., Binétruy P., Plagnol E., eds., p. 103
  * Sesana (2016) Sesana A., 2016, Phys. Rev. Lett., 116, 231102
  * Shapiro et al. (2010) Shapiro C., Bacon D. J., Hendry M., Hoyle B., 2010, MNRAS, 404, 858
  * Skinner & Wise (2020) Skinner D., Wise J. H., 2020, MNRAS, 492, 4386
  * Smith et al. (2017) Smith B. D. et al., 2017, MNRAS, 466, 2217
  * Smith et al. (2018) Smith B. D., Regan J. A., Downes T. P., Norman M. L., O’Shea B. W., Wise J. H., 2018, MNRAS, 480, 3762
  * Stacy, Bromm & Lee (2016) Stacy A., Bromm V., Lee A. T., 2016, MNRAS, 462, 1307
  * Stahler, Palla & Salpeter (1986) Stahler S. W., Palla F., Salpeter E. E., 1986, ApJ, 302, 590
  * Storey & Hummer (1995) Storey P. J., Hummer D. G., 1995, MNRAS, 272, 41
  * Susa, Hasegawa & Tominaga (2014) Susa H., Hasegawa K., Tominaga N., 2014, ApJ, 792, 32
  * Toyouchi et al. (2020) Toyouchi D., Hosokawa T., Sugimura K., Kuiper R., 2020, MNRAS
  * Tremmel et al. (2018) Tremmel M., Governato F., Volonteri M., Pontzen A., Quinn T. R., 2018, ApJ, 857, L22
  * Turk, Abel & O’Shea (2009) Turk M. J., Abel T., O’Shea B., 2009, Science, 325, 601
  * Turk et al. (2012) Turk M. J., Oishi J. S., Abel T., Bryan G. L., 2012, ApJ, 745, 154
  * Turk et al. (2011) Turk M. J., Smith B. D., Oishi J. S., Skory S., Skillman S. W., Abel T., Norman M. L., 2011, ApJS, 192, 9
  * Valiante et al. (2020) Valiante R. et al., 2020, MNRAS
  * Whalen, Abel & Norman (2004) Whalen D., Abel T., Norman M. L., 2004, ApJ, 610, 14
  * Wise & Abel (2011) Wise J. H., Abel T., 2011, MNRAS, 414, 3458
  * Wise et al. (2014) Wise J. H., Demchenko V. G., Halicek M. T., Norman M. L., Turk M. J., Abel T., Smith B. D., 2014, MNRAS, 442, 2560
  * 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
  * Wise et al. (2012) Wise J. H., Turk M. J., Norman M. L., Abel T., 2012, ApJ, 745, 50
  * Wolcott-Green, Haiman & Bryan (2011) Wolcott-Green J., Haiman Z., Bryan G. L., 2011, MNRAS, 418, 838
  * Woods et al. (2019) Woods T. E. et al., 2019, PASA, 36, e027
  * Woods, Heger & Haemmerlé (2020) Woods T. E., Heger A., Haemmerlé L., 2020, MNRAS, 494, 2236
  * Woods et al. (2017) Woods T. E., Heger A., Whalen D. J., Haemmerlé L., Klessen R. S., 2017, ApJ, 842, L6
  * Xu et al. (2014) Xu H., Ahn K., Wise J. H., Norman M. L., O’Shea B. W., 2014, ApJ, 791, 110
  * Xu, Wise & Norman (2013) Xu H., Wise J. H., Norman M. L., 2013, ApJ, 773, 83