RELATIVISTIC COLLAPSE AND EXPLOSION OF ROTATING SUPERMASSIVE STARS WITH THERMONUCLEAR EFFECTS

2.1.1. BSSN Formulation

We follow the 3+1 formulation in which the spacetime is foliated into a set of non-intersecting spacelike hypersurfaces. In this approach, the line element is written in the form

$$\begin{equation} ds^2 = -(\alpha ^{2} -\beta _{i}\beta ^{i}) dt^2 + 2 \beta _{i} dx^{i} dt +\gamma _{ij} dx^{i} dx^{j}, \end{equation} \tag{ 1 }$$

where α, βⁱ, and γ_ij are the lapse function, the shift three-vector, and the three-metric, respectively. The latter is defined by

$$\begin{equation} \gamma _{\mu \nu }=g_{\mu \nu }+n_{\mu }n_{\nu }, \end{equation} \tag{ 2 }$$

where n^μ is a timelike unit-normal vector orthogonal to a spacelike hypersurface.

We make use of the BSSN formulation (Nakamura et al. 1987; Shibata & Nakamura 1995; Baumgarte & Shapiro 1999c) to solve the Einstein equations. Initially, a conformal factor ϕ is introduced, and the conformally related metric is written as

$$\begin{equation} {\tilde{\gamma }}_{ij}=e^{-4\phi }\gamma _{ij} \end{equation} \tag{ 3 }$$

such that the determinant of the conformal metric, $\tilde{\gamma }_{ij}$, is unity and ϕ = ln (γ)/12, where $\gamma =\det (\gamma _{ij})$. We also define the conformally related traceless part of the extrinsic curvature K_ij,

$$\begin{equation} {\tilde{A}}_{ij}=e^{-4\phi }A_{ij}=e^{-4\phi }\left(K_{ij}-\frac{1}{3}\gamma _{ij}K\right), \end{equation} \tag{ 4 }$$

where K is the trace of the extrinsic curvature. We evolve the conformal factor, defined as χ ≡ e^−4ϕ (Campanelli et al. 2006), and the auxiliary variables $\tilde{\Gamma }^{i}$, known as the conformal connection functions, defined as

$$\begin{equation} \tilde{\Gamma }^{i}\equiv \tilde{\gamma }^{jk}\tilde{\Gamma }^{i}_{\hspace{5.69046pt}jk}=-\partial _{j}\tilde{\gamma }^{ij}, \end{equation} \tag{ 5 }$$

where $\tilde{\Gamma }^{i}_{\hspace{5.69046pt}jk}$ are the connection coefficients associated with $\tilde{\gamma }_{ij}$.

During the evolution, we also enforce the constraints ${\rm Tr} (\tilde{A}_{ij})=0$ and ${\rm det} (\tilde{\gamma }_{ij})=1$ at every time step.

We use the Cartoon method (Alcubierre et al. 2001) to impose axisymmetry while using Cartesian coordinates.

2.1.2. Gauge Choices

In addition to the BSSN spacetime variables (${\tilde{\gamma }_{ij}},{\tilde{A}_{ij}},K, \chi,\tilde{\Gamma }^{i}$), there are two more quantities left undetermined, the lapse, α, and the shift vector, βⁱ. We used the so-called non-advective 1+log slicing (Bona et al. 1997) by dropping the advective term in the "1+log" slicing condition. In this case, the slicing condition takes the form

$$\begin{equation} \partial _{t}\alpha = -2 \alpha K. \end{equation} \tag{ 6 }$$

For the shift vector, we choose the "Gamma-freezing condition" (Alcubierre et al. 2003), written as

$$\begin{equation} \hspace{-13.1pt}\partial _{t}\beta ^{i} = \frac{3}{4}B^{i},\hspace{13.1pt} \end{equation} \tag{ 7 }$$

$$\begin{equation} \partial _{t}B^{i} = \partial _{t}\tilde{\Gamma }^{i}-\eta B^i, \end{equation} \tag{ 8 }$$

where η is a constant that acts as a damping term, originally introduced both to prevent long-term drift of the metric functions and to prevent oscillations of the shift vector.

The general relativistic hydrodynamic equations, expressed through the conservation equations for the stress-energy tensor T^μν and the continuity equation, are

$$\begin{equation} \nabla _\mu T^{\mu \nu } = 0\;,\;\;\;\;\;\; \nabla _\mu \left(\rho u^{\mu }\right) = 0, \end{equation} \tag{ 9 }$$

where ρ is the rest-mass density, u^μ is the fluid four-velocity, and ∇ is the covariant derivative with respect to the spacetime metric. Following Shibata (2003), the general relativistic hydrodynamic equations are written in a conservative form in cylindrical coordinates. Since the Einstein equations are solved only in the y = 0 plane with Cartesian coordinates (two-dimensional), the hydrodynamic equations are rewritten in Cartesian coordinates for y = 0. The following definitions for the hydrodynamical variables are used:

$$\begin{equation} \rho _{*}\equiv \rho W e^{6\phi }, \end{equation} \tag{ 10 }$$

$$\begin{equation} v^{i}\equiv \frac{u^{i}}{u^{t}}=-\beta ^{i}+\alpha \gamma ^{ij}\frac{\mbox{\^{u}}_{j}}{hW}, \end{equation} \tag{ 11 }$$

$$\begin{equation} \mbox{\^{u}}_{i}\equiv h u_{i}, \end{equation} \tag{ 12 }$$

$$\begin{equation} \mbox{\^{e}}\equiv \frac{e^{6\phi }}{\rho _{*}}T_{\mu \nu }n^{\mu }n^{\nu }=hW-\frac{P}{\rho W}, \end{equation} \tag{ 13 }$$

$$\begin{equation} W\equiv \alpha u^{t}, \end{equation} \tag{ 14 }$$

where W and h are the Lorentz factor and the specific fluid enthalpy, respectively, and P is the pressure. The conserved variables are ρ_*, $J_{i}=\rho _{*}\mbox{\^{u}}_{i}$, and $E_{*}= \rho _{*}\mbox{\^{e}}$. We refer to Shibata (2003) for further details.

Isentropic SMSs are self-gravitating equilibrium configurations of masses in the range of 10⁴–10⁸ M_☉, which are mainly supported by radiation pressure, while the pressure of electron–positron pairs and of the baryon gas is only a minor contribution to the EOS. Such configurations are well described by Newtonian polytropes with polytropic index n = 3 (adiabatic index Γ = 4/3). The ratio of the gas pressure to the total pressure (β) for spherical SMSs can be written as (Fowler & Hoyle 1966)

$$\begin{equation} \beta = \frac{P_g}{P_{{\rm tot}}} \approx \frac{4.3}{\mu } \left(\frac{M_{\odot }}{M}\right)^{1/2}, \end{equation} \tag{ 15 }$$

where μ is the mean molecular weight. Thus, β  ≈  10⁻² for M  ≈  10⁶ M_☉.

Since nuclear burning timescales are too long for M  ≳  10⁴ M_☉, evolution of SMSs proceeds on the Kelvin–Helmholtz timescale and is driven by the loss of energy and entropy by radiation, as well as loss of angular momentum via mass shedding in the case of rotating configurations.

Although corrections due to the non-relativistic gas of baryons and electrons and general relativistic effects are small, they cannot be neglected for the evolution. First, gas corrections raise the adiabatic index slightly above 4/3:

$$\begin{equation} \Gamma \approx \frac{4}{3}+ \frac{\beta }{6}+0(\beta ^2). \end{equation} \tag{ 16 }$$

Second, general relativistic corrections lead to the existence of a maximum for the equilibrium mass as a function of the central density. For spherical SMSs this means that for a given mass the star evolves to a critical density beyond which it is dynamically unstable against radial perturbations (Chandrasekhar 1964):

$$\begin{equation} \rho _{{\rm crit}} = 1.994 \times 10^{18} \left(\frac{0.5}{\mu }\right)^3\left(\frac{M_\odot }{M}\right)^{7/2} \,{\rm g} \,{\rm cm^{-3}}. \end{equation} \tag{ 17 }$$

The onset of the instability also corresponds to a critical value of the adiabatic index Γ_crit, i.e., configurations become unstable when the adiabatic index drops below the critical value:

$$\begin{equation} \Gamma _{{\rm crit}} = \frac{4}{3}+1.12\frac{2GM}{Rc^2}. \end{equation} \tag{ 18 }$$

This happens when the stabilizing gas contribution to the EOS does not raise the adiabatic index above 4/3 to compensate for the destabilizing effect of general relativity expressed by the second term on the right-hand side of Equation (18).

Rotation can stabilize configurations against the radial instability. The stability of rotating SMSs with uniform rotation was analyzed by Baumgarte & Shapiro (1999a, 1999b). They found that stars at the onset of the instability have an equatorial radius R  ≈  640GM/c², a spin parameter q ≡ cJ/GM²  ≈  0.97, and a ratio of rotational kinetic energy to the gravitational binding energy of T/W  ≈  0.009.

To close the system of hydrodynamic equations (Equation (9)), we need to define the EOS. We follow a treatment that separately includes the baryon contribution on the one hand and photons and electron–positron pair contributions, in a tabulated form, on the other hand. The baryon contribution is given by the analytic expressions for the pressure and specific internal energy

$$\begin{equation} P_{b}=\frac{\mathcal {R}\rho T}{\mu _{b}},\hspace{8.5pt} \end{equation} \tag{ 19 }$$

$$\begin{equation} \epsilon _{b}=\frac{2}{3}\frac{\mathcal {R}\rho T}{\mu _b}, \end{equation} \tag{ 20 }$$

where $\mathcal {R}$ is the universal gas constant, T is the temperature, _b is the baryon specific internal energy, and μ_b is the mean molecular weight due to ions, which can be expressed as a function of the mass fractions of hydrogen (X), helium (Y), and heavier elements (metals; Z_CNO) as

$$\begin{equation} \frac{1}{\mu _b}\approx X + \frac{Y}{4}+\frac{Z_{{\rm CNO}}}{\langle A\rangle }, \end{equation} \tag{ 21 }$$

where 〈A〉 is the average atomic mass of the heavy elements. We assume that the composition of SMSs (approximately that of primordial gas) has a mass fraction of hydrogen X = 0.75 − Z_CNO and helium Y = 0.25, where the metallicity Z_CNO = 1 − X − Y is an initial parameter, typically of the order of Z_CNO  ∼  10⁻³ (see Table 1 for details). Thus, for the initial compositions that we consider the mean molecular weight of baryons is μ_b  ≈  1.23 (i.e., corresponding to a molecular weight for both ions and electrons of μ  ≈  0.59).

Effects associated with photons and the creation of electron–positron pairs are taken into account employing a tabulated EOS. At temperatures above 10⁹ K, not all the energy is used to increase the temperature and pressure, but part of the photon energy is used to create the rest mass of the electron–positron pairs. As a result of pair creation, the adiabatic index of the star decreases, which means that the stability of the star is reduced.

Given the specific internal energy, , and rest-mass density, ρ, as evolved by the hydrodynamic equations, it is possible to compute the temperature T by a Newton–Raphson algorithm that solves the equation *(ρ, T) = for T,

$$\begin{equation} T_{n+1} = T_{n} - (\epsilon ^{*}(\rho,T_{n})-\epsilon ) \left.\left(\frac{\partial \epsilon ^{*}(\rho,T)}{\partial T}\right)\right|_{T_{n}}^{-1}, \end{equation} \tag{ 22 }$$

where n is the iteration counter.

In order to avoid the small time steps and CPU-time demands connected with the solution of a nuclear reaction network coupled to the hydrodynamic evolution, we apply an approximate method to take into account the basic effects of nuclear burning on the dynamics of the collapsing SMSs. We compute the nuclear energy release rates by hydrogen burning (through the pp-chain, cold and hot CNO cycles, and their breakout by the rp-process) and helium burning (through the 3-α reaction) as a function of rest-mass density, temperature, and mass fractions of hydrogen X, helium Y, and CNO metallicity Z_CNO. These nuclear energy generation rates are added as a source term on the right-hand side of the evolution equation for the conserved quantity E_*.

The change rates of the energy density due to nuclear reactions, in the fluid frame, expressed in units of erg cm⁻³ s⁻¹ are given by

• 1.  
  pp-chain (Clayton 1983):

  $$\begin{eqnarray} \left(\frac{\partial e}{\partial t}\right)_{{\rm pp}} &=& \rho \big(2.38\times 10^{6}\rho g_{11}X^{2}T_{6}^{-0.6666} \nonumber \\ && e^{-33.80/T_{6}^{0.3333}}\big), \end{eqnarray} \tag{ 23 }$$

  where T₆ = T/10⁶ K and g₁₁ is given by

  $$\begin{eqnarray} g_{11} &=& 1+0.0123T_6^{0.3333}+0.0109T_6^{0.66666}\nonumber \\ && +0.0009T_6. \end{eqnarray} \tag{ 24 }$$
• 2.  
  3-α (Wiescher et al. 1999):

  $$\begin{eqnarray} \left(\frac{\partial e}{\partial t}\right)_{3\alpha } &= \rho \big(5.1\times 10^8\rho ^2Y^3T_9^{-3} e^{-4.4/T_9}\big), \end{eqnarray} \tag{ 25 }$$

  where T₉ = T/10⁹ K.
• 3.  
  Cold-CNO cycle (Shen & Bildsten 2007):

  $$\begin{eqnarray} \left(\frac{\partial e}{\partial t}\right)_{{\rm CCNO}} &=& 4.4\times 10^{25}\rho ^2 X Z_{{\rm CNO}} \big(T_9^{-2/3}e^{-15.231/T_9^{1/3}} \nonumber \\ && + 8.3\times 10^{-5}T_9^{-3/2}e^{-3.0057/T_9}\big). \end{eqnarray} \tag{ 26 }$$
• 4.  
  Hot-CNO cycle (Wiescher et al. 1999):

  $$\begin{eqnarray} \left(\frac{\partial e}{\partial t}\right)_{{\rm HCNO}} &= 4.6\times 10^{15}\rho Z_{{\rm CNO}}. \end{eqnarray} \tag{ 27 }$$
• 5.  
  rp-process (Wiescher et al. 1999):

  $$\begin{eqnarray} \left(\frac{\partial e}{\partial t}\right)_{{\rm rp}} &=& \rho \big(1.77\times 10^{16}\rho Y Z_{{\rm CNO}} \nonumber \\ &&\times 29.96 T_9^{-3/2} e^{-5.85/T_9}\big). \end{eqnarray} \tag{ 28 }$$

Since we follow a single fluid approach in which we solve only the hydrodynamic equations (Equation (9), i.e., we do not solve additional advection equations for the abundances of hydrogen, helium, and metals), the elemental abundances during the time evolution are fixed. Nevertheless, this assumption most possibly does not significantly affect the estimate of the threshold metallicity needed to produce a thermal bounce in collapsing SMSs. The average energy release through the 3-α reaction is about 7.275 MeV for each ¹²C nucleus formed. Since the total energy due to helium burning for exploding models is ∼10⁴⁵ erg (e.g., 9.0 × 10⁴⁴ erg for model S1.c) and even considering that this energy is released mostly in a central region of the SMS containing 10⁴ M_☉ of its rest mass (Fuller et al. 1986), it is easy to show that the change in the metallicity is of the order of 10⁻¹¹. Therefore, the increase of the metallicity in models experiencing a thermal bounce is much smaller than the critical metallicities needed to trigger the explosions. Similarly, the average change in the mass fraction of hydrogen due to the cold and hot CNO cycles is expected to be ∼10% for exploding models.

After each time iteration the conserved variables (i.e., ρ_*, J_x, J_y, J_z, E_*) are updated and the primitive hydrodynamical variables (i.e., ρ, v^x, v^y, v^z, ) have to be recovered. The recovery is done in such a way that it allows for the use of a general EOS of the form P = P(ρ, ). We calculate the function f(P*) = P(ρ*, *) − P*, where ρ* and * depend only on the conserved quantities and the pressure guess P*. The new pressure is then computed iteratively by a Newton–Raphson method until the desired convergence is achieved.

The EOS allows us to compute the neutrino losses due to the following processes, which become most relevant just before BH formation.

• 1.  
  Pair annihilation ($e^{+}+e^{-} \rightarrow \bar{\nu } + \nu$): the most important process above 10⁹ K. Due to the large mean free path of neutrinos in the stellar medium at the densities of SMSs, the energy loss by neutrinos can be significant. For a 10⁶ M_☉ SMS most of the energy release in the form of neutrinos originates from this process. The rates are computed using the fitting formula given by Itoh et al. (1996).
• 2.  
  Photo-neutrino emission ($\gamma +e^{\pm } \rightarrow e^{\pm }+ \bar{\nu }+ \nu$): dominates at low temperatures T  ≲  4 × 10⁸ K and densities ρ  ≲  10⁵ g cm⁻³ (Itoh et al. 1996).
• 3.  
  Plasmon decay ($\gamma \rightarrow \bar{\nu }+\nu$): this is the least relevant process for the conditions encountered by the models we have considered because its importance increases at higher densities than those present in SMSs. The rates are computed using the fitting formula given by Haft et al. (1994).

Since it is not possible to follow the gravitational collapse of an SMS from the early stages to the phase of BH formation with a uniform Cartesian grid (the necessary fine zoning would be computationally too demanding), we adopt a regridding procedure (Shibata & Shapiro 2002). During the initial phase of the collapse we rezone the computational domain by moving the outer boundary inward, decreasing the grid spacing while keeping the initial number of grid points fixed. Initially we use N × N = 400 × 400 grid points and place the outer boundary at L  ≈  1.5r_e, where r_e is the equatorial radius of the star. Rezoning onto the new grid is done using a polynomial interpolation. We repeat this procedure three to four times until the collapse timescale in the central region is much shorter than in the outer parts. At this point, we both decrease the grid spacing and also increase the number of grid points N in dependence of the lapse function typically as follows: N × N = 800 × 800 if 0.8 > α  >  0.6, N × N = 1200 × 1200 if 0.6 > α  >  0.4, and N × N = 1800 × 1800 if α < 0.4. This procedure ensures the error in the conservation of the total rest mass to be less than 2% on the finest computational domain.

To deal with the spacetime singularity from the newly formed BH, we use the method of excising the matter content in a region within the horizon as proposed by Hawke et al. (2005) once an AH is found. This excision is done only for the hydrodynamical variables, and the coordinate radius of the excised region is allowed to increase in time. On the other hand, we use neither excision nor artificial dissipation terms for the spacetime evolution and solely rely on the gauge conditions.

Here we define some of the quantities listed in Table 1. We compute the total rest mass M_* and the gravitational mass M as

$$\begin{equation} M_{*}= 4\pi \int _{0}^{L}{x}dx\int _{0}^{L}{\rho _{*}}dz, \end{equation} \tag{ 29 }$$

$$\begin{eqnarray} M&=& {-}2\int _{0}^{L}{x}dx\int _{0}^{L}{}dz\left[-2\pi E e^{5\phi }+\frac{e^{\phi }}{8}\tilde{R}\right. \nonumber \\ &&\left.-\,\frac{e^{5\phi }}{8}\left(\tilde{A}_{ij}\tilde{A}^{ij}-\frac{2}{3}K^{2}\right)\right], \end{eqnarray} \tag{ 30 }$$

where E = n_μn_νT^μν (n^μ being the unit normal to the hypersurface) and $\tilde{R}$ is the scalar curvature associated with the conformal metric $\tilde{\gamma }_{ij}$.

The rotational kinetic energy T_k and the gravitational potential energy W are given by

$$\begin{equation} T_k=2\pi \int _{0}^{L}{x^{2}}dx\int _{0}^{L}{\rho _{*}\mbox{\^{u}}_{y}\Omega }dz, \end{equation} \tag{ 31 }$$

where Ω is the angular velocity, and

$$\begin{equation} W=M-(M_{*}+T_k+E_{{\rm int}}), \end{equation} \tag{ 32 }$$

where the internal energy is computed as

$$\begin{equation} E_{{\rm int}}= 4\pi \int _{0}^{L}{x}dx\int _{0}^{L}{\rho _{*}\epsilon }dz. \end{equation} \tag{ 33 }$$

In axisymmetry the AH equation becomes a nonlinear ordinary differential equation for the AH shape function, h = h(θ) (Shibata 1997; Thornburg 2007). We employ an AH finder that solves this ordinary differential equation by a shooting method using ∂_θh(θ = 0) = 0 and ∂_θh(θ = π/2) = 0 as boundary conditions. We define the mass of the AH as

$$\begin{equation} M_{{\rm AH}}=\sqrt{\frac{\mathcal {A}}{16\pi }}, \end{equation} \tag{ 34 }$$

where $\mathcal {A}$ is the area of the AH.

Axisymmetric calculations without rotation (i.e., models S1 and S2) retain the spherical symmetry of the initial conditions. There are no physical phenomena like convective or overturn instabilities¹ to produce asphericity, i.e., we can directly compare our two-dimensional non-rotating models with those computed in spherical symmetry (one-dimensional calculations) by Fuller et al. (1986) and Linke et al. (2001).

The main differences with respect to the results obtained by Fuller et al. (1986) are most likely due to two reasons. First, we apply a fully general relativistic treatment while they used a post-Newtonian treatment of gravity, and second, there are differences in the treatment of nuclear burning (Fuller et al. 1986 solved the relevant nuclear network without the approximations adopted in our work; see Section 3.3 for further details). Despite these differences, the results agree fairly well. As discussed in detail in Section 7.1, the initial metallicities required to produce an explosion are similar, and a thermal bounce can be produced only if sufficient energy is liberated during the phase when the HCNO cycle is active.

The main difference with respect to the work of Linke et al. (2001) resides in the formulation of Einstein's field equations, in particular, in the foliation of the spacetime (foliation into a set spacelike hypersurfaces versus a foliation with outgoing null hypersurfaces). In order to compare with the results of Linke et al. (2001), we computed the redshifted total energy output for a model with the same rest mass, conducting the time integration until approximately the same evolutionary stage as in Linke et al. (2001). We find that the total energies released in neutrinos differ by less than 10% (for more details see Section 7.3).

Previous simulations of SMS collapse to a BH in general relativity have been performed with a Γ-law EOS with Γ = 4/3 (with the only exception being the work of Linke et al. 2001). In order to elucidate the influence of the EOS on the dynamics of collapsing SMSs, we performed three simulations of the same initial model (model R1.0, a marginally unstable uniformly rotating SMS with zero initial metallicity) without nuclear burning effects, and with three EOSs: a Γ-law EOS with Γ = 4/3 (i.e., a similar setup as in Shibata & Shapiro 2002) and the microphysical EOS with and without including electrons and the e^± pairs, i.e., for the last EOS case we consider Equations (19) and (20) for the baryons plus _γ = aT⁴ and P_γ = (1/3)_γ for the photons (where a is the radiation density constant). We denote the Γ-EOS as EOS-0, the full microphysical EOS, our canonical one for the studies of this work as EOS-1, and the reduced microphysical case as EOS-2.

In Figure 1 we show with a dotted line the time evolution of the central density of model R1.0 with EOS-0, and with a dashed (solid) line the time evolution of the central density with EOS-2 (EOS-1). The first thing to note is that the collapse timescale obtained with the Γ-law EOS is shorter than that obtained with the microphysical EOS (both EOS-1 and EOS-2) because the ion pressure contribution to the EOS raises the adiabatic index above 4/3 (see Equation (16)). This increase in the adiabatic index helps to stabilize the star against the gravitational instability and therefore delays the collapse.

Zoom In Zoom Out Reset image size

On the other hand, the effect of pair creation reduces the adiabatic index below 4/3 at T  ≳  10⁹ K. This explains the differences between the solid and dashed lines in Figure 1 at central densities ρ_c  ≳  10 g cm⁻³, which correspond to central temperatures T_c  ≳  10⁹ K. Once the collapse enters this regime, pair creation becomes relevant enough to reduce the adiabatic index below 4/3, which destabilizes the collapsing star. Compared to previous works (Shibata & Shapiro 2002), the use of a microphysical EOS instead of a Γ-law EOS delays the collapse (mostly due to the baryons while e^± destabilize) of an initially gravitationally unstable configuration.

We also performed two-dimensional axisymmetric simulations of differentially rotating SMSs. First, we investigated the influence of the EOS on gravitationally unstable stars using model D1 as a reference. This model corresponds, within the accuracy to which the initial conditions can be reproduced, to a differentially rotating unstable SMS discussed by Saijo & Hawke (2009, i.e., their model I). Results are displayed in Figure 2. In the upper panel of this figure, we show the time evolution of the central density for model D1 with EOS-0 (dashed line) and with the microphysical EOS-1 (solid line). Opposite to the behavior in the case of the uniformly rotating SMS R1.0, the collapse timescale is longer with the Γ-law than with the microphysical EOS. The reason for this difference is that the initial central temperature (T_c  ≈  1.4 × 10⁹ K) of model D1 is an order of magnitude higher than the initial central temperature in R1.0. Therefore, electron–positron pair creation already reduces the stability of the star (by reducing Γ) during the initial stages of the collapse. This behavior is also expected to be present in three dimensions, and since the collapse timescale is reduced when using the microphysical EOS, non-axisymmetric instabilities would have even less time to grow before the formation of a BH. It reinforces the conclusions of Saijo & Hawke (2009), who showed that the three-dimensional collapse of rotating stars proceeds in an approximately axisymmetric manner.

Zoom In Zoom Out Reset image size

The lower two panels of Figure 2 display the growth of the AH mass and the disk mass (defined as the rest mass outside the AH of the newly formed BH) as a function of time (in units of gravitational mass for comparison with Figures 9 and 10 of Saijo & Hawke 2009), respectively. The values of both quantities at the end of the simulation agree, within a 5% difference, with those obtained in three dimensions by Saijo & Hawke (2009). We note that there is also good agreement (less than 5% difference) regarding the time at which an AH is first detected. These observed small differences are likely due to differences in the initial models and numerical techniques rather than to the influence of non-axisymmetric effects. This suggests that our collapse simulation with the same treatment of physics yields good agreement with the three-dimensional simulations of Saijo & Hawke (2009).

We also investigated the influence of electron–positron pair creation on the evolution of gravitationally stable differentially rotating SMSs using model D2, which is similar to the stable differentially rotating model III of Saijo & Hawke (2009). We performed three simulations of model D2 varying the EOS. In Figure 3 we show the central rest-mass density as a function of time for a Γ-law (dashed line), and for the microphysical EOS-1 (solid line) and EOS-2 (dotted line). In agreement with the results obtained by Saijo & Hawke (2009), we find that model D2 represents a stable differentially rotating SMS when a Γ-law is used. A persistent series of oscillations is triggered by the initial perturbation in the pressure. This is also the case with the microphysical EOS-2 without the inclusion of electrons and the e^± pairs. However, the time evolution of D2 is completely different when e^± pairs are taken into account. The influence of pairs is large enough to destabilize model D2 against gravitational collapse. We note that unlike all other SMSs considered in this paper, which are Γ = 4/3 models initially unstable to gravitational collapse, model D2 is an initially stable Γ = 4/3 model that becomes gravitationally unstable only by the creation of electron–positron pairs at high temperatures. Hence, using a microphysical EOS with electron–positron pairs is crucial to determine the stability of differentially rotating SMSs.

Zoom In Zoom Out Reset image size

We note that the central temperature of the initial models D1 and D2 is of the order of ≈10⁹ K. At this temperature, the main source of thermonuclear energy is hydrogen burning via the rp-process. It is, however, expected that such SMSs would previously experience a phase of hydrogen burning via the cold and hot CNO cycles, which would significantly affect the evolution of the models such that configurations with high T_c as in models D1 and D2 might never be reached. Therefore, models D1 and D2 are not particularly well suited to investigate the existence of a thermal bounce during collapse (see Section 7). Exploring in detail the parameter space for the stability of differentially rotating SMSs with the microphysical EOS, as well as the existence of a thermal bounce during the collapse phase depending on the initial stellar metallicity, is a major task on its own, which is beyond the scope of this paper. For these reasons we do not consider differentially rotating SMSs in this work.

First we consider a gravitationally unstable spherically symmetric SMS with a gravitational mass of M = 5 × 10⁵ M_☉ (S1.a, S1.b, and S1.c), which corresponds to a model extensively discussed in Fuller et al. (1986) and therefore allows for a comparison with the results presented here. Fuller et al. (1986) found that unstable spherical SMSs with M = 5 × 10⁵ M_☉ and an initial metallicity Z_CNO = 2 × 10⁻³ collapse to a BH, while models with an initial metallicity Z_CNO = 5 × 10⁻³ explode due to the nuclear energy released by the hot CNO burning. They also found that the central density and temperature at thermal bounce (where the collapse is reversed to an explosion) are ρ_{c, b} = 3.16 g cm⁻³ and T_{c, b} = 2.6 × 10⁸ K, respectively.

The left panels in Figure 4 show the time evolution of the central rest-mass density (upper panel) and central temperature (lower panel) for models S1.a, S1.c, R1.a, and R1.d, i.e., non-rotating and rotating models with a mass of M = 5 × 10⁵ M_☉. In particular, the solid lines represent the time evolution of the central density and temperature for models S1.c (Z_CNO = 7 × 10⁻³) and R1.d (Z_CNO = 5 × 10⁻⁴). As the collapse proceeds, the central density and temperature rise rapidly, which increases the nuclear energy generation rate by hydrogen burning. Since the metallicity is sufficiently high, enough energy can be liberated to increase the pressure and to produce a thermal bounce. This is the case for model S1.c. In Figure 4 we show that a thermal bounce occurs (at approximately t  ∼  7 × 10⁵ s) entirely due to the hot CNO cycle, which is the main source of thermonuclear energy at temperatures in the range 2 × 10⁸ K  ⩽  T  ⩽  5 × 10⁸ K. The rest-mass density at bounce is ρ_{c, b} = 4.8 g cm⁻³, and the temperature is T_{c, b} = 3.05 × 10⁸ K. These values, as well as the threshold metallicity needed to trigger a thermonuclear explosion (Z_CNO = 7 × 10⁻³), are higher than those found by Fuller et al. (1986) (who found that a spherical non-rotating model with the same rest mass would explode if the initial metallicity was Z_CNO = 5 × 10⁻³).

Zoom In Zoom Out Reset image size

On the other hand, the dashed lines show the time evolution of the central density and temperature for model S1.a (Z_CNO = 5 × 10⁻³). In this case, as well as for model S1.b, the collapse is not halted by the energy release and continues until an AH is found, indicating the formation of a BH.

We note that the radial velocity profiles change continuously near the time where the collapse is reversed to an explosion due to the nuclear energy released by the hot CNO burning, and an expanding shock forms only near the surface of the star at a radius R  ≈  1.365 × 10¹³ cm (i.e., R/M  ≈  180), where the rest-mass density is ≈3.5 × 10⁻⁶ g cm⁻³. We show in Figure 5 the profiles of the x-component of the three-velocity v_x along the x-axis (in the equatorial plane) for the non-rotating spherical stars S1.a (dashed lines) and S1.c (solid lines) at three different time slices near the time at which model S1.c experiences a thermal bounce. Velocity profiles of model S1.c are displayed up to the radius where a shock forms at t  ≈  7.31 × 10⁵ s and begins to expand into the low-density outer layers of the SMS.

Zoom In Zoom Out Reset image size

The evolutionary tracks for the central density and temperature of the rotating models R1.a and R1.d are also shown in Figure 4. A dashed line corresponds to model R1.a, with an initial metallicity Z_CNO = 5 × 10⁻⁴, which collapses to a BH. A solid line denotes model R1.d with Z_CNO = 2 × 10⁻³, which explodes due to the energy released by the hot CNO cycle. We find that model R1.c, with a lower metallicity of Z_CNO = 1 × 10⁻³, also explodes when the central temperature is in the range dominated by the hot CNO cycle.

As a result of the kinetic energy stored in the rotation of models R1.c and R1.d, the critical metallicity needed to trigger an explosion decreases significantly relative to the non-rotating case. We observe that rotating models with initial metallicities up to Z_CNO = 8 × 10⁻⁴ do not explode even via the rp-process, which is dominant at temperatures above T  ≈  5 × 10⁸ K and increases the hydrogen burning rate by 200–300 times relative to the hot CNO cycle. We also note that the evolution timescales of the collapse and bounce phases are reduced because rotating models are more compact and have a higher initial central density and temperature than the spherical ones at the onset of the gravitational instability.

The right panels in Figure 4 show the time evolution of the central rest-mass density (upper panel) and the central temperature (lower panel) for models S2.a, S2.b, R2.a, and R2.d, i.e., of models with a mass of M = 10⁶ M_☉. We find that the critical metallicity for an explosion in the spherical case is Z_CNO = 5 × 10⁻² (model S2.b), while model S2.a with Z_CNO = 3 × 10⁻² collapses to a BH. We note that the critical metallicity leading to a thermonuclear explosion is higher than the critical value found by Fuller et al. (1986, i.e., Z_CNO = 1 × 10⁻²) for a spherical SMS with the same mass. The initial metallicity leading to an explosion in the rotating case (model R2.d) is more than an order of magnitude smaller than in the spherical case. As for the models with a smaller gravitational mass, the thermal bounce takes place when the physical conditions in the central region of the star allow for the release of energy by hydrogen burning through the hot CNO cycle. Overall, the dynamics of the more massive models indicate that the critical initial metallicity required to produce an explosion increases with the rest mass of the star.

Figure 6 shows the total nuclear energy generation rate in erg s⁻¹ for the exploding models as a function of time during the late stages of the collapse just before and after bounce. The main contribution to the nuclear energy generation is due to hydrogen burning by the hot CNO cycle. The peak values of the energy generation rate at bounce lie between several 10⁵¹ erg s⁻¹ for the rotating models (R1.d and R2.d) and ≈10⁵²–10⁵³ erg s⁻¹ for the spherical models (S1.c and S2.b). As expected, the maximum nuclear energy generation rate needed to produce an explosion is lower in the rotating models. Moreover, as the explosions are due to the energy release by hydrogen burning via the hot CNO cycle, the ejecta would mostly be composed of ⁴He.

Zoom In Zoom Out Reset image size

As a result of the thermal bounce, the kinetic energy rises until most of the energy of the explosion is in the form of kinetic energy. We list in the second to last column of Table 1 the radial kinetic energy after thermal bounce, which ranges between E_RK = 1.0 × 10⁵⁵ erg for the rotating star R1.c and E_RK = 3.5 × 10⁵⁷ erg for the spherical star S2.b.

Due to the lack of resolution at the surface of the star, it becomes difficult to compute accurately the photosphere and its effective temperature from the criterion that the optical depth is τ = 2/3. Therefore, in order to estimate the photon luminosity produced in association with the thermonuclear explosion, we make use of the fact that within the diffusion approximation the radiation flux is given by

$$\begin{equation} {F_{\gamma }} =-\frac{c}{3 \kappa _{\rm es}\rho }\nabla {U}, \end{equation} \tag{ 35 }$$

where U is the energy density of the radiation and κ_es is the opacity due to electron Thompson scattering, which is the main source of opacity in SMSs. The photon luminosity in terms of the temperature gradient and for the spherically symmetric case can be written as

$$\begin{equation} {L_{\gamma }} =-\frac{16 \pi a c r^2 T^3}{3 \kappa _{es}\rho }\frac{\partial {T}}{\partial r}, \end{equation} \tag{ 36 }$$

where a is the radiation constant and c is the speed of light. As can be seen in the last panel of Figure 7, the distribution of matter becomes spherically symmetric during the phase of expansion after the thermonuclear explosion. In this figure (Figure 7) we show the isodensity contours for the rotating model R1.d. The frames have been taken at the initial time (left panel), at t = 0.83 × 10⁵ s (central panel) just after the thermal bounce (at t = 0.78 × 10⁵ s), and at t = 2.0 × 10⁵ s when the radius of the expanding matter is roughly four times the radius of the star at the onset of the collapse.

Zoom In Zoom Out Reset image size

The photon luminosity computed using Equation (36) for model R1.d is displayed in Figure 8, where we also indicate with a dashed vertical line the time at which the thermal bounce takes place. The photon luminosity before the thermal bounce is computed at radii inside the star unaffected by the local dynamics of the low-density outer layers, which is caused by the initial pressure perturbation and by the interaction between the surface of the SMS and the artificial atmosphere. Once the expanding shock forms near the surface, the photon luminosity is computed near the surface of the star. The light curve shows that, during the initial phase, the luminosity is roughly equal to the Eddington luminosity ≈5 × 10⁴³ erg s⁻¹ until the thermal bounce. Then, the photon luminosity becomes super-Eddington when the expanding shock reaches the outer layers of the star and reaches a value of L_γ  ≈  1 × 10⁴⁵ erg s⁻¹. This value of the photon luminosity after the bounce is within a few percent difference with respect to the photon luminosity Fuller et al. (1986) found for a non-rotating SMS of the same rest mass. The photon luminosity remains super-Eddington during the phase of rapid expansion that follows the thermal bounce. We compute the photon luminosity until the surface of the star reaches the outer boundary of the computational domain ≈1.0 × 10⁵ after the bounce. Beyond that point, the luminosity is expected to decrease and then rise to a plateau of ∼10⁴⁵ erg s⁻¹ due to the recombination of hydrogen (see Fuller et al. 1986 for a non-rotating star).

Zoom In Zoom Out Reset image size

The outcome of the evolution of models that do not generate enough nuclear energy during the contraction phase to halt the collapse is the formation of a BH. The evolutionary tracks for the central density and temperature of some of these models are also shown in Figure 4. The central density typically increases up to ρ_c  ∼  10⁷ g cm⁻³ and the central temperature up to T_c  ∼  10¹⁰ K just before the formation of an AH.

Three isodensity contours for the rotating model R1.a collapsing to a BH are shown in Figure 9, which displays the flattening of the star as the collapse proceeds. The frames have been taken at the initial time (left panel), at t = 0.83 × 10⁵ s (central panel) approximately when model R1.d with higher metallicity experiences a thermal bounce, and at t = 1.127 × 10⁵s, where a BH has already formed and its AH has a mass of 50% of the total initial mass.

Zoom In Zoom Out Reset image size

At the temperatures reached during the late stages of the gravitational collapse (in fact at T  ⩾  5 × 10⁸ K), the most efficient process for hydrogen burning is the breakout from the hot CNO cycle via the ¹⁵O(α, γ)¹⁹Ne reaction. Nevertheless, we find that models, which do not release enough nuclear energy by the hot CNO cycle to halt their collapse to a BH, are not able to produce a thermal explosion due to the energy liberated by the ¹⁵O(α, γ)¹⁹Ne reaction. We note that above 10⁹ K, not all the liberated energy is used to increase the temperature and pressure, but it is partially used to create the rest mass of the electron–positron pairs. As a result of pair creation, the adiabatic index of the star decreases, which means the stability of the star is reduced. Moreover, due to the presence of e^± pairs, neutrino energy losses grow dramatically.

Figure 10 shows (solid lines) the time evolution of the redshifted neutrino luminosities of four models collapsing to a BH (S1.a, R1.a, S2.a, and R2.a) and (dashed lines) of four models experiencing a thermal bounce (S1.c, R1.d, S2.b, and R2.d). The change of the slope of the neutrino luminosities at ∼10⁴³ erg s⁻¹ denotes the transition from photo-neutrino emission to the pair-annihilation-dominated region. The peak luminosities in all forms of neutrino for models collapsing to a BH are L_ν  ∼  10⁵⁵ erg s⁻¹. Neutrino luminosities can be that important because the densities in the core prior to BH formation are ρ_c  ∼  10⁷ g cm⁻³, and therefore neutrinos can escape. The peak neutrino luminosities lie between the luminosities found by Linke et al. (2001) for the collapse of spherical SMSs and those found by Woosley et al. (1986, who only took into account the luminosity in the form of electron antineutrino). The maximum luminosity decreases slightly as the rest mass of the initial model increases, which was already observed by Linke et al. (2001). In addition, we find that the peak of the redshifted neutrino luminosity does not seem to be very sensitive to the initial rotation rate of the star. We also note that the luminosity of model R1.a reflects the effects of hydrogen burning at L_ν  ∼  10⁴³ erg s⁻¹.

Zoom In Zoom Out Reset image size

The total energy output in the form of neutrinos is listed in the last column of Table 1 for several models. The total radiated energies vary between E_ν  ∼  10⁵⁶ erg for models collapsing to a BH and E_ν  ∼  10⁴⁵–10⁴⁶ erg for exploding models. These results are in reasonable agreement with previous calculations. For instance, Woosley et al. (1986) obtained that the total energy output in the form of electron antineutrinos for a spherical SMS with a mass 5 × 10⁵ M_☉ and zero initial metallicity was 2.6 × 10⁵⁶ erg, although their simulations neglected general relativistic effects that are important to compute accurately the relativistic redshifts. On the other hand, Linke et al. (2001), by means of relativistic one-dimensional simulations, found a total radiated energy in the form of neutrinos of about 3 × 10⁵⁶ erg for the same initial model and about 1 × 10⁵⁶ erg when redshifts were taken into account. In order to compare with the results of Linke et al. (2001), we computed the redshifted total energy output for model S1.a, having the same rest mass, until approximately the same evolution stage as Linke et al. (2001) did (i.e., when the differential neutrino luminosity dL_ν/dr  ∼  4 × 10⁴⁵ erg s⁻¹ cm⁻¹). We find that the total energy released in neutrinos is 1.1 × 10⁵⁶ erg.

The neutrino luminosities for models experiencing a thermonuclear explosion (dashed lines in Figure 10) peak at much lower values L_ν  ∼  10⁴²–10⁴³ erg s⁻¹ and decrease due to the expansion and disruption of the star after the bounce.

The axisymmetric gravitational collapse of rotating SMSs with uniform rotation is expected to emit a burst of GWs (Saijo et al. 2002; Saijo & Hawke 2009) with a frequency within the LISA low-frequency band (10⁻⁴–10⁻¹ Hz). Although through the simulations presented here we could not investigate the development of non-axisymmetric features in our axisymmetric models that could also lead to the emission of GWs, Saijo & Hawke (2009) have shown that the three-dimensional collapse of rotating stars proceeds in an approximately axisymmetric manner.

In an axisymmetric spacetime, the ×-mode vanishes and the +-mode of GWs with l = 2 computed using the quadrupole formula is written as (Shibata & Sekiguchi 2003)

$$\begin{equation} h^{{\rm quad}}_{+} = \frac{\ddot{I}_{xx}(t_{{\rm ret}})-\ddot{I}_{zz}(t_{{\rm ret}})}{r}\sin ^2\theta, \end{equation} \tag{ 37 }$$

where $\ddot{I}_{ij}$ refers to the second time derivative of the quadrupole moment. The GW quadrupole amplitude is $A_2(t)=\ddot{I}_{xx}(t_{{\rm ret}})-\ddot{I}_{zz}(t_{{\rm ret}})$. Following Shibata & Sekiguchi (2003), we compute the second time derivative of the quadrupole moment by finite-differencing the numerical results for the first time derivative of I_ij obtained by

$$\begin{equation} \dot{I}_{ij}= \int \rho _{*}(v^ix^j+x^iv^j)d^3x. \end{equation} \tag{ 38 }$$

We calculate the characteristic GW strain (Flanagan & Hughes 1998) as

$$\begin{equation} h_{{\rm char}}(f) = \sqrt{\frac{2}{\pi }\frac{G}{c^3}\frac{1}{D^2}\frac{dE(f)}{df}}, \end{equation} \tag{ 39 }$$

where D is the distance of the source and dE(f)/df is the spectral energy density of the gravitational radiation given by

$$\begin{equation} \frac{dE(f)}{df} = \frac{c^3}{G}\frac{(2\pi f)^2}{16 \pi }|\tilde{A}_2(f)|^2, \end{equation} \tag{ 40 }$$

with

$$\begin{equation} \tilde{A}_2(f) = \int A_2(t)e^{2\pi ift}dt. \end{equation} \tag{ 41 }$$

We have calculated the quadrupole GW emission for the rotating model R1.a collapsing to a BH. We plot in Figure 11 the characteristic GW strain (Equation (39)) for this model assuming that the source is located at a distance of 50 Gpc (i.e., z  ≈  11), together with the design noise spectrum $h(f)=\sqrt{f S_h(f)}$ of the LISA detector (Larson et al. 2000). We find that, in agreement with Saijo et al. (2002), Saijo & Hawke (2009), and Fryer & New (2011), the burst of GWs due to the collapse of a rotating SMS could be detected at a distance of 50 Gpc and at a frequency that approximately takes the form (Saijo et al. 2002)

$$\begin{equation} f_{{\rm burst}}\sim 3\times 10^{-3} \left(\frac{10^{6}\,M_{\odot }}{M}\right)\left(\frac{5M}{R}\right)^{3/2} {\rm [Hz]}, \end{equation} \tag{ 42 }$$

where R/M is a characteristic mean radius during BH formation (typically set to R/M = 5).

Zoom In Zoom Out Reset image size

Furthermore, Kiuchi et al. (2011) have recently investigated, by means of three-dimensional general relativistic numerical simulations of equilibrium tori orbiting BHs, the development of the non-axisymmetric Papaloizou–Pringle instability (PPI) in such systems (Papaloizou & Pringle 1984) and have found that a non-axisymmetric instability associated with the m = 1 mode grows for a wide range of self-gravitating tori orbiting BHs, leading to the emission of quasi-periodic GWs. In particular, Kiuchi et al. (2011) have pointed out that the emission of quasi-periodic GWs from the torus resulting after the formation of an SMBH via the collapse of an SMS could be well above the noise sensitivity curve of LISA for sources located at a distance of 10 Gpc. Such instability appears for tori whose angular velocity in the equatorial plane expressed as $\bar{\Omega }(r)\propto r^{q}$ has q < q_kep, where q_kep corresponds to the Keplerian limit, i.e., q = −1.5 in Newtonian gravity.

We find that the torus (defined as the rest mass outside the AH) that forms after the collapse to a BH of the uniformly rotating model R1.a (when the mass of the AH exceeds 50% of the gravitational mass) does not fulfill the above condition for the development of the PPI. However, we find that the torus that forms when the differentially rotating model D1 collapses to a BH has a distribution of angular momentum such that $\bar{\Omega }(r)\propto r^{q}$ with q  ≈  − 1.62. This suggests that the torus may be prone to the development of the non-axisymmetric PPI, which would lead to the emission of quasi-periodic GWs with peak amplitude ∼10⁻¹⁸–10⁻¹⁹ and frequency ∼10⁻³ Hz during an accretion timescale ∼10⁵ s.

We have presented results of general relativistic simulations of collapsing supermassive stars using the two-dimensional general relativistic numerical code Nada, which solves the Einstein equations written in the BSSN formalism and the general relativistic hydrodynamic equations with high-resolution shock-capturing schemes. These numerical simulations have used an EOS that includes the effects of gas pressure and tabulated those associated with radiation pressure and electron–positron pairs. We have also taken into account the effects of thermonuclear energy release by hydrogen and helium burning. In particular, we have investigated the effects of hydrogen burning by the β-limited hot CNO cycle and its breakout via the ¹⁵O(α, γ)¹⁹Ne reaction (rp-process) on the gravitational collapse of non-rotating and rotating SMSs with non-zero metallicity.

We have presented a comparison with previous studies and investigated the influence of the EOS on the collapse. We emphasize that axisymmetric calculations without rotation (i.e., models S1 and S2) retain the spherical symmetry of the initial configurations as there are no physical phenomena to produce asphericity and numerical artifacts associated with the use of Cartesian coordinates are negligibly small. Overall, our collapse simulations yield good agreement with previous works when using the same treatment of physics. We have also found that the collapse timescale depends on the ion contributions to the EOS, and electron–positron pair creation affects the stability of SMSs. Interestingly, differentially rotating stars that are gravitationally stable with a Γ = 4/3 EOS can become unstable against gravitational collapse when the calculation is performed with the microphysical EOS including pair creation.

We have found that objects with a mass of ≈5 × 10⁵ M_☉ and an initial metallicity greater than Z_CNO  ≈  0.007 explode if non-rotating, while the threshold metallicity for an explosion is reduced to Z_CNO  ≈  0.001 for objects that are uniformly rotating. The critical initial metallicity for a thermal explosion increases for stars with a mass of ≈10⁶ M_☉. The most important contribution to the nuclear energy generation is due to the hot CNO cycle. The peak values of the nuclear energy generation rate at bounce range from ∼10⁵¹ erg s⁻¹ for rotating models (R1.d and R2.d) to ∼10⁵²–10⁵³ erg s⁻¹ for spherical models (S1.c and S2.b). After the thermal bounce, the radial kinetic energy of the explosion rises until most of the energy is kinetic, with values ranging from E_K  ∼  10⁵⁶ erg for rotating stars up to E_K  ∼  10⁵⁷ erg for the spherical star S2.b. The neutrino luminosities for models experiencing a thermal bounce peak at L_ν  ∼  10⁴² erg s⁻¹.

The photon luminosity roughly equals the Eddington luminosity during the initial phase of contraction. Then, after the thermal bounce, the photon luminosity becomes super-Eddington with a value of about L_γ  ≈  1 × 10⁴⁵ erg s⁻¹ during the phase of rapid expansion that follows the thermal bounce. For those stars that do not explode we have followed the evolution beyond the phase of BH formation and computed the neutrino energy loss. The peak neutrino luminosities are L_ν  ∼  10⁵⁵ erg s⁻¹.

SMSs with masses less than ≈10⁶ M_☉ could have formed in massive halos with T_vir  ≳  10⁴ K. Although the amount of metals that was present in such environments at the time when SMSs might have formed is unclear, it seems possible that the metallicities could have been smaller than the critical metallicities required to reverse the gravitational collapse of an SMS into an explosion. If so, the final fate of the gravitational collapse of rotating SMSs would be the formation of a SMBH and a torus. In a follow-up paper, we aim to investigate in detail the dynamics of such systems (collapsing of an SMS to a BH-torus system) in three dimensions, focusing on the post-BH evolution, and the development of non-axisymmetric features that could emit detectable gravitational radiation.