On the Maximum Mass of Accreting Primordial Supermassive Stars

The possible existence of "supermassive" stars, with M ≳ 10⁴ M_⊙, has been suggested since the early 1960s (e.g., Iben 1963; Fowler 1964). Only recently, however, have they been suspected to be necessary to explain the formation of at least the most massive quasars found at z ≳ 6 (with M ∼ 10⁹ M_⊙; e.g., Mortlock et al. 2011; Wu et al. 2015). In particular, lower-mass black holes formed from more typical Population III (Pop III) stars could not have sustained the high accretion rates needed to grow to such masses by this time (Whalen et al. 2004; Park & Ricotti 2011; Whalen & Fryer 2012; though hyper-Eddington accretion rates may be possible, see, e.g., Pezzulli et al. 2016).

In the supermassive star scenario, a primordial halo grows to masses of 10⁷–10⁸ M_⊙ without ever having formed a star, most likely because it is exposed to a strong Lyman–Werner UV field from nearby star-forming regions (Agarwal et al. 2012; Dijkstra et al. 2014). This destroys H₂ and prevents the early collapse and fragmentation of the cloud into lower-mass objects. In the absence of cooling by molecular hydrogen lines, the gas reaches temperatures of ≈8000 K, at which point cooling by line emission from collisionally excited atomic hydrogen becomes activated. As the cloud collapses isothermally, the accretion rate scales with the cube of the sound speed, with atomically cooled halos permitting catastrophic infall rates of 0.01–10 M_⊙ yr⁻¹ (Wise et al. 2008; Regan & Haehnelt 2009; Shang et al. 2010; Wolcott-Green et al. 2011; Choi et al. 2013; Latif et al. 2013). Other scenarios have also been proposed that could lead to such halos (e.g., Inayoshi & Omukai 2012; Inayoshi et al. 2015). A candidate direct collapse black hole (DCBH) has now been discovered: CR7, a Lyα emitter at z = 6.6 (Sobral et al. 2015; Agarwal et al. 2016; Hartwig et al. 2016).

In the simplest case, this produces a 10⁴–10⁵ M_⊙ star at the center of the halo, which will collapse directly to a black hole via the general relativistic (GR) instability (Iben 1963; Chandrasekhar 1964). This arises because general relativity requires a slightly greater adiabatic exponent than that of the radiation-pressure-dominated gas with Γ₁ = 4/3 in order for pressure support to stabilize the star against radial pulsations. For an n = 3 polytrope, this occurs when

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{1}-\displaystyle \frac{4}{3}\lesssim 1.12\,\displaystyle \frac{{R}_{{\rm{S}}}}{R},\end{eqnarray} \tag{ 1 }$$

where R_S = 2 GMc⁻² is the Schwarzschild radius of the stars and ${{\rm{\Gamma }}}_{1}={(\partial \mathrm{ln}P/\partial \mathrm{ln}\rho )}_{\mathrm{ad}}$ is the adiabatic exponent (Fowler 1964). For radiation-pressure-dominated stars (β = P_gas/P_tot ≪ 1), we have the approximate relation

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{1}\approx \displaystyle \frac{4}{3}+\displaystyle \frac{\beta }{6},\end{eqnarray} \tag{ 2 }$$

where β ≈ 4 k_B/(μ s), μ is the mean molecular weight, and s/k_B is the entropy per baryon.

Fuller et al. (1986) were the first to simulate the full evolution of a supermassive star through the hydrodynamic collapse due to the GR instability, but considered only the case of monolithically formed stars, that is, beginning their calculations with initially supermassive models. More recently, three studies have examined the more realistic case, following the growth of a supermassive Pop III protostellar core given the accretion rates expected in line-cooled halos. Hosokawa et al. (2013) followed the growth of supermassive Pop III stars up to ∼10⁵ M_⊙ at several constant accretion rates and found that they remain red and cool until they reach a few 10⁴ M_⊙. Sakurai et al. (2015) studied the evolution of such stars in clumpy accretion scenarios and found that the protostar could become intermittently blue and hot at low masses but eventually evolved onto a redder, cooler track. Finally, Umeda et al. (2016) included post-Newtonian corrections in their stellar evolution calculations, considering a sparse sample of accretion rates. All three studies employed stellar structure codes that were not hydrodynamical and therefore could not follow the collapse of such stars by the GR instability. Thus, they could not independently determine the final mass of the star, in order to estimate the supermassive black hole (SMBH) seed mass at birth, but rather compared the structure of their models with the above criterion, which is strictly valid only for polytropes.

We have now modeled the birth, evolution, and collapse of supermassive Pop III stars to DCBHs with the one-dimensional (1D) implicit hydrodynamics and stellar evolution code Kepler, which includes post-Newtonian corrections to the structure of the star allowing it to capture the GR instability. Kepler can transition to resolving hydrodynamic timescales should the stellar model under consideration encounter a dynamical instability. This allows us to model the complete evolution of the star until the end of its life and through the onset of its collapse to a black hole. Our approach includes accelerated nuclear burning driven by collapse that might be capable of slowing or reversing it. We obtain final masses for these stars over the range of central collapse rates in atomically cooled halos found in high-resolution cosmological simulations. Our Kepler models are described in Section 2, and the evolution and final masses of Pop III supermassive stars (SMSs) are discussed in Section 3. We consider the implications of our simulations for the properties of high-redshift quasars in Section 4.

2.1. Kepler

2.2. Initial Conditions

The evolution of Pop III SMSs depends sensitively on the balance between the accretion of new material, the nuclear evolution of the core, and the thermal relaxation of the star. In Figure 1, we show the evolution of the stellar structure and energy generation with time in four Kippenhahn diagrams for the 0.01, 0.1, 1, and 10 M_⊙ yr⁻¹ cases. Each model develops a distinct inner convective core and outer, high-entropy envelope, as expected (e.g., Begelman 2010). There is some uncertainty in the precise outer boundary of the convective core, becoming most pronounced for the lowest accretion rates, where the difference in entropy between core and envelope is lowest. This in turn is a consequence of the greater thermal relaxation of the envelope at lower accretion rates. This ambiguity is an inescapable consequence of the simplified 1D treatment of convection utilized in 1D stellar evolution codes. At relatively low accretion rates (≲10^−1.5 M_⊙ yr⁻¹), the fate of the star is principally set by the hydrogen nuclear-burning lifetime of the star. At higher accretion rates, the star encounters the GR instability while still on the hydrogen-burning main sequence, due to the growth of the convective core.

Zoom In Zoom Out Reset image size

As a representative example, we consider a star that grows at a rate of 1 M_⊙ yr⁻¹. It reaches a central temperature of 10⁸ K after ∼1485 yr at a density of 20 g cm⁻³. This allows the core to produce a rapid spike in CNO elements through the triple-α process, catalyzing hydrogen burning until the core is stabilized against further contraction (see Figure 2). The star develops a small convective core shortly after this (≈2000 yr), with an extended, high-entropy envelope. The core grows approximately linearly with time until encountering the GR instability when the total mass of the star reaches ∼300,000 M_⊙ (that is, after ∼300,000 yr). The inner convective core at this time has reached ≈60,000 M_⊙ and a helium fraction of ∼50%. After the onset of collapse, infall velocities quickly reach a few percent of the speed of light. We halt our calculation here, before the collapse becomes strongly relativistic, as our models include only post-Newtonian corrections to gravity but no relativistic hydrodynamics.

Zoom In Zoom Out Reset image size

For all of our models that accrete above ≳10^−1.5 M_⊙ yr⁻¹, the evolution follows similarly to that described above, except that their longer lifetimes prior to collapse via the GR instability allow them to reach core helium fractions of up to 99%. For accretion rates of 0.01 M_⊙ yr⁻¹ and 10^−1.5 M_⊙ yr⁻¹, we find that our stellar models reach total masses of only 32,000 M_⊙ and 84,000 M_⊙, respectively, before exhausting hydrogen in their cores. Our 10^−1.5 M_⊙ yr⁻¹ model collapses immediately upon core hydrogen-exhaustion, whereas our 0.01 M_⊙ yr⁻¹ model survives to silicon burning before the collapse of the core.

Turning to the highest accretion rates, the most strikingly different feature is the emergence of additional convective zones in the envelope, driven by instabilities arising near the surface and accompanied by radial pulsations. These convective zones are also seen in Umeda et al. (2016), as is evident from the small flat plateaus in the entropy profiles in the envelope of their 10 M_⊙ yr⁻¹ model, as seen in their Figure 3 (top left). Therefore, this cannot be a source of any discrepancy between our results and theirs regarding the evolution and fate of the star. Our model encounters the GR instability at a convective core mass ≈42,000 M_⊙, similar to the point at which they cease their calculations after entering the pair-unstable regime. This, however, occurs at a total mass of only ≈330,000 M_⊙ in our models, contrasted with their 800,000 M_⊙. Indeed, evidently the 10 M_⊙ yr⁻¹ model of Umeda et al. (2016) is far less evolved than our 10 M_⊙ yr⁻¹ model at any fixed total mass. The primary reason for this is unclear, but must arise through some difference in the entropy profile and convective stability of the star. This emphasizes the uncertainty in the treatment of convection and the cooling of the outer envelope. Haemmerlé et al. (2017) use a third stellar evolution code and find final masses similar to ours at low accretion rates (as do Umeda et al. 2016), but intermediate between those presented here and that of Umeda et al. (2016) for 10 M_⊙ yr⁻¹.

In Figure 3, we visualize the criterion for the onset of the GR instability as derived for an n = 3 polytrope (Equation (1); see Chandrasekhar 1964), plotting β/6 and 1.12 r_S/r against mass coordinate, m = m(r). Here, r_S = r_S(m) = 2 Gmc⁻² is the "local Schwarzschild radius" of the enclosed mass m. As found by Umeda et al. (2016), long before the onset of collapse β/6 > 1.12 r_S/r throughout the convective core and, indeed, almost the entirety of the envelope, whereas near the point of collapse the region where this condition is violated has passed into the isentropic core. The moment when this is satisfied at the surface of the core, however, does not mark the onset of collapse in our models, nor did it mark the endpoint in the calculations of Umeda et al. (2016), as in the final panel of their Figure 3 the intersection of the two curves is already well within the core boundary. This is reasonable, as the entirety of the star is clearly not well described by an n = 3 polytrope (given the accreted envelope), nor is the core alone reasonably approximated as such, although β ≪ 1. We then conclude that, even when applied only to the isentropic core, the Chandrasekhar (1964) criterion does not adequately describe the mass needed to trigger collapse in realistic supermassive stars. A modified polytropic approximation (Begelman 2010) may be able to better capture the critical pressure gradient needed, but this should also depend on the accretion and evolutionary history of the star, rendering it less useful. Following the hydrodynamic response of the star is necessary then in order to be certain of the moment the GR instability sets in.

Zoom In Zoom Out Reset image size

Our results are summarized in Figure 4, where the final masses and evolutionary characteristics of supermassive stars are given as a function of the accretion rate. For accretion rates ≲10^−1.5 M_⊙ yr⁻¹, stars collapse only after the onset of He burning, while for greater accretion rates, GR collapse occurs during core hydrogen burning. For accretion rates ≳0.1 M_⊙ yr⁻¹, the final mass at collapse varies as

$$\begin{eqnarray}&&{M}_{\mathrm{SMS},\mathrm{final}}\approx \left[0.83\,{\mathrm{log}}_{10}\left(\displaystyle \frac{\dot{M}}{{M}_{\odot }\,{\mathrm{yr}}^{-1}}\right)+2.48\right]\times {10}^{5}\,{M}_{\odot },\end{eqnarray} \tag{ 3 }$$

fitting our results to within ≈3%. Also shown for reference in Figure 4 are the maximum masses for hydrostatic hydrogen and helium burning for monolithically formed (as opposed to accreting) supermassive stars, found from trial computations beginning with supermassive, n = 3 polytropes.

Zoom In Zoom Out Reset image size

We find that, in general, supermassive stars survive well past the point where the isentropic core satisfies the Chandrasekhar (1964) criterion, before finally encountering the hydrodynamic GR instability. This is a natural consequence of the much different structure of an accreting supermassive star compared with a simple n = 3 polytrope. Our models predict that the general relativistic instability imposes a characteristic final total mass of ≲300,000 M_⊙ on accreting, primordial supermassive stars, as for accretion rates above ≈0.1 M_⊙ yr⁻¹ the final mass at collapse grows only logarithmically with the accretion rate. Such infall rates are thought to be typical of atomically cooled primordial composition clouds, providing the seeds of SMBHs formed via direct collapse including high-redshift, massive quasars. Since BHs built up in runaway collisions are generally an order of magnitude lower in mass (Devecchi & Volonteri 2009), these results are likely an upper limit on the masses of SMBHs that may have formed via direct collapse.

This picture could still change if stellar rotation is taken into account (Shibata et al. 2016). Rotation rates of up to 50% of the critical, or breakup, velocity at birth have been found for much lower-mass Pop III stars in numerical simulations (Stacy et al. 2011, 2013; L. Haemmerle et al. 2017, in preparation). Rotational mixing could alter the structure of the star and support it against collapse up to somewhat higher masses. This potentially has observational implications as the collapse of such stars could be accompanied by strong gravitational-wave emission (e.g., Fryer et al. 2001). Rotation may or may not also lead to significant mass loss driven by outflows during the final collapse of the star, depending on the internal redistribution of angular momentum during collapse (Fiacconi & Rossi 2017; Uchida et al. 2017).

The collapse of massive, atomically cooled primordial halos allows for extreme accretion rates sustained for a remarkable duration, far surpassing any star formation environment encountered later in the evolution of the universe. The resulting supermassive objects must then have been among the most massive stars to have ever existed and are strong candidates for the progenitors of some of the first and most luminous quasars.

A.H. was supported by an Australian Research Council (ARC) Future Fellowship (FT120100363) and NSF grant PHY-1430152 (JINA-CEE). D.J.W. was supported by STFC New Applicant Grant ST/P000509/1. L.H., D.J.W., and R.S.K. were supported by the European Research Council under the European Community's Seventh Framework Programme (FP7/2007-2013) via the ERC Advanced Grant "STARLIGHT: Formation of the First Stars" (project number 339177). Part of this work was supported by the Swiss National Science Foundation. Work at LANL was done under the auspices of the National Nuclear Security Administration of the U.S. Department of Energy at Los Alamos National Laboratory under Contract No. DE-AC52-06NA25396.