The Final Fate of Supermassive $M\sim 5\times 10^{4}\;M_{\odot}$ Pop III Stars: Explosion or Collapse?

A 1D stellar evolution code, HOngo Stellar Hydrodynamics Investigator (HOSHI) is used to evolve supermassive population III stars from hydrogen burning to the onset of gravitational instability (Takahashi et al., 2016, 2018, 2019; Yoshida et al., 2019). In this work, the code uses a 49 isotope nuclear network and does not include mass loss. We assume the direct collapse scenario as the formation channel, and assume no accretion.

In this paper, $M_{\rm r}$ refers to the mass inside a given radius $r$, $\rho$ is the density, $\mu$ is the mean molecular mass, $s$ the entropy per baryon, $T$ the temperature, and $X(i)$ the mass fraction of an isotope $i$. Quantities with a subscript c such as $\rho_{\rm c}$ refer to the central value of that quantity.

We set rotation at ZAMS to be $\Omega\approx 10^{-15}$ s^-1. The code then treats rotation self consistently (e.g. Heger et al., 2000, 2005; Meynet & Maeder, 2000; Maeder & Meynet, 2004), but we do not expect $\Omega>10^{-10}$ s^-1 so this is effectively a non rotating case. We investigate the effects of including slow rotation in Sec. 3.4. More detailed descriptions about the progenitor evolution and faster rotating cases are given elsewhere (Umeda et al. 2020 in preparation).

We have also added the first order post Newtonian correction to general relativity in the form of the Tolman-Oppenheimer-Volkoff equation (Oppenheimer & Volkoff, 1939; Tolman, 1939).

where the second line is the post Newtonian correction and includes only terms of order $c^{-2}$. Our results are dependent on including this correction. In the Newtonian case we found no evidence for an explosion.

After the HOSHI calculation, we switch to a hydrodynamics code with nuclear reactions (HYDnuc). This code is based on a 1D spherical, Lagrangian, general relativistic, neutrino radiation hydrodynamics code (nuRADHYD) developed by Yamada (1997) and modified in Sumiyoshi et al. (2005). Nuclear reaction networks were added in Takahashi et al. (2016). In this work, the code uses the same 49 isotope nuclear network as the stellar evolution code and it has an NSE solver for high temperature, as well as neutrino cooling (Itoh et al., 1996) which is the same as in the stellar evolution code. Since the HYDnuc code does not include radiative transfer of photons, the envelope will collapse if the model is static. For this reason, we can only switch to the hydrodynamics code once dynamical collapse of the core has begun. Switching earlier than this yields collapse starting in the envelope instead of core collapse.

We report the results of this code partially in terms of the energies of the star. Specifically, the internal energy is

where u is the internal energy per unit mass, the gravitational energy is

where $g_{\rm effective}$ is the Post Newtonian approximation of the TOV equation (Eq. 1), the kinetic energy is

where $v$ is the radial velocity. The total energy is the sum of these three energies

Furthermore, we include the time integral of the energy generation rate

from the start of the HYDnuc calculation

At some point during the HOSHI calculation, the star becomes unstable due to the general relativistic instability (Chandrasekhar, 1964) and it is possible for the entire star to experience dynamical collapse. This collapse induces rapid ⁴He burning in the central region and may cause an explosion (Chen et al., 2014). Our results are sensitive to the condition for switching between codes; explosion or non-explosion is mainly determined by the ⁴He mass fraction of the central convective core when we switch to the HYDnuc code. If we switch too early, the ⁴He mass fraction is larger and the model tends to explode.

We will now explain our criterion for the switch between codes; to first approximation, the GR instability coincides with the onset of dynamical collapse. There are a few ways of calculating the condition for the instability, but here we will assume that the star is radiation dominated and that the total entropy per baryon is constant in the isentropic core. Then, the entropy is approximately the radiation entropy:

(Shapiro & Teukolsky, 1983). This means the star has the mass distribution of an $n=3$ polytrope. Then, the critical central density above which the star is unstable is (e.g. Fuller et al., 1986)

where $M_{\rm core}$ is the mass of the isentropic core.

Using $\rho_{\rm c}>\rho_{\rm crit}$ as a criterion for switching codes, however, does not give a realistic result. This is because up to this point we have neglected the effects of nuclear burning. In reality, strong nuclear burning supports the core for some time before dynamical collapse. In our case, strong helium burning prevents rapid collapse even after the condition $\rho_{\rm c}>\rho_{\rm crit}$ is satisfied. As shown in Fig. 1, after the star satisfies $\rho_{\rm c}>\rho_{\rm crit}$, it takes more than $10^{12}$ seconds until collapse and during this period, there are changes in the physical quantities. Specifically, the central ⁴He mass fraction and the temperature change significantly (Fig. 1).

However, after $10^{9}$ seconds until collapse, the ⁴He mass fraction is mostly unchanged. In practice, we switch codes around $10^{6-7}$ seconds before collapse, when the timestep is $\Delta t_{\rm HOSHI}\sim 10^{2-3}$ s, though we expect the results to be the same for switching at any time after $10^{9}$ seconds until collapse; note that the time step is essentially monotonically decreasing in the flat section of Fig. 1. The Kelvin-Helmholtz time scale at this stage, $\tau_{\rm KH}$, is on the order of $10^{9}$ s.

In summary, we switch codes when $\Delta t_{\rm HOSHI}\ll\tau_{\rm KH}$ which occurs after $\rho_{\rm c}>\rho_{\rm crit}$ and in between these two times, the composition of the star can change significantly.

For most of our models, the star collapses too quickly for the internal energy to become significantly larger than the gravitational energy (Fig. 3). During the helium burning period, central temperature increases and nuclear reactions produce a large amount of energy (see internal energy in Fig. 3), which in some cases is enough to cause the total energy to become positive. However, the total energy does not stay positive for long enough to cause an explosion, and the total energy never becomes significantly larger than the kinetic energy. 42 seconds before the simulation terminates, the total energy once again becomes negative because of a rapid increase in the negative of the gravitational energy, $-E_{\rm grav}$. The calculation terminates due to instability in calculating the NSE around $g_{00}=0.35$ (the time component of the metric) for the central mesh.

We calculate the neutrino cooling luminosity, which peaks at $L_{\rm\nu,cooling}=1.40\times 10^{58}$ ergs s^-1. Because HYDnuc only contains neutrino cooling, this should be viewed as an upper limit for the true luminosity, and we plan to perform a realistic calculation of the neutrino light curve in the future using nuRADHYD.

It is also possible that the collapse process ejects some heavy elements. In Uchida et al. (2017), which simulated the collapse of a rotating SMS in the helium burning phase, a torus of material was ejected from the collapsing supermassive star; such a torus could contain ⁵⁶Ni or other heavy elements if it is ejected in the 10 seconds before black hole formation. However, because the heavy elements are located in the core and because of the slow rotation rate, we do not expect this effect to apply to our models.

We found one realistic exploding model (Fig. 4). This model had an initial rotation which caused the ⁴He mass fraction at dynamical collapse to be higher than in the non rotating case for this mass (See Sec. 3.4).

Even though this model has a lower helium mass fraction than the model in the previous section ($X_{\rm c}(^{4}$He$)=0.190$ compared to 0.199), it is able to explode because it has a lower mass (Fig. 2) and because of additional weight from the helium envelope (see Sec. 3.4). The initial temperature and density are similar to Fig. 3, but the kinetic energy is lower (Fig. 4). This allows the internal energy to overcome the kinetic energy and cause an explosion with explosion energy $E_{\rm tot}=5.45\times 10^{54}$ ergs. The final outward velocity exceeds the escape velocity in the outer regions of the star, with a maximum value of $v_{\rm max}=3.11\;v_{\rm esc}$. However, the inner eighty percent of the star ($M_{\rm r}<40200$ $M_{\odot}$) does not exceed the escape velocity. Thus the final fate of this object after the explosion is unknown, but it is clear that a large compact object does not immediately form.

The large explosion energy is caused by nuclear burning in a significant proportion of the star during the explosion. In Fig. 4, helium burning is sustained for 10000 seconds and during this time more than $1500\;M_{\odot}$ of ⁴He and $3300\;M_{\odot}$ of ¹⁶O is burnt (Table 1). The processes responsible for the latter is alpha capture of ¹⁶O, then of ²⁰Ne, and— in the inner ten percent of the star— of ²⁴Mg (Fig. 5). The explosion energy mostly comes from nuclear burning with $\Delta E=5.10\times 10^{54}$ ergs.

In Fig. 5, we can see the initial and final isotope distributions. Since the temperature never gets higher than Log $T_{\rm c}$ [K] $=8.7$, no elements beyond ²⁸Si are produced. The ⁴He-¹⁶O core is initially convective, as dynamical collapse occurs and central temperature increases, alpha process reactions occur in the center, disrupting the convection.

We should note that the timescale for this explosion is very long ($\sim 10000$ s compared to $\tau_{\rm dynamical}\sim 500$ s for the isentropic core.). In other types of supernovae, such as core collapse or pair instability supernovae, most of the burning takes place on time scales several orders of magnitude shorter. This explosion is caused by the quantity of material burnt, as opposed to reactions producing large amounts of energy quickly.

In the above subsections, we have examined two models, one of which collapses and one of which explodes via rapid helium and oxygen burning. In our simulations, the primary determinant of the fate of the model is the central ⁴He mass fraction when we switch codes. Models with $X_{\rm c}(^{4}$He$)$ over a certain threshold value explode and those with $X_{\rm c}(^{4}$He$)$ below this value collapse. We have observed that the threshold value tends to increase with mass (Fig. 2), although this is not strictly true due to effects from the envelope which introduce some randomness into our calculations (Sec. 3.4). For the mass range considered in this paper, the threshold value is between $0.2<X_{\rm c}(^{4}$He$)<0.4$. As demonstrated in Fig. 1, toward the end of the calculation in HOSHI, the helium mass fraction in the core decreases due to helium burning, so the later we switch codes, the less helium will be present in the core and thus the star will be more likely to collapse.

We do not find any evidence for an explosion among non rotating stars in the mass range $M=35000-60000\;M_{\odot}$. However, we do notice a peak in the central ⁴He mass fraction when changing codes (Fig. 6) which occurs at $M=55000\;M_{\odot}$. The location of this peak is near to the exploding model found by Chen et al. (2014) at $M=55500\;M_{\odot}$, though their model had a higher helium mass fraction.

The existence of the peak at $M=55000\;M_{\odot}$ is interesting and we will now attempt to explain the decrease in helium mass fraction on either side of $55000\;M_{\odot}$. For $M<55000\;M_{\odot}$, Eq. 9 indicates that higher mass means a lower critical density, so that, for increasing mass, the star has a lower central density and— assuming no large changes in entropy– central temperature when $\rho_{\rm c}=\rho_{\rm crit}$. Lower central temperature in turn means that less nuclear burning has occurred. Thus as we increase the mass, we also expect to increase $X_{\rm c}(^{4}$He$)$ at the time of collapse. This can be seen directly in the top panel of Fig. 1, where the initial temperature decreases with increasing mass. Because the rate of helium burning scales with temperature, any increase in temperature is mirrored by a decrease in the central helium mass fraction (the bottom panel of Fig. 1). Thus to have a large helium mass fraction and a higher chance of exploding, a star needs to have a low temperature at the time of collapse. If this $M<55000\;M_{\odot}$ trend continued to higher mass, we would eventually find models that exploded; however, around $M=55000\;M_{\odot}$, it is supplanted by a different trend.

For $M>55000\;M_{\odot}$, when the instability condition is satisfied, $X_{\rm c}(^{4}$He$)\approx 1$ and ⁴He burning has not yet fully turned on. Without full helium burning, the core cannot support itself against the GR instability (Eq. 9) as efficiently and it rapidly contracts to a higher temperature. In the top panel of Fig. 1, the temperature of the $M=56000\;M_{\odot}$ model starts lower, but the star is stable for longer than the other models. This is because the sudden jump in temperature after the GR instability releases a large amount of energy which supports the star and allows it to burn helium for significantly longer than the lower mass models (Fig. 1). During this period of stability, helium burning continues and this reduces the central helium mass fraction from $X(^{4}$He$)\approx 1$ down to $X(^{4}$He$)\approx.02$

This effect explains why the dropoff of helium mass fraction at dynamical collapse in Fig. 6 is so sharp for masses greater than $55000\;M_{\odot}$. Stars in this mass range have almost zero helium available by the time of dynamical collapse, because their helium was burnt during the stable period which is— almost paradoxically— caused by the jump in temperature after the instability equation is satisfied.

As mentioned above, Chen et al. (2014) found an explosion for a model with mass $M=55500\;M_{\odot}$, which is close to the peak we identify in this paper. For comparison to the model of Chen et al. (2014), we have created an artificial explosion by switching codes much earlier for $55000\;M_{\odot}$; this mass is the peak discussed above and it comes the closest to exploding of any of our non rotating models. The artificial model is slightly more energetic then the explosion of Chen et al. (2014), but shares the basic characteristics (Table 2). We are not sure of the cause of the discrepancy in results, but we note that the exploding model in Chen et al. (2014) seems to have $X_{\rm c}(^{4}$He$)\sim 0.4$; such a model would likely explode using our HYDnuc code, but we saw no evidence from the stellar evolution calculation that a model with such a high helium mass fraction is possible.

Although rotation cannot be included in the calculation during dynamical collapse, it can be implemented in the stellar evolution calculation. It is reasonable to assume that rotation would stabilize the star and thus prolong the period of ⁴He burning before the onset of collapse. This in turn would greatly reduce the chances of an explosion. However, we find that increased rotation affects the onset of collapse in a non standard way (Fig. 7). This figure shows that slow rotation may in some cases lead to larger $X_{\rm c}(^{4}$He$)$ than in Fig. 6 (the leftmost point in Fig. 7 is the non rotating case, see Sec. 2.1). The increased helium mass fraction can then result in an explosion, such as the case of $\Omega=10^{-10.4}$ s^-1 (Fig. 4).

The behavior in Fig. 7 is likely due to the effect of convection in the hydrogen and helium envelope; larger convective regions can form a heavy shell just above the core (Fig. 8). In Fig. 8, the two models shown have similar central helium mass fractions, which is the quantity that we have analyzed up to this point. However, the model which collapses has $M_{\rm final}(^{4}$He$)=3568\;M_{\odot}$ in the 5000 $M_{\odot}$ outside the core, whereas the explosion model has $M_{\rm final}(^{4}$He$)=4057\;M_{\odot}$. The extra mass in this shell increases the effective weight of the core and can induce collapse earlier than if only the weight of the core is considered. This early collapse can then result in a higher helium fraction and an explosion.

As mentioned above, the size of convective regions in the envelope is stochastic, and this explains the randomness present in Fig. 7. Since our investigation of rotation was not very fine grained in the mass or rotation spaces, it is likely that there are other models which explode when slow rotation is introduced, though we suspect they are rare. Finally, the decrease of the threshold $X_{\rm c}(^{4}$He$)$ value with decreasing mass (Fig. 2) suggests that it is easier for an explosion to occur at lower mass, perhaps even below those considered in this paper.