Modelling supermassive primordial stars with mesa

ABSTRACT

Supermassive stars forming at z ∼ 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 ultraviolet backgrounds or in highly supersonic baryon streaming flows. Atomic cooling triggers catastrophic baryon collapse capable of building up stars at rates of up to ∼1 M_⊙ yr⁻¹. Here, we examine the evolution of supermassive stars with a much larger and finer grid of accretion rates than in previous studies with the mesa stellar evolution code. We find that their final masses range from 3.5 × 10³ to 3.7 × 10⁵ M_⊙ at accretion rates of 0.001–1 M_⊙ yr⁻¹, respectively. We also find that supermassive star evolution diverges at accretion rates of 0.01–0.02 M_⊙ yr⁻¹, 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.

1 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 z > 6 (Fan et al. 2003, 2006), including eight at z > 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 z ∼ 10–20 in primordial haloes immersed in high Lyman–Werner ultraviolet (UV) backgrounds (Latif et al. 2014b; Schauer et al. 2015; Agarwal et al. 2016; Schauer et al. 2017a) or in highly supersonic baryon streaming motions (BSMs; Latif, Niemeyer & Schleicher 2014a; Hirano et al. 2017; Schauer et al. 2017b) that suppress normal star formation until they reach masses of 10⁷ M_⊙ and virial temperatures of 10⁴ K. These temperatures trigger atomic cooling that causes gas to collapse to the centre of the halo at rates of up to ∼1 M_⊙ yr⁻¹ (Bromm & Loeb 2003; Lodato & Natarajan 2006; Regan & Haehnelt 2009; Latif & Volonteri 2015). Such flows create massive, hot accretion discs that build up either a single SMS or binaries and small multiples (Wise et al. 2019; Latif, Khochfar & Whalen 2020; Patrick et al. 2020; Regan et al. 2020; Patrick et al., in preparation). However, it has now been found that the rare haloes capable of forming quasars by z ∼ 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 1964, 1966) and numerical simulations since the 1970s (e.g. Appenzeller & Fricke 1972; Shapiro & Teukolsky 1979; Fuller, Woosley & Weaver 1986; Baumgarte & Shapiro 1999; Butler et al. 2018; Sun, Ruiz & Shapiro 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, 2021a, b; Haemmerlé et al. 2018a). These models indicate that rapidly accreting SMSs can reach masses of a few 10⁵ M_⊙ and lifetimes of 1–2 Myr before collapsing to direct-collapse 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 haloes (Kitayama et al. 2004; Whalen, Abel & Norman 2004; Alvarez, Wise & Abel 2009; Whalen & Fryer 2012; Smith et al. 2018). DCBHs are born with much larger masses and in much higher densities in haloes that retain their fuel supply, even when it is heated by X-rays (Johnson et al. 2013a).

Previous studies have mostly found that rapidly accreting primordial stars evolve as cool, red hypergiants along the Hayashi limit, with surface temperatures of 5000–10 000 K due to H⁻ opacity in their atmospheres, at least until they reach ∼10⁵ M_⊙ (Hosokawa et al. 2013). Haemmerlé et al. (2018a) found that SMSs can remain cool even above these masses, reaching luminosities greater than 10¹⁰ L_⊙. In other studies, SMSs evolving from similar initial conditions soon settle on to 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⁻¹) also evolve along blue tracks, as can stars with clumpy accretion due to fragmentation or turbulence in the disc (Sakurai et al. 2015, but see Sakurai et al. 2016). It is not clear whether 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 SMS 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 per cent 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, Turk & Abel 2008). Another is the rise of magnetic fields in the accretion disc 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 on to the star without forming the usual hot, high-pressure bubbles that drive all the baryons out of less massive Pop III star-forming haloes (Susa, Umemura & Hasegawa 2009; Susa 2013; Hirano et al. 2014, 2015; Sugimura et al. 2020; Latif et al. 2021; Latif, Whalen & Khochfar 2022b).

Pop III SMSs accreting at rates of 0.1–1 M_⊙ yr⁻¹ would have extremely large luminosities that could, in spite of their low surface temperatures, be detected in the near-infrared (NIR) by the JWST and large ground-based telescopes in the coming decade (Johnson et al. 2012; Hosokawa et al. 2013; Surace et al. 2018, 2019; Whalen et al. 2020a). Their BHs could also be found in the NIR at z ∼ 15–20 by JWST (Pacucci et al. 2015; Natarajan et al. 2017; Barrow, Aykutalp & Wise 2018; Whalen et al. 2020b) and at z ∼ 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 z ∼ 15 (Vikaeus, Whalen & Zackrisson 2022). DCBHs will only be marginally visible to the Square Kilometre Array or next-generation Very Large Array at z ≳ 6–8 (Whalen et al. 2020a, 2021) but would become more luminous after growing to larger masses at later times (Whalen & Mezcua, in preparation). They could also be detected out to z ∼ 10 by future X-ray missions such as the Advanced Telescope for High-Energy Astrophysics 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–10 M_⊙ yr⁻¹. 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 set-ups 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.

2 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. 2011, 2013, 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 data sets such as the OPAL/SCVH tables (Saumon, Chabrier & van Horn 1995; Rogers & Nayfonov 2002), 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 fig. 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 correction to the equation of hydrostatic equilibrium:

which is implemented as a post-Newtonian correction to the gravitational constant, G,

G_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 G at the beginning of each run.

mesa adaptively rezones the stars as they evolve. Mesh refinement is triggered when gradients in temperature, pressure, and ⁴He abundances between adjacent zones exceed preset values: δlog P/P > 1/30, δlog T/T > 1/80, and δlog (χ + 0.01) > 1/20 log χ, where χ is the ⁴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 centre of the star to resolve nuclear burning and convective processes, and the stars are typically partitioned into about 1400 mass zones. Accretion flows on to 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 (Baraffe, Heger & Woosley 2001; Vink, de Koter & Lamers 2001). We initialize all the stars as 20 M_⊙ fully convective, n = 1.5 polytropes with temperatures below 10⁶ 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 three decades: 10⁻³, 10⁻², and 10⁻¹ M_⊙ yr⁻¹.

Models with accretion rates at or above 0.02 M_⊙ yr⁻¹ 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 general relativistic instability (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

where μ is the mean molecular weight, and s/k_B is the entropy per baryon (Woods et al. 2017). We run both cases for each star to determine whether 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⁻¹ 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⁻¹ 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⁷ 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  per cent (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 per cent 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 T_c > 8.15, by which point the central hydrogen fraction in most of the stars has fallen to below 30  per cent.

Pressure must be imposed on the outer boundaries of these stars to prevent superluminal velocities from arising in their low-density outermost zones and ensure the stability of the star (Woods et al. 2017; Woods, Heger & Haemmerlé 2020). This pressure, P, is parametrized by |$P_{\mathrm{extra\_factor}}$|⁠,

where

and g = GM/R², the surface gravity of the star, κ and τ are the opacity and optical depth, and M and L are the mass interior to those layers and their luminosity, respectively. We set |$P_{\mathrm{extra\_factor}} = 1.36$| for our low accretion rate models, 0.01–0.001 M_⊙ yr⁻¹.

3 SMS EVOLUTION

We show Hertzsprung–Russell (HR) plots of the evolution of our stars at accretion rates of 0.001–1 M_⊙ yr⁻¹ in Figs 1 and 2 and Kippenhahn diagrams of their structures over time in Figs 3–5. How the stars evolve primarily depends on how accretion time-scales, |$t_{\textrm {acc}} = M/\dot{M}$|⁠, compare to Kelvin–Helmholz (KH) contraction times, t_KH ≈ GM²/RL, early in their lives, where |$\dot{M}$| and L 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 Fig. 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 reconfiguration that, as discussed earlier, follows either a cool red track or a hot blue one. As shown in Figs 1 and 2, stars growing at rates above 0.01–0.02 M_⊙ yr⁻¹ 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⁻¹, 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⁴ K even as they reach luminosities of 10⁹–10¹⁰ L_⊙ as the stars evolve along the Hayashi limit. Below 0.01 M_⊙ yr⁻¹, the stars evolve on to the zero-age main sequence (ZAMS) and reach temperatures of up to 10^5.5 K and luminosities of 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 Fig. 4. Above 0.02 M_⊙ yr⁻¹, 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 Fig. 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 on to the ZAMS.

3.1 |$\dot{M} =$| 0.1–1 M_⊙ yr⁻¹

As shown in the top panel of Fig. 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 10 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 Fig. 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 |$\dot{M} =$| 0.02–0.1 M_⊙ yr⁻¹

As shown in the middle panel of Fig. 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⁻¹ 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 time-steps upon excursion from one effective temperature to another.

The large swings in the structure of the 0.02 M_⊙ yr⁻¹ 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 on to a blue track and becoming hot. As shown in Fig. 4, a significant fraction of the mass of the 0.02 M_⊙ yr⁻¹ 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 Fig. 4 that the 0.02 M_⊙ yr⁻¹ star is evolving towards the low-accretion rate regime, as this model is comparatively much older once it reaches 10⁴ M_⊙.

3.3 |$\dot{M} =$| 0.001–0.01 M_⊙ yr⁻¹

HR tracks for the SMSs in the lowest accretion range are shown in the bottom panel of Fig. 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 post-main 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 Fig. 5. Most of the mass of these compact blue stars resides in their convective cores. As at high accretion rates, the staircase-like 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 per cent, 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 Fig. 6, such ingestion events only occur near the end of the life of the star at 0.001 M_⊙ yr⁻¹. 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.

3.4 Collapse

As discussed in greater detail in the next section, SMSs growing at 0.02–1 M_⊙ yr⁻¹ 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⁻¹ as shown in Fig. 7. Stars growing at 0.02–0.1 M_⊙ yr⁻¹ evolve further along the main sequence before collapse, with the 0.02 M_⊙ yr⁻¹ SMS reaching final central hydrogen fractions of 0.12. The 0.01 M_⊙ yr⁻¹ SMS marks the transition to low accretion rate evolution in which the star reaches the central hydrogen fraction cut-off 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⁻¹ 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–0.02 M_⊙ yr⁻¹ thus mark a key divide in SMS evolution: whether they evolve as cool red hypergiants or compact blue supergiants and whether they collapse via the GRI or because of core fuel depletion.

3.5 Hot CNO cycle

If core temperatures in SMSs reach 2 × 10⁸ 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⁸ 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 Fig. 8 core temperatures calculated with the 21-isotope APPROX network never exceed 2 × 10⁸ 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 and 1.0 M_⊙ yr⁻¹ 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, and CNO and hCNO cycles to compare core temperatures and energy production rates, which are plotted in Fig. 9 (we also ran the 1.0 M_⊙ yr⁻¹ star with the eight-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.

3.6 0.02 M_⊙yr SMS atmosphere effects

As shown in Fig. 2, the 0.02 M_⊙ yr⁻¹ 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 mesaT(τ) 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 |$P_{\rm {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). Fig. 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 stabilizes 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 and 0.03 M_⊙ yr⁻¹ stars and few if any in the others, confirming that 0.02 M_⊙ yr⁻¹ is the transitional accretion rate for rapidly accreting SMSs.

During the protostellar phase, the 0.02 M_⊙ yr⁻¹ SMS initially contracts on a shorter time-scale 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 T–ρ space that favour 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 favourable to H⁻ formation. As the star expanded, its outer layers migrated into regions of T–ρ space 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 T–ρ space that are less favourable to H⁻ formation and cause the star to settle on to a blue track at earlier times. Excursions between red and blue states are dampened in the varying case and the star settles on to 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 on to 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.

4 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⁻¹. 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⁻¹. 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 time-scales that were shorter than the time-steps taken by our simulations. Further study is required to determine whether they occur and what impact they have on the evolution of the star, such as whether they can lead to collapse at lower masses than those due to the GRI.

We show final masses for the stars in Fig. 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⁻¹ they fall in between those of Haemmerlé et al. (2018a) and Woods et al. (2017). We modelled stars growing at above 0.01 M_⊙ yr⁻¹ with and without hydrodynamics to test its effects on their evolution and final masses. From 0.01 M_⊙ yr⁻¹ to 0.1 M_⊙ yr⁻¹, they converge to essentially the same mass. Stars growing at ≥0.1 M_⊙ yr⁻¹ 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⁻¹ stars begin to collapse because of hydrogen depletion, while at rates of 0.02 M_⊙ yr⁻¹ 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⁻¹ 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 velocities we would associate with collapse. Likewise, when we ran the 0.01 M_⊙ yr⁻¹ 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–0.02 M_⊙ yr⁻¹.

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. 2021a, b).

As mentioned in the Introduction, SMSs 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. However, 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, Janka & Müller 2012; Johnson et al. 2013b; Whalen et al. 2013a, b, 2014; Chen et al. 2014; Moriya et al. 2021; Nagele et al. 2021, 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, Fuller & Halzen 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⁴ M_⊙ SMS initially encloses 40–50 M_⊙ of the core of the star, evoking the possibility of the birth of a quasi-star (e.g. Begelman, Rossi & Armitage 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⁶ 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⁻¹, as shown in Fig. 7. Consequently, DCBHs are born with the mass at which their progenitors die.

5 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 three 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 |$\le \dot{M} \le$| 0.02 M_⊙ yr⁻¹, 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 destabilize 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 quasi-star that could create DCBHs of up to 10⁶ 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.

Acknowledgement

We thank Bill Paxton and the mesa community for valuable discussions that made this work possible. NPH acknowledges funding from the European Research Council for the Horizon 2020 ERC Consolidator Grant project ICYBOB, grant number 818940. DJW was supported by the Ida Pfeiffer Professorship at the Institute of Astrophysics at the University of Vienna and by STFC New Applicant Grant ST/P000509/1. TEW acknowledges support from the NRC-Canada Plaskett Fellowship. Software used in this study are mesa (Paxton et al. 2015, 2019), mesasdk 20.12.1 (Townsend 2019), and py_mesa_reader (Wolf & Schwab 2017).

DATA AVAILABILITY

The data in this study will be made available upon request to the corresponding author.

References