\section{Introduction}
\label{introduction}


For much of the history of astronomy, the post recombination early universe has been inaccessible to observation. While this state of affairs still holds broadly, inroads are being made. Observers have detected a star at redshift 6 \citep{welch2022}, galaxies as early as redshift 11 \citep{oesch2016}, a long gamma ray burst (GRB) at redshift 8 \citep{tanvir2009}, and quasars at redshifts 6 and 7 \citep{mortlock2011,wu2015,banados2018,matsuoka2019,wang2021}. These high redshift observations have greatly increased our knowledge of the early universe, but they also raise questions--- most notably, where did the high redshift quasars and their supermassive black hole (SMBH) engines come from?



Several theories have been put forward to explain the existence of SMBHs so soon after the big bang \citep[e.g.][]{rees1984,inayoshi2020}. These include, but are not limited to, the direct collapse scenario \citep{bromm2003}, super-Eddington accretion onto solar mass black holes \citep{haiman2001}, and rapid mergers of either primordial \citep{bean2002} or astrophysical black holes \citep{omukai2008}. Of these scenarios, the direct collapse scenario may be the easiest to test observationally; Population III (Pop III) stars are the first generation of stars and because of the lack of metals in the primordial gas out of which they form, a small fraction of these stars may be supermassive stars (SMS). Some of these Pop III SMSs may undergo a general relativistic instability supernova (GRSN)\citep{chen2014,nagele2020}, and recently \citet{moriya2021} showed that the GRSN should be visible to JWST \citep{gardner2006}, even at high redshift. Pop III SMSs may also produce gamma ray bursts (GRBs) \citep{sun2017} and gravitational waves \citep{shibata2016,li2018} visible to future detectors such as THESEUS \citep{amati2018}, LISA \citep{amaro-Seoane2017}, and DECIGO \citep{kawamura2011}.



Cosmological simulations of the early universe often find massive hot gas clouds, which may, under special circumstances \citep{shang2010,inayoshi2015,hirano2017,regan2017}, form a single SMS with total mass up to $10^6\;{\rm M}_\odot$ \citep[e.g.][]{latif2013}. Two subsequent scenarios are studied, the first where all of the accretion occurs early during the protostellar phase, that is before hydrogen burning, and thus referred to as the non accreting scenario \citep{fuller1986,montero2012,chen2014,nagele2020}. The second, where the accretion occurs at a constant rate, is known as the accreting scenario \citep{hosokawa2012,hosokawa2013,schleicher2013,umeda2016,woods2017,hammerle2018a}. It should be noted that some combination of these two scenarios is also possible. One subtlety is that if the early accretion occurs too rapidly, hydrogen burning cannot ignite and the star will become rotationally supported, eventually collapsing to a black hole \citep{shibata2002}. Even if the star does not become rotationally supported, it is thought to rotate near the $\Omega-\Gamma$ limit \citep{hammerle2018}. 


Considering the non accreting scenario, \citet{chen2014} discovered that for a small mass range around 55000 ${\rm M}_\odot$ the general relativstic instability \citep{chandrasekhar1964} occurs when the SMS has a large reserve of helium, and the subsequent collapse triggers explosive helium burning which disrupts the star in a GRSN. Previously, we investigated this phenomenon using a post Newtonian stellar evolution code (as in \citealt{chen2014}) followed by a general relativistic hydrodynamical code adapted from pair instability and core collapse supernovae simulations \citep{nagele2020}. We found slightly different results from \citet{chen2014}, but did confirm that a GRSN was possible. In \citet{nagele2020}, our main difficulty had been determining when to connect the two codes. In this paper, we perform a general relativistic stability analysis for a better determination, and using this we find a significantly altered mass range for the GRSN.

This paper is organized as follows. In Sec. \ref{methods_num} we discuss the stellar and hydrodynamical codes as well as numerical models. In Sec. \ref{methods_stab} we outline the stability analysis. In Sec. \ref{results_comp} we compare the results of the stability analysis to other methods. In Sec. \ref{results_explo} we present simulations of GRSNe using the results of the stability analysis and compare the resulting ejecta to observations of metal poor stars. In Sec. \ref{results_col}, we discuss how this change might effect the neutrino emission of collapsing SMSs. Finally, we conclude with a discussion in Sec. \ref{discussion}.





\begin{figure}
    \centering
    \includegraphics[width=\columnwidth]{figures/schem2.pdf}

    \caption{Schematic showing the mechanism of the explosion. The horizontal axis is time since hydrogen ignition (not to scale) while the vertical axis is mass (also shown by the size of the stars; not to scale). If the GR instability (grey arrows) triggers while the star is in the helium burning phase, it will explode or pulsate. On the other hand if the star is in hydrogen burning or further developed (C+ stands for carbon and heavier), rapid helium burning does not occur and a black hole is formed.}
    \label{fig:schem}
\end{figure}


\section{Methods}
\label{methods}

The GRSN occurs when a non accreting SMS experiences the general relativstic (GR) instability during helium burning (Fig. \ref{fig:schem}). The star then contracts before rapidly burning a fraction of it's helium and then exploding. In order to properly model the phenomenon, three elements are required, a stellar evolution code, a determination of when the GR instability occurs, and a dynamical code for the explosion.  

\subsection{Numerical Models and Codes}
\label{methods_num}

The post Newtonian (PN) stellar evolution code (HOSHI) and the GR hydrodynamics code (HYDnuc) used in this paper are largely the same as in our previous works \citep{yamada1997,sumiyoshi2005,takahashi2016,takahashi2018,takahashi2019,yoshida2019,nagele2020,nagele2021}. Both are 1D hydrodynamics codes which include nuclear burning and neutrino cooling. In HOSHI, we use a 52 isotope nuclear network, which is the 49 isotope network from our previous paper with the addition of $^{14}$O, $^{18}$Ne, and $^{19}$Ne. These isotopes are required for the hot CNO cycle \citep{fuller1986} which is an ingredient of high mass SMS evolution.


HOSHI uses the first order PN approximation to the Tolman Oppenheimer Volkoff (TOV) equation. However, unlike in \citet{chen2014,nagele2020}, we include the correction of energy to density:
\begin{equation}
    \rho = \rho_0 \bigg(1+\frac{\epsilon}{c^2}\bigg)
    \label{eq:epsilon}
\end{equation}
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 it's inclusion is necessary to correctly model the SMS envelope. Fig. \ref{fig:HOSHI_r} shows the results of the stellar evolution calculation compared to \citet{nagele2020}.   

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 (e.g. Fig. \ref{fig:schem}). 

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. \ref{fig:exp_nconv} shows the explosion energy as a function of mesh number, total isotope number, and $\mathcal{V}$, the limit on the maximum variation of the independent variables 
\begin{equation}
    \mathcal{V} > \rm{max}_{i,j} \; \Bigg|  \frac{x_{i}(j,k)}{x_{i}(j,k+1)} \Bigg|^{\pm 1}
    \label{eq_varlim}
\end{equation}
where $x_i$ is one of the independent variables of HYDnuc \citep{yamada1997}, j is the mesh number, and k is the time-step. In this paper, we use $\mathcal{V}=10^{-5}$, a 61 isotope network, and 767 meshes. 


The explosion energy depends on $\mathcal{V}$ and mesh 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 $^{20}$Ne and $^{24}$Mg; we use 1.5 $ \times $ \citet{caughlan1988}, 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 burnt via carbon alpha capture on the timescale of the explosion. However, if nucleons are present, catalysis enhances the carbon alpha capture rate with 
\begin{equation}
    \rm ^{12}C(p,\gamma) ^{13}N, \;\; ^{13}N(\alpha,p) ^{16}O.
\end{equation}


In the networks with isotope number less than 61 (Fig. \ref{fig:cx_t}, 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 
\begin{equation}
    \rm ^{24}Mg(p,\gamma) ^{25}Al.
\end{equation}
Indeed, aluminium is of particular importance (Fig. \ref{fig:cx_t}, right panel) because the inclusion of it's isotopes is the only difference between the 58 isotope network and the 61 isotope network, and from Fig. \ref{fig:exp_nconv}, 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. \citet{nagele2020} used a 49 isotope network and \citet{chen2014} 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. If we were to use a less accurate, smaller network, we would find more exploding models, some of which would have explosion energies comparable to those in previous works. Thus the network size not only determines the explosion energy, but also the existence of an explosion and so must be treated with great care. 


\subsection{GR Stability Analysis}
\label{methods_stab}

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 \citet{nagele2020}, we relied on the PN stellar evolution code to determine when collapse would occur, and our results were in rough agreement with those of \citet{chen2014}. In this work we perform a stability analysis on the normal modes of radial perturbations of a star in GR \citep{chandrasekhar1964}. 

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 \citep{shapiro1983}

        \begin{equation}
            \label{eq_N}
            \dv{}{r} \; \bigg[ \Gamma_1 \frac{P}{r^2} \dv{}{r} \; (r^2\xi)\bigg] - \frac{4}{r}\dv{P}{r} \xi +\omega^2\rho_0\xi = 0
        \end{equation}
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):
\begin{equation}
    \Gamma_1 = \pdv{\ln P}{\ln \rho_0}  \bigg\rvert_s.
\end{equation}


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 
\begin{equation}
    \omega^2_0<\omega^2_1<\omega^2_2<...
\end{equation}
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

\begin{equation}
    \omega^2_0 < 0.
\end{equation}
The corresponding equation in GR is \citep{chandrasekhar1964}
        \begin{equation}
            \nonumber
            e^{-2a-b}\dv{}{r} \; \Big[ e^{3a+b}\Gamma_1 \frac{P}{r^2} \dv{}{r} \; (e^{-a}r^2\xi)\Big] - \frac{4}{r}\dv{P}{r} \xi +e^{-2a-2b}\omega^2(P+\rho c^2)\xi 
        \end{equation}     
        \begin{equation}
            \label{eq_GR}
            - \frac{8\pi G}{c^4}e^{2b}P(P+\rho c^2)\xi - \frac{1}{P+\rho c^2} \bigg(\dv{P}{r}\bigg)^2 \xi = 0
        \end{equation}
where $a,b$ are the metric coefficients as in \citet{haemmerle2020} and the density is defined in Eq. $\ref{eq:epsilon}$. 


In order to solve Eqs. \ref{eq_N},\ref{eq_GR} for $\omega^2_0$, we adopt an iterative method similar to that outlined in exercise 6.11 of \citet{shapiro1983}. For a detailed discussion, see Appendix A. 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 \citep{chandrasekhar1964}. Fig. \ref{fig:nconv_inst} 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. \ref{fig:insta_hoshi}), the calculation is mapped to the hydrodynamical code. 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. \ref{fig:insta_hydro} 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. \ref{fig:insta_hydro}), 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 burnt most of its helium (Table \ref{tab:models}).








\begin{figure}
    \centering
    \includegraphics[width=\columnwidth]{figures/figk.pdf}

    \caption{Comparison of radial time snapshots for the 2 $\times 10^4$ ${\rm M}_\odot$ model in the HOSHI code used in this work (dashed lines) and the one used in \citet{nagele2020} (dotted lines). The inclusion of internal energy to the GR density when evaluating the PN TOV equation causes a more compact inner envelope, while the outer envelope remains largely unchanged.}
    \label{fig:HOSHI_r}
\end{figure}


\begin{figure*}
    \centering
    \includegraphics[width=2\columnwidth]{figures/figp.pdf}

    \caption{Numerical convergence of explosion energy for three parameters of HYDnuc for the 3 $\times 10^4$ ${\rm M}_\odot$ model. Left panel --- dependence on $\mathcal{V}$ (Eq. \ref{eq_varlim}) with total mesh number = 255 and isotope number = 52. Middle panel ---  dependence on number of isotopes in the nuclear network (52, 58, 61, 79, 89, 153, 300; see Appendix B, Table \ref{tab:networks}) with total mesh number = 255 and $\mathcal{V} =5\times 10^{-4}$. Right panel ---  dependence of explosion energy on number of mesh grids, with $\mathcal{V} = 5\times 10^{-4}$ and isotope number = 52.}
    \label{fig:exp_nconv}
\end{figure*}




\begin{figure}
    \centering
    \includegraphics[width=\columnwidth]{figures/figj.pdf}

    \caption{Numerical convergence of the stability analysis code for numerical polytropes with spacing d$\xi$. These polytropes are evaluated at $\Gamma_1=4/3$ for the Newtonian analysis and $\Gamma_1=4/3+\kappa GM/(rc^2) $ for the GR analysis so that we should find $\omega_0^2 = 0$. Thus, the y axis can be regarded as the error associated with our method, and its amplitude decreases for increasing numerical resolution. The number of meshes in HOSHI is variable, but is usually between 1000 and 2000, translating to somewhere on the left side of this plot. }
    \label{fig:nconv_inst}
\end{figure}



\begin{figure}
    \centering
    \includegraphics[width=\columnwidth]{figures/figo.pdf}

    \caption{GR stability analysis (Sec. \ref{methods_stab}) applied to the results of the 3 $\times 10^4$ ${\rm M}_\odot$ model in HOSHI. Upper panel --- amplitude of the fundamental mode frequency as a function of time for the General Relativistic (blue) and Newtonian (red) analyses. Unstable models, having a negative value of frequency squared, are denoted by os while stable models are denoted by crosses. Lower panel --- Central mass fraction of various isotopes as a function of time. The first instability occurs during helium burning. }
    \label{fig:insta_hoshi}
\end{figure}





\begin{figure}
    \centering
    \includegraphics[width=\columnwidth]{figures/figl.pdf}

    \caption{Stability analysis performed on the results of the HYDnuc calculation for the first unstable model of 2 $\times 10^4$ ${\rm M}_\odot$. Upper panel --- stability analysis. Lower panel --- central temperature. }
    \label{fig:insta_hydro}
\end{figure}





\begin{figure}
    \centering
    \includegraphics[width=\columnwidth]{figures/figm.pdf}

    \caption{Illustration of the polytropic criterion for the models in this paper. Upper panel --- mass of the helium core from HOSHI, determined by the mass outside which X($^1$H) $>$ 1e-5. Middle panel --- critical density (Sec. \ref{results_comp}) of a pure helium core with the mass from the top panel. Lower panel --- central density of the HOSHI models at the maximum central helium mass fraction. A comparison of the middle and lower panels shows that these models are not yet unstable according to the polytropic criterion.}
    \label{fig:polytropic}
\end{figure}

\begin{figure}
    \centering
    \includegraphics[width=\columnwidth]{figures/fign.pdf}

    \caption{Comparison of the first instability reached by two methods, Sec. \ref{methods_stab} and that of \citet{haemmerle2020}. Note that the values in this figure differ from Table \ref{tab:models} for lower masses because the first instability does not always collapse in HYDnuc (Sec. \ref{methods_stab}).} 
    \label{fig:haemmerle}
\end{figure}



\begin{figure*}
    \centering
    \includegraphics[width=2\columnwidth]{figures/fige.pdf}

    \caption{Left panel --- local energy as a function of mass coordinate for the 2.6 $\times 10^4$ ${\rm M}_\odot$ model (pulsation). The ejection criterion is $e_{\rm tot}>0$ for all $m_r > \rm M - M_{\rm ej}$ at the time of shock breakout. Right panel --- illustration of the stability of the ejection criterion. We measure the ejected mass at shock breakout and this quantity is roughly consistent with the post Newtonian escape velocity criterion near the end of the simulation. .}
    \label{fig:E_pulse}
\end{figure*}




\begin{figure*}
    \centering
    \includegraphics[width=2\columnwidth]{figures/figa.pdf}

    \caption{Time evolution of central temperature (1st panel), density (2nd panel), rate of change of specific internal energy (3rd panel) and entropy relative to the initial value (4th panel). The legend groups models by outcome, whereas colors vary with mass. The left column shows the exploding and pulsating models as a function of time. The right column shows the temperature as a function of time, while the other three panels are functions of temperature. }
    \label{fig:central}
\end{figure*}


\begin{figure*}
    \centering
    \includegraphics[width=2\columnwidth]{figures/figf.pdf}

    \caption{Upper panel --- comparison of central isotope mass fractions for the M = 3 $\times 10^4$ ${\rm M}_\odot$ model with a 58 isotope model (left) and a 61 isotope model (right). The high mass fraction of protons and $^{13}$N in the left panel facilitates artificially high carbon burning. Lower panel ---  time evolution of central $\dot{\epsilon}$.}
    \label{fig:cx_t}
\end{figure*}


\begin{figure}
    \centering
    \includegraphics[width=\columnwidth]{figures/figb.pdf}

    \caption{Total energy, kinetic energy, gravitational energy, and thermal energy (as defined in \citet{nagele2020}) in HYDnuc as a function of time for the explosion of the 3 $\times 10^4$ ${\rm M}_\odot$ model. Unless otherwise specified, subsequent figure will show the 3 $\times 10^4$ ${\rm M}_\odot$ model.}
    \label{fig:E_explo}
\end{figure}





\begin{figure*}
    \centering
    \includegraphics[width=2\columnwidth]{figures/figg.pdf}

    \caption{Initial (upper panel) and final (lower panel) isotope mass fractions as a function of the mass coordinate for the pulsation, M = 2.6 $\times 10^4$ ${\rm M}_\odot$ (left) and the explosion, M = 3 $\times 10^4$ ${\rm M}_\odot$ (right). }
    \label{fig:cx_mr}
\end{figure*}

\begin{figure}
    \centering
    \includegraphics[width=\columnwidth]{figures/figc.pdf}

    \caption{Log velocity as a function of time and mass coordinate. Shock breakout occurs about $10^5$ seconds after peak temperature.}
    \label{fig:v}
\end{figure}

\begin{figure}
    \centering
    \includegraphics[width=\columnwidth]{figures/figd.pdf}

    \caption{Log radius as a function of time and mass coordinate. }
    \label{fig:r_mr}
\end{figure}





\begin{figure}
    \centering
    \includegraphics[width=1\columnwidth]{figures/figr.pdf}

    \caption{Explosion yield (lines) compared with observed metal poor stars (crosses are observations, triangles upper limits). The line shows the fraction relative to a minimum mass of hydrogen with which it would have to mix \citet{magg2020}, and can thus be regarded as an upper limit. }
    \label{fig:XH}
\end{figure}





\begin{table*}
	\centering
	\caption{Summary table for all models. The columns are total mass, outcome of HYDnuc, mass of the isentropic core, central helium mass fraction at the start of HYDnuc, change in helium mass fraction, explosion energy, maximum central temperature, and maximum velocity of the outermost mesh, denoted $v_R$.}
	\label{tab:models}
	\begin{tabular}{|c|l|l|l|l|l|l|l|} 
		
    		M [$10^4$ ${\rm M}_\odot$] & Outcome & $M_{\rm core}$ [${\rm M}_\odot$] & $X_{\rm c}(^4$He) & $\Delta \; X_{\rm c}(^4$He) & $E_{\rm exp}$ [ergs] & max $T_{ \rm c} $ [K]& max $v_R$/c  \\
    		\hline
2 & Collapse & 10926 & 1.37e-3 &  ---  &  ---  &  ---  &  ---  \\ 
2.1 & Collapse & 11368 & 2.23e-4 &  ---  &  ---  &  ---  &  ---  \\ 
2.2 & Collapse & 11729 & 1.22e-4 &  ---  &  ---  &  ---  &  ---  \\ 
2.3 & Collapse & 12595 & 3.17e-2 &  ---  &  ---  &  ---  &  ---  \\ 
2.4 & Collapse & 13180 & 3.44e-18 &  ---  &  ---  &  ---  &  ---  \\ 
2.5 & Collapse & 13798 & 2.69e-3 &  ---  &  ---  &  ---  &  ---  \\ 
2.6 & Pulsation & 14772 & 0.104 & 0.104 & 4.32e53 & 7.58e8 & 0.032 \\ 
2.7 & Pulsation & 14964 & 0.222 & 0.147 & 4.70e52 & 6.62e8 & 0.021 \\ 
2.8 & Collapse & 15596 & 0.713 &  ---  &  ---  &  ---  &  ---  \\ 
2.9 & Pulsation & 16183 & 0.589 & 0.153 & 7.56e53 & 7.33e8 & 0.041 \\ 
2.95 & Explosion & 16504 & 0.599 & 0.168 & 1.23e54 & 7.69e8 & 0.046 \\ 
3 & Explosion & 16817 & 0.652 & 0.152 & 1.43e54 & 8.06e8 & 0.048 \\ 
3.05 & Collapse & 17144 & 0.734 &  ---  &  ---  &  ---  &  ---  \\ 
3.1 & Collapse & 17516 & 0.794 &  ---  &  ---  &  ---  &  ---  \\ 
3.15 & Collapse & 17793 & 0.815 &  ---  &  ---  &  ---  &  ---  \\ 
3.2 & Collapse & 18091 & 0.815 &  ---  &  ---  &  ---  &  ---  \\ 
3.3 & Collapse & 18888 & 1.000 &  ---  &  ---  &  ---  &  ---  \\ 
3.4 & Collapse & 19460 & 1.000 &  ---  &  ---  &  ---  &  ---  \\ 
3.5 & Collapse & 19933 & 0.950 &  ---  &  ---  &  ---  &  ---  \\ 
4 & Collapse & 23891 & 0.960 &  ---  &  ---  &  ---  &  ---  \\ 

		\hline
	\end{tabular}
\end{table*}

\begin{table*}
	\centering
	\caption{Mass ejecta by isotope for the explosions and the pulsations. Except for the first column which is consistent with Table \ref{tab:models}, values are recorded in units of ${\rm M}_\odot$.}
	\label{tab:ejecta}
	\begin{tabular}{|c|l|l|l|l|l|l|l|l|l|} 
		
    		M [$10^4$ ${\rm M}_\odot$] &M$_{\rm ej}$ & M($^1$H)  & M($^4$He) & M($^{12}$C) & M($^{16}$O) & M($^{20}$Ne) & M($^{24}$Mg) & M($^{28}$Si) & M($^{32}$S)  \\
    		\hline
2.6 & 2808 & 1877 & 974 & < 0.1 & < 0.1 & < 0.1 & < 0.1 & < 0.1 & < 0.1 \\
2.7 & 2299 & 1584 & 759 & < 0.1 & < 0.1 & < 0.1 & < 0.1 & < 0.1 & < 0.1 \\
2.9 & 2078 & 1465 & 651 & < 0.1 & < 0.1 & < 0.1 & < 0.1 & < 0.1 & < 0.1 \\
2.95 & 29500 & 5441 & 16946 & 3006 & 2505 & 481 & 812 & 306 & 1.2 \\
3 & 30000 & 5537 & 18077 & 2986 & 1829 & 367 & 702 & 497 & 5.1 \\

		\hline
	\end{tabular}
\end{table*}





\section{Results}
\label{results}

The primary advantage of this paper is the GR stability analysis, so before discussing applications, we will compare this analysis to two previous methods. 

\subsection{Comparison to previous works}
\label{results_comp}


First, and most common in the literature \citep[e.g.][]{fuller1986,umeda2016,woods2017} 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}$ \citep{shapiro1983}. From here, the polytropic constant is determined
\begin{equation}
    K=\frac{a}{3} \bigg( \frac{3s_r}{4m_pa} \bigg)^{4/3}.
\end{equation}
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$, 
\begin{equation}
     \frac{4}{3} + \frac{\beta}{6} + \mathcal{O}(\beta^2) = \frac{4}{3} + \frac{2GM_{\rm He}\kappa}{R_{\rm crit} c^2}
    \label{eq:rcrit}
\end{equation}
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 \citet{chandrasekhar1964} by a factor of 2). Eq. \ref{eq:rcrit} 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 \citep{shapiro1983}
\begin{equation}
    \rho_{\rm crit} = \Bigg[\frac{R_{\rm crit}}{\xi_1} \bigg( \frac{K}{\pi G}\bigg)^{-1/2} \Bigg]^{-3} \propto M_{\rm He} ^ {-7/2}
\end{equation}
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. \ref{methods_stab}. Fig. \ref{fig:polytropic} 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. \ref{methods_stab}.

Next we compare our method to the results of \citet{haemmerle2020}, who also evaluate Eq. \ref{eq_GR}, 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. \ref{fig:haemmerle} 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 \ref{tab:models}. This would mean that the explosion would occur at larger mass, if we were to use the method of \citet{haemmerle2020}, and this is because their method is merely a sufficient condition for instability. 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 \citet{chen2014,nagele2020}.



\subsection{Application to GRSNe}
\label{results_explo}

\subsubsection{HYDnuc explosion}

As hinted at in previous sections, we find two GRSNe as well as several pulsations (Table \ref{tab:models}). The pulsations are less energetic events which eject a portion of the envelope (Table \ref{tab:ejecta}). The relationship between the GRSNe and these pulsations is likely similar to the relationship between pair instability supernovae and the pulsational pair instability process \citep{woosley2017}. Previous papers did not find any pulsations because the models with mass lower than the explosion mass collapsed from entering the pair unstable region, as opposed to being triggered by the GR instability. Thus, the discovery of pulsations is a straightforward consequence of using the more accurate GR instability analysis. For the pulsations, we determine the ejecta mass in Table \ref{tab:models} using the local energy, e(r) (Fig. \ref{fig:E_pulse}, left panel), which is the integrand of the global energies defined in \citet{nagele2020}. We measure this quantity at shock breakout (right panel), when some of the energy is still in the form of thermal energy. The energy evolution after shock breakout may not be completely accurate. We confirm that the escape velocity criterion (right panel) converges roughly to this same value, validating the use of this ejection criterion. Note that $2.9 \times 10^4$ ${\rm M}_\odot$ is a marginal case. Although it is not an explosion, it does not have the steady behaviour of Fig. \ref{fig:E_pulse} after shock breakout, and we expect the value in Table \ref{tab:ejecta} to underestimate the ejecta mass for this case. 


The mass range for the explosion follows from straightforward considerations. Models which experience the instability before helium burning has reduced the binding energy of the star cannot explode. On the other hand, models with helium mass fractions less than ten percent also do not explode or pulsate because they lack fuel for alpha capture reactions sufficient to halt the collapse.

The pulsating models all have lower mass than the exploding models and the most massive model being the most energetic before a sharp drop-off is reminiscent of Fig. 6 of \citet{nagele2020}. This is a likely characteristic of the GRSNe even if the mass range found by the current analysis is not completely correct.


Fig. \ref{fig:central} shows the behavior of the central temperature, baryonic density, nuclear energy generation rate, and entropy as a function of time. For the exploding and pulsating models (left column), temperature and density increase and decrease smoothly, as was the case for the GRSNe in \citet{chen2014,nagele2020}, and the nuclear energy generation rate is also smooth. The models which collapse (right column) also show relatively smooth increases in temperature and density, but the energy generation rate and change in entropy are more complicated, as the stars enter different phases of burning and photodissociation. Roughly speaking, the first three peaks in $\dot{\epsilon}_c$ (first peak at $10^9$ K) correspond to carbon burning, sulfur burning, and calcium to iron peak element burning. Stars with a more evolved core, such as 2.4 $\times $ $10^4$ ${\rm M}_\odot$ do not have much carbon remaining, and the first peak corresponds to neon burning. At the other extreme, stars with reserves of hydrogen, such as 4 $\times $ $10^4$ ${\rm M}_\odot$ do not exhibit peaks because the presence of protons increases the number of possible reactions. After the core becomes iron/nickel, the next major reactions are assisted by free nucleons created by photodissociation (one peak), which is finally followed by photodissociation of nickel (one peak) and photodissociation of helium (two peaks). In reality, neutrino reaction would then begin to dominate, but HYDnuc does not include most of the relevant neutrino reactions. 




For the explosion of 3 $\times 10^4$ ${\rm M}_\odot$--- which we will use as an example in the figures--- Fig. \ref{fig:E_explo} shows the time evolution of the energy quantities, including the total energy which eventually determines the explosion energy (Table \ref{tab:models}). There are three phases, the initial contracting phases with $E_{\rm tot}<0$, the pre-shock breakout phase with $E_{\rm tot}>E_{\rm kin}>0$, and the post-shock breakout phase with $E_{\rm tot}=E_{\rm kin}>0$. Fig. \ref{fig:cx_mr} shows the isotope mass fraction as a function of mass coordinate for the initial (upper) and final (lower) time steps. The SMS cores are initially isentropic and shell hydrogen burning often occurs. The envelope is divided into many convective layers, the exact layout of which can alter the stability of the star \citep{nagele2020}. The final isotope distributions show that the majority of the nuclear burning takes place within the inner 5000 ${\rm M}_\odot$ for the pulsating model and the inner 10000 ${\rm M}_\odot$ for the exploding model. For the exploding model, the star is totally disrupted (e.g. Fig. \ref{fig:v}), so these elements will be ejected into the inter stellar medium (ISM) (Table \ref{tab:ejecta}).


After the explosive nuclear burning, the inwards velocity rapidly reverses and the shock propagates towards the surface of the SMS with a typical velocity of a few percent the speed of light (Table \ref{tab:models}); shock breakout occurs on a timescale of $10^5$ s (Fig. \ref{fig:v}). The initial inwards velocity is largest in the envelope, and we emphasize that the GRSN involves the collapse of the entire star, not just the core. This is another reason why the analysis in Sec. \ref{methods_stab} is necessary, as an analysis of the core alone will not always capture the instability. Fig. \ref{fig:r_mr} shows the time evolution of radius as a function of mass coordinate. The path of the shock can be seen in the upper part of the figure.  

\subsubsection{Comparison to observed metal poor stars}

We find a wider explosion window than previous works, as well as pulsations which could later explode. However, the explosion window is still narrow compared to the feasible mass range of SMSs. This narrowness, in combination with the probable paucity of SMSs themselves, means that the likelihood of a GRSN having occurred in the Milky Way volume ($r = 1.4 \pm 0.4$ Mpc at z = 127 \citealt{griffen2016}) is low. With that disclaimer in place, however, we will now compare the GRSN ejecta to observed metal poor stars. 

In this section, we will consider only the mono-enrichment scenario (see e.g. \citealt{hartwig2018a,hartwig2019}), so that the metals in the metal poor star have come exclusively from the GRSN ejecta and thus we ignore potential pollution from ISM accretion. we determine the minimum mass of ISM material with which the ejecta could mix via Eq. 2 of \citet{magg2020} where the explosion energy is taken from Table \ref{tab:models} and we use the fiducial value of the number density from \citet{magg2020}, $n_0 = 1$ cm$^{-3}$. We combine this mass with our ejecta and plot the inferred abundances relative to hydrogen, relative to solar (Fig. \ref{fig:XH}). Since the mixing mass is a lower limit, the abundance ratios are an upper limit. 

Also shown are the inferred abundances of the twenty stars in the SAGA database \citep{suda2008} with the lowest values of Fe / H: \citet{keller2014,ezzeddine2019,aguado2018,aoki2006,frebel2008,christlieb2004,bonifacio2018,bonifacio2015,frebel2015,caffau2016,norris2007,caffau2011,hansen2014,hansen2015,starkenburg2018,roederer2014a,roederer2014b}. We select these stars because the GRSN produces negligible amounts of iron. Fig. \ref{fig:XH} shows that our inferred yields do not match any of the observed stars. Even for the star of \citet{keller2014} which has strict upper limits on [Fe/H], our yield misses the observed value of calcium by several orders of magnitude. 

While we can definitively rule out the mono-enrichment scenario for metal poor stars in the vicinity of our galaxy, the multi-enrichment scenario is more challenging to completely rule out. The best we can do at current is to note that although many metal poor stars are carbon enriched, they are not generally silicon and magnesium enriched, which is strong evidence against an SMS being involved in the multi-enrichment scenario. Thus we can say that it is likely the case that there did not exist a non-accreting SMS in the explosion mass range, in the Milky Way volume.


\subsection{Application to SMS Collapse and neutrino emission}
\label{results_col}

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 \citep{nagele2021}. 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 \citet{nagele2021}, 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. 



\section{Discussion}
\label{discussion}

Using the general relativistic stability analysis in Sec. \ref{methods_stab}, we can more accurately predict when a SMS will be unstable and will collapse explode or pulse. This is necessary because the timescale of the collapse is many orders of magnitude shorter than the evolutionary timescale, meaning that our stellar evolution code often misses the GR instabilities. For some models, the instability is countered by increased nuclear burning, but for masses around 3 $\times 10^4$ ${\rm M}_\odot$, this is not the case. The 2.95 $\times 10^4$ ${\rm M}_\odot$ and 3 $\times 10^4$ ${\rm M}_\odot$ models explode in GRSNe, while several lower mass models pulsate and eject portions of their envelope. The final fate of the pulsating models is unclear, as they will reenter the evolutionary track with different properties, most notably the chemical composition of the core and total energy. If multiple pulsations were to occur, it could cause a collisional supernova. 

In comparison to the GRSNe from \citet{chen2014,nagele2020}, our GRSNe have a much lower explosion energy, though if we were to use an unphysically small nuclear network, we would find a GRSN with comparably large energies. The composition of our explosion ejecta is similar to the previous works (Table \ref{tab:ejecta}). On the other hand, the mass range of the explosion is wider because our more accurate determination of the GR instability means that these models are not near the pair instability region. 

\citet{moriya2021} found that the observational duration of a GRSN will be on the order of $10^{2-3}$ seconds for the peak emission and $10^{2-3}$ days for the plateau, which shares some similarities with Type IIP SNe. Because the GRSN would have to occur in the high redshift universe, this means that the observer duration is longer by about a factor of ten. \citet{moriya2021} demonstrated that the GRSN plateau may be differentiated from other persistent sources if it is observed in multiple bands. In the future, we plan to investigate if the lower energy of our GRSN would signifigantly alter any of their findings.

We determined that none of the observed metal poor stars match the inferred abundance pattern from the GRSNe, and from this we conclude that none of these stars were singly enriched by a GRSN. However, we note that the multi-enrichment scenario cannot be completley ruled out.

As in all numerical studies, there are numerous uncertainties. In HOSHI, along with the usual sources of numerical error, we also need to consider the error due to using the post Newtonian approximation. After the inclusion of internal energy to density (Eq. \ref{eq:epsilon}) in the TOV equation, the post Newtonian pressure gradient matches the TOV pressure gradient to one part in $10^5$. Although this level of accuracy likely supersedes numerical error, we have to keep in mind that the TOV equation assumes a hydrostatic configuration, which, for instance, is not true as the star contracts towards the end of hydrogen burning. Finally, and perhaps most importantly, our evolutionary models are not rotating, whereas real SMSs are expected to initially be medium rotators \citep{hammerle2018} which may spin up over the course of their lifetimes \citep{maeder2001}. We intend to more fully investigate the effects of rotation in the future.  

Regarding the GR stability analysis, we were able to quantify the error of $|\omega_0^2|$ for polytropes (Fig. \ref{fig:nconv_inst}) as being around $10^{-7,-8}$. As seen in Fig. \ref{fig:insta_hoshi}, typical values of $\omega_0^2$ are greater than this error. Although we demonstrated that the error decreases for increasing resolution (Fig. \ref{fig:nconv_inst}) the gain in accuracy is low compared to the gain in computational time that an increased mesh number of an order of magnitude or two would require. 

The largest source of uncertainty in HYDnuc may be the coordinate perturbation at the start of the calculation. We sought to minimize this perturbation via Eq. \ref{eq:epsilon}, but the possibility remains that the perturbation size may affect the explosion energy. Indeed in Fig. \ref{fig:exp_nconv} left panel we show how sensitive the explosion energy can be to the collapse timescale, which would vary with a different initial perturbation. That being said, a change in the initial perturbation would likely only shift the mass range of the explosion. As far as other sources of error go, \citet{chen2014} showed that multidimensional effects may not play a big role in the GRSN, and they are also not thought to contribute to the instability \citep{chandrasekhar1965}. We verified that radiative and convective energy transport do not effect the explosion outcome and tested that the explosion energy does not depend on the choice of the carbon alpha capture rate. 


In the future, we plan to apply the GR stability analysis to rotating SMSs, as well as verifying our current results with multidimensional simulations. We also intend to asses the observability of the GRSN using methods similar to those in \citet{moriya2021}. Finally, we will investigate the possibility of multiple pulsations. 



\section*{Data Availability}

The data underlying this article will be shared on reasonable request to the corresponding author.

\section*{Acknowledgements}

This study was supported in part by the Grant-in-Aid for the Scientific Research of Japan Society for the Promotion of Science (JSPS, Nos. JP19K03837, JP20H01905, JP20H00158, JP21H01123) and by Grant-in-Aid for Scientific Research on Innovative areas (JP17H06357, JP17H06365) from the Ministry of Education, Culture, Sports, Science and Technology (MEXT), Japan. For providing high performance computing resources, YITP, Kyoto University is acknowledged. K.S. would also like to acknowledge computing resources at KEK and RCNP Osaka University. C.N. would like to thank Tilman Hartwig for discussions of metal poor stars. 











\section{Introduction}
\label{introduction}


For much of the history of astronomy, the post recombination early universe has been inaccessible to observation. While this state of affairs still holds broadly, inroads are being made. Observers have detected a star at redshift 6 \citep{welch2022}, galaxies as early as redshift 11 \citep{oesch2016}, a long gamma ray burst (GRB) at redshift 8 \citep{tanvir2009}, and quasars at redshifts 6 and 7 \citep{mortlock2011,wu2015,banados2018,matsuoka2019,wang2021}. These high redshift observations have greatly increased our knowledge of the early universe, but they also raise questions--- most notably, where did the high redshift quasars and their supermassive black hole (SMBH) engines come from?



Several theories have been put forward to explain the existence of SMBHs so soon after the big bang \citep[e.g.][]{rees1984,inayoshi2020}. These include, but are not limited to, the direct collapse scenario \citep{bromm2003}, super-Eddington accretion onto solar mass black holes \citep{haiman2001}, and rapid mergers of either primordial \citep{bean2002} or astrophysical black holes \citep{omukai2008}. Of these scenarios, the direct collapse scenario may be the easiest to test observationally; Population III (Pop III) stars are the first generation of stars and because of the lack of metals in the primordial gas out of which they form, a small fraction of these stars may be supermassive stars (SMS). Some of these Pop III SMSs may undergo a general relativistic instability supernova (GRSN)\citep{chen2014,nagele2020}, and recently \citet{moriya2021} showed that the GRSN should be visible to JWST \citep{gardner2006}, even at high redshift. Pop III SMSs may also produce gamma ray bursts (GRBs) \citep{sun2017} and gravitational waves \citep{shibata2016,li2018} visible to future detectors such as THESEUS \citep{amati2018}, LISA \citep{amaro-Seoane2017}, and DECIGO \citep{kawamura2011}.



Cosmological simulations of the early universe often find massive hot gas clouds, which may, under special circumstances \citep{shang2010,inayoshi2015,hirano2017,regan2017}, form a single SMS with total mass up to $10^6\;{\rm M}_\odot$ \citep[e.g.][]{latif2013}. Two subsequent scenarios are studied, the first where all of the accretion occurs early during the protostellar phase, that is before hydrogen burning, and thus referred to as the non accreting scenario \citep{fuller1986,montero2012,chen2014,nagele2020}. The second, where the accretion occurs at a constant rate, is known as the accreting scenario \citep{hosokawa2012,hosokawa2013,schleicher2013,umeda2016,woods2017,hammerle2018a}. It should be noted that some combination of these two scenarios is also possible. One subtlety is that if the early accretion occurs too rapidly, hydrogen burning cannot ignite and the star will become rotationally supported, eventually collapsing to a black hole \citep{shibata2002}. Even if the star does not become rotationally supported, it is thought to rotate near the $\Omega-\Gamma$ limit \citep{hammerle2018}. 


Considering the non accreting scenario, \citet{chen2014} discovered that for a small mass range around 55000 ${\rm M}_\odot$ the general relativstic instability \citep{chandrasekhar1964} occurs when the SMS has a large reserve of helium, and the subsequent collapse triggers explosive helium burning which disrupts the star in a GRSN. Previously, we investigated this phenomenon using a post Newtonian stellar evolution code (as in \citealt{chen2014}) followed by a general relativistic hydrodynamical code adapted from pair instability and core collapse supernovae simulations \citep{nagele2020}. We found slightly different results from \citet{chen2014}, but did confirm that a GRSN was possible. In \citet{nagele2020}, our main difficulty had been determining when to connect the two codes. In this paper, we perform a general relativistic stability analysis for a better determination, and using this we find a significantly altered mass range for the GRSN.

This paper is organized as follows. In Sec. \ref{methods_num} we discuss the stellar and hydrodynamical codes as well as numerical models. In Sec. \ref{methods_stab} we outline the stability analysis. In Sec. \ref{results_comp} we compare the results of the stability analysis to other methods. In Sec. \ref{results_explo} we present simulations of GRSNe using the results of the stability analysis and compare the resulting ejecta to observations of metal poor stars. In Sec. \ref{results_col}, we discuss how this change might effect the neutrino emission of collapsing SMSs. Finally, we conclude with a discussion in Sec. \ref{discussion}.





\begin{figure}
    \centering
    \includegraphics[width=\columnwidth]{figures/schem2.pdf}

    \caption{Schematic showing the mechanism of the explosion. The horizontal axis is time since hydrogen ignition (not to scale) while the vertical axis is mass (also shown by the size of the stars; not to scale). If the GR instability (grey arrows) triggers while the star is in the helium burning phase, it will explode or pulsate. On the other hand if the star is in hydrogen burning or further developed (C+ stands for carbon and heavier), rapid helium burning does not occur and a black hole is formed.}
    \label{fig:schem}
\end{figure}


\section{Methods}
\label{methods}

The GRSN occurs when a non accreting SMS experiences the general relativstic (GR) instability during helium burning (Fig. \ref{fig:schem}). The star then contracts before rapidly burning a fraction of it's helium and then exploding. In order to properly model the phenomenon, three elements are required, a stellar evolution code, a determination of when the GR instability occurs, and a dynamical code for the explosion.  

\subsection{Numerical Models and Codes}
\label{methods_num}

The post Newtonian (PN) stellar evolution code (HOSHI) and the GR hydrodynamics code (HYDnuc) used in this paper are largely the same as in our previous works \citep{yamada1997,sumiyoshi2005,takahashi2016,takahashi2018,takahashi2019,yoshida2019,nagele2020,nagele2021}. Both are 1D hydrodynamics codes which include nuclear burning and neutrino cooling. In HOSHI, we use a 52 isotope nuclear network, which is the 49 isotope network from our previous paper with the addition of $^{14}$O, $^{18}$Ne, and $^{19}$Ne. These isotopes are required for the hot CNO cycle \citep{fuller1986} which is an ingredient of high mass SMS evolution.


HOSHI uses the first order PN approximation to the Tolman Oppenheimer Volkoff (TOV) equation. However, unlike in \citet{chen2014,nagele2020}, we include the correction of energy to density:
\begin{equation}
    \rho = \rho_0 \bigg(1+\frac{\epsilon}{c^2}\bigg)
    \label{eq:epsilon}
\end{equation}
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 it's inclusion is necessary to correctly model the SMS envelope. Fig. \ref{fig:HOSHI_r} shows the results of the stellar evolution calculation compared to \citet{nagele2020}.   

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 (e.g. Fig. \ref{fig:schem}). 

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. \ref{fig:exp_nconv} shows the explosion energy as a function of mesh number, total isotope number, and $\mathcal{V}$, the limit on the maximum variation of the independent variables 
\begin{equation}
    \mathcal{V} > \rm{max}_{i,j} \; \Bigg|  \frac{x_{i}(j,k)}{x_{i}(j,k+1)} \Bigg|^{\pm 1}
    \label{eq_varlim}
\end{equation}
where $x_i$ is one of the independent variables of HYDnuc \citep{yamada1997}, j is the mesh number, and k is the time-step. In this paper, we use $\mathcal{V}=10^{-5}$, a 61 isotope network, and 767 meshes. 


The explosion energy depends on $\mathcal{V}$ and mesh 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 $^{20}$Ne and $^{24}$Mg; we use 1.5 $ \times $ \citet{caughlan1988}, 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 burnt via carbon alpha capture on the timescale of the explosion. However, if nucleons are present, catalysis enhances the carbon alpha capture rate with 
\begin{equation}
    \rm ^{12}C(p,\gamma) ^{13}N, \;\; ^{13}N(\alpha,p) ^{16}O.
\end{equation}


In the networks with isotope number less than 61 (Fig. \ref{fig:cx_t}, 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 
\begin{equation}
    \rm ^{24}Mg(p,\gamma) ^{25}Al.
\end{equation}
Indeed, aluminium is of particular importance (Fig. \ref{fig:cx_t}, right panel) because the inclusion of it's isotopes is the only difference between the 58 isotope network and the 61 isotope network, and from Fig. \ref{fig:exp_nconv}, 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. \citet{nagele2020} used a 49 isotope network and \citet{chen2014} 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. If we were to use a less accurate, smaller network, we would find more exploding models, some of which would have explosion energies comparable to those in previous works. Thus the network size not only determines the explosion energy, but also the existence of an explosion and so must be treated with great care. 


\subsection{GR Stability Analysis}
\label{methods_stab}

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 \citet{nagele2020}, we relied on the PN stellar evolution code to determine when collapse would occur, and our results were in rough agreement with those of \citet{chen2014}. In this work we perform a stability analysis on the normal modes of radial perturbations of a star in GR \citep{chandrasekhar1964}. 

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 \citep{shapiro1983}

        \begin{equation}
            \label{eq_N}
            \dv{}{r} \; \bigg[ \Gamma_1 \frac{P}{r^2} \dv{}{r} \; (r^2\xi)\bigg] - \frac{4}{r}\dv{P}{r} \xi +\omega^2\rho_0\xi = 0
        \end{equation}
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):
\begin{equation}
    \Gamma_1 = \pdv{\ln P}{\ln \rho_0}  \bigg\rvert_s.
\end{equation}


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 
\begin{equation}
    \omega^2_0<\omega^2_1<\omega^2_2<...
\end{equation}
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

\begin{equation}
    \omega^2_0 < 0.
\end{equation}
The corresponding equation in GR is \citep{chandrasekhar1964}
        \begin{equation}
            \nonumber
            e^{-2a-b}\dv{}{r} \; \Big[ e^{3a+b}\Gamma_1 \frac{P}{r^2} \dv{}{r} \; (e^{-a}r^2\xi)\Big] - \frac{4}{r}\dv{P}{r} \xi +e^{-2a-2b}\omega^2(P+\rho c^2)\xi 
        \end{equation}     
        \begin{equation}
            \label{eq_GR}
            - \frac{8\pi G}{c^4}e^{2b}P(P+\rho c^2)\xi - \frac{1}{P+\rho c^2} \bigg(\dv{P}{r}\bigg)^2 \xi = 0
        \end{equation}
where $a,b$ are the metric coefficients as in \citet{haemmerle2020} and the density is defined in Eq. $\ref{eq:epsilon}$. 


In order to solve Eqs. \ref{eq_N},\ref{eq_GR} for $\omega^2_0$, we adopt an iterative method similar to that outlined in exercise 6.11 of \citet{shapiro1983}. For a detailed discussion, see Appendix A. 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 \citep{chandrasekhar1964}. Fig. \ref{fig:nconv_inst} 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. \ref{fig:insta_hoshi}), the calculation is mapped to the hydrodynamical code. 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. \ref{fig:insta_hydro} 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. \ref{fig:insta_hydro}), 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 burnt most of its helium (Table \ref{tab:models}).








\begin{figure}
    \centering
    \includegraphics[width=\columnwidth]{figures/figk.pdf}

    \caption{Comparison of radial time snapshots for the 2 $\times 10^4$ ${\rm M}_\odot$ model in the HOSHI code used in this work (dashed lines) and the one used in \citet{nagele2020} (dotted lines). The inclusion of internal energy to the GR density when evaluating the PN TOV equation causes a more compact inner envelope, while the outer envelope remains largely unchanged.}
    \label{fig:HOSHI_r}
\end{figure}


\begin{figure*}
    \centering
    \includegraphics[width=2\columnwidth]{figures/figp.pdf}

    \caption{Numerical convergence of explosion energy for three parameters of HYDnuc for the 3 $\times 10^4$ ${\rm M}_\odot$ model. Left panel --- dependence on $\mathcal{V}$ (Eq. \ref{eq_varlim}) with total mesh number = 255 and isotope number = 52. Middle panel ---  dependence on number of isotopes in the nuclear network (52, 58, 61, 79, 89, 153, 300; see Appendix B, Table \ref{tab:networks}) with total mesh number = 255 and $\mathcal{V} =5\times 10^{-4}$. Right panel ---  dependence of explosion energy on number of mesh grids, with $\mathcal{V} = 5\times 10^{-4}$ and isotope number = 52.}
    \label{fig:exp_nconv}
\end{figure*}




\begin{figure}
    \centering
    \includegraphics[width=\columnwidth]{figures/figj.pdf}

    \caption{Numerical convergence of the stability analysis code for numerical polytropes with spacing d$\xi$. These polytropes are evaluated at $\Gamma_1=4/3$ for the Newtonian analysis and $\Gamma_1=4/3+\kappa GM/(rc^2) $ for the GR analysis so that we should find $\omega_0^2 = 0$. Thus, the y axis can be regarded as the error associated with our method, and its amplitude decreases for increasing numerical resolution. The number of meshes in HOSHI is variable, but is usually between 1000 and 2000, translating to somewhere on the left side of this plot. }
    \label{fig:nconv_inst}
\end{figure}



\begin{figure}
    \centering
    \includegraphics[width=\columnwidth]{figures/figo.pdf}

    \caption{GR stability analysis (Sec. \ref{methods_stab}) applied to the results of the 3 $\times 10^4$ ${\rm M}_\odot$ model in HOSHI. Upper panel --- amplitude of the fundamental mode frequency as a function of time for the General Relativistic (blue) and Newtonian (red) analyses. Unstable models, having a negative value of frequency squared, are denoted by os while stable models are denoted by crosses. Lower panel --- Central mass fraction of various isotopes as a function of time. The first instability occurs during helium burning. }
    \label{fig:insta_hoshi}
\end{figure}





\begin{figure}
    \centering
    \includegraphics[width=\columnwidth]{figures/figl.pdf}

    \caption{Stability analysis performed on the results of the HYDnuc calculation for the first unstable model of 2 $\times 10^4$ ${\rm M}_\odot$. Upper panel --- stability analysis. Lower panel --- central temperature. }
    \label{fig:insta_hydro}
\end{figure}





\begin{figure}
    \centering
    \includegraphics[width=\columnwidth]{figures/figm.pdf}

    \caption{Illustration of the polytropic criterion for the models in this paper. Upper panel --- mass of the helium core from HOSHI, determined by the mass outside which X($^1$H) $>$ 1e-5. Middle panel --- critical density (Sec. \ref{results_comp}) of a pure helium core with the mass from the top panel. Lower panel --- central density of the HOSHI models at the maximum central helium mass fraction. A comparison of the middle and lower panels shows that these models are not yet unstable according to the polytropic criterion.}
    \label{fig:polytropic}
\end{figure}

\begin{figure}
    \centering
    \includegraphics[width=\columnwidth]{figures/fign.pdf}

    \caption{Comparison of the first instability reached by two methods, Sec. \ref{methods_stab} and that of \citet{haemmerle2020}. Note that the values in this figure differ from Table \ref{tab:models} for lower masses because the first instability does not always collapse in HYDnuc (Sec. \ref{methods_stab}).} 
    \label{fig:haemmerle}
\end{figure}



\begin{figure*}
    \centering
    \includegraphics[width=2\columnwidth]{figures/fige.pdf}

    \caption{Left panel --- local energy as a function of mass coordinate for the 2.6 $\times 10^4$ ${\rm M}_\odot$ model (pulsation). The ejection criterion is $e_{\rm tot}>0$ for all $m_r > \rm M - M_{\rm ej}$ at the time of shock breakout. Right panel --- illustration of the stability of the ejection criterion. We measure the ejected mass at shock breakout and this quantity is roughly consistent with the post Newtonian escape velocity criterion near the end of the simulation. .}
    \label{fig:E_pulse}
\end{figure*}




\begin{figure*}
    \centering
    \includegraphics[width=2\columnwidth]{figures/figa.pdf}

    \caption{Time evolution of central temperature (1st panel), density (2nd panel), rate of change of specific internal energy (3rd panel) and entropy relative to the initial value (4th panel). The legend groups models by outcome, whereas colors vary with mass. The left column shows the exploding and pulsating models as a function of time. The right column shows the temperature as a function of time, while the other three panels are functions of temperature. }
    \label{fig:central}
\end{figure*}


\begin{figure*}
    \centering
    \includegraphics[width=2\columnwidth]{figures/figf.pdf}

    \caption{Upper panel --- comparison of central isotope mass fractions for the M = 3 $\times 10^4$ ${\rm M}_\odot$ model with a 58 isotope model (left) and a 61 isotope model (right). The high mass fraction of protons and $^{13}$N in the left panel facilitates artificially high carbon burning. Lower panel ---  time evolution of central $\dot{\epsilon}$.}
    \label{fig:cx_t}
\end{figure*}


\begin{figure}
    \centering
    \includegraphics[width=\columnwidth]{figures/figb.pdf}

    \caption{Total energy, kinetic energy, gravitational energy, and thermal energy (as defined in \citet{nagele2020}) in HYDnuc as a function of time for the explosion of the 3 $\times 10^4$ ${\rm M}_\odot$ model. Unless otherwise specified, subsequent figure will show the 3 $\times 10^4$ ${\rm M}_\odot$ model.}
    \label{fig:E_explo}
\end{figure}





\begin{figure*}
    \centering
    \includegraphics[width=2\columnwidth]{figures/figg.pdf}

    \caption{Initial (upper panel) and final (lower panel) isotope mass fractions as a function of the mass coordinate for the pulsation, M = 2.6 $\times 10^4$ ${\rm M}_\odot$ (left) and the explosion, M = 3 $\times 10^4$ ${\rm M}_\odot$ (right). }
    \label{fig:cx_mr}
\end{figure*}

\begin{figure}
    \centering
    \includegraphics[width=\columnwidth]{figures/figc.pdf}

    \caption{Log velocity as a function of time and mass coordinate. Shock breakout occurs about $10^5$ seconds after peak temperature.}
    \label{fig:v}
\end{figure}

\begin{figure}
    \centering
    \includegraphics[width=\columnwidth]{figures/figd.pdf}

    \caption{Log radius as a function of time and mass coordinate. }
    \label{fig:r_mr}
\end{figure}





\begin{figure}
    \centering
    \includegraphics[width=1\columnwidth]{figures/figr.pdf}

    \caption{Explosion yield (lines) compared with observed metal poor stars (crosses are observations, triangles upper limits). The line shows the fraction relative to a minimum mass of hydrogen with which it would have to mix \citet{magg2020}, and can thus be regarded as an upper limit. }
    \label{fig:XH}
\end{figure}





\begin{table*}
	\centering
	\caption{Summary table for all models. The columns are total mass, outcome of HYDnuc, mass of the isentropic core, central helium mass fraction at the start of HYDnuc, change in helium mass fraction, explosion energy, maximum central temperature, and maximum velocity of the outermost mesh, denoted $v_R$.}
	\label{tab:models}
	\begin{tabular}{|c|l|l|l|l|l|l|l|} 
		
    		M [$10^4$ ${\rm M}_\odot$] & Outcome & $M_{\rm core}$ [${\rm M}_\odot$] & $X_{\rm c}(^4$He) & $\Delta \; X_{\rm c}(^4$He) & $E_{\rm exp}$ [ergs] & max $T_{ \rm c} $ [K]& max $v_R$/c  \\
    		\hline
2 & Collapse & 10926 & 1.37e-3 &  ---  &  ---  &  ---  &  ---  \\ 
2.1 & Collapse & 11368 & 2.23e-4 &  ---  &  ---  &  ---  &  ---  \\ 
2.2 & Collapse & 11729 & 1.22e-4 &  ---  &  ---  &  ---  &  ---  \\ 
2.3 & Collapse & 12595 & 3.17e-2 &  ---  &  ---  &  ---  &  ---  \\ 
2.4 & Collapse & 13180 & 3.44e-18 &  ---  &  ---  &  ---  &  ---  \\ 
2.5 & Collapse & 13798 & 2.69e-3 &  ---  &  ---  &  ---  &  ---  \\ 
2.6 & Pulsation & 14772 & 0.104 & 0.104 & 4.32e53 & 7.58e8 & 0.032 \\ 
2.7 & Pulsation & 14964 & 0.222 & 0.147 & 4.70e52 & 6.62e8 & 0.021 \\ 
2.8 & Collapse & 15596 & 0.713 &  ---  &  ---  &  ---  &  ---  \\ 
2.9 & Pulsation & 16183 & 0.589 & 0.153 & 7.56e53 & 7.33e8 & 0.041 \\ 
2.95 & Explosion & 16504 & 0.599 & 0.168 & 1.23e54 & 7.69e8 & 0.046 \\ 
3 & Explosion & 16817 & 0.652 & 0.152 & 1.43e54 & 8.06e8 & 0.048 \\ 
3.05 & Collapse & 17144 & 0.734 &  ---  &  ---  &  ---  &  ---  \\ 
3.1 & Collapse & 17516 & 0.794 &  ---  &  ---  &  ---  &  ---  \\ 
3.15 & Collapse & 17793 & 0.815 &  ---  &  ---  &  ---  &  ---  \\ 
3.2 & Collapse & 18091 & 0.815 &  ---  &  ---  &  ---  &  ---  \\ 
3.3 & Collapse & 18888 & 1.000 &  ---  &  ---  &  ---  &  ---  \\ 
3.4 & Collapse & 19460 & 1.000 &  ---  &  ---  &  ---  &  ---  \\ 
3.5 & Collapse & 19933 & 0.950 &  ---  &  ---  &  ---  &  ---  \\ 
4 & Collapse & 23891 & 0.960 &  ---  &  ---  &  ---  &  ---  \\ 

		\hline
	\end{tabular}
\end{table*}

\begin{table*}
	\centering
	\caption{Mass ejecta by isotope for the explosions and the pulsations. Except for the first column which is consistent with Table \ref{tab:models}, values are recorded in units of ${\rm M}_\odot$.}
	\label{tab:ejecta}
	\begin{tabular}{|c|l|l|l|l|l|l|l|l|l|} 
		
    		M [$10^4$ ${\rm M}_\odot$] &M$_{\rm ej}$ & M($^1$H)  & M($^4$He) & M($^{12}$C) & M($^{16}$O) & M($^{20}$Ne) & M($^{24}$Mg) & M($^{28}$Si) & M($^{32}$S)  \\
    		\hline
2.6 & 2808 & 1877 & 974 & < 0.1 & < 0.1 & < 0.1 & < 0.1 & < 0.1 & < 0.1 \\
2.7 & 2299 & 1584 & 759 & < 0.1 & < 0.1 & < 0.1 & < 0.1 & < 0.1 & < 0.1 \\
2.9 & 2078 & 1465 & 651 & < 0.1 & < 0.1 & < 0.1 & < 0.1 & < 0.1 & < 0.1 \\
2.95 & 29500 & 5441 & 16946 & 3006 & 2505 & 481 & 812 & 306 & 1.2 \\
3 & 30000 & 5537 & 18077 & 2986 & 1829 & 367 & 702 & 497 & 5.1 \\

		\hline
	\end{tabular}
\end{table*}





\section{Results}
\label{results}

The primary advantage of this paper is the GR stability analysis, so before discussing applications, we will compare this analysis to two previous methods. 

\subsection{Comparison to previous works}
\label{results_comp}


First, and most common in the literature \citep[e.g.][]{fuller1986,umeda2016,woods2017} 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}$ \citep{shapiro1983}. From here, the polytropic constant is determined
\begin{equation}
    K=\frac{a}{3} \bigg( \frac{3s_r}{4m_pa} \bigg)^{4/3}.
\end{equation}
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$, 
\begin{equation}
     \frac{4}{3} + \frac{\beta}{6} + \mathcal{O}(\beta^2) = \frac{4}{3} + \frac{2GM_{\rm He}\kappa}{R_{\rm crit} c^2}
    \label{eq:rcrit}
\end{equation}
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 \citet{chandrasekhar1964} by a factor of 2). Eq. \ref{eq:rcrit} 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 \citep{shapiro1983}
\begin{equation}
    \rho_{\rm crit} = \Bigg[\frac{R_{\rm crit}}{\xi_1} \bigg( \frac{K}{\pi G}\bigg)^{-1/2} \Bigg]^{-3} \propto M_{\rm He} ^ {-7/2}
\end{equation}
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. \ref{methods_stab}. Fig. \ref{fig:polytropic} 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. \ref{methods_stab}.

Next we compare our method to the results of \citet{haemmerle2020}, who also evaluate Eq. \ref{eq_GR}, 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. \ref{fig:haemmerle} 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 \ref{tab:models}. This would mean that the explosion would occur at larger mass, if we were to use the method of \citet{haemmerle2020}, and this is because their method is merely a sufficient condition for instability. 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 \citet{chen2014,nagele2020}.



\subsection{Application to GRSNe}
\label{results_explo}

\subsubsection{HYDnuc explosion}

As hinted at in previous sections, we find two GRSNe as well as several pulsations (Table \ref{tab:models}). The pulsations are less energetic events which eject a portion of the envelope (Table \ref{tab:ejecta}). The relationship between the GRSNe and these pulsations is likely similar to the relationship between pair instability supernovae and the pulsational pair instability process \citep{woosley2017}. Previous papers did not find any pulsations because the models with mass lower than the explosion mass collapsed from entering the pair unstable region, as opposed to being triggered by the GR instability. Thus, the discovery of pulsations is a straightforward consequence of using the more accurate GR instability analysis. For the pulsations, we determine the ejecta mass in Table \ref{tab:models} using the local energy, e(r) (Fig. \ref{fig:E_pulse}, left panel), which is the integrand of the global energies defined in \citet{nagele2020}. We measure this quantity at shock breakout (right panel), when some of the energy is still in the form of thermal energy. The energy evolution after shock breakout may not be completely accurate. We confirm that the escape velocity criterion (right panel) converges roughly to this same value, validating the use of this ejection criterion. Note that $2.9 \times 10^4$ ${\rm M}_\odot$ is a marginal case. Although it is not an explosion, it does not have the steady behaviour of Fig. \ref{fig:E_pulse} after shock breakout, and we expect the value in Table \ref{tab:ejecta} to underestimate the ejecta mass for this case. 


The mass range for the explosion follows from straightforward considerations. Models which experience the instability before helium burning has reduced the binding energy of the star cannot explode. On the other hand, models with helium mass fractions less than ten percent also do not explode or pulsate because they lack fuel for alpha capture reactions sufficient to halt the collapse.

The pulsating models all have lower mass than the exploding models and the most massive model being the most energetic before a sharp drop-off is reminiscent of Fig. 6 of \citet{nagele2020}. This is a likely characteristic of the GRSNe even if the mass range found by the current analysis is not completely correct.


Fig. \ref{fig:central} shows the behavior of the central temperature, baryonic density, nuclear energy generation rate, and entropy as a function of time. For the exploding and pulsating models (left column), temperature and density increase and decrease smoothly, as was the case for the GRSNe in \citet{chen2014,nagele2020}, and the nuclear energy generation rate is also smooth. The models which collapse (right column) also show relatively smooth increases in temperature and density, but the energy generation rate and change in entropy are more complicated, as the stars enter different phases of burning and photodissociation. Roughly speaking, the first three peaks in $\dot{\epsilon}_c$ (first peak at $10^9$ K) correspond to carbon burning, sulfur burning, and calcium to iron peak element burning. Stars with a more evolved core, such as 2.4 $\times $ $10^4$ ${\rm M}_\odot$ do not have much carbon remaining, and the first peak corresponds to neon burning. At the other extreme, stars with reserves of hydrogen, such as 4 $\times $ $10^4$ ${\rm M}_\odot$ do not exhibit peaks because the presence of protons increases the number of possible reactions. After the core becomes iron/nickel, the next major reactions are assisted by free nucleons created by photodissociation (one peak), which is finally followed by photodissociation of nickel (one peak) and photodissociation of helium (two peaks). In reality, neutrino reaction would then begin to dominate, but HYDnuc does not include most of the relevant neutrino reactions. 




For the explosion of 3 $\times 10^4$ ${\rm M}_\odot$--- which we will use as an example in the figures--- Fig. \ref{fig:E_explo} shows the time evolution of the energy quantities, including the total energy which eventually determines the explosion energy (Table \ref{tab:models}). There are three phases, the initial contracting phases with $E_{\rm tot}<0$, the pre-shock breakout phase with $E_{\rm tot}>E_{\rm kin}>0$, and the post-shock breakout phase with $E_{\rm tot}=E_{\rm kin}>0$. Fig. \ref{fig:cx_mr} shows the isotope mass fraction as a function of mass coordinate for the initial (upper) and final (lower) time steps. The SMS cores are initially isentropic and shell hydrogen burning often occurs. The envelope is divided into many convective layers, the exact layout of which can alter the stability of the star \citep{nagele2020}. The final isotope distributions show that the majority of the nuclear burning takes place within the inner 5000 ${\rm M}_\odot$ for the pulsating model and the inner 10000 ${\rm M}_\odot$ for the exploding model. For the exploding model, the star is totally disrupted (e.g. Fig. \ref{fig:v}), so these elements will be ejected into the inter stellar medium (ISM) (Table \ref{tab:ejecta}).


After the explosive nuclear burning, the inwards velocity rapidly reverses and the shock propagates towards the surface of the SMS with a typical velocity of a few percent the speed of light (Table \ref{tab:models}); shock breakout occurs on a timescale of $10^5$ s (Fig. \ref{fig:v}). The initial inwards velocity is largest in the envelope, and we emphasize that the GRSN involves the collapse of the entire star, not just the core. This is another reason why the analysis in Sec. \ref{methods_stab} is necessary, as an analysis of the core alone will not always capture the instability. Fig. \ref{fig:r_mr} shows the time evolution of radius as a function of mass coordinate. The path of the shock can be seen in the upper part of the figure.  

\subsubsection{Comparison to observed metal poor stars}

We find a wider explosion window than previous works, as well as pulsations which could later explode. However, the explosion window is still narrow compared to the feasible mass range of SMSs. This narrowness, in combination with the probable paucity of SMSs themselves, means that the likelihood of a GRSN having occurred in the Milky Way volume ($r = 1.4 \pm 0.4$ Mpc at z = 127 \citealt{griffen2016}) is low. With that disclaimer in place, however, we will now compare the GRSN ejecta to observed metal poor stars. 

In this section, we will consider only the mono-enrichment scenario (see e.g. \citealt{hartwig2018a,hartwig2019}), so that the metals in the metal poor star have come exclusively from the GRSN ejecta and thus we ignore potential pollution from ISM accretion. we determine the minimum mass of ISM material with which the ejecta could mix via Eq. 2 of \citet{magg2020} where the explosion energy is taken from Table \ref{tab:models} and we use the fiducial value of the number density from \citet{magg2020}, $n_0 = 1$ cm$^{-3}$. We combine this mass with our ejecta and plot the inferred abundances relative to hydrogen, relative to solar (Fig. \ref{fig:XH}). Since the mixing mass is a lower limit, the abundance ratios are an upper limit. 

Also shown are the inferred abundances of the twenty stars in the SAGA database \citep{suda2008} with the lowest values of Fe / H: \citet{keller2014,ezzeddine2019,aguado2018,aoki2006,frebel2008,christlieb2004,bonifacio2018,bonifacio2015,frebel2015,caffau2016,norris2007,caffau2011,hansen2014,hansen2015,starkenburg2018,roederer2014a,roederer2014b}. We select these stars because the GRSN produces negligible amounts of iron. Fig. \ref{fig:XH} shows that our inferred yields do not match any of the observed stars. Even for the star of \citet{keller2014} which has strict upper limits on [Fe/H], our yield misses the observed value of calcium by several orders of magnitude. 

While we can definitively rule out the mono-enrichment scenario for metal poor stars in the vicinity of our galaxy, the multi-enrichment scenario is more challenging to completely rule out. The best we can do at current is to note that although many metal poor stars are carbon enriched, they are not generally silicon and magnesium enriched, which is strong evidence against an SMS being involved in the multi-enrichment scenario. Thus we can say that it is likely the case that there did not exist a non-accreting SMS in the explosion mass range, in the Milky Way volume.


\subsection{Application to SMS Collapse and neutrino emission}
\label{results_col}

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 \citep{nagele2021}. 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 \citet{nagele2021}, 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. 



\section{Discussion}
\label{discussion}

Using the general relativistic stability analysis in Sec. \ref{methods_stab}, we can more accurately predict when a SMS will be unstable and will collapse explode or pulse. This is necessary because the timescale of the collapse is many orders of magnitude shorter than the evolutionary timescale, meaning that our stellar evolution code often misses the GR instabilities. For some models, the instability is countered by increased nuclear burning, but for masses around 3 $\times 10^4$ ${\rm M}_\odot$, this is not the case. The 2.95 $\times 10^4$ ${\rm M}_\odot$ and 3 $\times 10^4$ ${\rm M}_\odot$ models explode in GRSNe, while several lower mass models pulsate and eject portions of their envelope. The final fate of the pulsating models is unclear, as they will reenter the evolutionary track with different properties, most notably the chemical composition of the core and total energy. If multiple pulsations were to occur, it could cause a collisional supernova. 

In comparison to the GRSNe from \citet{chen2014,nagele2020}, our GRSNe have a much lower explosion energy, though if we were to use an unphysically small nuclear network, we would find a GRSN with comparably large energies. The composition of our explosion ejecta is similar to the previous works (Table \ref{tab:ejecta}). On the other hand, the mass range of the explosion is wider because our more accurate determination of the GR instability means that these models are not near the pair instability region. 

\citet{moriya2021} found that the observational duration of a GRSN will be on the order of $10^{2-3}$ seconds for the peak emission and $10^{2-3}$ days for the plateau, which shares some similarities with Type IIP SNe. Because the GRSN would have to occur in the high redshift universe, this means that the observer duration is longer by about a factor of ten. \citet{moriya2021} demonstrated that the GRSN plateau may be differentiated from other persistent sources if it is observed in multiple bands. In the future, we plan to investigate if the lower energy of our GRSN would signifigantly alter any of their findings.

We determined that none of the observed metal poor stars match the inferred abundance pattern from the GRSNe, and from this we conclude that none of these stars were singly enriched by a GRSN. However, we note that the multi-enrichment scenario cannot be completley ruled out.

As in all numerical studies, there are numerous uncertainties. In HOSHI, along with the usual sources of numerical error, we also need to consider the error due to using the post Newtonian approximation. After the inclusion of internal energy to density (Eq. \ref{eq:epsilon}) in the TOV equation, the post Newtonian pressure gradient matches the TOV pressure gradient to one part in $10^5$. Although this level of accuracy likely supersedes numerical error, we have to keep in mind that the TOV equation assumes a hydrostatic configuration, which, for instance, is not true as the star contracts towards the end of hydrogen burning. Finally, and perhaps most importantly, our evolutionary models are not rotating, whereas real SMSs are expected to initially be medium rotators \citep{hammerle2018} which may spin up over the course of their lifetimes \citep{maeder2001}. We intend to more fully investigate the effects of rotation in the future.  

Regarding the GR stability analysis, we were able to quantify the error of $|\omega_0^2|$ for polytropes (Fig. \ref{fig:nconv_inst}) as being around $10^{-7,-8}$. As seen in Fig. \ref{fig:insta_hoshi}, typical values of $\omega_0^2$ are greater than this error. Although we demonstrated that the error decreases for increasing resolution (Fig. \ref{fig:nconv_inst}) the gain in accuracy is low compared to the gain in computational time that an increased mesh number of an order of magnitude or two would require. 

The largest source of uncertainty in HYDnuc may be the coordinate perturbation at the start of the calculation. We sought to minimize this perturbation via Eq. \ref{eq:epsilon}, but the possibility remains that the perturbation size may affect the explosion energy. Indeed in Fig. \ref{fig:exp_nconv} left panel we show how sensitive the explosion energy can be to the collapse timescale, which would vary with a different initial perturbation. That being said, a change in the initial perturbation would likely only shift the mass range of the explosion. As far as other sources of error go, \citet{chen2014} showed that multidimensional effects may not play a big role in the GRSN, and they are also not thought to contribute to the instability \citep{chandrasekhar1965}. We verified that radiative and convective energy transport do not effect the explosion outcome and tested that the explosion energy does not depend on the choice of the carbon alpha capture rate. 


In the future, we plan to apply the GR stability analysis to rotating SMSs, as well as verifying our current results with multidimensional simulations. We also intend to asses the observability of the GRSN using methods similar to those in \citet{moriya2021}. Finally, we will investigate the possibility of multiple pulsations. 



\section*{Data Availability}

The data underlying this article will be shared on reasonable request to the corresponding author.

\section*{Acknowledgements}

This study was supported in part by the Grant-in-Aid for the Scientific Research of Japan Society for the Promotion of Science (JSPS, Nos. JP19K03837, JP20H01905, JP20H00158, JP21H01123) and by Grant-in-Aid for Scientific Research on Innovative areas (JP17H06357, JP17H06365) from the Ministry of Education, Culture, Sports, Science and Technology (MEXT), Japan. For providing high performance computing resources, YITP, Kyoto University is acknowledged. K.S. would also like to acknowledge computing resources at KEK and RCNP Osaka University. C.N. would like to thank Tilman Hartwig for discussions of metal poor stars. 









\bibliographystyle{mnras}
