Stability analysis of supermassive primordial stars: a new mass range for general relativistic instability supernovae.

The HOSHI code is a stellar evolution code which, in this work, solves the hydrodynamical equations in post Newtonian (PN) gravity using a Henyey type implicit method (Newton-Raphson method), taking the density, entropy, radius, and luminosity as the independent variables (Takahashi et al., 2016, 2018, 2019; Yoshida et al., 2019; Nagele et al., 2020). Although we describe HOSHI as a stellar evolution code (it solves the stellar structure equations, including energy transport), the code is able to follow hydrodynamical evolution to a degree, though it does not include a shock capture scheme. This, along with the lack of full GR, is why we require HYDnuc to model the GRSN. HOSHI includes a nuclear reaction network (52 isotopes), neutrino cooling, mass loss (though it is minimal for non rotating primordial stars), and rotation. The equation of state includes contributions from photons, averaged nuclei, electrons, and positrons. We use the Rosseland mean opacity of the OPAL project (Iglesias & Rogers, 1996) and solve the Saha equation to determine the ionization of hydrogen, helium, carbon, nitrogen, and oxygen. Several effects of convection are modeled including chemical mixing through diffusion in convective and semi-convective regions. Finally, energy transfer due to convection is treated according to 1D mixing length theory.

The star is initiated with Log $T_{c}$ $<$ 7.7 in a high entropy state relative to ZAMS and has a primordial composition except for deuterium which we have removed to avoid the proto-stellar burning phase, which is not important for the evolution (see e.g. Hosokawa & Omukai 2009). Because we initialize the star as a supermassive protostar, its structure is very nearly that of an $n=3$ polytrope. The star will undergo a moderate period ($\sim 10^{12}$ s) of contraction. The p-p chain is insufficient to stop this contraction and so the star must reach a central temperature around Log $T_{c}$ = 8.2 at which point the triple alpha reaction produces enough carbon to ignite the CNO cycle, which stabilizes the star. The star continues to be supported by the CNO cycle until hydrogen is exhausted, at around one million years. The star then contracts, with the central temperature increasing to Log $T_{c}$ = 8.5 before helium burning stabilizes the star. The stars in this paper become unstable to the GR radial instability (Sec. 2.2) in late hydrogen burning or in helium burning. Fig. 1 shows a Kippenhahn diagram for the M=$3\times 10^{4}$ ${\rm M}_{\odot}$ model. At the beginning of hydrogen burning, the star is nearly fully convective (diagonal hatches), though towards the end of this period, a core and envelope form. The core remains convective until the instability while the envelope has several convective layers separated by areas of non convection (dots) or semi-convection (crosses). Hydrogen shell burning proceeds in the non convective layer just above the core. The onset of the instability is sudden, as is the case with the pair instability, and is not precipitated by major changes in the composition (see color-bar) as in the case with core collapse.

HOSHI uses the first order PN approximation to the Tolman Oppenheimer Volkoff (TOV) equation. However, unlike in Chen et al. (2014); Nagele et al. (2020), we include the correction of energy to density:

where $\rho_{0}$ is the baryonic density and $\epsilon$ is the specific internal energy in units of ergs g^-1. We use the convention that the rest mass energy due to the mass excess of isotopes is included in the internal energy. $\epsilon/c^{2}$ is between 0.01 and 0.001 throughout the star and its inclusion is necessary to correctly model the SMS envelope. Fig. 2 shows the results of the stellar evolution calculation compared to Nagele et al. (2020).

Beyond correctly modeling the SMS envelope, the inclusion of internal energy is necessary for consistency between HOSHI and HYDnuc. If we only use the baryonic density when calculating the TOV PN terms in HOSHI, then the $\mathcal{O}(1\%)$ difference in the gravity will perturb the star too strongly in HYDnuc, causing most models— even some stable configurations— to collapse to black holes.

This work contains twenty supermassive stars with different masses, where we have chosen those masses to be centered around the explosion window. The explosion window is heavily dependent on when the GR instability triggers, as an explosion is only possible if the instability occurs during helium burning.

The question of when a SMS collapses due to the GR instability is a challenging one because the collapse takes place on timescales far smaller than the typical timestep of an evolutionary calculation. In addition, the TOV equation lacks the dynamical GR terms which make the collapse so rapid. It would likely require a relativistic stellar evolution code to fully address the problem of SMS collapse. In Nagele et al. (2020), we relied on the PN stellar evolution code to determine when collapse would occur, and our results were in rough agreement with those of Chen et al. (2014). In this work we perform a stability analysis on the normal modes of radial perturbations of a star in GR (Chandrasekhar, 1964).

Consider an infinitesimal, radial, Lagrangian perturbation which varies in time ($t$) as $\xi\propto e^{i\omega t}$ for $\omega^{2}\in\mathbb{R}$. In Newtonian gravity, this perturbation obeys the equation (Shapiro & Teukolsky, 1983)

where $P$ is the pressure, $r$ is the radius, $\rho_{0}$ is the baryonic density and $\Gamma_{1}$ is the local adiabatic index at constant entropy (s):

Solving this differential equation for $\xi$ and $\omega$ is a Sturm–Liouville eigenvalue problem. Finding any solution with $\omega^{2}<0$ is a sufficient condition for instability, as the motion of the perturbation will be exponential. Sturm–Liouville equations have the property that a sequence of solutions exist

corresponding to $\xi_{i}$s where $i$ is the number of nodes in the perturbation. Because of the above property, a necessary condition for instability is

The corresponding equation in GR is (Chandrasekhar, 1964)

where $a,b$ are the metric coefficients as in Haemmerlé (2021) and the density is defined in Eq. 1.

In order to solve Eqs. 2,6 for $\omega^{2}_{0}$, we adopt an iterative method similar to that outlined in exercise 6.11 of Shapiro & Teukolsky (1983) (for a detailed discussion, see Appendix A). For a given stellar profile ($r,m_{r}(r),\rho(r),P(r)$), we choose a value of $\omega^{2}$ and integrate Eq. 6 (or Eq. 2, the procedure is identical) to find $\epsilon$. This is done twice, once starting from the center and once starting from the stellar surface. If the two solutions match at some test radius, then $\omega^{2}$ is a solution to Eq. 6. If not, we repeat the procedure, extrapolating new values of $\omega^{2}$ based on the Wronskian of the two solutions at the test radius. Once a solution is found, we check if it is the fundamental mode by determining the number of nodes in $\xi$. If the number of nodes is zero, we have found $\omega_{0}^{2}$, and if it is not, we repeat the procedure for a lower initial value of $\omega^{2}$ (see Eq. 4).

In order to test our method, we construct numerical Lane-Emden Polytropes. We check that the polytropes satisfy $\omega^{2}=0$ at $\Gamma_{1}=4/3$ in the Newtonian case and $\Gamma_{1}=4/3+\kappa\frac{2GM}{Rc^{2}}$ in the relativistic case, where $M$ is the mass of the star, R the radius, and $\kappa$ is a constant determined numerically (Chandrasekhar, 1964). Fig. 3 shows the accuracy of these relations for increasing numerical resolution of the polytropes. We also verify $\xi_{0}\propto r$ for these values of $\Gamma_{1}$, a condition that should hold for all polytropes.

Once an unstable model is found in the stellar code (e.g. Fig. 4), the calculation is mapped to the hydrodynamical code (HYDnuc, Sec. 2.3). There, the star will either begin to collapse or stabilize due to nuclear burning. It is important to note that the stability analysis considers the stellar structure at a single moment, and cannot account for the energy generated by nuclear burning.

Thus, from the start of the HYDnuc calculation, there is a competition between the growing perturbation of the unstable model and nuclear energy generation. For lower mass models ($M\leq 2.7\times 10^{4}\;{\rm M}_{\odot}$), energy generation can sometimes stabilize the star. Fig. 5 shows the stability analysis applied to a HYDnuc model which stabilizes (M = 2 $\times 10^{4}$ ${\rm M}_{\odot}$, first instability). The stability analysis assumes zero velocity, so we cannot rely on it too heavily in a dynamical scenario, but it can be illustrative. The model is initially unstable and begins to contract; the star continues to contract, but the stability analysis fluctuates between ’stability’ and ’instability’. Here, ’instability’ and ’stability’ refer to whether the contraction is growing exponentially or not. In order for the star to become stable in the usual sense, nuclear burning must increase to counteract the perturbation, which occurs around $10^{5}$ s. Afterwards, the temperature decreases and reaches a new equilibrium which is higher than that of the initial model.

For the high mass models (and later time low mass models), the perturbation overcomes the energy generation, and the SMS moves onto a collapsing phase that triggers rapid alpha capture burning and may lead to a GRSN.

If the result of the hydrodynamical code is that the star stabilizes (as in Fig. 5), we map the next instability in the stellar code to HYDnuc and repeat until a model either explodes, pulsates, or collapses. In the case of e.g. 2 $\times 10^{4}$ ${\rm M}_{\odot}$, the first such model has burned most of its helium (Table 1).

HYDnuc is a 1D GR implicit Lagrangian code (Takahashi et al., 2016, 2018, 2019; Yoshida et al., 2019; Nagele et al., 2020, 2021) based on the nuRADHYD code (Yamada, 1997; Sumiyoshi et al., 2005; Nagele et al., 2021) which includes a Boltzmann solver for neutrino transport not present in HYDnuc. HYDnuc uses a Roe-type approximate linearized Riemann solver. The independent variables are density, velocity, internal energy, entropy, electron fraction, radius, baryonic mass, enthalpy, and the two components of the Misner Sharp metric (Misner & Sharp, 1964). Internal energy changes primarily via energy generation from a nuclear network or cooling from thermal neutrino reactions. The equation of state is the same as in HOSHI.

In order to transport a model from HOSHI to HYDnuc, we first determine the radial values of the HYDnuc mesh according to the frequency function in Appendix B of Takahashi et al. (2019). This method is used to provide appropriate resolution to both the core and envelope of the SMS. Other variables are then inferred using linear radial interpolation of the log values of those variables, which we have found is the most accurate method of ensuring the fidelity (as a function of enclosed mass) of the remapping. Care must be taken at this step so as not to introduce unphysical perturbations. This is the same reason that we introduced the relativistic energy density to HOSHI (Sec. 2.1). However, an unavoidable coordinate perturbation due to numerical error will always be present when changing mesh definitions, and this is a potential source of error in our simulation.

We continue the calculations for the explosions until there are convergence problems due to poor spatial resolution in the outer regions of the star. This typically occurs when the shock reaches $10^{15-16}$ cm. For the pulsations, a similar stopping condition is reached, but in this case we excise the ejected material from the simulation and verify that the remaining material becomes hydrostatic. For the models which collapse, we continue the calculation to a central temperature of 1 MeV, at which point neutrino heating begins to play a role.

In order to set the numerical parameters for HYDnuc, we perform several numerical convergence tests using the the explosion energy, which is defined as the total energy when the shock reaches the stellar surface. Fig. 6 shows the explosion energy as a function of mesh point number, total isotope number, and $\mathcal{V}$, the limit on the maximum variation of the independent variables

where $x_{i}$ is one of the independent variables of HYDnuc (Yamada, 1997), j is the mesh point number, and k is the time-step. If the limit in Eq. 7 is violated for a time-step k, that time-step is repeated with a reduced value of dt. Thus, a smaller $\mathcal{V}$ will require more time-steps which increases the time resolution of the simulation, making it more physically accurate. The simulations with greater resolution tend to reach slightly higher temperature, which correlates to an increase in $E_{\rm exp}$. In this paper, we use $\mathcal{V}=10^{-5}$, a 61 isotope network, and 767 mesh points. The 61 isotope network includes the p-p chain, CNO and hot CNO cycles, triple alpha, detailed alpha process until silicon, a more basic alpha process up to nickel, and photodissociation of heavy elements.

The explosion energy depends on $\mathcal{V}$ and mesh point number in straightforward ways, but the dependence on isotope number requires explanation. The main driver of the explosion is alpha capture reactions, but not all of these reactions proceed at the same rate. In particular, the carbon alpha capture rate is lower than the reactions involving ²⁰Ne and ²⁴Mg; we use 1.5 $\times$ the rate of Caughlan & Fowler (1988), but have verified that the explosion energy depends very weakly on this reaction rate. $\mathcal{O}(10\%)$ of the mass of the star is carbon, and very little of this is burned via carbon alpha capture on the timescale of the explosion. However, if nucleons are present, catalysis enhances the carbon alpha capture rate with

In the networks with isotope number less than 61 (Fig. 11, left panel), a reservoir of free nucleons is built up during the explosion (Appendix B), and these then serve as the catalyst for carbon burning. In the networks with higher isotope numbers, however, the nucleons are absorbed in reactions such as

Indeed, aluminium is of particular importance (Fig. 11, right panel) because the inclusion of its isotopes is the only difference between the 58 isotope network and the 61 isotope network, and from Fig. 6, we can see that this is the isotope number where the explosion energy converges. Appendix B contains a steady state calculation verifying this explanation.

Usually, the inclusion of a larger network increases the nuclear energy generation, but in this case, a threshold number of elements is required to properly follow the nucleonic reactions. Nagele et al. (2020) used a 49 isotope network and Chen et al. (2014) used a 19 isotope network, so it is possible that the explosion energies found in those works are overestimated. In this work, the maximum explosion energy is a factor of 3-4 smaller than in the previous works. One reason for this is that the lower mass means there is less fusion material, but the effect of the nuclear network also plays a role. For a more detailed discussion of how this work compares to previous ones, see Sec. 4.

First, and most common in the literature (e.g. Fuller et al., 1986; Umeda et al., 2016; Woods et al., 2017) is the polytropic criterion. SMSs have high entropy and are supported mostly by radiation pressure and this invites analytic approximations. In particular, it is often assumed that the SMS core is very nearly an $n=3$ polytrope.

The explosion in the next section involves explosive helium burning, so we will calculate the instability condition for a polytrope consisting of pure helium, as a function of the mass of the helium core ($M_{\rm He}$) which is taken from the stellar evolution calculation. Once the mass of the helium core is known, the radiation entropy may be approximated as $s_{r}\propto M_{\rm He}^{1/2}$ (Shapiro & Teukolsky, 1983). From here, the polytropic constant is determined

Next, we calculate the outer radius at which the star will be unstable, $R_{\rm crit}$ by setting the SMS $\Gamma$ equal to the general relativistic $\Gamma_{1}$,

where $\beta\approx 4.3/\mu(M/{\rm M}_{\odot})^{-1/2}$ is the ratio of the gas pressure to total pressure, and $\kappa=2.249$ for $n=3$ (this form of $\kappa$ differs from Chandrasekhar (1964) by a factor of 2). Eq. 11 can be solved for $R_{\rm crit}$ as a function of $M_{\rm He}$, so that we finally arrive at an expression for the critical density (Shapiro & Teukolsky, 1983)

as a function of $M_{\rm He}$, where $\xi_{1}=6.897$ is determined numerically.

Thus, we have an expression for the critical density of a purely helium SMS core as a function of $M_{\rm He}$. By construction, the SMS models in this paper are near to this point according to the GR stability analysis from Sec. 2.2. Fig. 7 shows that, for these models, the central density when the star is pure helium is more than an order of magnitude below the critical density. So, the polytropic criterion underestimates the GR instability in comparison to Sec. 2.2.

Next we compare our method to the results of Haemmerlé (2021), who also evaluate Eq. 6, though they make the simplifying assumption $\xi\propto r$. This assumption is valid for polytropic stars and possibly for higher mass hydrogen SMSs, but for our models, the perturbation is not always proportional to the radius (see Appendix A).

Fig. 8 shows the helium mass fraction at the first instability — that is, the first model in the HOSHI calculation which is unstable — determined by either method. The central helium mass faction for an explosion or pulsation is roughly 0.1 < $X_{\rm c}(^{4}$He) < 0.7, c.f. Table 1. Our method uses a necessary condition for instability while the method of Haemmerlé (2021) uses a sufficient condition. This means that there are models which would be stable according to the method of Haemmerlé (2021), but not according to our method. Furthermore, any model which is unstable according to the method of Haemmerlé (2021) will also be unstable according to our method. Thus, our method will find an instability earlier in the evolutionary calculation than the method of Haemmerlé (2021). The longer the time before the instability in the evolutionary calculation, the higher the mass range for the explosion, because lower mass models which would explode if the instability occurred earlier instead burn all of their helium ($M<3\times 10^{4}$ in Fig. 8). This means that the explosion would occur at larger mass, if we were to use the method of Haemmerlé (2021), although that being said, the explosion mass range found by their method would still be significantly smaller than the case of $M\sim 5\times 10^{4}\;{\rm M}_{\odot}$ discussed in Chen et al. (2014); Nagele et al. (2020).

We compare the neutrino light-curves of the lowest (2 $\times 10^{4}$ ${\rm M}_{\odot}$) and highest (4 $\times 10^{4}$ ${\rm M}_{\odot}$) mass models in this study with the results from our previous work (Nagele et al., 2021). The nuRADHYD code is unchanged so the only difference is the amount of time the star spends in the evolutionary stage. Both stars have higher entropy than in our previous work because there is less time for neutrino cooling during the evolutionary stage, and even though the models in this work have different chemical compositions at the start of collapse, they will eventually undergo the same reactions, namely alpha capture until sulfur, then production of calcium through nickel, followed by photodisassociation, first into helium and then into nucleons.

In our previous paper, we identified that many physical quantities scaled with the entropy at fixed density. Thus, the stars in this paper having higher entropy would suggest that they should also have higher temperature, and neutrino luminosity. Both of these turn out to be true, but while the hydrodynamical quantities match the trends in Fig. 8 of Nagele et al. (2021), the neutrino quantities do not. Thus, the neutrino luminosity, and number flux are all increased relative to the previous work, but not by as much as we expected. Of particular interest, the total neutrino number increased by 13$\%$ for 2 $\times 10^{4}$ ${\rm M}_{\odot}$ and 42$\%$ for 4 $\times 10^{4}$ ${\rm M}_{\odot}$. Furthermore, although we would expect the average neutrino energy to decrease because of the higher entropy, they increase by small fractions, 2$\%$ for 2 $\times 10^{4}$ ${\rm M}_{\odot}$ and 11$\%$ for 4 $\times 10^{4}$ ${\rm M}_{\odot}$.

Although these both trend in the right direction regarding detection of the diffuse SMS neutrino background, they are not large enough increases to alter our previous conclusion that if SMSs collapse in this mass range, then the detection of this background is not feasible using current methods.