Modeling Supermassive Primordial Stars with MESA

Supermassive stars forming at $z \sim$ 15 - 20 are one of the leading contenders for the origin of the first quasars, over 200 of which have now been discovered at $z>$ 6. These stars likely form in pristine, atomically cooled haloes immersed in strong Lyman-Werner UV backgrounds or in highly supersonic baryon streaming flows. Atomic cooling triggers catastrophic baryon collapse capable of building up stars at rates of up to $\sim$1 M$_{\odot}$ yr$^{-1}$. Here we examine the evolution of supermassive stars with a much larger and finer grid of accretion rates than in previous studies with the \texttt{MESA} stellar evolution code. We find that their final masses range from 3.5 $\times$ 10$^3$ M$_{\odot}$ - 3.7 $\times$ 10$^5$ M$_{\odot}$ at accretion rates of 0.001 M$_{\odot}$ yr$^{-1}$ - 1 M$_{\odot}$ yr$^{-1}$, respectively. We also find that supermassive star evolution diverges at accretion rates of 0.01 M$_{\odot}$ yr$^{-1}$ - 0.02 M$_{\odot}$ yr$^{-1}$, above which they evolve as cool red hypergiants along the Hayashi track and collapse via the general relativistic instability during central hydrogen burning, and below which they evolve as hot blue supergiants and collapse at the end of their nuclear burning lifetimes after exiting the main sequence.


INTRODUCTION
Supermassive stars (SMSs) are one of the leading candidates for the origin of the first quasars, more than 200 of which have been discovered at > 6 ( Fan et al. 2003Fan et al. , 2006, including eight at > 7 ( Mortlock et al. 2011;Wu et al. 2015;Bañados et al. 2018;Smidt et al. 2018;Matsuoka et al. 2019;Zhu et al. 2020;Wang et al. 2021). Until recently, they were thought to form at ∼ 10 -20 in primordial halos immersed in high Lyman-Werner UV backgrounds (Latif et al. 2014b;Agarwal et al. 2016;Schauer et al. 2015Schauer et al. , 2017a or in highly supersonic baryon streaming motions (BSMs; Latif et al. 2014a;Hirano et al. 2017;Schauer et al. 2017b) that suppress normal star formation until they reach masses of 10 7 M and virial temperatures of 10 4 K. These temperatures trigger atomic cooling that causes gas to collapse to the center of the halo at rates of up to ∼ 1 M yr −1 (Bromm & Loeb 2003;Lodato & Natarajan 2006;Regan & Haehnelt 2009;Latif & Volonteri 2015). Such flows create massive, hot accretion disks that build up either a single SMS or binaries and small multiples (Wise et al. 2019;Regan et al. 2020;Latif et al. 2020;Patrick et al. 2021Patrick et al. , 2020. However, it has now been found that the rare haloes capable of forming quasars by ∼ ★ E-mail: nh478@exeter.ac.uk 6 can form SMSs without the need for UV backgrounds, BSMs, or even atomic cooling (Latif et al. 2022a) SMSs have been the subject of analytical studies since the 1960s (e.g., Iben 1963;Chandrasekhar 1964;Fowler 1964Fowler , 1966 and numerical simulations since the 1970s (e.g., Appenzeller & Fricke 1972;Shapiro & Teukolsky 1979;Fuller et al. 1986;Baumgarte & Shapiro 1999;Sun et al. 2018;Butler et al. 2018), but it has only been recently that they have been studied in the extreme flows in which they form (Hosokawa et al. 2013;Umeda et al. 2016;Woods et al. 2017;Haemmerlé et al. 2018a;Woods et al. 2021a,b). These models indicate that rapidly accreting SMSs can reach masses of a few 10 5 M and lifetimes of 1 -2 Myr before collapsing to directcollapse black holes (DCBHs). SMSs are the leading contenders for the seeds of the first quasars because it is difficult for ordinary Population III (Pop III) star BHs to grow rapidly after birth because they form in low densities and in some cases are ejected from their host halos (Whalen et al. 2004;Kitayama et al. 2004;Alvarez et al. 2009;Smith et al. 2018). DCBHs are born with much larger masses and in much higher densities in halos that retain their fuel supply, even when it is heated by X-rays (Johnson et al. 2013a).
with surface temperatures of 5,000 -10,000 K due to H − opacity in their atmospheres, at least until they reach ∼ 10 5 M (Hosokawa et al. 2013). Haemmerlé et al. (2018a) found that SMSs can remain cool even above these masses, reaching luminosities greater than 10 10 L . In other studies, SMSs evolving from similar initial conditions soon settle onto hotter and bluer tracks with temperatures of 20,000 -40,000 K (Woods et al. 2017). Haemmerlé et al. (2018a) found that stars growing at low rates ( 0.005 M yr −1 ) also evolve along blue tracks, as can stars with clumpy accretion due to fragmentation or turbulence in the disk (Sakurai et al. 2015 -but see Sakurai et al. 2016). It is not clear if these differences arise from the opacities used in the models, code physics (such as the numerical treatment of convection), or resolution (see also Bear & Soker 2020 for additional discussion of convergence issues and Woods et al. 2019 for further discussion of supermassive star models).
Recent work suggests the ΩΓ limit, in which radiation pressure and centrifugal forces equal gravity in the star, restricted SMS rotational velocities to at most a few percent of Keplerian (Haemmerlé et al. 2018b). Consequently, the flows that created such stars in the early Universe must have had efficient mechanisms for shedding angular momentum. One possibility is the rise of 'bars within bars' instabilities on small scales (e.g., Wise et al. 2008). Another is the rise of magnetic fields in the accretion disk of the star, first on small scales due to turbulent subgrid dynamos (Schober et al. 2012) and then their subsequent amplification via the Γ instability on larger scales (Latif & Schleicher 2016).
The relatively low surface temperatures of SMSs in most models to date imply that radiation from the star cannot overcome the enormous ram pressures of catastrophic infall and that they form, evolve and die in the accretion flows that create them. This has been corroborated at the initial stages of SMS formation by radiation hydrodynamical simulations by Ardaneh et al. (2018) and Luo et al. (2018) (see also Smith et al. 2017). Any ionizing UV from the star would also only heat the gas to at most a few times the temperatures at which it already falls onto the star without forming the usual hot, high-pressure bubbles that drive all the baryons out of less massive Pop III star-forming halos (Susa et al. 2009;Susa 2013;Hirano et al. 2014Hirano et al. , 2015Sugimura et al. 2020;Latif et al. 2021Latif et al. , 2022b. Pop III SMSs accreting at rates of 0.1 -1 M yr −1 would have extremely large luminosities that could, in spite of their low surface temperatures, be detected in the near infrared (NIR) by the James Webb Space Telescope (JWST) and large ground-based telescopes in the coming decade (Johnson et al. 2012;Hosokawa et al. 2013;Surace et al. 2018Surace et al. , 2019Whalen et al. 2020a). Their BHs could also be found in the NIR at ∼ 15 − 20 by JWST (Pacucci et al. 2015;Natarajan et al. 2017;Barrow et al. 2018;Whalen et al. 2020b) and at ∼ 6 -8 by Euclid and the Roman Space Telescope, although lensing by massive galaxies and galaxy clusters in their wide fields could extend these detections up to ∼ 15 (Vikaeus et al. 2022). DCBHs will only be marginally visible to the Square Kilometre Array or next-generation Very Large Array at 6 -8 (Whalen et al. 2020a but would become more luminous after growing to larger masses at later times (Whalen & Mezcua 2022). They could also be detected out to ∼ 10 by future X-ray missions such as the Advanced Telescope for High-Energy Astrophysics (ATHENA) and Lynx (Aird et al. 2013;The Lynx Team 2018).
Previous studies of supermassive primordial star evolution considered small representative samples of 4 -7 accretion rates in the range of 0.001 M yr −1 -10 M yr −1 . Here we revisit the evolution of supermassive primordial stars with the Modules for Experiments in Stellar Astrophysics (MESA) code by examining a much larger and finer grid of accretion rates than previously explored. We describe the physics in our runs and their setups in Section 2. The evolution of the stars is examined in detail in Section 3. We tally final masses for SMSs as a function of accretion rate in Section 4 and conclude in Section 5.

NUMERICAL METHOD
MESA is a one-dimensional (1D) Lagrangian stellar evolution code that solves equations of stellar structure that are implicitly coupled to convective mixing and nuclear burning (version 12778; Paxton et al. 2011Paxton et al. , 2013Paxton et al. , 2018. MESA has implicit hydrodynamics that can solve for the velocity of each zone in the model at the onset of collapse. We use the 21-isotope APPROX network, which includes the pp chain, triple-alpha burning, and the CNO cycle. Our models have an equation of state (EOS) that is a composite of several datasets such as the OPAL/SCVH tables (Rogers & Nayfonov 2002;Saumon et al. 1995), which are used at lower densities and temperatures like those in the outer regions of the star and its atmosphere, and the HELM and PC EOSs (Timmes & Swesty 2000;Potekhin & Chabrier 2010), which are applied to high densities and temperatures like those in the core of the star (see Figure 1 and Section 4.2 of Paxton et al. 2011). We use the Henyey method for convective mixing length theory, with the mixing length parameter MLT = 2, and an overshooting parameter set to 0.1. These two values are the defaults for massive stars in MESA and are consistent with those used in previous SMS simulations (e.g., Woods et al. 2017). We also include superadiabatic convection in radiative regions when the default criteria for it in MESA are met. We also use the Ledoux criterion for convection along with the MESA prescription for smoothing composition gradients in non-convective regions that trail retreating convection zones. This measure reduces unphysical steps in entropy that can develop in the interior of the star over time.
The effects of general relativity (GR) on the structure of the star are approximated by the Tolman-Oppenheimer-Volkoff (TOV) correction to the equation of hydrostatic equilibrium, which is implemented as a post-Newtonian correction to the gravitational constant, , (2) rel is computed in every zone of the star every time step in the user-defined run_star_extras subroutine in the MESA_star subroutine and then used to update Equation 1. We turn on post-Newtonian corrections to at the beginning of each run.
MESA adaptively rezones the stars as they evolve. Mesh refinement is triggered when gradients in temperature, pressure, and 4 He abundances between adjacent zones exceed preset values: log / > 1/30, log / > 1/80, and log( + 0.01) > 1/20 log , where is the 4 He mass fraction (see section 6.4 of Paxton et al. 2011). These criteria result in a larger number of zones usually being allocated near the center of the star to resolve nuclear burning and convective processes, and the stars are typically partitioned into about 1400 mass zones. Accretion flows onto the star have a primordial composition that does not change as the star evolves. Their entropy is matched to that of the surface of the star so it is assumed that the accretion luminosity is radiated away rather than deposited at the surface. We exclude mass loss due to stellar winds because they are thought to be negligible at primordial compositions (Vink et al. 2001;Baraffe et al. 2001). We initialise all the stars as 20 M fully-convective, = 1.5 polytropes with temperatures below 10 6 K to preclude nuclear burning (Section 6.1 of Paxton et al. 2011). They have primordial compositions, with mass fractions H = 0.7516 and He = 0.2484. Our 16 models are uniformly partitioned in log accretion rate over 3 decades: 10 −3 M yr −1 , 10 −2 M yr −1 and 10 −1 M yr −1 .
Models with accretion rates at or above 0.02 M yr −1 are first evolved up to 10,000 M without hydrodynamics because it can lead to numerical instabilities if the core is still contracting or burning is about to begin. Each model is then branched into two parallel runs: with and without hydrodynamics. Stars evolved with hydrodynamics eventually collapse because of pulsations induced by the GRI. Stars without hydrodynamics cannot actually collapse because there are no fluid motions so it is instead inferred from the 'softening' of the EOS in the core of the star, when the adiabatic exponent falls below the critical value for radiation-dominated stars when GR is taken into account: where is the mean molecular weight, and / B is the entropy per baryon (Woods et al. 2017). We run both cases for each star to determine if the use of the softening of the EOS as a proxy for collapse produced accurate final stellar masses in previous studies. Stars with hydrodynamics are evolved up to the mass at which their cores collapse irreversibly at velocities of 500 -1000 km s −1 and those without hydrodynamics are evolved until Γ falls below the critical value in the core. GR is turned on in these runs at all times.
Stars with accretion rates below 0.02 M yr −1 never reach masses at which they encounter post-Newtonian instabilities so GR is turned off in these models, but hydrodynamics is turned on from the beginning of the run in zones whose temperatures exceed 10 7 K. Numerical instabilities that arise during the contraction of the core and the expansion of the star when it becomes fully convective after the depletion of central hydrogen make it difficult to evolve the stars to more advanced stages of burning so we halt their evolution at central H fractions of 1% (He burning due to the triple process begins well before this point). However, our main objective is to determine the final mass of the star, and at these low accretion rates the star cannot gain much mass in the short times over which late stages of burning proceed. We get to within a few percent of the true final mass of the star even though we exclude these stages. We also activated a control that adds more resolution to the models when log > 8.15, by which point the central hydrogen fraction in most of the stars has fallen to below 30%.
Pressure must be imposed on the outer boundaries of these stars to prevent superluminal velocities from arising in their lowdensity outermost zones and ensure the stability of the star (Woods et al. 2017. This pressure, , is parametrised by extra_factor , and = / 2 , the surface gravity of the star, and are the opacity and optical depth, and and are the mass interior to those layers and their luminosity, respectively. We set extra_factor = 1.36 for our low accretion rate models, 0.01 -0.001 M yr −1 .

SMS EVOLUTION
We show Hertzsprung-Russel (HR) plots of the evolution of our stars at accretion rates of 0.001 M yr −1 -1 M yr −1 in Figures 1  and 2 and Kippenhahn diagrams of their structures over time in Figures 3 -5. How the stars evolve primarily depends on how accretion timescales, acc = / , compare to Kelvin-Helmholz (KH) contraction times, KH ≈ 2 / , early in their lives where and are the accretion rate and luminosity of the star, respectively. If the protostar gains enough mass during contraction it will have sufficient opacity in the outer layers to expand as a cool red SMS. If not, it will remain compact and hot and evolve along bluer tracks. At first, our protostars are fully convective, with central temperatures below those required for deuterium burning. They contract under accretion until their temperatures are high enough to ignite proton-proton (PP) burning.
However, as shown by the early dips in all the HR tracks in Figure 1, PP burning does not generate enough energy to halt contraction, which continues until the stars begin CNO burning, which is quickly followed by the triple alpha process. The onset of the CNO cycle thus marks the beginning of stable nuclear burning in the star. The PP chain raises central temperatures that drive changes in internal opacities that create a short-lived radiative core. The mass at which the star forms a radiative core increases with accretion rate because the core acquires more mass during initial contraction. However, CNO burning dramatically increases energy production so the core becomes convective to transport this energy outward more efficiently. The stars eventually develop large convective cores and high-entropy outer envelopes like those in Woods et al. (2017) and Haemmerlé et al. (2018a).
As in Haemmerlé et al. (2018a), each star converges to a monotonic mass-luminosity evolution after an early period of reconfigu- ration that, as discussed above, follows either a cool red track or a hot blue one. As shown in Figures 1 and 2, stars growing at rates above 0.01 M yr −1 -0.02 M yr −1 branch off to evolve as red hypergiants while those below these rates become compact blue supergiants. The transition from one regime to the other occurs at ∼ 0.02 M yr −1 , as this SMS oscillates between red and blue phases. As in Hosokawa et al. (2013) and Haemmerlé et al. (2018a), effective temperatures above this rate remain at ∼ 10 4 K even as they reach luminosities of 10 9 − 10 10 L as the stars evolve along the Hayashi limit. Below 0.01 M yr −1 the stars evolve onto the zero-age main sequence (ZAMS) and reach temperatures of up to 10 5.5 K and luminosities of 10 8 − 10 10 L . The two tracks bifurcate during initial contraction, prior to any nuclear burning. This transition can also be seen in the interiors of the stars shown in Figure 4. Above 0.02 M yr −1 there are more radiation dominated structures in which less of the interior resides within the convective core. At lower accretion rates the stars are almost fully convective, and most of the star is a burning core. It is clear from Figure 1 that high accretion rate stars grow significantly in radius, and this is due to the fact that they gain mass more quickly than they can contract. Low accretion rate stars contract more quickly than they gain mass so their radii decrease as they settle onto the ZAMS.

3.1
= 0.1 -1 M yr −1 As shown in the top panel of Figure 2, these cool red stars evolve very quickly over a narrow range of effective temperatures. As they grow in mass, they encounter opacity bumps that lead to excursions in the HR diagram until they converge to a direct track along the Hayashi limit. Energy production by the PP, CNO and triple alpha chains causes the stars to grow in radius by a factor of ten during main sequence burning. This expansion is driven by temperature sensitive H − opacity in their atmospheres, which leads to efficient transport of energy from the interior to the outer layers. In Figure 3, the internal structures of these models change with decreasing accretion rate. At lower accretion rates more of the stellar interior is subsumed by the convective core. The staircase-like growth of the central convective zone is due to the finite resolution of our models and the limitations of our 1D prescription for convective instability.
3.2 = 0.02 -0.1 M yr −1 As shown in the middle panel of Figure 2, most of the stars in this accretion range evolve in the same manner as our highest accretion rate stars, with cool red tracks and radial expansion after converging to the Hayashi limit. The transition from red to blue SMSs occurs at the lower end of this range, as we discuss in greater detail in Section 3.6. The star accreting at 0.02 M yr −1 alternates between red and blue tracks because densities and temperatures in its atmosphere at times are like those in high accretion rate models, but the star does not gain sufficient mass over time to sustain them and periodically contracts back down to a bluer state. These transitions also cause MESA to readjust its timesteps upon excursion from one effective temperature to another. The large swings in the structure of the 0.02 M yr −1 star prematurely end its run but adjustments to the pressure at its surface allow it to be evolved to the end of main sequence burning, as discussed in Section 3.6. There we show that the star continues to cycle between red and blue states before finally settling onto a blue track and becoming hot. As shown in Figure 4, a significant fraction of the mass of the 0.02 M yr −1 star resides within its convective core. At lower accretion rates the percentage of the star that is convective rises. Nuclear energy generation in the convective core is mostly due to CNO, triple alpha and nitrogen burning. Outside the convective zone, weak energy generation proceeds by the PP chain but is dwarfed by core burning. It is clear from Figure 4 that the 0.02 M yr −1 star is evolving towards the low accretion rate regime, as this model is comparatively much older once it reaches 10 4 M .

3.3
= 0.001 -0.01 M yr −1 HR tracks for the SMSs in the lowest accretion range are shown in the bottom panel of Figure 2. They evolve as hot blue stars that can be more than a hundred times smaller than red SMSs of comparable mass. At early stages in their evolution these stars have similar radii as they approach the ZAMS. Just before reaching the postmain sequence, they begin to expand as they deplete their central hydrogen supply, with those with the highest rates having the largest radii and luminosities. We show the internal structures of these stars in Figure 5. Most of the mass of these compact blue stars resides in their convective cores. As at high accretion rates, the staircaselike growth of the central convective zone is again due to the finite resolution of our models and the limitations of our 1D prescription for convective instability.
In principle, convection in these stars could periodically replenish their cores with H from higher mass coordinates. If such an event were to occur after core H mass fractions fall below 1% our models would underestimate the true lifetime and final mass of the star because the ingestion of fresh H into the core would have allowed the star to continue to burn after we terminated the run. However, as we show in Figure 6, such ingestion events only occur near the end of the life of the star at 0.001 M yr −1 . It is not possible at present to evolve this star further because of numerical instabilities that arise in the run during the contraction of the core and the expansion of the star as it becomes fully convective at this stage, so the final mass of this star should be taken to be a lower limit.

Collapse
As discussed in greater detail in the next section, SMSs growing at 0.02 -1 M yr −1 become unstable during hydrogen burning because of changes in central temperatures and densities induced by pulsations due to the post-Newtonian instability and collapse before the end of the main sequence. Once collapse proceeds, infall rates in the core rapidly reach 1000 km s −1 as shown in Figure 7. Stars growing at 0.02 -0.1 M yr −1 evolve further along the main sequence before collapse, with the 0.02 M yr −1 SMS reaching final central hydrogen fractions of 0.12. The 0.01 M yr −1 SMS marks the transition to low accretion rate evolution in which the star reaches the central hydrogen fraction cutoff of 0.01 and develops a helium core as it advances to later stages of burning. Stars that grow at 0.001 -0.01 M yr −1 reach the post-main sequence, evolving in a similar manner as massive Pop III stars. They never encounter the post-Newtonian instability and, although we do not model it here, are expected to collapse during He, C or O burning. Accretion rates of 0.01 M yr −1 -0.02 M yr −1 thus mark a key divide in SMS evolution: whether they evolve as cool red hypergiants or compact blue supergiants and if they collapse via the GRI or because of core fuel depletion.

Hot CNO cycle
If core temperatures in SMSs reach 2 × 10 8 K while still at high central hydrogen fractions, energy generation due rapid proton (rp) captures can rival that of the CNO cycle and become several hundred times greater at temperatures above 5 × 10 8 K because CNO reaction rates are limited by the half-life of one of its beta decays ( -limited CNO; Fuller et al. 1986). However, as shown in Figure 8 core temperatures calculated with the 21-isotope APPROX network never exceed 2 × 10 8 K in our stars except for brief episodes late in   their lives when central hydrogen fractions are low. In principle, the inclusion of more isotopes and rp captures (the hot CNO, or hCNO, cycle) could have produced larger core temperatures than those in our models. We ran the 0.001 M yr −1 and 1.0 M yr −1 stars with the 44-isotope HBURN network with one hCNO cycle and the reduced HCNO network with just the p-p chain, triple alpha chain, CNO and hCNO cycles to compare core temperatures and energy production rates, which are plotted in Figure 9 (we also ran the 1.0 M yr −1 star with the 8-isotope BASIC network, which has the p-p chain, triple alpha chain, and CNO cycle, for comparison). Unlike non-accreting SMS models, which exhibit core temperatures at which the hCNO cycle becomes important (e.g., Fuller et al. 1986;Nagele et al. 2022), we find that the inclusion of the hCNO cycle in our models does not produce central temperatures at which it becomes important. Indeed, there is little difference in core temperatures with the three networks. The inclusion of more isotopes stabilizes energy production in the core at earlier times but results in the same evolution and final mass for the star.

0.02 M / yr SMS Atmosphere Effects
As shown in Figure 2, the 0.02 M yr −1 star alternates between red and blue tracks, suggesting that this is the accretion rate at which SMSs transition from blue to red. Its evolution is sensitive to the choice of atmosphere and surface pressure. We evolved the star with two versions of the standard MESA T( ) atmosphere, one in which the opacity is constant throughout the atmosphere ('fixed', which was used for all the models in our study) and one in which it varies in a manner that is consistent with the local pressure and temperature ('varying'). For the fixed case we tested two extra pressures at the surface of the star that are set by extra_factor = 1.2 and 1.3, as discussed at the end of Section 2. We also consider an alternate formulation for the energy equation, De/Dt, that better conserves energy (Paxton et al. 2019). Figure 10 shows that changing only the pressure boundary for this star can lead to either a hot blue or cool red SMS after the onset of main sequence burning. The varying atmosphere option stabilises the thermal oscillations and leads to a track that is intermediate to the cool red and hot blue regimes. Testing these options in our other models produced much smaller deviations in evolution in the 0.01 M yr −1 and 0.03 M yr −1 stars and few if any in the others, confirming that 0.02 M yr −1 is the transitional accretion rate for rapidly accreting SMSs.
During the protostellar phase the 0.02 M yr −1 SMS initially contracts on a shorter timescale than it can accrete and evolves as a hot blue star. After central nuclear burning begins, the choice of atmosphere leads to deviations in evolution by sending the outer layers of the star into regions of -space that favor or suppress H − formation, whose opacity intercepts energy from the core of the star and can cause the outer layers to expand along the Hayashi limit. In the original model, the surface layers of the star were at temperatures and densities that were marginally favorable to H − formation. As the star expanded its outer layers migrated into regions ofspace in which H − tended to be destroyed (likely because falling densities decreased the two-body H − formation rate) and the star again began to contract to a hotter, bluer phase. Increasing the pressure on the outer boundary moves these layers into regions in -space that are less favorable to H − formation and cause the star to settle onto a blue track at earlier times. Excursions between red and blue states are dampened in the varying case and the star settles onto an intermediate track at early times because the temperature structure of the atmosphere can respond quickly to changes in opacity due to those in density. The run with De/Dt settles onto the hot blue track at early times because the energy equation produces consistently higher surface temperatures that prevent H − from forming and expanding the star. The varying case suggests that the oscillations of this star in the HR diagram are probably not in reality as large as the other test runs suggest but, as noted earlier, this and the De/Dt option had little effect of the evolution of the other stars.

DISCUSSION
High accretion rate stars encounter pulsations due to the GRI at masses above a few tens of thousands of solar masses. These pulsations are initially mild, and can drive shocks into the core with speeds of a few tens of km s −1 . The shocks only induce small changes in central densities and temperatures from which the star can easily recover, unless code time steps that are shorter than nuclear burning times but too long for the hydrodynamics to return accurate velocities lead to unphysically large infall speeds, as discussed earlier. However, as the star grows in mass the pulsations become more violent and drive shocks into the core with velocities of hundreds of km s −1 . These shocks elevate central temperatures and densities that in turn exacerbate the GRI and lead to even stronger pulsations.
The stars reach their final masses when a pulsation finally triggers collapse.
In tests with no hydrodynamics, there are no shocks to induce changes in central densities or temperatures so the model ends when the GRI causes enough softening of the EOS in the core, usually at significantly larger masses as discussed below. We note that both high and low accretion rate stars could be prone to other types of pulsations that did not appear in our models because they had characteristic timescales that were shorter than the time steps taken by our simulations. Further study is required to determine if they occur and what impact they have on the evolution of the star, such as if they can lead to collapse at lower masses than those due to the GRI.
We show final masses for the stars in Figure 11 with those from previous studies at the same accretion rates. As in earlier work, our final masses rise monotonically with accretion rate, and above 0.01 M yr −1 they fall in between those of Haemmerlé et al. (2018a) and Woods et al. (2017). We modeled stars growing at above 0.01 M yr −1 with and without hydrodynamics to test its effects on their evolution and final masses. From 0.01 M yr −1 -0.1 Ms yr −1 they converge to essentially the same mass. Stars growing at ≥ 0.1 M yr −1 without hydrodynamics have consistently larger final masses that are closer to those of Umeda et al. (2016) because they do not manifest the pulsations that would have collapsed the star at earlier times and because the code can take larger time steps that can mitigate other effects of the GRI.
As noted earlier, at accretion rates of 0.01 M yr −1 stars begin to collapse because of hydrogen depletion while at rates of 0.02 M yr −1 they collapse via the GRI. What are the relative roles of the two processes in collapse over this range? Either one can lead to excessive changes in central densities and temperatures that drive catastrophic reductions in time step that signify the death of the star. When we ran the 0.02 M yr −1 star without post-Newtonian corrections it evolved further into hydrogen depletion and encountered numerical difficulties as it entered the post main sequence, but it did not develop large infall velocitiies we would associate with collapse. Likewise, when we ran the 0.01 M yr −1 star without post-Newtonian corrections, the star struggled through hydrogen depletion and then collapsed. These two tests indicate that the GRI was primarily to blame for the collapse of the first star and depletion of core hydrogen caused the collapse of the second. Both likely contribute to collapse at 0.01 M yr −1 -0.02 M yr −1 .
Stars accreting at the lowest decade in rate have high temperatures and large ionizing UV fluxes that could in principle disperse the flows that create them, raising the question of whether such stars could grow at these rates in the first place. Latif et al. (2021) studied the growth of such stars and included the effects of their radiation on their inflows. They found that they could continue to grow at rates in the lowest decade we considered. In reality, no SMS grows at uniform accretion rates, and even cool, red stars whose radiation has no effect on accretion are subject to time-dependent cosmological flows. However, studies of SMSs have begun to consider their evolution in such flows and find that they evolve in much the same way as at constant accretion rates (Woods et al. 2021b,a).
As mentioned in the Introduction, SMS are not expected to lose mass to line-driven winds because of the low opacity of primordial gas. However, the stars soon become highly convective and heavier elements produced by the core could be dredged up and drive winds later in the life of the star. But it is unlikely that these winds would lead to significant mass loss because they would be modest and probably be overcome by the ram pressure of the infalling gas. Consequently, the final mass of the star will just be the mass it accrued over its lifetime. SMSs in a narrow range of masses around 55,000 M in which accretion has subsided and the star has thermally relaxed (e.g., Woods et al. 2020) have been found to explode in previous studies (Montero et al. 2012;Whalen et al. 2013a;Johnson et al. 2013b;Whalen et al. 2013bWhalen et al. , 2014Chen et al. 2014;Nagele et al. 2021;Moriya et al. 2021;Nagele et al. 2022). However, there could be other observational signatures of SMS collapse even if they do not explode. For example, if the star rotates and core collapse drives a jet that pierces its outer layers, it could imprint distinctive nucleosynthetic patterns on surrounding gas that could later be found in the atmospheres of low mass stars forming in it (e.g., Joggerst et al. 2010). Collapse could also produce a strong neutrino signal that is far brighter than that unleashed by conventional core-collapse SNe, but it would still only be detectable at Mpc distances (Shi et al. 1998;Muñoz et al. 2021;Nagele et al. 2021).
Although we do not follow the collapse of our stars to late times, Nagele et al. (2021) found that event horizon formation during the collapse of a 10 4 M SMS intially encloses 40 -50 M of the core of the star, evoking the possibility of the birth of a quasi-star (e.g., Begelman et al. 2008;Volonteri & Begelman 2010). In this picture, X-rays from the BH support the rest of the star from prompt collapse and form a stable envelope that would appear to be a cool, red giant star to an external observer. Such stars can grow to ∼ 10 6 M before the BH becomes so massive that a hydrostatic envelope is no longer possible. However, even if the BH initially had similar masses in our stars it is unlikely that X-rays could halt the collapse of the star because so much of its mass is falling inward at velocities above several hundred km s −1 , as shown in Figure 7. Consequently, DCBHs are born with the mass at which their progenitors die.

CONCLUSION
We have carried out a systematic study of the evolution of rapidly accreting SMS with the publicly-available stellar evolution code MESA. Our grid of models spans 3 decades in accretion rate that bracket the range expected for primordial environments conducive to SMS formation. We find that SMSs evolve along two different pathways, as cool red hypergiants or compact blue supergiants, at accretion rates above and below 0.01 ≤ ≤ 0.02 M yr −1 , respectively. This range also marks the transition between stars collapsing because of the depletion of core fuel after the end of the main sequence at low rates and collapse via the GRI during the main sequence at higher rates.
We also find that hydrodynamics is crucial to capturing how the GRI causes the death of the star at higher accretion rates: triggering pulsations that eventually lead to its collapse. Without hydrodynamics, the GRI still leads to the collapse of the star but at significantly higher masses by softening the EOS in the core and triggering ingestion events that raise central densities and temperatures and destabilise the core. High accretion rate models with hydrodynamics encounter fatal unstabilities at lower final masses before ingestion events occur. When our SMS models reach collapse, large infall velocities enclose most of the mass of the star so it will go into the BH soon after birth, preventing the formation of a quasistar that could create DCBHs of up to 10 6 M (e.g., Begelman et al. 2008). Our results are broadly consistent with previous, more sparsely-sampled accretion rates (Woods et al. 2017;Haemmerlé et al. 2018a) and, critically, provide a framework for future SMS modelling efforts with open-source tools.