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

Halos with masses in excess of the atomic limit are believed to be ideal environments in which to form heavy black hole seeds with masses above 10^3 Msun. In cases where the H_2 fraction is suppressed this is expected to lead to reduced fragmentation of the gas and the generation of a top heavy initial mass function. In extreme cases this can result in the formation of massive black hole seeds. Resolving the initial fragmentation scale and the resulting protostellar masses has, until now, not been robustly tested. Cosmological simulations were performed with the moving mesh code Arepo using a primordial chemistry network until z = 11. Three haloes with masses in excess of the atomic cooling mass were then selected for detailed examination via zoom-ins. The highest resolution simulations resolve densities up to 10^-6 g cm^-3 (10^18 cm^-3) and capture a further 100 yr of fragmentation behaviour at the center of the halo. Our simulations show intense fragmentation in the central region of the halos, leading to a large number of near-solar mass protostars. Despite the increased fragmentation the halos produce a protostellar mass spectrum that peaks at higher masses relative to standard Population III star forming halos. The most massive protostars have accretion rates of 10^-3-10^-1 Msun yr^-1 after the first 100 years of evolution, while the total mass of the central region grows at 1 Msun yr^-1. Lower resolution zoom-ins show that the total mass of the system continues to accrete at 1 Msun yr^-1 for at least 10^4 yr, although how this mass is distributed amongst the rapidly growing number of protostars is unclear. However, assuming that a fraction of stars can continue to accrete rapidly the formation of a sub-population of stars with masses in excess of 10^3 Msun is likely in these halos.


Introduction
Quasi-stellar radio sources (quasars) are thought to be powered by supermassive black holes (SMBHs) which populate the centers of most, if not all, massive galaxies.Quasars have been detected out to redshifts of z > 7 with black hole (BH) masses believed to be in excess of 10 9 M ⊙ (e.g.Mortlock et al. 2011;Matsuoka et al. 2019), implying that the seeds for these SMBHs formed in the early Universe.The existence of SMBHs with masses in excess of 10 9 M ⊙ within the first billion years of the Universe poses a challenge to our understanding of both BH formation and BH accretion.How could such massive objects appear so early in cosmic history?Two mainstream pathways have emerged over the last four decades; SMBHs may originate from so-called light seeds with masses less than 10 3 M ⊙ , or from heavy BH seeds with masses significantly in excess of 10 3 M ⊙ .
which only need to last for brief periods of time (Lupi et al. 2016;Inayoshi et al. 2016).However, in both Eddington limited and super-Eddington cases, radiative feedback from the accretion of matter onto BHs heats the surrounding gas and lowers the accretion rate (Johnson & Bromm 2007;Milosavljević et al. 2009;Alvarez et al. 2009;Smith et al. 2018) -making maintaining either Eddington accretion and/or super-Eddington accretion unlikely over sustained periods (Regan et al. 2019;Su et al. 2023;Massonneau et al. 2023).Finally, both Eddington and super-Eddington limited growth requires that the light seed sit at the centre of a powerful gas inflow and that the embryonic BH can readily accrete the surrounding dense gas.However, dynamical studies of BHs have shown that this is also a challenge, with light seeds tending to walk random trajectories around the host halo centres (Beckmann et al. 2019;Pfister et al. 2019).As a result of these obstacles to light seed growth, the possibility of much heavier BH seeds has also been studied, beginning with Rees (1978).
Forming heavy seeds has hinged on two separate but not necessarily distinct pathways.On the one hand dynamical processes have been invoked to explain heavy BH seed formation; runaway collisions in dense young star clusters could produce massive black holes (MBHs) of ∼ 10 3 M ⊙ (e.g.Portegies Zwart et al. 2004;Glebbeek et al. 2009;Katz et al. 2015;Reinoso et al. 2023) that are candidates for the ultraluminous X-ray sources observed in young star forming regions (Ptak & Colbert 2004).These MBHs may grow into SMBHs through binary mergers and/or gas accretion (Micic et al. 2007).Additionally, BH mergers within a dense BH cluster may achieve the same outcome (Stone et al. 2017;Schleicher et al. 2022) through either collisions and/or gas accretion, although the growth prospects are far from certain (Arca Sedda et al. 2023).
On the other hand conditions may exist in the early Universe conducive to the formation of truly massive stars with final masses well in excess of 10 3 M ⊙ (Regan et al. 2020;Latif et al. 2022;Regan et al. 2023) and possibly up as high as a few times 10 5 M ⊙ (Woods et al. 2017).The formation of these massive primordial Pop III stars requires high inflow rates onto the stellar surface (Haemmerlé et al. 2018;Woods et al. 2020) but repeated numerical experiments have shown the stars to be stable for at least 2 Myr until their inevitable direct collapse into a MBH (e.g.Hosokawa & Omukai 2009;Hosokawa et al. 2012a).Finally, it may also be that the processes that lead to a dense stellar cluster or a massive Pop III star form part of a continuum with a massive star forming at the very centre of a dense stellar cluster where collisions drive the formation of a very massive star (Boekholt et al. 2018;Chon & Omukai 2020;Schleicher et al. 2023).The goal of this paper will be to test whether so-called atomic cooling halos can lead to the formation of a heavy seed Pop III star using state-of-the-art, high resolution, hydrodynamic simulations.
Pristine atomic cooling halos have long been suggested as the ideal environment in which to seed MBHs (Loeb & Rasio 1994;Spaans & Silk 2006;Prieto et al. 2013).If the H 2 abundance inside the massive halo can be suppressed then the gas must cool predominatly through atomic hydrogen line emission and H − free-bound emission.The thermal pathway then taken by the gas inside an atomic cooling halo in the absence of effective H 2 cooling deviates significantly from the standard Pop III star formation scenario (e.g.Omukai 2000;Klessen & Glover 2023).
Pop III star formation in minihalos is facilitated by molecular hydrogen.The DM halo potential well pulls in the gas and shock heats it up to ∼ 1000 K, At these temperatures, the H 2 abundance increases to ∼ 10 −4 (e.g.Tegmark et al. 1997;Greif et al. 2008) where rotational transitions can occur via electrical quadrupole radiation, which allows the gas to cool down to minimum temperatures of ∼200 K (e.g.Abel et al. 1997;Bromm et al. 1999;Glover & Abel 2008) and collapse, decoupling from the DM halo.
One potential way in which H 2 abundances can be reduced is via a nearby source of Lyman-Werner (LW) radiation.These farultraviolet (FUV) photons in the Lyman and Werner bands of H 2 (11.2 -13.6 eV) can dissociate H 2 via the two-step Solomon process (Field et al. 1966;Stecher & Williams 1967).Additionally, photons of above 0.76eV can photodissociate H − , disrupting the primary H 2 formation channel (e.g.Chuzhoy et al. 2007).Before the Strömgren spheres of Pop III stars overlap, the UV background below the ionization threshold is able to penetrate large clouds and suppress their H 2 abundance (Haiman et al. 1997).This photodissociation of H 2 suppresses further star formation inside small halos and delays reionization until larger halos form (Haiman et al. 2000).While other physical mechanism can have a similar effect, we are focused here in studying the gas collapse inside atomically cooling halos, which are both pristine and have had their H 2 cooling efficiency suppressed.
In halos with virial temperatures below T vir ∼ 8000 K, star formation is suppressed entirely if LW radiation reduces the H 2 abundance below the level at which gas can cool within a Hubble time.For example, an intense burst of LW radiation from a neighbouring star-bursting protogalaxy just before the gas cloud undergoes gravitational collapse is proposed to prevent the cloud from collapsing or forming stars (Regan et al. 2017).However, the halo will continue to grow through hierarchical mergers (e.g.Chon et al. 2016;Dong et al. 2022).Once the virial temperature reaches ∼ 8000 K, it becomes possible for the gas to cool via Lyman-α emission.The required virial temperature is related to the virial mass through the relation (Fernandez et al. 2014) giving a virial mass of ∼ 3 × 10 7 M ⊙ at z = 12 for a virial temperature of 8000 K.If a halo can grow to this mass it will become hot enough to cool via atomic line emission and begin to collapse (Oh & Haiman 2002;Bromm & Loeb 2003;Bromm & Yoshida 2011).In this scenario, the collapse occurs almost isothermally, and fragmentation of the gas is thought to be suppressed throughout.
The maximum density that can be reached in simulations of this process is related to the resolution of the simulation.To avoid artificial fragmentation, it is necessary to resolve the Jeans length, which progressively shrinks as the gas collapses to higher densities (Truelove et al. 1997).Therefore, the better the resolution of the simulation, the higher the density that it can reach.We know from simulations of the standard Pop III star formation scenario that the formation of the primordial protostar occurs at densities of 10 −6 -10 −4 g cm −3 where the gas becomes adiabatic (Omukai 2000;Machida & Nakamura 2015).The Jeans length at this stage is 0.01-0.1 au, and so the required resolution is roughly a factor of ten smaller than this.Despite this, the resolution of most atomically cooled halo simulations is relatively poor in comparison, with most studies not resolving their gas past densities of 10 −17 g cm −3 (e.g.Shang et al. 2010;Sugimura et al. 2014;Regan et al. 2014;Hartwig et al. 2015;Agarwal & Khochfar 2015;Glover 2015;Agarwal et al. 2016;Regan et al. 2017;Dunn et al. 2018).This resolution is sufficient to determine whether or not H 2 cooling is important during the initial collapse of the gas, but does not allow one to draw conclusions about the later stages of the collapse.Some studies find evidence for smallscale fragmentation, even in the absence of effective H 2 cooling (e.g.Becerra et al. 2015Becerra et al. , 2018a;;Chon et al. 2018;Latif et al. 2020;Patrick et al. 2023).The number of fragments formed is generally much smaller than the number found in recent simulations of the standard Pop III star formation scenario.However, this may be a consequence of the limited peak density: in most of these studies, the gas density never exceeds ρ ∼ 10 −10 g cm −3 , four orders of magnitude smaller than the point at which we expect the collapse to become adiabatic (Becerra et al. 2018b).
This study aims to provide the most accurate picture of atomically cooled halo collapse at high densities to date, answering whether atomic cooling halos do experience reduced fragmentation and higher stellar/BH seed masses compared to the Pop III minihalo scenario, or whether the fragmentation at high densities produces protostellar masses similar to what we have seen in simulations of Pop III star forming minihalos.To that end, we simulate the collapse of atomically cooled halos from cosmological initial conditions with zoom-in simulations running up to the lower limit of protostellar formation at 10 −6 g cm −3 and capturing a further ∼100 yr of disc fragmentation.We directly compare our results to a recent study, Prole et al. (2023) (hereafter LP23) which examined H 2 cooling minihalos with the same simulation code, chemical set-up and maximum resolution as the simulations presented in this work.
The format of this paper is as follows: in §2 we outline the numerical technique used including the simulation code and chemical network.In §3 we discuss the collapse of the gas up to point immediately prior to sink formation, while in §4 we discuss the fragmentation of the gas after the insertion of sink particles.In §5 we analyse the growth of the stellar system and compare, via a convergence study, the evolution of the star particles into main sequence stars and discuss their eventual evolution into (massive) BHs.In §6 we discuss some caveats before concluding in §7.

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

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

Simulation Setup
As discussed in the Introduction the goal of this study is to investigate the formation of primordial Pop III stars at the centre of a metal-free atomically cooling halo.Such halos have previously been investigated as ideal sites for heavy seed formation (Haiman et al. 1996;Regan & Haehnelt 2009;Latif et al. 2013b;Latif & Volonteri 2015;Latif et al. 2015Latif et al. , 2020;;Wise et al. 2019).
To generate the appropriate initial conditions we generated cosmological initial conditions using MUSIC (Hahn & Abel 2011).An initial cosmological simulation was performed within a comoving box of side length 1 h −1 Mpc using a ΛCDM cosmology with parameters h = 0.6774, Ω 0 = 0.3089, Ω b = 0.04864, Ω Λ = 0.6911, n = 0.96 and σ 8 = 0.8159 (Planck-Collaboration et al. 2020).The simulations were initialized at z = 127 with an initial dark matter distribution using the transfer functions of Eisenstein & Hu (1998).The gas distribution was set within Arepo to initially follow the dark matter (i.e.GENERATE_GAS_IN_ICS = 1).We model the dark matter with 512 3 particles and the gas was modelled with 512 3 grid cells (prior to refinement).During the simulation, an additional Jeans refinement criterion was applied such that the Jeans length of the gas is always resolved with at least 4 grid cells.Hence, the gas is able to dynamically refine during the simulation allowing maximum resolution where re- During this initial phase, we disabled the molecular chemistry functionality and hence utilised a simpler 6 species model such that H 2 and HD abundances remained at their initial value.This eliminated the need for a LW background radiation field at this initial stage of the calculation since the goal of this study is to investigate the idealised case of Pop III formation inside of a pristine, atomically cooling halo at a virial temperature of approximately 8000 K.We note that applying a super-critical LW field has the same result on the H 2 abundance.As the dominant HD formation pathway relies on the H 2 abundance (Nakamura & Umemura 2002), a super-critical LW field also suppresses the HD abundance.
At z ∼ 13.3 the first atomic halo begins to collapse within our box.Two more, physically distinct halos collapse at z = 11.5.These first three halos to begin atomically cooling and collapsing were then extracted.The selection criteria was simple; while H 2 rich halos begin cooling down from 1000 K from ∼10 −23 g cm −3 , reaching 200 K by ∼10 −22 g cm −3 , we have disabled H 2 formation, hence the only way for the gas to collapse to these densities is via atomic cooling.A density threshold for selecting a collapsing halo was therefore chosen as slightly higher than this density at 10 −21 g cm −3 .Only halos containing gas cells with a density exceeding this threshold were selected for extraction and resimulation.The central coordinate of the halos was found using a friends of friends (FoF) algorithm and a new box length of 2 kpc (physical units) was cut around it.Using the new box cut from the parent box, the simulations were restarted as "zoom-in" simulations.The units were converted into physical units for the remainder of the zoom-in calculation and the simulation essentially run as an isolated galaxy simulation with periodic boundary conditions.Table 1 shows the virial mass, radius and temperature of each halo as calculated by the FoF algorithm as well as the redshift when it was extracted, while Figure 1 shows the temperature-density profiles at that stage and Figure 2 shows the radial profiles of the enclosed mass for both the gas and the DM.
For the zoom-in simulations, the full chemistry model (12 species) is again enabled to get an accurate picture of H 2 formation at very high densities inside the core of the collapsing halo.We invoke a super-critical LW radiation field of J LW = 10 5 J 21 (where J 21 is in units of 10 −21 erg s −1 cm −2 Hz −1 sr −1 ) to suppress H 2 formation.We model the effects of H 2 self-shielding using the TreeCol algorithm (Clark et al. 2012).
We turn-on the 12 species model to examine the impact of H 2 production which is inevitable at the highest densities.In our highest resolution simulations, we evolve the collapse up to a density of 10 −6 g cm −3 before inserting sink particles (see §2.5) and capturing a further ∼100 yr of disc fragmentation and accretion behaviour, while in the lowest resolution simulation we capture ∼10 4 yr of accretion.See Table 2 for a list of simulation realisations and resolutions employed.The structure of Halo 3 at the end of the high resolution simulation is visualised in Figure 3, which shows the density field along with H 2 fraction and temperatures at various zoom-in scales.

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

Sink particles
The simulation mesh must be refined during a gravitational collapse to ensure the local Jeans length is resolved.Once the mesh refines down to its minimum cell volume, sink particles must be introduced to represent the dense gas, preventing artificial instability in cells larger than their Jeans length.Our sink particle implementation was introduced in Wollenberg et al. ( 2020) and Tress et al. (2020).A cell is converted into a sink particle if it satisfies three criteria: 1.The cell reaches a threshold density.2. It is sufficiently far away from pre-existing sink particles so that their accretion radii do not overlap.3. The gas occupying the region inside the sink is gravitationally bound and collapsing.
Likewise, for the sink particle to accrete mass from surrounding cells, the cell must meet two criteria: 1.The cell lies within the accretion radius.2. It is gravitationally bound to the sink particle.
A sink particle can accrete up to 90% of a cell's mass, above which the cell is removed and the total cell mass is transferred to the sink.Increasing the threshold density for sink particle creation drastically increases the degree of fragmentation, reducing the masses of subsequent secondary protostars (Prole et al. 2022a).However, increasing the sink particle threshold density also increases the computational challenge beyond which it is currently intractable.Ideally, sink particles would be introduced when the gas becomes adiabatic at ∼ 10 −4 g cm −3 (Omukai 2000).However, the zero metallicity protostellar model of Machida & Nakamura (2015) suggests that stellar feedback kicks in to halt collapse at ∼ 10 −6 g cm −3 (10 18 cm −3 ), so we choose this as our sink particle creation density for our highest density simulations.
The accretion radius of a sink particle R sink is chosen to be the Jeans length λ J corresponding to the sink particle creation Fig. 1: 2D histograms of the temperature -density profiles for the halos at the point at which they are extracted from the initial cosmological simulation, weighted by total gas mass within each 2D bin.The halos have begun to gravitationally collapse via atomic cooling as seen from the close to isothermal temperature profile.They differ from the regular H 2 minihalo case by the absence of a sharp drop in temperature to ∼200 K beginning at ∼ 10 −23 g cm −3 .Fig. 2: Radial profiles of enclosed gas and DM mass for the three atomic halos at the point where they are extracted from the initial cosmological simulation.The baryonic component becomes comparable to the DM within the central ∼ 100 pc where it begins to decouple from the DM.In Halo 3, the baryonic component is dominant over DM on these scales.density and corresponding temperature.Taking our high resolution simulation as an example, at a density of 10 −6 g cm −3 , we take the temperature value from Prole et al. (2022a) of 4460 K to give a Jeans length of 1.67 × 10 12 cm (0.11 au).We set the minimum cell length to fit 8 cells across the sink particle accretion radius in compliance with the Truelove condition, by setting a minimum cell volume V min = (R sink /4) 3 .The minimum gravitational softening length for cells and sink particles L soft is set to R sink /2.The simulation parameters for the low, medium and high resolution simulations are summarised in Table 2.
The sink particle treatment also includes the accretion luminosity feedback from Smith et al. (2011), as implemented in Arepo by Wollenberg et al. (2020).Stellar internal luminosity is not included in this work, which is not a problem in our high resolution simulations because the Kelvin-Helmholtz times of the protostars formed are much longer than the period simulated, meaning that no stars will have yet begun nuclear burning.This however would affect our lower resolution simulations, which run for significantly longer periods, although the sink particles in those simulations (potentially) represent groups of protostars rather than individual stars.A comprehensive treatment of protostellar feedback is nonetheless outside the scope of this study.
Lastly, as protostellar mergers are often reported in primordial star forming simulations (e.g.Greif et al. 2012;Hirano & Bromm 2017; Susa 2019), we include the treatment of sink par-  Since sink particles carry no thermal data, the last criteria simply require that their gravitational potential exceeds the kinetic energy of the system.When these criteria are met, the larger of the sinks gains the mass and linear momentum of a smaller sink, and its position is shifted to the centre of mass of the system.We allow multiple mergers per time-step.For example, if sink A is flagged to merge into sink B, and sink B is flagged to merge into sink C, then both A and B will be merged into sink C simultaneously.

Initial collapse
Figure 4 summarises the state of the collapse at a point just before the formation of the first sink particle.The temperaturedensity relationship in the top panel shows that the initial contraction up to densities of 10 −15 g cm −3 follows a near isothermal collapse as expected (e.g.Omukai 2000;Klessen & Glover 2023).Above this density, the temperature drop steepens, but remains close to isothermal.Here the temperature drops by around a factor of 2 over a density range of 10 4 , giving an effective gamma of ∼ 0.95.This steepening of the temperature profile was first reported on in the one-zone calculations of Omukai (2001); at this density, the H 2 abundance is still too low to provide significant cooling, instead the temperature drop corresponds to the point where cooling becomes dominated by H − free-bound cooling.The near-isothermal collapse of the gas is further evidenced in the bottom panel of Figure 4 where we clearly see a ρ ∝ R −2 scaling of the density versus radius over more than six orders of magnitude in radius.
From the very bottom panel of Figure 4, the H 2 abundance begins to build up, albeit from very low abundance levels, within the inner 10 −2 pc (∼ 2000 au) of the halo.While our LW field strength of J 21 = 10 5 is already quite extreme, we run a second realisation of Halo 3 with an extremely high value of 10 10 (as shown in the top panel of Figure 5) to demonstrate that above roughly 10 −15 g cm −3 LW photo-dissociation can not prevent an increase in the H 2 abundance (nor probably can any other physical process).The three-body H 2 formation rate per unit volume is proportional to n 3 (where n is number density), whereas the corresponding scalings for the photodissociation rate and the H 2 collisional dissociation rate are n and n 2 , respectively.Therefore H 2 formation will inevitably overcome its destruction at high densities.In practice then the formation of H 2 at high densities is inevitable.
Furthermore, a build in the H 2 fraction at high densities owes to the exponential nature of the self-shielding.To demonstrate this, we calculate the H 2 shielding factor given in Wolcott-Green et al. (2011) as The shielding factor is given as a function of density in the bottom panel of Figure 5, which quickly falls to 0 above densities of 10 −15 g cm −3 , independently of the LW intensity.This shows that the core will always be shielded from LW radiation above these densities regardless of the value of J 21 .A build up of H 2 fraction within the dense core has been seen to various degrees in the radial profiles of previous studies.For example, the one-zone calculations of Omukai (2001) show a H 2 fraction of ∼ 0.01 at a density of ∼ 10 −8 g cm −3 , while simulations by Becerra et al. (2015) reach a H 2 fraction of ∼ 0.1 by ∼ 10 −6 g cm −3 and Regan et al. ( 2017) reaches a H 2 fraction of 10 −4 by their maximum density of 10 −15 g cm −3 .What the results here show is that the very central regions of atomic cooling halos are highly likely to have small pockets of fully molecular gas.

Fragmentation and the Initial Mass Function
We now move onto to discuss the build of the initial mass function and the resulting protostellar masses that are found in our simulations at the very highest densities when we introduce sink particles.In Figure 6 we show the mass flux into concentric shells just before the formation of the first sink particle, comparing the standard Pop III star forming minihalos (∼ 10 6 M ⊙ ) of LP23 with the more massive, atomically cooled halos simulated here.Within the inner ∼100 pc, the accretion rate into the atomic halos exceeds the minihalo case by a factor of between 10-1000 due to the combined effects of a stronger DM gravitational well and a larger available reservoir of baryonic gas.In principle this should lead to the formation of more massive sinks given the higher accretion rates possible.
Figure 7 shows how the total mass in sink particles, the mass of the most massive sink particle and the total number of sinks evolves in time.The total mass in sinks of the system exceeds the minihalo case by almost 2 order of magnitude due to the higher mass flux.While the initial accretion onto the most massive sink particle far exceeds the same quantity in the minihalos, the most massive sink is ejected from the system in all three atomic halos within the first ∼50 yr after its formation.This results in the most massive sink particle aligning with the upper limit seen in the minihalos.The most surprising result is that the fragmentation of the gas in the center of the atomic halos far exceeds the minihalo cases.
As discussed earlier, H − formation becomes the dominant cooling process above 10 −15 g cm −3 , leading to a slightly steeper decline in temperature with density.However, the temperature is still higher in the atomically cooled halos when compared to the minihalo case at the same density (see e.g.Omukai 2000;Yoshida et al. 2006;Prole et al. 2022a), corresponding to a higher Jeans mass.The increase in fragmentation is therefore attributable to the higher mass infall rate.For example, assuming fragmentation is perfectly efficient, the maximum number of fragments that can form is the number of Jeans masses present within the disc.We can therefore get an estimate of the fragmentation from Figure 8, which shows the number of enclosed Jeans masses of gas at or above a given density.At high densities the number of Jeans masses present in the atomically cooled halos exceeds the minihalo case by an order of magnitude and roughly matches the number of sink particles formed in each halo (see Figure 7), which explains the increase in fragmentation when compared to the minihalos.
The combined sink particle mass functions from all three halos are shown in Figure 9 at ∼ 100 yr after the formation of the first sink particle.We also show the mass functions from the minihalos at the end of the LP23 simulations (∼300 yr).Despite increased fragmentation and a factor of 3 less in time to accrete, the atomic halo mass function peaks at a higher protostellar mass of 3 M ⊙ compared to the 0.3 M ⊙ peak in the minihalos.We note that these protostars will continue to accrete for roughly 10 4 yr before the end of the pre-main sequence.While the zero- age main sequence (ZAMS) masses of these stars are unknown (see §5), we have confirmed that atomically cooled halos do produce a population of higher mass protostars when compared to regular Pop III star forming minihalos, at least after the first 100 yr after the formation of the first protostar.These protostars can then accrete the available gas in competition with further star formation.

Final stellar and black hole masses
In order to estimate the ZAMS masses of the sinks we run lower resolution simulations of Halo 3 which allows us to evolve the system over longer timescales.The top panel of Figure 10 shows how the total mass of the system of sink particles grows as a function of time across the different resolutions tested.We show the formation of new sink particles in the 10 −13 and 10 −10 g cm −3 simulations as star shaped markers.We do not show the forma-10 22 10 19 10 16 10 13 10 10 10 7 [g cm 3 ] The Jeans mass at each density bin is calculated using the mass weighted average temperature and density within the bin.
The calculated Jeans mass is then compared to the total mass of gas at or above the density of the bin.tion of sinks for the 10 −6 g cm −3 cases as fragmentation occurs almost immediately.Clearly the higher the resolution used, the earlier fragmentation occurs.If we take t = 0 to be the time at which the first sink forms in each case then the second sink forms almost immediately in the highest resolution case, after ∼ 10 yr in the 10 −10 g cm −3 case and only after ∼ 1000 yr in the 10 −13 g cm −3 case.As the resolution used does not significantly affect the growth of the total system, the lowest resolution run shows that the system will continue to grow linearly through the pre-main sequence phase to reach a mass of ∼ 10 4 M ⊙ by 10 4 yr i.e. there will be ∼ 10 4 M ⊙ available for star formation within the central 264 au.The increased fragmentation in the higher resolution simulations complicates how this mass will be distributed amongst the protostars.The bottom panel of Figure 10 shows the accretion rate onto the most massive surviving sink particle i.e. the most massive non-ejected protostar.At the end of the high resolution simulations, the largest survivors The times at which new sink particles formed are indicated with star-shaped markers in the medium and low resolution simulations.We do not show the formation of sink particles for the high resolution simulations as fragmentation occurs almost instantly.Bottom -accretion rates onto the most massive surviving (non-ejected) sink particle.Also shown are the masses of the most massive surviving sink at the end of the simulations.
in the three halos have accretion rates in the range 10 −3 -10 −1 M ⊙ yr −1 and have masses in the range of 7.5-12 M ⊙ .Among the two lower resolution simulations the final masses of the protostars are approximately 1500 M ⊙ (after 1000 years) and 25000 M ⊙ (after 20,000 years) respectively.If the most massive sink particle from the highest resolution simulations can accrete or grow at or near the rates found in the lower resolution simulations then super-massive star formation may be realised.
The zero-age main sequence (ZAMS) mass of a star depends on the pathway between the formation of the protostar and the eventual beginning of core hydrogen burning, which in turn depends strongly on the evolution of the protostellar accretion rate.When a protostar's Kelvin-Helmholtz (KH) timescale (time to radiate away its own gravitational energy) is shorter than its accretion timescale (time to double its mass by accretion), it will radiate away its energy and begin to contract, which typically occurs at masses ∼10 M ⊙ (Palla & Stahler 1991;Omukai & Palla 2003) though this depends sensitively on the assumed accretion rate (e.g.Nandal et al. 2023).The contraction causes an increased extreme ultraviolet (EUV) luminosity and surface temperature, ionizing infalling gas and accelerating it outwards.This triggers a runaway collapse as the decreased accretion rate causes further contraction until hydrogen burning begins and the protostar reaches the ZAMS.If the accretion rate onto the largest sink particle in our high resolution simulations falls and remains below ∼ 4 × 10 −3 M ⊙ yr −1 , the KH contraction will begin at ∼10 M ⊙ and the growth of the protostar will be self-regulated by these radiative feedback effects, limiting the final mass to a few tens of M ⊙ (Hosokawa et al. 2011(Hosokawa et al. , 2012b)).However, numerous studies have examined the impact that accretion has on the contraction of a Pop III star (e.g.Hirano et al. 2014).A key quantity here is the critical accretion rate -this is the accretion rate onto the stellar surface that prevents contraction of the star.Recently Nandal et al. (2023) investigated in detail this quantity in terms of episodic accretion rates using the Geneva Stellar Evolution code (Genec: Eggenberger et al. 2008) .They found that during the pre-main sequence, which is of most relevance here, the critical accretion rate is ∼ 2.5 × 10 −2 M ⊙ yr −1 .In this case the effective surface temperature remains below 10 4 K, hence feedback can not form a HII region and accretion is not prevented, leading to the formation of a super-giant protostar, which can grow up to several thousand solar masses depending on the details of the accretion (Umeda et al. 2016;Woods et al. 2017;Haemmerlé et al. 2018).The end of accretion onto a super-giant protostar is caused by a fast contraction when the accretion rate falls below ∼ 7 × 10 −3 M ⊙ yr −1 .
At the end of the high resolution simulations, the accretion rates onto the largest surviving sink particles vary significantly in time.As the accretion rates onto ∼10 M ⊙ protostars appear to be a good indicator of future contraction/accretion evolution (Hosokawa et al. 2011(Hosokawa et al. , 2012b;;Hirano & Bromm 2017) and the largest survivor in Halos 2 and 3 have accretion rates between 10 −2 -10 −1 M ⊙ yr −1 , they could feasibly go on to form (supermassive) Pop III stars with masses anywhere in the range of 10 M ⊙ up to 10 4 M ⊙ within their main sequence (MS) lifetime.If they grow in excess of 25 M ⊙ the type II supernova explosion would be too weak to eject much of the star and the subsequent fallback of material causes the resulting neutron star to collapse into a BH (MacFadyen et al. 2001), while above 40 M ⊙ the neutron star is unable to launch a supernova shock and collapses directly to form a BH with no mass loss (Fryer 1999).If these stars avoid the disruptive pair instability supernova mass range of 140-260 M ⊙ , they will result in a massive BH equal in mass to the stellar progenitor.However, whether the protostars will maintain their accretion rates is uncertain.In the optimistic case that they maintain accretion rates of 10 −1 M ⊙ yr −1 , these results are in line with many previous, lower resolution simulations of heavy seed black hole formation (e.g.Johnson et al. 2011;Latif et al. 2013a,b;Regan et al. 2014;Latif & Volonteri 2015;Choi et al. 2015;Regan et al. 2017;Smidt et al. 2017;Ardaneh et al. 2018) despite the increased gas fragmentation.The largest survivor in Halo 1 ends with an accretion rate of ∼10 −3 M ⊙ yr −1 , which if maintained will likely lead to an early KH contraction and limit the final mass to 10-30 M ⊙ .It is unclear if the accretion rate onto Halo 1 will remain below 10 −3 or if it will experience an increase similar to that of Halo 3. Certainly the formation of massive BH seeds of 10 4 M ⊙ needed to explain z ∼ 7 quasar observations would rely on frequent mergers with other protostars.
The sink particle masses in the lower resolution simulations represent whole groups of protostars in the high resolution simulations.Frequent mergers of secondary protostars into the main protostar would push the stellar and later BH masses up to close to what is achieved in the low resolution simulations, with an upper limit of 10 4 M ⊙ within a 10 kyr period.Recent studies have pointed out the importance of super competitive accretion, in which a central few stars grow supermassive while a large number of other stars are competing for the gas reservoir (Chon & Omukai 2020), with central objects growing significantly through mergers (Sassano et al. 2021;Vergara et al. 2021;Trinca et al. 2022;Schleicher et al. 2022;Zwick et   where M 0 is the initial mass of a sink particle.Bottom -ratio of cumulative number of sink mergers against the current number of surviving sink particles.

2023).
Here the mass growth by collisions can be comparable to its growth via accretion (Schleicher et al. 2023).The supercompetitive accretion scenario was initially invoked to explore the 'low metallicity' regime (Z ∼ 10 −6 − 10 −3 Z ⊙ ) where fragmentation is expected to be more active compared to the metalfree case.What we find here is that fragmentation is already vigorous within the central 2000 au of the halo.It therefore appears that for all metallicities below ∼ 10 −3 Z ⊙ we can expect a scenario where vigorous fragmentation competes with a rapidly growing protostar.Differentiating the environments which produce a dense stellar cluster versus those that produce a central massive Pop III star may be the next frontier.
To investigate the merger behaviour further, we plot the total number of sink mergers, ratio of total merged mass to accreted sink mass and the ratio of the number of mergers to the number of surviving sink particles as a function of time in Figure 11.From this we see that the mass gained through mergers constitutes between 1-10% of the total sink particle mass.The frequency of mergers is high, with the number of mergers totalling between 10-50% of the total number of surviving sink particles.If this holds throughout the following 10 4 yr, the total mass gained through mergers alone could be as high as 10 3 M ⊙ , although it is unclear if this mass would be shared between the growing number of protostars randomly or be preferentially received by the largest growing BH seed.

Caveats
We have not included magnetic fields in these simulations.While studies of primordial magnetic fields suggest that they can increase the mass of protostars (e.g.Saad et al. 2022;Hirano & Machida 2022;Stacy et al. 2022), the fields have no effect when they are properly resolved, distributing the magnetic energy from the small-scale turbulent dynamo to smaller spatial scales (Prole et al. 2022b).
We have also not included radiative feedback from our protostars, which has the effect of heating the surrounding gas and lowering accretion rates.While our high resolution simulations end at a point when the protostars would not produce significant levels of feedback.The low resolution simulations however run for much longer and would be subject to feedback effects.
As our high resolution simulations have only captured the first ∼ 100 yr of accretion after the formation of the first protostar and the lower resolution simulations lack the capacity to resolve individual star foratmion, to what extent competitive accretion and mergers affect the subsequent final stellar/BH masses is unclear.
Ideally we would resolve up to the protostellar formation density of 10 −4 g cm −3 , which is currently computationally unfeasible.Failure to resolve up to that density means we have not achieved numerical convergence.Despite this, our implemented threshold density of 10 −6 g cm −3 and post-sink particle run time of 100 yr represents the current state of the art in this regime.
As mentioned in §2.2, although we account self-consistently for the effects of opacity when computing the H 2 cooling rate (due to line cooling and collision-induced emission), we do not yet account for the continuum opacity of the gas when computing the cooling rate due to other radiative processes, e.g.H − freebound cooling.At most densities in our models, the optical depth of the gas in the continuum is small and this simplification has little impact -for example, it should be valid at all of the densities traced in our medium and low resolution models.However, above a density of ∼ 10 −8 g cm −3 , we expect the optical depth of the gas in the continuum to become significant, and hence our current simulations over-estimate the cooling rate of the gas above this density.This will likely give us slightly more effective fragmentation at the highest densities than we would find in reality.However, we do not expect this to significantly impact the main conclusions of our study.As 8 demonstrates, the number of Jeans masses present in the gas at ρ = 10 −8 g cm −3 in the atomic cooling halos already substantially exceeds the corresponding value in H 2 -cooled minihalos, demonstrating that this result is not a consequence of our simplified treatment of very high density cooling.Our finding that the gas fragments extensively in the centre of the atomic-cooled halos should therefore be a robust result, although some of the details (e.g. the precise number of fragments and their formation masses) will depend on the treatment of the cooling at very high densities.Our finding that the fragments quickly grow to become more massive than their Pop III counterparts should also be a robust result: properly accounting for the continuum opacity will likely give us fragments with slightly larger initial masses, which if anything will increase the rate at which they later gain mass, exacerbating the difference between the atomic-cooled and H 2 -cooled results.
We have assumed that the background radiation field from Pop III stars takes the form of a black body spectrum with an effective temperature of 10 5 K, peaking in the UV.However, this relies on the stars having masses of ∼ 300 M ⊙ (Schaerer 2002).The masses of Pop III stars are currently uncertain, so our choice of effective temperature affects the photodissociation of the molecules present in the gas.For example, the effective temperature for a 1 M ⊙ Pop III star is ∼ 7000 K (Larkin et al. 2022), which peaks at lower frequencies, resulting in significantly less H 2 photodissociating LW radiation but more infrared radiation, which can still photodissociate H − to disrupt H 2 formation (Chuzhoy et al. 2007).
The methodology here was idealised in that we used a six species chemical model to construct a pristine atomic cooling halo before switching to a full 12 species chemical model and employing an intense LW background.In future work we will apply the same methodology to rapidly assembling halos without the reliance on super-critical LW values.

Conclusions
We have performed high resolution simulations of three pristine atomically cooled halos of mass ∼ 10 8 M ⊙ which begin to collapse at z ∼ 11.5 − 13.3.The goal of our simulation suite was to examine the fragmentation and protostellar formation within the core of such systems as they reach protostellar densities.For the highest resolution we could tractably simulate (ρ ∼ 10 −6 g cm −3 ) we followed the protostar formation for approximately 100 years by introducing sink particles representing individual protostars.The main results are as follows: 1. Resolving protstellar densities reveals that the center of (pristine) atomically cooled halos is subject to intense fragmentation, even more-so than in the canonical Pop III minihalo case.This is driven by H − free-bound cooling and an increased number of Jeans masses within the core, owing to enhanced accretion rates.2. Despite increased fragmentation, the atomically cooled halos formed a population of higher mass protostars compared to Pop III star forming minihalos.3. The initial accretion rates onto the most massive surviving protostars indicates that their zero-age main sequence masses could range from 10 2 − 10 4 M ⊙ depending on subsequent accretion and/or mergers.4. Our coarser resolution simulations show that the total mass of the system continues to grow steadily for at least 10 4 years after the formation of the first protostar (achieving central sink particle masses in excess of 10 4 M ⊙ ), although how this is distributed amongst the rapidly growing number of protostars is unknown and requires further study.5.The formation of a massive Pop III star (M * ≳ 1000 M ⊙ ) is therefore realistic and achievable in a high-z setting but relies on (super-competitive) accretion and frequent mergers with secondary protostars 6.The H 2 fraction within the inner ∼2000 au of all three halos was able to build up independently of the strength of the external LW radiation field.This strongly indicates that all collapsing halos (even those subject to strong LW fields and likely other physical processes which could suppress star formation up to the atomic limit e.g.streaming velocities between DM and baryonic gas) contain small pockets of cold, H 2 rich gas at their centre.
The findings here show that the formation of massive Pop III stars (i.e.heavy seed MBHs) is entirely plausible in sufficiently massive halos but will depend sensitively on the future evolution of the protostars and on the balance between stellar accretion and mergers in the dense cluster that forms within the center of the atomic cooling halo.While we model the impact of a supercritical LW radiation field here, any physical process which results in the initial suppression of star formation up to the atomic cooling limit should have a similar impact and allow for the formation of more massive Pop III stars.

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

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

Fig. 5 :
Fig. 5: Top -Comparison of the build up of H 2 in high density regions when a J 21 = 10 5 and 10 10 LW radiation field are used.Bottom -H 2 shielding factor calculated using equation 12 from Wolcott-Green et al. (2011).

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

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

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

Fig. 9 :
Fig.9: Comparison of the sink particle mass function from the 3 atomic cooling halos at ∼100 yr versus the 15 H 2 cooling minihalos from LP23 at ∼300 yr.Power laws of M 0.85 and M −2 are superimposed to give the reader an idea of the slopes involved.

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

Fig. 11 :
Fig.11: Time evolution of the merger data.Top -cumulative number of sink particle mergers.Middle -ratio of cumulative merged mass to the total accreted sink particle mass (M tot −M 0 ) where M 0 is the initial mass of a sink particle.Bottom -ratio of cumulative number of sink mergers against the current number of surviving sink particles.

Table 1 :
Halo virial masses, radii and temperatures along with redshifts at the point where they are extracted from the initial cosmological simulation.

Table 2 :
Simulation parameters.Left to right -sink particle creation density, accretion radius, gas minimum cell volume and minimum gravitational softening length, DM mass resolution and gravitational softening length.resolution ρ sink [g cm −3 ] R sink [au] V min [au 3 ] L soft [au] M DM [M ⊙ ] L DM [pc] al.