THE GENERAL RELATIVISTIC INSTABILITY SUPERNOVA OF A SUPERMASSIVE POPULATION III STAR

The existence of supermassive black holes (SMBHs) in most massive galaxies today poses significant challenges to the paradigm of hierarchical structure formation (Kormendy & Richstone 1995; Ferrarese & Merritt 2000; Ferrarese & Ford 2005; Gebhardt et al. 2000; Di Matteo et al. 2005; Beifiori et al. 2012; McConnell & Ma 2013). In particular, it is not known how some SMBHs reached masses of 10⁹ M_☉ by z ∼ 7, less than a Gyr after the big bang (Fan et al. 2002, 2006; Mortlock et al. 2011). The two leading contenders for the origins of SMBHs are 100–500 M_☉ Population III (Pop III) stars at z ∼ 20 and supermassive (10⁴–10⁵ M_☉) Pop III stars at z ∼ 10–15 (e.g., Rees 1984; Madau & Rees 2001; Volonteri 2010, 2012; Karlsson et al. 2013; Inayoshi et al. 2013).

The main argument against conventional Pop III stars as SMBH seeds is that their black holes (BHs) must accrete at the Eddington limit continuously down to z ∼ 7 to become supermassive. How they could sustain such growth is unclear because they form in low-density relic H ii regions that delay initial accretion (Whalen et al. 2004; Johnson & Bromm 2007). Later, when they do accrete, the relatively shallow dark matter potential wells in which they reside cannot retain their fuel supply because of radiative feedback (Alvarez et al. 2009), so they grow at well below the Eddington limit (Park & Ricotti 2011, 2012, 2013). Furthermore, many low-mass Pop III BHs are ejected from their halos (and thus their fuel supply) at birth by asymmetries in their explosion engines (Whalen & Fryer 2012).

For these reasons, there has been growing interest in supermassive stars (SMSs) as candidates for SMBH seeds (e.g., Loeb & Rasio 1994; Bromm & Loeb 2003; Begelman et al. 2006; Bromm & Yoshida 2011; Johnson et al. 2013c). Such stars could form in atomically cooled halos at z ∼ 15. In this scenario, a primeval galaxy forms in a strong Lyman–Werner ultraviolet (UV) background that sterilizes its constituent halos of H₂ prior to assembly (Agarwal et al. 2012; Johnson et al. 2013a), preventing them from forming stars. When the protogalaxy reaches a mass of ∼10⁸ M_☉ and virial temperature of ∼10⁴ K, atomic cooling triggers the rapid collapse of gas at its center at rates of up to 1 M_☉ yr⁻¹ (Lodato & Natarajan 2006; Spaans & Silk 2006; Wise et al. 2008; Regan & Haehnelt 2009; Shang et al. 2010; Latif et al. 2013a, 2013b; Schleicher et al. 2013). In this manner, a 10⁴–10⁵ M_☉ Pop III star could be formed in less than a Myr.

The evolution of supermassive protostars and stars has been studied with both analytical and numerical techniques for more than 40 yr (Fowler 1966; Wheeler 1977; Bond et al. 1984; Carr et al. 1984; Fuller et al. 1986; Bisnovatyi-Kogan 1998; Baumgarte & Shapiro 1999; Begelman 2010; Hosokawa et al. 2013; Schleicher et al. 2013). Most SMSs are thought to become fully convective and then later collapse to massive BHs (e.g., Shapiro & Teukolsky 1979; Shibata & Shapiro 2002; Reisswig et al. 2013). However, previous calculations of SMS evolution were done at low resolution; considered a relatively coarse grid of models in mass; and in some cases, ignored first-order, post-Newtonian general relativistic corrections to gravity. The latter could play an important role in the evolution of 10³–10⁵ M_☉ stars, especially during their pre-explosion phase.

We have revisited the evolution of Pop III SMSs with high-resolution one-dimensional (1D) simulations with updated nuclear reaction rates (Heger et al. 2001; Heger & Woosley 2002) and post-Newtonian corrections to gravity. We have found a 55,500 M_☉ star that explodes as a highly energetic thermonuclear supernova (SN) instead of collapsing to a BH (see also Montero et al. 2012). Our discovery suggests that there may be a range of high-mass stars that can die as SNe. In Section 2, we describe our stellar evolution model. The explosion of the star is examined in one dimension (1D) in Section 3 and in two dimensions (2D) in Section 4. We discuss the potential impact of this explosion on its host galaxy and its visibility at high redshift in Section 5.

The SMS is evolved from the beginning of the main sequence to the onset of collapse and then explosion in 1D in KEPLER. To verify these processes in a separate code, and to capture the violent fluid instabilities that can occur during the SN, we also model the collapse and explosion of the SMS in 2D in CASTRO with initial conditions taken from KEPLER. Combining 1D and 2D simulations in this manner allows us to realistically simulate all phases of the SMS within a practical computational budget.

2.1. KEPLER

KEPLER (Weaver et al. 1978) is a 1D Lagrangian hydrodynamics and stellar evolution code. It includes nuclear burning and mixing due to convection. For a given stellar mass, the initial profiles of density and temperature of a star are determined by solving the Lane–Emden equation (Chandrasekhar 1939). We use the 19 isotope APPROX nuclear reaction network (Weaver et al. 1978; Timmes 1999), which includes heavy ion reactions, alpha chain reactions, hydrogen-burning cycles, photo-disintegration of heavy nuclei, and energy loss through thermal neutrinos. Nuclear burning is self-consistently coupled to hydrodynamics, and we also account for energy deposition due to radioactive decay of ⁵⁶Ni → ⁵⁶Co → ⁵⁶Fe. KEPLER uses 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 1201 zones in mass for an effective resolution of 46.2 M_☉. Consistent with the usual convention that massive Pop III stars lose little of their mass over their lives (Kudritzki 2000; Vink et al. 2001; Baraffe et al. 2001; Krtička & Kubát 2006; Ekström et al. 2008), we turn off mass loss in our models.

When stars become extremely massive, ≳1000 M_☉, general relativity (GR) must be taken into account in gravity in stellar evolution models. Instead of full GR, we apply first-order, post-Newtonian GR corrections to gravity in our calculations. We use the Tolman–Oppenheimer–Volkoff equation for hydrostatic equilibrium in GR to calculate the pressure P (Zeldovich & Novikov 1971; Kippenhahn & Weigert 1990):

$$\begin{eqnarray} &&\frac{dP}{dr} = -\frac{Gm}{r^{2}}\varrho \left(1+\frac{P}{\varrho c^{2}} \right) \left(1+\frac{4\pi r^3 P}{m c ^{2}} \right) \left(1-\frac{2Gm}{r c^{2}} \right)^{-1}, \nonumber\\ \end{eqnarray} \tag{ 1 }$$

where r is the radius of star, c is the speed of light, G is the gravitational constant, and m is the enclosed relativistic mass at r, which is the rest mass plus total energy (internal plus gravitational) divided by c². The relativistic density ϱ is therefore ϱ₀ + U/c², where ϱ₀ and U are the rest mass and total energy densities. When c² → ∞, Equation (1) reduces to Newton's Law:

$$\begin{equation} \frac{dP}{dr} = -\frac{Gm}{r^{2}}\varrho. \quad \end{equation} \tag{ 2 }$$

The effective gravity $\tilde{g}$ is calculated from

$$\begin{equation} \tilde{g} = -\frac{Gm}{r^{2}} \left(1+\frac{P}{\varrho c^{2}} \right) \left(1+\frac{4\pi r^3 P}{m c^{2}} \right) \left(1 - \frac{2Gm}{r c^{2}} \right)^{-1}. \quad \end{equation} \tag{ 3 }$$

The U/c² term in ϱ is negligible in our models. It only becomes comparable to the rest energy when T becomes ∼10¹³ K (kT ∼ m_pc², where m_p is the proton mass). Equation (3) therefore becomes

$$\begin{equation} \tilde{g} = g \left(1+\frac{P}{\varrho _0 c^{2}}+\frac{4\pi r^3 P}{m c^{2}}+\frac{2Gm}{r c^{2}} \right), \quad \end{equation} \tag{ 4 }$$

where g is the Newtonian gravity field, which is determined from the gravity equation in KEPLER and from the monopole approximation to gravity for nearly spherical matter distributions in CASTRO.

2.2. CASTRO

The SMS has a radius of 1.73 × 10¹³ cm, an effective surface temperature of about 7.19 × 10⁴ K and a luminosity of 5.70 × 10⁴² erg s⁻¹ during its pre-collapse phase. The star evolves about 1.69 Myr before beginning to collapse. When it begins helium burning, the central density and temperature are ∼10 g cm⁻³ and 2.0 × 10⁸ K, respectively. Radiation dominates the pressure in the core, and the adiabatic index, γ_ad, is slightly above 4/3. At this density and temperature, even the most energetic photons in the Maxwellian distribution cannot create electron–positron pairs and cause γ_ad to fall below 4/3. Pair production therefore does not trigger collapse in the SMS. Instead, the energy density of the thermal photons in the high temperature and very low density of the core begins to affect the gravitational field by contributing to the source term in Einstein's field equation, in effect becoming an additional source of gravity. The ratio of the radiative pressure to the relativistic density P/(ρc²) becomes ∼1–3× 10⁻³. This causes γ_ad in the core to fall below 4/3, and it begins to contract.

Although the temperature and density in the core now rise rapidly and accelerate nuclear burning, they are not high enough to ignite carbon burning, which is the next stage after helium burning. Helium instead begins to burn explosively in the core, burning to carbon first through the triple α reaction and then to ¹⁶O, ²⁰Ne, ²⁴Mg, and ²⁸Si through α capture reactions. These reactions release enough energy to reverse collapse. We show the evolution of radial velocities in the star in Figure 1. The onset and reversal of collapse is evident, together with the formation of the outgoing shock, which propagates outward at several thousand kilometers per second. The shock breaks through the surface of the star 11,880 s after core bounce.

Zoom In Zoom Out Reset image size

The binding energy of the star prior to collapse is 5.76 × 10⁵³ erg, but explosive helium burning releases 6.52 × 10⁵⁴ erg. We show the evolution of density and temperature at the center of the star in Figure 2. They are initially ∼20 g cm⁻³ and 3.6 × 10⁸ K, and then rise to peak values of 360 g cm⁻³ and 8.26 × 10⁸ K at core bounce before beginning to fall. The relative magnitudes of the binding and explosion energies together with the rapidly falling density at the center of the star strongly suggest that the SN completely unbinds the star, with no BH formation.

Zoom In Zoom Out Reset image size

We initialize CASTRO with the KEPLER profile of the collapsing star 600 s prior to maximum core compression. In the beginning of the CASTRO run, the core temperature T_c and density ρ_c are 8.21 × 10⁸ K and 356 g cm⁻³, respectively. Element masses as a function of radius are shown in Figure 3. The mass of the helium core is ∼31,000 M_☉, and it also contains ²⁸Si, ²⁴Mg, ²⁰Ne, ¹⁶O, and ¹²C. As in KEPLER, the collapse is reversed by nuclear energy release, but violent fluid instabilities now arise in the core during the explosion. Helium burning creates pressure gradients that are opposite to the gradients of density and mass fraction, which in turn gives rise to Rayleigh–Taylor instabilities. These instabilities are more violent than those in non-reactive flows because they mix the fuel with hot ash that then quickly burns the fuel, which in turn enhances the pressure gradient and drives more mixing. This explosion released 8.82 × 10⁵⁴ erg, slightly more than in the 1D model because of the additional burning due to mixing.

Zoom In Zoom Out Reset image size

Burning continues until the shock reaches the hydrogen envelope. Table 1 shows elemental abundances before and after the explosion. Energy release in this SN is primarily due to the burning of helium, oxygen, and neon, which consumes 757 M_☉, 1116 M_☉, and 291 M_☉, respectively. Burning in turn yields 1195 M_☉ of ²⁴Mg and 970 M_☉ of ²⁸Si. Few isotopes beyond ²⁸Si are synthesized, with ≪1 M_☉ of ⁵⁶Ni being formed, not enough to be detected. Since the hydrogen envelope of the SMS is not as extended as in red supergiants, the shock does not plow up much mass when it crashes into it. The reverse shock is therefore weak and does not drive further mixing. We show the degree to which ¹⁶O, ²⁴Mg, ²⁸Si, and ³²S have mixed by the time the shock has broken out of the surface in Figure 4. Most mixing has now ceased. The major drivers of mixing in this model are fluid instabilities that emerge during burning. Mixing can cause heavy elements to be dredged up from deeper layers and appear in the SN spectra at early times.

Zoom In Zoom Out Reset image size

In KEPLER models that exclude GR, the 55, 500 M_☉ star collapses to a BH instead of exploding because its core continues evolving through helium burning, then encounters the pair instability before it ignites carbon burning. Once the collapse of core triggered by pair instabilities happens, the nuclear burning is not able to halt the collapse at this time and the entire star eventually dies as a BH. Likewise, if the star is 56,000 M_☉  it collapses to a BH, even when GR is included. This suggests that there is a narrow window in mass around 55,500 M_☉ for the progenitor masses of general relativistic supernovae (GSNe). More masses need to be tested and more extensive nuclear reaction networks may be needed to capture all the nucleosynthetic pathways in GSNe. Our KEPLER and CASTRO simulations show that GSNe are only weakly dependent on dimensionality.

GSNe are different from pair-instability SNe (PSNe; Heger & Woosley 2002; Scannapieco et al. 2005; Chen et al. 2014) in every aspect: what triggers the collapse, what drives the explosion, the degree of mixing, and ⁵⁶Ni production. We compare properties of PSNe and GSNe in Table 2. In GSNe, only a trace amount of ⁵⁶Ni (≪1 M_☉) is produced at the edge of the oxygen-burning shell by α capture. PSNe can synthesize up to 50 M_☉ of ⁵⁶Ni through explosive ²⁸Si burning. The light curves of GSNe are mainly powered by the thermal emission by hot ejecta and the conversion of kinetic energy into entropy by the shock, not ⁵⁶Ni decay. Large masses of elements with atomic masses between ¹²C and ²⁸Si are synthesized during the explosion, and they are all dispersed into the surrounding medium.

We have discovered that Pop III stars with masses of 55,500 M_☉ may explode as SNe instead of collapsing to BHs. GR effects, rather than the pair instability, trigger the explosion, which at ∼9 × 10⁵⁴ erg is the most energetic thermonuclear SN known. Mixing in 2D enhances burning and thermonuclear yields during the explosion, which completely unbinds the star and leaves no compact remnant. Energy release is primarily due to explosive helium burning after the onset of central helium burning. The explosion yields mostly elements between carbon and silicon, with almost no iron group elements. Besides enhancing burning, mixing during the SN can dredge up heavy elements from deeper layers and cause them to appear at earlier times in spectra. Montero et al. (2012) recently studied the relativistic collapse of SMS and discovered that explosions of SMS occurred even for stellar masses up to 5 × 10⁵ M_☉, which are powered by the hydrogen burning due to hot CNO cycle. In their simulations, they use simplified stellar models with a single fluid of fixed hydrogen, helium, and CNO metal abundance. More importantly, all of their exploded models are associated non-zero metallicity; however, there is no metal inside the Pop III stars when they are born. Hence the hot CNO cycle may not occur for the Pop III SMS early on. Instead, our SMS model evolves from the main sequence through collapse to explosion with the detailed stellar physics that self-consistently produces CNO from the primordial elements. This provides a more realistic pre-supernova model. The differences in pre-supernova models as well as other microphysics used can significantly alter the explosion mechanics that may explain the discrepancies between their findings and ours.

Our simulations are approximate for several reasons. First, we do not model the evolution of the protostar from much lower masses. Instead, the star is initialized at the beginning of the main sequence in KEPLER. Next, we did not evolve the star under the heavy ongoing accretion that gave birth to the star. We also did not consider stellar rotation, which could lower the mass at which the star explodes (e.g., Chatzopoulos & Wheeler 2012; Chatzopoulos et al. 2013). There is growing evidence that many Pop III stars are born with high rates of spin (Stacy et al. 2011; Greif et al. 2012). Finally, we did not include UV feedback from the star, which could regulate its accretion rates. However, in some cases, luminosity from the star is thought to terminate accretion (Johnson et al. 2012), so our assumption that the star has a constant mass is plausible. Efforts are now underway to survey Pop III protostellar evolution under a variety of accretion rates and SMS evolution with much larger and finer grids in mass.

What effect would such an explosion have on the protogalaxy that gave birth to the star? The answer depends on the ambient density of the SN. If the SMS grows by accretion through a disk, its UV radiation will break free of the disk and ionize the surrounding envelope. In these circumstances, the star would explode in low densities, and Johnson et al. (2013b) show that the SN would drive all the gas from the protogalaxy, perhaps engulfing nearby protogalaxies with metals. Much of this material would later fall back to the halo on timescales of 50–100 Myr and form stars. If accretion is spherical, then the star will not ionize its envelope and it will explode in very high densities (Johnson et al. 2012; Hosokawa et al. 2013; Schleicher et al. 2013). In this scenario, Whalen et al. (2013b, 2013c) find that much of the energy of the explosion is promptly radiated away by bremsstrahlung X-rays and inverse Compton scattering of cosmic microwave background photons, but that the SN remnant still expands to roughly the virial radius of the halo, ∼1 kpc, and then collapses back into the halo. Collapse thoroughly mixes baryons in the halo with metals and may trigger a starburst that would easily distinguish this protogalaxy from its dimmer and less rapidly evolving neighbors.

Radiation hydrodynamical simulations show that the GSN would be visible at z ≲ 20 to future all-sky NIR surveys by Euclid, the Wide-Field Infrared Survey Telescope, and the Wide-field Imaging Surveyor for High-Redshift (Whalen et al. 2013a). The wide survey fields and high sensitivities of these missions would enable them to detect such an event anywhere in the universe, even if their numbers are small. Furthermore, a single GSN can create 100 times the chemical yield of a PSN. However, unlike PSNe, which synthesize more iron group elements, GSNe would mainly enrich primordial gas with elements from ¹²C to ²⁸Si. Traces of GSNe might therefore be found in early galaxies that are ⁵⁶Fe deficient but enhanced with ¹²C and ¹⁶O (Keller et al. 2014). Whether in future NIR campaigns or in the fossil abundance record, the detection of GSNe in the early universe may soon signal the birth of SMBH seeds and the first quasars.

The authors thank the anonymous referee for reviewing this manuscript and providing insightful comments, and the members of CCSE at LBNL for help with CASTRO. We also thank Volker Bromm, Dan Kasen, Lars Bildsten, John Bell, and Adam Burrows for many useful discussions. K.C. was supported by an IAU Gruber Fellowship, a Stanwood Johnston Fellowship, and a KITP Graduate Fellowship. A.H. was supported by a Future Fellowship from the Australian Research Council (ARC FT 120100363). D.J.W. was supported by the Baden–Württemberg–Stiftung by contract research under the programme Internationale Spitzenforschung II (grant P-LS-SPII/18). All numerical simulations were performed at the University of Minnesota Supercomputing Institute and the National Energy Research Scientific Computing Center. This work was supported by the DOE grants DE-SC0010676, DE-AC02-05CH11231, DE-GF02-87ER40328, and DE-FC02-09ER41618, and by NSF grants AST-1109394 and PHY02-16783. 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.