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

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 M$_{\odot}$. The most massive star has a stellar mass of over 6000 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 M$_{\odot}$ and 10,000 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.


INTRODUCTION
Supermassive stars (SMSs) with masses between 10 4 and 10 5 M 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) 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 there-fore 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 M 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;Woods et al. 2017). However, the SMS must continue to accrete above this threshold rate. If the accretion rate falls below the critical rate 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 M yr −1 , however, a sufficient baryon reservoir is also required and furthermore the metallicity of the gas being accreted should be below 10 −3 Z (Chon & Omukai 2020). For these reasons metal-poor (i.e., Z 10 −3 Z ) 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 Z ) 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 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 M 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 J LW ∼ 2 J 21 existed). In HaloB the most massive star in the halo has a mass of approximately 173 M . The impact of ionising radiation is not considered in these simulations but postprocessing 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.

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. 2012Wise et al. , 2014Regan, 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 Abel et al. 1997). In particular the photo-dissociation of 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 Brummel-Smith et al. 2019) 1 . 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 × 10 4 M . 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 halo ∼ 2 × 10 6 M .
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, ∆x ∼ 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 some showed almost ideal isothermal collapse 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.
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 found (z = 15.6 for HaloA and z = 15.0 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 DM ∼ 170 M inside the target halo and its surrounding Lagrangian volume.

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. 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 simulation and allows us to examine the impact of what happens when no LW field is present.

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 ∆x ∼ 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, Mc-Kee & 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. The cell is at the highest refinement level 2. The cell exceeds the Jeans density 3. The flow around the cell is converging 4. The cooling time of the cell is less than the freefall time 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) whereṀ is the mass accretion rate, Ω is the accretion zone over which we integrate, ρ 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 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 M min,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).
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. 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).
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 M 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 radia-tion 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 M (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 T eff = 5500 K ). 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 M 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 H 2 for photons with energies within the Lyman-Werner band and the photo-detachment of H − and 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 H − , H + 2 and 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 H 2 against LW radiation we use the prescription of Wolcott-Green, Haiman & Bryan (2011). 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 ∼ 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.
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 formation.

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 singleelectron 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 × 10 −3 M yr −1 , which will evolve on the Hayashi track due to H − opacity in their inflated envelopes (e.g., , 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 ≈ 10 5 K (Schaerer 2002;Woods, Heger & Haemmerlé 2020). We further take their luminosities to be approximately Eddington (L ≈ 1.3 × 10 38 × (M/M ) 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 (∼ 3 × 10 −4 pc) and terminate our calculations at the outer boundary of each nebula once either the gas tem- perature 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 §3.2.

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 resimulation. 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 J LW ∼ 2 J 21 until a redshift of z ∼ 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 radiation 2 ) in the direction of HaloA. The LW field reaches its zenith at a redshift of z ∼ 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 ∼ 16.2 and lasts until z ∼ 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 HaloA = 9.3 × 10 7 M ) and the halo has remained metal-free. HaloB, on the other hand, collapses early with a mass of M HaloB = 3.7×10 6 M due to the lack of a LW background and thus can be classi-fied 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.

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.
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 sphericallysymmetric 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Ṅ 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: where n H (r, t) is the spherically-averaged density profile centered on the most massive star in the Halo at time t (see above), and α B (H 0 , T) is the case B recombination coefficient for hydrogen, which we take as ≈ 1.4 × 10 −13 cm 3 s −1 for primordial gas at a temperature typical of plasmas photoionised by Population III stars, T ∼ 2 × 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 vol-ume (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 ∼ 500 and ∼ 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 M and crucially it moves to a lower density region away from the central overdensity. 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). 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 density 3 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 M and is 28 kyr old, the second star has a mass of 6 M and is 4 kyr old. Stars are coloured in red if the accretion rate exceeds the SMS critical rate of 0.04 M 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 su- However, it should be noted that due to our finite resolution (∆x ∼ 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.

The case for super-massive star formation
In Figure 4 we saw that the mass of the most massive star in HaloA was 6127 M . 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 ∼ 10 5 M . 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, 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 H 2 fraction is very different between the two haloes, with HaloA showing a steep decline in H 2 towards the centre due to the combination of the LW background and the rapid assembly of mass. HaloB on the otherhand shows a more typical 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 M 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 M yr −1 at 1 pc down to approximately 10 −3 M yr −1 in the very centre. As we will see this result is reflected in the 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 the simulation itself. In Figure 7 we show the mass evolution and the accretion rate onto the most massive star in each halo. In the left panel we show the mass evolution as a function of the stellar age. The most massive star in HaloA is denoted by the blue line, the most massive star in HaloB by the red line. The dashed line in the panel shows the mass evolution of a star with the median 4 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 M within approximately 20 kyr after which it stops growing. The mass evolution of the most massive star in HaloB is a little different. In this case the star accretes slower at approximately 10 −3 M yr −1 (see right panel). What then drives the star to increase its mass is a stellar merger with another star. This occurs when the main progenitor star is approximately 200 kyr old and results in a star of 173 M .
In the right hand panel of Figure 7 we see the accretion rates onto the most massive star and onto the median star for each simulation. Initially each star follows 10 3 10 2 10 1 10 0 10 1 10 2  Figure 6. We examine the radial profiles for both HaloA (blue line) and HaloB (red line). The four panels are from bottom left going clockwise: temperature, 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 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 H 2 fractions. The lower 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 0.1 pc) of HaloA are above 0.1 M yr −1 compared to below 0.001 M 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).

Radius [pc]
approximately the mass accretion rate as found in the radial profiles (see Figure 6). However, in both cases 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 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.
By the end of each simulation the total stellar mass in HaloA is 90,000 M and 1300 M in HaloB. The star formation efficiency in the inner 20 pc (5 pc) of HaloA (HaloB) is 28% (18%). 20 pc is the approximate Jeans length of the gas in HaloA while 5 pc is the approximate Jeans length of the gas in HaloB. In order for the formation of a truly supermassive star in HaloA all of the stellar mass (i.e. ∼90,000 M ) 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 subclouds 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 300,000 M inside 20 pc, approximately 2% of the baryons went into the most massive star. To form a SMS with a mass of close to 100,000 M we would need a third 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. What appears more likely, for this halo at least, is that the most massive star(s) in the halo will form of population of intermediate mass black holes (IMBHs) with masses ranging from 300 M to several thousand M . The subsequent growth of these IMBHs within the halo then through mergers and accretion will then slowly build the black hole mass. Nonetheless, the key issue with forming 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.

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 8 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, t coll , for a stellar population consisting of two populations can be written as (Freitag 2008): (3) where n * is the stellar density, σ vel is the one dimensional velocity dispersion, R 1 is the stellar radius of the first population of stars and R 2 is the stellar radius of the second population. We choose the median stellar mass, 683 M , for M 1 and we choose our maximum stellar mass, 6127 M , for M 2 . 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 re-quired for the system to contract due to two body interactions in a populous cluster (N 10). The relaxation time is given by: where λ = ln(0.02N ) and N is the number of cluster stars and m is the median stellar mass. The collisional time, t coll , is greater than the Hubble time at all radii indicating that collisions are very unlikely for the stars in HaloA. The relaxation time, t relax , 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). This result is backed up by more recent investigations of black hole merging in similar environments as we now discuss. At the end of the simulation HaloA has 22 stars with masses greater than 1000 M with the most massive star having a mass of 6127 M . The simulations show that accretion onto the stars is focused on the early stages in the life of the star with most of the mass accreted within the first 100 kyr. Therefore, if we assume that the proto-galaxy gets populated with a plethora of weakly accreting black holes early in the proto-galactic evolution then it is 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 (∆x ∼ 10 pc) with black hole masses in the range 10 4 M to 10 5 M in highz 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 highz 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 M 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. Toyouchi et al. (2020) investigated a more idealised scenario where they modelled the dynamics of a single 10 4 M mass black hole incorporating the effects of dynamical friction and radiative feedback. They found that if the gas density is low that the black hole does not sink to the centre confirming the results of Pfister et al. (2019). However, in the case where the black hole encounters dense gas 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 much shorter than the dynamical time. 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 ∼ 10 4 M ) 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 which may help to supply the necessary ram pressure and extract orbital energy. 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 M (upper limit of pair instability supernovae) up to 6000 M (and likely higher). In Figure 9 we plot the characteristic strain expected from the merger of two black holes with masses of 5000 M each at z = 15. 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 M at z = 15 is approximately 32 and hence well within LISA's detection parameters 5 . The merger of a 10,000 M binary black hole system will enter the LISA band several 1 month before merger and complete thousands of orbits before the final plunge. With this number of cycles LISA will be able to detect the redshifted mass with a precision of close to 1% (Sesana 2013). 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 dom- inant 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 (Arridge et al. in prep).

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 (∆x ∼ 1000 au) and the mass resolution by a factor of 169 (M DM ∼ 170 M ). Additionally, self-shielding of 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 M . The most massive star in HaloA has a mass of 6127 M . HaloB, resimulated 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 M . 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 M its accretion rate had fallen to less than 10 −6 M 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 M 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 subclouds 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 M 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 minihaloes (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×10 3 M -similar to what we find here. In the other case they found much less fragmentation leading to an almost spherical collapse and the formation of a handful of stars with masses in the range  . The LISA sensitivity curve with the strain due to two massive black holes each with a mass of 5000 M . The merger of two 5000 M 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 × 10 −4 Hz, corresponding to a time to merger of approximately 1 month. of a few ×10 4 M . 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 starbursting 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 M (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 9. 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\R1\191132.
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 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.