On the Maximum Mass of Accreting Primordial Supermassive Stars

TL;DR

This study models the birth, evolution, and collapse of accreting primordial supermassive stars, predicting their maximum masses and implications for early quasars using advanced stellar evolution simulations.

Contribution

It introduces a comprehensive simulation approach including post-Newtonian corrections to predict the maximum mass of supermassive stars before collapse.

Findings

Stars collapse at 150,000-330,000 solar masses depending on accretion rate.

Final star mass scales logarithmically with accretion rate.

Structural stability is sensitive to convection treatment and outer envelope heat content.

Abstract

Supermassive primordial stars are suspected to be the progenitors of the most massive quasars at z~6. Previous studies of such stars were either unable to resolve hydrodynamical timescales or considered stars in isolation, not in the extreme accretion flows in which they actually form. Therefore, they could not self-consistently predict their final masses at collapse, or those of the resulting supermassive black hole seeds, but rather invoked comparison to simple polytropic models. Here, we systematically examine the birth, evolution and collapse of accreting non-rotating supermassive stars under accretion rates of 0.01-10 solar masses per year, using the stellar evolution code KEPLER. Our approach includes post-Newtonian corrections to the stellar structure and an adaptive nuclear network, and can transition to following the hydrodynamic evolution of supermassive stars after they encounter the general relativistic instability. We find that this instability triggers the collapse of the star at masses of 150,000-330,000 solar masses for accretion rates of 0.1-10 solar masses per year, and that the final mass of the star scales roughly logarithmically with the rate. The structure of the star, and thus its stability against collapse, is sensitive to the treatment of convection, and the heat content of the outer accreted envelope. Comparison with other codes suggests differences here may lead to small deviations in the evolutionary state of the star as a function of time, that worsen with accretion rate. Since the general relativistic instability leads to the immediate death of these stars, our models place an upper limit on the masses of the first quasars at birth.

Equations

$${\Gamma }_1-\frac{4}{3}\lesssim 1.12\frac{R_{\mathrm{S}}}{R},$$

$${\Gamma }_1\approx \frac{4}{3}+\frac{\beta }{6}$$

$$M_{\mathrm{SMS},\mathrm{final}}\approx \left[0.83{\log }_{10}\left(\frac{\dot{M}}{{\mathrm{M}}_{\odot }{\mathrm{yr}}^{-1}}\right)+2.48\right]\times 10^5{\mathrm{M}}_{\odot }$$

Full text

1. Introduction

The possible existence of “supermassive” stars, with $M\gtrsim 10^{4}\,{\mathrm{M}_{\odot}}$, 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\gtrsim 6$ (with $M\sim 10^{9}\,{\mathrm{M}_{\odot}}$, 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^{7}-10^{8}\,{\mathrm{M}_{\odot}}$ 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 $\rm{H}_{2}$ 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 $\approx 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\,{{\mathrm{M}_{\odot}}\,{\mathrm{yr}}^{-1}}$ (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 which could lead to such halos (e.g., Inayoshi & Omukai, 2012; Inayoshi et al., 2015, Yoshida et al., submitted). A candidate direct collapse black hole (DCBH) has now been discovered: CR7, a Ly-$\alpha$ emitter at $z=6.6$ (Sobral et al., 2015; Hartwig et al., 2016; Agarwal et al., 2016).

In the simplest case, this produces a $10^{4}-10^{5}\,{\mathrm{M}_{\odot}}$ 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 $\Gamma_{1}=4/3$ in order for pressure support to stabilize the star against radial pulsations. For an $n=3$ polytrope, this occurs when

[TABLE]

where ${R_{\mathrm{S}}}=2\,GMc^{-2}$ is the Schwarzschild radius of the stars, and $\Gamma_{1}=\left(\partial\ln P/\partial\ln\rho\right)_{\mathrm{ad}}$ is the adiabatic exponent (Fowler, 1964). For radiation-pressure dominated stars ($\beta=P_{\mathrm{gas}}/P_{\mathrm{tot}}\ll 1$), we have the approximate relation:

[TABLE]

where $\beta\approx 4\,k_{\mathrm{B}}/\left(\mu\,s\right)$, $\mu$ is the mean molecular weight, and $s/k_{\mathrm{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 $\sim 10^{5}\,{\mathrm{M}_{\odot}}$ at several constant accretion rates and found that they remain red and cool until they reach a few $10^{4}\,{\mathrm{M}_{\odot}}$. 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 SMS are discussed in Section 3. We consider the implications of our simulations for the properties of high redshift quasars in Section 4.

2. Numerical Method

2.1. Kepler

Kepler (Weaver et al., 1978; Fuller et al., 1986; Woosley et al., 2002) is a one-dimensional (1D) Lagrangian hydrodynamics and stellar evolution code with nuclear burning and mixing due to convection. Kepler can use a 19-isotope “APPROX” nuclear reaction network (Weaver et al., 1978) which is adequate for many stellar evolution situations. However, for the studies presented here, we use an adaptive full nuclear reaction network implicitly coupled to the hydrodynamics (Woosley et al., 2004). Kepler uses an equation of state similar to, and compatible with, the Helmholtz equation of state (EOS; Timmes & Swesty, 2000), which includes contributions from degenerate and non-degenerate relativistic and non-relativistic electrons, electron-positron pair production, and radiation.

The star is partitioned into a maximum of 1982 zones in mass. Consistent with the general consensus that massive Pop III stars do not lose much mass over their lives, we turn off mass loss in our models, assuming that any mass loss will be negligible relative to the rapid accretion rates considered here. In particular, pulsational mass loss is also unlikely to be significant for accreting supermassive stars (Hosokawa et al., 2013). In order to capture the physics of the GR instability, we use the 1st order post-Newtonian approximation to the Tolman-Oppenheimer-Volkoff correction to the structure of the star (Zeldovich & Novikov, 1971; Kippenhahn et al., 2012) as described in Section 2.1 of Chen et al. (2014).

The accreted material is assumed to be of the same primordial composition as the cloud from which the original seed of the star has formed. We assume cold accretion (following, e.g., Hosokawa et al., 2013), which means material is accreted with the surface entropy. In order to accurately follow the transfer of heat through the accreted envelope, we switch from a wholly Lagrangian to a semi-Eulerian treatment in the outermost part of the star. We make the transition at a specified total optical depth measured inward from the surface, here chosen to be $10^{6}$ in order to ensure the outer surface layers remain well-resolved. In practice, this amounts to adding a homologous term to the gravitational energy release in the outer layers of the star, in order to account for the work done by compression of the accreted material. In principle, some fraction of the accretion luminosity should be released at the accretion shock, a part of which should contribute to an additional heating term at the surface of the star. This was found to have a negligible effect on the evolution of supermassive stars, assuming any fraction $0\leq\eta\leq 1$ of this luminosity is added to the bottom of the photosphere (Hosokawa et al., 2013). A more detailed treatment may change this picture slightly, however given the current uncertainties in the formation of supermassive stars, we ignore this term in the present work.

2.2. Initial Conditions

Due to the difficulty in achieving numerically stable results when beginning at lower stellar masses, we initialize all of our models as $10\,{\mathrm{M}_{\odot}}$, chemically homogeneous, $n=3$ polytropes. The initial central density, $\rho_{\mathrm{c}}$, is set to $10^{-3}\,{\mathrm{g}\,\mathrm{cm}^{-3}}$, giving a central temperature $T_{\mathrm{c}}=1.2\times 10^{6}\,{\mathrm{K}}$ capable of sustaining deuterium-burning. We assume a baryon-to-photon ratio, $\eta$, equal to $6\times 10^{-10}$ in the primordial gas (Cyburt et al., 2001, 2002) which is consistent with the results of Big Bang Nucleosynthesis, and is consistent with initial mass fractions of hydrogen, 4He, 3He, 2H, and Li of $\approx 0.75$, $\approx 0.25$, $\approx 2.1\times 10^{-5}$, $\approx 4.3\times 10^{-5}$, and $\approx 1.9\times 10^{-9}$ respectively (B. D. Fields, priv. comm.).

We consider seven constant accretion rates uniformly spaced in $\log\dot{M}$ from $0.01\,{{\mathrm{M}_{\odot}}\,{\mathrm{yr}}^{-1}}$ to $10\,{{\mathrm{M}_{\odot}}\,{\mathrm{yr}}^{-1}}$ that span the infall rates expected in Lyman-cooled halos based on cosmological simulations. The lowest, $0.01\,{{\mathrm{M}_{\odot}}\,{\mathrm{yr}}^{-1}}$, corresponds to the maximum accretion rate due to H2 cooling in less massive halos.

3. The Evolution of Accreting Supermassive Stars

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 Fig. 1, we show the evolution of the stellar structure and energy generation with time in four Kippenhahn diagrams for the $0.01\,{{\mathrm{M}_{\odot}}\,{\mathrm{yr}}^{-1}}$, $0.1\,{{\mathrm{M}_{\odot}}\,{\mathrm{yr}}^{-1}}$, $1\,{{\mathrm{M}_{\odot}}\,{\mathrm{yr}}^{-1}}$, and $10\,{{\mathrm{M}_{\odot}}\,{\mathrm{yr}}^{-1}}$ 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 ($\lesssim 10^{-1.5}\,{{\mathrm{M}_{\odot}}\,{\mathrm{yr}}^{-1}}$), 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.

As a representative example, we consider a star that grows at a rate of $1\,{{\mathrm{M}_{\odot}}\,{\mathrm{yr}}^{-1}}$. It reaches a central temperature of $10^{8}\,{\mathrm{K}}$ after $\sim 1,485$ yr at a density of $20\,{\mathrm{g}\,\mathrm{cm}^{-3}}$. This allows the core to produce a rapid spike in CNO elements through the triple-$\alpha$ process, catalyzing hydrogen-burning until the core is stabilized against further contraction (see Fig. 2). The star develops a small convective core shortly after this ($\approx 2,000\,{\mathrm{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 $\sim 300,000\,{\mathrm{M}_{\odot}}$ (that is, after $\sim 300,000\,{\mathrm{yr}}$). The inner convective core at this time has reached $\approx 60,000\,{\mathrm{M}_{\odot}}$ and a helium fraction of $\sim\,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.

For all of our models which accrete above $\gtrsim 10^{-1.5}\,{{\mathrm{M}_{\odot}}\,{\mathrm{yr}}^{-1}}$, 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\,{{\mathrm{M}_{\odot}}\,{\mathrm{yr}}^{-1}}$ and $10^{-1.5}\,{{\mathrm{M}_{\odot}}\,{\mathrm{yr}}^{-1}}$, we find that our stellar models reach total masses of only $32,000\,{\mathrm{M}_{\odot}}$ and $84,000\,{\mathrm{M}_{\odot}}$, respectively, before exhausting hydrogen in their cores. Our $10^{-1.5}\,{{\mathrm{M}_{\odot}}\,{\mathrm{yr}}^{-1}}$ model collapses immediately upon core hydrogen-exhaustion, whereas our $0.01\,{{\mathrm{M}_{\odot}}\,{\mathrm{yr}}^{-1}}$ 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\,{{\mathrm{M}_{\odot}}\,{\mathrm{yr}}^{-1}}$ model, as seen in their Fig. 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 $\approx 42,000\,{\mathrm{M}_{\odot}}$, 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 $\approx 330,000\,{\mathrm{M}_{\odot}}$ in our models, contrasted with their $800,000\,{\mathrm{M}_{\odot}}$. Indeed, evidently the $10\,{{\mathrm{M}_{\odot}}\,{\mathrm{yr}}^{-1}}$ model of Umeda et al. (2016), is far less evolved than our $10\,{{\mathrm{M}_{\odot}}\,{\mathrm{yr}}^{-1}}$ 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\,{{\mathrm{M}_{\odot}}\,{\mathrm{yr}}^{-1}}$.

In Fig. 3, we visualize the criterion for the onset of the GR instability as derived for an $n=3$ polytrope (Eq. 1, see Chandrasekhar, 1964), plotting $\beta/6$ and $1.12\,{r_{\mathrm{S}}}/r$ against mass coordinate, $m=m\left(r\right)$. Here ${r_{\mathrm{S}}}={r_{\mathrm{S}}}\left(m\right)=2\,Gmc^{-2}$ is the “local Schwarzschild radius” of the enclosed mass $m$. As found by Umeda et al. (2016), long before the onset of collapse $\beta/6>1.12\,{r_{\mathrm{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 Fig. 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 a $n=3$ polytrope (given the accreted envelope), nor is the core alone reasonably approximated as such, although $\beta\ll 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.

Our results are summarized in Fig. 4, where the final masses and evolutionary characteristics of supermassive stars are given as a function of the accretion rate. For accretion rates $\lesssim 10^{-1.5}\,{{\mathrm{M}_{\odot}}\,{\mathrm{yr}}^{-1}}$, stars collapse only after the onset of He-burning, while for greater accretion rates, GR collapse occurs during core hydrogen-burning. For accretion rates $\gtrsim 0.1\,{{\mathrm{M}_{\odot}}\,{\mathrm{yr}}^{-1}}$, the final mass at collapse varies as:

[TABLE]

fitting our results to within $\approx 3\,\%$. Also shown for reference in Fig. 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.

4. Discussion and Conclusion

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 $\lesssim 300,000\,{\mathrm{M}_{\odot}}$ on accreting, primordial supermassive stars, as for accretion rates above $\approx 0.1\,{{\mathrm{M}_{\odot}}\,{\mathrm{yr}}^{-1}}$ 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 which 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, Haemmerle et al., in prep.). 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.

AH 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. LH,DJW, and RSK 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.